Neodymium isotopes as a paleo-water mass tracer: A model-data reassessment

Proxy reconstructions from deep ocean sediments have helped to shape our understanding of the role of the global overturning circulation in past climate change. Neodymium (Nd) isotopes have contributed to this knowledge, as a tracer of past bottom water provenance and mixing. Here, we extend the imple- mentation of Nd isotopes in the physical-biogeochemical Bern3D model by revising a number of critical parameterizations, which result in an improved description of the marine Nd cycle. We exploit the dynamically consistent framework of the model, which allows us to assess the processes driving non-conservative Nd isotope behavior with a particular focus on the Last Glacial Maximum (LGM) and its substantially different climatic, oceanic, and biogeochemical boundary conditions. We show that the more radiogenic Nd isotopic compositions found throughout the glacial ocean can be explained by changes in the weathering input ﬂ uxes and do not require large reorganizations of the deep circulation. Our ﬁ ndings further highlight that the Nd isotopic composition of a water mass can not only be signi ﬁ cantly affected by a benthic Nd ﬂ ux, but also be modi ﬁ ed by the vertical downward transport of Nd via reversible scavenging. While these non-conservative processes only have a limited impact in the modern ocean, they were substantially more pronounced during the LGM and mostly independent of the circulation state, with their contributions being non-linear, partially opposing, and spatially variable. During the transiently simulated deglaciation Nd isotope variations induced by major circulation weakenings and resumptions are found to be most pronounced in the South Atlantic, while they are increasingly muted towards the north. Hence, it emerges that the interpretation of authigenic Nd isotope records requires more spatially speci ﬁ c considerations of non-conservative processes in order to more reliably infer basin-scale ocean circulation and water mass mixing of the past. (http://creativecommons.org/licenses/by/4.0/).


Introduction
Over the past decades, neodymium (Nd) isotopes emerged as a valuable geochemical tracer of water mass provenance (e.g., Piotrowski et al., 2004;P€ oppelmeier et al., 2020b;Roberts et al., 2010). Since its radiogenic isotope composition is independent of environmental fractionation and archived in marine sediments (Frank, 2002;van de Flierdt et al., 2016), Nd isotopes provide a unique opportunity to assess the origin of water masses today and to reconstruct the physical state of the ocean in the past. In combination with other paleoceanographic tracers, this proxy promises to be a key element in disentangling biologically-mediated from physical processes taking part in the transfer of carbon from the surface to the deep ocean (Brovkin et al., 2007;Sigman and Boyle, 2000). Hence, Nd isotopes may help to better understand the individual contributions of both these processes to CO 2 draw-down or release under changing climatic boundary conditions, not only during past glacial-interglacial cycles, but also for the storage of anthropogenic carbon.
Neodymium is a trace metal that is supplied to the ocean through weathering of continental rocks in different forms, such as mineral dust dissolution, dissolved and particulate Nd by rivers (Frank, 2002), and from sediments as a benthic flux (Abbott et al., 2015). In the water column Nd adsorbs to sinking particles, is partly released back into dissolved form in the deep, but is ultimately buried in the sediments (Siddall et al., 2008). Neodymium has a mean ocean residence time of 500e800 years (P€ oppelmeier et al., 2020c;Tachikawa et al., 2003), which is shorter than that of the global ocean circulation. The radiogenic Nd isotopic composition, denoted as εNd (which is the 143 Nd/ 144 Nd deviation from the Chondritic Uniform Reservoir per 10,000; Jacobsen and Wasserburg, 1980), of continental rocks is a function of age and initial Sm/Nd ratio, and thus εNd shows pronounced global variability (Robinson et al., 2021). The old cratons of North America and Greenland provide Nd with low isotopic signatures (termed unradiogenic) to the North Atlantic while input into the Pacific is dominated by young volcanic material with high (termed radiogenic) Nd isotopic compositions. This leads to a distinct inter-basin gradient between the Atlantic and Pacific with intermediate values in the Indian and Southern Ocean (Tachikawa et al., 2017;van de Flierdt et al., 2016) ( Fig. 1). These characteristics of the marine Nd (isotope) cycle are fundamental for the utility of εNd as a powerful tool to investigate global-scale ocean circulation. Furthermore, its isotopic composition of past bottom water can be extracted from sedimentary archives opening the door for many paleoceanographic applications (e.g., Blaser et al., 2016;Gutjahr et al., 2007;Wilson et al., 2013). However, critical knowledge gaps remain in the understanding of modern marine Nd cycling, such as the driving processes and the magnitude of the benthic Nd flux (Abbott, 2019;Abbott et al., 2015;Haley et al., 2017), Nd release from particulate phases transported by rivers (Rousseau et al., 2015), and the nature of reversible scavenging (Siddall et al., 2008;Stichel et al., 2020). Additional uncertainties further arise for the application in the paleoceanographic context. These include a limited understanding of the variations of input fluxes from all sources and the impact of substantially different seawater particle concentrations on the scavenging efficiency. As a consequence, the evolution of the εNd endmembers in the North Atlantic, Southern Ocean, and North Pacific remain insufficiently constrained even though recent studies have made important progress on this front (P€ oppelmeier et al., 2020b, 2020cZhao et al., 2019). Yet more importantly, no information on past seawater Nd concentrations is available to date due to missing suitable archives. Thus, for the interpretation of εNd-based water mass mixing calculations past Nd concentrations are commonly assumed constant and equal to modern (e.g., Howe et al., 2016a,b;P€ oppelmeier et al., 2020b). While this assumption is reasonable in the absence of any constraining information, it is unlikely to be correct. In particular, on glacial-interglacial timescales it is expected that the dramatic changes in the climatic conditions also changed the continental weathering flux, very likely equally affecting the runoff of Nd, and as a consequence the ocean's Nd inventory. For instance, during the Last Glacial Maximum (LGM) the large continental northern hemispheric ice sheets are thought to be associated with limited Nd weathering and transport to the North Atlantic, which is one of the leading hypotheses for a more radiogenic North Atlantic εNd end-member during that time (Blaser et al., 2020;P€ oppelmeier et al., 2019, 2020cZhao et al., 2019). The evidence for striking changes in the global dust deposition (Lambert et al., 2008(Lambert et al., , 2012Mahowald et al., 2006) as well as biogenic export fluxes during the last glacial cycle (Anderson et al., 2014;Kumar et al., 1995) must have further led to changes in seawater particle concentrations and consequently in Nd scavenging. The issue of unknown past Nd concentrations is further exacerbated by the fact that major modern subsurface water masses exhibit Nd concentrations that differ by a factor of up to 4 (~10 pmol/kg for Antarctic Intermediate Water versus~40 pmol/kg for Pacific Deep Water; van de Flierdt et al., 2016). This leads to a large impact of changing Nd concentrations on the resulting Nd isotopic signatures during water mass mixing.
It thus becomes clear that the interpretation of reconstructed εNd requires careful consideration of the respective climate background states and their influence on the marine Nd cycle. For the quantification of these climatic impacts on Nd cycling dedicated models are critically necessary. Approaches with simple box models give first-order estimates on these impacts (e.g., Du et al., 2018Du et al., , 2020P€ oppelmeier et al., 2020b), but they generally lack important aspects such as dynamical biogeochemical cycling and realistic physical descriptions of ocean circulation and climate interactions. Neodymium has also been implemented in a number of Earth system models of intermediate complexity (EMIC) that capture a more complete view of the Nd cycle (e.g., Arsouze et al., 2009;Gu et al., 2019;Rempfer et al., 2011). Yet, most of these implementations were limited by assumptions on the Nd cycle that do not reflect the current knowledge about dominating processes (most notably the likely incorrect assumption that the benthic fluxes are restricted to shallow margins, Abbott et al., 2015). Sensitivity tests with these EMICs were further limited to simplified experiments such as globally altered Nd fluxes or isotopic signatures (e.g., Gu et al., 2019;P€ oppelmeier et al., 2020c;Rempfer et al., 2012). As a result, none of these EMICs was successful in simulating the glacial εNd distribution in convincing agreement with the reconstructed data on global scale (see Du et al., 2020 for a review).
Here, we use the Bern3D EMIC with an updated Nd implementation that alleviates many earlier applied simplifications and allows for a more comprehensive interpretation of reconstructed Nd isotope records in particular from times with considerably different climate boundary conditions. After validating our Nd implementation against seawater observations, we diagnose the contributions of conservative versus non-conservative effects to the simulated modern Atlantic εNd distribution. With this more complete interpretative framework at hand we then simulate the global steady-state Nd cycle under LGM boundary conditions also considering changes in biogeochemistry and Nd fluxes. This enables us to test the agreement between various circulation regimes and reconstructed εNd of the LGM and ultimately also the capability of the Nd isotope tracer as a proxy for water mass provenance in a vastly different climate regime. Finally, we perform transient simulations from the LGM to present day, with particular focus on the deglacial perturbations of the Atlantic Meridional Overturning Circulation (AMOC) associated with the iceberg discharge events of Heinrich Stadial 1 and the Younger Dryas. These simulations thus not only inform us on the influence of changing Nd fluxes but also the impact of transient water mass changes and their profound biogeochemical implications on the Nd cycle.
This study therefore aims at elucidating modern and past marine Nd cycling and building a consistent and comprehensive framework toward more robust interpretations of Nd isotope records.

Methods
We here combine published reconstructions of past seawater εNd from marine sediments covering the past 25 kyr with simulations of the Nd-enabled Bern3D model. The Bern3D EMIC (v2.0) consists of a dynamic geostrophic-frictional balance ocean (Müller et al., 2006) coupled to a thermodynamic sea-ice component, and a single-layer energy-moisture balance atmosphere (Ritz et al., 2011). All model components have a spatial resolution of 41 Â 40 grid cells with the ocean having 32 logarithmically scaled depth layers. The physical ocean features an isopycnal diffusion scheme and a Gent-McWilliams parameterization for eddy-induced transport (Griffies, 1998). Wind stress and cloud cover are prescribed as monthly climatologies (Kalnay et al., 1996). We further use an artificial dye tracer in order to unambiguously diagnose the expansion of North Atlantic Deep Water (NADW). The dye tracer is restored to 1 in the surface of the North Atlantic (40e70 N) and is reset to 0 at all other surfaces.
The prognostic biogeochemistry module calculates export production and remineralization in the water column of particulate organic matter (POM), calcium carbonate (CaCO 3 ), and biogenic opal (Battaglia and Joos, 2018;Parekh et al., 2008;Tschumi et al., 2011). Primary productivity is limited by light availability, sea-ice extent, and phosphate and iron concentrations, and for biogenic opal also by silicic acid concentrations. The new production is split into two carbon pools, with two thirds transferred to dissolved organic carbon and one third exported as POM. The dust flux is prescribed after Mahowald et al. (2006) (see Figure S1). For details on the biogeochemical implementation we refer to Parekh et al. (2008) and Tschumi et al. (2011).
The model was integrated over 35,000 years to a pre-industrial steady-state corresponding to 1765 CE boundary conditions. The additional initialization to the LGM state is described in section 2.2.

Updated Nd parameterization in the Bern3D model
The Nd cycle has been fully implemented in the Bern3D model (P€ oppelmeier et al., 2020c;Rempfer et al., 2011), which includes the three sources of dust, rivers, and marine sediments and the internal cycling and removal by reversible scavenging. Since the last update of Nd isotopes in the Bern3D model (P€ oppelmeier et al., 2020c) a number of new observational studies were published that investigated different aspects of the Nd cycle in detail.
Most recently, Robinson et al. (2021) published a new highresolution Nd isotope map of the global seafloor that is Sketch of a section from the North to South Pacific through the Southern Ocean (SO) further extending from the South to the North Atlantic. Black arrows roughly indicate modern water mass propagation of the meridional overturning. Small colored arrows depict Nd inputs to the ocean with colors indicating their respective isotopic composition (lighter and darker colors indicating more and less radiogenic Nd input, respectively). Red symbols mark the locations of Nd isotope sediment records the model outputs of section 5 are compared to. Distributions of (b) simulated Nd isotopic compositions and (c) Nd concentrations of the pre-industrial control. Colored circles depict measurements of seawater that are close to the section (±15 ) on the same color scale as the background. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) constrained by more than 5100 detrital εNd data points and the lithology of seafloor sediments ( Figure S2a). We now use this map as the basis for the isotopic signatures of the benthic Nd flux but slightly adjusted it in the eastern and northern Pacific, where the bulk detritus is thought to be not representative of the Nd isotopic signatures of the benthic flux (Du et al., 2020). In that region we increased the detrital Nd isotopic signatures by 4 ε-units. Further, based on the rationale that authigenic and detrital phases both exchange their Nd in the pore water during early diagenesis (e.g., Du et al., 2016;Blaser et al., 2019), we adopted a new parameterization for the calculation of the εNd signatures of the benthic flux. For this, the weighted averages of the above mentioned detrital map, and the bottom water signatures are dynamically calculated at all bottommost grid cells and at every time step in a 2:1 ratio (a 1:2 ratio has also been tested and is discussed in section 4, see also Figure S10). Finally, we use the regional benthic flux scaling as described in P€ oppelmeier et al. (2020c) (their Fig. 3d), but with the exception that we here set the benthic Nd flux in the entire Arctic to zero, because Nd concentrations were far too high in this region in the previous implementation compared to observations (P€ oppelmeier et al., 2020c), which was also found for a recent implementation of chromium in the Bern3D model that uses a similar description for the benthic flux (P€ oppelmeier et al., 2021b).
New constraints also became available for the riverine Nd source of the Congo River (Rahlf et al., 2021) that indicate large Nd release from particulate matter in the estuary in addition to already very high dissolved Nd concentrations. Thus, the total Nd source of the Congo River is estimated to be almost twice as high as the one of the Amazon River despite its discharge rate being only one fifth (Rahlf et al., 2021).
We further take this opportunity to additionally revise the Nd isotope signatures of the dust source that were previously based on first-order assumptions by Tachikawa et al. (2003) assigning one isotope signature to each ocean basin (see Rempfer et al., 2011). For the revised εNd dust map presented here ( Figure S2d) we take into account the dust plume expansion as constrained by Mahowald et al. (2006) and the mean Nd isotopic signatures of the dust sources as constrained from detrital data in the respective dust source areas (Blanchet, 2019). The polar regions are linearly extrapolated without any additional constraints, but this does not impact the effective Nd dust source since the dust fluxes in these regions are virtually absent.
Additionally, we implemented a new parameterization of the particle sinking velocity responsible for the downward transport of Nd via reversible scavenging. Previously, particle sinking was prescribed with a globally uniform velocity of 1000 m/yr. This simplified description was identified as the main driver for the too low simulated surface concentrations compared to observations not only in the Bern3D model (P€ oppelmeier et al., 2020c;Rempfer et al., 2011) but also in CESM (Gu et al., 2019), which uses the same description as Bern3D. As such, we developed a new parameterization that takes into account the lower sinking velocities within the mixed layer due to enhanced turbulence (e.g., Chamecki et al., 2019;Noh et al., 2006) and complexation with organic ligands (e.g., Byrne and Kim, 1990). We describe these two processes combined as a first-order approximation by introducing a density dependent scaling of the particle sinking velocity as follows: À r min Þ Â ð1 À aÞ ðr max À r min Þ where W s,const is the globally uniform and constant particle sinking velocity of 1000 m/yr used previously, a is the minimum of scaling here set to 0.2, and r(k), r min , and r max are the densities at depth layer k, the minimum and maximum in the water column, respectively. Hence, particle settling velocity vary between 200 and 1000 m/yr throughout the water column. Finally, we performed a re-tuning of the Nd module analogous to P€ oppelmeier et al. (2020c). For this, we updated the database of observational seawater Nd concentrations and isotopic signatures with data from most recent studies ( Figure S3a; Behrens et al., 2018;Dausmann et al., 2020;Fuhr et al., 2021;Rahlf et al., 2021;Stichel et al., 2018;Wu, 2019). The details of the tuning and a comparison to P€ oppelmeier et al. (2020c) are described in supplementary Text S1 and Table S1.

