Effect of the Atlantic Meridional Overturning Circulation on atmospheric pCO₂ variations

Proxy records show large variability of atmospheric pCO₂ on different timescales. Most often such variations are attributed to a forced response of the carbon cycle to changes in external conditions. Here, we address the problem of internally generated variations in pCO₂ due to pure carbon cycle dynamics. We focus on the effect of the strength of Atlantic Meridional Overturning Circulation (AMOC) on such internal variability. Using the Simple Carbon Project Model v1.0 (SCP-M), which we have extended to represent a suite of nonlinear carbon cycle feedbacks, we efficiently explore the multi-dimensional parameter space to address the AMOC–pCO₂ relationship. We find that climatic boundary conditions and the representation of biological production in the model are most important for this relationship. When climate sensitivity in our model is increased, we find intrinsic oscillations due to Hopf bifurcations with multi-millennial periods. The mechanism behind these oscillations is clarified and related to the coupling of atmospheric pCO₂ and the alkalinity cycle, via the river influx and the sediment outflux. This mechanism is thought to be relevant for explaining atmospheric pCO₂ variability during glacial cycles.

1 Introduction

Atmospheric pCO₂ values show large variations on many different timescales. Over the Cenozoic, pCO₂ values gradually decreased from values of up to 2500 ppmv in the Eocene to 300 ppmv at the end of the Pliocene. When considering the Pleistocene glacial–interglacial cycles, one of the remarkable results is the strong correlation between pCO₂ and temperature, with dominant variations of about 100 ppmv in 100 000 years, as reconstructed from ice cores (Petit et al., 1999). Over the industrial period, pCO₂ values have increased by 130 ppmv due to human activities (Friedlingstein et al., 2020). This forced trend is superposed on natural variability associated with the seasonal cycle and longer-timescale climate variability (Gruber et al., 2019). The effect of the natural variability is much lower than the forced trend on such relatively short timescales. For example, the El Niño–Southern Oscillation (ENSO), a dominant mode of interannual climate variability, induces atmospheric pCO₂ variations of only 1–2 ppmv (Jiang and Yung, 2019). Most studies seek to explain such variations in pCO₂ as a forced response of the carbon cycle to changes in external conditions. For example, glacial cycles are thought to be caused by orbital variations in insolation, possibly amplified by physical processes in the climate system (Muller and MacDonald, 2000). Such variations in temperature (and other quantities, e.g., precipitation) then affect the carbon cycle, leading to changes in pCO₂. On the other hand, changes in pCO₂ will affect global mean temperature and hence may amplify any temperature anomaly. Hence it is questionable whether the pCO₂ response to orbital insolation changes can be considered as a solely forced response, with no internal dynamics of the carbon being involved (Rothman, 2015).

The carbon cycle is comprised of an extremely complex entangled set of processes which act in the different components of the climate system (e.g., land, ocean) on many different timescales. The marine carbon cycle, with its three main carbon pumps, is a major player in this cycle, at present day resulting in the uptake of about 25 % of human-released emissions (Sabine et al., 2004). The carbon pumps involve physical processes, biological processes, and processes in ocean sediments. Many carbon cycle feedbacks exist, either between only physical quantities or between biological and physical quantities. An example of such a feedback is the solubility feedback: for higher atmospheric pCO₂, solubility of CO₂ decreases due to higher ocean temperatures, resulting in relatively less CO₂ uptake by the ocean and thus relatively higher atmospheric pCO₂. Given this strongly nonlinear system, it would be strange if it would not show strong internal variability, i.e., variability which would exist even if the carbon cycle system were driven by a time-independent external forcing. There are indeed examples (Rothman, 2019), where oscillatory behavior in the carbon cycle has been attributed to internal carbon cycle dynamics.

The physical context of all carbon pumps is the three-dimensional ocean circulation, which can be roughly decomposed in a wind-driven and an overturning component, the latter strongly related to the deep-ocean circulation. The Atlantic Meridional Overturning Circulation (AMOC) is a major component of the global overturning circulation because of its associated meridional transport of heat, salt, and nutrients.

The relation between the AMOC and atmospheric pCO₂ is complicated. A direct effect of a changing AMOC is a change in the distribution of tracers such as temperature, dissolved inorganic carbon (DIC), alkalinity (Alk), and nutrients. For example, after an AMOC weakening, the distributions of these tracers affect biological export production via reduced nutrient upwelling (Marchal et al., 1998; Menviel et al., 2008; Mariotti et al., 2012; Nielsen et al., 2019) and gas exchange via changing solubility of CO₂ in the ocean (Menviel et al., 2014). Besides these direct effects, the AMOC also influences mixing in the Southern Ocean. Changes in this mixing due to a weaker AMOC can result in a higher outgoing flux of carbon to the atmosphere (e.g., Schmittner et al., 2007; Huiskamp and Meissner, 2012; Menviel et al., 2014). Furthermore, changes in the AMOC also influence the general sensitivity of the marine carbon cycle to, for example, changes in the wind field (Munday et al., 2014). These processes form a complex puzzle where the sign of atmospheric pCO₂ change following an AMOC strength change is difficult to determine. Currently, different models produce different results with respect to the sign of the atmospheric pCO₂ change, which can be attributed to the assessed timescale, model used, and what climatic boundary conditions are used (Gottschalk et al., 2019).

On the other hand, pCO₂ also influences the AMOC (Toggweiler and Russell, 2008), and present-day climate models forced with anthropogenic emissions simulate a weaker AMOC for larger atmospheric pCO₂ (Gregory et al., 2005; Weijer et al., 2020). By contrast, proxy data suggest that in the Last Glacial Maximum both atmospheric pCO₂ and the strength of the AMOC were lower (Duplessy et al., 1988). This shows that there is also a sensitivity to climatic boundary conditions in the relation (Zhu et al., 2015) between the AMOC and atmospheric pCO₂. The AMOC can also display tipping behavior (Weijer et al., 2019) under an increase in pCO₂, which can have large effects on climate. Examples of these effects are disrupted heat transport (Ganachaud and Wunsch, 2000), changing precipitation patterns (Vellinga and Wood, 2002), and a different distribution of important tracers in the ocean. Such tipping can hence have strong consequences on the carbon cycle, and hence on atmospheric pCO₂.

In this paper, we perform a systematic study of internal carbon cycle variability and the relation AMOC–pCO₂ connection, using the Simple Carbon Project Model v1.0 (SCP-M). This model (O'Neill et al., 2019) simulates the most important carbon cycle processes in a simple global ocean box structure. The simple box setup enables us to efficiently scan the parameter space of the carbon cycle model using parameter continuation methods. With this approach we aim to answer the following three questions: (i) how does atmospheric pCO₂ respond to changes in the strength of a constant (in time) AMOC? (ii) Does the pCO₂–AMOC feedback lead to new variability phenomena? And (iii) are there tipping points and internal oscillations in the carbon cycle?

When answering these questions, we pay special attention to different (non-linear) carbon cycle feedbacks. We will also use two different model configurations to take account of different climatic boundary conditions, the pre-industrial (PI) configuration and the Last Glacial Maximum (LGM) configuration. The SCP-M, its configurations, the different additional feedbacks implemented, and the parameter continuation approach are described in Sect. 2. In Sect. 3 we present the results of the different cases considered, and we conclude the paper with a summary and discussion in Sect. 4.

2 Methods

2.1 SCP-M

The SCP-M is a carbon cycle box model focused on the marine carbon cycle. Because of its simple structure, it is well suited to test high-level concepts in both modern and past configurations. In the ocean several tracers are resolved. In this study we will only use the three most important tracers, DIC, Alk, and phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$), to reduce the problem size. In the original model there is also a terrestrial biosphere component, and several sources of CO₂ to the atmosphere. We will not use these, since our focus is on the marine carbon cycle. The processes which are resolved are ocean overturning circulation, sea–air gas exchange, biological production, calcium carbonate (CaCO₃) production and dissolution, and river and sediment fluxes.

