The Bern Simple Climate Model (BernSCM) v1.0: an extensible and fully documented open-source re-implementation of the Bern reduced-form model for global carbon cycle–climate simulations

The Bern Simple Climate Model (BernSCM) is a free open-source re-implementation of a reduced-form carbon cycle–climate model which has been used widely in previous scientific work and IPCC assessments. BernSCM represents the carbon cycle and climate system with a small set of equations for the heat and carbon budget, the parametrization of major nonlinearities, and the substitution of complex component systems with impulse response functions (IRFs). The IRF approach allows cost-efficient yet accurate substitution of detailed parent models of climate system components with near-linear behavior. Illustrative simulations of scenarios from previous multimodel studies show that BernSCM is broadly representative of the range of the climate–carbon cycle response simulated by more complex and detailed models. Model code (in Fortran) was written from scratch with transparency and extensibility in mind, and is provided open source. BernSCM makes scientifically sound carbon cycle–climate modeling available for many applications. Supporting up to decadal time steps with high accuracy, it is suitable for studies with high computational load and for coupling with integrated assessment models (IAMs), for example. Further applications include climate risk assessment in a business, public, or educational context and the estimation of CO₂ and climate benefits of emission mitigation options.

1 Introduction

Simple climate models (SCMs) consist of a small number of equations, which describe the climate system in a spatially and temporally highly aggregated form. SCMs have been used since the pioneering days of computational climate science to analyze the planetary heat balance (Budyko, 1969; Sellers, 1969) and to clarify the role of the ocean and land compartments in the climate response to anthropogenic forcing through carbon and heat uptake (e.g., Oeschger et al., 1975; Siegenthaler and Oeschger, 1984; Hansen et al., 1984). Due to their modest computational demands, SCMs enabled pioneering research using the limited computational resources of the time and continue to play a useful role in the hierarchy of climate models today.

Recent applications of SCMs are often found in research in which computational resources are still limiting. Examples include probabilistic or optimization studies involving a large number of simulations, or the use of a climate component as part of a detailed interdisciplinary model. SCMs are also much easier to understand and handle than large climate models, which makes them useful as practical tools that can be used by non-climate experts for applications for which detailed spatiotemporal physical modeling is not essential. This applies to interdisciplinary research, educational applications, or the quantification of the impact of emission reductions on climate change.

Figure 1 BernSCM as a box-type model of the carbon cycle–climate system based on impulse response functions. Heat and carbon taken up by the mixed ocean surface layer and the land biosphere, respectively, is allocated to a series of boxes with characteristic timescales for surface-to-deep ocean transport (τ_O) and of terrestrial carbon overturning (τ_L). The total perturbations in land and surface ocean carbon inventory and in surface temperature are the sums over the corresponding individual perturbations in each box ($m_{\mathrm{S}k},\mathrm{\Delta }{T}_k,m_{\mathrm{L}k}$). Using pattern scaling, the response in SAT can be translated to regional climate change for fields v(x,t) of variables such as SAT or precipitation.

An important application of SCMs is in integrated assessment models (IAMs). IAMs are interdisciplinary models that couple a climate component with an energy-economy model to simulate emissions and their climate consequences. Another application of simple models (e.g., Boucher and Reddy, 2008; Bruckner et al., 2003; Enting et al., 1994; Good et al., 2011; Hooss et al., 2001; Huntingford et al., 2010; Joos and Bruno, 1996; Oeschger et al., 1975; Raupach, 2013; Siegenthaler and Oeschger, 1984; Smith et al., 2017; Tanaka et al., 2007; Urban and Keller, 2010; Wigley and Raper, 1992) is to compare, analyze, or emulate more complex models (Geoffroy et al., 2012a, b; Meinshausen et al., 2011; Raper et al., 2001; Thompson and Randerson, 1999). Simple models also play a significant role in previous assessments of the Intergovernmental Panel on Climate Change (e.g., Harvey et al., 1997). The comprehensive scope and interdisciplinarity of such models raise the challenge of maintaining a high and balanced scientific standard across all model components, especially when human resources are limited. This may apply particularly to the climate component, as IAMs are mostly used within the economic and engineering disciplines. Climate and carbon cycle representation are central parts of an IAM and have been critically assessed in the literature (Joos et al., 1999a; Schultz and Kasting, 1997; van Vuuren et al., 2011).

BernSCM is a zero-dimensional global carbon cycle–climate model built around impulse-response representations of the ocean and land compartments, as described previously in Joos et al. (1996) and Meyer et al. (1999). The linear response of more complex ocean and land biosphere models with detailed process descriptions is captured using impulse-response functions (IRFs). These IRF-based substitute models are combined with nonlinear parametrizations of carbon uptake by the surface ocean and the terrestrial biosphere as a function of atmospheric CO₂ concentration and global mean surface temperature. Pulse response models have been shown to accurately emulate spatially resolved, complex models (Joos et al., 1996; Joos and Bruno, 1996; Meyer et al., 1999; Joos et al., 2001; Hooss et al., 2001).

BernSCM (Fig. 1) is designed to compute decadal- to millennial-scale perturbations in atmospheric CO₂, in climate and in fluxes of carbon and heat relative to a reference state, typically preindustrial conditions. The uptake of excess, anthropogenic carbon from the atmosphere is described as a purely physicochemical process (Prentice et al., 2001). As in pioneering modeling approaches with box-type (Oeschger et al., 1975; Revelle and Suess, 1957) and general ocean circulation models (Maier-Reimer and Hasselmann, 1987; Sarmiento et al., 1992), modification of the natural carbon cycle through potential changes in circulation and the marine biological cycle (Heinze et al., 2015) are not explicitly considered. While such modifications and their potential socioeconomic consequences are vividly discussed in the literature (Gattuso et al., 2015), associated climate–CO₂ feedbacks are likely of secondary importance. Estimated uncertainties in the marine carbon uptake due to climate change, including warming-driven changes in CO₂ solubility, are found to be smaller in magnitude than uncertainties arising from imperfect knowledge of surface-to-deep physical transport (see Fig. 2d and e in Friedlingstein et al., 2006). The exchange of CO₂ between the atmosphere and the surface ocean is described by two-way fluxes, from the atmosphere to the surface ocean and vice versa, and the net flux of CO₂ into the ocean is proportional to the air–sea partial pressure difference. CO₂ reacts with water to form carbon and bicarbonate ions (Dickson et al., 2007; Orr et al., 2015), and acid–base equilibria are described here using the well-established Revelle factor formalism (Siegenthaler and Joos, 1992; Zeebe and Wolf-Gladrow, 2001). The first-order climate–carbon feedback of a decreasing solubility in warming water is considered. Surface-to-deep exchange, the rate limiting step of ocean carbon and heat uptake, is described using an IRF. On timescales of up to a few millennia, processes associated with ocean sediments and weathering can be neglected. In such a closed ocean–atmosphere–land biosphere system, excess CO₂ is partitioned between the ocean and the atmosphere and a substantial fraction of the emitted CO₂ remains in the atmosphere and in the surface ocean in a new equilibrium (Joos et al., 2013). This corresponds to a constant term (infinitely long removal timescale) in the IRF representing surface-to-deep mixing. On multimillennial timescales, excess anthropogenic CO₂ is removed from the ocean–atmosphere–land system by ocean–sediment interactions and changes in the weathering cycle (Archer et al., 1999; Lord et al., 2016), and the IRF is readily adjusted to account for these processes, important for simulations extending over many millennia.

BernSCM simulates global mean surface temperature and the heat uptake by the planet. The latter is equivalent to the net top-of-the-atmosphere energy flux. Changes in the Earth's heat storage in response to anthropogenic forcing are dominated by warming of the surface ocean and the interior ocean (T. F. Stocker et al., 2013) due to their large heat capacity in comparison with that of the atmosphere and their large thermal conductivity in comparison to that of the land surface. Consequently, the atmospheric and land surface heat capacity is formally lumped with the heat capacity of the surface ocean in the BernSCM. The uptake of heat by the ocean (or planet) is, as for carbon, formulated as a two-way exchange flux. The flux of heat from the atmosphere into the surface ocean is taken to be proportional to the radiative forcing resulting from changes in CO₂ and other agents (Etminan et al., 2016). The upward loss of heat from the surface is proportional to the product of the simulated surface temperature perturbation and the (prescribed) climate sensitivity λ (Siegenthaler and Oeschger, 1984; Winton et al., 2010).

As with carbon, surface-to-deep transport is the rate-limiting step for ocean heat uptake and thus for the adjustment of surface temperature to radiative forcing. This transport is key to determine the lag between realized warming and equilibrium warming (Frölicher and Paynter, 2015). Again, this transport is described using an IRF. This IRF encapsulates the finite volume of the entire ocean. It also represents the range of transport timescales associated with advection, diffusion, and convection ranging from decades for the ventilation of thermocline to more than a millennium for deep Pacific ventilation as evidenced by transient tracers such as chlorofluorocarbons and radiocarbon (Olsen et al., 2016). The simulated surface ocean temperature perturbation, taken as a measure of global mean surface air temperature (SAT) change, may be combined with spatial patterns of change in temperature, precipitation, or any other variable of interest to compute regionally explicit changes (Hooss et al., 2001; Joos et al., 2001; B. D. Stocker et al., 2013) (Fig. 1).

Non-CO₂ radiative forcing may be prescribed, e.g., following estimates from complex climate–chemistry models (Myhre et al., 2013) or from simple emission-driven non-CO₂ modules of radiative forcing related to CO₂ chemistry (Joos et al., 2001; Smith et al., 2017) and reconstructions of solar and volcanic forcing (Eby et al., 2012; Jungclaus et al., 2017) and considering the forcing efficacy of non-CO₂ agents relative to CO₂ forcing (Hansen et al., 2005). Climate sensitivity characterizing the response to radiative forcing is a free parameter in BernSCM. Climate sensitivity may change under increasing warming, particularly in high-emission scenarios (Geoffroy et al., 2012a; Gregory et al., 2015; Pfister and Stocker, 2017). Here, climate sensitivity is assumed to be time invariant and a potential state dependency of climate sensitivity is not considered. This may be changed when more solid information on state dependency becomes available or for the purpose of sensitivity analyses. Similarly, ocean heat uptake efficacy (Winton et al., 2010), influencing the atmospheric temperature response to ocean heat uptake forcing, is set to 1 here.

