# Removing Climatic Overprints in Sedimentary Cosmogenic Beryllium Records: Potentials and Limits

Tatiana Savranskaia, Ramon Egli, Quentin Simon, Jean-pierre Valet, Franck

Bassinot, Nicolas Thouveny

► To cite this version:

Tatiana Savranskaia, Ramon Egli, Quentin Simon, Jean-pierre Valet, Franck Bassinot, et al. Re-moving Climatic Overprints in Sedimentary Cosmogenic Beryllium Records: Potentials and Limits. Geochemistry, Geophysics, Geosystems, 2024, 25 (12), pp.e2024GC011761. 10.1029/2024gc011761. hal-04842285

HAL Id: hal-04842285

https://hal.science/hal-04842285v1

Submitted on 17 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2024GC011761

**Key Points:**

* A new technique for removing climatic overprints from geomagnetic records is presented

* A worst-case example based on a single cosmogenic $^{10}$Be record yields a reduction of environmental overprints by a factor of 2

* Cosmogenic $^{10}$Be is a better geomagnetic field proxy for the North Atlantic Ocean than $^{10}$Be normalized by $^{9}$Be

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

T. Savranskaia,
tatiana.savranskaia@lmu.de

**Citation:**

Savranskaia, T., Egli, R., Simon, Q., Valet, J.-P., Bassinot, F., & Thouveny, N. (2024). Removing climatic overprints in sedimentary cosmogenic beryllium records: Potentials and limits.

*Geochemistry, Geophysics, Geosystems*, 25, e2024GC011761. https://doi.org/10.1029/2024GC011761

Received 24 JUL 2024

Accepted 16 OCT 2024

**Author Contributions:**

**Conceptualization:** Tatiana Savranskaia

**Formal analysis:** Tatiana Savranskaia, Ramon Egli

**Funding acquisition:** Jean-Pierre Valet

**Investigation:** Tatiana Savranskaia, Ramon Egli, Quentin Simon, Franck Bassinot

**Methodology:** Tatiana Savranskaia, Ramon Egli

**Resources:** Quentin Simon, Jean-Pierre Valet, Franck Bassinot, Nicolas Thouveny

**Validation:** Ramon Egli

© 2024 The Author(s). Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Removing Climatic Overprints in Sedimentary Cosmogenic Beryllium Records: Potentials and Limits

Tatiana Savranskaia¹,² ID, Ramon Egli³,⁴, Quentin Simon⁵ ID, Jean-Pierre Valet³ ID, Franck Bassinot⁶, and Nicolas Thouveny⁵

¹GeoForschungsZentrum (GFZ) Potsdam, Potsdam, Germany, ²Now at Department of Earth and Environmental Sciences, Ludwig Maximilians University, Munich, Germany, ³Institut de Physique du Globe de Paris, CNRS, Université de Paris, Paris, France, ⁴Department of Geophysics, GeoSphere Austria, Vienna, Austria, ⁵CEREGE UM34, CNRS, IRD, INRAE, Coll France, Aix Marseille University, Aix en Provence, France, ⁶Laboratoire des Sciences Du Climat et de l'Environnement (CEA-CNRS-UVSQ), Domaine Du CNRS, Gif-sur-Yvette, France

**Abstract** Continuous reconstructions of past variations of the Earth's magnetic field are based mainly on paleomagnetic and cosmogenic $^{10}$Be records in marine sediments. In both cases, the recording mechanisms can be affected by environmental processes. Climatic overprints are only partially removed by normalization procedures, so that stacking is used to further remove site-specific effects. Regionally or globally correlated artifacts, however, cannot be removed by stacking. Here we present a modified approach where geomagnetic records are complemented by environmental proxies representing processes that might affect the field recording mechanism. Geomagnetic and environmental records are jointly processed with principal component analysis to obtain a set of components supposed to represent true variations of the geomagnetic field and climatic overprints, respectively. After discussing the theoretical background of this new approach and its underlying assumptions, a practical example is presented, using a worst-case scenario based on a single $^{10}$Be record from the North Atlantic with strong climatic overprints, covering the last 600 ka. The first two principal components, which represent the modulation of $^{10}$Be by global climatic variations and by the geomagnetic field, respectively, explain 66.3% of the signal variance. Comparison of the geomagnetic principal component with global relative paleointensity stacks shows that the original climatic overprint can be reduced by a factor of 2, outperforming a $^{10}$Be/$^{9}$Be stack obtained from two sites with little glacial-interglacial variability. The proposed method for removing climatic overprints can be applied to multiple sites more efficiently than conventional stacking.

**Plain Language Summary** Continuous records of the Earth's magnetic field rely on measurements of magnetic minerals or cosmogenic isotopes in sediments. Both types of records are also sensitive to environmental conditions and are thus affected by past climatic variations. These unwanted climatic overprints are difficult to remove: one strategy consists in stacking records from different sites; however, regionally or globally correlated artifacts cannot be completely removed by this technique. Here we present a new method for separating the geomagnetic signal from unwanted climatic overprint, which is based on the principal component analysis (PCA). The efficiency of this new method is tested with a worst-case example based on a single site located in the North Atlantic, which is characterized by strong glacial-interglacial variability. The first two principal components obtained from PCA represent the modulation of $^{10}$Be by global climatic variations and by the geomagnetic field, respectively. Comparisons of the geomagnetic field component with reference data show that the original climatic overprint has been reduced by a factor of 2, outperforming a $^{10}$Be/$^{9}$Be stack obtained from two sites with little glacial-interglacial variability. The proposed method for removing climatic overprints can be applied to multiple sites more efficiently than conventional stacking.

## 1. Introduction

Sedimentary paleomagnetic records depend on the partial alignment of magnetic minerals with the Earth's magnetic field during and after sediment deposition. This alignment produces a natural remanent magnetization (NRM) which is proportional to the local field intensity and to the concentration of magnetic carriers in the sediment. Changes in magnetic mineral concentrations are compensated by normalizing the NRM with a laboratory-induced remanent magnetization, typically the anhysteretic remanent magnetization (ARM) or the isothermal remanent magnetization (IRM). Ideally, the RPI record resulting from this normalization procedure is proportional to the local field vector intensity. In reality, ARM and IRM display different sensitivities to the size,
**Visualization:** Tatiana Savranskaia,
Ramon Egli

**Writing – original draft:**
Tatiana Savranskaia

**Writing – review & editing:**

Tatiana Savranskaia, Ramon Egli, Jean-Pierre Valet, Nicolas Thouveny

shape, and composition of the remanence carriers, none of which is exactly proportional to the original magnetic moments, so that normalization artifacts can occur in case of downcore variations of the magnetic properties (Roberts et al., 2013). The magnetic alignment mechanism during NRM acquisition also depends on the interplay between sediment lithology and the size and shape of remanence carriers (e.g., Franke et al., 2004). Consequently, RPI reconstructions tend to be affected by environmental factors such as change of deep oceanic currents circulation (Channell et al., 1998; Kissel et al., 2009; Xuan & Channell, 2008) or enhanced fluvial contributions on the continental margins (Moreno et al., 2002). Such environmental contaminations are particularly evident at the level of glacial-interglacial cycles.

The cosmogenic isotope $^{10}$Be is produced in the stratosphere and the upper troposphere by cosmic rays through spallation reactions involving collisions of high-energy cosmic particles with nitrogen and oxygen atoms (Lal & Peters, 1967). The $^{10}$Be production rate is primarily controlled by the Earth's magnetic field through its shielding action against incoming charged particles, and, to a lesser extent, by variations of the incoming flux of primary galactic cosmic rays and by the solar activity (Beer et al., 2012; Blinov, 1988; Kocharov et al., 1989). Accordingly, $^{10}$Be production is maximized during periods of low geomagnetic dipole intensity, such as during geomagnetic excursions and reversals. Atmospheric $^{10}$Be is eventually removed by aerosols (Heikkilä et al., 2008) before being deposited on the continents and into the oceans, where it is scavenged by suspended matter, entering finally into the sedimentary record (von Blanckenburg & Bouchez, 2014). The lithological dependence of the $^{10}$Be scavenging efficiency is eliminated by normalizing the concentration of $^{10}$Be in sediment with the authigenic fraction of the terrestrial isotope $^{9}$Be (Bourlès et al., 1989) released by rock weathering (Measures & Edmond, 1983; von Blanckenburg et al., 2012). Both isotopes are adsorbed to the same type of sediment particles with the same efficiency (Sharma et al., 1987), so that $^{10}$Be/$^{9}$Be is free of artifacts caused by variations of the scavenging efficiency. However, normalization by $^{9}$Be can introduce additional environmental overprints to the geomagnetic response of $^{10}$Be, owing to the different source location, transport path and sedimentation cycling of the two isotopes (McHargue & Donahue, 2005; Simon, Thouveny, Bourles, Nuttin, et al., 2016; von Blanckenburg & Bouchez, 2014; von Blanckenburg & Igel, 1999; von Blanckenburg et al., 2012).

In summary, because of their distinct recording mechanisms and sensitivity to global and local geomagnetic field features, RPI and $^{10}$Be records represent independent proxies for the local and global field intensity, respectively, providing different constraints to global paleomagnetic field reconstructions (e.g., Korte et al., 2011; Nilsson et al., 2014; Ziegler et al., 2011). First comparisons of cosmogenic $^{10}$Be with relative paleointensity stacks revealed strong similarities, as expected from their common geomagnetic modulation, but also significant temporal offsets during geomagnetic dipole lows when considering individual cores (Carcaillet et al., 2004; Thouveny et al., 2008). Contamination of these records by environmental factors thus represents a persistent problem, which is only partially eliminated by stacking different sites. While site-specific environmental overprints can be averaged out by stacking distributed records from different basins, global climatic artifacts, such as those induced by glacial-interglacial cycles, cannot be removed in this way. Significant improvements to the fidelity of geomagnetic records can be brought by a method that is capable of separating the geomagnetic signal from environmental overprints using the limited number of available RPI stacks and $^{10}$Be records.

Here we present a new method for removing paleoclimatic overprints from RPI and $^{10}$Be records using information derived from climatic proxies. For this purpose, we test the worst-case scenario of a single high-latitude $^{10}$Be record from the North Atlantic site MD95-2016 covering the last 600 ka, in which glacial-interglacial cycles represent the main environmental forcing. Principal component analysis (PCA) is used to isolate $^{10}$Be variations due to the geomagnetic field and to global and regional climate changes, respectively, based on authigenic Be isotope measurements integrated with X-ray fluorescence (XRF) elemental ratios, planktic $\delta^{18}$O, and magnetic susceptibility. Residual artifacts in the isolated geomagnetic component of the $^{10}$Be record are shown to be smaller than those of a global $^{10}$Be/$^{9}$Be stack (by Simon, Thouveny, Bourlès, Valet, et al. (2016)), demonstrating the power of the proposed approach.

## 2. Sampling Site, Lithology, and Environmental Setting

The sedimentary core MD95-2016 (57°42.46′N; 29°25.44′W) was retrieved in the North Atlantic Ocean at the Reykijanes Ridges (2,318 m water depth) during the IMAGES MD101 cruise of R/V Marion Dufresne (Figure 1). The site is located in the Gardar Drift area (L. D. Labeyrie, 1996), which is a flat contourite sheeted drift deposited
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Bathymetric map showing the location of studied site MD95-2016, along with the nearby ODP sites 982 and ODP 983 (Jansen et al., 1996). The red and blue thick lines depict the pathways of North Atlantic Deep Water and the Iceland Scotland Overflow Water, respectively.</figcaption></figure>

by a broad current (Faugères et al., 1999). The intermediate and deep waters flowing at the site, which stands on the pathways of the North Atlantic Current (NAC) (Petit et al., 2018), include the Labrador Sea Water, formed by open ocean convection in the Labrador Sea, the warm North East Atlantic water, the Lower Deep Water (Müller-Michaelis & Uenzelmann-Neben, 2014), and the Iceland Scotland Overflow Water (ISOW), respectively. The sediment of core MD95-2016 is mainly composed of clayey/sandy nanno ooze and silt. The sedimentation rate is exceptionally high over the last 100ka, reaching up to 22 cm/ka (Figure S1e in Supporting Information S1). Data from two neighboring cores ODP 982 (57°52′N; 15°87′W) and 983 (60°40′N; 23°64′W), retrieved during the Ocean Drilling Program campaigns at 1,135 and 1,194 m water depth, respectively, were included in the present study for comparison purposes.

## 3. Measurements

### 3.1. Beryllium Isotopes

Beryllium measurements of marine sediments are performed on the so-called authigenic fraction, which contains ⁹Be and ¹⁰Be atoms scavenged from the water column by adsorption onto suspended matter (Bourlès et al., 1989). About 50% of authigenic ¹⁰Be and 30% of authigenic ⁹Be are carried by iron and manganese oxyhydroxides, ~20% and ~10% by the exchangeable and carbonate fractions, respectively, and the remaining part by organic matter (Bourlès et al., 1989). Core MD95-2016 was sampled every 10 cm for beryllium measurements. A total of 247 samples covering the topmost 25.2 m (0–615 ka, see Section 3.3 for the age model), were processed for Be isotope analyses. For each sample, ~1 g of dry sediment was chemically treated for authigenic Be extraction at the CEREGE National Cosmogenic Nuclides Laboratory (France), following the protocol summarized in Simon, Thouveny, Bourlès, Valet, et al. (2016). Authigenic ⁹Be was measured using a graphite-furnace atomic absorption spectrophotometer with double-beam correction. Authigenic ¹⁰Be was measured after chemical preparation at the French accelerator mass spectrometer (AMS) national facility ASTER. Chemistry blank ratios range from 10⁻¹⁴ to 10⁻¹⁵, which is at least 3 orders of magnitude lower than ¹⁰Be/⁹Be ratios in the sediment. Concentration values of ¹⁰Be were corrected for radioactive decay using a half-life of 1.387 ±0.012 Ma (Chmeleff et al., 2010; Korschinek et al., 2010).

### 3.2. XRF Scanning

Elemental composition records provide essential information about sediment sources, transport, primary productivity, and atmospheric/oceanic circulation (Calvert & Pedersen, 2007). Elemental ratios have been obtained from high-resolution micro-core XRF spectrometry using the XRF core scanner AVAATECH (EPOC, Université de Bordeaux). The XRF spectrum is derived from the fluorescence emitted by electrons excited by a primary X-ray source, in form of photon counts per photon energy interval (Oyedotun, 2018; Richter et al., 2006). Each element is characterized by a series of peaks and the integrated intensity is proportional to element concentration. The procedure for converting raw XRF outputs into absolute elemental concentrations requires non-trivial calibration procedures affected by large uncertainties (Weltje & Tjallingii, 2008). For this reason, results are usually expressed as elemental ratios. Elemental ratios are widely used to track multiple provenance sources (Richter et al., 2006). XRF measurements of 15 elements were conducted at 1 cm intervals on U-Channels using a 10 keV source with 0.5 mA current. A total of 2004 measurements covering the topmost 25.2 m have been obtained (Figure 2). Correlation coefficients between elemental ratios and δ¹⁸O, ¹⁰Be, and ¹⁰Be/⁹Be are listed in Table S1 of Supporting Information S1.

### 3.3. Oxygen Isotope and Age Model