Figure 1 The box structure and fluxes for the SCP-M based on O'Neill et al. (2019). ψ₁ (red) is the GOC, ψ₂ (orange) is the AMOC, and γ₁ and γ₂ (blue) represent bidirectional mixing. Biological fluxes are represented by the green arrows, calcifier fluxes by the light gray arrows, and general dissolution of calcium carbonate by the gray wiggles in the boxes. k_w (gray) represents the gas exchange between the ocean and the atmosphere. Lastly, there is an influx in box 1 via the rivers (yellow) and an outflux to the sediments (light gray).

The model consists of eight boxes: one atmospheric box and seven oceanic boxes (Fig. 1). This means that the sediment stock is not explicitly solved for in the model. In the model, the ocean boxes are differentiated on latitude and depth. Consequently, there is no longitudinal variation, and no differentiation between ocean basins. The used boxes are (1) a low-latitude surface box, (2) a northern high-latitude surface box, (3) an intermediate ocean box, (4) a deep-ocean box, (5) a southern high-latitude surface box, (6) an abyssal ocean box, and (7) a sub-polar surface box. This division in the ocean is based on regions in the ocean where the water masses have similar characteristics. The different boxes are connected via ocean circulation and mixing, which is based upon a conceptual view of the ocean circulation (Talley, 2013). The largest circulation is the Global Overturning Circulation (GOC; ψ₁). This circulation connects boxes 4–7 and represents the formation of Antarctic Bottom Water. Next to the GOC, the other major circulation is the AMOC (ψ₂) which connects boxes 2–4 and 7. Lastly, there is bidirectional (vertical) mixing between boxes 4 and 6 (γ₁) and boxes 1 and 3 (γ₂).

To be able to solve for several fluxes, such as the air–sea gas exchange, the pH in the ocean needs to be determined. Unfortunately, pH is not a conservative tracer, which means that we need a carbonate chemistry module to solve for pH. In the SCP-M, a direct solver is used where the pH value of the previous time step is used as an estimate for the new step (Follows et al., 2006). Using this carbonate chemistry, the model is able to determine the carbonate (${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$) concentration and oceanic pCO₂. This latter quantity is used to model the exchange of CO₂ between the atmosphere and ocean. For each surface box, the flux is proportional to the pCO₂ difference between atmosphere and ocean and a constant piston velocity (k_w).

Biological production is constant in the SCP-M. Per surface box, a constant value is used to denote the biological export production at 100 m depth. The organic flux is remineralized in the subsurface boxes following the power law of Martin et al. (1987). The biological export production is also important for the carbonate pump. Via a constant rain ratio, the biological production is linked to the production of calcifiers. Besides organic growth via photosynthesis, calcifiers also take up DIC and Alk to form shells (CaCO₃). Upon death, these calcifiers sink to deeper boxes where the shells are dissolved. The dissolution of the shells is dependent on a constant dissolution rate and a saturation-dependent dissolution. If the total dissolution of CaCO₃ in the ocean is smaller (larger) than the production in the surface ocean, there is an outflux (influx) of DIC and Alk to (from) the sediments. The river flux for ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ is constant in the SCP-M and balanced by a constant outflux into the sediments. Influx of DIC and Alk via the rivers is variable and proportional to atmospheric pCO₂ according to

$$\begin{array}{}\text{(1)}& C_{\mathrm{river}}=W_{\mathrm{SC}}+\left(W_{\mathrm{SV}}+W_{\mathrm{CV}}\right)\times {p{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{atm}},\end{array}$$

where W_SC is a parameter reflecting constant silicate weathering, W_SV a parameter representing variable silicate weathering, W_CV a parameter representing variable carbonate weathering, and pCO₂^atm the partial pressure of CO₂ in the atmosphere. This relation was already used in the SCP-M and is directly taken from Toggweiler and Russell (2008), meaning no adaptations were made to the relation or to the parameter values. This expression is adapted from a formulation used in Walker and Kasting (1992). The two different processes parameterized here, carbonate and silicate weathering, are important on different timescales. Terrestrial carbonate weathering is typically important on timescales of 10³ to 10⁴ years, whereas silicate weathering balances volcanic carbon input on timescales of 10⁵ to 10⁶ years (Archer et al., 1997).

The difference between the influx of DIC and Alk and the outflux into the sediments determines the change in total carbon and Alk in the system. The outflux of DIC and Alk via the sediments i, is related to CaCO₃ burial, which is the difference between CaCO₃ production in the surface ocean and CaCO₃ dissolution in the ocean and sediments. CaCO₃ dissolution consists of two parts: a saturation-driven part and a constant part following

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{C}_{\mathrm{diss}}& =\left(\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]\right)\times k_{\mathrm{Ca}}\\ & \times \left(\mathrm{1}-\left(\min \left(\left({\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}{K_{\mathrm{sp}}}}\right),\mathrm{1}\right)\right)\right)+D_{\mathrm{C}},\end{array}\end{array}$$

where the first part is related to the saturation state: $\frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{Ca}}^{\mathrm{2}+}\right]}{K_{\mathrm{sp}}}$. If the saturation state is larger than 1, the saturation-dependent dissolution is zero, and only the constant term remains (D_C). CaCO₃ production is linearly dependent on biological production using a rain ratio.

A big advantage of the SCP-M is that it has two configurations: a PI configuration and an LGM configuration. The parameter values in both configurations have been determined via extensive tuning of the model to observations and proxies in O'Neill et al. (2019). The configurations are differentiated on surface ocean temperature and salinity, ocean circulation, sea ice cover in box 5, and total volume of the ocean. The parameter values of the two different SCP-M configurations can be found in Table 1.

Table 1 The parameter values that are different between the two configurations (PI and LGM). In columns 1 and 2 the parameter symbol and description are given. In column 3 the PI value is given, and in column 4 the LGM value is given.

2.2 Representation of carbon cycle processes and feedbacks

The carbon cycle has many (non-linear) feedbacks which are not represented in the original SCP-M version to keep the model as simple as possible. The absence of these feedbacks can lead to non-physical behavior (e.g., negative concentrations) when parameter values, such as the AMOC strength, are changed. We have implemented several additional feedbacks which can be divided into two categories: those that mostly concern physical processes and those associated with biological processes. The feedbacks are included through parameters λ; when such a parameter is zero, the feedback is not active in the SCP-M and the original version applies. For all feedbacks, except the feedback on the rain ratio (Eq. 11 below), the sign of the feedback (positive or negative) is unclear beforehand as multiple (carbon cycle) processes are involved.

2.2.1 Physical processes