The present version 1.0 of BernSCM is fundamentally analogous to the Bern model as used already in the IPCC Second Assessment Report, Bern-SAR (whereas different versions of the Bern model family were used in the more recent IPCC reports). BernSCM represents the relevant processes more completely than Bern-SAR, thanks to additional alternative representations of the land and ocean components, which contain a more complete set of relevant sensitivities to temperature and atmospheric CO₂.

Here, BernSCM model simulations are compared to previous multimodel studies. The model is run for an idealized atmospheric pulse CO₂ emission experiment of Joos et al. (2013), for an idealized CO₂ forcing experiment similar to simulations from the Climate Model Intercomparison Project 5 (CMIP5), and for the SRES A2 emission scenario used in the C4MIP study (Friedlingstein et al., 2006).

Together with this publication, BernSCM v1.0 is provided as an open-source Fortran code for free use. The code was also rewritten from scratch, with flexibility and transparency in mind. The model is comprehensively documented, and easily extensible. New alternative model components can be added using the existing ones as a template. A range of numerical solution schemes is implemented. Up to decadal time steps are supported with high accuracy, suitable for the coupling with emission models of coarse time resolution, for example. However, the published code is a ready-to-run stand-alone model, which may also be useful in its own right.

BernSCM offers a physically sound carbon cycle–climate representation, but it is small enough for use in IAMs and other computationally tasking applications. In particular, the support of long time steps is ideally suited to the application of BernSCM as an IAM component, as these complex models often use time steps on the order of 10 years.

BernSCM also offers a tool to realistically assess the climate impact of carbon emissions or emission reductions and sinks, for example in aviation, forestry (Landry et al., 2016), blue carbon management, peat development (Mathijssen et al., 2017), life cycle assessments (Levasseur et al., 2016), or to assess the interaction of climate engineering interventions such as terrestrial carbon dioxide removal with the natural carbon cycle (Heck et al., 2016).

In this paper, we describe the model equations (Sect. 2 and Appendix), illustrative simulations in comparison with previous multimodel studies, and uncertainty assessment (Sect. 3), followed by a discussion (Sect. 4) and conclusions (Sect. 5).

2 The BernSCM model framework and equations

BernSCM simulates the relation among CO₂ emissions, atmospheric CO₂, radiative forcing (RF), and global mean SAT by budgeting carbon and heat fluxes globally among the atmosphere, the (abiotic) ocean, and the land biosphere compartments. Given CO₂ emissions and non-CO₂ RF, the model solves for atmospheric CO₂ and SAT (e.g., in the examples of Sect. 3), but can also solve for carbon emissions (or residual uptake) when atmospheric CO₂ (or SAT and non-CO₂ RF) is prescribed, or for RF when SAT is prescribed.

The transport of carbon and heat to the deep ocean, as well as the decay of land carbon, results from complex but linear to first order behavior of the ocean and land compartments. These are represented in BernSCM using IRFs (Green's function). The IRF describes the evolution of a system variable after an initial perturbation, e.g., the pulse-like addition of carbon to a reservoir. It fully captures linear dynamics without representing the underlying physical processes (Joos et al., 1996). More illustratively, the ocean and land models can be considered to consist of systems of uncoupled first-order ordinary differential equations or “box models”, which are an equivalent representation of the IRF model components (Fig. 1).

The net primary production (NPP) of the land biosphere and the surface ocean carbon uptake depend on atmospheric CO₂ and surface temperature in a nonlinear way. These essential nonlinearities are described by parametrizations linking the linear model components.

2.1 Carbon cycle component

The budget equation for atmospheric carbon is

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{A}}}{\mathrm{d}t}}=e-f_{\mathrm{O}}-{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{L}}}{\mathrm{d}t}},\end{array}$$

where m_A denotes the atmospheric carbon stored in CO₂, e CO₂ emissions, f_O the flux to the ocean, m_L the land biosphere carbon stock, and t time. Here, m_L refers to the (potential) natural biosphere. Human impacts on the land biosphere exchange including land use and land use changes are not simulated in the present version and are treated as exogenous emissions (e). These emissions may be prescribed based on results from spatially explicit terrestrial models. An overview of the model variables and parameters is given in Tables A1 and A2.

The change in land carbon is given by the balance of NPP and decay of assimilated terrestrial carbon,

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{L}}}{\mathrm{d}t}}=f_{\text{NPP}}-f_{\text{decay}}.\end{array}$$

Decay includes heterotrophic respiration (RH), fire, and other disturbances due to natural processes.

Carbon is taken up by the ocean through the air–sea interface (f_O) and distributed to the mixed surface layer (m_S) and the deep ocean interior (f_deep):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{f}_{\mathrm{O}}={\displaystyle \frac{\mathrm{d}{m}_{\mathrm{S}}}{\mathrm{d}t}}+f_{\text{deep}}.\end{array}$$

Global NPP (f_NPP) is assumed to be a function of the partial pressure of atmospheric CO₂ ($p^{{\mathrm{CO}}_{\mathrm{2}}}$) and the SAT deviation from preindustrial equilibrium (functions for the implemented land components are given in Appendix A),

The net flux of carbon into the ocean is proportional to the gas transfer velocity (k_g) and the CO₂ partial pressure difference between surface air and seawater:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{f}_{\mathrm{O}}=k_{\mathrm{g}}{A}_{\mathrm{O}}\mathit{\epsilon }(p_{\mathrm{A}}^{{\mathrm{CO}}_{\mathrm{2}}}-p_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}),\end{array}$$

where A_O is ocean surface area and ε the atmospheric mass of C per mixing ratio of CO₂.

The global average perturbation in surface water $\mathrm{\Delta }{p}_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}$ is a function of dissolved inorganic carbon change (ΔDIC) in the surface ocean at constant alkalinity (Joos et al., 1996) and SAT (Takahashi et al., 1993); ΔDIC and $p_{\mathrm{A}}^{{\mathrm{CO}}_{\mathrm{2}}}$ are related to model variables (see Appendix A),

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\text{DIC}={\displaystyle \frac{m_{\mathrm{S}}}{H_{\text{mix}}{A}_{\mathrm{O}}\mathit{\varrho }{M}_{\mathrm{\mu }\mathrm{mol}}{\mathrm{10}}^{-\mathrm{15}}\mathrm{Gt}{\mathrm{g}}^{-\mathrm{1}}}},\text{(6)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{A}}^{{\mathrm{CO}}_{\mathrm{2}}}=m_{\mathrm{A}}{\mathit{\epsilon }}^{-\mathrm{1}}.\end{array}$$

The carbon cycle equation set is closed by the specification of f_decay and f_deep (Sect. 2.3), as well as ΔT, i.e., the coupling to the climate component (Sect. 2.2).

2.2 Climate component

BernSCM simulates the deviation in global mean SAT from the preindustrial state. SAT is approximated by the temperature perturbation of the surface ocean ΔT, which is calculated from heat uptake by the budget equation

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}\mathrm{\Delta }T}{\mathrm{d}t}}{c}_{\mathrm{S}}=f_{\mathrm{O}}^{\mathrm{H}}-f_{\text{deep}}^{\mathrm{H}},\end{array}$$

where c_S is the heat capacity of the surface layer, $f_{\mathrm{O}}^{\mathrm{H}}$ is ocean heat uptake, and $f_{\text{deep}}^{\mathrm{H}}$ is heat uptake by the deep ocean (and accounts for the bulk of the effective heat capacity of the ocean). Continental heat uptake is neglected due to the much higher heat conductivity of the ocean in comparison to the continent.

$f_{\mathrm{O}}^{\mathrm{H}}$ is taken to be proportional to RF (Forster et al., 2007) and the deviation of SAT from radiative equilibrium (ΔT=ΔT^eq(RF); see Table A2 for parameter definitions),

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{f}_{\mathrm{O}}^{\mathrm{H}}=\text{RF}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{\Delta }{T}^{\text{eq}}}}\right){\displaystyle \frac{A_{\mathrm{O}}}{a_{\mathrm{O}}}}.\end{array}$$

This relation follows from the assumption that feedbacks are linear in ΔT (e.g., Hansen et al., 1984). ΔT^eq is given by

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{T}^{\text{eq}}=\text{RF}{\displaystyle \frac{\mathrm{\Delta }{T}_{\mathrm{2}\times }}{{\text{RF}}_{\mathrm{2}\times }}},\end{array}$$

where ΔT_2× is climate sensitivity (defined as the equilibrium temperature change corresponding to twice the preindustrial CO₂ concentration). Equation (8) describes ocean heat uptake as the difference between RF and the climate system's response, λ⋅ΔT, with $\mathit{\lambda }=\text{RF}/\mathrm{\Delta }{T}^{\text{eq}}$ the climate sensitivity expressed in $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$.

Climate sensitivity is an external parameter, as the model does not represent the processes determining equilibrium climate response. RF of CO₂ is calculated as (Myhre et al., 1998)

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\text{RF}}_{{\mathrm{CO}}_{\mathrm{2}}}=\ln \left({\displaystyle \frac{p_{\mathrm{A}}^{{\mathrm{CO}}_{\mathrm{2}}}}{p_{\text{A0}}^{{\mathrm{CO}}_{\mathrm{2}}}}}\right){\displaystyle \frac{{\text{RF}}_{\mathrm{2}\times }}{\ln \left(\mathrm{2}\right)}},\end{array}$$

where $p_{\text{A0}}^{{\mathrm{CO}}_{\mathrm{2}}}$ is the preindustrial reference concentration of atmospheric CO₂, and RF_2× is the RF at twice the preindustrial CO₂ concentration. RF of other greenhouse gases (GHGs), aerosols, etc. can be parametrized in similar expressions involving GHG and pollutant emissions and concentrations (Prather et al., 2001). In the provided BernSCM code, non-CO₂ RF is treated as an exogenous boundary condition. Total RF is then

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\text{RF}={\text{RF}}_{{\mathrm{CO}}_{\mathrm{2}}}+{\text{RF}}_{\text{non}{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

The calculation of $f_{\text{deep}}^{\mathrm{H}}$ (Sect. 2.3) completes the climate model.

2.3 Impulse response model components

The response of a time-invariant linear system to a time-dependent forcing f can be expressed by

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }m\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}f\left(t^{\prime }\right)r(t-t^{\prime })\mathrm{d}{t}^{\prime }.\end{array}$$

