Inclusion of MyAMI-derived Mg ∕ Ca corrections to the marine carbonate system in the cGENIE.cookie Earth system model (v.0.91)

The concentrations of the major cations (esp., Ca²⁺, Mg²⁺) in Earth's oceans have undergone large-scale fluctuations in the geological past. This is important because the key geochemical properties of the marine environment that underpin the global carbon cycle – the aqueous carbonate system equilibria and solubility of solid calcium carbonate (CaCO₃) – are heavily influenced by ion-pairing, which in turn depends on the activity of the major cations and anions. An appropriate interpretation of marine proxies as well as the reconstruction of past states of ocean geochemistry and carbon cycle across geologic events requires that these effects are considered. However, most current global carbon cycle models use empirical carbonate system dissociation constants fitted to laboratory experiments with present-day seawater major cation and anion concentrations and in the few simulations of global carbon cycling that have accounted for paleo-seawater composition, only relatively simplified empirical adjustments of the equilibrium constants (from Ben-Yaakov and Goldhaber, 1973; Tyrrell and Zeebe, 2004) have been implemented (e.g., Panchuk et al., 2008).

Here we develop and evaluate a new scheme in the cGENIE Earth system model for correcting carbonate system equilibrium constants and accounting for variations in the dissolved calcium and magnesium concentrations in the ocean. We base our new parameterization on the MyAMI specific ion interaction model of Hain et al. (2015) and implement this in cGENIE by means of linear interpolation within a 4-dimensional parameter look-up table of pre-calculated carbonate system equilibrium constant values. For modern seawater composition, our implementation of MyAMI-based equilibrium constants yields no meaningful deviation from model results using empirically-based equilibrium constants, validating our look-up/interpolation approach. However, for simulations conducted under non-modern $\mathrm{Mg}/\mathrm{Ca}$, we find substantial differences in carbon chemistry and CaCO₃ saturation when using our new MyAMI-based equilibrium constants as compared to the existing (default) correction scheme in cGENIE. Specifically, our new MyAMI-based correction scheme exhibits a much lower sensitivity to an instantaneous change in $\mathrm{Mg}/\mathrm{Ca}$ from modern to Eocene and an approximate doubling of Ca²⁺ concentration, in both surface ocean pH (0.01) and calcite saturation state (9.41), which were overestimated by 0.05 and 4.04, respectively, with the previous correction scheme. Any bias in carbonate chemistry and CaCO₃ saturation will also affect the preservation and burial of CaCO₃ in deep-sea sediments and hence potentially impact model-data comparisons. We illustrate this by contrasting the ocean carbon inventory arising under Eocene $\mathrm{Mg}/\mathrm{Ca}$ with the same total weathering (and hence CaCO₃ burial) flux for the different possible equilibrium constant corrections. We find that the new MyAMI-based and previous default corrections give rise to a dissolved inorganic ocean carbon inventory 20 µmol kg⁻¹ (348 PgC) higher and 59 µmol kg⁻¹ (950 PgC) lower, respectively, relative to the same experiment conducted using empirical equilibrium constants without any $\mathrm{Mg}/\mathrm{Ca}$ correction. Applying no correction at all for a different-from-modern $\mathrm{Mg}/\mathrm{Ca}$ ratio in the ocean would appear to be better than applying an overly-approximated correction but explicitly accounting for past dissolved calcium concentrations remains of fundamental importance. We provide this new carbonate system equilibria correction as an option in cGENIE.muffin version 0.9.76, and as standard in a forthcoming new cGENIE code release – cGENIE.cookie v.0.91.

1 Introduction

The major ion (Na⁺, K⁺, Ca²⁺, Mg²⁺, Cl⁻, SO${}_{\mathrm{4}}^{\mathrm{2}-}$) composition of the ocean has changed substantially through Earth history (e.g., Horita et al., 1991; Hardie, 1996, Fig. 1). For example, the Ca²⁺ concentration ([Ca²⁺]) was higher during the Cretaceous (up to ∼ 33 mmol kg⁻¹ compared to 10 mmol kg⁻¹ today) while [Mg²⁺] was lower (35–40 mmol kg⁻¹ compared to 52 mmol kg⁻¹ today, Timofeeff et al., 2006; Brennan et al., 2013; Hain et al., 2015). The first order consequence of higher Cretaceous [Ca²⁺] (as we will illustrate later) is that assuming a relatively similar global weathering (and hence carbonate burial) rate to the present day, the ocean carbonate ion (CO${}_{\mathrm{3}}^{\mathrm{2}-}$) concentration would have been lower. This is a simple and direct consequence of carbonate saturation involving the product of [Ca²⁺] and [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] – for the same saturation state (the primary control on marine carbonate burial), if one goes up, the other must go down (Ridgwell, 2005; Ridgwell and Zeebe, 2005). This trade-off between higher Cretaceous [Ca²⁺] and lower [CO${}_{\mathrm{3}}^{\mathrm{2}}$] at relatively invariant saturation then allows for elevated atmospheric CO₂ and reduced ocean surface pH (e.g., Henehan et al., 2019) to be reconciled with an abundance of geological carbonate deposits at that time (Hönisch et al., 2012). At the same time, lower [CO${}_{\mathrm{3}}^{\mathrm{2}}$] also creates a seawater carbonate system that would have been less buffered against a rise in CO₂ (Hain et al., 2015), i.e. CO₂ absorption leads to larger pH changes. Changing proportions of [Mg²⁺] relative to [Ca²⁺] (hereafter: $\mathrm{Mg}/\mathrm{Ca}$) has also been linked to major changes in the dominant mineralogy of carbonate-precipitating organisms (Hardie, 1996; Stanley and Hardie, 1998; Stanley et al., 2005; Stanley, 2006), with calcite precipitation favoured in low-[Mg²⁺] seawater (such as during the Cretaceous) but calcite precipitation inhibited in favour of the aragonite polymorph in higher $\mathrm{Mg}/\mathrm{Ca}$ modern seawater (e.g., Berner, 1975; Morse et al., 2007).

Figure 1 Changes of seawater [Mg²⁺], [Ca²⁺], [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] and $\mathrm{Mg}/\mathrm{Ca}$ since the Cretaceous as reconstructed from halite inclusions (Timofeeff et al., 2006; Brennan et al., 2013).

[Ca²⁺] and [Mg²⁺] have additional but much less appreciated impacts on the aqueous equilibria of carbon: both contribute to the ionic strength of seawater and both interact with carbonate species and borate (Millero and Thurmond, 1983; Harvie et al., 1984) as a result of their ability to form complexes with carbonate, bicarbonate and hydroxide ion (Larson et al., 1973; Harvie et al., 1984; Pitzer, 1991; Stefánsson et al., 2017). As a result, changing [Ca²⁺] and [Mg²⁺] substantially affects the activity coefficients of free anions, the proportion of free and complexed anions, as well as the value of the equilibrium constants determining carbon speciation, CO₂ fugacity (i.e. the effective pressure of CO₂, accounting for molecular interactions), and the CaCO₃ saturation state of seawater. While environmental (temperature, salinity, pressure) controls on the carbonate system equilibrium are accounted for in empirically-determined modern seawater equilibrium constants, equilibrium constants for paleo-seawater also need to be adjusted to account for different [Ca²⁺] and [Mg²⁺] (Hain et al., 2015). The necessary corrections can be computed using Pitzer-type models based on potentiometric data (Pitzer, 1991; Millero and Pierrot, 1998; Hain et al., 2015; Clegg et al., 2023). This is relevant in reconstructions of atmospheric CO₂ from B isotope-based seawater pH, when solving the carbonate system in Earth System models, and especially when using reconstructed pH to validate simulated past global carbon cycle states and carbon-climate dynamics occurring on geological timescales.

The “muffin” release of the Earth system model of intermediate complexity “cGENIE” (Ridgwell et al., 2007) has been widely used to generate realizations of marine carbon cycling and atmospheric CO₂ for a variety of geological intervals and events as well as for carrying out model-data comparison against carbon cycle-related paleoceanographic proxies (e.g., Gutjahr et al., 2017; Greene et al., 2019). The-atmosphere-ocean-sediment biogeochemistry originally developed in cGENIE was calibrated against spatial observations of modern ocean geochemistry (Ridgwell et al., 2007) and surface sediment composition (Ridgwell and Hargreaves, 2007) and based on carbonate chemistry understanding rooted in modern seawater composition (Ridgwell, 2001). However, in the first paleo carbon cycle application of directly contrasting simulated and observed deep-sea sedimentary calcium carbonate (CaCO₃) contents across the Paleocene-Eocene Thermal Maximum warming event (“PETM”, ca. 55 Ma) (Panchuk et al., 2008), the potential importance of higher early Eocene ocean [Ca²⁺] (18.2 mmol kg⁻¹) and lower [Mg²⁺] (29.9 mmol kg⁻¹) in modifying carbonate preservation and burial in marine sediments was recognized. Following Tyrrell and Zeebe (2004), two adjustments were made to the aqueous carbonate chemistry scheme of Ridgwell et al. (2007) (described in the Supplement of Panchuk et al., 2008 but reproduced here for completeness):