An important feedback is the coupling of temperature to atmospheric pCO₂. There are several ways temperature effects the carbon cycle. For example, decreasing temperatures increase the solubility of CO₂, which results in more uptake of CO₂ by the ocean. For this feedback, we make a distinction between box 5 and boxes 1, 2, and 7. Box 5, the southern high-latitude surface box, is more isolated than the other boxes due to the Antarctic Circumpolar Circulation (ACC). Therefore, we have included the option in the model to use a different sensitivity in box 5. The temperature in the boxes is calculated as follows.

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{T}_i=T_{i,\mathrm{base}}+\mathrm{\Delta }{T}_i,i=\mathrm{1},\mathrm{2},\mathrm{5},\mathrm{7}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{T}_i={\mathit{\lambda }}_T\times \mathrm{0.54}\times \mathrm{5.35}\times \ln {\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}}{{\mathrm{CO}}_{{\mathrm{2}}_{\mathrm{base}}}}},i=\mathrm{1},\mathrm{2},\mathrm{7}\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{T}_{\mathrm{5}}={\mathit{\lambda }}_{T_{\mathrm{5}}}\times \mathrm{0.54}\times \mathrm{5.35}\times \ln {\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}}{{\mathrm{CO}}_{{\mathrm{2}}_{\mathrm{base}}}}}\end{array}$$

Here T_i,base is the base temperature in the SCP-M. The change in temperature is dependent on atmospheric pCO₂ and a base value of atmospheric pCO₂. This base value is the steady-state solution in the SCP-M without feedbacks (Table 1). Climate sensitivity can be changed via the λ parameters. For a λ of 1, sea surface temperatures increase 2 K per CO₂ doubling. As a reference, a 2 K warming for surface air temperatures is at the lower end of the range found in CMIP6 models (Zelinka et al., 2020).

Besides an effect on solubility, temperature can also affect the piston velocity. In the often used (Wanninkhof, 1992) formulation, the piston velocity is dependent on temperature via the Schmidt number (Eqs. 6 and 7). In our model, we use this dependency on the Schmidt number, which causes the piston velocity to increase for warmer temperatures. Hence

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{k}_{w,i}=& \left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{P}}\right)\times k_{w,i\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{P}}\times k_{w,i\mathrm{base}}\times {\left({\displaystyle \frac{S{c}_i}{\mathrm{660}}}\right)}^{-\mathrm{0.5}}\\ & i=\mathrm{1},\mathrm{2},\mathrm{5},\mathrm{7},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(7)}& \begin{array}{rl}S{c}_i& =\mathrm{2116.8}-\mathrm{136.25}{T}_i+\mathrm{4.7353}{T}_i^{\mathrm{2}}-\mathrm{0.092307}{T}_i^{\mathrm{3}}\\ & +\mathrm{0.0007555}{T}_i^{\mathrm{4}},i=\mathrm{1},\mathrm{2},\mathrm{5},\mathrm{7}.\end{array}\end{array}$$

In these equations, k_w is the used piston velocity, k_w,base is the piston velocity in the SCP-M (3 m d⁻¹), and T is the temperature of the box in degrees Celsius. The λ parameter needs to be either 0 (constant piston velocity, as in SCP-M) or 1 (variable piston velocity). Notice that if the temperature feedback is used (λ_T>0), the Schmidt number depends on atmospheric pCO₂.

2.2.2 Biological processes

A large limitation in the original SCP-M is the constant biological production. Nutrient availability introduces a large constraint on biological production, but this process is completely absent in the original SCP-M. This process is introduced in the model here by adopting the expression used in the Long-term Ocean-atmosphere-Sediment Carbon cycle Reservoir model (LOSCAR) (Zeebe, 2012). In LOSCAR, production is dependent on the upwelling of nutrients, which in our model translates to the following expressions.

$$\begin{array}{}\text{(8a)}& {\displaystyle }\begin{array}{rl}{Z}_{\mathrm{1}}=& \left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{BI}}\right)\times Z_{\mathrm{1},\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{BI}}\times \\ & \left({\mathit{\gamma }}_{\mathrm{2}}\times {\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}_{\mathrm{3}}+R_{{\mathrm{PO}}_{\mathrm{4}}}\right)\times {\mathit{ϵ}}_{\mathrm{1}}\end{array}\text{(8b)}& {\displaystyle }\begin{array}{rl}& Z_{\mathrm{2}}=\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{BI}}\right)\times Z_{\mathrm{2},\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{BI}}\times {\mathit{\psi }}_{\mathrm{2}}\times {\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}_{\mathrm{3}}\times {\mathit{ϵ}}_{\mathrm{2}}\end{array}\text{(8c)}& {\displaystyle }\begin{array}{rl}& Z_{\mathrm{5}}=\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{BI}}\right)\times Z_{\mathrm{5},\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{BI}}\times \mathit{\alpha }\times {\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}_{\mathrm{7}}\times {\mathit{ϵ}}_{\mathrm{5}}\end{array}\text{(8d)}& {\displaystyle }\begin{array}{rl}{Z}_{\mathrm{7}}& =\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{BI}}\right)\times Z_{\mathrm{7},\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{BI}}\times \left(\mathit{\alpha }\times {\mathit{\psi }}_{\mathrm{1}}+{\mathit{\psi }}_{\mathrm{2}}\right)\\ & \times {\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}_{\mathrm{4}}\times {\mathit{ϵ}}_{\mathrm{7}}\end{array}\end{array}$$

In these equations Z represents the production in the surface box, and Z_base the value used in the original SCP-M. Furthermore, α is the fraction of ψ₁ that moves from box 4 to box 7, and ϵ is the biological efficiency in the box. As with the piston velocity, λ_BI is either 0 (SCP-M) or 1. Notice that the current branch represented by ψ₁, which flows from box 4 to box 5, does not influence the production in box 5. We do not use this branch, since it is assumed to flow into box 5 below the euphotic zone.

In Eq. (8a–d) the biological efficiency (ϵ) is also introduced. There are studies (e.g., Cael et al., 2017) where they relate biological efficiency to temperature. We have adopted a simple linear relation to represent the influence of temperature on biological efficiency, i.e.,

$$\begin{array}{}\text{(9)}& {\mathit{ϵ}}_i={\mathit{\lambda }}_{\mathit{ϵ}}\times \left(-\mathrm{0.1}\mathrm{\Delta }{T}_i\right)+{\mathit{ϵ}}_{i,\mathrm{base}},i=\mathrm{1},\mathrm{2},\mathrm{5},\mathrm{7}.\end{array}$$

In this equation, λ_ϵ controls how strong the relation is between the efficiency and temperature change (ΔT). In addition, ϵ_base is the base value of the biological efficiency. These values have been fitted so that Z is equal to Z_base under the original parameter values in the SCP-M.

In the SCP-M, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ is the only nutrient. In the real ocean, additional nutrients play a role in biological production, one of them being nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$). During photosynthesis, organisms take up nitrate, and thereby increase Alk. This biological influence on Alk is not incorporated in the SCP-M but is present in many other models (e.g., Kwon and Primeau, 2008). We have included this influence as follows:

$$\begin{array}{}\text{(10)}& A_{\mathrm{Bio},i}={\mathit{\lambda }}_{\mathrm{BA}}\times \left(-{\displaystyle \frac{\mathrm{16}}{\mathrm{106}}}\right)\times C_{\mathrm{Bio},i},i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5},\mathrm{6},\mathrm{7}.\end{array}$$

In this equation A_Bio is the biological flux affecting Alk. This flux is related to the DIC biological flux (C_Bio) and the N:C Redfield ratio ($\frac{\mathrm{16}}{\mathrm{106}}$). For this relation, the λ_BA parameter can be 0 (not included, original SCP-M) or 1 (included).

Finally, we have also included a feedback for the rain ratio, which is the fraction of calcifiers in the total biological production. In the original SCP-M this is a constant value for all boxes. Calcifiers can be limited in growth when ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations are too low. Ridgwell et al. (2007) model this limitation via the saturation state of CaCO₃ as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{F}_{\mathrm{Ca},i}& =\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{F}}\right)\times F_{\mathrm{Ca},\mathrm{base}}+{\mathit{\lambda }}_{\mathrm{F}}\\ & \times \mathrm{0.022}{\left({\displaystyle \frac{[\mathrm{Ca}{]}_i{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}_i}{K_{{\mathrm{sp}}_i}}}-\mathrm{1}\right)}^{\mathrm{0.81}},i=\mathrm{1},\mathrm{2},\mathrm{5},\mathrm{7}.\end{array}\end{array}$$

Here, F_Ca is the used rain ratio, and F_Ca,base is the value used in the original SCP-M (0.07). The saturation state is determined via the concentrations of calcium ([Ca]), the carbonate ion concentration [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$], and an equilibrium constant K_sp. In this feedback, λ_F is either 0 (SCP-M) or 1. The rain ratio feedback is a negative feedback. When carbonate concentrations increase in the surface layer, the rain ratio increases, and therefore more calcium carbonate is removed from the surface layer, effectively lowering the carbonate concentration.

2.3 Parameter continuation methodology

The SCP-M, including our representations of the additional feedbacks, leads to a system of ordinary differential equations of the form

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}\mathit{u}}{\mathrm{d}t}}=\mathit{f}\left(\mathit{u}\right(t),\mathit{p}),\end{array}$$

where u is the state vector (containing all the dependent quantities in all boxes), f contains the right-hand side of the equations, and p is the parameter vector. Usually, such models are integrated in time from a certain initial condition, and the equilibrium behavior is determined for different values of the parameters. However, this is not very efficient to scan the parameter space, and, moreover, it is difficult to detect tipping behavior. A much more efficient approach is to determine the equilibrium solutions directly versus parameters, avoiding time integration, using continuation methods.

Here, we use the continuation and bifurcation software program AUTO to scan the parameter space and detect bifurcations efficiently (Doedel et al., 2007). The SCP-M is very suitable to be implemented into AUTO and to easily compute branches of equilibrium solutions, such as steady states of Eq. (12), versus parameters. The equations of the SCP-M turn out to have a singular Jacobian matrix (because carbon, alkalinity, and phosphate quantities are determined up to an additive constant), which requires adding integral conservation equations. We have added such integral conservation equations for carbon (DIC and atmospheric pCO₂), Alk, and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ to the model equations to replace the equations for box 4. The conservation law for ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ is straightforward and already present in the model equations. The constant influx of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ via the rivers is equal to the constant outflux via the sediments.

In the original SCP-M model, carbon and Alk are conserved in the ocean, atmosphere, continents, and sediments. However, the continental and sediment stocks are not explicitly represented in the version of the SCP-M we use. However, we can describe the change of total carbon and total Alk in the combined atmosphere and ocean stocks over time as

$$\begin{array}{}\text{(13a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{dTC}}{\mathrm{dt}}}& =C_{\mathrm{river}}\times V_{\mathrm{1}}+\sum _{n=\mathrm{1}}^{\mathrm{7}}\left(C_{\mathrm{carb},n}\times V_n\right)\\ & +\sum _{n=\mathrm{1}}^{\mathrm{7}}\left(C_{\mathrm{bio},n}\times V_n\right),\end{array}\text{(13b)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\mathrm{dTAlk}}{\mathrm{d}t}}=A_{\mathrm{river}}\times V_{\mathrm{1}}+\sum _{n=\mathrm{1}}^{\mathrm{7}}\left(A_{\mathrm{carb},n}\times V_n\right).\end{array}\end{array}$$

In these equations “TC” and “TAlk” are the total carbon and alkalinity in the system. As with ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, total carbon and Alk change due to influx via the rivers (C_river and A_river) and outflux via the sediments. The carbon outflux via the sediments is determined by the sum of carbonate (C_carb) and biological (C_bio) fluxes in the system. For Alk, the biological influence is absent. Model simulations with the original SCP-M have shown that the influence of the biological fluxes is negligible; i.e., all biologically produced organic matter is respired in the ocean itself. Therefore, this term can be set to zero in Eq. (13a). This makes Eq. (13b) proportional to Eq. (13a), and hence we include only the latter and use it to determine the change in total Alk in the model.

We also changed the carbonate chemistry in the model. The original SCP-M uses the algorithm of Follows et al. (2006), which solves the carbonate chemistry by using hydrogen ion concentrations from a previous time step. Therefore, the algorithm is inherently transient and, since we directly solve for steady-state solutions, not suitable. We therefore adopted a simple “textbook” carbonate chemistry based on carbonate alkalinity (Williams and Follows, 2011; Munhoven, 2013). This method approximates oceanic pCO₂ by assuming that Alk is equal to carbonate alkalinity ($A_{\mathrm{C}}=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}]$). A disadvantage of this method is that pH values are generally a bit higher (0.15–0.2) than using more complicated algorithms (Munhoven, 2013). These higher pH values are one of the reasons our atmospheric pCO₂ values are lower than in the original SCP-M (approximately 60 ppm for case P-CTL described in Sect. 3).

Eventually, by including Eq. (13a) and the overall conservation equations, the version of SCP-M used is a dynamical system with a state vector of dimension d=20. There is one equation for atmospheric pCO₂; there are six for DIC, Alk, and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ in the ocean; and there is one equation for the total carbon content. Except for the new carbonate chemistry, the necessary changes made to the SCP-M do not change the outcome of the model compared to the original model. When the original model is fitted with the same carbonate chemistry based on carbonate alkalinity, the AUTO implementation and the original code produce the same results.

3 Results

In Sect. 3.1 we present the general sensitivity of atmospheric pCO₂ to variations in the AMOC strength. We extend these results in Sect. 3.2 by adding a coupling between the AMOC strength and atmospheric pCO₂. Internal variability found in the model will be presented in Sect. 3.3. An overview of all cases considered is given in Table 2. Our control experiment uses the original model, which is tuned to accurately represent the pre-industrial and last glacial maximum carbon cycle. From this “realistic” model we investigate the sensitivity of the carbon cycle to specific carbon cycle feedbacks that can be found in more detailed models. By gradually increasing the number of feedbacks in the model, we can assess the effects of the (combined) feedbacks.

