Inclusion of a suite of weathering tracers in the cGENIE Earth System Model-muffin release v.0.9.10

The metals strontium (Sr), lithium (Li), osmium (Os) and calcium (Ca) and their isotopes are important tracers in the study of changes in weathering rates and volcanism, two main processes which shape the long-term cycling of carbon and other biogeochemically important elements at the Earth’s surface. Traditionally, isotopic shifts of these four elements in the geologic record are interpreted with isotope-mixing, tracer-specific box models because of their long residence times in the ocean. However, these models often lack mechanistic links between the cycling of the four metals and other geochemically 5 relevant elements, particularly carbon. Here we develop and evaluate the implementation of Sr, Li, Os and Ca isotopes into the Earth system model cGENIE. The model has the potential to study these metal systems at equilibrium and under perturbations alongside other biogeochemical cycles. We provide examples of how to apply this new model to investigate Sr, Li, Os and Ca isotope dynamics and responses to environmental change.

1 Introduction 10 The evolution of life and climate on Earth is intrinsically linked to the dynamics of carbon, oxygen and nutrients at the Earth's surface. While complex interactions on a range of spatial and temporal scales determine the distribution of these elements between surficial (non-lithological) reservoirs like oceans, atmosphere, biomass and soils, the total amount and isotopic distribution of these elements are ultimately governed by the balance between crust weathering, marine sediment deposition and mantle inputs during oceanic crust formation. Geological records show that the abundance of these elements at Earth's surface 15 varied on time scales of thousands to millions of years. For example, the changing composition of mantle fluxes into the ocean and nutrient supply during weathering of large mountain ranges were potentially drivers of the accumulation of surficial oxygen during the Archean and Proterozoic (e.g. Kump and Barley, 2007;Campbell and Allen, 2008). Likewise, the supply of partic-ulate iron from airborne grounded rock to the oceans of the Quaternary is invoked as a mechanism to drive glacial-interglacial carbon cycle dynamics (e.g Martin, 1990;Martínez-Garcia et al., 2011;Loveley et al., 2017;Hooper et al., 2019). The emplacement of large igneous provinces (LIPs) repeatedly led to climatic and biotic crises by increasing the supply of mantle-derived carbon and nutrients to the oceans, e.g. during oceanic anoxic events (e.g. Erba et al., 2010;Jenkyns, 2010;Bottini et al., 2012;Monteiro et al., 2012;Percival et al., 2015). Increased basalt weathering in the Cenozoic could have led to a continuous decline 5 in CO 2 concentrations, contributing to the reconstructed long-term cooling trend (e.g. Kent and Muttoni, 2013;Mills et al., 2014). Mass fluxes between the lithosphere and Earth's dynamic surface are thus important drivers of environmental processes on the Earth's surface over geological time scales, but our ability to measure these fluxes directly is limited by their long evolution times and large spatial extents. Instead, records of variations in the content and isotopic composition of trace metals like strontium (Sr), osmium (Os), lithium (Li) and calcium (Ca), preserved in marine sediments, are used to study lithological 10 processes. The dynamics of these metals are controlled by processes that also shape the long-term carbon and nutrient cycles, namely continental and oceanic weathering and direct mantle emissions, but variations in their isotopic compositions can be linked more directly to changes in these processes since radioactive decay and mass-dependent fractionation form isotopically distinct lithological metal reservoirs. The residence time of these metals in the ocean is generally assumed to be greater than the mixing time of the ocean and therefore their distribution is comparably homogeneous (e.g. Faure and Mensing, 2005). 15 While laboratory methods to measure and reconstruct the evolving trace metal composition of seawater have become increasingly precise, our ability to interpret these records is limited by our incomplete mechanistic understanding of the cycling of these metals and their isotopes. Quantifying reservoir sizes and fluxes constitute particular challenges under varying chemical conditions and marine productivity over geological time. Mass balances and box models have been most commonly used to interpret geological records of trace metal variations (e.g. Tejada et al., 2009;Misra and Froelich, 2012;Kristall et al., 2017;20 Them et al., 2017). They however rely on assumed fluxes and residence times based on today's ocean which might not always be applicable throughout the geological record.
Mechanistic Earth system models are useful tools to test our conceptual and mathematical understanding of processes against observations, and to explore characteristics and processes of biogeochemical cycling beyond modern observational data. The Earth system model of intermediate complexity cGENIE has been used successfully to study marine biogeochemical cycles 25 (including C and Ca) under various boundary conditions and external forcings (e.g. Ridgwell and Zeebe, 2005;John et al., 2014;Death et al., 2014;Turner and Ridgwell, 2016;Hülse et al., 2019). Implementing isotope enabled cycling of Sr, Os, Li and Ca in cGENIE allows us to test the influence of different assumptions about source and sink mechanisms, and efficiencies, on metal distributions in the ocean and their behaviour under environmental perturbations.
Here, we present our implementation of the marine cycling of Sr, Os, Li and Ca in cGENIE. We evaluate the model's 30 performance by comparing simulations of pre-industrial trace metal distributions in the oceans to seawater measurements, and study their equilibration times to imposed geochemical perturbations compared to estimated seawater residence times for these elements.
We first provide a brief review of the systematics of Sr, Os, Li and Ca cycling in the ocean, as well as their observed concentrations and isotopic compositions. This discussion provides the theoretical basis for our model implementation, and the data compilations for the model evaluation at the end of this study.