The function r is the system's IRF, as can be shown by evaluating the integral for a Dirac impulse ($f\left(t^{\prime }\right)=\mathit{\delta }\left(t^{\prime }\right)$). The IRF indicates the fraction remaining in the system at time t of a pulse input at a previous time t^′. Because of linearity of the integral, any physically meaningful integrand f can be represented as a sequence of such impulses of varying size.

In BernSCM, an IRF is used to calculate the perturbation of heat and carbon in the mixed surface ocean layer (mixed layer IRF; Joos et al., 1996). For carbon,

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{m}_{\mathrm{S}}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\mathrm{O}}\left(t^{\prime }\right)r_{\mathrm{O}}(t-t^{\prime })\mathrm{d}{t}^{\prime },\end{array}$$

and similarly, for heat

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }T\left(t\right)c_{\mathrm{S}}=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\mathrm{O}}^{\mathrm{H}}\left(t^{\prime }\right)r_{\mathrm{O}}(t-t^{\prime })\mathrm{d}{t}^{\prime }.\end{array}$$

This approach has been shown to faithfully reproduce atmospheric CO₂ and SAT as simulated with the models from which the IRF is derived (Joos et al., 1996). For temperature, the linear approach works since relatively small and homogeneous perturbations of ocean temperatures do not affect the circulation strongly and can be treated as a passive tracer (Hansen et al., 2010). Note that for compatibility with commonly used units, carbon fluxes are expressed in Gt yr⁻¹, while heat fluxes are expressed in joules per second (watt) in Eqs. (13) and (14), respectively.

Equation (13) closes the ocean C budget equation (Eq. 3), as can be seen by taking the derivative with respect to time (using r(0)=1),

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{S}}}{\mathrm{d}t}}=f_{\mathrm{O}}\left(t\right)-\underset{f_{\text{deep}}}{\underbrace{\left(-\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\mathrm{O}}\left(t^{\prime }\right){\displaystyle \frac{\mathrm{d}{r}_{\mathrm{O}}}{\mathrm{d}t}}(t-t^{\prime })\mathrm{d}{t}^{\prime }\right)}},\end{array}$$

where f_deep is the flux to the deep ocean. Similarly, Eq. (14) closes the heat budget equation (Eq. 7) for the surface ocean,

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}\mathrm{\Delta }T}{\mathrm{d}t}}{c}_{\mathrm{S}}=f_{\mathrm{O}}^{\mathrm{H}}\left(t\right)-\underset{f_{\text{deep}}^{\mathrm{H}}}{\underbrace{\left(-\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\mathrm{O}}^{\mathrm{H}}\left(t^{\prime }\right){\displaystyle \frac{\mathrm{d}{r}_{\mathrm{O}}}{\mathrm{d}t}}(t-t^{\prime })\mathrm{d}{t}^{\prime }\right)}}.\end{array}$$

Another IRF is used for the carbon m_L in living or dead biomass reservoirs of the terrestrial biosphere,

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{m}_{\mathrm{L}}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\text{NPP}}\left(t^{\prime }\right)r_{\mathrm{L}}(t-t^{\prime })\mathrm{d}{t}^{\prime }.\end{array}$$

Again, Eq. (17) closes the budget equation for the land biosphere (Eq. 2), as shown by the derivative with respect to time,

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{L}}}{\mathrm{d}t}}=f_{\text{NPP}}\left(t\right)-\underset{f_{\text{decay}}}{\underbrace{\left(-\underset{-\mathrm{\infty }}{\overset{t}{\int }}{f}_{\text{NPP}}\left(t^{\prime }\right){\displaystyle \frac{\mathrm{d}{r}_{\mathrm{L}}}{\mathrm{d}t}}(t-t^{\prime })\mathrm{d}{t}^{\prime }\right)}}.\end{array}$$

The time derivative of the land IRF is also known as the decay response function Joos et al. (e.g., 1996).

The IRFs above can be expressed as a sum of exponentials,

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }r\left(t\right)=a_{\mathrm{\infty }}+\sum _k{a}_k{e}^{-t/{\mathit{\tau }}_k},\end{array}$$

where the constant term a_∞ corresponds to an infinite decay timescale.

The ocean IRF contains a positive constant coefficient a_∞, indicating a fraction of the perturbation that will remain indefinitely (implied by carbon conservation in the ocean model). CaCO₃ compensation by sediment dissolution and weathering (Archer et al., 1999) are not considered here, but could be described using analogous elimination processes with timescales on the order of 10⁴ to 10⁵ years (Joos et al., 2004). We emphasize that the implementation considering only the partitioning of excess carbon among atmosphere, land, and ocean (hence a_∞≠ 0), neglecting ocean sediment interactions and weathering flux perturbations, is only valid for timescales shorter than about 2000 years. In land biosphere models, in contrast, organic carbon is lost to the atmosphere by oxidation to CO₂ at nonzero rates, and consequently all timescales are finite (i.e., a_∞=0), and the IRF tends to zero (Fig. 2).

Figure 2 IRFs of ocean (blue) and land (green) model components (without temperature dependence). Ocean components are normalized to a common mixed-layer depth of 50 m (multiplied by H_mix∕50 m), causing initial response to deviate from 1.

Inserting the formula (Eq. 19) in the pulse response equation (Eq. 12) yields (f is a perturbation flux when a_∞≠ 0)

$$\begin{array}{}\text{(20)}& {\displaystyle }{\displaystyle }m\left(t\right)=\sum _k\underset{-\mathrm{\infty }}{\overset{t}{\int }}f\left(t^{\prime }\right)a_k{e}^{-(t-t^{\prime })/{\mathit{\tau }}_k}\mathrm{d}{t}^{\prime }+\underset{-\mathrm{\infty }}{\overset{t}{\int }}f\left(t^{\prime }\right)a_{\mathrm{\infty }}\mathrm{d}{t}^{\prime }.\end{array}$$

Thus the expression (Eq. 12) separates into a set of independent integrals m_k corresponding to the number of timescales of the response. Taking the time derivative of the expression (Eq. 20) reveals the equivalence to a system of uncoupled first-order ordinary differential equations.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{m}_k}{\mathrm{d}t}}=f\left(t\right)a_k-m_k/{\mathit{\tau }}_k;{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{\infty }}}{\mathrm{d}t}}=f\left(t\right)a_{\mathrm{\infty }}\\ \text{(21)}& {\displaystyle }& {\displaystyle }m=\sum _k{m}_k+m_{\mathrm{\infty }}\end{array}$$

The direct numerical evaluation of Eq. (12) involves integrating over all previous times at each time step. The differential form Eq. (21) allows a recursive solution, which is much more efficient, especially for long simulations (the recursive solution implemented in BernSCM is described in Appendix B).

The differential equation system (Eq. 21) can be considered to consist of several boxes, whereby each box m_k receives a fraction a_k of the input f and has a characteristic turnover time τ_k (Fig. 1). In the following this is referred to as a box model. For the mixed ocean surface layer the carbon content of box k is given by

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{m}_{{\mathrm{S}}_k}}{\mathrm{d}t}}=f_{\mathrm{O}}\left(t\right)a_{{\mathrm{O}}_k}-m_{{\mathrm{S}}_k}/{\mathit{\tau }}_{{\mathrm{O}}_k};{\displaystyle \frac{\mathrm{d}{m}_{{\mathrm{S}}_{\mathrm{\infty }}}}{\mathrm{d}t}}=f_{\mathrm{O}}\left(t\right)a_{{\mathrm{O}}_{\mathrm{\infty }}}\end{array}$$

and the change in total carbon content in the mixed layer is

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }{m}_{\mathrm{S}}=\sum _k{m}_{{\mathrm{S}}_k}+m_{{\mathrm{S}}_{\mathrm{\infty }}}.\end{array}$$

Similar equations describe the heat content in the ocean surface layer, as well as the carbon stored in the land biosphere (Fig. 1).

The timescales of an IRF describing a linear system are equivalent to the inverse eigenvalues of the model matrix of that system and may also be interpreted in the context of the Laplace transformation (Enting, 2007; Raupach, 2013). For example, the timescales of the mixing layer IRF are the inverse eigenvalues of a matrix describing a diffusive multilayer ocean model (Hooss et al., 2001). A large model matrix yields a spectrum of many eigenvalues and timescales and corresponding model boxes. In practice, IRFs are approximated with fewer fitting parameters and, equivalently, timescales (four to six in the case of BernSCM). Joos et al. (1996) used IRFs combined from two or more functions to minimize the number of parameters needed for an accurate representation. In BernSCM, simple IRFs of the form (Eq. 19) are used exclusively. This allows adequate accuracy and a consistent interpretation as a multibox model.

Thinking of IRF components as box models is conceptually meaningful. The simple Bern 4-box biosphere model (Siegenthaler and Joos, 1992), for example, contains boxes corresponding to ground vegetation, wood, detritus, and soil (Appendix A). The High-Resolution Biosphere Model (HRBM) land component (Meyer et al., 1999), however, is abstractly defined by an IRF, but corresponds to boxes which correlate with biospheric reservoirs. However, since different box models may show a similar response, in practice the coefficients a_k and timescales τ_k may not be uniquely defined by the IRF and should be interpreted primarily as abstract fitting parameters (Enting, 2007; Li et al., 2009).

The IRF representation is, strictly speaking, only valid if the described subsystem is linear and the timescales of the system are time invariant. Then, the response function r does not depend on time and on state variables. In the BernSCM, major nonlinearities in the carbon cycle, namely air–sea gas exchange and the nonlinear carbonate chemistry, and changes in NPP in response to changes in environmental conditions are treated by separate nonlinear equations (Eqs. 4 and 5), while surface-to-deep ocean transport of carbon and heat and respiration of carbon in litter and soils are viewed as approximately linear processes using IRFs. However ocean circulation and the respiration of carbon from soil and litter is likely to change under global warming, violating the assumption of linearity. In practice, the IRF representation remains a useful approximation as long as the impact of associated nonlinearities on simulated atmospheric CO₂ and temperature remain moderate.