Table 2 Overview of the cases considered and their notation. The left column displays the used feedback. The other columns show the notation and what feedbacks are included in the specific case. The “x” in the notation is replaced with either P for the PI configuration or L for the LGM configuration. Bold rows indicate that this combination of feedbacks is also used for cases with a coupling between the AMOC and atmospheric pCO₂ (Sect. 3.2). For these cases, “C” is added to either “P” or “L=” to denote the coupling. The last column represents the feedback combinations used in Sect. 3.3. Case x-CTL is the original SCP-M.

3.1 Sensitivity of atmospheric pCO₂ to the AMOC

In this section the AMOC strength is used as a control parameter, and steady states are calculated versus this parameter. For each configuration (PI and LGM) we use three reference cases (x-CTL, the original SCP-M configuration; x-BIO, with a different parameterization for biological production; and x-ALL, with all feedbacks included, in Table 2, where x is either P for the PI or L for the LGM configuration). The steady-state value of atmospheric pCO₂ versus AMOC is shown for the reference cases in Fig. 2, where all branches represent stable fixed points. For the cases where the biological feedback is not included, the solutions for smaller values of AMOC ($<\sim \mathrm{12}$ Sv) display negative ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations in box 2 and hence are not allowed. Such boundaries can be automatically monitored in AUTO, and the continuation is stopped once a boundary is exceeded.

Figure 2 Atmospheric pCO₂ in parts per million under varying AMOC strength in sverdrups for three reference cases (blue: no feedback; red: with biological feedback; black: all feedbacks) in two configurations (solid: PI; dashed: LGM). All branches represent stable fixed points.

For the PI configuration, Fig. 2 shows that, whereas pCO₂ increases for larger AMOC strengths in the P-CTL case, it remains fairly constant in P-BIO and P-ALL. Atmospheric pCO₂ in the P-BIO case peaks around 5 Sv and then decreases until approximately 20 Sv, after which it increases slightly again. This different behavior occurs because, in the P-BIO case, the AMOC has competing influences on DIC concentrations of the surface ocean. A first effect of an increasing AMOC is to increase the ventilation of the deep ocean, which also increases DIC concentrations in the surface layer. This promotes outgassing to the atmosphere. However, by increasing the AMOC strength, biological production in boxes 2 and 7 is also increased. As a result, DIC and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ are transported from the surface layer to the deep ocean. The first effect is dominant after 20 Sv and the second effect in the range of 5 to 20 Sv. The absence of the second effect in P-CTL explains the difference in sensitivity between P-CTL and P-BIO. P-ALL behaves fairly similarly to P-BIO, except in the regime with a weak AMOC strength ($<\sim \mathrm{4}$ Sv). This behavior is caused by the saturation-dependent rain ratio.

When we look at the other cases (Fig. 3a), we see that they either behave qualitatively like P-CTL (the cases without “BIO”) or behave like P-BIO (cases with “BIO”). Looking in more detail, we can see that when we include the rain ratio feedback (cases P-RAIN, P-BIO-RAIN) atmospheric pCO₂ is higher, and when we include the biological influence on alkalinity, atmospheric pCO₂ is lower (cases P-BALK, P-BIO-BALK). The results in Fig. 3a show that the biological feedback (λ_BI=1) is the most dominant feedback in the PI configuration; i.e., including this feedback leads to a completely different sensitivity of the carbon cycle to changes in the AMOC strength.

Figure 3 Atmospheric pCO₂ in parts per million (color shading) under varying AMOC strength in sverdrups for the cases considered in Table 2. (a) Pre-industrial configuration. (b) Last Glacial Maximum configuration. Note that the range of the color shading differs between the two configurations and that some CO₂ concentrations fall outside the displayed range. The AMOC ranges of the bars differ, because for some cases the steady solution becomes nonphysical (e.g., negative concentrations or large subzero temperature). The vertical black lines represent the AMOC strength in the original SCP-M.

For the LGM configuration (Fig. 2), two important differences with respect to the PI configuration appear: (1) atmospheric pCO₂ is approximately 80 ppm lower, and (2) the L-BIO and L-ALL cases have a different sensitivity than the P-BIO and P-ALL cases for lower AMOC values. Where in P-BIO atmospheric pCO₂ decreases for an increasing AMOC between 5 and 20 Sv, L-BIO shows a monotonous increase in atmospheric pCO₂ from 3 Sv onward. We see this different relation because in the LGM configuration, deep-ocean ventilation, which can be seen as the sum of the GOC and AMOC, is lower due to a weaker GOC. Consequently, deep-ocean ventilation is more sensitive to changes in the AMOC. This eventually causes the different response of the L-BIO and L-ALL cases with respect to the P-BIO and P-ALL cases. The L-TEMP to L-BIO-TEMP(A)-EFF cases (Fig. 3b) relate to the L-CTL and L-BIO cases as in the PI configuration. Figure 3b shows that in the LGM configuration, as is the case in the PI configuration, the biological feedback is most dominant. The other feedbacks only influence the offset of CO₂ concentrations, but they do not result in large changes to the relation between the AMOC and atmospheric pCO₂.

3.2 Coupling AMOC–carbon cycle

The AMOC strength also depends on atmospheric pCO₂, and below we will discuss the steady-state model solutions when a coupling between the AMOC and atmospheric pCO₂ is applied. This coupling is based on how the AMOC responds to increasing atmospheric pCO₂ in CMIP6 models (e.g., Bakker et al., 2016) and given by

$$\begin{array}{}\text{(14)}& {\mathit{\psi }}_{\mathrm{2}}={\mathit{\psi }}_{\mathrm{2},\mathrm{base}}\times \left(\mathrm{1}-{\mathit{\lambda }}_A\times \mathrm{0.1}\times \mathrm{0.54}\times \mathrm{5.35}\times \ln {\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}}{{\mathrm{CO}}_{{\mathrm{2}}_{\mathrm{base}}}}}\right).\end{array}$$

In this equation ψ_2,base is a base value of the AMOC taken from the uncoupled case (where the AMOC is prescribed), ψ₂ is the actual AMOC strength in m³ s⁻¹ in the coupled case, and λ_A determines the strength of the coupling. We use three different values of λ_A in this section: (1) 0 (no coupling), 1 (20 % decrease for a CO₂ doubling), and 4 (80 % decrease for a CO₂ doubling). As the AMOC strength ψ₂ is now part of the state vector, we need other quantities as control parameters. We will use three different parameters here: (1) the rain ratio (F_Ca), (2) the biological production (Z), and (3) the piston velocity (k_w). We have chosen these three parameters since they (approximately) represent the three carbon pumps: the carbonate pump, the soft tissue pump, and the solubility pump We follow the steady-state solution in these parameters, where possible, between 0.1 and 10 times the reference value (indicated by the multiplier in Fig. 4). This large, though not necessarily realistic, range is used to test the sensitivity of the carbon cycle to the parameters, and to see whether bifurcations can arise in the carbon cycle.

Figure 4 Atmospheric pCO₂ in parts per million (color shading) under varying parameter values. Panels (a, c, e) represent cases of the PI configuration and (b, d, f) of the LGM configuration. Panels (a, b) show cases where the strength of the rain ratio is varied between 0.1 and 10 times the original value. Panels (c)–(f) show the same but for cases where biological production (c, d) and the piston velocity (e, f) are varied. In total seven feedback combinations are used, denoted by the text within the graph. For each case, three different coupling strengths have been used: (1) λ_A=0, (2) λ_A=1, and (3) λ_A=4.