Stable oxygen isotope measurements of core MD95-2016 were performed at ~10 cm intervals on picks of the planktonic foraminifera *Globigerina bulloides* in the 250–315 µm size fraction, using VG-Optima and Elementar Isoprime dual-inlet mass spectrometers at the Laboratoire des Sciences du Climat et de l’Environnement (LSCE, Gif-Sur-Yvette, France). All results are expressed as δ¹⁸O with respect to the V-PDB standard. The external
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Downcore variations of Be records, δ<sup>18</sup>O, X-ray fluorescence elemental ratios, and magnetic susceptibility. Shaded bands highlight interval of elevated <sup>9</sup>Be concentrations. Notice that some scales are reversed to maintain a positive correlation between curves.</figcaption></figure>

analytical reproducibility of ±0.05‰ (1σ) was determined from replicate measurements of a laboratory carbonate standard. The standard approach for depth-age model construction relies on fitting the δ<sup>18</sup>O record to a reference, commonly a stack curve, and following astronomical tuning if more precise orbitally paced signals are available. The global benthic δ<sup>18</sup>O stack LR04 (Lisiecki & Raymo, 2005) is the most widely used target curve because it integrates 57 globally distributed benthic records spanning the past 5.3 Ma. The HMM stack used here is a recent update of LR04 based on a hidden Markov model that includes a larger number of records (Ahn et al., 2017). Tuning a single-location δ<sup>18</sup>O record to the global δ<sup>18</sup>O stack might bias the age-depth model due to the incomplete global synchronization of δ<sup>18</sup>O variability (L. Labeyrie et al., 2005; Waelbroeck et al., 2011). Age lags were addressed by constructing regional benthic δ<sup>18</sup>O stacks covering the past 150 ka (Lisiecki & Stern, 2016). Regional stacks enable a better synchronization of geographically close records.