LGM forcings
For the initialization of the LGM state we first spun up the model analogous to the pre-industrial except for a closed Bering Strait (see P€ oppelmeier et al., 2021a). Then continuing from this pre-industrial state we ran the model for 20 kyr under constant LGM forcing. This forcing includes orbital parameters set to 20 kyr BP (Berger, 1978) and radiative forcing corresponding to greenhouse gas concentrations of CO 2 ¼ 191 ppm, CH 4 ¼ 370 ppb, and N 2 O ¼ 208 ppb (K€ ohler et al., 2017). Further, the continental ice sheet extent and the related changes in albedo are set to the reconstructions by Peltier (2004). The corresponding freshwater relocation responsible for the increase in salinity by about 1 PSU is performed during the first 2 kyr of the LGM spin-up, which gives the ocean enough time to reach a new steady-state. We additionally apply LGM wind stress anomalies (Muglia and Schmittner, 2015) and account for the increased vertical mixing due to lower sea-level increasing tidal dissipation (Wilmes et al., 2019) as described by P€ oppelmeier et al. (2021a).
These forcings alone lower the Global Mean Surface Temperature (GMST) and Mean Ocean Temperature (MOT) by only 3.2 C and 1.7 C, respectively (P€ oppelmeier et al., 2021a), while recent reconstruction suggest coolings of~6.1 C for GMST (Tierney et al., 2020) and 2.6 C for MOT (Bereiter et al., 2018). We therefore apply an additional top-of-the-atmosphere aerosol forcing of À2.5 W/m 2 to all LGM simulations and increase the continental ice sheet albedo (only affecting the Laurentide and Fennoscandian ice sheets) from 0.7 to 0.85. We acknowledge that such a strong change in the aerosol forcing is at the very upper limit of previous estimates (Hopcroft et al., 2015), and we therefore consider this here only as a tuning mechanism for achieving colder LGM temperatures.
At the LGM many aspects of the biogeochemical cycling experienced profound changes as well, yet some of these are not dynamically simulated in the Bern3D model and need to be prescribed externally. Therefore, we reduced the POM remineralization in the upper water column as a consequence of lower microbial activity due to the lower ocean temperatures ( Figure S4; Laufk€ otter et al., 2017). Further, we here exchanged the prescribed modern dust field with the LGM one from the same study by Mahowald et al. (2006) (Figure S1d). This has a substantial direct impact on the Nd dust flux, which is increased by a factor of three relative to the pre-industrial. We additionally slightly reduced the global riverine flux by 10% to reflect the drier conditions of the LGM (e.g., Russell et al., 2014;Wang et al., 2017). Yet, we strongly reduced the riverine Nd flux into the North Atlantic (90% reduction north of 30 N) in order to represent the strongly limited weathering at the northern high latitudes that was restricted by the large continental ice sheets (Zhao et al., 2019). Since the North Atlantic benthic flux is thought to be controlled to a significant extent by the reactivity of the detrital material primarily derived from the North American continent (Howe et al., 2016a,b;P€ oppelmeier et al., 2019), which was locked under the ice sheet during the LGM we also reduced the elevated benthic flux in the North Atlantic to the background level of 5 pmol/cm 2 /yr (see P€ oppelmeier et al., 2020c). The choices of these parameterizations of the past Nd input fluxes are substantiated by geochemical data indicating reduced weathering input from the North American continent to the North Atlantic during the last glacial. For instance, the consistently more radiogenic Nd isotope signatures found throughout the entire glacial North Atlantic indicate a reduced influence from the old cratons of North America and Greenland (Howe et al., 2016a,b;P€ oppelmeier et al., 2019). More directly, lead isotope data from around the western North Atlantic show the first increase of continental weathering only during the initiation of the last deglaciation (e.g., Gutjahr et al., 2009;Crocket et al., 2012). A detailed list of all changes in boundary conditions can be found in supplementary Table S2.