When we look at the effect of increasing λ_A, i.e., the coupling, we see that the general sensitivity of the solution to changes in model parameters decreases. This effect is best seen in the LC-CTL case, but also present in the other cases, though less pronounced.

In Fig. 4a and b we plot the results when we use the rain ratio as a control parameter in the continuation. There are no large differences between the different cases and configurations. Generally, we see two regimes. For low rain ratios, the solution is quite sensitive to changes in the rain ratio. Where the coloring in Fig. 4a and b is yellow (around 230 ppm for the PI and around 140 ppm for the LGM configuration), we see a shift: the solution becomes less sensitive to changes in the rain ratio. To explain the regimes of sensitivity, we note that the CaCO₃ production is linearly related to the rain ratio. The production minus the dissolution of CaCO₃ in the water column determines the outflux of Alk and DIC via the sediments. The different regimes can be explained by the amount of CaCO₃ dissolution in the deep and abyssal ocean. For low rain ratios, the saturation state in the ocean is larger than 1, which means there is no saturation-driven dissolution and only constant dissolution. This makes the outflux of Alk and DIC linearly proportional to the production: if the rain ratio is low, the outflux is also low. Because we are looking at a steady-state solution, this decrease in burial has to be compensated for by a weaker influx, i.e., a lower river influx. This is only possible when atmospheric pCO₂ is lower. For larger rain ratios, we have both saturation-dependent and constant dissolution in the subsurface boxes, i.e., more dissolution in the water column. Due to the variable dissolution, the outflux of Alk and DIC is no longer fully determined by CaCO₃ production. This results in a lower sensitivity of the outflux to changes in the rain ratio. Therefore, atmospheric pCO₂ is also less sensitive to the rain ratio.

For biological production as a control parameter (Fig. 4c and d) again all cases show comparable behavior. We can see that the parameter range for higher biological production is short. This is because ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations become negative at this point, even when we include the biological feedback. All cases have a maximum in atmospheric pCO₂ around 0.7–0.8 times the original value. When the multiplier is lower than this value, we see a positive relation (higher biological production, higher atmospheric pCO₂). For values larger than the maximum, we see an opposite relation, i.e., lower atmospheric pCO₂ for higher biological production. This second regime is generally what we would expect when biological production is increased; i.e., when biological production removes more carbon from the surface layer, more carbon can be taken out of the atmosphere by the surface ocean, which reduces atmospheric pCO₂. The first regime is not what we would expect at first, but this can be explained by the same mechanism as for the rain ratio: reduced biological production leads to low production of CaCO₃, leading to low burial rates of CaCO₃. Lower burial rates lead to lower river influx because the sources and sinks of alkalinity to the ocean balance, which can only be achieved by decreasing atmospheric pCO₂. Again, increasing the AMOC coupling only reduces the sensitivity of the solutions.

Figure 5 (a) Bifurcation diagram for case L-HB in atmospheric pCO₂–AMOC space. Blue solid lines denote stable steady states (or fixed points, FPs), red dashed lines indicate unstable steady states, black solid lines indicate a stable limit cycle (LC), and the black square denotes the location of a (supercritical) Hopf bifurcation (HB). (b) The oscillation of atmospheric pCO₂ in parts per million versus time in years for the limit cycle at 15 Sv. The period is 5814 years.

In Fig. 4e and f, we plot the results when we use the piston velocity parameter (k_w) as a control parameter. By the gradually changing colors, we can see a logarithmic relation with higher sensitivities for lower piston velocities. The different feedbacks, configurations, and coupling strengths have the same effect as for the other two control parameters discussed above.

3.3 Internal oscillation

The feedback strengths we have used so far have been quite modest. The continuation methodology enables us to efficiently look at cases with different feedback strengths and to see whether different combinations can induce bifurcations in the carbon cycle, and by extensively scanning the parameter space we found such bifurcations. Especially in the LGM configuration, when climate sensitivity (λ_T) and the biological efficiency feedback (λ_ϵ) are increased, bifurcations arise on the branches of steady solutions. With case L-HB (for parameter values, see Table 2), we present an example where we find a supercritical Hopf bifurcation (HB) around 13 Sv (Fig. 5a) in the uncoupled case (λ_A=0, so the AMOC strength is a control parameter again). The HB produces a stable limit cycle extending to larger AMOC strengths with a period between 5000 and 6000 years where all state variables oscillate. In this section we look at the internal oscillation at 15 Sv (Fig. 5b). The oscillation has a period of 5814 years, and atmospheric pCO₂ has a range of 72 ppm.

The HB described in this section exists for a large range of parameter values and is thus robust. One important constraint on the existence of the bifurcation is the coupling strength between atmospheric pCO₂ and biological production. This coupling comes down to the effect of atmospheric pCO₂ on the biological efficiency (ϵ), which can be increased by increasing the temperature feedback (λ_T) and/or the efficiency feedback (λ_ϵ). We do not find this bifurcation in the PI configuration, because when the biology feedback (λ_BI=1) is included, atmospheric pCO₂ is insensitive to changes in the AMOC strength (P-BIO case, Fig. 2). Because of this low sensitivity, surface ocean temperature and biological efficiency are also insensitive to changes in the AMOC strength in the PI configuration. Therefore, the coupling between the two is less effective in this configuration, and we do not find a HB.

To explain the mechanism behind the oscillation, we have to look at the time-dependent solution of the model. What is important for this oscillation is the coupling between atmospheric pCO₂ and the alkalinity cycle. Alkalinity influences the gas exchange between the ocean and the atmosphere via the carbonate chemistry and is, in turn, influenced by atmospheric pCO₂ because the source and sink of alkalinity are coupled to pCO₂. The source, the river influx, is directly proportional to atmospheric pCO₂ according to Eq. (1).

The sink, i.e., outfluxing via the sediments, is related to CaCO₃ burial, which is the difference between CaCO₃ production in the surface ocean and CaCO₃ dissolution in the ocean and sediments (Eq. 2). In the oscillation, the saturation state of CaCO₃ in the ocean is larger than 1 everywhere. This happens when the river influx is larger than the biogenic flux in the surface ocean (Zeebe and Westbroek, 2003), which is plausible for the past oceans. Therefore, total dissolution in the ocean is constant and does not vary. This means that CaCO₃ burial becomes a function of CaCO₃ formation and thus of biological production. Since this production is dependent on the biological efficiency, which is directly proportionate to atmospheric pCO₂, the sink is also influenced by atmospheric pCO₂. However, the effect of atmospheric pCO₂ on the source and sink is opposite. When atmospheric pCO₂ is high, the river influx is high, while the sediment outflux is low. This is key to the general mechanism sketched in Fig. 6.

Figure 6 Schematic representation of the mechanism of the internal oscillation. The rectangles represent state variables, while the pointed blocks represent fluxes or model parameters. Some boxes have thicker outlining. These boxes cause a chain of events. The chains corresponding to the box are the boxes with the same color shading. The gray and blue rectangles in the background represent a quarter period. In the second half of the period, processes are replaced with a dashed line. These processes are the opposite of what happens in the first half of the period.