The interpretation of the IRF representation as a box model provides a starting point for considering nonlinearities in the response. To account for nonlinearities, the response timescales τ_k and the coefficients a_k may be gradually adjusted as a function of state variables such as temperature. As the integral form (Eq. 12) involves integration over the whole history at each time step, changing parameters along the way would result in inconsistencies. In contrast, the differential or box model form (Eq. 21) does not depend on previous time steps. Changing the model parameters from one step to the next thus equates to applying a slightly different model at each time step. Within each time step, the parameters remain constant, and the solution for the linear case applies. As time steps are small compared to the whole simulation, this discretization yields accurate results, which is confirmed by the close agreement between the different time resolutions (Table B1).

Varying coefficients have been successfully implemented and tested for the HRBM land component and its decay IRF (Meyer et al., 1999). In this way, the enhancement of biomass decay by global warming is captured (see the Appendix A and Sect. 3.1). In such a modification, the advantage of the IRF and the equivalent box model representation – the faithful representation of the characteristic response timescale of a model system – is largely maintained, while at the same time the impact of time- and state-dependent system responses on simulated outcomes is approximated.

3 Illustrative simulations with the BernSCM

3.1 Model setup for sensitivity analyses and uncertainty assessment

The carbon cycle–climate uncertainty of simulations with BernSCM can be assessed in two ways. First, to assess structural uncertainty, different substitute models for the ocean and land components can be used. Currently, this approach is quite limited by the set of available substitute models (see Appendix A). Second, parameter uncertainty can be assessed by varying the temperature and CO₂ sensitivities of the model, based on a standard set of components that represent the key dependencies as completely as possible (here, the IRF substitutes for the high-latitude exchange/interior diffusion–advection (HILDA) ocean model (Joos et al., 1996) and for the HRBM land biosphere model (Meyer et al., 1999) are used in the standard setup).

The uncertainties of the global carbon cycle concern the sensitivity of the modeled fluxes of carbon and heat to changing atmospheric CO₂ and climate. Key uncertainties strongly affecting the overall climate response are associated with land C storage: the dependency of NPP on CO₂ (CO₂ fertilization), and the dependency of land C on temperature (f_decay increases with warming). This gives rise to large and opposed carbon flux perturbations which are both very uncertain in magnitude (Le Quéré et al., 2016). While all substitute land models available for BernSCM include CO₂ fertilization, only the HRBM substitute model represents temperature sensitivity of biomass decay (Appendix A2).

As for the ocean, the uncertainty of heat uptake into the surface ocean is treated in terms of climate sensitivity (Eq. 8). The efficiency of the uptake of heat ($f_{\text{deep}}^{\mathrm{H}}$) and carbon (f_deep) into the deep ocean is not sensitive to temperature, as the currently available substitute models all represent a fixed circulation pattern (IRF or box model parameters are not temperature dependent; Appendix A1). The nonlinear chemistry of CO₂ dissolution in the surface ocean (Eq. 4), which determines the sensitivity of ocean C uptake to atmospheric CO₂, is scientifically well established (Dickson et al., 2007; Orr and Epitalon, 2015) and is not treated as an uncertainty in BernSCM. The temperature sensitivities of NPP and CO₂ dissolution in the surface ocean are treated as uncertain here, but have secondary influence on the climate response.

Figure 3 Fraction of realized warming (temperature divided by the equilibrium temperature for the current RF) for idealized experiments with prescribed atmospheric CO₂ concentration increase from preindustrial levels; panel (a) shows an exponential CO₂ increase by 1 % yr⁻¹ over 140 years to approximately 4 times the preindustrial concentration (and linear increase in RF); panel (b) shows an abrupt increase to 4-fold CO₂ concentration. BernSCM simulations are shown for climate sensitivities of 2, 3, and 4.5 K and the three available ocean model substitutes as indicated in the legend. Arrows in panel (a) indicate the corresponding warming fractions at year 99 compiled by Frölicher and Paynter (2015, Tables 1, 2 in their Supplement) for Earth system models (ESMs, right-pointing) and Earth system Models of Intermediate Complexity (EMICs, left-pointing); arrow colors indicate climate sensitivities below 2.5 K (green), between 2.5 and 3.5 K (black), and above 3.5 K (red).

Similar to previous studies using models from the Bern family (Plattner et al., 2008; Joos et al., 2001; Meehl et al., 2007; Van Vuuren et al., 2008), the parameter uncertainty range is assessed using the following setups.

Coupled.

All temperature and CO₂ sensitivities are set to their standard values.

Uncoupled.

All sensitivities are set to zero (except for the ocean CO₂ dissolution chemistry).

C-only.

Only CO₂ dependencies are considered (CO₂ fertilization).

T-only.

Only temperature dependencies are considered in the land module (NPP, decay).

We performed simulations with these different setups. In Sect. 4.2, we probe the timescales of the temperature response in simulations in which atmospheric CO₂ is abruptly (instantaneously) quadrupled or by increasing CO₂ radiative forcing linearly within 140 years. In Sect. 4.3, we probe the response of the coupled system to a pulse-like release of 100 Gt C into the atmosphere. Finally in Sect. 4.4, we analyze carbon cycle–climate feedbacks relying on simulations over the industrial period and for the SRES A2 scenario. BernSCM results are compared with the results from three multimodel intercomparison projects: the Climate Model Intercomparison Project 5 (CMIP5) with results as summarized by Frölicher and Paynter (2015), an analysis of carbon dioxide and climate impulse response functions (Joos et al., 2013, here referred to as IRFMIP), and the C4MIP climate–carbon cycle feedback analysis (Friedlingstein et al., 2006).

3.2 Fraction of realized warming and idealized forcing experiments

The climate response of BernSCM is illustrated using idealized simulations with prescribed forcing. One series of simulations (a) was run for CO₂ concentration increasing exponentially from the preindustrial value by 1 % yr⁻¹ over 140 years to approximately 4 times the preindustrial concentration, corresponding to a linear increase in RF (Fig. 3a); in a second series of simulations (b), CO₂ was abruptly increased to 4 times the preindustrial concentration (Fig. 3b).

Frölicher and Paynter (2015) compare similar simulations of Earth system models (ESMs) performed within the Coupled Model Intercomparison Project Phase 5 (CMIP5), and Earth System Models of Intermediate Complexity (EMICs) (Joos et al., 2013). As a model comparison metric sensitive to the long-term climate response, Frölicher and Paynter (2015) use the fraction of realized warming, defined by the ratio of the temperature response at a given year and the equilibrium temperature for the corresponding RF. They show that the smaller realized warming of ESMs in comparison to EMICs (Fig. 3) is connected to a higher long-term warming response; this implies an increase in the coefficient relating global warming to cumulative carbon emissions on multicentennial timescales and suggests a lower quota on allowed emissions for a given global warming target (Frölicher and Paynter, 2015). The realized warming fraction simulated with BernSCM is in good agreement with the responses of the ESMs (and lower on average than that of the EMICs). The validity of the IRF approach has also been shown by Good et al. (2011) using a SCM to reconstruct and interpret atmosphere–ocean general circulation model (AOGCM) projections. For the 150-year timescale of the CMIP5 experiments, Geoffroy et al. (2012a, b) show that the climate response of AOGCMs is well captured by a two-layer energy balance model with two effective response timescales.

Figure 4 IRFMIP pulse response range compared to BernSCM range for parameter uncertainty (colors according to legend) and structural uncertainty, with model versions HILDA–HRBM (solid lines), HILDA–4-box (dots), and Princeton–HRBM (dashed). Standard climate sensitivity is 3 ^∘C, and a climate sensitivity range of 2–4.5 ^∘C is shown by the white area (envelope of all BernSCM runs). Single-model ensemble ranges from IRFMIP are included as error bars indicating the 5–95 % range and dots indicating the median. The multimodel IRFMIP range is shown by box plots indicating median (bold black line), first quartiles (box), and extreme values (whiskers) excluding outliers deviating from the median by more than 1.2 times the interquartile distance (asterisks).

In BernSCM, the fraction of realized warming depends primarily on the choice of climate sensitivity and is qualitatively similar for the different model setups. Such a clear relationship is not seen in the EMS and EMICs. Thus the structural uncertainty and model differences of complex models are not fully represented in BernSCM. The BernSCM climate response to abrupt warming (Fig. 3b) is qualitatively similar, especially on multicentennial timescales.

3.3 Impulse response experiment

Coupled carbon cycle–climate models can be characterized and compared based on their response to a CO₂ emission pulse to the atmosphere (Joos et al., 2013). The airborne fraction (AF) denotes the fraction of emissions found in the atmosphere at a given time. In IRFMIP, the AF for a pulse of 100 Gt C, emitted on top of current (i.e., year 2010) atmospheric CO₂ concentrations, was simulated by a set of 15 carbon cycle–climate models of different complexity. For three of these models (Bern3D-LPJ, GENIE, MAGICC), ensembles sampling the parameter uncertainty of these models are included in IRFMIP. Thus, IRFMIP captures structural as well as parameter uncertainty.

The IRFMIP pulse experiment was repeated with BernSCM, exploring parameter uncertainty of the carbon cycle (Sect. 3.1), as well as structural uncertainty, using the ocean model IRFs HILDA and Princeton (Sarmiento et al., 1992) in various combinations with the land biosphere components HRBM and the Bern 4-box model (Fig. 4). Simulations were run for equilibrium climate sensitivities of 3 (standard setup), 2, and 4.5 ^∘C.

The AF simulated with BernSCM broadly agrees with the set of simulations from IRFMIP. At 100 years after the pulse, the AF is 0.40 (0.34–0.57) for a climate sensitivity of 3 ^∘C (for coupled setup with uncertainty range in brackets). Climate sensitivity uncertainty only slightly affects the upper end of this range (Fig. 4). For AF simulated with BernSCM, the standard coupled setup is close to the IRFMIP multimodel median. The BernSCM uncertainty range is asymmetric, like the IRFMIP multimodel range. For the MAGICC and GENIE ensembles, the medians also correspond with the BernSCM standard case, while the uncertainty ranges are more symmetric.