• 1.

  Firstly, the solubility coefficient for calcite was adjusted as a function of the deviation of ambient $\mathrm{Mg}/\mathrm{Ca}$ from a modern reference ratio (Tyrrell and Zeebe, 2004).

  $$\begin{array}{}\text{(1)}& K_{\mathrm{sp}}^{\ast }=K_{\mathrm{sp},\mathrm{modern}}^{\ast }-\mathit{\alpha }\cdot \left({\displaystyle \frac{{\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]}_{\mathrm{modern}}-\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]}{{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}_{\mathrm{modern}}-\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}}\right)\end{array}$$

  where $K_{\mathrm{sp}}^{*}$ is the apparent (see below) solubility coefficient of calcite for a modern seawater composition (Mucci, 1983), [Mg²⁺]_modern and [Ca²⁺]_modern are the modern mean seawater concentrations (52.82 and 10.25 mmol kg⁻¹, respectively), and $\mathit{\alpha }=\mathrm{3.655}\times {\mathrm{10}}^{-\mathrm{8}}$ is a scaling constant.

• 2.

  Secondly, because of the tendency of Mg²⁺ to form ion pairs in seawater (Zeebe and Wolf-Gladrow, 2001), the 1st and 2nd apparent equilibrium constants of carbonic acid ($K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$) were scaled as a function of the relative deviation of [Mg²⁺] from modern, following Ben-Yaakov and Goldhaber (1973).

  $$\begin{array}{}\text{(2)}& \begin{array}{rl}{K}_{\left(x\right)}^{\ast }& =\left(\mathrm{1}+s_{K^{\ast }\left(x\right)}\cdot {\displaystyle \frac{\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]-{\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]}_{\mathrm{modern}}}{{\left[{\mathrm{Mg}}^{\mathrm{2}+}\right]}_{\mathrm{modern}}}}\right)\\ & \cdot K_{\left(x\right),\mathrm{modern}}^{\ast }\end{array}\end{array}$$

  where $K_{\left(x\right),\mathrm{modern}}^{*}$ is the apparent equilibrium constant for modern seawater composition (in cGENIE – Mehrbach et al., 1973 as refit by Dickson and Millero, 1987), and s_K(x) is a sensitivity parameter. x is either 1 or 2 – corresponding to $K_{\mathrm{1}}^{*}$ or $K_{\mathrm{2}}^{*}$. The values of s_K(x) are 0.155 and 0.442 for $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$, respectively (Ben-Yaakov and Goldhaber, 1973).

Note that all K^*s in this manuscript are the apparent dissociation constants. In seawater, ion activities differ from ion concentrations due to ion interactions (Millero, 1979). Thermodynamic dissociation constants are defined using ion activities and thus explicitly account for ion interactions, while apparent dissociation constants incorporate those ion interaction effects so that the dissociation equilibria can be expressed with ion concentrations. Apparent dissociation constants therefore vary with seawater composition.

The carbonate system corrections for [Ca²⁺] and [Mg²⁺] have also been employed in several subsequent cGENIE.muffin paleo studies (e.g., Ridgwell and Schmidt, 2010; Kirtland Turner and Ridgwell, 2013; Jennions et al., 2015; Gutjahr et al., 2017; Greene et al., 2019; Henehan et al., 2019). However, the availability of recently developed seawater chemical speciation models now leads us to re-visit the existing scheme in cGENIE, and implement and fully evaluate an updated correction scheme for use in paleo studies.

The current state-of-the-art in seawater chemical speciation modelling combines the classical ion-pairing approach (Sillen, 1961; Garrels and Thompson, 1962; Morel and Morgan, 1972; Millero and Schreiber, 1982) with specific ion interaction modelling of the free ions (Pitzer, 1991). This approach was implemented as the “MIAMI” spreadsheet model and evaluated against empirical measurements of standard modern seawater by Millero and Pierrot (1998). MIAMI builds on a substantial empirical database of dissociation constants in artificial and, particularly for $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$, natural seawater (Millero and Pierrot, 1998; Millero et al., 2002). It is a comprehensive and often-used model to calculate chemical speciation specifically in seawater (Turner et al., 2016). A simplified version of MIAMI was re-implemented (MyAMI) in the Python language by Hain et al. (2015) to make the approach more accessible to paleo seawater studies and Earth System modelling. To reduce the computational effort required in MIAMI, MyAMI only contains Pitzer equations for the activities of the species most relevant for the carbonate system. Furthermore, these were shortened by removing the higher-order electrostatic terms that are more relevant for media with much more elevated ionic strength than open ocean seawater. MyAMI enables carbonate system equilibria to be solved for a wide range of temperatures, salinities, [Mg²⁺] and [Ca²⁺] (assuming the modern seawater relationship between salinity and ionic strength). MyAMI predicts the dissociation constants of the carbonate system in the present-day ocean with an error of just a few percent (Hain et al., 2015) and was, at the time calibrated against the best-practice modern seawater empirical equilibrium constants (Chap. 5, Sect. 7 in Dickson et al., 2007). Subsequently, Zeebe and Tyrrell (2018) recognized an inaccuracy in the K₁ predicted by MyAMI, and this issue was resolved by replacing the literature source for the calcium-bicarbonate interaction parameters (Hain et al., 2018; https://github.com/MathisHain/MyAMI, last access: 26 July 2026). MyAMI-derived corrections to the temperature- and salinity-dependencies of the equilibration constants are available in Seacarbx (Raitzsch et al., 2022), cbsyst (Branson et al., 2023), and Kgen (Whiteford et al., 2025) and are commonly used to reconstruct local carbonate system parameters for the late Cretaceous through the Cenozoic (e.g., Sosdian et al., 2018; Anagnostou et al., 2020; Rae et al., 2021; CenCO2PIP, 2023).

Carbonate system solvers that include MyAMI compute three effects of a change in seawater major ion composition (here: [Ca²⁺] and [Mg²⁺]) on seawater carbon chemistry that would be expected be realized instantaneously in a closed ocean-atmosphere system (without any addition or removal of carbon or alkalinity from the system) : (1) the saturation (Ω) of CaCO₃ changes proportionally with seawater [Ca²⁺] (e.g., Ridgwell, 2005), (2) the equilibrium constants respond with modest changes due to weak specific ion-ion interactions, affecting the equilibrium carbonate speciation between CO₂, bicarbonate, and carbonate ion and shifting pH, and (3) the degree of carbonate ion complexation changes proportionally to the change in total of divalent cations, as formulated by Millero and Pierrot (1998) and implemented in MyAMI (Hain et al., 2015, 2018). However, they cannot calculate changes in CO₂ and Ω due to carbon cycle feedbacks (Hain et al., 2024) and the adjustment of global carbon and alkalinity inventories to re-balance weathering on land and CaCO₃ burial in the ocean. Accounting for this third effect of [Ca²⁺] and [Mg²⁺] change requires global carbon cycle or Earth system models.

In this paper, we describe and fully evaluate an implementation of MyAMI-derived carbonate system constants including the effect of different-from-modern [Ca²⁺] and [Mg²⁺] in the cGENIE Earth system model.

2 Methods

We start by providing a full description of the existing aqueous carbonate chemistry scheme in the cGENIE Earth system model. We then describe how the major equilibrium constants are derived from MyAMI and implemented as an alternative option in cGENIE. Finally, we describe our testing and evaluation methodology for the MyAMI-enabled marine carbon cycle in cGENIE and how this compares both to the current $\mathrm{Mg}/\mathrm{Ca}$ correction scheme as well as to not accounting for different-from-modern ocean $\mathrm{Mg}/\mathrm{Ca}$ ratios at all. Note that “muffin” (original) and forthcoming “cookie” releases of cGENIE differ only in the criteria for [H⁺] convergence (described below) and whether or not the MyAMI-based correction scheme is employed by default (cookie) or only as an option that needs to be specified (muffin).

2.1 The cGENIE Earth system model

Aqueous carbonate chemistry in cGENIE is governed by two tracers in the ocean circulation model – (1) dissolved inorganic carbon (DIC) which is the sum of: CO_2(aq) (ignoring the contribution from carbonic acid H₂CO₃), ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (bicarbonate ions), and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (carbonate ions), and (2) alkalinity (ALK) which is defined following Dickson (1981) (see Eq. 13, below).

The default empirical fits employed in cGENIE (muffin release) for the 1st and 2nd dissociation constants of carbonic acid (treating the activity of H₂O as effectively constant),

$$\begin{array}{}\text{(3)}& {\displaystyle }{K}_{\mathrm{1}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]\cdot \left[{\mathrm{CO}}_{\mathrm{2}}\right]}}\text{(4)}& {\displaystyle }{K}_{\mathrm{2}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}}\end{array}$$

are those of Mehrbach et al. (1973) as refitted by Dickson and Millero (1987) and are functions of the ambient environmental conditions of temperature (T) and salinity (S).