The results show that atmospheric pCO₂ is affected by the amount of ingassing into box 1. Therefore, we start the explanation of the oscillation in Fig. 6 at this point. At the beginning of the oscillation (time t=0 in Fig. 6), ingassing in box 1 starts to decrease. As a result, atmospheric pCO₂ starts to increase approximately 200 years later. There is a delay, since atmospheric pCO₂ is not solely determined by the gas exchange with box 1. The increase in atmospheric pCO₂ has multiple effects. First of all, temperatures start to increase, which lowers biological efficiency. This in turn reduces CaCO₃ production, and thus the sink of alkalinity is also reduced. Another effect of increasing atmospheric pCO₂ is an increasing river flux, i.e., an increasing source of alkalinity into the ocean. After a quarter period (time $t=T/\mathrm{4}$ in Fig. 6), the source becomes larger than the sink, and total alkalinity in the ocean starts to increase. Meanwhile, atmospheric pCO₂ is still increasing. As a result, the river influx also keeps increasing, while the sediment outflux keeps decreasing. After half a period (time $t=T/\mathrm{2}$ in Fig. 6), oceanic pCO₂ in box 1 starts to decrease because alkalinity concentrations in box 1 have increased. The lower oceanic pCO₂ causes ingassing into box 1 to increase, which in turn decreases atmospheric pCO₂. The other half of the period is as explained above, but then the opposite.

The processes described above are important for driving the oscillation. However, these are not the only processes represented in the model. The concentrations of DIC, Alk, and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ in the ocean boxes are subtle balances of multiple larger fluxes where the sum of these fluxes can be more than 100 times smaller than the individual fluxes. It is therefore difficult to assess the effects of all the individual fluxes, since they also depend on each other. We do see that the DIC concentrations in the surface ocean boxes lag atmospheric pCO₂ by multiple centuries. The lag in DIC concentrations thereby increases oceanic pCO₂ after atmospheric pCO₂ has reached its maximum, which dampens the amplitude of the oscillation. The solubility constant (K₀) and dissociation constants (K₁ and K₂), which are also important for the air–sea gas exchange, oscillate due to the dependency on temperature and also dampen the amplitude of the oscillation. It is good to note that all these processes are responsible for the exact shape and amplitude of the oscillation. However, the coupling between atmospheric pCO₂ and the alkalinity cycle appears to be the driving mechanism.

Figure 7 (a) Atmospheric pCO₂ (red, left y axis) in parts per million, and total alkalinity (blue, right y axis) in the ocean in petamoles for one oscillation. Total alkalinity lags pCO₂ by approximately a quarter period. (b) The source (blue) and the absolute value of the sink (red) of alkalinity in the ocean. The source represents river influx, and the sink represents the sediment outflux. When the lines cross, i.e., around 1500 and 4400 years, total alkalinity in (a) has a minimum and a maximum, respectively.

In Fig. 7a, we can see that total alkalinity in the ocean lags atmospheric pCO₂ by approximately a quarter period. In Fig. 7b we can also see the anti-correlation between the source and sink of alkalinity to the ocean. Comparing the sink and source, we can clearly see a strong (anti-) correlation between atmospheric pCO₂ and the (sink) source of alkalinity. The anti-correlation between the source and sink is the driving mechanism behind the oscillatory behavior. It is good to point out that the amplitude of the sink of alkalinity is larger than that of the source. The timescale of the oscillation (∼6000 years) is related to the adjustment time of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ to an imbalance between the influx and outflux of alkalinity and DIC in the ocean. This process, termed the calcium carbonate homeostat (Sarmiento and Gruber, 2006), has a timescale between 5 and 10 kyr (Archer et al., 1997). The period of the internal oscillation corresponds well to this range. The river influx, which plays a role in the oscillation, is usually viewed as a slow process because of the long timescales of silicate weathering (10 kyr or more). In the oscillation, however, the river flux varies on shorter timescales. This is because in the model, carbonate weathering is more important than silicate weathering and acts on shorter timescales (1 to 10 kyr timescales). Furthermore, the system does not reach equilibrium and continuously oscillates, whereby the river flux responds directly to the oscillations in atmospheric pCO₂. It is also good to note that, even though it seems box 1 is a main driver in the oscillation, it is in fact a global process due to the role of CaCO₃ burial; the amplitude of CaCO₃ burial is more than 2 times larger than the amplitude of the river flux.

4 Summary and discussion