Transient deglacial forcing
Finally, we also performed transient simulations starting from the LGM steady-state that were run for 20 kyr covering the entire deglaciation and Holocene. The orbital forcing follows Berger (1978) and greenhouse gas evolutions are prescribed after K€ ohler et al. (2017). Continental ice sheets are scaled by the global benthic d 18 O stack of Lisiecki and Stern (2016), which is also used to scale the LGM wind stress and vertical mixing anomalies, the aerosol forcing, the dust field, and all biogeochemical variables that were externally prescribed for the LGM runs including the changes in Nd source fluxes ( Figure S11). Changes in the AMOC were forced by North Atlantic (45e70 N) freshwater fluxes following Menviel et al. (2012). While this is a simplification for some parameters, such as the changes in the POM remineralization, it is consistent with the intermediate complexity of the Bern3D model and should only add some small temporal uncertainty to the evolution of these parameters. All these parameters were updated every 50 years to reduce computational overhead.
At 12 kyr BP we open the Bering Strait following evidence from local sea-level reconstructions (Pico et al., 2020). Further, from 10.5 to 8 kyr BP we apply an additional Nd-specific forcing that increases the benthic Nd flux in the northern North Atlantic fivefold. As such, we attempt to capture the εNd anomalies reconstructed in many Atlantic records during that time, which were already investigated in detail with the Bern3D model by P€ oppelmeier et al. (2020c).

Simulated Nd isotopes in the modern ocean
The observed modern εNd distribution of the Atlantic agrees well with diagnosed water masses from conservative tracers on basin scale, while non-conservative processes appear to have a greater impact in the Pacific (Du et al., 2020;Horikawa et al., 2011;Tachikawa et al., 2017). It is thought that the benthic flux, altering the bottom water isotope signatures and Nd concentration (or [Nd]) independent of water mass mixing, has a larger influence on the sluggishly circulating Pacific than on the vigorously circulating Atlantic (Du et al., 2020). However, many unknowns remain regarding other non-conservative processes that may only have a limited impact today and are therefore poorly quantified, but may have been significant in the past. As it emerged from early on that Nd isotopes behave more conservative in the Atlantic than the Pacific, the majority of studies focused on the former for the reconstruction of water mass provenance (e.g., Howe et al., 2016a,b;Piotrowski et al., 2004;P€ oppelmeier et al., 2019;Roberts et al., 2010). Hence, here we also focus our investigations on the Atlantic basin.