The BernSCM SAT response also broadly agrees with IRFMIP. The standard coupled simulation is somewhat lower than the IRFMIP median, which is explained in part by the climate sensitivity (3 ^∘C) being slightly lower than the IRFMIP average (3.2 ^∘C). The short-term temperature response of BernSCM in particular is on the lower side of the IRFMIP range, suggesting stronger ocean mixing. The quickest initial temperature increase in the BernSCM simulations is obtained with the Princeton ocean model component (dashed lines), which shows a slower initial mixing to the deep ocean than the other implemented components (Fig. 2). The comparability of the SAT projections is limited, as the range of climate sensitivities considered in the BernSCM simulations (2–4.5 ^∘C) differ somewhat from those of the IRFMIP multimodel set (1.5–4.6 ^∘C) and the single model ensembles (1.9–5.7 ^∘C) and are compounded with RF differences resulting from the uncertainty in atmospheric CO₂.

Figure 5 Land, ocean, and airborne fractions of the 100 Gt C CO₂ pulse shown in Fig. 4 for the coupled (solid lines and colored areas), T-only (dashed), C-only (dotted), and uncoupled (dash-dotted) model setups. In the T-only case, the land biosphere exhibits a net release (light green shading), and the ocean uptake consists of the sum of this area and the area delimited by the dashed line below the line at 1; for the uncoupled case, land uptake is zero and ocean uptake extends from the dashed–dotted line to unity.

Figure 5 shows how the added carbon is redistributed within the Earth system. In the coupled setup, the fraction of the initial pulse sequestered by the land and by the ocean increases over the first century, while the airborne fraction decreases. After 100 years, slightly more than 20 % of the added carbon is stored in the land and about 40 % in the ocean. The ocean continues to sequester excess carbon in the following centuries to become the dominant sink for excess carbon. In contrast, the land returns part of the sequestered carbon back to the atmosphere and ocean as decreasing atmospheric CO₂ reduces the modeled CO₂ fertilization of the land biosphere. In the T-only setup, in which CO₂ fertilization is not operating, the land is a source of carbon to the atmosphere due to accelerated soil turnover in response to warming. The largest land sink is simulated in the C-only setup, in which soil turnover timescales remain invariant and CO₂ fertilization is on. The different BernSCM setups span a range of plausible land biosphere and ocean responses to continued anthropogenic CO₂ emissions as reflected in the simulated range in the airborne fraction (Figs. 4a and 5).

3.4 Carbon cycle–climate feedbacks

Climate models with explicit and detailed carbon cycle components exhibit a wide range of responses, as shown in the intercomparison studies of climate models with a detailed carbon cycle, C4MIP (Friedlingstein et al., 2006) and CMIP5 (Jones et al., 2013). The authors analyzed the feedback of carbon cycle–climate models using linearized sensitivity measures. These are derived from a simulation with temperature dependence (“coupled”) and one without (“uncoupled”; note that these names have a different meaning in BernSCM). Total CO₂ emissions for the coupled (left-hand side) and uncoupled (right-hand side) simulations can be expressed as

$$\begin{array}{}\text{(24)}& \mathrm{\Delta }{C}_{\mathrm{A}}^{\mathrm{c}}(\mathit{\epsilon }+{\mathit{\beta }}_{\mathrm{L}}+{\mathit{\beta }}_{\mathrm{O}}+\mathit{\alpha }({\mathit{\gamma }}_{\mathrm{L}}+{\mathit{\gamma }}_{\mathrm{O}}\left)\right)=\mathrm{\Delta }{C}_{\mathrm{A}}^{\mathrm{u}}(\mathit{\epsilon }+{\mathit{\beta }}_{\mathrm{L}}+{\mathit{\beta }}_{\mathrm{O}}),\end{array}$$

where ΔC_A is the cumulative change in atmospheric CO₂ (in parts per million) in the coupled (c) or uncoupled (u) cases, and the terms in parentheses represent the total sensitivity of C storage to ΔC_A; in particular, β is the sensitivity of carbon storage to atmospheric CO₂ (in Gt C ppm⁻¹) on land (β_L) or in the ocean (β_O). γ is similarly the sensitivity in carbon storage to climate change, and α is the linear transient climate sensitivity to CO₂ (^∘C ppm⁻¹) as in Friedlingstein et al. (2006); ε converts ppm to Gt C (see Table A2; the formula in the original paper implies identical units for atmospheric and stored carbon).

The climate–carbon cycle feedback is measured by the feedback metric g, defined by

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle \frac{\mathrm{\Delta }{C}_{\mathrm{A}}^{\mathrm{c}}}{\mathrm{\Delta }{C}_{\mathrm{A}}^{\mathrm{u}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-g}}\end{array}$$

and is thus estimated by

$$\begin{array}{}\text{(26)}& {\displaystyle }{\displaystyle }g=-{\displaystyle \frac{\mathit{\alpha }({\mathit{\gamma }}_{\mathrm{L}}+{\mathit{\gamma }}_{\mathrm{O}})}{\mathit{\epsilon }+{\mathit{\beta }}_{\mathrm{O}}+{\mathit{\beta }}_{\mathrm{L}}}}.\end{array}$$

Thus the feedback strength scales with the assumed climate sensitivity and the temperature sensitivities and is reduced by CO₂-induced sinks.

The C4MIP study used a SRES A2 emission scenario to compare the carbon cycle sensitivities of a range of models. As in the C4MIP exercise, BernSCM was run for SRES A2 without any non-CO₂ forcings (Fig. 6; prescribed historical and scenario emissions were smoothed with the R smooth.spline function (R Core Team, 2015) for 41 degrees of freedom for use with different time steps). Land use was treated as an exogenous CO₂ emission, while the land model simulates an undisturbed biosphere.

Figure 6 BernSCM simulations of the SRES A2 scenario used for C4MIP, with a climate sensitivity of 2.5 ^∘C and the HILDA–HRBM ocean–land components. Results for three numerical schemes are overlaid; (i) 0.1-year Euler forward time step (solid thin line), (ii) 1-year implicit time step (dashed bold line), (iii) 10-year implicit time step with piecewise linear approximation of fluxes (circles); the difference at this resolution is only visible in the C uptake. The C4MIP model range at 2100 is indicated by grey bars; numbers above or below the bars indicate values outside of the chart range.

The BernSCM sensitivity setups can be expressed in terms of the C4MIP sensitivity metrics: T-only corresponds to β_L=0, C-only to ${\mathit{\gamma }}_{\mathrm{L}}={\mathit{\gamma }}_{\mathrm{O}}=\mathrm{0}$, and uncoupled to ${\mathit{\beta }}_{\mathrm{L}}={\mathit{\gamma }}_{\mathrm{L}}={\mathit{\gamma }}_{\mathrm{O}}=\mathrm{0}$. This can be used to estimate climate–carbon cycle feedback g captured in BernSCM. The sensitivity metrics for the BernSCM standard simulation (HILDA-HRBM with coupled carbon cycle) lie within the C4MIP range (Table 1). The uncertainty range for BernSCM, however, is not congruent with the multimodel range of C4MIP. Maximum and standard sensitivity for BernSCM are practically identical. Notably, this sensitivity is smaller (absolutely) than the C4MIP average for the land carbon response to CO₂ increase and warming. The resulting gain g is also smaller, though this results in large part from the lower climate sensitivity in BernSCM (which corresponds to 2.5 ^∘C as used for the Bern-CC model contribution to C4MIP). The lower end (in absolute terms) of the BernSCM carbon cycle sensitivity range is, however, zero per definition for all but the ocean-CO₂ sensitivity β_O (see Sect. 3.1). As a consequence, the climate–carbon cycle feedback range also includes zero. In contrast, the C4MIP range does not include zero for all sensitivity parameters.

The land carbon uptake until 2100, under the different BernSCM configurations, varies over 500 Gt C (Fig. 6), more than 3 times the range of ocean uptake (180 Gt C). This partly reflects the limited coverage of the uncertainty in ocean mixing but also the fact that the land carbon sink is, together with the source related to land use, the most uncertain item in the budget (Le Quéré et al., 2009).

4 Discussion

We simulated illustrative scenarios from two recent multimodel studies, C4MIP and IRFMIP, to compare BernSCM to the literature of carbon cycle–climate models. The results show that BernSCM is broadly representative of the current understanding of the global carbon cycle–climate response to anthropogenic forcing (in a time-averaged sense that does not address internal variability). The BernSCM uncertainty range in CO₂ and SAT projections is broadly similar to the ranges spanned by probabilistic single-model ensembles and multimodel “ensembles of opportunity” such as the 15 IRFMIP models. The BernSCM uncertainty range shown consists mainly of parameter uncertainty and to a small extent of structural uncertainty. For the standard coupled model setup, the sensitivities of ocean and land carbon uptake to changing CO₂ and climate (Table 1) of BernSCM are within the range of the detailed carbon cycle models in C4MIP. However, as some C4MIP models show much higher sensitivities, the BernSCM range does not capture the full C4MIP multimodel range. However, the C4MIP set is unlikely to sample uncertainty exhaustively, as each model contributed only a single, most likely simulation. Thus it does not include zero (or weak) sensitivities, whereas the BernSCM range does.

As Fig. 6 shows, solutions with different time steps and numerical schemes as implemented in BernSCM are largely equivalent for a sufficiently smooth forcing. This offers the flexibility to opt for simplicity of implementation or maximum speed as required by the application (see also Appendix B).

BernSCM does not explicitly distinguish between surface atmosphere and surface ocean temperature to compute global mean SAT perturbation. This is in contrast to some energy balance calculations used to analyze results from state-of-the-art ESMs (e.g., Geoffroy et al., 2012b). The BernSCM approach follows earlier work of Siegenthaler and Oeschger (1984). It is further guided by the similarity in reconstructions of marine nighttime air and sea surface temperature perturbations (T. F. Stocker et al., 2013) that are consistent with the short, monthly relaxation timescale for air–sea heat exchange. The focus of the BernSCM is on the representation of the transport of heat from the surface into the thermocline and the deep ocean on decadal to multicentury timescales, while information on seasonal and spatial changes such as on land–sea air temperature differences or polar amplification may be obtained by applying suitable spatial perturbation patterns as derived from state-of-the-art models.