Strictly speaking, the empirically-determined equations for the various dissociation (and stability) constants are valid only over the range of experimental conditions from which they were derived. For $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$ from Mehrbach et al. (1973) (refit by Dickson and Millero, 1987) – the current default choice in cGENIE.muffin – T and S are limited to the range: $\mathrm{2}\le T\le \mathrm{35}$ °C, and $\mathrm{26}\le S\le \mathrm{43}$ PSU. Model-projected values of T and/or S lying outside of this range are simply truncated at the minimum or maximum empirical limits for the purpose of the carbonate chemistry calculation. For consistency, the same T and S range limitations as for $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$ are placed on all other dissociation (and carbonate stability) constants. In contrast, the default environmental limits associated with the Bunsen gas solubility coefficients are $\mathrm{2}\le T\le \mathrm{35}$ °C and $\mathrm{26}\le S\le \mathrm{43}$ PSU, while the temperature limits for calculating the Schmidt numbers used in air-sea gas transfer are $\mathrm{0}\le T\le \mathrm{30}$ °C, following Wanninkhof (1992). Implementing MyAMI-derived carbonate system constants allows the T and S ranges for $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$ to be expanded (although the Bunsen gas solubility coefficients and Schmidt numbers are left unchanged).

cGENIE solves for acidity on the seawater pH scale (pH_SWS), with all dissociated constants converted to pH_SWS if originally fitted on a different scale. To do so, the model first estimates the sulphate and flouride speciation using K^*s on the free scale, assuming the modern empirical relationship between salinity and total dissolved sulphate and flouride. These speciations are used to determine all the required conversion factors to the seawater scale. This includes the K^*s of sulphate and fluoride, but these adjusted, seawater-scale values are only used in the calculation of carbonate alkalinity. Numerical solution of the carbonate system is via an implicit iterative method to obtain the equilibrium hydrogen ion concentration ([H⁺]). Iteration n+1 of this calculation depends on the results of the previous iteration (n) via:

$$\begin{array}{}\text{(5)}& {\left[{\mathrm{H}}^{+}\right]}_{n+\mathrm{1}}={\left({\left[{\mathrm{H}}^{+}\right]}_n^{K_{\mathrm{1}}^{\mathrm{2}}}+{\left[{\mathrm{H}}^{+}\right]}_n^{K_{\mathrm{2}}^{\mathrm{2}}}\right)}^{\mathrm{0.5}}\end{array}$$

where:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\left[{\mathrm{H}}^{+}\right]}_n^{K_{\mathrm{1}}^{\ast }}=K_{\mathrm{1}}^{\ast }\cdot {\displaystyle \frac{{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_n}{{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}_n}}\text{(7)}& {\displaystyle }{\left[{\mathrm{H}}^{+}\right]}_n^{K_{\mathrm{2}}^{\ast }}=K_{\mathrm{2}}^{\ast }\cdot {\displaystyle \frac{{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}_n}{{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}_n}}\end{array}$$

In turn, the concentrations of: CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ in the nth iteration are estimated from Skirrow (1975):

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{\mathrm{CO}}_{\mathrm{2},n}& ={\mathrm{DIC}}_n-{\mathrm{ALK}}_{\mathrm{DIC},n}\\ & +{\displaystyle \frac{k\cdot \left({\mathrm{ALK}}_{\mathrm{DIC},n}-{\mathrm{DIC}}_n\right)-\mathrm{4.0}\cdot {\mathrm{ALK}}_{\mathrm{DIC},n}+z_n}{\mathrm{2.0}\cdot \left(k-\mathrm{4.0}\right)}}\end{array}\text{(9)}& {\displaystyle }{\mathrm{HCO}}_{\mathrm{3},n}^{-}={\displaystyle \frac{k\cdot {\mathrm{DIC}}_n-z_n}{k-\mathrm{4.0}}}\text{(10)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{3},n}^{\mathrm{2}-}={\displaystyle \frac{k\cdot \left({\mathrm{ALK}}_{\mathrm{DIC},n}-{\mathrm{DIC}}_n\right)-\mathrm{4.0}\cdot {\mathrm{ALK}}_{\mathrm{DIC},n}+z_n}{\mathrm{2.0}\cdot \left(k-\mathrm{4.0}\right)}}\end{array}$$

where:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{z}_n& ={\left(\left(\mathrm{4.0}-k\right)\cdot {\mathrm{ALK}}_{\mathrm{DIC},n}+k\cdot {\mathrm{DIC}}_n\right)}^{\mathrm{2}}\\ & +{\left(\mathrm{4}\cdot (k-\mathrm{4.0})\cdot ({\mathrm{ALK}}_{\mathrm{DIC},n}{)}^{\mathrm{2}}\right)}^{\mathrm{0.5}}\end{array}\end{array}$$

and

$$\begin{array}{}\text{(12)}& k={\displaystyle \frac{K_{\mathrm{1}}}{K_{\mathrm{2}}}}\end{array}$$

We define carbonate alkalinity, ALK_DIC, using the full alkalinity definition of Dickson (1990) (but excluding only the contribution from S²⁻, which is important only at very low values of pH; Zeebe and Wolf-Gladrow, 2001):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathrm{ALK}}_{\mathrm{DIC},n}& =\mathrm{ALK}-{\mathrm{H}}_{\mathrm{4}}{\mathrm{BO}}_{\mathrm{4}}^{-}-{\mathrm{OH}}^{-}-{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\\ & -\mathrm{2.0}\cdot {\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}-{\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}-{\mathrm{NH}}_{\mathrm{3}}-{\mathrm{HS}}^{-}\\ & +{\mathrm{H}}^{+}+{\mathrm{HSO}}_{\mathrm{4}}^{-}+\mathrm{HF}+{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\end{array}\end{array}$$

In calculating these components, we use the apparent ionization constant of boric acid ($K_{\mathrm{B}}^{*}$) from Dickson (1990);

$$\begin{array}{}\text{(14)}& K_{\mathrm{B}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{H}}_{\mathrm{4}}{\mathrm{BO}}_{\mathrm{4}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]\cdot \left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{BO}}_{\mathrm{3}}\right]}}\end{array}$$

converting from the original total pH scale to the seawater pH scale (Millero, 1995), and the apparent ionization constant of water ($K_{\mathrm{B}}^{*}$) from Millero (1992) (fitted on pH_SWS scale):

$$\begin{array}{}\text{(15)}& K_{\mathrm{H}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{OH}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]}}\end{array}$$

The 1st, 2nd, and 3rd dissociation constants of phosphoric acid, as well as the dissociation constants of silicic acid and ammonium, are all from Yao and Millero (1995) (and were all originally fitted on the pH_SWS scale).

$$\begin{array}{}\text{(16)}& {\displaystyle }{K}_{P_{\mathrm{1}}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\right]}}\text{(17)}& {\displaystyle }{K}_{P_{\mathrm{2}}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\right]}}\text{(18)}& {\displaystyle }{K}_{P_{\mathrm{3}}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}}\text{(19)}& {\displaystyle }{K}_{{\mathrm{NH}}_{\mathrm{4}}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{NH}}_{\mathrm{3}}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}}\text{(20)}& {\displaystyle }{K}_{\mathrm{Si}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{4}}{\mathrm{SiO}}_{\mathrm{4}}\right]}}\end{array}$$

The dissociation constant of hydrogen sulphide is from Hershey et al. (1988) (converting from total to seawater pH scale using estimated sulphur speciation and fluor speciation on the free scale);

$$\begin{array}{}\text{(21)}& K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{HS}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right]}}\end{array}$$

the dissociation constant of bisulfate is from Dickson (1990) (converting from free to seawater pH scale),

$$\begin{array}{}\text{(22)}& K_{{\mathrm{HSO}}_{\mathrm{4}}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}{\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]}}\end{array}$$

and for HF, we adopt the dissociation constant of hydrogen fluoride from Dickson and Riley (1979) and convert from the free to seawater pH scale:

$$\begin{array}{}\text{(23)}& K_{\mathrm{HF}}^{\ast }={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\cdot \left[{\mathrm{F}}^{-}\right]}{\left[\mathrm{HF}\right]}}\end{array}$$