Model-data agreement of updated Nd parameterization
Many new constraints of the modern Nd cycle have been considered here for the first time in an EMIC, and we therefore briefly evaluate the performance of the updated Nd module in representing the ever-growing dataset of seawater observations of Nd isotopes and concentrations (see section 2.1). In the preindustrial (PI) control simulation εNd exhibits a pronounced inter-basin gradient with unradiogenic signatures in the deep North Atlantic steadily increasing toward the North Pacific (Fig. 1b). The three main intermediate and deep water masses of the Atlantic, Antarctic Intermediate Water (AAIW), Antarctic Bottom Water (AABW), and North Atlantic Deep Water (NADW), can all be clearly identified (Fig. 1b). However, due to the coarse grid resolution and simplified large-scale dynamics the model cannot resolve the individual constituents of NADW, i.e., Labrador Sea Water and the Nordic Seas Overflow Waters that are characterized by distinct Nd isotope signatures, which Results in a more homogeneous Nd isotope signature of NADW than observed today . We therefore define this water mass simply as northernsourced water (NSW) as it is common in paleoceanographic studies.
Compared to the previous implementation by P€ oppelmeier et al. (2020c), we find the largest deviations in the top 1 km, where local signals, mostly from riverine input, are now more pronounced ( Figure S5a). This is mainly a consequence of the revised parameterization of the particle sinking velocity that is now attenuated towards shallower depths leading to reduced scavenging and thus a longer residence time in the upper water column. As a result, Nd concentrations are reduced less at shallow waters compared to P€ oppelmeier et al. (2020c) and hence are in better agreement with observations. Deep water Nd concentrations generally also agree well with observations, but are somewhat too high in the deepest North Atlantic (Fig. 2). In sum, the updated Nd parameterization slightly improves the agreement of both, the Nd isotopic compositions and concentrations, between observations and the model compared to previous iterations and thus the revised model represents their global distributions more accurately.

How conservative is εNd in the modern atlantic?
In order to evaluate the (non-)conservative properties of εNd in the model we make use of a fully conservative dye tracer that tracks the expansion of NSW. This dye tracer is implemented in the model such that the fraction of NSW is 100% at the northern North Atlantic surface ocean (40e70 N) and is diluted elsewhere due to water mass mixing. To assess how conservative εNd behaves in the Atlantic, we compute the zonally averaged water mass distribution based on binary mixing of southern-sourced water (SSW) and NSW. Since Nd end-member characteristics cannot be defined at the surface due to local inputs, we choose SSW and NSW end-members at 60 S e 3.1 km water depth and 50 N e 2.6 km water depth, respectively (see supplementary text S2 for details). As such, we renormalized the %NSW distribution of the dye tracer to these endmember locations to allow for a direct comparison between εNd and dye-based %NSW. Yet, we note that this renormalization only makes a marginal difference for the dye tracer distribution.
As is common for such binary mixing calculations, we define the depth of NSW as the 50% isoline (e.g., Howe et al., 2016a,b;Oppo et al., 2018;P€ oppelmeier et al., 2020a). This is different from the AMOC depth, often diagnosed in models from the zero-isoline of the stream function (Fig. 3), but we note that both diagnostic parameters are closely correlated. In the pre-industrial control NSW reaches down to the seafloor in the North Atlantic, while the water mass boundary is situated at about 3.5 km in the South Atlantic. The Nd-based %NSW distribution generally displays the same features as the dye-based %NSW. Yet, the South Atlantic water mass boundary is situated about 500 m deeper and the fraction of NSW is overestimated by up to 20% in the entire deep Atlantic. At the same time NSW is slightly underestimated at intermediate depths ( Fig. 3d), while even shallower depths of <1 km are dominated by local input sources and do not allow for meaningful inferences of water mass mixing, as was already noted by observational studies (e.g., Rahlf et al., 2021;van de Flierdt et al., 2016).
This pattern of underestimated shallow and overestimated deep NSW is primarily a consequence of the reversible scavenging that transports Nd downward in the water column. Thus, the radiogenic Nd of AAIW is adsorbed on particles but quickly desorbs again in the underlying water, which in this case is upper NSW that consequently becomes more radiogenic. This has recently also been observed at the GEOTRACES transect GA10 in the South Atlantic where the AAIW Nd isotope signal contributed to the underlying upper NSW (Wang et al., 2021). The same occurs for unradiogenic Nd of NSW being transported to the underlying SSW via reversible scavenging. Yet, since Nd concentrations increase and particle concentrations decrease with water depth, this effect is expected to diminish in the deeper water column. Instead, we here find the largest deviation from conservative mixing in the deepest part of the water column. This is attributed to the second major driver of  non-conservative Nd alteration, namely the benthic flux. In most of the modern Atlantic the detrital εNd signatures of sediments are less radiogenic than the overlying bottom water (see Figure S2a). As such, the bottom water is slowly modified towards less radiogenic signatures independent of water mass mixing. Generally, the overprinting is a function of benthic exposure time and the difference of bottom water to detrital εNd (Du et al., 2018). This thus also explains why the deep North Atlantic deviates most from conservative mixing, as sluggish SSW slowly accumulates increasingly unradiogenic Nd on its way north.
To further quantify the effect of both non-conservative processes on the Nd isotope tracer, we derived a hypothetical fully conservative Atlantic εNd distribution by normalizing the dye tracer distribution to the Nd end-members ( Fig. 4; see supplementary Text S2 for details). This normalization reveals that the deviation from fully conservative εNd is in the range of ±0.5 ε-units for most of the intermediate and mid-depth Atlantic. Only in the deep Atlantic, this deviation reaches values of about À1 ε-unit, which corresponds to an overestimation of the NSW fraction bỹ 20% (excluding shallow depths and the Nordic Seas, where this calculation is not applicable). This is in good agreement with observation-based calculation of non-conservative εNd (Tachikawa et al., 2017). The εNd-based over-and underestimations of the NSW distribution diagnosed here therefore are relatively minor. However, it cannot be assumed that they were equally small in the past under different boundary conditions, when particle fluxes and circulation regimes were substantially different.