5
Of the four metals modelled in this study, Sr has the second highest marine concentration (3 -4 times more abundant than Li and more than one billion times more abundant than Os, Broecker and Peng, 1982). The main source for marine Sr is continental weathering, replenishing the ocean reservoir through rivers and potentially groundwater discharge (Basu et al., 2001;Beck et al., 2013). Smaller sources are hydrothermal input of mantle-derived Sr and refluxes from diagenetic alteration of carbonates at the sea floor. The chemical and physical similarity of Sr to Ca leads to its substitution in aragonite and, to 10 a lesser degree, calcite (Fietzke and Eisenhauer, 2006;Rüggeberg et al., 2008;Böhm et al., 2012;Stevenson et al., 2014).  argue that biological uptake is the dominant driver of spatial gradients in Sr concentrations in today's oceans, and Krabbenhöft et al. (2010) predict that marine carbonate burial could be the most important sink of seawater Sr. The concentration of Sr incorporated into biogenic carbonates depends on the mineralogy and the growth rate with higher Sr/Ca ratios in aragonite precipitated at fast growth rates (Rickaby et al., 2002;Stevenson et al., 2014). By contrast, direct effects of 15 temperature on elemental fractionation during the formation of calcite, the major carbonate buried in marine sediments, are small (Tang et al., 2008). Assuming rivers to be the only supply of continental Sr to the oceans, Hodell et al. (1989) estimates the marine residence time of Sr to be 1.9 -3.45 Myr. Considering the potentially appreciable Sr influx from groundwater, the actual residence time might be at or below the lower end of this estimate.
Elemental fractionation during magmatic processes creates different Rb/Sr ratios that over time generate reservoirs with distinct Sr isotopic signatures with the continental crust being considerably more radiogenic than the mantle (Faure and Mensing, 25 2005). Dust and rainwater isotopic signatures, both radiogenic and stable, are generally lower than modern day seawater values (Pearce et al., 2015). Volcanic material brings unradiogenic particulate Sr into sediments and can over time significantly affect the radiogenic composition of dissolved Sr in the sediment column while having little effect on the isotopic composition of dissolved Sr in seawater (Elderfield and Gieskes, 1982;Pearce et al., 2015). Changes in isotope abundances caused by kinetic and mass-dependent fractionation during carbonate formation are much smaller, and one can thus assume that the 87 Sr/ 86 Sr 30 of seawater and deposited carbonates is only controlled by inputs to the ocean (Krabbenhöft, 2011). By contrast, variations in δ 88/86 Sr are only controlled by fractionation during chemical processes, biogenic carbonate formation being the most important one in the ocean (e.g. Krabbenhöft et al., 2010). There is no evidence for a generic dependence of δ 88/86 Sr fractionation during biogenic carbonate formation on environmental conditions, although temperature and growth rate effects were observed in some calcifying species (Fietzke and Eisenhauer, 2006;Böhm et al., 2012;Stevenson et al., 2014;Vollstaedt et al., 2014).