In this study we investigated steady states of an extended version of the simple carbon cycle box model (SCP-M), where additional feedbacks have been included. Focus was on the relation between the AMOC and atmospheric pCO₂ for these steady states, with a special attention to the effect of feedbacks and climatic boundary conditions on this relation. Although the model we use is a simple box model, the original SCP-M was shown to be quite capable of simulating present-day observations and proxy data (LGM) (O'Neill et al., 2019).

In Sect. 3.1 we looked into how the carbon cycle, and specifically atmospheric pCO₂, responds to changes in the AMOC. These cases include different combinations of additional feedbacks. Our results (Sect. 3.1) suggest that the most important feedback is the biological feedback, represented by Eq. (8a–d). In both the PI and the LGM configurations, this feedback leads to a different sensitivity of atmospheric pCO₂ to the AMOC (Fig. 2). Other feedbacks did not introduce large effects on the sensitivity (Fig. 3). This shows that biology can exert a large effect on atmospheric pCO₂, which supports studies with more detailed models where biological production plays a role in the response of atmospheric pCO₂ to changes in the AMOC (e.g., Nielsen et al., 2019). The results also show the importance of the climatic boundary conditions, as was already stated in Gottschalk et al. (2019). Generally, cases with the biological feedback (x-BIO, and other cases including “BIO”) respond differently in the LGM configuration than in the PI configuration. This is related to the difference in deep-ocean ventilation between the two configurations.

When a coupling between the AMOC and atmospheric pCO₂ is included (Section 3.2), the pCO₂ of the steady solutions becomes less sensitive to changes in model parameters (k_w, Z, F_Ca). This shows that the coupling works as a negative feedback in the carbon cycle dynamics. What is interesting to see is that the carbon cycle feedbacks do not have a large effect on the AMOC–pCO₂ relation. This implies that ocean circulation is very effective in damping changes in gas exchange (k_w) and biological (Z) and CaCO₃ (F_Ca) production.

When considering bifurcations of the steady solutions, an important result is what we did not find: saddle-node bifurcations. Hence, although quite nonlinear carbon cycle processes have been captured in this model, no multiple equilibrium regimes and associated hysteresis occur. As a consequence, any sharp transition in carbon cycle quantities cannot be easily linked to a transition between different steady states. However, we did find internal oscillations in the model, in particular with a period of 5000 to 6000 years related to the CaCO₃ homeostat (Fig. 6). Important for this oscillation is the process representation that CaCO₃ production reduces for increasing temperatures, which is supported by studies that suggest a decreased production under high-atmospheric-CO₂ concentrations (Barker and Elderfield, 2002). However, this assumption is under debate as there are studies that find an increased calcifier production for higher temperatures (Cole et al., 2018) in specific situations. Whether this internal oscillation also exists in a system where the AMOC strength and atmospheric pCO₂ are coupled (as in Sect. 3.2) is not known. The internal oscillations were found using the AMOC strength as a control parameter, which is not possible with a relation as in Sect. 3.2.

Linking this oscillation to proxy data is difficult, especially since the variation in atmospheric pCO₂ is relatively high (72 ppm) for reasonable AMOC values. If we look for example at the record of the last glacial period, pCO₂ variations are of the order of 20 ppm (Bauska et al., 2021). The variation found in our model is closer to that during the Pleistocene glacial cycles, but on a much shorter timescale. The timescale is actually closer to that of the Heinrich events. It is therefore hard to find an oscillation like this in the past record, but this does not mean the mechanism is not relevant. If we look at more fundamental work, our mechanism shares similarities to the internal oscillation found in a conceptual model where only Alk and DIC are resolved (Rothman, 2019). The mechanism in Rothman (2019) is based on the imbalance between the influx and outflux of DIC in the surface ocean and is thus comparable to our mechanism. The phase differences in our model between quantities in the carbonate system (i.e., DIC, Alk, pH, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and H₂CO₃) in the top 250 m compare well to those in Rothman (2019) (not shown). However, the responsible processes are different. In Rothman (2019) there is an important role for respiration of organic matter. In our model, this flux is implicitly modeled, and we can reconstruct a similar flux from the export production. This reconstructed flux has comparable phase differences with the carbonate content as in Rothman (2019), but the relative strength of the flux does not match the burial flux in our model. This means that the SCP-M captures a different internal oscillation. In Rothman (2019) there is an important role for the ballast feedback because it couples the sources and sinks of DIC using the carbonate-ion concentration. In our oscillation, it is not the ballast feedback that drives the oscillation, but the CaCO₃ homeostat, coupling the sources and sinks of alkalinity through atmospheric pCO₂.

The results in this study are achieved with a very simple framework with multiple assumptions and limitations. The main assumption we make is that the SCP-M is a well-performing model for the Last Glacial Maximum and pre-industrial periods. Comparison of the model results with observations in O'Neill et al. (2019) supports this assumption. Assumptions in the river flux parameterization that possibly affect the oscillation are the parameter values and the fact that there is no delay between the river influx and atmospheric pCO₂. The parameter values are important for the amplitude of the oscillation, and decreasing the parameters would result in a decrease in amplitude of the oscillation. However, the assumed parameters are fitted to represent estimated river influx in present-day conditions. The river flux parameterization assumes that changes in atmospheric pCO₂ result in changes in the river flux into the ocean mainly due to carbonate weathering. These changes are not relevant on short (annual to centennial) timescales but will affect model results on (multi) millennial timescales. In our results, this parameterization is relevant since the approximate quarter period delay between atmospheric pCO₂ and total alkalinity (and alkalinity in box 1) is important in the oscillation mechanism. The river influx plays a role in this by changing the alkalinity in box 1. It would be interesting to introduce a time delay between atmospheric pCO₂, continental weathering, and the river flux. The inclusion of such a time delay would make the carbon system more complicated, and moreover, additional oscillatory behavior is expected. However, this extension is outside the scope of this paper and will not affect our results regarding the existence of oscillatory behavior in the carbon cycle.

We also made assumptions for the features we added to the model. The first assumption is that refitting of the parameters is not necessary. For most changes we made, we do expect this assumption to be valid since for most features the elemental cycles remained the same, and constant parameter values were replaced by equations which keep the parameter values close to their original values. The addition of the biological alkalinity flux might make refitting of parameters necessary since a completely new process is added to the alkalinity cycle. Refitting of the parameters would be a large exercise and would also make comparison between the different cases difficult. However, cases with this feedback do not show divergent results compared to other cases. Maybe the most impactful change we made is the simplification of the carbonate chemistry. This change typically reduces pH by 0.15–0.2 (Munhoven, 2013) and changes equilibrium pCO₂ values by 20% (Munhoven, 2013), explaining the approximately 60 ppm lower atmospheric pCO₂ in our model. The assumption that biological efficiency is linearly related to change in temperature might not be valid while this assumption is important for the driving mechanism of the internal oscillation. However, what seems to be important for the oscillation is that the coupling between atmospheric pCO₂ and the biological efficiency is strong enough and not necessarily the exact formulation of the feedback. Limitations of the model include the incapability to discern between ocean boxes and strict, slightly arbitrary box boundaries. Due to these limitations, this model is not suitable for looking at regional processes. However, the original SCP-M simulates representative global values, making it suitable for the application in this study.

In this study we have scanned large ranges of parameter values, with some values outside realistic ranges. The parameters we have varied are AMOC strength, the rain ratio, biological production, the piston velocity, and climate sensitivity. By using such a wide range for certain parameters, we can be quite certain that there are no saddle-node bifurcations and therefore no multiple equilibria in realistic parameter ranges. We believe that most results are within a realistic range for the AMOC strength since present-day model simulations show maximum AMOC strengths of around 25 Sv (Weijer et al., 2019), while model simulations simulating an AMOC collapse show very weak AMOC strengths. In Sect. 3.2, we studied a large range of rain ratio, biological production and piston velocity values. The main purpose of this large range was to see whether bifurcations would occur, which did not. The climate sensitivity variations we used are all within CMIP6 ranges (1.8–5.6 K; Zelinka et al., 2020). Cases without the temperature feedback, however, do yield unrealistic results since ocean temperatures remain above freezing temperatures even for near-zero atmospheric pCO₂ values.

In conclusion, we have found that the relation between atmospheric pCO₂ and the AMOC strength relies mostly on biological processes and climatic boundary conditions. Therefore, we suggest that by comparing results of different models special attention should be given to the way biological production is represented. Our study also shows that atmospheric pCO₂ appears to be rather insensitive to changes in the AMOC strength, which suggests that projected weakening of the AMOC in the future does not lead to a large response in atmospheric pCO₂. In this study we also searched for saddle-node bifurcations, but we did not find any, suggesting that tipping points in the carbon cycle are unlikely to occur. Our most interesting result is the discovery of an internal oscillation in the carbon cycle, and we hope that the mechanism behind this oscillation will stimulate further model work and be useful for explaining past atmospheric pCO₂ variability.

Appendix A: Model parameters

In this appendix values and descriptions of the parameters in the extended SCP-M are given. In Tables A1 to A3 the parameter values used in our model are presented. The values presented here are for the pre-industrial configuration. The parameter values that are different in the Last Glacial Maximum configuration are presented in Table 1. All parameter values, except the biological efficiency (ϵ) parameters, are taken from the SCP-M. The biological efficiency parameters have been fitted such that Z in Eq. (8a–d) is equal to Z_base under the original parameter values in the SCP-M. When determining the value of ϵ_2,base, we also took the effect of the biological production in the North Pacific into account, which leads to a value of ${\mathit{ϵ}}_{\mathrm{2},\mathrm{base}}>\mathrm{1}$. In Table A4 we also present the literature where the expressions for the equilibrium constants were taken from.

Table A1 Symbol (column 1), description (column 2), value (column 3), and units (column 4) of the general parameters used in our model.

Table A2 Symbol (column 1), description (column 2), value (column 3), and units (column 4) of parameters concerning the dimensions of the boxes used in our model.

Table A3 Symbol (column 1), description (column 2), value (column 3), and units (column 4) of the other parameters used in our model.

Weiss (1974)Lueker et al. (2000)Lueker et al. (2000)Mucci (1983)Millero (1983)

Table A4 The symbols and description of the equilibrium constants are presented in the first two columns. The third column presents the source of the used expression.