In solving the carbonate system, the initial [H⁺] value ($[{\mathrm{H}}^{+}{]}_{n=\mathrm{1}}$) everywhere in the ocean is seeded at the start of a model experiment with an approximately representative bulk ocean value of 10^−7.8 (pH = 7.8). Thereafter, the initial [H⁺] value each time (and time-step) in which the carbonate system is solved is taken from the equilibrium value calculated at the previous time-step at the same ocean model grid point. In the original pH solution code (muffin), we judge the system to be sufficiently converged when the relative change in [H⁺] between iterations is below some threshold (which itself scales inversely with [H⁺]). Typically, this threshold equates to changes in [H⁺] of less than 0.01 % between iterations, equivalent to changes in pH and the fugacity of CO₂ (fCO₂) to within ±0.001 units (pH_SWS) and ±0.2 µatm, respectively (and compared to a fully converged solution. In cookie, we simplify this and adopt an explicit pH convergence threshold of 0.001 units (pH_SWS) (although this can be adjusted). Both variants of the same convergence scheme are generally stable for plausible (including future and much of geologically relevant) differences between DIC and ALK. Extreme ratios of DIC : ALK or ALK < ca. 500 µmol eq. kg⁻¹ (less than 25 % of modern) can lead to numerical instability if the two independent [H⁺] estimates (Eqs. 6 and 7) are sufficiently far apart, or one becomes negative.

By default, pH (including constants) is only recalculated every time-step in the ocean surface grid points in the model – a necessity for calculating air-sea gas exchange and saturation state. At the seafloor, if the sediment model (e.g., Ridgwell and Hargreaves, 2007) is used, pH is updated every sediment time-step (every year). Otherwise, to improve computational efficiency, pH is only solved in the ocean interior if required in providing model output, or if ocean interior carbonate speciation is required such as in the case of methanotrophy (Reinhard et al., 2020).

Dissolved Ca²⁺ and total S (HSO${}_{\mathrm{4}}^{-}$ + SO${}_{\mathrm{4}}^{\mathrm{2}-}$) are typically configured as prognostic tracers in cGENIE model experiments, in which case their oceanic distributions are simulated explicitly. If not selected, their concentrations are estimated from salinity following Millero, 1982, 1995):

$$\begin{array}{}\text{(24)}& {\displaystyle }\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]=\mathrm{0.01028}\cdot {\displaystyle \frac{S}{\mathrm{35.0}}}\text{(25)}& {\displaystyle }\left[{\mathrm{SO}}_{\mathrm{4}\left(\mathrm{tot}\right)}\right]=\mathrm{0.02824}\cdot {\displaystyle \frac{S}{\mathrm{35.0}}}\end{array}$$

Even if sulfate is carried as an explicit tracer in the model, for internally correcting stability constants between different pH scales, [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] is always derived using the modern relationship with salinity as above (Eq. 25). The justification for this is that the stability constants used in cGENIE, including the new MyAMI-derived ones, are all based on a modern seawater sulfate composition so as to avoid bias when calculating pH on the total scale. Note that the explicit tracer concentration of sulfate, if available, is used in the calculation of ALK_DIC (Eq. 13 above).

The concentrations of total boric acid and fluorine are typically not included as prognostic tracers in cGENIE model experiments, and are estimated from salinity following Millero (1982, 1995):

$$\begin{array}{}\text{(26)}& {\displaystyle }\left[{\mathrm{B}}_{\left(\mathrm{tot}\right)}\right]=\mathrm{0.000416}\cdot {\displaystyle \frac{S}{\mathrm{35.0}}}\text{(27)}& {\displaystyle }\left[{\mathrm{F}}_{\left(\mathrm{tot}\right)}\right]=\mathrm{0.00007}\cdot {\displaystyle \frac{S}{\mathrm{35.0}}}\end{array}$$

If a marine cycle of silica and hence the tracer silicic acid (H₄SiO₄) is not included in a given experiment configuration, a zero concentration is assumed throughout the ocean. For the modern ocean, the error in atmospheric CO₂ induced by this simplification compared to a carbon cycle utilizing observed H₄SiO₄ concentrations is <1 µatm (Ridgwell, 2001). Note that if the H₄SiO₄ tracer is included but used as a diagnostic tracer for tracking rates of silicate weathering and there is no corresponding opal sink (e.g., Hülse and Ridgwell, 2025), it can be omitted from the calculation of ALK_DIC so as to avoid unintended impacts on ocean pH.

Finally, the saturation index (or “saturation state”) with respect to CaCO₃, Ω, is defined as the product of calcium and carbonate ion concentrations divided by solubility product $K_{\mathrm{sp}}^{\ast }$, where seawater with Ω>1 is oversaturated and Ω<1 undersaturated with respect to CaCO₃:

$$\begin{array}{}\text{(28)}& \mathrm{\Omega }={\displaystyle \frac{\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{K_{\mathrm{sp}}^{\ast }}}\end{array}$$

The solubility products for calcite and aragonite, $K_{\mathrm{sp},\mathrm{cal}}^{*}$ and $K_{\mathrm{sp},\arg }^{*}$, respectively, are from Mucci (1983) and the corresponding saturation states are given the notation Ω_cal and Ω_arg, respectively.

All dissociation constants are corrected for pressure (P) following Millero (2005) and assuming that we can approximately relate depth below the ocean surface and pressure as 1 dbar m⁻¹ (which induces an error of no more than ∼ 3 % even at the deepest depths of the modern ocean, Ridgwell, 2001). Pressure corrections are applied to all dissociation constants used in cGENIE – $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$, $K_{\mathrm{B}}^{*}$, $K_{\mathrm{W}}^{*}$, $K_{\mathrm{Si}}^{*}$, $K_{\mathrm{HF}}^{*}$, $K_{{\mathrm{HSO}}_{\mathrm{4}}}^{*}$, $K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{*}$, $K_{{\mathrm{NH}}_{\mathrm{4}}}^{*}$, $K_{P_{\mathrm{1}}}^{*}$, $K_{P_{\mathrm{2}}}^{*}$, and $K_{P_{\mathrm{3}}}^{*}$, plus $K_{\mathrm{sp},\mathrm{cal}}^{*}$ and $K_{\mathrm{sp},\arg }^{*}$. Following Millero (1979), the general form of the pressure correction factor is:

$$\begin{array}{}\text{(29)}& \ln \left({\displaystyle \frac{K_{\left(P\right)}^{\ast }}{K_{\left(\mathrm{0}\right)}^{\ast }}}\right)=-\left({\displaystyle \frac{\mathrm{\Delta }V}{R\cdot T}}\right)\cdot P+\left({\displaystyle \frac{\mathrm{0.5}\cdot \mathrm{\Delta }\mathit{\kappa }}{R\cdot T}}\right)\cdot P^{\mathrm{2}}\end{array}$$

where P is the applied pressure in bars, T is the temperature (K), ΔV and Δκ are the molar volume and compressibility change for the dissociation reactions, respectively, and R is the gas constant (83.145 bar cm³ mol⁻¹ K⁻¹). For each dissociation reaction the values of ΔV and Δκ in seawater are approximated as function only of temperature (assuming a salinity of 35 PSU) by

$$\begin{array}{}\text{(30)}& {\displaystyle }\mathrm{\Delta }V=a_{\mathrm{0}}+a_{\mathrm{1}}\cdot T+a_{\mathrm{2}}\cdot T^{\mathrm{2}}\text{(31)}& {\displaystyle }{\mathrm{10}}^{\mathrm{3}}\cdot \mathrm{\Delta }\mathit{\kappa }=b_{\mathrm{0}}+b_{\mathrm{1}}\cdot T\end{array}$$

using the coefficients a₀…b₁. In the absence of an available relationship describing the effect of pressure on the dissociation of silicic acid, the same pressure effect as for boric acid is assumed.

The values of the various coefficients are corrected for historical typographical errors in the literature where necessary (see Lewis and Wallace, 1998 for an overview of some (but not all) of the typos prevalent in the literature, and Orr et al. (2015) for a more recent intercomparison of packages solving the aqueous carbonate system and corrected coefficient values). We detail all the dissociation constant coefficients used in cGENIE in Table S1 in the Supplement. Pressure corrections to the various dissociation constants are formulated based on Lewis and Wallace (1998) (also see Orr et al., 2015). Again, because of the occurrence of typographical errors in the literature, we also, for completeness, detail all the pressure correction coefficients used in cGENIE in Table S2.

2.2 Implementation of MyAMI-derived lookup tables in cGENIE

Since [Mg²⁺] and [Ca²⁺] are spatially and temporally variable in the ocean, a direct application of MyAMI in cGENIE would require running the MyAMI carbonate system solver for every grid cell and at every time step. This process would be excessively computational expensive, and hence we did not attempt to explicitly incorporate the MyAMI code (written in Python) itself into cGENIE (written in FORTRAN 77, F77, and Fortran 90, f90). Instead, we followed the approach of Ridgwell et al. (2003) (in that particular case, sediment dissolution fluxes were pre-calculated as a function of 5 different boundary conditions and substituted for a mechanistic 1D reaction-transport model for the seafloor flux) and generated look-up tables of carbonate system equilibrium constants under various boundary conditions, with the gridded values created through offline calculations with MyAMI. The gridded values are included in the cGENIE source code repository in the form of an ASCII format file for each equilibrium constant that is read in by the cGENIE model at runtime. The equilibrium constants are quadri-linearly interpolated from the look-up table using run-time environmental conditions (T, S, [Ca²⁺], [Mg²⁺]) in each ocean model grid-cell. We validate the accuracy of this implementation in Sect. 3.1 below.

We generated one look-up table for each of the equilibrium constants: $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$, $K_{\mathrm{W}}^{*}$, $K_{\mathrm{B}}^{*}$, $K_{{\mathrm{HSO}}_{\mathrm{4}}}^{*}$, $K_{\mathrm{sp},\arg }^{*}$, and $K_{\mathrm{sp},\mathrm{cal}}^{*}$. To create these values, the MyAMI model version 1.0 (Hain et al., 2015, 2018) was run for all combinations of seawater temperature (−2 to 50 °C in 1 °C steps), salinity (30–45 PSU in 1 PSU steps), [Ca²⁺] (1–60 mmol kg⁻¹ in 1 mmol⁻¹ steps) and [Mg²⁺] (1–60 mmol kg⁻¹ in 1 mmol kg⁻¹ steps). For environmental conditions at any location in the cGENIE ocean grid falling outside of these limits, no extrapolation is applied, and parameter values are capped at the limit value. Compared to the default Mehrbach et al. (1973) equilibrium constants, the permissible range in T and S in this new parameterization is now expanded to $-\mathrm{2}\le T\le \mathrm{50}$ °C, and $\mathrm{30}\le S\le \mathrm{45}$ PSU (instead of $\mathrm{2}\le T\le \mathrm{35}$ °C, and $\mathrm{26}\le S\le \mathrm{43}$ PSU). When MyAMI is selected in the cGENIE model, for consistency the permissible T and S ranges for the carbonate system solver need to be adjusted accordingly. The ranges for Bunsen gas solubility coefficients and Schmidt numbers are left at their defaults (see earlier).

To retain backwards compatibility, MyAMI carbonate dissociation constants are used in place of Mehrbach et al. (1973) if explicitly selected by the user in muffin but will become the defaults in cookie. When MyAMI is enabled, the original (muffin) carbonate dissociation constants, including the [Mg²⁺] and [Ca²⁺] corrections of Ben-Yaakov and Goldhaber (1973) and Tyrrell and Zeebe (2004), are simply replaced with ones derived from MyAMI. For compatibility with the pH units used in cGENIE, we convert the MyAMI-derived constants from the total to the seawater pH scale, and from molarities (mol per kg pure water solvent) to amount concentrations (mol per kg seawater). Further, we combine $K_{\mathrm{0}}^{*}$ and $K_{\mathrm{1}}^{*}$ from the MyAMI output as is done in the default carbonate system scheme in cGENIE to treat H₂CO₃ implicitly (Pines et al., 2016):

$$\begin{array}{}\text{(32)}& \log K_{\mathrm{cGENIE}}^{\ast }=\log K_{\mathrm{1},\mathrm{cGENIE}}^{\ast }+{\displaystyle \frac{\mathrm{1}}{\log K_{\mathrm{0},\mathrm{cGENIE}}^{\ast }}}\end{array}$$

The interpolated constants are then pressure-corrected consistent with other cGENIE carbonate constants with the pressure correction carried out via partial volumes and compressibility parameters taken from Millero (1979, 1983, 1995).

Note that we did not change [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] from modern in the current study. An expansion of the carbonate system solver to correct for varying [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] will be a focus of future model development.

2.3 Implementation of calcium carbonate diagenesis in cGENIE

The lookup table approach to calculating CaCO₃ dissolution in surface sediments of the deep ocean (and hence carbonate burial) (Ridgwell and Hargreaves, 2007), although able to account for changing ocean [Ca²⁺], is uncorrected for deviations in $\mathrm{Mg}/\mathrm{Ca}$ from modern (Ridgwell, 2001). For paleo seawater compositions, cGENIE.cookie explicitly employs the (F77 converted to f90) code of Archer (1991) in calculating the equilibrium rate of CaCO₃ dissolution (given mean wt % CaCO₃, bottom-water [O₂], organic carbon rain rate, seafloor depth, carbonate chemistry, etc.). Compared to Archer (1991), we propagate $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$, and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ from the main cGENIE carbonate chemistry solver, allowing us to hence also propagate a $\mathrm{Mg}/\mathrm{Ca}$ correction. The profiles of porosity and bioturbation rate in the original model are substituted with those of Ridgwell (2001) as described in Ridgwell (2007). Here (in comparison to e.g., Ridgwell and Hargreaves, 2007), because the model of Archer (1991) is an “oxic-only” approximation and omits the role of NO${}_{\mathrm{3}}^{-}$ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ reduction at the higher organic matter rain fluxes characteristic of continental margins, we limit the calculation of CaCO₃ dissolution (and hence burial) to depths greater than 1000 m and assume no carbonate preservation at shallower depths. Finally, the numerical solution of the 1-D reaction-transport equations for steady-state CaCO₃ dissolution (Archer, 1991) is not unconditionally stable. In the very few (<1 %) model sediment grid points where at some point in time a robust solution is not achieved, we substitute the explicit scheme with the lookup table of Ridgwell et al. (2003) and omit values for those particular grid points when making comparisons (Figs. 4, S5 and S8 in the Supplement).

2.4 Methodology for the evaluation of carbonate chemistry and steady-state marine carbon cycling in cGENIE.cookie

We perform two tests to evaluate the updated carbonate system (using the cGENIE.cookie codebase rather than the very slightly different pH convergence criteria of muffin). Firstly, we assess the accuracy of interpolating K^*s from our lookup tables rather than explicitly calculating them for each grid cell. We do this by simply contrasting the equilibrium carbonate constants derived internally in cGENIE by interpolating within the MyAMI-derived look-up tables vs. those explicitly calculated by MyAMI under pre-industrial [Mg²⁺] and [Ca²⁺] and the same environmental (boundary) conditions (simulation PI-C-ocean in Table 1). We create a configuration of the cGENIE model in which the only differences that can occur in marine carbon cycling and carbonate geochemistry are due to changing the carbonate dissolution constants. Specifically, this involves: (a) utilizing the non-seasonal ocean-atmosphere-only preindustrial (278 ppm CO₂) configuration of Ridgwell et al. (2007) so as to minimize any difference between climate states for the same imposed value of atmosphere CO₂ (which can arise in the seasonal model configuration of Cao et al. (2009) as a result of highly non-linear interactions between seasonal sea-ice extent and local ocean-atmosphere climate state), and (b) imposing a fixed (annual mean) pattern of CaCO₃ : POC export rain ratio derived from Ridgwell et al. (2007) to remove feedbacks on carbonate chemistry arising from a response of CaCO₃ export to carbonate saturation, (e.g., as per Ridgwell, 2007, 2009). We employ fixed and spatially uniform remineralization profiles of POC and CaCO₃ within the ocean interior (Ridgwell et al., 2007) and instantaneous remineralization of the residual fluxes occurs at the sea floor (a “reflective” boundary condition, Hülse et al., 2017). In order to explore parameter space in making the comparisons we take advantage of the fact that a single model experiment provides an array of combinations of environmental (T, S) conditions – a total of 934 surface grid cells with varying permutations of T and S spanning ∼ −2 to 37 °C and ∼ 33–39 PSU, respectively. This reveals the error introduced by interpolating from a lookup table rather than directly calculating the equilibrium constants at a wide range of seawater conditions.

Secondly, we investigate to what extend the new K^*s alter the marine carbonate system and carbon cycling in the model. For these purposes, we produce steady state carbon cycle simulations with three sets of K^*s (summarized in Table 1): (A) the default set of carbonate constants (Mehrbach et al., 1973) in cGENIE plus Ben-Yaakov and Goldhaber (1973) and Tyrrell and Zeebe (2004) $\mathrm{Mg}/\mathrm{Ca}$ correction scheme – the current cGENIE model default, (B) the new carbonate system corrections based on MyAMI (Hain et al., 2015) which are calibrated to the total pH scale (Dickson and Millero, 1987) and then converted to seawater pH scale assuming internally consistent modern seawater [SO${}_{\mathrm{4}}^{\mathrm{2}-}$]_T and [F⁻]_T total concentrations, and (C) the default carbonate constant set (Mehrbach et al., 1973) in cGENIE but no $\mathrm{Mg}/\mathrm{Ca}$ correction scheme.

Table 1 Summary of our simulation ensemble.

We apply these sets of K^*s in three different configurations of ocean-atmosphere-sediment carbon cycling: “ocean”, “closed” and “open”. The “ocean” configuration accounts only for ocean-atmosphere exchange (and ocean mixing and redistribution of dissolved carbon and alkalinity) and is used for the lookup table evaluation as described above. In this configuration, we are able to assess the direct effect of $\mathrm{Mg}/\mathrm{Ca}$ corrections on the dissolved carbonate system in a 3D ocean.

Next, we configure an ocean-atmosphere-sediment system as per Ridgwell and Hargreaves (2007). This is the same non-seasonally-forced ocean configuration as Ridgwell et al. (2007), but with the imposition of the same fixed spatial field of CaCO₃ : POC as used above. However, instead of dissolving all CaCO₃ that reaches the ocean-sediment interface, CaCO₃ now enters the upper sediment layers and is dissolved in situ or buried depending on local porewater chemistry and calculated by the 1D reaction-transport model of Archer (1991). Loss from the ocean of DIC, ALK, and Ca²⁺ through CaCO₃ burial is automatically balanced (tracked) by a continuously changing and equal input flux of Ca²⁺, DIC and ALK into the surface ocean routed via the runoff scheme (Colbourn et al., 2013). The result is that the total ocean inventories of ALK and Ca²⁺ do not change (hence “closed” configuration) although CO₂ can independently exchange with the atmosphere depending on the state of the biological pump (and hence the DIC inventory can evolve). In this configuration, we assess the impact of carbonate chemistry scheme and correction on CaCO₃ deposition and burial.

In the real Earth system, marine inputs would not instantly adjust to changes in marine burial but rather the marine carbonate system would adjust to re-equilibrate marine burial and inputs (the carbonate compensation feedback on 1–10 kyr time-scales – e.g., Ridgwell and Zeebe, 2005). We assess this full system response by configuring the model as an “open” system in which while the solute input (from weathering) is fixed, ocean chemistry and CaCO₃ burial dynamically adjusts to balance the (fixed) input flux. We set the fixed input flux of dissolved CaCO₃ to 10×10¹² mol C yr⁻¹ and hence close to modern open ocean burial (Ridgwell and Hargreaves, 2007). Note that in the “closed” and “open” configurations, we deviate from Ridgwell and Hargreaves (2007) by substituting the 1D reaction-transport model of Archer (1991) in place of the original sediment dissolution look-up tables (see Sect. 2.3). (It should be noted that having substituted the sediment model component (and carbonate chemistry scheme) we do not attempt to re-tune the global distribution of core-top sediment composition (or burial) – this will be the focus of a future paper.

Finally, we apply major ion concentration assumptions representing both modern ([Ca²⁺] = 10.2 mmol kg⁻¹, mean [Mg²⁺] = 52.8 mmol kg⁻¹) and idealized early Eocene ([Ca²⁺] = 20.0 mmol kg⁻¹, mean [Mg²⁺] = 30.0 mmol kg⁻¹, the same as in Hain et al., 2015), to all combinations of carbonate chemistry scheme and ocean/closed/open configuration for a total of 3 × 3 × 2 (18) permutations of experiments. Of these, we list in Table 1 the subset of 12 permutations that we focus on in the Results, plus one simulation (PI-C2-ocean) testing the importance of user-set temperature limits for the simulated carbonate system.

3 Results and Discussion

We start our presentation and discussion of the results with an assessment of the accuracy of our lookup-table implementation of MyAMI-derived interpolated K^* values (Sect. 3.1), followed by an evaluation of the impacts on the carbonate system parameters and marine carbon cycling of changing from the default set of carbonate constants to the MyAMI-derived constant set – both experiments conducted under modern (PI) seawater [Mg²⁺] and [Ca²⁺] and hence implicitly with no $\mathrm{Mg}/\mathrm{Ca}$ correction being applied (Sect. 3.2). In Sect. 3.3 for assumed Early Eocene (EE) seawater [Mg²⁺] and [Ca²⁺] conditions, we assess the implications for ocean carbonate chemistry as well as CaCO₃ preservation and burial in marine sediments, of applying the default cGENIE and MyAMI-derived $\mathrm{Mg}/\mathrm{Ca}$ correction schemes and vs. no applied correction. We end with a discussion of the implications for the interpretation and data assimilation of paleoenvironmental proxies (Sect. 3.4).

3.1 Evaluation of the interpolated lookup approximation of MyAMI carbonate constant characteristics

The values of $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ calculated directly using MyAMI are compared with the internally interpolated values in cGENIE for each of the model surface ocean grid cells (934 points) and for the specific combination of T, S, [Mg²⁺], and [Ca²⁺] of simulation PI-C-ocean with pre-industrial $\mathrm{Mg}/\mathrm{Ca}$ (Fig. 2). The differences are vanishingly small (∼ 0.02 % at most) and two to three orders of magnitude smaller than their respective spatial variability. Indeed, the interpolation errors are smaller than the precision of MyAMI (Hain et al., 2015) and the empirical uncertainty of the constants (Orr et al., 2018) and we thus conclude that the interpolated carbonate constants could not be distinguished from the directly calculated ones in the real world. Still, we analyse the errors to understand how the design of the discrete, evenly sampled look-up table and multi-dimensional linear interpolation affects the carbonate constants applied to cGENIE. At a given combination of T, S, [Mg²⁺], and [Ca²⁺] the total error is composed of the errors caused by linear interpolation in each dimension. The largest error is caused by the interpolation that least captures the functionality in the respective dimension, either because that functionality is poorly approximated by a linear fit or because the sampling in that dimension is too coarse to capture its complexity. For $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$, the interpolation of the salinity dependence causes the largest errors, evidenced by the sinusoidal distribution of the error on the salinity grid, with the largest errors occurring at half-distance between the sampled salinities (every 1 PSU). This suggests that the coarse sampling in the salinity-space dominates the interpolation error in $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$. In contrast, for $K_{\mathrm{2}}^{*}$, the interpolation along the temperature axis creates the largest errors, again distributed sinusoidally (Fig. 2k). The interpolation errors are generally larger for $K_{\mathrm{2}}^{*}$, which is more sensitive to [Mg²⁺] and [Ca²⁺] changes than $K_{\mathrm{1}}^{*}$ due to the high complexation potential of the carbonate ion (Hain et al., 2015, 2018). Overall, this reaffirms our assertion that use of an interpolation method to substitute for the underlying MyAMI code does not introduce substantial errors.

Figure 2 Errors due to interpolating rather than calculating in situ K^*s in the surface ocean in simulation PI-C-ocean. All errors are given as percentage of the exact MyAMI value.

3.2 Impact of changing carbonate system constants (cGENIE default vs. MyAMI-derived)

We next turn to evaluating any differences in the marine carbon cycle that might be induced by using the MyAMI-derived carbonate system parameters instead of the ones calculated based on Mehrbach et al. (1973) (the default parameters in cGENIE), assuming modern seawater (PI) and with cGENIE configured in an ocean(/atmosphere)-only configuration (PI-C-ocean vs. PI-B-ocean). In addition, MyAMI is calibrated to the standard “best practice” carbonate equilibrium constants from Dickson et al. (2007) instead of Mehrbach et al. (1973). The Dickson constant set, though originally based on Mehrbach et al. (1973)'s measurements, is calibrated to the total pH scale and intended for use in modern seawater chemistry. Thus, the shift from the Mehrbach-PI ocean to the MyAMI-PI ocean theoretically improves cGENIE's representation of preindustrial ocean chemistry, making it more directly comparable to empirical measurements that employ Dickson's standard practice K^*s. The equilibrium constants of Mehrbach et al. (1973) vs. Dickson et al. (2007) have slightly different temperature and salinity dependencies, causing differences between the default and new cGENIE PI state. Overall, the new dependence on local [Mg²⁺] and [Ca²⁺] and the change in parameter set introduce larger differences than those caused by the interpolation but they are still minimal compared to the spatial variability of these constants, as we will show next (Fig. 3). Though the differences are small, we describe their spatial pattern in the surface ocean to provide a sense of their effect on the simulated carbonate system.

Figure 3 Surface ocean values of $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$, and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ in default cGENIE (PI-B-ocean) and cGENIE + MyAMI (PI-C-ocean) solutions for the pre-industrial carbonate system (a–f) and difference between the two (simulations PI-C-ocean minus PI-B-ocean) for the surface and benthic ocean. The surface layer represents the uppermost 175 m of the water column and benthic values are from the deepest ocean grid box in every location. In the “ocean” set up, organic and carbonate particles are instantaneously respired/dissolved when reaching the sediment-water interface. In shallow locations (less than 600 m water depth) where the particle flux is large, this results in unusual local carbonate system states that appear as brown cells in (j)–(l).

Differences between default and MyAMI-derived K^*s are largest at the margins of the sampled parameter range, particularly in the cold and relatively fresh Arctic and salty Mediterranean (Fig. 3). In the Arctic, $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$ and $K_{\mathrm{B}}^{*}$ are all lower, which reduces [H⁺] and thus increases pH (Fig. S1). These changes cause a reduction in [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] and a small increase in [HCO${}_{\mathrm{3}}^{-}$] despite the lower $K_{\mathrm{1}}^{*}$, which drives a marginal overall DIC increase in Arctic waters. Reductions of $K_{\mathrm{sp},\mathrm{cal}}^{*}$ coincident with the $K_{\mathrm{2}}^{*}$ changes lessen the impact on Ω, which consequently shows little change. These high-latitude differences are strongly affected by the different temperature limits for the default Mehrbach scheme and the MyAMI-derived constants and are not present when we apply the same temperature limits in both cases (Fig. S3). The temperature limits do not affect the rest of the surface ocean as sea-surface temperatures outside the high-latitudes are within the temperature limits of both schemes.

Outside of the Arctic, MyAMI-derived surface ocean $K_{\mathrm{1}}^{*}$ is much closer to the default scheme values, with the next largest differences being in the Mediterranean, followed by the Atlantic. The difference in $K_{\mathrm{2}}^{*}$ shows a similar pattern, but with an increase rather than decrease in the Mediterranean. Thus, in the Mediterranean, the dissociation of carbonic acid is less favourable and the dissociation of [HCO${}_{\mathrm{3}}^{-}$] remains unchanged or is slightly increased, respectively, resulting in less [HCO${}_{\mathrm{3}}^{-}$] and more [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] (Fig. S1). Despite increased [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], Ω is lower in Mediterranean and Atlantic surface waters because of a higher $K_{\mathrm{sp},\mathrm{cal}}^{*}$ (Fig. 3i).

In the Western Pacific surface, the $K_{\mathrm{2}}^{*}$ decrease outcompetes the effect of reduced $K_{\mathrm{1}}^{*}$, like in the Arctic, and thus leads to increased [HCO${}_{\mathrm{3}}^{-}$] and reduced [CO${}_{\mathrm{3}}^{\mathrm{2}-}$].

DIC anomalies in the surface are entrained into deep water, impacting carbon speciation throughout the ocean interior in addition to changes arising from the use of different constants. For instance, in the deep Atlantic, higher [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] in North Atlantic deep water compensates for the slightly increased $K_{\mathrm{sp},\mathrm{cal}}^{*}$ to cause a small Ω increase (Fig. S1) while in most of the abyssal Indo-Pacific, which is ventilated by southern-sourced waters with no or negative [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] differences, the $K_{\mathrm{sp},\mathrm{cal}}^{*}$ change dominates and slightly reduces Ω.

Figure 4 Absolute distributions and differences in the carbonate content in surface sediments and burial rates with the default scheme (simulation PI-B-open) and cGENIE + MyAMI (PI-C-open) for pre-industrial $\mathrm{Mg}/\mathrm{Ca}$. Sediments shallower than 1000 m are masked in black and locations where the sedimentary carbonate system could not be solved dynamically and was instead replaced by a value derived from the lookup table of Ridgwell et al. (2003) is hatched in yellow.

When we include marine sediments in our model configuration (PI-C-closed vs. PI-B-closed), any adjustments made to the marine carbonate system constants are explicitly propagated into the 1D reaction-transport sediment model (Archer, 1991). The benthic Ω decrease across most of the Indo-Pacific leads to a small reduction of CaCO₃ accumulation, and ultimately burial, in sediments (Fig. 4). In the Atlantic, the slightly increased Ω results in increased CaCO₃ preservation. In a fully open system (PI-C-open vs. PI-B-open), the changed carbonate preservation pattern alters the flux of DIC, calcium and alkalinity back into the ocean, which shifts the marine carbonate system until total burial fluxes are re-equilibrated with the terrestrial ALK and DIC supply, which in our setup remains constant. These feedback relationships constitute the “carbonate compensation” dynamic of the open system carbon cycle (Ridgwell and Zeebe, 2005; Hain et al., 2025). In the case of modern seawater (PI), these shifts are small in the global average. The global DIC inventory of the ocean is 0.06 % (24 PgC) higher and mean weight percent of carbonate in surface sediments is −0.12 % lower with the MyAMI-derived K^*s in the open system. Locally, the differences can be slightly higher, especially in the transitory depth zone between full carbonate preservation and full carbonate dissolution, resulting in slightly increased burial in the Indian Ocean and open Pacific, decreased burial rates in the Atlantic, and altered sediment composition in a few isolated grid cells (Fig. 4).

3.3 Implications of applying $\mathrm{Mg}/\mathrm{Ca}$ carbonate system corrections for Early Eocene-like seawater composition

Next, we assess these three main impacts of a [Ca²⁺] and [Mg²⁺] change (CaCO₃ saturation change, equilibrium carbonate speciation change, carbonate ion complexation) in the context of both closed (i.e. constant carbon and ALK), and open (variable carbon and ALK) atmosphere-ocean system and under seawater $\mathrm{Mg}/\mathrm{Ca}$ representative of the Early Eocene.

Figure 5 Changes of the global mean surface values of selected dissociation constants and carbonate system metrics due to changing [Mg²⁺] and [Ca²⁺] from pre-industrial to Early Eocene-like in simulations without corrections (simulations EE-A-closed minus PI-A-closed, on the left side of each set of bars) default cGENIE correction (simulations EE-B-closed minus PI-B-closed, middle bar in each set) and cGENIE + MyAMI (EE-C-closed minus PI-C-closed, on the right side of each pair of bars). The changes are given as percentage of the PI values.

The impact of changing seawater major ion composition in a closed system under (i) default (EE-B-closed), (ii) MyAMI-based (EE-C-closed), and (iii) no $\mathrm{Mg}/\mathrm{Ca}$ correction (EE-A-closed) are shown in Fig. 5. Ignoring for now the response of Ω (discussed later), out of the three effects we find that the complexation of carbonate ions is the most significant. This is because the calcium and magnesium concentration changes effectively reduce divalent cation concentration by about 20 %, thereby significantly increasing the activity coefficient of total carbonate ion and reducing $K_{\mathrm{2}}^{*}$. For example, at modern surface DIC/ALK the decrease in carbonate complexation effectively raises pH and repartitions alkalinity from carbonate ion to borate ion, driving 5 %–10 % reductions in total carbonate ion, CO₂ and [H⁺]_total and corresponding increases in bicarbonate (Fig. 5) and borate. The differences between the K^*s in the uncorrected cGENIE simulation (EE-A-closed) and those in the cGENIE + MyAMI simulation (EE-C-closed) also vary spatially, with the largest $K_{\mathrm{1}}^{*}$ differences in cold waters of the deep ocean and polar surface oceans, the largest $K_{\mathrm{2}}^{*}$ differences in the warm tropical surface waters and the largest $K_{\mathrm{sp},\mathrm{cal}}^{*}$ differences in the deep ocean (Fig. 6).

Figure 6 Differences in $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$, and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ between uncorrected cGENIE and cGENIE + MyAMI solutions for the carbonate system of the surface and benthic ocean with Eocene-like $\mathrm{Mg}/\mathrm{Ca}$ (simulations EE-C-closed minus EE-A-closed). In (d)–(f), locations with sediments shallower than 1000 m are masked in black and locations where the sedimentary carbonate system could not be solved dynamically and was instead replaced by a value derived from the lookup table of Ridgwell et al. (2003) is hatched in yellow.

Because the magnitude of the equilibrium constant corrections is significant for Eocene seawater, the choice of correction scheme also becomes relevant. Hain et al. (2015) showed that the equilibrium constants adjusted with MyAMI differ from those adjusted with the cGENIE default schemes, especially for large deviations from PI conditions, because of different sensitivities to [Ca²⁺] and [Mg²⁺]. Specifically, $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ decrease more (PI to EE) with the standard corrections than with those calculated by MyAMI (Fig. 5; Fig. 2 of Hain et al., 2015). The coincident change in [Ca²⁺] and [Mg²⁺] only causes a minor (2 %) decrease in $K_{\mathrm{1}}^{*}$ with MyAMI, suggesting that the changes of hydrogen ion and bicarbonate ion activity almost cancel out. Carbonate ion activity, however, is sensitive to changing ionic strength (e.g., Garrels and Thompson, 1962; Pytkowicz and Hawley, 1974). Any reduction in the total divalent cation concentration of seawater – as in the change from modern to Eocene – will tend to reduce the fraction of total carbonate ion that is complexed into stable ion-pairs with the divalent cations. That is, in modern seawater 36 % of total carbonate ion molecules are free (64 % complexed), and with the lower Eocene total divalent cation concentration 39.5 % of total carbonate ion molecules are free – corresponding to a ∼ 10 % increase in the total carbonate ion activity coefficient, and hence a decrease in $K_{\mathrm{2}}^{*}$. Carbonate ion activity also dominates the response of the CaCO₃ solubility constants $K_{\mathrm{sp},\mathrm{cal}}^{*}$, and both $K_{\mathrm{2}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ need to decrease by the same amount to account for changes in carbonate complexation, as is the case for MyAMI (Hain et al., 2015) but in contrast to the combination of $K_{\mathrm{sp},\mathrm{cal}}^{*}$ of Tyrrell and Zeebe (2004) and $K_{\mathrm{2}}^{*}$ of Ben-Yaakov and Goldhaber (1973).

When considering the effects of changing speciation and seawater major ion composition on the CaCO₃ saturation index Ω there are additional factors to consider. In approximately doubling Eocene seawater [Ca²⁺] relative to modern and in the absence of any change in the concentration of total carbonate ion, we would expect Ω to increase proportionately. Indeed, for no applied $\mathrm{Mg}/\mathrm{Ca}$ correction and virtually no change in $K_{\mathrm{sp},\mathrm{cal}}^{*}$ or [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] (which in any case largely cancel out – see Eq. 28), this is what we observe in the model (+98 %, Fig. 5). Including a $\mathrm{Mg}/\mathrm{Ca}$ correction modifies this response. With the MyAMI-based carbonate constants, there is an additional +10 % increase in Ω due to reduced carbonate complexation (at constant total carbonate ion), which, coupled with a 7 % decrease resulting from the simulated total carbonate ion decline, gives a net 103 % increase in Ω in our simulations (Fig. 5). In contrast, with the $K_{\mathrm{sp},\mathrm{cal}}^{*}$ correction factor of Tyrrell and Zeebe (2004) as previously implemented in cGENIE (e.g., Panchuk et al., 2008) the Ω increase is ∼ 2.7-fold (+170 %), or about two thirds greater than computed with MyAMI (Hain et al., 2015; Fig. 5), and with only a small portion of this attributable to the difference in [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] decrease. Hence, using the Tyrrell and Zeebe (2004) correction factor for the CaCO₃ solubility constants $K_{\mathrm{sp},\mathrm{cal}}^{*}$ may lead to significant biases compared to the use of MyAMI which explicitly includes carbonate complexation by divalent cations.

Figure 7 Effect of changing $\mathrm{Mg}/\mathrm{Ca}$ from pre-industrial to Eocene-like for mean ocean DIC and ALK changes, mean surface pH and Ω, average CaCO₃ content of marine sediments and total CaCO₃ burial in simulations without correcting for major ion changes, default cGENIE and cGENIE + MyAMI.

In our closed system configuration, simulated Ω changes drive very large increases in CaCO₃ preservation (reflected in core-top wt % CaCO₃) and hence burial (three brown bars in Fig. 7b, e, f), with the global sediment accumulation rate of CaCO₃ increasing disproportionately in response to Ω – from 9.0 Tmol yr⁻¹ (PI-A-closed) to 27.6 Tmol yr⁻¹ (EE-A-closed, Fig. 7f shows the difference of the two simulations) Because the [Ca²⁺] effect on Ω is dominant (see above) we find that the different $\mathrm{Mg}/\mathrm{Ca}$ correction schemes exert only a relatively minor modulation of CaCO₃ wt % and burial. In the closed system experiments, carbonate system changes cause minor shifts in marine carbon storage (three orange bars in Fig. 7a), with a ∼ 13 µmol kg⁻¹ (−0.006 %) mean ocean DIC loss in default cGENIE but a ∼ 18 µmol kg⁻¹ (+0.008 %) mean ocean DIC gain in cGENIE + MyAMI because the applied restoring of the atmospheric CO₂ concentration allows for net C loss or gain in the atmosphere-ocean system at constant ALK.

In contrast, if we configure cGENIE with a constant global terrestrial weathering rate (and hence invariant fluxes of DIC and ALK to the ocean) in an “open” system configuration, the initial imbalance in (enhanced) CaCO₃ burial vs. (fixed) weathering drives the DIC and ALK composition of the ocean lower and away from that of the closed system (Fig. 7a, d). Once the system has re-balanced with CaCO₃ burial equal to weathering (we run our experiments for 100 kyr to achieve this), DIC, ALK, and pH are all lower than was the case for PI [Ca²⁺] and [Mg²⁺], while Ω is slightly higher (three blue bars, Fig. 7a–d). For no applied $\mathrm{Mg}/\mathrm{Ca}$ correction and hence a purely [Ca²⁺] induced reorganization of marine carbon cycling (EE-A-open minus PI-A-open), the changes are respectively: −465 µmol kg⁻¹ DIC (−0.21 %, Fig. 7a) and −538 µmol kg⁻¹ ALK (−0.22 %, Fig. 7d) (on a mean global ocean), −0.094 pH units (Fig. 7c) and +1.28 Ω (Fig. 7b, all on a global ocean surface mean basis). When applying the MyAMI-based $\mathrm{Mg}/\mathrm{Ca}$ correction (EE-C-open minus PI-C-open), we find that the decreases in DIC and ALK are slightly reduced (by 20 and 12 µmol kg⁻¹, respectively, Fig. 7a, d), resulting in a slightly greater increase in Ω but smaller decline in pH. In contrast, the default cGENIE $\mathrm{Mg}/\mathrm{Ca}$ correction scheme (EE-B-open minus PI-B-open) results in the carbonate saturation increase almost being doubled (Fig. 7b), and the pH decrease halved (Fig. 7c) – differences that largely occur because $K_{\mathrm{1}}^{*}$, $K_{\mathrm{2}}^{*}$ and $K_{\mathrm{sp},\mathrm{cal}}^{*}$ decrease less when corrected with MyAMI rather than the default scheme (Fig. 5). Spatially, differences between correction schemes appear across the whole non-Arctic surface ocean and are largest in warm waters for $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$ and saltier waters for $K_{\mathrm{sp},\mathrm{cal}}^{*}$ (Fig. S4). Consequently, pH is 0.04–0.06 lower and Ω reduced by 0.8–2.5 across most of the surface ocean when the MyAMI rather than the default scheme is used in our simulations (Fig. S5).

3.4 Implications for interpreting paleoenvironmental proxies

By updating the simulated carbonate system in cGENIE, the MyAMI-derived constants also adjust the representation of model variables that can be directly compared to paleoenvironmental proxies. Examples are the CaCO₃ fraction of marine sediments and marine CaCO₃ burial, which have been used to constrain past marine carbonate system states (e.g. Panchuk et al., 2008; Si et al., 2023; Li et al., 2024). For Eocene-like [Mg²⁺] and [Ca²⁺], simulated CaCO₃ fraction of marine sediments and marine CaCO₃ burial are slightly different with MyAMI-derived constants (simulation EE-C-open) than with the default correction scheme (simulation EE-B-open, Fig. S5e–f) due to the discussed changes in deep ocean Ω, highlighting the potential for a small error in model-data comparisons without the new scheme.

In addition to correcting the representation of major ion effects on marine carbon cycling, the new carbonate system correction also enables a more direct comparison of simulated and reconstructed pH (when both are reported on the total scale). This is because the $\mathrm{Mg}/\mathrm{Ca}$ correction of B alkalinity adds to the differences of surface ocean carbonate system solutions between MyAMI and the default cGENIE carbonate system. We demonstrate this by comparing the marine dissolved boron (B) speciation that is simulated in cGENIE with the default vs. new MyAMI-derived scheme.

Figure 8 Zonally-averaged profiles of differences in borate ion concentrations and δ¹¹B of borate derived from $K_{\mathrm{B}}^{*}$ and pH between cGENIE + MyAMI and default cGENIE (simulations EE-C-open minus EE-B-open) for the Pacific. The underlying pH differences are shown in Fig. S6c.

cGENIE accounts for borate concentrations in the total alkalinity, but the dissociation constant for boric acid $K_{\mathrm{B}}^{*}$ is not corrected for $\mathrm{Mg}/\mathrm{Ca}$ in the default scheme. The new scheme imports $K_{\mathrm{B}}^{*}$ for local conditions from MyAMI. In our simulation EE-C-open, this leads to a slightly lower surface $K_{\mathrm{B}}^{*}$ in the global average than under PI conditions (Fig. 5), and hence the carbonate system is solved with lower borate ion concentrations than in simulation EE-B-open with the default scheme. Figure 8 shows that the borate concentration differences are highest in the surface ocean where they reach up to −12 µmol kg⁻¹.

Assuming modern-day δ¹¹B of seawater (39.61 ‰) and a α value of 1.0272 (Klochko et al., 2006), we estimate δ¹¹B of borate from the difference between p$K_{\mathrm{B}}^{*}$ and pH (Fig. 8):

$$\begin{array}{rl}{\mathit{\delta }}^{\mathrm{11}}{\mathrm{B}}_{\mathrm{Borate}}& =(\mathrm{39.61}\mathrm{‰}\times (\mathrm{1.0}+{\mathrm{10}}^{\mathrm{p}{K}_{\mathrm{B}}^{\ast }-\mathrm{pH}})\\ & -\mathrm{27.2}\times {\mathrm{10}}^{\mathrm{p}{K}_{\mathrm{B}}^{\ast }-\mathrm{pH}})/(\mathrm{1.0}+\mathit{\alpha }\times {\mathrm{10}}^{\mathrm{p}{K}_{\mathrm{B}}^{\ast }-\mathrm{pH}})\end{array}$$

The differences in the local dissociation constants for boric acid translate into differences of up to 0.5 ‰ in the δ¹¹B of borate ions in the surface ocean, corresponding to up to 0.04 pH units. This highlights the relevance of $\mathrm{Mg}/\mathrm{Ca}$ corrections for B-based pH estimates and shows how internal inconsistencies can arise when different $\mathrm{Mg}/\mathrm{Ca}$ corrections are used in comparisons of local pH reconstructions and simulated carbonate systems in Earth system models.

4 Conclusions

We implemented a new carbonate chemistry scheme option in the cGENIE Earth system model to correct the simulated carbonate system equilibria based on the MyAMI model and numerically implemented this via look-up tables. We evaluated the accuracy of the scheme by comparing our interpolated carbonate system parameters to those directly calculated with MyAMI. We then assessed the effects of the new scheme on the simulated carbonate system in cGENIE. Using MyAMI-derived carbonate system parameters introduces small differences in the simulated pre-industrial carbonate system. When simulating a carbonate system which requires ion-pairing corrections, in our case by changing marine [Ca²⁺] and [Mg²⁺] to Early Eocene-like values, the new scheme results in lower surface ocean saturation state and pH increases than the default scheme in a closed system and less ALK and DIC loss and a larger surface ocean pH decline in an open system than the default scheme. These differences demonstrate the importance of updating the default cGENIE carbonate system corrections for major ion concentrations and exemplify the systematic bias that exists when comparing the carbonate system simulated with cGENIE's default scheme to observations.