Table 1 C4MIP sensitivity metrics. The BernSCM range covers the carbon cycle settings as discussed in Sect. 3.1, and different combinations of model components (HILDA–HRBM, HILDA–4-box, Princeton–HRBM); the C4MIP range covers all participating models.

Currently, a limited set of substitute models is available and included with BernSCM. The simple structure and open-source policy of BernSCM allows users to address these current limitations according to the needs of their applications. More components can be added using the existing ones as a template. This requires the specification of the IRF and the parametrization of gas exchange for the surface ocean or NPP for the land biosphere (as described in Joos et al., 1996; Meyer et al., 1999).

Ocean transport is known to vary under climate change with some consequences for heat and carbon uptake (Joos et al., 1999b). Here, we applied time-invariant ocean transport parameters ($a_{{\mathrm{O}}_k}$, ${\mathit{\tau }}_{{\mathrm{O}}_k}$). It is in principle possible to represent temperature dependency of ocean transport in a similar way as it is performed for the climate dependency of heterotrophic respiration for the HRBM land biosphere substitute model (Meyer et al., 1999). In the current BernSCM version, the same IRF parameters are applied for the transport of carbon and heat from the surface ocean to the interior ocean. Thereby, it is implicitly assumed that the spatial pattern of change is the same for temperature and carbon. This appears to be a reasonable first-order approximation on decadal to century timescales as perturbations in temperature and carbon show similar patterns with decreasing perturbations from the surface to depth. In future efforts, one may differentiate the ocean IRF for heat and carbon, in particular when more information from long-term multicentury to millennial-scale ESM simulations becomes available. The application of the same IRF for carbon and heat in individual model runs implies that modeled carbon and heat transport tend to be physically consistent. In contrast, some other simple models employ different transport parameters for heat and carbon and varied these parameters independently in probabilistic studies.

A distribution of timescales applies to ocean transport processes as evidenced by observations of transient and time-dependent tracers such as chlorofluorocarbons and bomb-produced and natural radiocarbon and biogeochemical tracers (Key et al., 2004; Olsen et al., 2016). This continuum is sometimes approximated by one timescale, also termed heat uptake efficiency (e.g., Gregory et al., 2009), and by two timescales, as in Geoffroy et al. (2012b). The one-to-two timescale approximations were used to analyze relatively short ESM simulations that do not yet reveal the multicentury response timescales of the deep ocean. We note that the equivalent ocean depth of the simple energy balance model of Geoffroy et al. (2012b) for their AOGCM ensemble is only 1182 m compared to a mean ocean depth of about 3800 m. The ocean IRFs used in the BernSCM are derived from long simulations with ocean-only or simplified models. The range of distinct timescales used to construct the IRF faithfully approximates the sub-annual to multicentury response continuum of the parent models as shown in earlier work (Joos et al., 1996). Further, the BernSCM IRF model represents the heat capacity of the entire ocean.

The BernSCM model may be extended to model perturbation in the signatures and exchange fluxes of the carbon isotopes ¹³C and ¹⁴C as demonstrated in earlier work (Joos et al., 1996). This was not implemented here to keep the code as simple as possible and as most potential users are likely concerned with the evolution of climate and atmospheric CO₂ .

A potential future application of BernSCM is to use it as an emulator of the global long-term response of complex climate–carbon cycle models by adding the corresponding substitute model components. Additionally, pattern scaling can be applied to transfer the global mean temperature signal into spatially resolved changes in surface temperature, precipitation, cloud cover, etc., exploiting the correlation of global SAT with regional and local changes (Hooss et al., 2001). This allows us to drive spatially explicit models, e.g., of terrestrial vegetation Joos et al. (as in 2001); Strassmann et al. (as in 2008) or impacts related to climate change (e.g., as in Hijioka et al., 2009). Patterns of change are generally similar across models for temperature, whereas patterns in precipitation are more uncertain and show greater variability among models (Knutti and Sedlacek, 2013) and are forcing dependent (Shine et al., 2015). We also note that natural variability strongly influences the space-time evolution of climate change (Deser et al., 2012). Patterns may be scaled with changes in global mean SAT as indicated in Fig. 1 or dependencies on radiative forcing may be considered (Shine et al., 2015)

The addition of further alternative model components will extend the structural uncertainty that can be represented with BernSCM. A sufficient coverage of structural uncertainty could allow the interpolation among alternative model components to represent uncertainty with scalable parameters (and removing the distinction between structural and parameter uncertainty). Such a parametrization of the uncertainty would enhance the possibilities for probabilistic applications of BernSCM, although more sophisticated models are available for observation-constrained probabilistic quantification of climate targets (Holden et al., 2010; Steinacher and Joos, 2016; Steinacher et al., 2013).

5 Conclusions

BernSCM is a reduced-form carbon cycle–climate model that captures the characteristics of the natural carbon cycle and the climate system essential for simulating the global long-term response to anthropogenic forcing. Simulated atmospheric CO₂ concentrations and SATs are in good agreement with results from two comprehensive multimodel ensembles. Process detail is minimal, due to the use of IRFs for system compartments that can be described linearly and nonlinear parametrizations governing the carbon fluxes into these compartments. This framework allows, in particular, the representation of the wide range of response timescales of the ocean and land biosphere and the nonlinear chemistry of CO₂ uptake in the surface ocean – both essential for reliably simulating the global climate response to arbitrary forcing scenarios.

Due to its structural simplicity and computational efficiency, BernSCM has many potential applications. In combination with pattern scaling, BernSCM can be used to project spatial fields of impact-relevant variables for applications such as climate change impact assessment, coupling with spatially explicit land biosphere models, etc. With alternative numerical solutions of varying complexity and stability to choose from, applications range from educational to computationally intensive integrated assessment modeling. BernSCM also offers a model-based alternative to global warming potentials for estimation of the climate impact of emissions and can be used to quantify climate benefits of mitigation options by applying emissions- or concentration-driven simulations.

The generic implementation of linear IRF components offers a transparent, extensible climate model framework. Current limitations concern the number of available substitute models (limiting the uncertainty range represented), and ocean transport not influenced by climate change. An addition of further alternative model components and more flexible representation of sensitivities in terms of continuously variable parameters would further increase the models' usefulness, for example for probabilistic applications.

Appendix A: Model parameters and parametrizations

A1 Ocean

Currently available ocean components include substitute models for the high-latitude exchange/interior diffusion–advection model (HILDA Joos et al., 1996), Bern2D (Stocker et al., 1992), and the Princeton general circulation model (GCM) (Sarmiento et al., 1992). Ocean model parameters of the equations described in the main text are listed in Table A3 for the mixed-layer IRF/box models and in Table A2 for other equations. The IRF/box model parameters given here are recalculated by fitting a sum of six exponential functions and one constant to the original response functions as given in Joos and Bruno (1996). The original functions treated the first few years separately; the approximation to a purely exponential form simplifies the equations and has a negligible effect on accuracy. The parametrization of ocean surface CO₂ pressure is the same for all available ocean components and is given below.

Table A1 Model variables.

Table A2 Model parameters.

^* The first value is used in the climate component equations, the value in parentheses in the C cycle component equations.

Table A3 Mixed-layer IRF/Box parameters.

Table A4 Land C stock IRF/Box parameters.

Ocean surface CO₂ pressure perturbations are fitted as a function of the globally averaged unperturbed surface temperature T^* and perturbations in dissolved inorganic carbon (DIC) by Joos et al. (1996) using carbonate chemistry coefficients summarized by Millero (1995):

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }{\left.p_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}\right|}_{T^{*}}=& {\displaystyle }(\mathrm{1.5568}-\mathrm{1.3993}\times {\mathrm{10}}^{-\mathrm{2}}{T}^{*})\mathrm{\Delta }\text{DIC}\\ {\displaystyle }& {\displaystyle }+(\mathrm{7.4706}-\mathrm{0.20207}{T}^{*}){\mathrm{10}}^{-\mathrm{3}}\mathrm{\Delta }{\text{DIC}}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-(\mathrm{1.2748}-\mathrm{0.12015}{T}^{*}){\mathrm{10}}^{-\mathrm{5}}\mathrm{\Delta }{\text{DIC}}^{\mathrm{3}}\\ {\displaystyle }& {\displaystyle }+(\mathrm{2.4491}-\mathrm{0.12639}{T}^{*}){\mathrm{10}}^{-\mathrm{7}}\mathrm{\Delta }{\text{DIC}}^{\mathrm{4}}\\ {\displaystyle }& {\displaystyle }-(\mathrm{1.5468}-\mathrm{0.15326}{T}^{*}){\mathrm{10}}^{-\mathrm{10}}\mathrm{\Delta }{\text{DIC}}^{\mathrm{5}}.\end{array}$$

The expression holds for unperturbed global average surface water temperature T^* between 17.7 and 18.3 ^∘C and for $\mathrm{\Delta }{p}_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}$ between 0 and 1320 ppm.

Ocean surface CO₂ pressure for global surface temperature perturbation ΔT (Takahashi et al., 1993):

$${\displaystyle }{\displaystyle }{p}_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}={\left.p_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}\right|}_{T^{*}}\cdot e^{\mathrm{0.0423}\mathrm{\Delta }T}.$$

A2 Land biosphere

Currently available land biosphere components include substitute models for the High-Resolution Biosphere Model (HRBM) (Meyer et al., 1999) and the 4Box biosphere model (Siegenthaler and Joos, 1992).