The age model was constructed in three steps. First, we defined age tie points by matching the peaks of Ice Rafted Debris (IRD) events in core MD95-2016 with those of U1308 (Hodell et al., 2008), using Si/Sr and Zr/Sr as IRD proxies (Figure S1 in Supporting Information S1). The choice of U1308 as a target core is determined by its accurate age model being constrained by radiocarbon dates and oxygen isotopes stratigraphy with Greenland ice cores (Svensson et al., 2006) for the last 42 ka, and with the benthic stack LR04 prior to 60 ka. Planktonic δ<sup>18</sup>O records can be affected by significant age offsets caused by foraminifera dissolution (e.g., Barker et al., 2007), and size segregation (Savranskaia et al., 2022) effects. Indeed, offsets by up to 10 ka affects the δ<sup>18</sup>O based age model of MD95-2016 during the last ~100 ka, as seen for instance in the <sup>10</sup>Be record of the Laschamps excursion (Figure S1 and Table S3 in Supporting Information S1). In order to correct these offsets, additional tie points were obtained by matching the <sup>9</sup>Be record of MD95-2016 with the Deep North Atlantic (DNA) δ<sup>18</sup>O stack (Lisiecki & Stern, 2016) (Figure S10 in Supporting Information S1). Authigenic <sup>9</sup>Be is a valid substitute for δ<sup>18</sup>O at this high-latitude site, owing to its strong association with glacial-interglacial cycles and the relatively short residence time in the Atlantic Ocean (von Blanckenburg & O'nions, 1999).

In the second and third steps, we aligned the upper part of the MD95-2016 δ<sup>18</sup>O record with the DNA δ<sup>18</sup>O age model. Due to the limited DNA age coverage prior to 0–135 ka, the older part of sedimentary record was matched
with the HMM stack using the software package BIGMACS, which is based on the Bayesian Inference Gaussian Process regression (Lee et al., 2022). The accuracy of the automated correlation by BIGMACS has been further evaluated by comparing the MD95-2016 record with the planktic δ¹⁸O variations from the neighbor core ODP982 (57°30.8′N, 152°55.5′W) on its independent age scale (Venz et al., 1999). There is an overall good agreement between the planktic δ¹⁸O records of MD95-2016 and ODP982, respectively, with two notable exceptions around 210–250 and 380–430 ka (Figure S1 in Supporting Information S1). Climatic and geomagnetic components extracted from the beryllium records (Section 5) appear to be well synchronized with global paleointensity and δ¹⁸O stacks over the entire time interval, supporting the validity of the chosen age model.

### 3.4. Magnetic Measurements and Lightness

NRM, low-field magnetic susceptibility ($\chi$) and anhysteretic remanent magnetization (ARM) of core MD95-2016 were measured on U-channels and single samples using an Agico KLY-3 susceptibility meter and a 2G Enterprises cryogenic magnetometer, inside the shielded room of the paleomagnetic laboratory at the Institut de Physique du Globe de Paris (Supporting Information S1). The Sediment lightness (L*) variations were quantified with a Minolta CV-2002 hand-held spectrophotometer during the Marion Dufresne cruise. Measurements were performed at 5 cm intervals along the whole core.

The ratio ARM/$\chi$ is a proxy for the grain size of low-coercivity minerals such as magnetite and Ti-magnetite, owing to the sensitivity of $\chi$ and ARM to the total magnetic mineral content and its finest fraction, respectively (King et al., 1982). The downcore variability of $\chi$ and ARM/$\chi$ has frequently been used to trace sediment provenance (Kissel, 2005), deep oceanic circulation strength (Ballini et al., 2006), the dominant transport path of magnetic particles (Kissel et al., 1999), and changes in the balance between detrital and authigenic minerals (Watkins & Maher, 2003) in different parts of the North Atlantic Ocean. Large ARM/$\chi$ variations in MD95-2016 are driven by glacial-interglacial cycles (Figure S3 in Supporting Information S1), as seen by the presence of large peaks corresponding to δ¹⁸O minima, making this core unsuited for RPI determinations. Indeed, the RPI record reconstructed from (NRM/ARM)$_{50-80 mT}$ is dominated by large-amplitude oscillations unrelated to the geomagnetic field (Figure S2 in Supporting Information S1). The RPI stacks SINT-2000 (Valet et al., 2005) and PADM2M (Ziegler et al., 2011) have been used instead for comparison with cosmogenic Be records, since they are less affected by local field variability, which is not recorded by cosmogenic isotopes, and by climatic overprints.

## 4. Methods and Models

### 4.1. Model for Cosmogenic Records of the Earth's Magnetic Field

The global ¹⁰Be production rate P is related to the global field intensity and to the galactic cosmic ray flux by the empirical expression

$$ \frac{P}{P_0} \approx Q\left(\frac{m}{m_0}\right)\Phi\left(\frac{\phi}{\phi_0}\right) \qquad (1) $$

derived from cosmogenic isotope production calculations based on a dipole model of the geomagnetic field (e.g., Poluianov et al., 2016). In this expression, $P_0 \approx 1.48 \times 10^9$ atoms/s is the modern value of P, m is the geomagnetic dipole moment with modern value $m_0 = 78$ ZAm², and $\phi$ the modulation potential of the cosmic rays, with a mean value of $\phi_0 = 650$ MV. Q and $\Phi$ are empirical functions expressing the dependence of P on the geomagnetic field intensity and the cosmic ray flux, respectively. Data from Lal (1992) and Masarik and Beer (2009) are well approximated by $Q(x) = ((\xi + 1)/(\xi + x))^q$ with best-fit parameters $\xi = 0.31$ and $q = 0.50$ (Figure 3). P is inversely related to the geomagnetic dipole moment, with maximum values corresponding to periods of low field intensity during excursions or reversals. The non-linearity of Q makes ¹⁰Be records particularly sensitive to low-field conditions and relatively insensitive to secular variations under normal field intensity conditions (Figures 3b and 3c).

Cosmogenic sedimentary records of the Earth's magnetic field can be written in their most general form as

$$ R(t) = A + B \cdot S(t), \qquad (2) $$
<figure><img src="image_5.png" /></figure>

**Figure 3.** Relation between the dipole moment of the Earth magnetic field and the cosmogenic $^{10}$Be production rate. (a) Dipole moment *m* normalized by its present-day value *m*₀ versus $^{10}$Be production rate *P*, normalized by the present-day value *P*₀. Blue dots: calculated production rates from Lal (1992), replotted by Ménabréaz et al. (2012), and from the integration of latitudinal data by Masarik and Beer (2009). Solid line: empirical field modulation function *Q*(*x*) = ((*ξ* + 1)/(*ξ* + *x*))⁽q* with *ξ* = 0.31 and *q* = 0.50. Brown dots with error bars: data obtained from the comparison of $^{10}$Be/$^{9}$Be records with geomagnetic field models (Simon et al., 2018; Simon, Thouveny, Bourlès, Valet, et al., 2016). (b) Simulated geomagnetic event (dipole moment drop during an excursion or a reversal), superimposed to a sinusoidal signal representing short-term field variations with constant amplitude. (c) $^{10}$Be production rate obtained using to the dipole variation in (b) as input of *Q*. Notice the enhanced sensitivity to small field variations during low-field conditions. (d) Same long-term variation of the $^{10}$Be production rate as in (c), combined with a sinusoidal climatic modulation after Equation 1. (e) Dipole moment variations reconstructed from the $^{10}$Be record in (d), using $Q^{-1}$. Notice the amplification of the climatic modulation during periods of regular field strength.

where *R* is the recorded parameter in sediment (e.g., $^{10}$Be/$^{9}$Be), *S* a function that describes the primary source of the recorded signal (e.g., the global production rate *P* in case of $^{10}$Be), *B* a proportionality factor, and *A* an offset accounting for field-independent sources of the recorded signal. In the case of Be records, *B* describes the efficiency with which $^{10}$Be is transported to the site and scavenged by the depositing sediment, while *A* > 0 represents unsynchronized Be inputs originating from $^{10}$Be reservoirs with long residence time, which might affect sites with >3,000 m water depth. This term is negligible for MD95-2016 (Savranskaia et al., 2021).

Ideal cosmogenic records are characterized by $A = 0$ and a site-specific, time-independent proportionality coefficient $B$, so that $R$ in Equation 2 can be converted to a global production rate, and, assuming a constant flux of cosmic rays, to the geomagnetic dipole moment. In the case of cosmogenic Be records, this conversion is based on $m/m_0 \propto Q^{-1}(R/R_0)$, with $R_0 = R(t_0)$ for a reference epoch $t_0$ where the dipole field intensity was close to the modern value $m_0$. In reality, $B$ is affected by environmental factors and is therefore time-dependent. The time dependency of $B$ introduces a climatic modulation of the geomagnetic signal. In the case of marine sediments, $B = T \cdot K$, where $K$ is the sediment scavenging efficiency, and $T$ a transport term expressing the time-integrated concentration of dissolved Be encountered along the settling path of suspended sediment—$^{10}$Be adsorbed on land practically does not contribute to the final authigenic concentration (Wittmann et al., 2017). The use of isotopic ratios (e.g., $^{10}$Be/$^{9}$Be) sets $K = 1$, thus eliminating the effect of scavenging efficiency variations, provided that both isotopes are characterized by the same $T$. This simplification, however, is not necessarily correct given the different geographic distribution of $^{10}$Be and $^{9}$Be sources, which typically generate distinct source-site connecting paths.

Equation 2 with time-dependent coefficients $A$ and $B$ has been used to model $^{10}$Be/$^{9}$Be records of the Matyuama-Brunhes reversals in five marine cores, including MD95-2016 (Savranskaia et al., 2021). Differences between records from the Indian and Pacific oceans relative to MD95-2016 were explained by a combination of (a) site-specific, almost time-invariant coefficients $A$ and (b) site-specific time-dependent coefficients $B(t)$ whose variations turned out to be well correlated with the $\delta^{18}$O record. These observations led to the conclusion that climatic
contaminations of ¹⁰Be/⁹Be records arise from a multiplicative modulation mechanism represented by the product of a field-dependent and an environment-dependent term in Equation 2. During periods of normal field intensity, small climatic contaminations of the ¹⁰Be record introduce large errors in geomagnetic moment reconstructions obtained from Q⁻¹(P/P₀), owing to the non-linearity of Q (Figures 3d and 3e). These errors can easily exceed the typical amplitude of secular variation.

## 4.2. Theory of Geomagnetic and Climatic Signal Separation

Principal Component Analysis (PCA) is a dimensionality reduction technique commonly used to decrease signal complexity and extract relevant information (Jolliffe, 2002). For this purpose, a given set of signals is represented as linear combination of orthogonal functions called principal components (PC). The PCs correspond to directions of maximum variance of the data in the original measurement space after successive subtraction of each PC from all signals. PCs are mathematical constructs which might or might not coincide with physically meaningful signals. In the context of Earth sciences, PCA can be used as a stacking analog to isolate a global signal from a set of geographically distributed records (e.g., Spratt & Lisiecki, 2016). In this case, the extracted global signal is an optimally weighted average of the individual records. PCA can therefore be used to reconstruct past variations of the global field intensity. In the case of RPI records, this approach sacrifices detailed field geometry information to remove noise, overprints, and other unwanted signal variations from the data. This is not expected to occur with cosmogenic records, which are intrinsically a global proxy of the dipole field intensity. PCA has been applied to ¹⁰Be records from the Atlantic Ocean after normalization with ²³⁹Th (Christl et al., 2010). The different affinities of Be and Th to specific minerals in the sediment make Th-normalized ¹⁰Be records particularly sensitive to lithological changes, and therefore susceptible to containing significant environmental overprints. In this case, however, PCA did not lead to the identification of climate-related components (Christl et al., 2010).

Equation 2 also describes RPI records: in this case, S(t) is the ratio between the NRM of sediment and a laboratory magnetization supposed to represent variations of the magnetic mineralogy, while A = 0 by definition. B represents the product of (a) the alignment efficiency of magnetic minerals during NRM acquisition and (b) the ratio between the magnetic moments of remanence carriers in their NRM state and in the magnetic state induced by the laboratory magnetization, respectively (Hofmann & Fabian, 2009). Conventional stacking is commonly used to reduce RPI overprints at a regional or global scale. A single study adopted a different approach to eliminate climatic overprints, using PCA on a set of records which included RPI and a set of environmental parameters retrieved from the same core (Valet et al., 2011). This approach exploits additional information contained in environmental proxies, which cannot be incorporated by conventional stacking techniques. As shown here, the joint PCA analysis of geomagnetic records and climatic proxies can be applied even to cases where only a single geomagnetic record is available.

The successful use of PCA for the purpose of isolating physical signals rely on the fulfillment of the following conditions (Jolliffe, 2002):

1. *Linearity.* The analyzed records are linear combinations of a common set of signals (e.g., geomagnetic and climatic variations), superimposed to record-specific variations which, in the context of PCA, are formally equivalent to noise.

2. *Orthogonality.* The signal of interest must be orthogonal to the other signals (i.e., pairwise correlations between this signal and all other signals are zero).

3. *Variability.* Time samples from the common set of signals build Gaussian-like distributions with distinct variances.

None of these conditions are fully met in real applications. Condition 1 is obviously not met in the context of climatically modulated records of the geomagnetic field, because of the multiplicative term in Equation 2. Condition 2 holds in the case of sufficiently long geomagnetic records, if climatic variations are statistically uncorrelated to the Earth's magnetic field. Connections between climate and magnetic field variations, for example, through Milankovic orbital cycles, have been proposed, but evidence is tenuous and highly debated (Courtillot et al., 2007). A recent claim (Zhou et al., 2023) was based on the separation of paleomagnetic and climatic signals in ¹⁰Be records from the Chinese Loess Plateau (Zhou et al., 2007). A critical aspect of all correlation analyses based on sedimentary records is represented by the cross-contamination of supposedly
independent proxies, notably, climate and RPI, as we show in this study. Condition 2 is also violated by coincident events, regardless of causality. For instance, a number of paleomagnetic excursions during the Brunhes chron (0–780 ka) was noted to occur during minima in precession and warmer interglacial episodes (Carcaillet et al., 2004). Finally, condition 3 is met to a good degree of approximation in the case of records with unimodal statistics.

Additional considerations are thus required to establish whether, and under which conditions, $B(t)$ and $S(t)$ in Equation 2 can be separated by PCA. We begin by expressing the climatic modulation as $B(t) = b_0 + b_1 \delta B(t)$ with site-specific coefficients $b_0$ and $b_1$, and a site-independent function $\delta B(t)$ with zero mean (i.e., $\langle \delta B(t) \rangle = 0$) representing global changes of the field recording efficiency. In this case, Equation 2 becomes

$$R(t) = b_0 S(t) + b_1 S(t) \delta B(t) + \rho(t), \quad (3)$$

where $\rho(t)$ is a noise term that includes $A(t)$ and other record-specific inputs. In Equation 3, $S(t)$ is separable from the other terms only if $\langle S(t) \rho(t) \rangle = 0$ and $\langle S(t) \delta B(t) \rangle = 0$ (condition 2). The latter criterion is satisfied for sufficiently long time series, provided that geomagnetic and climatic variations are driven by completely independent processes. If both criteria are met, the geomagnetic component $S(t)$ is captured by a corresponding principal component $PC_{geo}$, while a second principal component reproduces the mixed global signal $S(t) \delta B(t)$ corresponding to the second term of Equation 3. Additional PCs represent record-specific inputs contained in $\rho(t)$.

The finite duration of the analyzed records introduces random violations of the orthogonality condition, adding uncertainty to the PCA results. In mathematical terms, residual contaminations of $PC_{geo}$ are caused by a rotation of the orthogonal set of PCs away from their correct orientation in the measurement space. Accordingly, the inclusion of climatic/environmental proxies play an important role in fixing the direction of PCs that supposedly represents the time variations of $B(t)$. Climatic proxies can be expected to be linearly related to $B(t)$, but not to $S(t) \delta B(t)$. The peculiar properties of geomagnetic field variations, however, allow to overcome this problem. This can be shown by reformulating Equation 3 in a form that provides the desired linearity. For this purpose, we use the long-term average $S_0 = \langle S(t) \rangle$ of the geomagnetic signal to write $S(t) = S_0 + \Delta S(t)$ and obtain

$$R(t) = b_0 S_0 + b_0 \Delta S(t) + b_1 S_0 \delta B(t) + b_1 \Delta S(t) \delta B(t) + \rho(t). \quad (4)$$

Now, $b_1 S_0$ can be replaced by a new constant coefficient, while the offset $b_0 S_0$, which corresponds to $\langle R(t) \rangle$ by definition, is irrelevant to the PCA. Accordingly, the only mixed term remaining in Equation 4 is $b_1 \Delta S(t) \delta B(t)$. This term is the product of two quantities that represent deviations from the long-term means of $S(t)$ and $B(t)$. In the case of good-quality paleomagnetic records, climatic contaminations can be assumed to be relatively small, so that $b_1 \delta B(t) \ll b_0$. In order for the mixed term to be negligible, the condition $\Delta S(t) \ll S_0$ must hold as well. This condition is fulfilled during periods of stable field conditions—when $S_0$ and $\Delta S(t)$ correspond to the mean geomagnetic signal and its changes due to secular variation, respectively—but not during field excursions and reversals, when the cosmogenic production rate, and thus $\Delta S(t)$, almost doubles (Figure 3). Nevertheless, the limited duration of significant geomagnetic events ensures that $\Delta S(t) \ll S_0$ for most of the samples taken from sufficiently long records. The remaining samples have a limited effect on covariance matrices and thus on PCA results. Accordingly, geomagnetic ($R$) and climatic ($X$) records take the desired linear form

$$R(t) = \langle R \rangle + b_0 \Delta S(t) + b_1 \delta B(t) + \rho(t), \quad (5)$$

$$X(t) = \langle X \rangle + x_1 \delta B(t) + \chi(t),$$

where $\rho(t)$ is a residual signal containing also the linearization errors associated with invalid assumption $\Delta S(t) \ll S_0$ during excursions and reversals.

### 4.3. PCA Implementation

In order to test the efficiency of the proposed PCA approach in removing climatic contamination from geomagnetic records, we provide an example based on the worst-case scenario of a single geomagnetic record—namely $^{10}$Be or $^{10}$Be/$^{9}$Be from core MD95-2016. Additional records used to constrain climatic effects are obtained from environmental proxies measured on the same core. These proxies include $^{9}$Be, $\delta^{18}$O, elemental ratios,
and the magnetic susceptibility $\chi$. Significant reduction of the climatic overprint of $^{10}$Be in this example represents an important proof of principle for future applications based on multiple sites.

In order to deal with highly contrasting variations of signals characterized by different units, PCA has been performed on the standardized records

$$r(t) = \frac{R - \langle R \rangle}{\text{sd}(R)} \quad \text{and} \quad x_i(t) = \frac{X_i - \langle X_i \rangle}{\text{sd}(X_i)} \qquad (6)$$

PCA analysis of $n$ signals yields $n$ principal components $PC_1, \dots, PC_n$ in order of decreasing contributions to the signal variance. One of these PCs, designated as $PC_g = PC_{\text{geo}}$, is expected to represent common geomagnetic variations in the corresponding records (in our example, a single record). At least one of the remaining components carries a common climate signature, whereby the exact number of PCs needed to represent relevant climatic variations depends on the number of independent signals contained in the climate proxies. Accordingly, the global part of the standardized geomagnetic measurements is given by

$$\tilde{r}(t) = c_g PC_g(t) + \sum_{\substack{i=1 \\ i \neq g}}^{l<n} c_i PC_i(t), \qquad (7)$$

where the first and second terms on the right-hand side represents the field-driven component of geomagnetic record and its climatic modulation, respectively. Combination of Equations 6 and 7 yields the geomagnetic component

$$S_g(t) = \langle R \rangle - A + c_g \text{sd}(R) PC_g(t) \qquad (8)$$

of the record $R$, where $A$ is the (unknown) offset defined in Equation 2. Determination of $A$ and a proper normalization is required to transform $S_g$ into the desired field intensity reconstruction. In the case where $R$ is a cosmogenic Be record, this last step is performed by solving

$$\frac{P(t)}{P_0} = b_0 [\langle R \rangle - A + c_g \text{sd}(R) PC_g(t)] = Q\left(\frac{m(t)}{m_0}\right) \qquad (9)$$

with respect to the offset $A$ and the scaling parameter $b_0$. The two main assumptions underlying this solution are that (a) PCA properly identifies the geomagnetic component with a single PC, and that (b) cosmic ray flux variations are negligible during the time interval covered by the records. Assumption (b) might not be valid: for instance, variations by up to a factor 1.5 over the last 35 ka have been detected from in-situ produced $^{14}$C in an ice core from Greenland (Lal et al., 2005). Equation 9 contains two unknown coefficients, $A$ and $b_0$, which need to be determined before solving with respect to $m/m_0$. The determination of $A$ and $b_0$ requires a minimum of two calibration points corresponding to epochs with different geomagnetic dipole strengths obtained from a reference paleomagnetic model. Two-point calibrations are prone to errors originating from the uncertainty of the reference model and from possible age offsets.

A more robust approach to the determination of the unknown coefficients consists in using the solution of Equation 9 with respect to $PC_g$ as model function for a least squares regression of an existing geomagnetic model. Best-fit $A$ and $b_0$ values tend to be strongly correlated, and therefore poorly determined, if the curvature of $Q$ over the range of dipole moments expected during the time interval under consideration is small compared to the scatter of $PC_g$. This is clearly seen in the limit case where $Q$, and therefore also $PC_g$, are nearly linear, in which case Equation 9 yields infinite combinations of $A$ and $b_0$ values with same fitting performance. If $A = 0$ can be assumed, as in the specific case of MD95-2016, a stable solution can always be obtained for $b_0$, in analogy with the usual normalization procedure for cosmogenic Be records (e.g., Simon, Thouveny, Bourlès, Valet, et al., 2016). Otherwise, solutions for $A$ and $b_0$ can be stabilized by adding a regularization term to the mean squared residuals to be minimized, such that solutions with unrealistically large values of $A$ are penalized. Taking the residuals $r$ of a regression with $A = 0$ as a measure of the typical error of the reconstructed record, and, at the same time, as a threshold above which $A$ becomes significant, the function to minimize for obtaining $A$ and $b_0$ is given by
$$ \frac{1}{Nc_g^2 \text{var}(R)} \sum_{i=1}^{N} \left[ \frac{1}{b_0} Q \left( \frac{m_i}{m_0} \right) + A - \langle R \rangle - c_g \text{sd}(R) \text{PC}_g(t_i) \right]^2 + \frac{(b_0 A)^2}{\langle r^2 \rangle} \quad (10) $$

As shown later, Equation 10 yields very small values of A when applied to the Be record of MD95-2016 and the corresponding PCA results.

## 4.4. Uncertainty Estimation

Violations of the three conditions required to retrieve the geomagnetic signal from a set of geomagnetic and environmental records (Section 4.2) add uncertainties to the PCA results. These uncertainties can be regarded as the result of the stochastic nature of input data (Görtler et al., 2020), where the particular set of records used for PCA represents a single realization of a multivariate probability distribution. In our case, two main stochastic sources contribute to the uncertainty of PCA results: (a) orthogonality errors (e.g., $\langle S(t)B(t) \rangle \neq 0$) caused by the finite duration of the records, and (b) site- or proxy-specific contributions to the records ($\rho(t)$ and $\chi(t)$ in Equation 5). The first error source cannot be accounted for, as $S(t)$ and $B(t)$ are unknown. The second error source can be estimated with a bootstrap approach if site- or proxy-specific contributions can be assumed to be identically distributed (Babamoradi et al., 2013). Bootstrap methods require a large number of input records and cannot be used in our application. Therefore, we opted for a Monte Carlo approach where an a-posteriori error estimate is assigned to each record.

The a-posteriori error estimate $\hat{E}(t)$ of a record $X(t)$ is identified with the difference between $X(t)$ and its reconstruction based on the PCs assigned to $S(t)$ and $B(t)$. In the context of the above mentioned stochastic approach, $\hat{E}(t)$ is a random realization of a stochastic vector $E(t)$ whose autocovariance function is approximated by $E_{XX}(\tau) = \langle X(t)X(t + \tau) \rangle$. The autocovariance function can be used to generate random copies of the original error estimate. These copies are noise signals with power spectrum $E_{XX}^*(\nu)$ defined by the Fourier transform of $E_{XX}$. Accordingly, the inverse Fourier transform of

$$ E_X^*(\nu) = \sqrt{E_{XX}^*(\nu)} e^{2\pi i \epsilon(\nu)}, \quad (11) $$

is a random realization of $E(t)$, with $\epsilon(\nu)$ being an independent random number comprised between 0 and 1 for each discrete frequency $\nu$ (Timmer & König, 1995). The Monte Carlo simulation of PCA errors consists in generating copies of the original records through the addition of random realizations of the a-posteriori error estimates. Each set of copies is used as a new PCA input, obtaining PCs that can be identified with realizations of the (unknown) stochastic functions that build the input data. It is thereby implicitly assumed that errors from different records are uncorrelated. This assumption cannot be verified and might not be completely valid. However, avoidance of redundant records based on proxies that are obviously correlated (e.g., $^{10}$Be and $^{10}$Be/$^{9}$Be) is expected to minimize the impact of correlated errors.

# 5. Results and Discussion

## 5.1. Authigenic $^{10}$Be and $^{9}$Be Records in the North Atlantic

Authigenic $^{10}$Be concentrations in MD95-2016 are comprised between $2.55 \times 10^8$ and $12.45 \times 10^8$ at·g$^{-1}$ (Figure 4a). The average ($6.11 \times 10^8$ at·g$^{-1}$) is ~3 times larger than that of sediments from the Central Baffin Bay during the last 140 ka (Simon, Thouveny, Bourlès, Valet, et al., 2016), ~1.6 times larger than that of a record from the Central Arctic Ocean (Frank et al., 2008), and compatible with the range observed for the Arctic and Norwegian Seas (Eisenhauer et al., 1994). A comparison with the neighbor site ODP 983 is possible only over a short time interval during the Icelandic Basin geomagnetic excursion (~190 ka): in this case, the authigenic $^{10}$Be concentration in MD95-2016 is ~3.6 times larger than in ODP 983 (Knudsen et al., 2008).

Authigenic $^{9}$Be concentrations in MD95-2016 range from $1.02 \times 10^{16}$ to $3.43 \times 10^{16}$ at·g$^{-1}$ (Figure 4b). These values are twice as large as for a Central Arctic record (Frank et al., 2008), slightly higher than the mean of the central Baffin Bay and ~10 times larger than at the neighbor ODP site 983. The $^{9}$Be concentration is largest during glacial marine isotopic stages, exceeding the average concentration by a factor ~1.5. The effect of glacial-interglacial cycles on $^{9}$Be is clearly demonstrated by the positive correlation with Ti/Ca ($r = +0.65$, Figure 2),
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Authigenic Be records. (a, b) Measured $^{10}$Be and $^{9}$Be concentrations ("Be$_{meas}$"), long term linear trends of the measured data ("trend"), and detrended measurements ("Be$_{d}$"). (c) Authigenic $^{10}$Be/$^{9}$Be (gray, left scale) and the global RPI stack SINT-2000 of Valet et al. (2005) converted to a synthetic $^{10}$Be production rate P/P$_{0}$ using the function Q defined in Section 4.1 (orange, right scale). Vertical gray bands highlight periodic maxima in the $^{9}$Be record.</figcaption></figure>

which reflects the low affinity of Be toward calcium carbonate (Chase et al., 2002). The correlation between Ca and $^{10}$Be is weaker than with $^{9}$Be, as expected from the fact that $^{10}$Be variations are mainly driven by the geomagnetically controlled production rate, but in some intervals (e.g., 105–150, 242–300, 324–393 ka) $^{10}$Be appears to be negatively correlated with the Ca content, contrary to $^{9}$Be. This might be explained by a dilution effect caused by high fluxes of IRD, since the terrigenous material transported by IRD is poor in $^{10}$Be and rich in $^{9}$Be.

Authigenic $^{10}$Be/$^{9}$Be is comprised between $1.32 \times 10^{-8}$ and $4.46 \times 10^{-8}$, with a mean of $2.92 \times 10^{-8}$ (Figure 4c). This mean is lower than the ($4-6) \times 10^{-8}$ range reported for DNA waters (von Blanckenburg et al., 1996), and ($5-8) \times 10^{-8}$ for deep Arctic waters (Frank et al., 2009). The overall smaller $^{10}$Be/$^{9}$Be values in the North Atlantic, relative to the Pacific, are caused by larger terrigenous inputs, as indicated by the correlation with $^{9}$Be ($r^2$ = -0.42), Zr/Sr ($r^2$ = -0.53) and Ti/Ca ($r^2$ = -0.47). Despite evident climatic contaminations of the $^{10}$Be/$^{9}$Be record, several peaks (e.g., 100 ka, 190 ka, 380, and 550 ka) agree with $^{10}$Be maxima during major geomagnetic events recorded by RPI (Figure 4c). As discussed in Section 4.1, the climatic modulation of cosmogenic records is smallest during such events, while it easily masks secular variations during periods of normal field intensity (Figures 3d and 3e). This explains the overall poor agreement between $^{10}$Be/$^{9}$Be and RPI.

A long-term concentration trend is observed for both Be isotopes. This trend implies an almost linear increase of authigenic Be concentrations over the last 600 ka (Figures 4a and 4b). Absence of a similar trend in $^{10}$Be/$^{9}$Be means that both isotopes are equally affected (Figure 4c), excluding explanations calling for a faulty $^{10}$Be, as well as long-term trends affecting the strength of the geomagnetic field or the cosmic ray flux. Alternate explanations include (a) the remobilization of weakly adsorbed authigenic Be inside the sedimentary column and subsequent
expulsion during compaction, (b) a systematic increase of the Be scavenging efficiency over time, (c) an increase of eolian inputs due to changes in the atmospheric circulation, and (d) a long-term change in deep water circulation controlling delivery of Be to the site. Authigenic Be remobilization has been shown to be an important source contributing to the global oceanic beryllium budget (Deng et al., 2023; Savranskaia et al., 2021). Expulsion of Be-rich pore water leads to the formation of a deep-water Be pool that does not represent atmospheric and continental inputs. Scavenging of this pool by depositing sediment adds a nearly constant contribution to the Be record, yielding $A > 0$ in Equation 1 (Savranskaia et al., 2021). This contribution has been detected only at sites with water depths $>2.5$ km, where poorly mixed bottom waters are characterized by $^{10}$Be residence times $>1$ ka (von Blanckenburg & Igel, 1999). Its absence in core MD95-2016 implies that there is no significant accumulation of benthic beryllium released into the water column, without excluding diagenetic losses.

A systematic increase of the sediment scavenging efficiency over time can be excluded, as there is no long-term change in sediment lithology that is compatible with the observed trend. Long-term trends of the eolian input or a systematic change of the atmospheric circulation pattern during the last 600 ka have not been reported. In contrast, an increased supply of terrigenous material to the North Atlantic Ocean during the last 2.5 Ma was inferred from changes of $^{206}$Pb/$^{204}$Pb and $^{143}$Nd/$^{144}$Nd in ferromanganese crusts (Frank, 2002; von Blanckenburg & Nägler, 2001). This trend is compatible with an intensification of the North Atlantic Deep Water current within the same period of time, due to the onset of Quaternary glaciations (Raymo, 1994). Because Pb and Nd isotopes cannot be used to distinguish intensified erosion/weathering from changes of deep water circulation (Burton et al., 1997), measurements of $^{10}$Be/$^{9}$Be in ferromanganese crusts from the northwest Atlantic have been used to differentiate between the two processes (von Blanckenburg & O'nions, 1999). Ferromanganese $^{10}$Be/$^{9}$Be in the North Atlantic and Pacific Oceans did not change significantly over the last 7 Ma; however, this result does not preclude the possibility that similar relative changes in the concentration of the two isotopes occurred over the same time interval.

Because the presence of artifacts in long-term trends in the input data can bias the results of PCA analysis, the $^{10}$Be and $^{9}$Be records have been detrended prior to further analyses. The correction has been performed by multiplying $^{10}$Be and $^{9}$Be concentrations with $1/(1 - ct)$, where $1 - ct$, with constant coefficient c, is the empirical linear trend determined by a least squares fit of the data shown in Figure 4. The same value of c was used for $^{10}$Be and for $^{9}$Be.

## 5.2. Glacial-Interglacial Cycles

Bulk geochemical analyses (XRF) allow us to interpret past environmental changes and their impact on Be records in terms of (a) mineralogical, (b) grain size, and (c) sediment provenance variations. The most significant changes in XRF elemental ratios are correlated with intervals of elevated $^{9}$Be concentrations occurring during glacial marine isotopic stages (Figures 2 and 5). Glacial stages are generally characterized by a reduced carbonate content and an increased terrigenous flux into the Atlantic Ocean (Balsam & McCoy, 1987). The evolution of carbonate concentration is mirrored by Ca/Ti—a proxy for the ratio between biogenic (Ca) and detrital (Ti) inputs, and by Ca/Sr—a proxy for the ratio between biogenic (Ca) and detrital (Sr) carbonates.

The evident 100 ka periodicity of $^{9}$Be and $^{10}$Be/$^{9}$Be records correlates with variations of Zr/Sr (Figure 5a), which is an indirect proxy of ice-rafted debris (e.g., Hodell et al., 2010; Margari et al., 2020). IRD material was mainly deposited within the primary IRD belt given by the North Atlantic surface water circulation (Ruddiman et al., 1970). The evident correspondence between Zr/Sr (Figure 5a) and IRD records (Figure 5b) from two neighbor ODP sites 982 (Venz et al., 1999) and 983 (Barker et al., 2015) confirms that the MD95-2016 site is located within the IRD accumulation zone (Venz et al., 1999). Zr/Sr and IRD peaks during glacial periods indicate that detrital material is delivered to the site via massive glacier calving as the ice front reached its maximum seaward position, rather than during glacier melting (Bond et al., 1992; McManus et al., 1999). IRD layers are poor in detrital carbonate and rich in silicates (Hodell et al., 2008), as seen by Si/Sr (Figure 5c). Authigenic $^{9}$Be concentration maxima during periods of high silicate content and IRD inputs (Figures 5b–5d) point to a common detrital source, with authigenic $^{9}$Be released during ice melting being transported by oceanic currents (Measures & Edmond, 1983; von Blanckenburg et al., 2012). Large oceanic currents are thus expected to be the main factor controlling the $^{9}$Be concentration in sediment.

Some studies suggest that the North Atlantic thermohaline circulation did not change significantly across glacial-interglacial cycles (Foster et al., 2007; von Blanckenburg & Nägler, 2001), while recent findings revealed a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between beryllium records and environmental proxies. (a) Zr/Sr as proxy for bulk sediment grain size (left scale, dotted gray line), and global stack HMM of benthic $\delta^{18}$O (solid blue line, right scale). (b) Ice Rafted Debris (IRD) (predominantly by quartz and volcanic material) in ODP cores 982 and 983. (c) Detrended $^{10}$Be (dotted black line, left scale), and Si/Sr as proxy of IRD (solid gray line, right scale). (d) Detrended $^{9}$Be (dotted black line, left scale) and Ti/Ca, as proxy of the ratio between detrital (Ti) and biogenic (Ca) inputs (solid gray line, right scale). Vertical gray bands highlight periodic maxima in the $^{9}$Be record.</figcaption></figure>

persistent 100 ka periodicity of $^{143}$Nd/$^{144}$Nd, which reflects an increased inflow of water masses from the Southern Ocean into the DNA during glacial periods and weaker overturning circulation (Kim et al., 2021). Another source of glacial-interglacial cyclicity of $^{9}$Be is related to the scavenging efficiency of sediment, which is sensitive to mineral size and composition (Boschi & Willenbring, 2021; Bourlès et al., 1989). Clay content, which is representative of fine sediment fraction, is reflected by Rb/Sr. While there is strong enrichment of Rb in mica and detrital clay minerals (Croudace & Rothwell, 2015; Kylander et al., 2011), Sr is often associated with the presence of aragonite, which forms in shallow waters (Frijia & Parente, 2008). The strong positive correlation between $^{9}$Be and Ti/Ca and Rb/Sr (Figure 2, Table S1 in Supporting Information S1) suggests a modulation of the scavenging efficiency by sediment grain size. This is also confirmed by the decreasing Ti (Ti/Fe) content during glacial periods (Figure 2), which indicates sediment fining, and thus a more efficient Be adsorption (Simon, Thouveny, Bourles, Nuttin, et al., 2016).

Environmental factors specifically affecting $^{10}$Be, and thus also $^{10}$Be/$^{9}$Be, are more difficult to assess. The source distribution of $^{10}$Be is different from that of $^{9}$Be, yielding distinct responses to changes of the oceanic circulation
pattern. Site-specific processes affecting the ¹⁰Be concentration in MD95-2016 might include inputs of ¹⁰Be accumulated on ice sheets during glacial times and subsequently released by melting (Frank et al., 2008; Simon, Thouveny, Bourles, Nuttin, et al., 2016). The lack of ¹⁰Be peaks corresponding to IRD events (Figure 5c) suggests that melting of ice masses does not deliver ¹⁰Be isotopes to the MD95-2016 sediment to the same extent as terrestrial ⁹Be.

## 5.3. PCA

In order to disentangle the geomagnetic and environmental components of the Be record of core MD95-2016, we performed a PCA analysis according to the method discussed in Section 4, using the following 10 records: δ¹⁸O, detrended ⁹Be and ¹⁰Be (Section 5.1), χ, and the logarithms of Si/Sr, Ca/Sr, Sr/Rb, Ti/Ca, Ti/K, and Zr/Sr ratios of elemental XRF intensities, which are linearly related to concentration ratios. All records were standardized using Equation 6. The first two principal components, PC₁ and PC₂, explain 49.6% and 16.7% of the total signal variance, respectively (Figure S4 in Supporting Information S1). PC₁ evidently reflects detrital input variations, as seen by the good (r² > 0.9, Table S2 in Supporting Information S1) correlations with Ti/Ca (detrital input), Zr/Sr (bulk sediment grain size) and Rb/Sr (weathering). These variations mainly reflect glacial-interglacial cycles (Figure 6a, r² = 0.65), and correlate well with ⁹Be (r² = 0.77). Due to the multiplicity of climatic proxies used as inputs for PCA, PC₁ (Figure 6a) explains the largest variance of the data set.

PC₂ is related mainly to ¹⁰Be (r² = -0.78), explaining 75% of its variance: it is therefore identified with PCgeo. The global ¹⁰Be productivity obtained from PCgeo by Equations 9 and 10, referred to as ¹⁰Begeo in the following, is overall coherent with the composite global record of relative paleointensity SINT-2000 (r² = -0.51, Figure 6b), and, to a much lesser extent, with PISO-1500 (r² = -0.29, Figure 6c), despite the contrary might be expected from the fact that PISO-1500 is based on sites from the same oceanic basin as MD95-2016. Because SINT-2000 was not used as input for PCA analysis, this result can only be explained by assuming that PISO-1500 contains significant climatic overprints. PCgeo contributes also significantly to the variances of χ and Ti/K (48% and 43%, respectively): this is unexpected, since geomagnetic components should not contribute to climatic proxies if variations of the Earth's magnetic field (ΔS in Equation 5) and of climatic parameters (δB in Equation 5) are not correlated. Possible reasons are discussed in Section 5.5.

PC₃ and higher components explain the remaining 33.7% of the total signal variance, with individual contributions not exceeding 10%. The main climatic and geomagnetic signals are therefore identified with PC₁ and PC₂, respectively. The remaining principal components have been treated as noise and used to estimate the uncertainty of ¹⁰Begeo using the Monte Carlo approach described in Section 4.4. The confidence intervals of ¹⁰Begeo have been defined by the 25% and 75% quantiles of the ensemble of PCA solutions obtained from 100 Monte Carlo simulations of the input records. The contribution of PC₃ to the ¹⁰Be record (Figure S5a in Supporting Information S1) is mainly reflected by the sea surface water temperature (SST), as indicated by the good match with a SST record from the neighbor core ODP 982 (Herbert et al., 2016). The SST record of ODP 982 is driven by global and regional processes over short (e.g., direct forcing by insolation cycles, changes in the surface wind fields, growth and decay of sea ice) and long (e.g., growth of continental ice sheets) time scales (Lawrence et al., 2009). The total contribution of PC₄ and PC₅ to the ¹⁰Be record does not show obvious relations with climatic proxies (Table S2 in Supporting Information S1), nor with the geomagnetic field, except for the 390–470 ka interval (Figure S5b in Supporting Information S1), where variations of PC₄+₅ match the ¹⁰Be production rate expected from SINT-2000. This means that either SINT-2000 contains a climatic contamination that is captured by PC₄ and PC₅, or that part of the geomagnetic modulation of ¹⁰Be leaked into PC₄ and PC₅.

The majority of geomagnetic excursions seen on SINT and PISO stacks are characterized by a significant decrease of the field intensity. Accordingly, most of the known geomagnetic events during the last 600 ka are reflected by corresponding ¹⁰Begeo peaks (Figure 6). Apparent time offsets with respect to RPI stacks are mostly limited to 1–4 ka, except the Calabrian Ridge II event with 9 ka (Table S3 in Supporting Information S1). These offsets lie within the uncertainty of the age model (Figure S3 in Supporting Information S1); nevertheless, RPI lags ¹⁰Begeo in most cases. This is expected in sediments carrying a post-depositional remanent magnetization, since this type of magnetization is locked below the typical mixing depth of bioturbated sediment (Egli & Zhao, 2015). Some excursions are missing in ¹⁰Begeo, while other peaks (e.g., at ~170 ka) are not reflected by RPI records. Major discrepancies between ¹⁰Begeo and SINT-2000, in excess of the PCA confidence intervals, occur at
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Principal component analysis results. (a) Rescaled climate-driven variations of $^{10}$Be obtained from PC$_{1}$ (left), and HMM benthic $\delta^{18}$O stack (right, Ahn et al., 2017). Numbers correspond to marine isotopic stages. (b, c) Rescaled geomagnetically driven variations of $^{10}$Be obtained from PC$_{2}$ (left, error bars represent confidence intervals), compared to the paleointensity stacks SINT-2000 (Valet et al., 2005) and PISO-1500 (Channell et al., 2009) converted to $P/P_{0}$ (right). (d) Comparison between $^{10}$Be$_{geo}$ and the stacks SINT-2000, PISO-1500, and $^{10}$Be/$^{9}$Be from Simon, Thouveny, Bourlès, Valet, et al. (2016). Identified geomagnetic events are based on Roberts (2008), Laj and Channell (2015), and Simon, Thouveny, Bourlès, Valet, et al. (2016).</figcaption></figure>

four discrete time intervals: 150–190 ka, 280–320 ka, 350–420 ka, and 520–550 ka. These intervals are not systematically related to specific climatic conditions, but, as discussed in Section 5.5, differences between $^{10}$Be$_{geo}$ and RPI stacks correlate with the magnetic susceptibility record of MD95-2016.

5.4. Assessing the Effectiveness of PCA

The improvement of the global field reconstruction obtained with our PCA analysis, relative to original $^{10}$Be and $^{10}$Be/$^{9}$Be records, is evident when analyzing scatter plots of $P/P_{0}$ versus $m/m_{0}$ obtained from a reference model, chosen to be SINT-2000 (Figure 7). In these plots, the congruence between $^{10}$Be- and RPI-based records is measured by the data scatter around the theoretical relation $P/P_{0} = Q(m/m_{0})$ predicted by cosmogenic isotope
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Relation between $m/m_0$ obtained from the SINT-2000 RPI stack and $P/P_0$ obtained from (a) $^{10}$Be (this work), (b) $^{10}$Be$_{geo}$ (this work), (c) $^{10}$Be/$^{9}$Be (this work), and (d) $^{10}$Be/$^{9}$Be (from Simon, Thouveny, Bourlès, Valet, et al. (2016)). Solid black lines represent the geomagnetic modulation function $Q$ of $P/P_0$.</figcaption></figure>

production models (solid lines in Figure 7). Taking SINT-2000 or PADM2M as error-free references, the scatter is quantified by the root mean square (RMS) of residuals given by $\Delta = P/P_0 - Q(m/m_0)$. The poor performance of $^{9}$Be as normalizer of the cosmogenic $^{10}$Be record of core MD95-2016 is seen by comparing the residuals obtained from $^{10}$Be/$^{9}$Be (RMS: 0.245 for SINT-2000, Figure 7c, and 0.256 for PADM2M) and from $^{10}$Be (RMS: 0.199 for SINT-2000, Figure 7a, and 0.207 for PADM2M). Consideration of the unnormalized $^{10}$Be record reduces the scatter significantly, but potentially introduces new artifacts caused by temporal changes of the scavenging efficiency. Scavenging artifacts that are reproduced by at least one of the climatic proxies chosen for the analysis of core MD95-2016 are expected to be eliminated by PCA. Indeed, a further reduction of the residuals is obtained with $^{10}$Be$_{geo}$ (RMS: 0.163 for SINT-2000, Figure 7d, and 0.172 for PADM2M). Despite being derived from a single-core record from a site heavily affected by glacial-interglacial cycles, $^{10}$Be$_{geo}$ is compatible, in terms of quality, with a $^{10}$Be/$^{9}$Be stack (RMS: 0.146) chosen for comparison (Simon, Thouveny, Bourlès, Valet, et al., 2016).

Despite the good performance of our PCA analysis, $^{10}$Be$_{geo}$ is not completely free of artifacts, as seen by the temporal patterns of $\Delta$ (Figure 8). We identify three types of artifacts contributing significantly to $\Delta$, with examples highlighted in Figure 8b:

* *Type 1—Spikes.* Spikes tend to be associated with similar short-term variations of L\*, a proxy of the CaCO$_3$ content (Figure 8d), and Ti/K (Figure 8e), a proxy for sediment provenance in the North Atlantic which discriminates basaltic (Ti-rich, large $\chi$) silts transported from local volcanic sources from finer-grained felsic (K-rich, smaller $\chi$) sediments (Mirzaloo et al., 2019).

* *Type 2—Trends.* $\Delta$ contain short- and long-term trends lacking an obvious relation with geomagnetic excursions and with all environmental proxies, except for magnetic susceptibility (Figure 8c). Similar
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Analysis of principal component analysis residuals. (a) $^{10}$Be$_{geo}$, PADM2M, and SINT-2000 converted to $P/P_0$. (b) Difference $\Delta$ between $^{10}$Be$_{geo}$ and SINT-2000 in (a). Examples for three types of artifacts described in the text are highlighted in (a–b). (c–e) Comparison of $\Delta$ (gray, left scale) with magnetic susceptibility $\chi$, lightness L*, and Ti/K (black, right scale). All records on panels c–d were normalized for comparison.</figcaption></figure>

discrepancies are also observed among SINT-2000, PISO-1500, and the $^{10}$Be/$^{9}$Be stack of Simon, Thouveny, Bourlès, Valet, et al. (2016) over the same time intervals (i.e., 150–190, 280–320, and 350–420 ka), thus excluding explanations calling for significant changes of the cosmic ray flux. The origin of these artifacts is discussed later.
* *Type 3—Age model artifacts.* Age model differences between MD95-2016 and the magnetic field reference (in this case, SINT-2000 and PADM2M) contribute to $\Delta$ proportionally to the rate of field change. Accordingly, increased residual amplitudes are expected to coincide with geomagnetic excursions and reversals, as particularly evident for the Icelandic Basin excursion (185–193 ka) highlighted in Figures 8a and 8b. Detailed analysis of this excursion and the Laschamps one (42–43 ka), yield an age offset of +2.7 ka and +3.5 ka for the SINT-2000 records of the Laschamps and the Icelandic Basin excursions, relative to MD95-2016. Correction of these offsets (Figure S6 in Supporting Information S1) eliminates the corresponding artifacts in $\Delta$. The detected age offsets are compatible with the estimated time model uncertainty.

Among the three types of residual artifacts discussed above, Type 3 is the easiest to remove, since it requires a simple adjustment of the age model according to detected age offsets during excursions and reversals. An adjustment of the $^{10}$Be$_{geo}$ record with respect to SINT-2000 reduces the residual RMS from 0.156 to 0.151 (Figures S6 and S7 in Supporting Information S1). Most importantly, it eliminates incorrect interpretations of differences between paleomagnetic and cosmogenic isotope records of the Earth's magnetic field in terms of field geometry effects when the dipole component is not dominant.

Type 1 artifacts are evidently associated with glacial melting. As discussed in Section 5.2, glacial melting causes strong variations in the Be transport pattern. $^{9}$Be is mostly affected, through the ice-rafted delivery of detrital material. Delivery of Be to the MD95-2016 site is controlled by the interplay between the NAC coming from the South, and the Iceland-Scotland Overflow Water (ISOW) and Norwegian Sea Deep Water (NSDW) coming from the North. The lower $^{10}$Be concentration in deep Arctic waters feeding the ISOW and NSDW currents—350–750 at·g$^{-1}$ (Frank et al., 2009), compared to 770–860 at·g$^{-1}$ for mid-latitude deep waters along the NAC (Luo & Ku, 2003)—makes the $^{10}$Be record sensitive to current patterns. Because the source distribution of $^{10}$Be is distinct from that of sediment, mineralogic proxies such as Zr/Sr, Ti/Ca, Si/Sr, and Ti/K, do not fully represent the climatic modulation of $^{10}$Be, which is the reason why artifacts of the $^{10}$Be record are not completely removed by PCA.

Finally, we discuss the observed correlation between $\Delta$ and magnetic susceptibility $\chi$ ($r^2 = 0.63$, Figure 8c), and, to a lesser extent, Ti/K ($r^2 = 0.53$, Figure 8e). Field geometry-related differences between $^{10}$Be$_{geo}$ and RPI records are expected to occur mainly during periods of low dipole field intensity, and this is obviously not the case for $\Delta$, except for Type 3 residuals caused by age model differences. Differences with paleomagnetic records might also reflect a modulation of the $^{10}$Be production rate by cosmic rays. In both cases, however, a correlation with environmental parameters is not expected. We are therefore left with the assumption that either $^{10}$Be$_{geo}$, or the RPI stacks chosen as reference, or both, contain a climatic contamination that is best represented by $\chi$. Coincident spikes (e.g., at 250–260 ka) and stepwise offsets (e.g., 380–400 ka) in $\Delta$ and $\chi$ point to features originating from $^{10}$Be$_{geo}$, since these sudden changes are absent in SINT-2000 and PADM2M (Figures 8a and 8b). Differences between SINT-2000 and PADM2M are generally much smaller than between these records and $^{10}$Be$_{geo}$. $^{10}$Be$_{geo}$ must therefore contains residual artifacts that have not been properly removed by our PCA analysis. One of the largest discrepancies between PADM2M and SINT-2000 occurs at 390–420 ka and coincides with a large Type 2 artifact of the $^{10}$Be$_{geo}$ record (Figures 8a and 8c).

We speculate that the reason for $\chi$ being the environmental parameter that most closely reproduces the residual artifacts of $^{10}$Be$_{geo}$ originates from variations of the scavenging efficiency. The partition coefficient $K_d$ for Be, which is a measure for the scavenging efficiency, is controlled by the particle size and mineralogy of sediment. For instance, calcium carbonate has a much lower affinity toward dissolved Be ($K_d \approx 2 \times 10^4$ ml/g) than biogenic silica ($K_d \approx 1 \times 10^6$ ml/g) and lithogenic fractions ($K_d \approx 2 \times 10^6$ ml/g) (Li, 2005; Luo & Ku, 2004). Since magnetic susceptibility is mainly controlled by the concentration magnetic iron oxides such as magnetite, which share similar $K_d$ values with other detrital components such as TiO$_2$ and silicates (Katsuta et al., 2017), we expect variations of $\chi$ to be representative for the lithogenic fraction, which in turn modulates the scavenging efficiency of the bulk sediment. The fact that $\Delta$ is best represented by $\chi$, rather than $^{9}$Be, might be explained by the different provenance of sediment and dissolved terrestrial Be.

The persistence of an overprint proportional to $\chi$, despite this parameter being one of the environmental proxies included in the set of input records, represents an evident limitation of PCA. Changing the weight of $\chi$ in the PCA input does not modify the results for $^{10}$Be$_{geo}$ (Figure S9 in Supporting Information S1). The unremoved overprint
must therefore be attributed to a violation of the orthogonality condition discussed in Section 4.2, for instance through the preferred occurrence of geomagnetic dipole lows during interglacial or transitional paleoclimatic conditions (Thouveny et al., 2008; Zhou et al., 2023), which was shown by independent observations based on RPI, $^{10}$Be, and records that cannot be affected by climatic overprints, such as seafloor magnetization.

## 5.5. Final Dipole Moment Reconstruction

As discussed in Section 5.4, the PCA-derived geomagnetic component of the $^{10}$Be record contains residual climatic contaminations which correlate well with magnetic susceptibility, as seen by the by the agreement between $\Delta$ and its best-fit analog $\Delta^* = -0.42 + 0.02\chi$ obtained from a linear transformation of the susceptibility signal (Figure 8c). Accordingly, an improved reconstruction of the cosmogenic Be production rate is given by $P*/P_0 = (P/P_0)(10\text{Be}_{\text{geo}}) - \Delta^*$. Even if the coefficients used to calculate $\Delta^*$ are based on the comparison with a reference geomagnetic field model (in this case, SINT-2000), the new model (Figure 9) does not use more external information than the standard processing of cosmogenic Be records, since linear rescaling is always required to convert these records into variations of the geomagnetic dipole moment (e.g., Simon, Thouveny, Bourlès, Valet, et al., 2016).

As expected, discrepancies between the final model and RPI stacks, expressed by the RMS of $P*/P_0 - Q(m/m_0)$, are further reduced in comparison to the PCA result of Figure 8 (0.121 instead of 0.163 for SINT-2000 and 0.126 instead of 0.172 for PADM2M), becoming smaller than those of the $^{10}$Be/$^{9}$Be stack (RMS: 0.146) chosen for comparison. This improvement is mainly due to a general reduction of scatter and the elimination of a systematic offset between ~390 and 460 ka (Type 1 and Type 2 artifacts, respectively, see Section 5.4). Two spikes (asterisks in Figure 9c) introduced by $\Delta^*$ are evident artifacts related to $\chi$ and have been excluded from further processing. Finally, the model of Figure 9c has been converted to a geomagnetic dipole moment (Figure 9d) using $m = m_0 Q^{-1}(P*/P_0)$, where $Q^{-1}$ is the inverse function of $Q$. Our model and the $^{10}$Be/$^{9}$Be stack of Simon, Thouveny, Bourlès, Valet, et al. (2016) feature significant discrepancies relative to each other and relative to the RPI stacks. While differences between cosmogenic and paleomagnetic records can be expected from an unsteady cosmic ray flux over geological times (e.g., Lal et al., 2005), no difference is expected among cosmogenic records because of their global nature and the fact that the cosmic ray flux can be considered constant during the <3 ka mixing time of Be in the oceans (von Blanckenburg & Igel, 1999). In fact, differences between the two cosmogenic records and the RPI stacks are completely unrelated ($r^2 = 0.016$), meaning that they are dominated by recording artifacts, rather than field geometry or the cosmic ray flux.

Overall, the best agreement between the three geomagnetic dipole moment reconstructions of Figure 9d is observed over the first ~120 ka and for the four clearly identifiable excursions at 42–43 ka (Laschamps), 62–64 ka (Norwegian Greenland Sea), 120 ka (Blake), and 185–193 ka (Icelandic Basin). At older ages, the $^{10}$Be/$^{9}$Be stack is affected by large oscillations that are not reflected by the other records. An overall good agreement between our model and RPI stacks is also observed at ~200–260 ka and 540–570 ka, including the Big Lost excursion. Assuming the remaining differences with RPI stacks to be entirely caused by the cosmogenic record, a further reduction of residual overprints by a factor ~4 (from $\delta P/P_0 \approx 0.2$, e.g. at 160–180 ka, to 0.05, which is the amplitude of short-time wiggles in SINT-2000 converted to $P/P_0$) would be required to resolve true differences with SINT-2000 or PADM2M. This result might be attained by stacking ~16 cosmogenic Be records carrying independent climatic overprints of similar amplitude as in MD95-2016. Less records would be needed in case of conditions more favorable than those of North Atlantic sites.

## 6. Conclusions and Outlook

Sedimentary records of the Earth's magnetic field are unavoidably affected by environmental overprints. Overprints with strong regional or global correlation, such as those driven by the global climate, are difficult or impossible to remove by conventional stacking. Furthermore, high-resolution cosmogenic records, for example, of $^{10}$Be, are limited in number, so that other techniques are required to disentangle geomagnetic and climatic signals in such records. PCA has been used in the past to separate different common contributions within a set of environmental or paleomagnetic records; however, the success of this approach depends heavily on the number of available records (usually $\gg 10$), and on how well critical assumptions for the interpretation of principal components are fulfilled. For instance, the multiplicative climatic modulation of geomagnetic signals formally violates the linearity requirements of PCA. Even after overcoming linearity issues, the limited number of
<figure><img src="image_11.png" /></figure>

**Figure 9.** Final model for the cosmogenic $^{10}$Be record of MD95-2016. (a) Normalized geomagnetic dipole moment from SINT-2000 versus normalized cosmogenic $^{10}$Be production rate for $^{10}$Be$_{geo}$ − Δ$^{*}$ (dots). The geomagnetic modulation function *Q* is shown for comparison (solid line). (b) Same as (a) for the $^{10}$Be/$^{9}$Be stack of Simon, Thouveny, Bourlès, Valet, et al. (2016). (c) Comparison between normalized cosmogenic $^{10}$Be production rates obtained from SINT-2000, PADM2M, $^{10}$Be$_{geo}$, and $^{10}$Be$_{geo}$ − Δ$^{*}$. Asterisks highlight two spikes produced by Δ$^{*}$, which are evidently an artifact related to $\chi$. (d) Comparison between normalized geomagnetic dipole moments obtained from SINT-2000, the $^{10}$Be/$^{9}$Be stack of Simon, Thouveny, Bourlès, Valet, et al. (2016), and $^{10}$Be$_{geo}$ − Δ$^{*}$. Artifacts marked by asterisks in (c) have been removed.

geomagnetic events during the last millions of years and the scarcity of $^{10}$Be records poses statistical limits that
adds uncertainties to the PCA results.

In this work, we introduced a new approach for dealing with these limitations, which is based on the integration of
geomagnetic records with a set of relevant environmental proxies expected to represent the climatic modulation to
be removed. Inexpensive environmental proxies can be obtained for instance from elemental XRF measurements,
colorimetric analyses, and magnetic measurements. As a proof of concept, we applied this new approach to the
worst-case scenario represented by a single authigenic $^{10}$Be record from the North Atlantic core MD95-2016
covering the last 600 ka. In this core, $^{10}$Be and, to a greater extent, $^{10}$Be/$^{9}$Be, are dominated by glacial-
interglacial cycles. The benefit of removing lithology-driven variations of the scavenging efficiency through normalization with ⁹Be is obscured by even stronger variations of ⁹Be delivery, for example, by ice-rafting debris.

Auxiliary measurements of δ¹⁸O, magnetic susceptibility χ, sediment lightness L*, and elemental ratios Si/Sr, Ca/Sr, Sr/Rb, Ti/Ca, Ti/K, and Zr/Sr from the same core have been therefore used as climatic proxies expected to represent—individually or as linear combination—the different climatic overprints of ¹⁰Be and ¹⁰Be/⁹Be. PCA analysis of the complete set of records yielded 10 principal components: the first two components, which explain 64.4% of the total signal variance, could be identified with a global climatic signal and a geomagnetic signal, respectively, by comparison with δ¹⁸O and with global RPI stacks. The third principal component reflects regional environmental overprints as revealed by comparison with SST variations obtained from a neighbor site. The remaining principal components do not display obvious relations with climatic proxies, nor with the RPI stack. Accordingly, ¹⁰Be variations related to the geomagnetic field have been identified with the contribution ¹⁰Begeo of the corresponding principal component to the ¹⁰Be record.

The geomagnetic component ¹⁰Begeo yields a reconstruction of the normalized cosmogenic ¹⁰Be production rate *P*, which can be compared with an independent reconstruction of *P* obtained from paleomagnetic records. The geomagnetic modulation function *Q* derived from cosmogenic isotope production models (e.g., Poluianov et al., 2016) serves as key for transforming *P* into the geomagnetic dipole moment *m* and vice-versa, through proper normalization of the corresponding records with modern-time values *P*₀ and *m*₀, respectively. Comparisons have been performed with three RPI stacks: two global ones (SINT-2000 and PADM2M), which includes records from the Atlantic, Indian, and Pacific Oceans, and a regional one, which is based mainly on sediment cores from the North Atlantic. Contrary to what can be expected from the common geographical provenance of our site and those underlying PISO-1500, ¹⁰Begeo agrees better with SINT-2000 and PADM2M (Figure 9c), testifying the presence of significant climatic overprints in PISO-1500, and the reduction or elimination of these overprints in the ¹⁰Be record. The effectivity of PCA in removing climatic overprints from the ¹⁰Be record is also demonstrated by the reduction of scatter in plots of *m* versus *P* (Figure 7) obtained from ¹⁰Be/⁹Be (0.245), ¹⁰Be (0.204), and ¹⁰Begeo (0.156).

Residual differences Δ between ¹⁰Begeo and SINT-2000 correlate surprisingly well with χ (Figure 8c). Careful analysis of relevant temporal features of Δ excludes a climatic contamination of SINT-2000 or PADM2M as possible source, attributing the origin of this signal to an unremoved ¹⁰Be overprint caused by variations of the beryllium scavenging efficiency. The persistence of this overprint in the geomagnetic PC, despite χ being one of the PCA input records, is attributed to a violation of the orthogonality condition by which the signals to be separated—in this case the geomagnetic modulation of ¹⁰Be and the climatic overprints—must be uncorrelated. The preferred occurrence of geomagnetic dipole lows during interglacial or transitional paleoclimatic conditions (Thouveny et al., 2008; Zhou et al., 2023) is most likely at the origin of this problem. After proving the origin of Δ and its relation with lithological variations, it was possible to reduce the residual climatic overprint by decomposing Δ into a component proportional to χ and a residual signal, and subtracting the former from ¹⁰Begeo. The final reconstruction of *P* yield a further reduction of the RMS difference with SINT-2000, from 0.156 to 0.121. The final result (Figure 9d) contains significantly less artifacts than a ¹⁰Be/⁹Be stack obtained from two sites in the Pacific Ocean with little glacial-interglacial variability. This demonstrates the effectivity of the proposed approach, which reduced climatic overprints of a single cosmogenic record by a factor of ~2. This improvement is similar to what can be expected from the simple stack of four independent records that do not share a common climatic modulation.

Our analysis of core MD95-2016 provides important insights into the mechanisms that control ¹⁰Be records in the North Atlantic, which can be summarized as follows:

* ¹⁰Be and ⁹Be display a long-term, quasi-linear increase in concentration with time. The absence of a similar trend in ¹⁰Be/⁹Be and the lack of a plausible mechanism that would explain a steady concentration increase of both isotopes over the last 600 ka points to a diagenetic release of authigenic Be, as postulated in previous studies.

* The lack of a constant term in the ¹⁰Be record of this core (corresponding to the coefficient A in Equation 1) indicates that, contrary to sites with a much larger water depth, diagenetically released Be does not accumulate at the ocean floor.
* ⁹Be is much more strongly modulated by glacial-interglacial cycles than ¹⁰Be. This modulation is also reflected by ¹⁰Be/⁹Be, which, for this site, is a poorer geomagnetic field recorder than unnormalized ¹⁰Be. This means that isotope-specific differences in the climatic modulation of dissolved Be play a more important role than variations of the scavenging efficiency driven by the sediment lithology.

* The main climatic modulation of the ¹⁰Be record identified by PCA matches the global δ¹⁸O stack. The second most important overprint matches regional environmental changes, as shown by comparisons with the sea surface temperature record derived from a neighbor site.

* The unnormalized ¹⁰Be record contains artifacts caused by temporal changes of the scavenging efficiency. While variations in sediment grain size are controlled mainly by glacial-interglacial cycles and captured by the first principal component, additional changes of the scavenging efficiency caused by Be-adsorbing minerals associated with detrital iron oxides are well reflected by magnetic susceptibility.

Our approach to the analysis of the cosmogenic ¹⁰Be record of a single site can be generalized and extended to analyze simultaneously several sites and different types of records, using the following processing steps:

1. For each site, collect, where possible, paleomagnetic and cosmogenic records of the Earth's magnetic field, along with a set of environmental proxies (e.g., δ¹⁸O, XRF elemental ratios, rock-magnetic parameters).

2. For each site, use PCA to obtain the geomagnetic component(s) from the original measurements.

3. For each site, compare the obtained geomagnetic components with reference records of the Earth's magnetic field in order to identify possible residual climatic overprints. In case of a good correlation between residual overprints and one of the environmental proxies, subtract the contribution of this proxy to obtain final, site-specific geomagnetic records.

4. Use PCA to identify the common geomagnetic signal among the site-specific records obtained previously.

This approach is expected to be much more efficient than traditional stacking.

## Data Availability Statement

The ¹⁰Be and ⁹Be (Savranskaia et al., 2024a), δ¹⁸O (Savranskaia et al., 2024b), and XRF (Savranskaia et al., 2024c) data supporting this study is openly available in the Zenodo repository.

## References

Ahn, S., Khider, D., Lisiecki, L. E., & Lawrence, C. E. (2017). A probabilistic Pliocene-Pleistocene stack of benthic δ¹⁸O using a profile hidden Markov model. *Dynamics and Statistics of the Climate System*, 2(1), dxz002. https://doi.org/10.1093/climsys/dzx002

Babamoradi, H., van den Berg, F., & Rinnan, A. (2013). Bootstrap based confidence limits in principal component analysis – A case study. *Chemometrics and Intelligent Laboratory Systems*, 120, 97–105. https://doi.org/10.1016/j.chemolab.2012.10.007

Ballini, M., Kissel, C., Colin, C., & Richter, T. (2006). Deep-water mass source and dynamic associated with rapid climatic variations during the last glacial stage in the North Atlantic: A multiproxy investigation of the detrital fraction of deep-sea sediments. *Geochemistry, Geophysics, Geosystems*, 7(2), Q02N01. https://doi.org/10.1029/2005gc001070

Balsam, W. L., & McCoy, F. W., Jr. (1987). Atlantic sediments: Glacial/Interglacial comparisons. *Paleoceanography*, 2(5), 531–542. https://doi.org/10.1029/pa002i005p00531

Barker, S., Broecker, W., Clark, E., & Hajdas, I. (2007). Radiocarbon age offsets of foraminifera resulting from differential dissolution and fragmentation within the sedimentary bioturbated zone. *Paleoceanography*, 22(2), PA2205. https://doi.org/10.1029/2006pa001354

Barker, S., Chen, J., Gong, X., Jonkers, L., Knorr, G., & Thornalley, D. (2015). Icebergs not the trigger for North Atlantic cold events. *Nature*, 520(7547), 333–336. https://doi.org/10.1038/nature14330

Beer, J., McCracken, K., & Steiger, R. (2012). *Cosmogenic radionuclides: Theory and applications in the terrestrial and space environments*. Springer Science & Business Media.

Blinov, A. (1988). The dependence of cosmogenic isotope production rate on solar activity and geomagnetic field variations. In *Secular Solar and Geomagnetic Variations in the Last 10,000 Years* (pp. 329–340). Springer Netherlands. https://doi.org/10.1007/978-94-009-3011-7_20

Bond, G., Heinrich, H., Broecker, W., Labeyrie, L., McManus, J., Andrews, J., et al. (1992). Evidence for massive discharges of icebergs into the North Atlantic ocean during the last glacial period. *Nature*, 360(6401), 245–249. https://doi.org/10.1038/360245a0

Boschi, V., & Willenbring, J. K. (2021). Chemical and physical drivers of beryllium retention in two soil endmembers. *Science of the Total Environment*, 754, 141591. https://doi.org/10.1016/j.scitotenv.2020.141591

Bourlès, D., Raisbeck, G., & Yiou, F. (1989). ¹⁰Be and ⁹Be in marine sediments and their potential for dating. *Geochimica et Cosmochimica Acta*, 53(2), 443–452. https://doi.org/10.1016/0016-7037(89)90395-5

Burton, K. W., Ling, H.-F., & O’Nions, R. K. (1997). Closure of the Central American Isthmus and its effect on deep-water formation in the North Atlantic. *Nature*, 386(6623), 382–385. https://doi.org/10.1038/386382a0

Calvert, S., & Pedersen, T. (2007). Chapter fourteen elemental proxies for palaeoclimatic and palaeoceanographic variability in marine sediments: Interpretation and application. *Developments in Marine Geology*, 1, 567–644. https://doi.org/10.1016/s1572-5480(07)01019-6

Carcaillet, J., Bourlès, D. L., Thouveny, N., & Arnold, M. (2004). A high resolution authigenic ¹⁰Be/⁹Be record of geomagnetic moment variations over the last 300 ka from sedimentary cores of the Portuguese margin. *Earth and Planetary Science Letters*, 219(3–4), 397–412. https://doi.org/10.1016/s0012-821x(03)00702-7

## Acknowledgments

This study was supported by Alexander von Humboldt Research Fellowships for Postdoctoral Researchers to Tatiana Savranskaia and the ERC advanced grant to Jean-Pierre Valet (1953–2024). Grant agreement ID: 339899–“EDIFICE”, funded under the ERC’s 7th Framework Program FP7-IDEAS-ERC (European Union). The ASTER AMS national facility (CEREGE, Aix en Provence, France) is supported by INSU/CNRS, IRD and by the ANR through the EQUIPEX “ASTER-CEREGE” action. Didier Bourlès, who passed away in 2021, directed measurements of beryllium isotopes and contributed important insights to the discussion of beryllium geochemistry. We also would like to express our gratitude to Naoufel Haddam and Isabelle Billy for conducting the XRF measurements, and to Anojh Thevarasan for magnetic measurements at IPGP. We acknowledge the Associated Editor, Lisa Tauxe, and an anonymous reviewer for their constructive comments, which helped improving the manuscript. Open Access funding enabled and organized by Projekt DEAL.
Channell, J. E. T., Hodell, D. A., McManus, J., & Lehman, B. (1998). Orbital modulation of the Earth's magnetic field intensity. *Nature*, **394**(6692), 464–468. https://doi.org/10.1038/28833

Channell, J. E. T., Xuan, C., & Hodell, D. A. (2009). Stacking paleointensity and oxygen isotope data for the last 1.5 myr (PISO-1500). *Earth and Planetary Science Letters*, **283**(1–4), 14–23. https://doi.org/10.1016/j.epsl.2009.03.012

Chase, Z., Anderson, R. F., Fleisher, M. Q., & Kubik, P. W. (2002). The influence of particle composition and particle flux on scavenging of Th, Pa and Be in the ocean. *Earth and Planetary Science Letters*, **204**(1–2), 215–229. https://doi.org/10.1016/s0012-821x(02)00984-6

Chmeleff, J., von Blanckenburg, F., Kossert, K., & Jakob, D. (2010). Determination of the ¹⁰Be half-life by multicollector ICP-MS and liquid scintillation counting. *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms*, **268**(2), 192–199. https://doi.org/10.1016/j.nimb.2009.09.012

Christl, M., Lippold, J., Steinhilber, F., Bernsdorff, F., & Mangini, A. (2010). Reconstruction of global ¹⁰Be production over the past 250 ka from highly accumulating Atlantic drift sediments. *Quaternary Science Reviews*, **29**(19–20), 2663–2672. https://doi.org/10.1016/j.quascirev.2010.06.017

Courtillot, V., Gallet, Y., Le Mouël, J.-L., Fluteau, F., & Genevey, A. (2007). Are there connections between the Earth's magnetic field and climate? *Earth and Planetary Science Letters*, **253**(3–4), 328–339. https://doi.org/10.1016/j.epsl.2006.10.032

Croudace, I. W., & Rothwell, R. G. (2015). *Micro-XRF studies of sediment cores: Applications of a non-destructive tool for the environmental sciences* (Vol. 17). Springer.

Deng, K., Rickli, J., Suhrhoff, T. J., Du, J., Scholz, F., Severmann, S., et al. (2023). Dominance of benthic fluxes in the oceanic Beryllium budget and implications for paleo-denudation records. *Science Advances*, **9**(23), eadg3702. https://doi.org/10.1126/sciadv.adg3702

Egli, R., & Zhao, X. (2015). Natural remanent magnetization acquisition in bioturbated sediment: General theory and implications for relative paleointensity reconstructions. *Geochemistry, Geophysics, Geosystems*, **16**(4), 995–1016. https://doi.org/10.1002/2014gc005672

Eisenhauer, A., Spielhagen, R., Frank, M., Hentzschel, G., Mangini, A., Kubik, P., et al. (1994). 10Be records of sediment cores from high northern latitudes: Implications for environmental and climatic changes. *Earth and Planetary Science Letters*, **124**(1–4), 171–184. https://doi.org/10.1016/0012-821x(94)00069-7

Faugères, J.-C., Stow, D. A., Imbert, P., & Viana, A. (1999). Seismic features diagnostic of contourite drifts. *Marine Geology*, **162**(1), 1–38. https://doi.org/10.1016/s0025-3227(99)00068-7

Foster, G. L., Vance, D., & Prytulak, J. (2007). No change in the neodymium isotope composition of deep water exported from the North Atlantic on glacial-interglacial time scales. *Geology*, **35**(1), 37–40. https://doi.org/10.1130/g23204a.1

Frank, M. (2002). Radiogenic isotopes: Tracers of past ocean circulation and erosional input. *Reviews of Geophysics*, **40**(1), 1–1–1–38. https://doi.org/10.1029/2000rg000094

Frank, M., Backman, J., Jakobsson, M., Moran, K., O'Regan, M., King, J., et al. (2008). Beryllium isotopes in central Arctic Ocean sediments over the past 12.3 million years: Stratigraphic and paleoclimatic implications. *Paleoceanography*, **23**(1), PA1S02. https://doi.org/10.1029/2007pa001478

Frank, M., Porcelli, D., Andersson, P., Baskaran, M., Björk, G., Kubik, P. W., et al. (2009). The dissolved Beryllium isotope composition of the Arctic Ocean. *Geochimica et Cosmochimica Acta*, **73**(20), 6114–6133. https://doi.org/10.1016/j.gca.2009.07.010

Franke, C., Hofmann, D., & von Dobeneck, T. (2004). Does lithology influence relative paleointensity records? A statistical analysis on South Atlantic pelagic sediments. *Physics of the Earth and Planetary Interiors*, **147**(2–3), 285–296. https://doi.org/10.1016/j.pepi.2004.07.004

Frijia, G., & Parente, M. (2008). Strontium isotope stratigraphy in the upper Cenomanian shallow-water carbonates of the southern Apennines: Short-term perturbations of marine ⁸⁷Sr/⁸⁶Sr during the oceanic anoxic event 2. *Palaeogeography, Palaeoclimatology, Palaeoecology*, **261**(1–2), 15–29. https://doi.org/10.1016/j.palaeo.2008.01.003

Görtler, J., Spinner, T., Streeb, D., Weiskopf, D., & Deussen, O. (2020). Uncertainty-aware principal component analysis. *IEEE Transactions on Visualization and Computer Graphics*, **26**(1), 822–831. https://doi.org/10.1109/tvcg.2019.2934812

Heikkilä, U., Beer, J., & Feichter, J. (2008). Modeling cosmogenic radionuclides ¹⁰Be and ⁷Be during the Maunder Minimum using the ECHAM5-HAM general circulation model. *Atmospheric Chemistry and Physics*, **8**(10), 2797–2809. https://doi.org/10.5194/acp-8-2797-2008

Herbert, T. D., Lawrence, K. T., Tzanova, A., Peterson, L. C., Caballero-Gill, R., & Kelly, C. S. (2016). Late Miocene global cooling and the rise of modern ecosystems. *Nature Geoscience*, **9**(11), 843–847. https://doi.org/10.1038/ngeo2813

Hodell, D. A., Channell, J. E., Curtis, J. H., Romero, O. E., & Röhl, U. (2008). Onset of “Hudson Strait” Heinrich events in the eastern North Atlantic at the end of the middle Pleistocene transition (~640 ka)? *Paleoceanography*, **23**(4), PA4218. https://doi.org/10.1029/2008pa001591

Hodell, D. A., Evans, H. F., Channell, J. E., & Curtis, J. H. (2010). Phase relationships of North Atlantic ice-rafted debris and surface-deep climate proxies during the last glacial period. *Quaternary Science Reviews*, **29**(27–28), 3875–3886. https://doi.org/10.1016/j.quascirev.2010.09.006

Hofmann, D. I., & Fabian, K. (2009). Correcting relative paleointensity records for variations in sediment composition: Results from a South Atlantic stratigraphic network. *Earth and Planetary Science Letters*, **284**(1–2), 34–43. https://doi.org/10.1016/j.epsl.2009.03.043

Jansen, E., Raymo, M. E., & Blum, P. (1996). 1. Leg 162: New frontiers on past climates. In *Proceedings ODP, Initial Reports* (Vol. 162, pp. 5–20).

Jolliffe, I. T. (2002). *Principal component analysis*. Springer.

Katsuta, S., Kanaya, N., Bessho, K., & Monjushiro, H. (2017). Adsorption behavior of trace Beryllium (II) onto metal oxide nanoparticles dispersed in water. *International Journal of Chemistry*, **9**(4), 62–70. https://doi.org/10.5539/ijc.v9n4p62

Kim, J., Goldstein, S. L., Pena, L. D., Jaume-Seguí, M., Knudson, K. P., Yehudai, M., & Bolge, L. (2021). North Atlantic deep water during Pleistocene interglacials and glacials. *Quaternary Science Reviews*, **269**, 107146. https://doi.org/10.1016/j.quascirev.2021.107146

King, J., Banerjee, S. K., Marvin, J., & Özdemir, Ö. (1982). A comparison of different magnetic methods for determining the relative grain size of magnetite in natural materials: Some results from lake sediments. *Earth and Planetary Science Letters*, **59**(2), 404–419. https://doi.org/10.1016/0012-821x(82)90142-x

Kissel, C. (2005). Magnetic signature of rapid climatic variations in glacial North Atlantic, a review. *Comptes Rendus Geoscience*, **337**(10–11), 908–918. https://doi.org/10.1016/j.crte.2005.04.009

Kissel, C., Laj, C., Labeyrie, L., Dokken, T., Voelker, A., & Blamart, D. (1999). Rapid climatic variations during marine isotopic stage 3: Magnetic analysis of sediments from Nordic Seas and North Atlantic. *Earth and Planetary Science Letters*, **171**(3), 489–502. https://doi.org/10.1016/s0012-821x(99)00162-4

Kissel, C., Laj, C., Mulder, T., Wandres, C., & Cremer, M. (2009). The magnetic fraction: A tracer of deep water circulation in the North Atlantic. *Earth and Planetary Science Letters*, **288**(3–4), 444–454. https://doi.org/10.1016/j.epsl.2009.10.005

Knudsen, M. F., Henderson, G. M., Frank, M., Mac Niocaill, C., & Kubik, P. W. (2008). In-phase anomalies in Beryllium-10 production and palaeomagnetic field behaviour during the Iceland Basin geomagnetic excursion. *Earth and Planetary Science Letters*, **265**(3–4), 588–599. https://doi.org/10.1016/j.epsl.2007.10.051
Kocharov, G., Blinov, A., Konstantinov, A., & Levchenko, V. (1989). Temporal $^{10}$Be and $^{14}$C variations: A tool for paleomagnetic research. *Radiocarbon, 31*(2), 163-168. https://doi.org/10.1017/s0033822200044829

Korschinek, G., Bergmaier, A., Faestermann, T., Gerstmann, U., Knie, K., Rugel, G., et al. (2010). A new value for the half-life of $^{10}$Be by heavy-ion elastic recoil detection and liquid scintillation counting. *Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 268*(2), 187-191. https://doi.org/10.1016/j.nimb.2009.09.020

Korte, M., Constable, C., Donadini, F., & Holme, R. (2011). Reconstructing the holocene geomagnetic field. *Earth and Planetary Science Letters, 312*(3-4), 497-505. https://doi.org/10.1016/j.epsl.2011.10.031

Kylander, M. E., Ampel, L., Wohlfarth, B., & Veres, D. (2011). High-resolution X-ray fluorescence core scanning analysis of Les Echets (France) sedimentary sequence: New insights from chemical proxies. *Journal of Quaternary Science, 26*(1), 109-117. https://doi.org/10.1002/jqs.1438

Labeyrie, L., Waelbroeck, C., Cortijo, E., Michel, E., & Duplessy, J.-C. (2005). Changes in deep water hydrology during the Last Deglaciation. *Comptes Rendus Geoscience, 337*(10-11), 919-927. https://doi.org/10.1016/j.crte.2005.05.010

Labeyrie, L. D. (1996). IMAGES-MD 101 à bord du Marion-Dufresne du 29 Mai au 11 Juillet 1995. A coring cruise of the R/V Marion Dufresne in the North Atlantic Ocean and Norwegian Sea. *Les Publications de l'Institut Français Pour la Recherche et la Technologie Polaires. Les Rapports des Campagnes à la Mer, 96*(1).

Laj, C., & Channell, J. (2015). 5.10 - Geomagnetic Excursions. In G. Schubert (Ed.), *Treatise on geophysics* (2nd ed., pp. 343-383). Elsevier. https://doi.org/10.1016/B978-0-444-53802-4-00104-4

Lal, D. (1992). Expected secular variations in the global terrestrial production rate of radiocarbon. In E. Bard & W. S. Broecker (Eds.), *The last deglaciation: Absolute and radiocarbon chronologies* (pp. 114-126). Springer.

Lal, D., Jull, A. J. T., Pollard, D., & Vacher, L. (2005). Evidence for large century time-scale changes in solar activity in the past 32 kyr, based on in-situ cosmogenic $^{14}$C in ice at Summit, Greenland. *Earth and Planetary Science Letters, 234*(3-4), 335-349. https://doi.org/10.1016/j.epsl.2005.02.011

Lal, D., & Peters, B. (1967). *Cosmic ray produced radioactivity on the Earth, Kosmische Strahlung II/Cosmic Rays II* (pp. 551-612). Springer Berlin Heidelberg.

Lawrence, K. T., Herbert, T. D., Brown, C. M., Raymo, M. E., & Haywood, A. M. (2009). High-amplitude variations in North Atlantic sea surface temperature during the early Pliocene warm period. *Paleoceanography, 24*(2), PA2218. https://doi.org/10.1029/2008pa001669

Lee, T., Rand, D., Lisiecki, L. E., Gebbie, G., & Lawrence, C. E. (2022). Bayesian age models and stacks: Combining age inferences from radiocarbon and benthic $^{δ^{18}O}$stratigraphic alignment. *EGUsphere, 2022*, 1-29.

Li, Y.-H. (2005). Controversy over the relationship between major components of sediment-trap materials and the bulk distribution coefficients of $^{230}$Th, $^{231}$Pa, and $^{10}$Be. *Earth and Planetary Science Letters, 233*(1-2), 1-7. https://doi.org/10.1016/j.epsl.2005.02.023

Lisiecki, L. E., & Raymo, M. E. (2005). A Pliocene-Pleistocene stack of 57 globally distributed benthic $^{δ^{18}O}$records. *Paleoceanography, 20*(1), PA1003. https://doi.org/10.1029/2004pa001071

Lisiecki, L. E., & Stern, J. V. (2016). Regional and global benthic $^{δ^{18}O}$stacks for the last glacial cycle. *Paleoceanography, 31*(10), 1368-1394. https://doi.org/10.1002/2016pa003002

Luo, S., & Ku, T. L. (2003). Constraints on deep-water formation from the oceanic distribution of $^{10}$Be. *Journal of Geophysical Research, 108*(C5), 3137. https://doi.org/10.1029/2002jc001670

Luo, S., & Ku, T.-L. (2004). On the importance of opal, carbonate, and lithogenic clays in scavenging and fractionating $^{230}$Th, $^{231}$Pa and $^{10}$Be in the ocean. *Earth and Planetary Science Letters, 220*(1-2), 201-211. https://doi.org/10.1016/s0012-821x(04)00027-5

Margari, V., Skinner, L. C., Menviel, L., Capron, E., Rhodes, R. H., Mleneck-Vautravers, M. J., et al. (2020). Fast and slow components of interstadial warming in the North Atlantic during the last glacial. *Communications Earth & Environment, 1*(1), 1-9. https://doi.org/10.1038/s43247-020-0006-x

Masarik, J., & Beer, J. (2009). An updated simulation of particle fluxes and cosmogenic nuclide production in the Earth's atmosphere. *Journal of Geophysical Research, 114*(D11), D11103. https://doi.org/10.1029/2008jd010557

McHargue, L., & Donahue, D. (2005). Effects of climate and the cosmic-ray flux on the $^{10}$Be content of marine sediments. *Earth and Planetary Science Letters, 232*(1-2), 193-207. https://doi.org/10.1016/j.epsl.2004.11.018

McManus, J. F., Oppo, D. W., & Cullen, J. L. (1999). A 0.5-million-year record of millennial-scale climate variability in the North Atlantic. *Science, 283*(5404), 971-975. https://doi.org/10.1126/science.283.5404.971

Measures, C., & Edmond, J. (1983). The geochemical cycle of $^{9}$Be: A reconnaissance. *Earth and Planetary Science Letters, 66*, 101-110. https://doi.org/10.1016/0012-821x(83)90129-2

Ménabréaz, L., Bourlès, D. L., & Thouveny, N. (2012). Amplitude and timing of the Laschamp geomagnetic dipole low from the global atmospheric $^{10}$Be overproduction: Contribution of authigenic $^{10}$Be/$^{9}$Be ratios in west equatorial Pacific sediments. *Journal of Geophysical Research, 117*(B11), B11101. https://doi.org/10.1029/2012jb009256

Mirzaloo, M., Nürnberg, D., Kienast, M., & van der Lubbe, H. J. L. (2019). Synchronous changes in sediment transport and provenance at the Icelad-Faroe Ridge linked to millennial climate variability from 55 to 6 ka BP. *Geochemistry, Geophysics, Geosystems, 20*(8), 4184-4201. https://doi.org/10.1029/2019gc008298

Moreno, E., Thouveny, N., Delanghe, D., McCave, I. N., & Shackleton, N. J. (2002). Climatic and oceanographic changes in the Northeast Atlantic reflected by magnetic properties of sediments deposited on the Portuguese Margin during the last 340 ka. *Earth and Planetary Science Letters, 202*(2), 465-480. https://doi.org/10.1016/s0012-821x(02)00787-2

Müller-Michaelis, A., & Uenzelmann-Neben, G. (2014). Development of the Western Boundary Undercurrent at Eirik Drift related to changing climate since the early Miocene. *Deep Sea Research Part I: Oceanographic Research Papers, 93*, 21-34. https://doi.org/10.1016/j.dsr.2014.07.010

Nilsson, A., Holme, R., Korte, M., Suttie, N., & Hill, M. (2014). Reconstructing Holocene geomagnetic field variation: New methods, models and implications. *Geophysical Journal International, 198*(1), 229-248. https://doi.org/10.1093/gji/ggu120

Oyedotun, T. D. T. (2018). X-ray fluorescence (XRF) in the investigation of the composition of Earth materials: A review and an overview. *Geology, Ecology, and Landscapes, 2*(2), 148-154. https://doi.org/10.1080/24749508.2018.1452459

Petit, T., Mercier, H., & Thierry, V. (2018). Firavranskaiat direct estimates of volume and water mass transports across the Reykjanes Ridge. *Journal of Geophysical Research: Oceans, 123*(9), 6703-6719. https://doi.org/10.1029/2018jc013999

Poluianov, S. V., Kovaltsov, G. A., Mishev, A. L., & Usokin, I. G. (2016). Production of cosmogenic isotopes $^{7}$Be, $^{10}$Be, $^{14}$C, $^{22}$Na, and $^{36}$Cl in the atmosphere: Altitudinal profiles of yield functions. *Journal of Geophysical Research: Atmospheres, 121*(13), 8125-8136. https://doi.org/10.1002/2016jd025034

Raymo, M. E. (1994). The initiation of Northern Hemisphere glaciation. *Annual Review of Earth and Planetary Sciences, 22*(1), 353-383. https://doi.org/10.1146/annurev-earth.22.1.353
Richter, T. O., Van der Gaast, S., Koster, B., Vaars, A., Gieles, R., de Stigter, H. C., et al. (2006). The Avaatech XRF core scanner: Technical description and applications to NE Atlantic sediments. *Geological Society, London, Special Publications*, **267**(1), 39–50. https://doi.org/10.1144/gsl.sp.2006.267.01.03

Roberts, A. P. (2008). Geomagnetic excursions: Knowns and unknowns. *Geophysical Research Letters*, **35**(17), L17307. https://doi.org/10.1029/2008gl034719

Roberts, A. P., Tauxe, L., & Heslop, D. (2013). Magnetic paleointensity stratigraphy and high-resolution geochronology: Successes and future challenges. *Quaternary Science Reviews*, **61**, 1–43.

Ruddiman, W. F., Tolderlund, D. S., & Bé, A. (1970). Foraminiferal evidence of a modern warming of the North Atlantic Ocean. In *Deep sea research and oceanographic abstracts* (Vol. 17, pp. 141–155). Elsevier. https://doi.org/10.1016/0011-7471(70)90093-8

Savranskaia, T., Egli, R., & Valet, J.-P. (2022). Multiscale Brazil nut effects in bioturbated sediment. *Scientific Reports*, **12**(1), 11450. https://doi.org/10.1038/s41598-022-14276-w

Savranskaia, T., Egli, R., Valet, J.-P., Bassinot, F., Meynadier, L., Bourlès, D. L., et al. (2021). Disentangling magnetic and environmental signatures of sedimentary $^{10}$Be/$^{9}$Be records. *Quaternary Science Reviews*, **257**, 106809. https://doi.org/10.1016/j.quascirev.2021.106809

Savranskaia, T., Egli, R., Valet, J.-P., Simon, Q., Bassinot, F., & Thouveny, N. (2024a). Authigenic $^{10}$Be and $^{9}$Be concentrations and $^{10}$Be/$^{9}$Be ratios of sediment core MD95-2016 [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.14011681

Savranskaia, T., Egli, R., Valet, J.-P., Simon, Q., Bassinot, F., & Thouveny, N. (2024b). Stable isotopic oxygen measurements of sediment core MD95-2016 [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.14011714

Savranskaia, T., Egli, R., Valet, J.-P., Simon, Q., Bassinot, F., & Thouveny, N. (2024c). X-ray fluorescence (XRF) scans of sediment core MD95-2016 [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.14011735

Sharma, P., Mahannah, R., Moore, W., Ku, T., & Southon, J. (1987). Transport of $^{10}$Be and $^{9}$Be in the ocean. *Earth and Planetary Science Letters*, **86**(1), 69–76. https://doi.org/10.1016/0012-821x(87)90189-0

Simon, Q., Thouveny, N., Bourlès, D. L., Bassinot, F., Savranskaia, T., Valet, J.-P., & Team, A. (2018). Increased production of cosmogenic $^{10}$Be recorded in oceanic sediment sequences: Information on the age, duration, and amplitude of the geomagnetic dipole moment minimum over the Matuyama-Brunhes transition. *Earth and Planetary Science Letters*, **489**, 191–202. https://doi.org/10.1016/j.epsl.2018.02.036

Simon, Q., Thouveny, N., Bourles, D. L., Nuttin, L., Hillaire-Marcel, C., & St-Onge, G. (2016). Authigenic $^{10}$Be/$^{9}$Be ratios and $^{10}$Be-fluxes ($^{230}$Th$_{xs}$-normalized) in central Baffin Bay sediments during the last glacial cycle: Paleoenvironmental implications. *Quaternary Science Reviews*, **140**, 142–162. https://doi.org/10.1016/j.quascirev.2016.03.027

Simon, Q., Thouveny, N., Bourlès, D. L., Valet, J.-P., Bassinot, F., Ménabréaz, L., et al. (2016). Authigenic $^{10}$Be/$^{9}$Be ratio signatures of the cosmogenic nuclide production linked to geomagnetic dipole moment variation since the Brunhes/Matuyama boundary. *Journal of Geophysical Research: Solid Earth*, **121**(11), 7716–7741. https://doi.org/10.1002/2016jb013335

Spratt, R. M., & Lisiecki, L. E. (2016). A late pleistocene sea level stack. *Climate of the Past*, **12**(4), 1079–1092. https://doi.org/10.5194/cp-12-1079-2016

Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D., Davies, S. M., et al. (2006). The Greenland ice core chronology 2005, 15–42 ka. Part 2: Comparison to other records. *Quaternary Science Reviews*, **25**(23–24), 3258–3267. https://doi.org/10.1016/j.quascirev.2006.08.003

Thouveny, N., Bourlès, D. L., Saracco, G., Carcaillet, J. T., & Bassinot, F. (2008). Paleoclimatic context of geomagnetic dipole lows and excursions in the Brunhes, clue for an orbital influence on the geodynamo? *Earth and Planetary Science Letters*, **275**(3–4), 269–284. https://doi.org/10.1016/j.epsl.2008.08.020

Timmer, J., & König, M. (1995). On generating power law noise. *Astronomy & Astrophysics*, **300**, 707–710.

Valet, J.-P., Meynadier, L., & Guyodo, Y. (2005). Geomagnetic dipole strength and reversal rate over the past two million years. *Nature*, **435**(7043), 802–805. https://doi.org/10.1038/nature03674

Valet, J.-P., Moreno, E., Bassinot, F., Johannes, L., Dewilde, F., Bastos, T., et al. (2011). Isolating climatic and paleomagnetic imbricated signals in two marine cores using principal component analysis. *Geochemistry, Geophysics, Geosystems*, **12**(8), Q08012. https://doi.org/10.1029/2011gc003697

Venz, K. A., Hodell, D. A., Stanton, C., & Warnke, D. A. (1999). A 1.0 myr record of Glacial North Atlantic Intermediate Water variability from ODP site 982 in the northeast Atlantic. *Paleoceanography*, **14**(1), 42–52. https://doi.org/10.1029/1998pa900013

von Blanckenburg, F., & Bouchez, J. (2014). River fluxes to the sea from the ocean's $^{10}$Be/$^{9}$Be ratio. *Earth and Planetary Science Letters*, **387**, 34–43. https://doi.org/10.1016/j.epsl.2013.11.004

von Blanckenburg, F., Bouchez, J., & Wittmann, H. (2012). Earth surface erosion and weathering from the $^{10}$Be (meteoric)/$^{9}$Be ratio. *Earth and Planetary Science Letters*, **351**, 295–305. https://doi.org/10.1016/j.epsl.2012.07.022

von Blanckenburg, F., & Igel, H. (1999). Lateral mixing and advection of reactive isotope tracers in ocean basins: Observations and mechanisms. *Earth and Planetary Science Letters*, **169**(1–2), 113–128. https://doi.org/10.1016/s0012-821x(99)00070-9

von Blanckenburg, F., & Nägler, T. F. (2001). Weathering versus circulation-controlled changes in radiogenic isotope tracer composition of the Labrador Sea and North Atlantic Deep Water. *Paleoceanography*, **16**(4), 424–434. https://doi.org/10.1029/2000pa000550

von Blanckenburg, F., & O'nions, R. (1999). Response of Beryllium and radiogenic isotope ratios in Northern Atlantic Deep Water to the onset of northern hemisphere glaciation. *Earth and Planetary Science Letters*, **167**(3–4), 175–182. https://doi.org/10.1016/s0012-821x(99)00028-x

von Blanckenburg, F., O'nions, R., Belshaw, N., Gibb, A., & Hein, J. (1996). Global distribution of Beryllium isotopes in deep ocean water as derived from Fe-Mn crusts. *Earth and Planetary Science Letters*, **141**(1–4), 213–226. https://doi.org/10.1016/0012-821x(96)00059-3

Waelbroeck, C., Skinner, L. C., Labeyrie, L., Duplessy, J.-C., Michel, E., Vazquez Riveiros, N., et al. (2011). The timing of deglacial circulation changes in the Atlantic. *Paleoceanography*, **26**(3), PA3213. https://doi.org/10.1029/2010pa002007

Watkins, S., & Maher, B. (2003). Magnetic characterisation of present-day deep-sea sediments and sources in the North Atlantic. *Earth and Planetary Science Letters*, **214**(3–4), 379–394. https://doi.org/10.1016/s0012-821x(03)00422-9

Weltje, G. J., & Tjallingii, R. (2008). Calibration of XRF core scanners for quantitative geochemical logging of sediment cores: Theory and application. *Earth and Planetary Science Letters*, **274**(3–4), 423–438. https://doi.org/10.1016/j.epsl.2008.07.054

Wittmann, H., von Blanckenburg, F., Mohtadi, M., Christl, M., & Bernhardt, A. (2017). The competition between coastal trace metal fluxes and oceanic mixing from the $^{10}$Be/$^{9}$Be ratio: Implications for sedimentary records. *Geophysical Research Letters*, **44**(16), 8443–8452. https://doi.org/10.1002/2017gl074259

Xuan, C., & Channell, J. E. (2008). Origin of orbital periods in the sedimentary relative paleointensity records. *Physics of the Earth and Planetary Interiors*, **169**(1–4), 140–151. https://doi.org/10.1016/j.pepi.2008.07.017

Zhou, W., Kong, X., Paterson, G. A., Sun, Y., Wu, Y., Ao, H., et al. (2023). Eccentricity-paced geomagnetic field and monsoon rainfall variations over the last 870 kyr. *Proceedings of the National Academy of Sciences of the United States*, **120**(17), e2211495120. https://doi.org/10.1073/pnas.2211495120
Zhou, W., Priller, A., Beck, J. W., Zhengkun, W., Maobai, C., Zhisheng, A., et al. (2007). Disentangling geomagnetic and precipitation signals in an 80-hyr Chinese Loess record of $^{10}$Be. *Radiocarbon*, **49**(1), 137–158. https://doi.org/10.1017/s0033822200041977

Ziegler, L., Constable, C., Johnson, C., & Tauxe, L. (2011). PADM2M: A penalized maximum likelihood model of the 0–2 ma palaeomagnetic axial dipole moment. *Geophysical Journal International*, **184**(3), 1069–1089. https://doi.org/10.1111/j.1365-246x.2010.04905.x

## References From the Supporting Information

Banerjee, S. K., King, J., & Marvin, J. (1981). A rapid method for magnetic granulometry with applications to environmental studies. *Geophysical Research Letters*, **8**(4), 333–336. https://doi.org/10.1029/gl008i004p00333

Channell, J. E., Singer, B. S., & Jicha, B. R. (2020). Timing of quaternary geomagnetic reversals and excursions in volcanic and sedimentary archives. *Quaternary Science Reviews*, **228**, 106114. https://doi.org/10.1016/j.quascirev.2019.106114

Coe, R. S., Grommé, S., & Mankinen, E. A. (1978). Geomagnetic paleointensities from radiocarbon-dated lava flows on Hawaii and the question of the Pacific nondipole low. *Journal of Geophysical Research*, **83**(B4), 1740–1756. https://doi.org/10.1029/jb083ib04p01740

Herrero-Bervera, E., Helsley, C., Sarna-Wojcicki, A., Lajoie, K., Meyer, C., McWilliams, M., et al. (1994). Age and correlation of a paleomagnetic episode in the western United States by $^{40}$Ar/$^{39}$Ar dating and tephrochronology: The Jamaica, Blake, or a new polarity episode? *Journal of Geophysical Research*, **99**(B12), 24091–24103. https://doi.org/10.1029/94jb01546

Laj, C., & Channell, J. E. (2007). Geomagnetic excursions. *Geomagnetism*, **5**, 373–416. https://doi.org/10.1016/b978-044452748-6.00095-x

Lanphere, M. A. (2000). Comparison of conventional K-Ar and $^{40}$Ar/$^{39}$Ar dating of young mafic volcanic rocks. *Quaternary Research*, **53**(3), 294–301. https://doi.org/10.1006/qres.1999.2122

Liu, J., Nowaczyk, N. R., Frank, U., & Arz, H. W. (2018). A 20–15 ka high-resolution paleomagnetic secular variation record from Black Sea sediments–no evidence for the ‘Hilina Pali excursion’? *Earth and Planetary Science Letters*, **492**, 174–185. https://doi.org/10.1016/j.epsl.2018.04.014

Lund, S. (2022). A summary of paleomagnetic secular variation and excursions for the last 380 ky of the Brunhes normal polarity Chron. *Physics of the Earth and Planetary Interiors*, **333**, 106923. https://doi.org/10.1016/j.pepi.2022.106923

Lund, S., Stoner, J. S., Channell, J. E., & Acton, G. (2006). A summary of Brunhes paleomagnetic field variability recorded in Ocean Drilling Program cores. *Physics of the Earth and Planetary Interiors*, **156**(3–4), 194–204. https://doi.org/10.1016/j.pepi.2005.10.009

Maher, B. A. (1988). Magnetic properties of some synthetic sub-micron magnetites. *Geophysical Journal International*, **94**(1), 83–96. https://doi.org/10.1111/j.1365-246x.1988.tb03429.x

Robinson, S. G., Maslin, M. A., & McCave, I. N. (1995). Magnetic susceptibility variations in Upper Pleistocene deep-sea sediments of the NE Atlantic: Implications for ice rafting and paleocirculation at the last glacial maximum. *Paleoceanography*, **10**(2), 221–250. https://doi.org/10.1029/94pa02683

Singer, B., Relle, M., Hoffman, K., Battle, A., Laj, C., Guillou, H., & Carracedo, J. (2002). Ar/Ar ages from transitionally magnetized lavas on La Palma, Canary Islands, and the geomagnetic instability timescale. *Journal of Geophysical Research*, **107**(B11). EPM 7-1-EPM 7-20. https://doi.org/10.1029/2001jb001613

Singer, B. S., Guillou, H., Jicha, B. R., Laj, C., Kissel, C., Beard, B. L., & Johnson, C. M. (2009). $^{40}$Ar/$^{39}$Ar, K-Ar and $^{230}$Th-$^{238}$U dating of the Laschamp excursion: A radioisotopic tie-point for ice core and climate chronologies. *Earth and Planetary Science Letters*, **286**(1–2), 80–88. https://doi.org/10.1016/j.epsl.2009.06.030

Singer, B. S., Guillou, H., Jicha, B. R., Zanella, E., & Camps, P. (2014). Refining the quaternary geomagnetic instability time scale (GITS): Lava flow recordings of the Blake and Post-Blake excursions. *Quaternary Geochronology*, **21**, 16–28. https://doi.org/10.1016/j.quageo.2012.12.005

Singer, B. S., Jicha, B. R., He, H., & Zhu, R. (2014). Geomagnetic field excursion recorded 17 ka at Tianchi Volcano, China: New $^{40}$Ar/$^{39}$Ar age and significance. *Geophysical Research Letters*, **41**(8), 2794–2802. https://doi.org/10.1002/2014gl059439

Tauxe, L., Pick, T., & Kok, Y. (1995). Relative paleointensity in sediments: A pseudo-Thellier approach. *Geophysical Research Letters*, **22**(21), 2885–2888. https://doi.org/10.1029/95gl03166

Thouveny, N., Carcaillet, J., Moreno, E., Leduc, G., & Nérini, D. (2004). Geomagnetic moment variation and paleomagnetic excursions since 400 kyr BP: A stacked record from sedimentary sequences of the Portuguese margin. *Earth and Planetary Science Letters*, **219**(3–4), 377–396. https://doi.org/10.1016/s0012-821x(03)00701-5