Implications of changes in Nd cycling at the LGM
The LGM control simulation with boundary conditions as described in section 2.2 produces a cooling of the GMST of 6.1 C and a MOT cooling of 2.4 C; both in good agreement with recent reconstructions (GMST: 6.1 C, Tierney et al., 2020 andMOT: 2.6 C, Bereiter et al., 2018). Compared to a recent study that investigated the glacial AMOC in the Bern3D model (P€ oppelmeier et al., 2021a), we here apply an additional aerosol forcing leading to the colder temperatures and thus also a weaker AMOC (18.5 Sv here versus 21.3 Sv by P€ oppelmeier et al., 2021a). The total Nd source to the oceans is comparable to the modern, yet with a substantially shifted apportionment of the individual contributions. The greatly increased glacial dust flux (Lambert et al., 2012;Mahowald et al., 2006) leads to nearly a tripling of the Nd dust source. In contrast, the global benthic source is slightly lower than in the PI control due to the curtailed North Atlantic flux (see Table S1). Similarly, the global riverine source is reduced by about 20% in the LGM runs. Thus, the glacial dust and riverine source contributions combined are of about the same size as the sedimentary source (22%, 30%, and 48%, respectively and in contrast to 9%, 32%, and 60%, respectively for the PI, P€ oppelmeier et al., 2020c).
In order to investigate the influence of different AMOC configurations on the Atlantic εNd distribution, we force the surface ocean with additional salt and freshwater fluxes that are continuously applied over the entire 20 kyr of the LGM spinups so that the ocean can fully equilibrate. By applying a freshwater flux to the North Atlantic, that is compensated for in the North Pacific with an equivalent salt flux, we are able to produce weakened and shoaled AMOC states (see P€ oppelmeier et al., 2021a for details). These freshwater fluxes represent continental melt water fluxes and hydrological processes not explicitly or only poorly represented in the model. In the same manner, we apply additional salt fluxes (i.e., negative freshwater fluxes) to the Weddell Sea, compensated for in the entire Southern Ocean, to stimulate a stronger formation of deep SSW. The latter mostly affects the AMOC depth and only to a lesser degree its strength. By combining and systematically varying both these forcing parameters we obtain 24 different circulation states, that span a range of AMOC strengths between 9.5 and 18.5 Sv and NSW depths between 2.7 and 4.1 km at 20 S ( Figure S6). Circulation states weaker than~9 Sv are unstable and lead to a full collapse of the AMOC. We here emphasize that these 24 simulations are forced by artificial freshwater fluxes and hence may not represent fully realistic circulation states.
This wide array of different LGM simulations primarily serves to provide further insights into the impact of coupled physical and biogeochemical changes on the Nd cycle and the ability of εNd to faithfully trace the water mass provenance of the past. We evaluate these 24 LGM circulation states by comparing the simulated Nd isotope distributions to reconstructed εNd and their anomalies with respect to the pre-industrial/late Holocene. For this we complemented the data compilation by Du et al. (2020) with the newest reconstructions (Bergh et al., 2021;P€ oppelmeier et al., 2020b, 2020cTachikawa et al., 2021). The model performance was evaluated by minimizing the cost function of the Mean Absolute Error (MAE) weighted by the grid cell volumes. The MAEs only span relatively small ranges of 0.9e1.2 and 0.8 to 1.1 ε-units for absolute values and anomalies, respectively (Fig. 5). A closer look at the global Nd isotope distribution of these simulations reveals that changes in εNd relative to the PI are dominated by a shift to more radiogenic values throughout the global deep ocean for all LGM simulations (Fig. 6). The greatest deviation towards more radiogenic signatures is found in the domain of NSW, which is a direct consequence of the strongly reduced North Atlantic weathering flux (riverine and sedimentary). This corroborates the hypothesis proposed by Zhao et al. (2019) that a reduction of weathering contribution caused by the North American ice sheet played a critical role in setting the North Atlantic εNd end-member. However, due to our model's inability to resolve the distinct signatures of the upper and lower NSW branches, we cannot rule out the possibility that shifts in the location of deep water formation contributed to the more radiogenic end-member as suggested by Du et al. (2020). A combination of both, changes in weathering fluxes and deep water formation, most likely contributed to the observed shift in the NSW end member. Further constraints from paleoceanographic reconstructions and (regional) climate models are required for a better understanding of this question. Yet, for the interpretation of glacial εNd variations downstream of the northern North Atlantic it is only of secondary importance how exactly the end-member signature was altered, and our interpretations thus remain valid independent of the ultimate cause of the more radiogenic glacial NSW.
At intermediate depths, changes are dominated by the increased release of dust-derived Nd. The greatest intensification of the dust flux occurs in the North Pacific and South Atlantic, where signals of the unradiogenic Asian dust and radiogenic Patagonian dust, respectively, are clearly imprinted into the surface to intermediate ocean Nd isotope signatures. Less radiogenic Nd isotope signatures are, however, not reconstructed for the North Pacific during the LGM suggesting that the influence of Asian dust may be overestimated in our model. Despite these general features found in all simulations, a trend in the MAE cost functions is apparent, that suggests a better model-data agreement with shallower and weaker AMOC states (Fig. 5). This trend is slightly more pronounced in the MAEs of the anomalies that show a clearer minimum at NSW depths at 20 S of around 2.8 km and AMOC strengths between 10 and 13 Sv (Fig. 5b) than the MAEs of the absolute values. In some key regions such as the North Atlantic and North Pacific critical differences between core top and bottom water εNd signatures have been observed (e.g., Blaser et al., 2016;Hu et al., 2016), but these offsets may be temporally constant if the pore water chemistry was rather invariable as well. In this case, the anomalies have more interpretive meaning by retaining more information on the past evolution than the absolute Nd isotope signatures, also in the context of model-data comparisons. Similar findings of better model-data agreement with shallower NSW were previously also made based on stable carbon isotopes (e.g., Gu et al., 2020;Muglia and Schmittner, 2021).
Yet, taking into account the uncertainties associated with reconstructing εNd from sedimentary archives (e.g., age and analytical uncertainties) and the relatively small difference between best and worst fits, we consider this trend as not decidedly robust. Moreover, the Bern3D model suffers from a number of shortcomings due to its intermediate complexity that may have biased these Results, such as the limited spatial resolution and its direct impact on deep water formation (strength and location), or the simplified biogenic particle remineralization that is globally uniform. Nevertheless, the question arises why all these different circulation states only yield relatively minor differences in their Nd isotope distributions (at least relative to the large differences in circulation). To better understand this, we further investigate three simulations with substantially different circulation states. These are the LGM control simulation with a strong and deep AMOC (in the following referred to as LGM_18Sv), a simulation with a nearly equally strong but much shallower AMOC (LGM_16Sv), and a simulation with a very weak and shallow AMOC (LGM_09Sv) (Figs. 5a and 7).
The only major differences between the εNd distributions of these three simulations are found in the South Atlantic (in fact for all 24 simulations), while changes are somewhat smaller in the North Atlantic and virtually absent elsewhere (Fig. 6). In the strong and deep AMOC regime (LGM_18Sv), vigorous NSW feeds the formation of AAIW that acquires highly radiogenic signatures by the dissolution of South American dust and voluminously penetrates northward. The tight coupling of NSW and AAIW formation is thus also responsible for the shoaled and nearly absent AAIW in the weak AMOC state (LGM_09Sv). As a consequence, DεNd (i.e., the εNd anomaly) in the intermediate depth South Atlantic is at a maximum in the former, close to zero in the latter case, and intermediate for LGM_16Sv. At first glance, this appears to be counter-intuitive, because the AMOC of the PI control shows more resemblance to the one of the LGM_18Sv and LGM_16Sv than the LGM_09Sv. However, the general shifts toward more radiogenic signatures of both, NSW and AAIW, need to be considered. Thus, for simulation LGM_09Sv the end-member change of NSW nearly fully compensates for the retracted AAIW so that the absolute Nd isotope signatures at intermediate depths are only marginally different from the PI. Reconstructions from the southern Brazil Margin (30 S) also show nearly constant εNd signatures at intermediate depths over the past 25 kyr (Howe et al., 2018;P€ oppelmeier et al., 2020a), which therefore support a circulation regime in which AAIW formation is strongly attenuated as simulated here in LGM_09Sv. Yet, the simulated εNd anomaly at intermediate depth is to an important degree also determined by the nature of the South Atlantic Nd dust source that supplies highly radiogenic Nd to AAIW. As such, a smaller Nd concentration of dust, less dissolution, or a shift in the source region away from a Patagonian source (see Basile et al., 1997) could strongly diminish the supply of radiogenic Nd, which would also result in less radiogenic AAIW. Consequently, stronger AAIW production than in LGM_09Sv could also fit the reconstructed εNd from the Brazil Margin equally well under such conditions.
In the deep South Atlantic the difference between the simulations is smaller, but εNd signatures are generally more radiogenic   under the shallow AMOC states of LGM_16Sv and LGM_09Sv due to the greater intrusion of radiogenic SSW. In order to diagnose the true expansion of SSW and how this relates to the εNd distributions in these simulations we again calculate the dye-based NSW expansion analogously to the PI (Fig. 7def). For LGM_09Sv, this reveals a strong shoaling of NSW to 2.7 km and 3.2 km water depth at 20 S and 20 N, respectively (Fig. 8) and slightly less shoaling for LGM_16Sv. In contrast, the εNd-based %NSW distribution paints a different picture that also includes a strong shoaling in the South but instead NSW still dominating the entire North Atlantic at all depths (i.e., fraction greater than 50%) for both LGM_16Sv and LGM_09Sv. Here, the overestimation of NSW in the deep grows from about 10% at 30 S to more than 40% north of 30 N (Fig. 7h) compared to a maximum of 20% in the PI control. Similar is true for LGM_18Sv, only that the fraction of NSW is already overestimated by 40% at the equator, and additionally, %NSW is strongly underestimated in the upper water column. The other 21 LGM simulations exhibit the same behavior with vertical εNd and water mass gradients plotting between LGM_18Sv and LGM_09Sv (Fig. 8).
Hence, this indicates that in the Atlantic εNd behaves substantially less conservative under LGM than PI boundary conditions, even with vigorous overturning circulations. In the past, the investigation of non-conservative behavior of εNd mainly focused on the impact of the benthic flux and was generally assessed with simple box models (e.g., Du et al., 2018;Haley et al., 2017;P€ oppelmeier et al., 2020b). Yet, the Results of the Fig. 7. (a,b,c) Zonally integrated Atlantic stream function, (d,e,f) fraction of northern-sourced water (NSW) derived from the dye tracer, (g,h,i) fraction of NSW based on εNd, and (j,k,l) difference between panels c and e, and f and g, respectively (indicating where εNd over-or underestimates the NSW fraction) for simulations LGM_18Sv (left), LGM_16Sv, and LGM_09Sv (right). LGM simulations performed here suggest a more complex interplay of processes being responsible for the non-conservative behavior, because substantial non-conservative effects also occur under strong circulation regimes with low benthic exposure times, are not limited to the deep ocean, and can also lead to a large underestimation of %NSW at intermediate depth. Furthermore, simulated end-member shifts of SSW and NSW towards more radiogenic values are not equal in magnitude here, which results in a reduction of the maximum meridional εNd gradient and thus increases the impact of non-conservative processes on binary mixing for these LGM simulations. In contrast, previous studies assumed a constant (P€ oppelmeier et al., 2020b) or even greater (Howe et al., 2016a,b) glacial εNd end-member difference, based on highly radiogenic signatures reconstructed from the deep South Atlantic (À4.0 to À5.5; Howe et al., 2016a,b;Huang et al., 2020;Skinner et al., 2013). Here we are, however, unable to simulate such radiogenic Nd isotope signatures in the South Atlantic, since the vigorous mixing of the Antarctic circumpolar current (ACC) strongly attenuates any radiogenic input signal (even though the ACC strength is substantially reduced in all LGM simulations compared to PI due to the northward shifted westerlies and sea ice edge). Alternatively, we may have underestimated the contributions by benthic Nd fluxes in the Southern Ocean. The underestimation of upper NSW in the South Atlantic was already featured in the PI control (Fig. 3d), but is substantially more pronounced in simulation LGM_18Sv. This is the consequence of two aspects acting together. First, particle fluxes that transport the radiogenic Nd of AAIW to the underlying NSW by reversible scavenging are higher because of the elevated dust flux and reduced biogenic particle remineralization. These also increase Nd burial and in combination with the downsized Nd weathering flux into the North Atlantic cause Nd concentrations to be reduced Atlanticwide (Fig. 9b), which in turn gives the downward Nd transport by reversible scavenging more leverage for overprinting upper NSW. This process is further amplified by the second aspect, which is the strongly elevated Nd dust flux to the South Atlantic raising the εNd signature of AAIW and thus also the isotopic difference to NSW. The enhanced and more impactful downward transport of Nd by reversible scavenging also markedly affects the deep South Atlantic, where the radiogenic signature of SSW is more effectively overprinted. In addition to this, SSW also becomes more sensitive to overprinting by the benthic flux in these LGM simulations as the difference between bottom water and detrital εNd and thus its leverage is larger. Together, both processes increase the nonconservative behavior of Nd isotopes despite comparable advection rates of LGM_18Sv and the PI control.
Turning to the sluggish LGM_09Sv circulation regime, we find the effect of underestimating upper NSW to be much smaller, largely because the expansion of AAIW and thus the source of radiogenic Nd is greatly reduced. Yet, fundamental changes in the biogeochemical cycling that arise from such sluggish overturning contribute as well to the better representation of true water mass mixing by Nd isotopes in the upper water column, while having the opposite effect in the deep. The weakened circulation strongly limits North Atlantic nutrient upwelling and hence export production. Consequently, particle concentrations and Nd scavenging are critically reduced there. Along with the restricted advective Nd export, this leads to an accumulation of Nd in the North Atlantic with Nd concentrations of NSW exceeding 35 pmol/kg, which is more than a doubling compared to LGM_18Sv (Fig. 9b). This effect is essentially limited to the Atlantic (Fig. 9d), because circulation changes are largest there (as by design of the experiments). Moreover, circulation-induced changes in export production are much smaller in the Southern Ocean, and substantial parts of the Pacific since these regions are rather iron than nutrient limited today (Behrenfeld and Kolber, 1999;Parekh et al., 2008), and iron is readily supplied by the large dust flux during the LGM. Therefore, NSW retains its εNd signature, because non-conservatively changing its Nd isotope signature requires an amount of Nd proportional to its high Nd concentration, which is simply not available from the Nd-depleted AAIW. On the other hand, only a small fraction of the strongly Nd-enriched NSW needs to be transported downward by reversible scavenging to leave a large imprint in the underlying SSW, which itself exhibits a Nd concentration virtually independent of the circulation regime (Fig. 9a). As such, SSW is not only notably overprinted by the benthic flux, its impact is also enhanced due to the longer benthic exposure time, but also by reversible scavenging. Hence, it emerges that the interplay of these two processes, benthic exposure and reversible scavenging, attenuates the vertical Nd isotope gradients induced by the water masses to such a degree that the reconstruction of North Atlantic water mass provenance becomes challenging under sluggish deep circulation. Here, it is important to note that this deep ocean overprinting is rather independent of the AMOC strength (that is defined at intermediate depth) as indicated by simulation LGM_16Sv that also displays large overprinting in the deep ocean. This is the result of the deep Atlantic ocean residence time being largely determined by the water mass provenance and not the advection rate ( Figure S8), because the difference in ventilation age of SSW and NSW is much larger than any changes induced by different deep water formation rates in the North Atlantic (Muglia and Schmittner, 2021). Hence, vertical Nd isotope gradients in the North Atlantic are rather difficult to distinguish between states of strong and weak AMOC (Fig. 8b), somewhat limiting the use of εNd as a faithful water mass tracer in this region. This effect does not appear to be as strong in the South Atlantic, where a considerable, yet still attenuated, vertical εNd gradient develops (>2 ε-units; Fig. 8a). This general trend of decreasing vertical εNd gradients under more sluggish circulation states, particularly in the North Atlantic, also appears to be a robust result in the Bern3D model given that we find no meaningful change in this behavior in a set of sensitivity tests for which we varied the scavenging efficiency of POM and the dust source flux ( Figure S9). Only when changing the parameterization of the benthic flux isotopic composition to weight the detrital and bottom water signatures from a 2:1 ratio (see section 2.1) to a 1:2 ratio (i.e., assuming that 67% of the benthic Nd flux is "recycled" from authigenic phases) the overprinting in the deep Atlantic slightly decreases ( Figure S10, only~30% overprinting relative to the~40% in Fig. 7).
In sum, the here performed simulations therefore serve to emphasize that non-conservative behavior of εNd can result from nonlinear combinations of partially opposing, and spatially variable, contributing factors.