For the HRBM model, temperature-dependent IRF/box model parameters as given by Meyer et al. (1999) are implemented:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{a}}_k={\displaystyle \frac{a_k{e}^{s_{a_k}T}}{{\sum }_j{a}_j{e}^{s_{a_j}T}}},\\ {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_k={\mathit{\tau }}_k{e}^{-s_{{\mathit{\tau }}_k}T},\end{array}$$

where ${\stackrel{\mathrm{̃}}{a}}_k$, ${\stackrel{\mathrm{̃}}{\mathit{\tau }}}_k$ are the adjusted and a_k, τ_k the unperturbed parameters. The IRF/box model parameter values for HRBM and the 4Box model are listed in Table A4. The temperature sensitivities of the HRBM IRF are parametrized for a warming of up to 5 ^∘C.

Net primary production for HRBM is given by (Meyer et al., 1999)

$$\begin{array}{ll}{\displaystyle }{\left.\text{NPP}\left(p\right)\right|}_{\mathrm{\Delta }T=\mathrm{0}}=& {\displaystyle }-e^{\mathrm{3.672801}}+e^{-\mathrm{0.430818}}\cdot p\\ {\displaystyle }& {\displaystyle }-e^{-\mathrm{6.145559}}\cdot p^{\mathrm{2}}+e^{-\mathrm{12.353878}}\cdot p^{\mathrm{3}}\\ {\displaystyle }& {\displaystyle }-e^{-\mathrm{19.010800}}\cdot p^{\mathrm{4}}+e^{-\mathrm{26.183752}}\cdot p^{\mathrm{5}}\\ {\displaystyle }& {\displaystyle }-e^{-\mathrm{34.317488}}\cdot p^{\mathrm{6}}-e^{-\mathrm{41.553715}}\cdot p^{\mathrm{7}}\\ {\displaystyle }& {\displaystyle }+e^{-\mathrm{48.265138}}\cdot p^{\mathrm{8}}-e^{-\mathrm{56.056095}}\cdot p^{\mathrm{9}}\\ {\displaystyle }& {\displaystyle }+e^{-\mathrm{64.818185}}\cdot p^{\mathrm{10}},\end{array}$$

where p is atmospheric CO₂ pressure. This expression holds up to a CO₂ concentration of 1274 ppm and is capped at that value. The model includes growth enhancement by SAT increase (but without a dynamical vegetation):

$$\begin{array}{ll}{\displaystyle }\text{NPP}(p,\mathrm{\Delta }T)=& {\displaystyle }{\text{NPP}}_{\mathrm{0}}\\ {\displaystyle }& {\displaystyle }\cdot (\mathrm{1}+\mathrm{0.11780208}\tanh (\mathrm{\Delta }T/\mathrm{50.9312421})\\ {\displaystyle }& {\displaystyle }+\mathrm{0.002430513}\cdot \tanh (\mathrm{\Delta }T/\mathrm{8.85326739})).\end{array}$$

This expression holds up to a SAT increase of 5 ^∘C.

Net primary production for the 4Box model is described after (Enting et al., 1994; Schimel et al., 1996):

$${\displaystyle }{\displaystyle }\text{NPP}={\text{NPP}}_{\mathrm{0}}+{\text{NPP}}_{\mathrm{0}}\cdot \mathit{\beta }\cdot \log \left(p^{{\mathrm{CO}}_{\mathrm{2}}}/p_{\mathrm{0}}^{{\mathrm{CO}}_{\mathrm{2}}}\right),$$

where NPP₀ is undisturbed NPP.

Appendix B: Implementation of the pulse-response model

B1 Discretization

For the solution of the pulse-response equation (Eq. 12), two discrete approximations are implemented, using the separation by timescales in Eq. (20) or, equivalently, in the differential equation system (Eq. 21). The recursive solution for a time step Δt can be obtained from Eq. (20) by substituting $t=t_n=t_{n-\mathrm{1}}+\mathrm{\Delta }t$, and $s=t^{\prime }-t_{n-\mathrm{1}}$,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m}_n=m_{\mathrm{\infty }n}+\sum _k{m}_{kn}\\ {\displaystyle }& {\displaystyle }{m}_{kn}=m_{kn-\mathrm{1}}{e}^{-\mathrm{\Delta }t/{\mathit{\tau }}_k}+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}f(t_{n-\mathrm{1}}+s)a_k{e}^{-(\mathrm{\Delta }t-s)/{\mathit{\tau }}_k}\mathrm{d}s\\ \text{(B1)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{\infty }n}=m_{\mathrm{\infty }n-\mathrm{1}}+\underset{\mathrm{0}}{\overset{\mathrm{\Delta }t}{\int }}f(t_{n-\mathrm{1}}+s)a_{\mathrm{\infty }}\mathrm{d}s,\end{array}$$

where $m_n=m\left(t_n\right)=m(t_{n-\mathrm{1}}+\mathrm{\Delta }t)$.

First, f can be taken as constant over a sufficiently short time step $\mathrm{\Delta }t=t_i-t_{i-\mathrm{1}}$. Evaluating equations (Eq. B1) yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m}_{kn}=m_{kn-\mathrm{1}}{e}^{-\mathrm{\Delta }t/{\mathit{\tau }}_k}+f\left(t^{*}\right)a_k{\mathit{\tau }}_k(\mathrm{1}-e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k})\\ \text{(B2)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{\infty }n}=m_{\mathrm{\infty }n-\mathrm{1}}+f\left(t^{*}\right)a_{\mathrm{\infty }}\mathrm{\Delta }t,\end{array}$$

where t^* is chosen to be t_n−1 (for explicit forward solution) or t_n (for implicit backward solution).

Second, for longer time steps, a better approximation is obtained by assuming linear variation in f over each time step. This yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m}_{kn}=m_{kn-\mathrm{1}}{e}^{-\mathrm{\Delta }t/{\mathit{\tau }}_k}\\ {\displaystyle }& {\displaystyle }+f_{n-\mathrm{1}}{a}_k{\mathit{\tau }}_k\left({\displaystyle \frac{{\mathit{\tau }}_k}{\mathrm{\Delta }t}}(\mathrm{1}-e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k})-e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k}\right)\\ {\displaystyle }& {\displaystyle }+f_n{a}_k{\mathit{\tau }}_k\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\tau }}_k}{\mathrm{\Delta }t}}(\mathrm{1}-e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k})\right)\\ \text{(B3)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{\infty }n}=m_{\mathrm{\infty }n-\mathrm{1}}+{\displaystyle \frac{f_{n-\mathrm{1}}+f_n}{\mathrm{2}}}{a}_{\mathrm{\infty }}\mathrm{\Delta }t.\end{array}$$

B2 Numerical schemes

For the solution of the BernSCM model equations, both explicit and implicit time stepping is implemented.

The stability requirement for the numerical solution depends on the equilibration time for the ocean surface CO₂ pressure $p_{\mathrm{S}}^{{\mathrm{CO}}_{\mathrm{2}}}$. Due to the buffering of the carbonate chemistry, the CO₂ equilibration time is smaller than the gas diffusion timescale (∼10 years) by a ratio given by the buffer factor. For undisturbed conditions (buffer factor ≃10) the equilibration time is about 1 year. With increasing DIC, the buffer factor increases and the equilibration time shortens, making the equation system stiffer. Accordingly, when the model is solved explicitly with a time step of 1 year, instability typically occurs after sustained carbon uptake by the ocean, which can occur in many realistic scenarios.

For the tested scenario range, the explicit solution is stable at a time step on the order of 0.1 year, for which the piecewise constant approximation is accurate. For a larger step size, an implicit solution is required to guarantee stability.

The piecewise constant approximation is adequate for time steps up to 1 year, and the piecewise linear approximation is adequate for up to decadal time steps. An overview of the performance of three representative settings (set at compile time) for the C4MIP A2 scenario is given in Table B1.

Table B1 Performance and accuracy for time steps of 1–10 years relative to a reference with a time step of 0.1 year. The reference simulation is solved explicitly; otherwise an implicit solution was used. The average execution time of the time integration loop is given as a fraction of the explicit case. For atmospheric CO₂ and SAT, the root mean square difference to the explicit case divided by the value range over the simulation is given. All values are for the C4MIP A2 scenario (years 1700–2100), using the HILDA ocean component and the HRBM land component with standard temperature and carbon cycle sensitivities (coupled).

The explicit solution is only implemented for the piecewise constant approximation (Eq. B2) and the implicit solution for both the piecewise constant (Eq. B2) and the piecewise linear approximation (Eq. B3). Equations (B2) and (B3) are expressed in a common equation by substituting

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\displaystyle }{m_k}_n={m_k}_{n-\mathrm{1}}{p}_{mk}+f_n{p}_{fk}+f_{n-\mathrm{1}}{p}_{fk}^{\text{old}}.\end{array}$$

In the following, the implicit solution for the piecewise constant discretization is derived. Here, the fully implicit scheme for land and ocean exchange is discussed, but for stability it is only crucial to treat ocean uptake implicitly. The parameters of Eq. (B4) for this case are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{p}_{mk}=e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k}\\ {\displaystyle }& {\displaystyle }{p}_{fk}=a_k{\mathit{\tau }}_k(\mathrm{1}-e^{-\mathrm{\Delta }t/{\mathit{\tau }}_k})\\ \text{(B5)}& {\displaystyle }& {\displaystyle }{p}_{fk}^{\text{old}}=\mathrm{0}.\end{array}$$

First consider the equation system for carbon, assuming temperature to be known (or neglecting temperature dependence of model coefficients). Equation (B4) is applied to land carbon exchange for the constant approximation (Eq. B5),

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m_{\mathrm{L}}}_n=m_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}+\mathrm{\Delta }{f}_{\text{NPP}}\sum _k{p_{fk}}_{\mathrm{L}}\\ \text{(B6)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}=\sum _k{m_{\mathrm{L}k}}_{n-\mathrm{1}}{p_{mk}}_{\mathrm{L}}+{f_{\text{NPP}}}_{n-\mathrm{1}}\sum _k{p_{fk}}_{\mathrm{L}},\end{array}$$

where $m_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}$ is the land carbon stock obtained after one time step if NPP remained constant (“constant flux commitment”), and $\mathrm{\Delta }{f}_{\text{NPP}}=({f_{\text{NPP}}}_n-{f_{\text{NPP}}}_{n-\mathrm{1}})$ is the change in NPP over one time step. For ocean carbon uptake,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{m_{\mathrm{S}}}_n=m_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}+{f_{\mathrm{O}}}_n\sum _k{p_{fk}}_{\mathrm{O}}\\ \text{(B7)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}=\sum _k{m_{\mathrm{S}k}}_{n-\mathrm{1}}{p_{mk}}_{\mathrm{O}},\end{array}$$

where $m_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}$ is the value of m_S after one time step if ${f_{\mathrm{O}}}_n=\mathrm{0}$ (“zero-flux commitment”).

To solve the implicit system, the nonlinear parametrizations need to be linearized around t_n−1. Linearizing ocean surface CO₂ pressure as a function of surface ocean carbon and inserting in Eq. (4) yields