Osmium
Osmium is a siderophile and chalcophile element which led to its accumulation in the Earth's core (Goldschmidt, 1922) and is also compatible during mantel melting. As a result, it is one of the rarest elements in the Earth's crust and the oceans. Because 5 of its low concentrations, the ability to measure Os in seawater only developed in the past few decades, limiting the number of observations of marine Os concentrations and isotopic compositions.
The cycling of Os at the Earth's surface is similar to that of radiogenic Sr in many aspects. Oxidative weathering of sediments is currently the most important natural source of Os to the oceans (Peucker-Ehrenbrink and Ravizza, 2000;Lu et al., 2017).
The estimated annual input of such Os to the oceans via aeolian and riverine transport and groundwater discharge is about 5 10 times larger than all other natural inputs combined (Lu et al., 2017). Hydrothermal inputs constitute the next biggest source of oceanic Os with high and low temperature systems being of roughly equal importance (Georg et al., 2013). The high temperature hydrothermal source can be split into a basaltic and a peridotitic source (Burton et al., 2010). Unlike for Sr, cosmic and terrestrial dust are significant sources of marine Os, as well as volcanic aerosols, though field observations do not suggest that the latter significantly contributes to the oceanic inventory (Sharma et al., 2007). Os is removed from the ocean by deposition at the seafloor, predominantly under suboxic water masses and in ferro-manganese nodules (Lu et al., 2017). Os incorporation into biogenic carbonates constitutes an additional minor sink (Burton et al., 2010). The marine residence time of Os is inherently uncertain given the uncertainty in Os fluxes, although estimates range from 3 -50 kyr (Sharma et al., 1997;Levasseur et al., 1998;Oxburgh, 2001 (Faure and Mensing, 2005). The isotopic composition of Os in marine sediments provides insight into changing fluxes between different Os reservoirs, particularly weathering related Os fluxes from the continents and mantle derived fluxes from volcanic activity. During partial melt in the upper mantle Os is compatible 10 and rhenium (Re) is incompatible, leading to an increased Re/Os ratio in continental crust relative to the mantle (Dąbek and Halas, 2007), and therefore elevated 187 Os/ 188 Os ratios in the continental crust compared to the mantle. No mass-dependent fractionation has been observed during Os uptake into macroalgae (Racionero-Gómez et al., 2017) and is not assumed to occur during other fluxes between surface reservoirs. Lithium is the 25th most abundant element in the Earth's crust, and it is more concentrated in continental crust than in oceanic crust (Baskaran, 2011). Weathering, in particular of silicate rocks, releases dissolved Li to rivers and soils where it is partially bound during clay formation (Kısakűrek et al., 2005;Dellinger et al., 2015;Pogge von Strandmann et al., 2017). The remaining dissolved Li is transported to the oceans via rivers and groundwater. Aeolian transport only contributes a minor flux of Li to the oceans (Baskaran, 2011). Li is also added to the oceans by hydrothermal vents and submarine weathering, but the size of 20 this flux is more uncertain (Chan et al., 1992;Hathorne and James, 2006). Removal from the ocean happens predominantly via Li adsorption onto clay minerals, with a minor proportion buried with Li-containing biogenic calcite (Hathorne and James, 2006). These removal fluxes are dependent on the abundance of inorganic carbon and pH (Hall and Chan, 2004;Marriott et al., 2004), and are not uniform (Hathorne and James, 2006). The residence time of Li in the ocean is estimated to be 0.3 -3 Ma (Stoffyn-Egli and Mackenzie, 1984) with more recent estimates closer to 1 Ma (Vigier and Goddéris, 2015).
Li has two stable isotopes (7.52% 6 Li and 92.48% 7 Li, Penniston-Dorland et al., 2017). Changes in the abundance of Li isotopes on Earth are related to mass-dependent fractionation processes. The isotopic composition of Li is expressed as δ 7 Li, which is the ‰ deviation of the 7 Li/ 6 Li ratio from the L-SVEC standard (δ 7 Li = ( earlier studies suggest considerable variation in seawater δ 7 Li (Carignan et al., 2004), recent studies suggest that seawater is remarkably homogeneous in its δ 7 Li (Hall et al., 2005;Rosner et al., 2007;Penniston-Dorland et al., 2017) consistent with its long residence time. Mass-dependent isotope fractionation has been observed during clay formation, adsorption onto minerals and incorporation of Li into calcite shells (Pistiner and Henderson, 2003;Rudnick et al., 2004;Dellinger et al., 2015;Hindshaw et al., 2019), while kinetic fractionation occurs in magmatic systems due to diffusion (Parkinson et al., 2007;Penniston-10 Dorland et al., 2017), and potentially in aqueous solutions (Richter et al., 2006). Isotopic differences in weathered lithologies and post-weathering formation of secondary minerals in rivers and soils result in a large spatial and temporal variability of the δ 7 Li of continental run-off (Huh et al., 1998;Pistiner and Henderson, 2003;Pogge von Strandmann and Henderson, 2015).
Temporal variations in the amount of secondary mineral formation on land have the potential to drive shifts in seawater δ 7 Li over geological time (Misra and Froelich, 2012;Pogge von Strandmann and Henderson, 2015). Fractionation during biogenic 15 carbonate formation results in carbonate δ 7 Li values that are a few ‰ lower than seawater, but this offset seems to be carbonate producer-dependent Hathorne and James (2006).

Calcium
Calcium cycling is closely linked to the C cycle, both shaping and shaped by the size of C reservoirs and C fluxes at the Earth's surface. Similar to Sr, Os and Li, the dominant Ca source for today's oceans is weathering-derived dissolved and particulate 20 Ca in continental runoff. Input through hydrothermal vents near ocean ridges is on the order of 20% of the riverine flux (e.g. Milliman, 1993;DePaolo, 2004). Unlike Sr, Os, and Li, Ca plays an important role in many biological systems, predominantly as an electrolyte and building block for biogenic minerals (e.g. shells, exoskeletons, bones and teeth). In the ocean, Ca ions are 6 https://doi.org/10.5194/gmd-2020-233 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
incorporated into biogenic minerals (e.g. foraminiferal tests and calcareous nannoplankton), or form hydrogenetic or authigenic minerals (Fantle et al., 2020) if waters are highly saturated (Ω=a Ca 2+ · a CO 2− 3 /K sp ). The resulting minerals are dissolved in undersaturated conditions, or buried, compacted, and lithified. Carbonate formation creates the biggest long-term Ca and C sinks in today's oceans, and marine carbonate accumulation and dissolution constitute a significant buffer mechanism to stabilize marine and atmospheric pCO 2 during periods of enhanced exogenic C input. The residence time of Ca in seawater is 5 estimated to be 0.5-1.3 million years (Milliman, 1993;Sime et al., 2007;Griffith et al., 2008). -1) · 1000). We use the δ 44/40 Ca notation in this study, from now on shortened to δ 44 Ca, with NIST SRM-915a as standard. Dissolution of carbonates and silicates, as well as the precipitation of secondary silicate minerals controls the Ca isotopic composition in soil pore fluids, lakes and rivers (Farkaš et al., 2007;Tipper et al., 2008;Hindshaw et al., 2013;Fantle and Tipper, 2014;Kasemann et al., 15 2014;Perez-Fernandez et al., 2017). Both, biotic and abiotic precipitation of carbonates fractionate Ca isotopically, generating minerals with low δ 44 Ca values relative to aqueous Ca 2+ . The isotopic fractionation factor is most strongly a function of precipitation rate (and thus and solution chemistry, e.g. DePaolo, 2011) and is close to one (∆ ∼ 0) in the marine sedimentary column (see reviews in Blättler et al., 2012;Fantle and Tipper, 2014). In the ocean, species-dependent fractionation has been observed for several groups of calcifiers, with a small dependence on temperature (e.g. Nägler et al., 2000). Dissolved Sr, Os, Li and Ca are largely homogeneous in sea water, which is illustrated in figure 1 by sorting all seawater measurements -independent of location and depth -by their measured value. Measurements of a perfectly homogeneous seawater property would appear as a horizontal line in this plot, since the same value would be measured everywhere in the ocean. Most trace metal measurements have a difference of less than one standard deviation from the respective mean, suggesting very small spatial heterogeneities. In particular, most measured isotopic compositions are analytically indistinguishable from the standard deviation of the overall data population. Larger differences between lowest and highest measured metal concentrations indicate heterogeneity (potentially horizontal and/or vertical gradients) in metal abundances or are an artefact of the small number of sub-surface measurements. One possible example is Sr, which is reportedly less abundant in surface waters and waters in the 5 North Atlantic than in deeper waters and sites outside the North Atlantic (see fig. A1). Higher Li concentrations in the surface ocean compared with deeper waters are only reported in one study, , while no vertical gradients are apparent in Fabricand et al. (1967) and Hall. These differences suggest that the spread in Li concentrations in fig. 1 could be the result of analytical uncertainty rather than real gradients in dissolved Li. Os concentrations vary between sites, but only a few sub-surface measurements are available and there is considerable variation in Os concentration between measurement 10 techniques ( Gannoun and Burton, 2014). Ca profiles show no apparent spatial patterns.

Application of Sr, Os, Li and Ca isotopes as proxies of weathering and mantle emissions
The isotopic differences between Sr, Os, Li and Ca in the ocean and exogenic reservoirs allows these metals to be used as proxies for mass exchange between the surficial Earth system and the mantle, continental crust (all four) and extra-terrestrial material (Os). However, they record different geochemical pathways because they behave differently in water and have different 15 predominant host lithologies.
Sr and Os are used as proxies for the balance between the weathering of old continental crust and juvenile basalt or direct mantle emission (e.g. Hodell et al., 1990;Goddéris and François, 1995;Tejada et al., 2009;Finlay et al., 2010;Bottini et al., 2012;Dickson et al., 2015). The fact that carbonates are the biggest source of continental Sr but not of Os (which resides predominantly in shales and evaporites), means differences in these records can be used to infer changes in weathered lithology 20 (e.g. Peucker-Ehrenbrink et al., 1995). Li is orders of magnitude more abundant in silicates than in carbonates and thus is discussed as the most direct proxy for silicate weathering (Kısakűrek et al., 2005;Millot et al., 2010). . Ca isotopes have also been used to examine weathering processes in the geological record (Kasemann et al., 2005;Farkaš et al., 2007;Kasemann et al., 2008;Blättler et al., 2011;Holmden et al., 2012;Fantle and Tipper, 2014;Kasemann et al., 2014), but are also crucial proxies for the quantification of carbonate precipitation rates (Pogge von Strandmann et al., 2019). Because of its tight links with the marine carbonate system, Ca also serves as a constraint on CO 2 concentration. Furthermore, Ca isotopes have been considered as a potential temperature proxy (e.g. Nägler et al., 2000), but this application might be complicated by the strong control of aqueous chemistry on Ca fractionation (DePaolo, 2011). Likewise, changes in the Sr/Ca ratio in calcifiers reflects ecosystem shifts, carbonate mineralogy and calcification rate (Stoll and Schrag, 2001), and stable Sr isotopes in calcifiers correlate with growth rate, mostly influenced by temperature and pCO 2 (Stevenson et al., 2014;Müller et al., 2018). While each system individually gives valuable insight into Earth system 5 dynamics, in concert, these trace metal systems can provide information on feedbacks and event durations, as well as improve the accuracy of our reconstructions (e.g. Kasemann et al., 2008).

Model implementation
We implemented the marine cycling of Sr, Os and Li and Ca isotopes in the Earth system model 'cGENIE'. cGENIE is a modular model, comprising modules for ocean physics ('GOLDSTEIN' Edwards and Marsh, 2005), marine biogeochemistry 10 ('BIOGEM , continental weathering and run-off ('ROKGEM' Colbourn et al., 2013), sea-floor sediment formation ('SEDGEM' , atmospheric chemistry ('ATCHEM') and atmospheric energy balance ('EMBM . As a whole, cGENIE includes a 3D ocean combined with a 1D atmosphere which capture the cycling of carbon and a range of other elements relevant for biogeochemical studies of the ocean water column and at the sea-sediment and sea-air interfaces, as well as climate-sensitive continental run-off and explicit sedimentary carbonate burial. Hashed arrows indicate processes during which isotopic fractionation occurs in the model. and SEDGEM that are described in the following sections.

Continental weathering and run-off
Continental weathering is the main source for marine Sr, Os, Li and Ca today and changes in the weathered lithology or weath-5 ering intensity are invoked as major drivers of the evolution of trace metal concentrations and their isotopic compositions in seawater over time (e.g. Misra and Froelich, 2012;Kristall et al., 2017). ROKGEM, the weathering module of cGENIE, provides a framework for calculating climate-and CO 2 -dependent additions of Ca, Mg and alkalinity from carbonate weathering (following Berner, 1994) and silicate weathering (following Brady, 1991) to the ocean (see Colbourn et al., 2013, for a full description of the weathering module). The module calculates a drainage map based on a prescribed continental topography 10 and determines the coastal locations where this input is added to the ocean. Depending on the primary source rock, we tied 11 https://doi.org/10.5194/gmd-2020-233 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
the rate M of Sr, Os or Li ion delivery to the ocean to Ca and magnesium (Mg) input rates from weathering of carbonates (Ca CaCO3 ) and/or silicates (Ca CaSiO3 ), related by a constant ratio k (see table 5): ROKGEM allows for different CaSiO 3 weathering schemes, including one which separates the contribution from CaSiO 3 weathering into contributions from granite and basalt weathering ('Wallmann'). For this case, we imposed individual parame-5 ters for Sr weathering for each lithology. The model does not consider secondary mineral formation on land, which should be taken into account by setting the parameters for terrestrial Li input based on run-off data rather than source rock composition.
The amount and isotopic composition of weathering-derived Os is mostly dependent on organic matter content (Jaffe et al., 2002;Georg et al., 2013;Dubin and Peucker-Ehrenbrink, 2015). ROKGEM does not contain an explicit representation of organic matter weathering, and a differentiation between carbonate and silicate derived Os would not capture the lithology 10 dependence of weathering derived Os appropriately. Hence, we decided not to differentiate the source lithology of Os in the model, but to prescribe the net abundance and isotopic composition of Os in continental run-off.

Hydrothermal inputs
Mantle-derived metal input via hydrothermal vents is parameterised in the SEDGEM module. To represent the vast range of chemical conditions and reaction rates at the sediment-water interface in the ocean efficiently, three biogeochemically distinct CaCO 3 deposition and dissolution as well as elemental fluxes across the sediment-water interface are parameterised differently in each of these environments, using masks to distinguish the different environments. Hydrothermal inputs of elements into bottom waters occurs only in grid cells which represent 'deep sea'. We parameterised this flux such that a global total input 20 flux can be prescribed, which is distributed equally across 'deep sea' grid cells. There is no differentiation between seafloor areas with predominant high-or low-temperature hydrothermal venting in our implementation. Hence, we suggest to prescribe total hydrothermal metal fluxes with an average isotopic composition in cGENIE to account for both temperature regimes.

Aeolian inputs
To account for inputs at the air-sea interface (only assumed to be relevant for Os, as described in the scientific background section), we added the option to prescribe a uniform or spatially explicit field of annual input into the surface ocean (or technically anywhere into the water column). This implementation follows the same equations of flux forcings for other marine tracers (e.g. oxygen) and is calculated as part of the BIOGEM module.

Metal export in carbonates
Ca export from the surface ocean is linked to the export of particulate organic carbon (POC) by a constant CaCO3:POC ratio (rain ratio) (see , for a description of POC export simulation in cGENIE). We parameterise the incorporation of Os, Li and Sr into carbonates by scaling their export in carbonates to the Ca export. This can be done by setting either a constant Sr/Li/Os-to-Ca ratio (r M/Ca ) or a constant scaling factor α which links the local Sr/Li/Os-to-Ca ratio 10 in seawater to the Sr/Li/Os-to-Ca ratio in the precipitate (following e.g. Tang et al., 2008): The Sr/Li/Os-to-Ca ratio of the precipitate is then used to calculate the metal consumption during carbonate formation, as 15 well as metal release during carbonate dissolution. To reflect the different Sr/Li/Os-to-Ca ratios in pelagic compared to reef carbonates, a separate factor can be set for reef production. Reef CaCO 3 immediately contributes to surface sediments, while pelagically formed CaCO 3 sinks through the water-column, and might be dissolved before reaching the 'deep sea' sedimentwater interface as a result of local chemical equilibria. cGENIE also calculates Li export by abiotically precipitated carbonates.

Scavenging from the water-column
To test implications of Os burial with particulate organic carbon, we included an option for [O 2 ]-sensitive scavenging from the water column. In the surface ocean, organic carbon is produced as a function of the concentration of inorganic carbon, light, nutrient availability (PO 3− 4 in our set up), sea ice cover and temperature (Ridgwell, 2001). The produced organic carbon is split into dissolved and particulate organic carbon (DOC and POC) at an adjustable ratio. DOC is advected and diffused in 5 the ocean and remineralized depending on a prescribed lifetime, while POC is instantly exported to deeper water layers and remineralized following a decay function tuned to POC flux measurements (Ridgwell, 2001). O 2 is exchanged between the atmosphere and the ocean at the air-water interface, depending on the solubility of O 2 in water. In the ocean, O 2 is released during primary production in the surface ocean and consumed during remineralization of DOC and POC under oxic conditions.
Following the example of existing scavenging parameterization schemes in cGENIE, we implemented this process by scaling 10 the flux of scavenged Os (f scav ) to the local concentrations of Os and P OC: Since there is evidence that Os needs to be reduced in order to be burried, we include a switch to use the scavenging code only where [O 2 ] is below an adjustable threshold.

Clay formation
The model calculates Li burial during clay formation f Li clay locally at the sea-sediment interface, which is scaled to the concentrations of Li in the deepest grid box of the water column and the detrital flux into the sediments: Recrystallization of SrCO 3 is an additional source of Sr to the water column (see above). Its parameterisation is similar to that of hydrothermal inputs, except that elemental fluxes from recrystallization only occur in bottom waters of grid cells labelled as 'reef'. The size and rate of this Sr flux can be either set as a total global or an area-weighted value.

Isotopes
The isotopic composition of Sr, Os, Li and Ca sources to the ocean is prescribed, in delta notation for stable isotopes and ratios for radiogenic isotopes. Once in the ocean, cGENIE tracks isotopes in total molar abundances rather than delta notation or ratios so that they can be advected and diffused like any other tracer. For this approach, prescribed isotopic compositions of marine inputs are converted into molar isotope abundances which requires a standard for each isotope. The hard-coded standards are 5 listed in table 12, and are either international standards or observations of modern-day sea water. Li isotopes are implemented similarly to stable C isotopes (Ridgwell, 2001). Sr, Os and Ca are different, because they have more than two principal stable isotopes. We chose to reduce the number of traced isotopes to three for Sr and Os and two for Ca, because these subsets contain the most abundant isotopes of the respective element (Sr) and/or are most relevant for their application as seawater proxies (Sr, Os and Ca). Calcium isotopes are thus treated similarly to Li isotopes. To track three isotopes rather than two, we 10 track two isotopes explicitly and one implicitly with the bulk abundance of the trace metal. In the case of Sr, abundances of 87 Sr and 88 Sr are tracked explicitly. The abundance of 86 Sr is then taken as the difference between the abundances of the two explicitly tracked isotopes and the bulk Sr abundance, since 84 Sr can be neglected. Os has more than two stable isotopes, but only two of them, 187 Os and 188 Os, are currently used as a proxy system and so these two isotopes are tracked explicitly and the bulk Os abundance contains all remaining stable isotopes. To account for isotopic fractionation during chemical reactions, 15 the amount f I of an isotope I transferred to the product-side of the reaction is calculated by multiplying the amount of Sr, Li or Ca involved in the reaction, f M , with the relative abundance k of I in the reaction product: k is derived from the prescribed relative abundance of I in the reaction product in standard sea water, usually the observed isotopic signature of the respective reaction product in today's ocean, and the difference between the local, simulated relative is still debated. Growth rate-dependent fractionation factors could be implemented in the future by using the ecosystem module ECOGEM (Ward et al., 2018), but this will significantly increase computational costs.

Pre-industrial configuration
We initialised cGENIE with a modern-day geography on a 36x36 grid with up to 8 depth layers in the ocean (following 10 Ridgwell and Hargreaves, 2007) and with a pre-industrial CO 2 concentration of 278 ppm. Marine biological productivity is simulated using a single nutrient scheme as outlined in  but with a constant CaCO 3 :POC rain ratio set to 0.043 (Panchuk et al., 2008). Burial and dissolution of CaCO 3 in deep-sea sediments follows .
This configuration results in a pelagic CaCO 3 burial rate of 0.125 GtC yr −1 , close to the observationally calibrated CaCO 3 sink of . In addition to the pelagic environment, covering 350.6 million km 2 , we simulated reef 15 deposition on shelfs, which covers 5.5 million km 2 in our set up. The sediment model bathymetry was derived from ETOPO5 1 .
Reefal deposition was simulated in grid cells representing marine sediments shallower than 1000 m and not further pole-wards than 41.8 • N/S, on the basis that reefal carbonate deposition is predominantly a tropical and subtropical process. We prescribe a reefal CaCO 3 sink of 0.04 GtC yr −1 . Temperature dependent terrestrial silicate and carbonate weathering (Lord et al., 2016) were included with the baseline temperature and rates of CaCO 3 and CaSiO 3 weathering set to 8.48 mid-ocean ridges and consumption of 1.5 PgC yr −1 with δ 13 C = 2 ‰ during seafloor weathering (Cocker et al., 1982). Table   14 summarizes the trace metal fluxes we prescribe to simulate the pre-industrial trace metal cycling. Note that the prescribed hydrothermal metal fluxes are net fluxes of high-and low-temperature hydrothermal activity.
By applying the chosen model parameters for our spin-up with modern boundary conditions and running the model to steady state, we assumed that metal fluxes are currently in equilibrium. This is a common assumption when modelling weathering 30 tracer isotopes and their perturbations in the geologic record (e.g. Misra and Froelich, 2012;Bauer et al., 2017). However, given the long residence time estimates in today's ocean, it is likely that at least the Sr, Li and Ca cycles are not fully in steady state today (i.e. Derry, 2009). Constant Os isotopic compositions of seawater over the last millennia suggest that Os has reached steady state since the last deglaciation, but residence time estimates of >20 kyr based on modern day Os fluxes put this into question (Oxburgh, 2001). The newly implemented tracers in cGENIE can be used to study different equilibria and transient adaptations of the marine metal reservoirs to new boundary conditions. As an example, we simulated two different steady states for Sr under pre-industrial boundary conditions: One with a best estimate of pre-industrial Sr fluxes (hereafter referred to as 5 FLUXES) and one tuned to best fit the spatial mean of observations (TUNED).
20 https://doi.org/10.5194/gmd-2020-233 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License. The spin up was done in three stages to improve computational efficiency. The first stage (20 kyr) was run with a closed marine carbonate system, where the marine C and alkalinity reservoirs are artificially restored by balancing losses through

Comparison between simulated and observed trace metal contents of seawater
One advantage of simulating trace metal cycles within a 3D earth system model is that we can test, for the first time, simulated metal concentrations and isotopes against observed spatial patterns. In the following sections, the simulated pre-industrial distribution of each metal is compared against observational data shown in section 2.5.  to be smaller in the model simulations than in observations. This is particularly true for Sr concentrations which are more homogeneous in the simulations than in observations. Measured isotopic compositions of Sr in seawater are very homogeneous, which is matched by the model simulations. The FLUXES simulation, in which we constrain fluxes rather than the marine 5 reservoir with observations, results in higher than observed Sr concentrations in seawater, and a more radiogenic and isotopically heavy composition. Uncertainties in size and isotopic composition of pre-industrial Sr fluxes are big enough to allow for realistic combinations that lead to equilibrated Sr concentrations close to observations (TUNED simulation) but continental input needs to be less radiogenic than the average of today's rivers to yield a steady state marine 87 Sr/ 86 Sr close to obser-23 https://doi.org/10.5194/gmd-2020-233 Preprint. Discussion started: 10 August 2020 c Author(s) 2020. CC BY 4.0 License.
vations (in agreement with e.g. Pearce et al., 2015). The observed fractionation factor for stable Sr isotopes during biogenic CaCO 3 formation is too big to equilibrate the model at the observed low δ 88 Sr, even if we assume that all Sr sources only provide mantle-like light Sr (δ 88 Sr = 0.256‰, TUNED simulation). Corroborating the conclusion of e.g. Vance et al. (2009), this suggests that the isotopic composition of marine dissolved Sr is currently not at equilibrium.
The simulations capture the homogeneous distribution of Sr concentrations observed at most site (see figures A1 and A2 in 5 the SI). However modelled Sr concentrations are higher than observed in some studies of North Atlantic water. Given that the low North Atlantic values come from one of the first studies measuring Sr concentrations in seawater , and that other studies of North Atlantic seawater yielded higher Sr concentrations (e.g. De Villiers, 1999), this difference could be partially the result of analytical errors or seasonal/inter-annual variability.  Woodhouse et al., 1999;Gannoun and Burton, 2014). Since the depletion of Os in Pacific surface waters is a robust feature against measurement methods, this indicates that 15 some processes affecting Os distribution in that region might be missing in our model. For example, Woodhouse et al. (1999) suggested that Os could bind onto organic matter in the water column under low oxygen conditions. cGENIE provides the necessary framework to investigate this and other assumptions, which should be the subject of a separate study. Our aim here is to focus on reproducing the basin-scale distribution in trace metals. Observations from other sites show no vertical Os concentration gradients, which is well-reproduced by our simulation set-up. The isotopic Os composition shows a uniform vertical distribution in the simulations, which is consistent with observations at all sites (see figures A3b). Given relatively large analytical uncertainty, it is more likely that this results from measurement inconsistencies rather than unknown processes affecting the isotopic composition of dissolved Li.

Osmium
The high concentrations which are underrepresented in our simulation were mostly measured in the upper water column 10 (see fig. A4). As shown in fig. A1, these measurements all come from one study in the Caribbean Sea , a marginal ocean which is not representative of the open ocean and which can barely be resolved in the coarse grid of cGENIE. Additionally,  was one of the first studies measuring Li concentrations in seawater, and their measurements may be associated with much larger uncertainties than reported. Were these measurements excluded, the spatial mean of observed Li concentrations would be lower, and the model could be re-tuned by prescribing slightly lower Li 15 fluxes in-and out of the ocean which are permissible given the current range of flux estimates (see table 3).
Only one measured vertical profile of Li isotopes is available in the literature (Hall), and it shows no vertical δ 7 Li variation which is reproduced by our simulation. The simulated 7 Li is offset from the profile reported by Hall since we tuned the model to the average of reported seawater 7 Li measurements but the profile in Hall contains some of the most 7 Li enriched published seawater values. Simulated concentrations of dissolved Ca and its homogeneous isotopic composition closely resemble the observations. Similar to the average Ca concentration in seawater, the uniform vertical distribution is well captured by the model simulation (see fig.   A5). Likewise, no vertical Sr/Ca gradients are observed or simulated. Observed Sr/Ca ratios are in between simulation results with tuned and un-tuned marine Sr reservoir. This is an artefact of the tuning to mean observed Sr concentrations which is lower than the mean Sr concentrations of the sites shown in figure A5 because it includes the unusually depleted samples from 10 the North Atlantic. This mismatch is another indication that the Sr measurements reported in  might be associated with a larger uncertainty than reported.

Transient perturbation experiments
Another advantage of the new metal cycles in cGENIE is that we can study their transient behaviour under external forcings in a fully coupled system, including effects from ocean circulation and primary productivity changes, climate-sensitive weathering 15 and marine carbonate accumulation and dissolution. Here we present one example, where we instantly release either 1000 Pg C or 5000 Pg C into the atmosphere to produce temporary weathering responses which re-equilibrate the model (similar to Lord et al., 2016). These scenarios allow us to showcase the sensitivity of the simulated metal systems and to discuss the Earth system processes that shape their transient evolution. Since our experiment are set up to understand model behaviour, it is of less concern here whether these scenarios cause isotopic perturbations beyond the analytical uncertainty of real-world measurements. Os being the most responsive system. spheric CO 2 and temperature to an exogenic CO 2 pulse. After the instantaneous carbon release, atmospheric CO 2 concentration increases and the δ 13 C of seawater decreases on the timescale of years. The sudden increase in atmospheric CO 2 concentrations replenishes the marine carbon reservoir, which decreases biogenic carbonate formation in the surface ocean and leads to the dissolution of pelagic carbonates. Over the next ∼10,000 yrs, this negative carbon cycle feedback removes the majority of 5 the initial excess carbon (Lord et al., 2016). Increased atmospheric CO 2 also enhances the greenhouse effect and thus leads to a rise in global mean air temperature, altering the hydrological cycle and speeding up reaction rates which result in increased rates of continental crust weathering. After a few ten thousands of years, increased continental weathering supplies enough alkalinity to the ocean to slow down carbonate dissolution and increase carbonate burial. Through enhanced silicate weathering and carbonate burial, the remaining excess carbon is sequestered over the timescale of ∼100,000 yrs after the emission. 10 Over the same timescale, temperature, which due to Earth system feedbacks initially recovers more slowly than atmospheric CO 2 concentrations, returns to pre-emissions values and ends the phase of enhanced weathering. Metal input from continental weathering (and marine carbonate dissolution for Sr and Ca) leads to a transient growth of marine metal reservoirs and isotopic excursions. Increased continental weathering increases the influx of radiogenic Sr and Os, driving the respective marine reservoirs to more radiogenic isotopic compositions. Continental crust is also depleted in heavy stable isotopes compared to 15 seawater, and the increased weathering influx therefore reduces the δ 88 Sr, δ 7 Li and δ 44 Ca values of seawater. The isotopic composition of carbonates and silicates in continental crust is fixed in our set-up, but since carbonate and silicate weathering rates vary differently with climate, the isotopic composition of weathering-derived marine input varies for metals with different isotopic compositions in these lithologies. For example, the radiogenic Sr excursion in the ocean is amplified by a transient increase in 87 Sr/ 86 Sr of weathering-derived Sr. Since more alkalinity than Ca 2+ is derived from continental weathering in our 20 set-up, the climate-driven increase in weathering results in a transient imbalance between alkalinity and Ca 2+ supplies which does not appear in C cycle simulations without an igneous Ca 2+ source at the seafloor (e.g. Lord et al., 2016;Vervoort et al., 2019). The imbalance is small (10% of the Ca input at its peak) but cumulatively the unbalanced alkalinity supply is large enough to causes a slight reduction of the steady-state marine Ca and Sr reservoirs (by ∼0.5% following the emission of 5000 Pg C). Our assumption of a climate-insensitive Ca 2+ source at the seafloor is a simplification since there is evidence that this 25 flux depends on bottom water properties, in particular temperature (e.g. Krissansen-Totton and Catling, 2017). However, the climate sensitivity of this Ca 2+ flux is lower than that from continental weathering (Brady and Gíslason, 1997) and temperature change is dampened in the deep ocean. In our transient simulations this would mean that Ca 2+ supply from the seafloor would increase less than from the continent, were this climate sensitivity to be included in the model. Without this additional feedback, our set-up thus produces an upper limit for a transient imbalance between alkalinity and Ca 2+ supplies which is the 30 consequence of decoupled alkalinity and Ca 2+ sources with different climate sensitivities.
In the presented simulations, only the radiogenic 87 Sr/ 86 Sr excursion amplitude is larger than measurement uncertainty in present-day seawater (see figures 3 -6). Stronger or prolonged climate change, or a more sensitive weathering regime, would be required to produce detectable weathering changes in the other isotope systems. However, our simulation results still illustrate the different responses of all five isotope systems to a sudden carbon injection. The timings of excursion peaks and minima differ strongly between the isotope systems. The radiogenic 187 Os/ 188 Os and negative δ 88 Sr and δ 44 Ca excursions peak a few tens of thousands of years after the C injection, while the radiogenic 87 Sr/ 86 Sr and negative δ 7 Li excursions take hundreds of thousands of years to reach their full amplitude. These time lags result from different residence times and fractionation processes: With the shortest residence time, estimated both from present-day inventories and in our simulation (see table 15), Os is the 5 first metal to re-equilibrate after the perturbation. The excursions in Os concentration and isotopes are only driven by enhanced weathering rates, and thus both recover within 300 kyrs once atmospheric temperature, and with it continental weathering rates, decrease sufficiently. The radiogenic 87 Sr/ 86 Sr excursion, however, continues to grow until weathering rates have fully returned to pre-event levels because concentration-dependent Sr sinks adapt more slowly to enhanced continental Sr input as a result of the longer residence time. The re-equilibration of the marine Sr reservoir is substantially slower than that of the Os reservoir 10 for the same reason. However, the removal of marine Sr is also dependent on carbonate chemistry, since marine carbonate preservation constitutes the largest marine Sr sink. Biogenic carbonate preservation initially decreases because of the pH drop following the C injection, but it increases during the C cycle recovery as a result of enhanced weathering-related alkalinity input.

Conclusions
Sr, Os, Li and Ca isotope records can help identify mantle activity and lithological responses to climate change in the geological record if processes governing their distribution in the ocean and their response to geological perturbations were well 30 understood. Our implementation of the marine cycling of Sr, Os, Li and Ca in the Earth system model cGENIE allows us to investigate these processes and their consistency with observations. Simulating pre-industrial Sr, Os, Li and Ca distributions in the ocean, the model achieves slightly more homogeneous fields than observed, but consistent with the widely used assumption of homogeneous metal isotope distributions in a fully equilibrated state (e.g. Elderfield, 1986;Chan and Edmond, 1988). It reproduces the mostly homogeneous mean observed depth profiles of concentrations and isotopic compositions well, but partially deviates from local observations. This could be an artefact of measurement uncertainties, non-equilibrated cycles in today's ocean or local processes not resolved or parameterised in cGENIE. We showed how the new cGENIE tracers thus provide an opportunity to investigate the consistency between reported metal compositions in the ocean and marine sources and sinks, 5 and between different ocean basins. Further investigations of the differences between cGENIE simulations and observations have thus the potential to improve our understanding of present-day marine trace metal cycling. Furthermore, the new implementation can be used to study the response of these metal cycles to transient environmental change. As an example, we show how these systems react to a sudden CO 2 release through perturbation to their input and output fluxes on different timescales.
cGENIE thus becomes a valuable tool to improve our understanding of Sr, Os, Li and Ca isotopes and distributions in today's 10 ocean, and to quantitatively interrogate the environmental signals of these isotope systems as preserved in sedimentary records.
Code availability. The model source code, as well as the simulation configuration files, are publicly available and can be accessed on github.
8 Figure A2. Comparison between measured and simulated Sr profiles in seawater. Shown are composites of all available measured Sr concentration and isotope ratio profiles and all profiles in the simulation. Data is taken from , Fabricand et al. (1967), Bernat et al. (1972, Brass andTurekian (1974), De Villiers (1999) and Mokadem et al. (2015).
9 Figure A3. Comparison between measured and simulated Os vertical profiles in seawater. Shown are composites of all available measured Os concentration and isotope ratio profiles and all profiles in the simulation. Data is taken from Levasseur et al. (1998), Woodhouse et al. (1999) and Gannoun and Burton (2014). Figure A4. Comparison between measured and simulated Li profiles in seawater. Shown are composites of all available measured Li concentration and isotope ratio profiles and all profiles in the simulation. Data is taken from , Fabricand et al. (1967), Chan (1987), Chan and Edmond (1988), You and Chan (1996), Moriguti and Nakamura (1998), Tomascak et al. (1999), James and Palmer (2000), Košler et al. (2001), Nishio and Nakai (2002), Hall, Bryant et al. (2003), Pistiner and Henderson (2003)  Implications for the osmium geochemical cycle and the use of osmium as a paleoceanographic tracer, Geochimica et Cosmochimica Acta,