Deglacial evolution under transient boundary conditions
The last deglaciation was punctuated by a sequence of abrupt climate changes that were tightly coupled to the evolution of the AMOC ( Fig. 10; Clark et al., 2012;McManus et al., 2004). Here, we simulated this transition from the LGM to the present day with a focus on the rapid major weakenings and resumptions of the AMOC and their transient impacts on the Nd isotope evolution. These simulations are continued from the two most distinct simulations, LGM_18Sv and LGM_09Sv, in order to investigate whether the markedly different LGM circulation states prime a distinct Nd isotope response during the early deglaciation. Four exemplary reconstructed εNd time series are chosen here for a direct modeldata comparison (Fig. 10e), which represent the general trends of the major ocean basins. These are situated in the North Atlantic (Zhao et al., 2019), the South Atlantic (Skinner et al., 2013), the South Pacific (Basak et al., 2018), and the North Pacific (Du et al., 2018) (see also red symbols of Fig. 1a), and essentially show a diminishing deglacial Nd isotope trend from the North Atlantic (DεNd (LGM-Holocene) ¼ 3.8) to the North Pacific (1.0 ε-units).
Over the simulated deglaciation all Nd source fluxes are slowly scaled back from the LGM to the present day values based on the benthic d 18 O evolution as described in section 2.3. Since benthic d 18 O mainly represents the evolution of the continental ice sheets, we also consider this a reasonable approximation for the weathering regime and thus for the Nd source fluxes. In addition, millennial-scale AMOC variability is induced by North Atlantic freshwater fluxes ( Figure S11; Menviel et al., 2012). As a response to the latter, the AMOC comes to a nearly complete halt during late Heinrich Stadial 1 (HS1) and the Younger Dryas (Fig. 11a). While such a strong AMOC weakening is widely accepted for HS1 (e.g., Adkins, 2013;McManus et al., 2004), proxy reconstructions and model simulations are more ambiguous about the AMOC state during the Younger Dryas (Keigwin et al., 1991;McManus et al., 2004;Meissner, 2007). The simulations therefore rather serve to illustrate the εNd response to strong AMOC variability, while we acknowledge that the true evolution of the Atlantic overturning may have been somewhat different.
The general patterns of the simulated Nd isotope evolution at the four sites representing the major ocean basins are consistent with the reconstructed time series (Figs. 10e and 11e). Yet, absolute εNd values partly deviate, which is most pronounced for the North Pacific site where sedimentary εNd signatures are about 3 ε-units more radiogenic than simulated. This, however, is the result of an equally sized mismatch between sedimentary and bottom water Nd isotope signatures at this location (Du et al., 2018), which is also observed extensively throughout the entire Pacific (Hu et al., 2016). Despite this offset, the North Pacific εNd record was interpreted to represent bottom water Nd isotope variability by the original authors (Du et al., 2018). Therefore, three key findings emerge from these transient simulations.
The first finding is that the εNd evolution at all four sites appears to be only weakly dependent on whether the transient simulations are continued from LGM_18Sv or LGM_09Sv. As mentioned in section 4, the meridional εNd gradient of the Atlantic is smaller under a strong circulation regime, since unradiogenic NSW penetrates farther southward. In contrast, a sluggish AMOC with expansive SSW in fact leads to less radiogenic Nd isotope signatures in the North Atlantic as was also reported for the CESM Earth system model (Gu et al., 2019). These differences in the Atlantic meridional εNd gradient disappear during mid-HS1, when the AMOC strength reaches its minimum. At the same time, Nd concentrations and their early deglacial evolution at these locations indeed show substantial differences. This is most pronounced in the North Atlantic, where under sluggish circulation conditions the Nd concentration is exceedingly high (section 4), but even further increases until the AMOC comes to a full halt and then abruptly decreases again, even though the circulation continues to be at its minimum strength. The variability in Nd concentration is even greater for the simulation that exhibits a strong AMOC and consequently a low North Atlantic Nd concentration at the LGM. During the HS1 AMOC slow-down the concentration doubles within about 1 kyr and then equally fast falls back to its initial value. These changes in Nd concentration are inversely correlated to the POM export production in the North Atlantic, which itself rather appears to be a function of the change in AMOC strength than to be dependent on its absolute value, due to the nature of the nutrient cycle (Mariotti et al., 2012;Schmittner, 2005). After HS1 the differences in Nd concentrations become minor and are negligible at the time of the Younger Dryas.
The second key finding is, that the Nd isotope response to the pronounced deglacial AMOC variability is only minor at these four locations. This is mechanistically closely related to the first finding that already indicated limited differences in the εNd signatures linked to on the initial LGM AMOC state. During the deglaciation this becomes even more evident as Nd isotope responses at the four sites are only on the order of 0.5 ε-units despite nearly full cessations of the AMOC. We described the processes responsible for this insensitivity in detail in section 4 (less reversible scavenging leading to higher Nd concentrations and greater impact of the benthic flux), but here it emerges that they can in fact also critically impact the marine Nd cycling on millennial timescales. As such, the non-conservative behavior of Nd can quickly become dominant in large parts of the Atlantic when the AMOC slows down substantially. Yet, in the deep South Atlantic the Nd isotope time series shows more pronounced short-term excursions to unradiogenic signatures during the reinvigorations of the AMOC after HS1 and the Younger Dryas (Fig. 11e). These are caused by the flushing of the  (Andersen et al., 2004). (b) 231 Pa/ 230 Th from the Bermuda Rise (34 N, 58 W,~4.6 km; Lippold et al., 2019;McManus et al., 2004). (c) Mean ocean temperatures (MOT, orange line) derived from noble gas measurements of EPICA Dome C (Bereiter et al., 2018). (d) Global benthic d 18 O stack of Lisiecki and Stern (2016). (e) Nd isotope records from the North Pacific (site 87JC/U1418, Du et al., 2018), the South Pacific (site PS75/073, Basak et al., 2018), the South Atlantic (site MD07-3076Q, Skinner et al., 2013), and the North Atlantic (site KNR198 GGC35, Zhao et al., 2019). Locations of these sites are also marked as red symbols in Fig. 1a North Atlantic that had accumulated highly unradiogenic Nd during the previous AMOC shut-downs. However, because of their very brief nature (100e200 years) they are unlikely to be preserved in the sediments due to bioturbation ( Figure S12) and even less likely to appear in sedimentary records due to the usually limited temporal resolutions of these archives. In contrast to the εNd evolution at the four exemplary sites, millennial scale variability is indeed pronounced at the lower mid-latitudes (Fig. 12). Particularly in the South Atlantic, Nd isotope signatures vary by about 3 ε-units at 2e4 km water depth over the deglaciation, yet HS1 can only be resolved if the preceding LGM circulation was vigorous (i.e., LGM_18Sv). Previously, Howe et al. (2018) inferred the evolution of water mass mixing based on Nd isotopes from the southern Brazil Margin (sites between 20 S and 30 S) at 1e2 km water depth. Yet, here we find that this depth range is relatively insensitive to changes in water mass mixing (albeit very short excursions appear after the Younger Dryas that are however unlikely to be resolvable in sediment cores) as εNd signatures vary between À9 and À11 throughout the past 20 kyr and independently of the LGM circulation state. P€ oppelmeier et al. (2020a) expanded the dataset by Howe et al. (2018) to greater water depths (2.8 km), but also found no distinct millennial scale variability in Nd isotopes at the southern Brazil Margin. Combining these datasets with our simulations suggests first, that the LGM exhibited a relatively weak circulation so that the transition to HS1 was not particularly pronounced (or instead both, the LGM and HS1, exhibited a moderately weak circulation), and second, that during the Younger Dryas the AMOC experienced only a minor weakening, which did not lead to a strong A third important finding is that the general glacial-interglacial transition of Nd isotopes is almost exclusively driven by changes in Nd source fluxes and not by major reorganizations of the water mass geometry. In order to demonstrate this further, we performed two additional transient simulations that are identical to the previous ones except that all Nd source fluxes were held constant at the LGM level (Figures S15 and S16). In these two simulations, starting again from LGM_18Sv and LGM_09Sv, Nd isotope trends are only punctuated by the small millennial scale variability during the deglaciation, that are further attenuated compared to the previous transient simulations, but else remain globally nearly constant throughout the entire past 20 kyr. As such, we confirm the hypothesis of P€ oppelmeier et al. (2019, 2020b) and Zhao et al. (2019) that changes in the weathering regime have been the main drivers of the glacial-interglacial εNd variability. While the previous observation-based studies strongly focused on the changes in Nd weathering inputs associated with the large northern hemispheric ice sheets that undoubtedly dominated the global εNd evolution, we further emphasize here the important role of changes in the dust fluxes on surface and intermediate depth Nd isotope signatures. During the LGM these introduced a large amount of radiogenic Nd into the South Atlantic leaving a substantial imprint on the intermediate depth water masses and thus also affecting the underlying NSW via reversible scavenging.
In summary, our simulations of the last deglaciation suggest that the Nd cycle was substantially affected by dynamically changing input fluxes and non-conservative processes that can mute the changes in water mass mixing and therefore require a more careful interpretation of Nd isotope reconstructions, in particular in the North Atlantic. Yet, while we consider these processes to be robust features of the Nd cycle, we note that the extent and spatial expansion are dependent on the parameterization of a number of processes as well as the model resolution. Hence, it still remains to be determined how substantial the impact of these nonconservative effects truly were. Moreover, here we are not able to simulate a range of observations, such as a highly unradiogenic glacial deep Labrador Sea (Blaser et al., 2020) and Nordic Seas (Struve et al., 2019), or highly radiogenic signatures in the Northeast (Roberts and Piotrowski, 2015) and South Atlantic (Huang et al., 2020;Skinner et al., 2013), which either hint to missing processes or lacking resolution and general shortcomings of the model.