$$\begin{array}{ll}{\displaystyle }{f_{\mathrm{O}}}_n& {\displaystyle }\simeq k_{\mathrm{g}}{A}_{\mathrm{O}}\left({m_{\mathrm{A}}}_n-\mathit{\epsilon }{p}_{\mathrm{S},n-\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}\right)\\ \text{(B8)}& {\displaystyle }& {\displaystyle }+k_{\mathrm{g}}{A}_{\mathrm{O}}{\left.\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{p}_S^{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{d}{m}_{\mathrm{S}}}}\right|}_{n-\mathrm{1}}({m_{\mathrm{S}}}_{n-\mathrm{1}}-{m_{\mathrm{S}}}_n),\end{array}$$

where Eqs. (5) and (6) were used. Similarly, NPP as a function of atmospheric carbon is linearized,

$$\begin{array}{}\text{(B9)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{f_{\text{NPP}}}_n\simeq {\left.{\displaystyle \frac{\mathrm{d}{f}_{\text{NPP}}}{\mathrm{d}{m}_{\mathrm{A}}}}\right|}_{n-\mathrm{1}}({m_{\mathrm{A}}}_n-{m_{\mathrm{A}}}_{n-\mathrm{1}}),\end{array}$$

using Eq. (6).

The system is completed with the discretized budget equation (Eq. 1).

$$\begin{array}{}\text{(B10)}& {\displaystyle }{\displaystyle }{m_{\mathrm{A}}}_n={m_{\mathrm{A}}}_{n-\mathrm{1}}+\left(e_{n-\frac{\mathrm{1}}{\mathrm{2}}}-{f_{\mathrm{O}}}_n\right)\mathrm{\Delta }t-({m_{\mathrm{L}}}_n-{m_{\mathrm{L}}}_{n-\mathrm{1}})\end{array}$$

Here, e_n−½ is assumed to be known (though this only applies to the “forward” solution for atmospheric CO₂ from emissions; solving for emissions from CO₂ is also implemented in the model code).

After calculating the “committed” values ${m_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}}_n,{m_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}}_n$ from the model state at t_n−1, Eqs. (B7) through (B10) are solved:

$$\begin{array}{ll}\text{(B11)}& {\displaystyle }\mathrm{\Delta }{f}_{\text{NPP}}=& {\displaystyle }{\displaystyle \frac{{\left.\frac{\mathrm{d}{f}_{\text{NPP}}}{\mathrm{d}{m}_{\mathrm{A}}}\right|}_{n-\mathrm{1}}}{UV+W}}{\displaystyle }& {\displaystyle }\times ({m_{\mathrm{L}}}_{n-\mathrm{1}}-m_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}+\mathrm{\Delta }t{e}_{n-\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{\Delta }t{k}_{\mathrm{g}}{A}_{\mathrm{O}}\\ {\displaystyle }& {\displaystyle }\times (\mathit{\epsilon }{p}_{\mathrm{S},n-\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}-{m_{\mathrm{A}}}_{n-\mathrm{1}}+{\left.\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{p}_S^{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{d}{m}_{\mathrm{S}}}}\right|}_{n-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }\times [m_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}-{m_{\mathrm{S}}}_{n-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }+\sum _k{p_{fk}}_{\mathrm{O}}({\displaystyle \frac{{m_{\mathrm{L}}}_{n-\mathrm{1}}-m_{\mathrm{L}}^{\mathrm{c}\mathrm{\Delta }}}{\mathrm{\Delta }t}}+e_{n-\frac{\mathrm{1}}{\mathrm{2}}})\left]\right)),\end{array}$$

with the auxiliary variables

$$\begin{array}{}\text{(B12)}& {\displaystyle }& {\displaystyle }U=k_{\mathrm{g}}{A}_{\mathrm{O}}{\left.\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{p}_S^{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{d}{m}_{\mathrm{S}}}}\right|}_{n-\mathrm{1}}\sum _k{p_{fk}}_{\mathrm{O}}+\mathrm{1},\text{(B13)}& {\displaystyle }& {\displaystyle }V={\left.{\displaystyle \frac{\mathrm{d}{f}_{\text{NPP}}}{\mathrm{d}{m}_{\mathrm{A}}}}\right|}_{n-\mathrm{1}}\sum _k{p_{fk}}_{\mathrm{L}}+\mathrm{1},\text{(B14)}& {\displaystyle }& {\displaystyle }W=\mathrm{\Delta }t{k}_{\mathrm{g}}{A}_{\mathrm{O}},\end{array}$$

and, after inserting into Eq. (B6),

$$\begin{array}{ll}\text{(B15)}& {\displaystyle }{f_{\mathrm{O}}}_n=& {\displaystyle }{\displaystyle \frac{k_{\mathrm{g}}{A}_{\mathrm{O}}}{U+W}}({m_{\mathrm{A}}}_{n-\mathrm{1}}-\mathit{\epsilon }{p}_{\mathrm{S},n-\mathrm{1}}^{{\mathrm{CO}}_{\mathrm{2}}}-{\left.\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{p}_S^{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{d}{m}_{\mathrm{S}}}}\right|}_{n-\mathrm{1}}{\displaystyle }& {\displaystyle }\times (m_{\mathrm{S}}^{\mathrm{c}\mathrm{0}}-{m_{\mathrm{S}}}_{n-\mathrm{1}})-({m_{\mathrm{L}}}_n-{m_{\mathrm{L}}}_{n-\mathrm{1}})+\mathrm{\Delta }t{e}_{n-\frac{\mathrm{1}}{\mathrm{2}}}).\end{array}$$

The remaining variables are then calculated using Eqs. (B7) and (B10), whereby first the components ${m_k}_n$ are calculated as in Eq. (B4) and then summed. Finally, the nonlinear parametrizations are recalculated with the updated model state.

The order of these equations matters, as the updated variables are successively inserted into the following equations. The land part is solved first, and can be substituted by an explicit step or a separate model, while keeping the ocean step implicit.

An implicit time step is also implemented for calculating SAT from RF (again, solving RF from SAT is also implemented but not discussed here). RF(t_n) can be assumed as known, as atmospheric CO₂ is calculated first (i.e., no linearization necessary). Applying Eq. (B4) to temperature,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\Delta }{T}_n{c}_{\mathrm{S}}=\mathrm{\Delta }{T}^{\mathrm{c}\mathrm{\Delta }}{c}_{\mathrm{S}}+\mathrm{\Delta }{f}_{\mathrm{O}}^{\mathrm{H}}\sum _k{p_{fk}}_{\mathrm{O}}\\ \text{(B16)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }{T}^{\mathrm{c}\mathrm{\Delta }}=\sum _k\mathrm{\Delta }{T_k}_{n-\mathrm{1}}{p_{mk}}_{\mathrm{O}}+{f_{\mathrm{O}}^{\mathrm{H}}}_{n-\mathrm{1}}/c_{\mathrm{S}}\sum _k{p_{fk}}_{\mathrm{O}},\end{array}$$

where ΔT^cΔ is the committed temperature for constant heat flux to the ocean, and $\mathrm{\Delta }{f}_{\mathrm{O}}^{\mathrm{H}}={f_{\mathrm{O}}^{\mathrm{H}}}_n-{f_{\mathrm{O}}^{\mathrm{H}}}_{n-\mathrm{1}}$ is the change in heat flux over one time step. Equations (8), (9), and (B16) are solved for $f_{\mathrm{O}}^{\mathrm{H}}$,

$$\begin{array}{}\text{(B17)}& {\displaystyle }{\displaystyle }{f_{\mathrm{O}}^{\mathrm{H}}}_n={\displaystyle \frac{{\text{RF}}_n-\frac{{\text{RF}}_{\mathrm{2}\times }}{\mathrm{\Delta }{T}_{\mathrm{2}\times }}\mathrm{\Delta }{T}^{\mathrm{c}\mathrm{\Delta }}+f_{Hn-\mathrm{1}}{\sum }_k{p_{fk}}_{\mathrm{O}}\frac{{\text{RF}}_{\mathrm{2}\times }}{\mathrm{\Delta }{T}_{\mathrm{2}\times }{c}_{\mathrm{S}}}}{\frac{{\text{RF}}_{\mathrm{2}\times }}{\mathrm{\Delta }{T}_{\mathrm{2}\times }{c}_{\mathrm{S}}}{\sum }_k{p_{fk}}_{\mathrm{O}}+a_{\mathrm{O}}/A_{\mathrm{O}}}}.\end{array}$$

Temperature change ΔT_n then follows from Eq. (B16).

The case of piecewise linear approximation (Eq. B3) differs from the piecewise constant one (Eq. B2) only in a nonzero contribution of f_n−1 and a slightly different budget equation,

$$\begin{array}{ll}\text{(B18)}& {\displaystyle }{m_{\mathrm{A}}}_n=& {\displaystyle }{m_{\mathrm{A}}}_{n-\mathrm{1}}{\displaystyle }& {\displaystyle }+\left(e_{n-\frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{{f_{\mathrm{O}}}_n+{f_{\mathrm{O}}}_{n-\mathrm{1}}}{\mathrm{2}}}\right)\mathrm{\Delta }t-({m_{\mathrm{L}}}_n-{m_{\mathrm{L}}}_{n-\mathrm{1}}).\end{array}$$

The first difference merely changes the calculation of committed changes, and only the second difference affects the solution of the implicit time step. In practice, however, this can be neglected without loss of accuracy, and thus Eqs. (B11)–(B15) and (B17) are also used to solve the piecewise linear system (while Eq. B18 is used to close the budget).

B3 Temperature-dependent parameters

BernSCM allows for temperature-dependent model parameters for IRF-based substitute models. This generalization of the IRF approach is possible using a box model form (Sect. 2.3). Currently, temperature-dependent coefficients and timescales are implemented for the HRBM land biosphere substitute model (Appendix A2).

BernSCM updates any temperature-dependent model parameters by approximating the current temperature ΔT_n by the committed temperature ΔT^cΔ as defined in Eq. (B16). Accuracy is further improved by substituting ΔT^cΔ for ΔT_n in evaluating Eq. (B8) with temperature-dependent parametrizations.