Conclusions
Here, we extended the implementation of Nd isotopes in the Bern3D model by revising a number of critical parameterizations, which provides us with an improved description of the marine Nd cycle. The coherent framework of the model, which includes tools to diagnose the true water mass distribution, allows us to assess the processes driving non-conservative Nd isotope behavior. Our findings highlight that the Nd isotopic composition of water masses is not only overprinted by the benthic flux, which was the main focus of previous investigations (Abbott et al., 2015;Du et al., 2020), but also by the vertical downward transport of Nd by reversible scavenging (Siddall et al., 2008).
In the modern Atlantic, both processes only have a limited impact, because benthic exposure times and deep ocean particle concentrations are low, and isotopic gradients between water masses as well as the differences to the sedimentary detritus are relatively small. These findings are consistent with global assessments of non-conservative behavior of modern seawater Nd isotopes (Du et al., 2020;Tachikawa et al., 2017) and therefore corroborate that the major processes influencing Nd cycling are adequately implemented in the Bern3D model.
During the LGM, the dramatic changes of the climatic boundary conditions were not only responsible for a substantial shift in the apportionment of Nd input fluxes, but also impacted many biogeochemical processes with important implications for the marine Nd cycle. We explored how these changes in Nd input fluxes and biogeochemistry and their complex interplay with the ocean circulation impacted the capacity of Nd isotopes to faithfully trace water mass mixing during glacial periods.
Our simulations indicate that non-conservative processes generally played a larger role under LGM boundary conditions, but their impacts are nonlinear, partially opposing, spatially variable, and distinctly dependent on the circulation state. In the Atlantic, this enhanced non-conservative behavior of εNd was the result of the interplay of a number of processes (see Fig. 13): First, the changed Nd input fluxes increased the isotopic differences between the Nd reservoirs (different water masses and sediments), which increased the effectiveness of the respective non-conservative exchange processes, i.e., reversible scavenging and benthic fluxes; second, as a result of the adjusted POM remineralization profile due to lower microbial activity under colder temperatures, the redistribution of biogenic particles in the water column from top to bottom gave rise to a stronger response of export production and hence particle concentrations to changes in circulation state, which in turn had a profound effect on scavenging and consequently on the Nd concentrations that nonlinearly increased in the North Atlantic with decreasing AMOC strength; finally, the increase of benthic exposure times with more sluggish bottom water circulation caused the benthic Nd flux to more strongly overprint the bottom water Nd isotope signatures. As a consequence, we find that εNd-based calculations of the Atlantic-wide NSW fraction lead to a large under-and overestimation at the intermediate and deep South Atlantic, respectively, under strong and deep glacial AMOC conditions. In contrast, a very sluggish and shallow AMOC causes the underestimation of the NSW fraction at shallow depths to be reduced while the εNd signature of deep SSW is heavily overprinted by sedimentary contributions along its northward pathway leading to a considerable overestimation of NSW in the deep North Atlantic. The vertical εNd gradient in the North Atlantic therefore exhibits only relatively small differences as a function of the true water mass configuration of the LGM derived from dye tracers, which renders the interpretation of Nd isotope records from this region more challenging than previously assumed. These findings are further corroborated by our transient simulations of the last deglaciation, for which we find a muted response in North Atlantic εNd to the sequence of rapid major weakenings and resumptions of the AMOC. In contrast, this simulated millennial-scale variability is indeed reflected most pronouncedly in the South Atlantic where the Nd isotope proxy is a more reliable water mass tracer than in the North Atlantic.
The ramifications of these findings for the general interpretative Fig. 13. Conceptual scheme of the interplay between ocean circulation, biogeochemistry, and Nd cycling and their impact on the Nd isotope signatures at mid-depth and in the deep ocean that is in contact with the seafloor. Blue and red represent positive and negative correlations, respectively, with biogeochemical and Nd cycling that are discussed in the text.
(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) framework of Nd isotopes are substantial, yet, we have to note that our model also has some shortcomings, such as low spatial resolution and simplifications in the biogeochemistry and Nd implementation, that may have biased the Results presented here.
Further, due to model limitations we have not fully explored all AMOC configurations previously proposed, including a very strong and shallow overturning circulation (Kurahashi-Nakamura et al., 2017;Lippold et al., 2012) or multiple distinct northern-sourced water masses sandwiching SSW (Keigwin and Swift, 2017). As such, the non-conservative effects described here may differ under consideration of a more complete description of the biogeochemical and Nd cycle or under other potential alternative circulation regimes. Hence, this study serves to illustrate that important assumptions commonly made for the interpretation of εNd records in the North Atlantic most likely do not hold due to the complex interactions of biogeochemistry, ocean circulation, and Nd cycling. By further investigating the processes that govern the Nd cycle and in particular the nature and spatial variability of the benthic flux, a more comprehensive understanding can be achieved, which will bolster our ability to predict and thus reliably correct for nonconservative processes in the framework of integrated modeldata assessments.

Data availability
Results of the simulations presented will be uploaded to the ZENODO repository (https://zenodo.org/) upon acceptance of the manuscript and are available upon request from the corresponding author.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.