SURFER v2.0: a flexible and simple model linking anthropogenic CO₂ emissions and solar radiation modification to ocean acidification and sea level rise

We present SURFER, a novel reduced model for estimating the impact of CO₂ emissions and solar radiation modification options on sea level rise and ocean acidification over timescales of several thousands of years. SURFER has been designed for the analysis of CO₂ emission and solar radiation modification policies, for supporting the computation of optimal (CO₂ emission and solar radiation modification) policies and for the study of commitment and responsibility under uncertainty. The model is based on a combination of conservation laws for the masses of atmospheric and oceanic carbon and for the oceanic temperature anomalies, and of ad-hoc parameterisations for the different sea level rise contributors: ice sheets, glaciers and ocean thermal expansion. It consists of 9 loosely coupled ordinary differential equations, is understandable, fast and easy to modify and calibrate. It reproduces the results of more sophisticated, high-dimensional earth system models on timescales up to millennia.

1 Introduction

Carbon emissions in the following decades will have a significant impact on the sea level and on the acidity of the oceans for millennia (e.g. Clark et al., 2016; Van Breedam et al., 2020), mainly because of the long residence time of CO₂ in the atmosphere–ocean system (Archer, 2005; Archer et al., 2009). While reducing CO₂ emissions is necessary, different geoengineering methods have been proposed to reduce some of the negative effects of climate change and complement the efforts of CO₂ emission reduction (Lawrence et al., 2018). One approach to such climate intervention is that of solar radiation modification, which proposes to reflect some of the incoming solar radiation to cool Earth’s surface. Solar radiation modification, in particular stratospheric aerosol intervention, acts on much shorter timescales than the residence time of CO₂ in the atmosphere.

Therefore, assigning a value judgement to climate mitigation policies for the next decades necessarily has to account for impacts (e.g. in terms of benefits and damages due to sea level rise, crossing planetary boundaries (Rockström et al., 2009; Steffen et al., 2015) and perhaps solar radiation modification and CO₂ emission reduction costs) over much longer timescales than the next decades.

Failing to do so may lead to short-sighted policies that commit upcoming generations to unmanageable impacts or that severely shrinks their range of viable CO₂ emission and/or solar radiation modification options (Clark et al., 2016; Mengel et al., 2018; Nauels et al., 2019).

Assessing the value of climate mitigation policies for the next decades over millennia requires accounting for epistemic (Shepherd et al., 2018; Shepherd, 2019) model uncertainties but also, and perhaps more importantly, for uncertainties about how these policies will actually be implemented: we know that decisions in the matter of greenhouse gas abatements may be implemented with delays and that solar radiation modification options for mitigating the impacts of high CO₂ concentrations in the atmosphere are not free from risks (see Gardiner, 2010; Moreno-Cruz and Keith, 2013; Robock, 2016, 2020; Helwegen et al., 2019; Zarnetske et al., 2021). Also, assessing the value of climate mitigation policies for the next decades based on sea level rise and ocean acidification necessarily has to account for uncertainties about the (anthropogenic) forcings to be expected after, say, 2100.

Depending on the methodology applied and on the level of guarantees (correctness) required, computing “best” climate policies under uncertainty and estimating measures of responsibility for specific climate decisions (Botta et al., 2021) can be more or less computationally expensive than assessing the value of a given policy under the same uncertainties. As a consequence, in the presence of uncertainty, various interesting assessments require models that are easy to understand and modify, and fast to apply. Such assessments include (among others) the analysis of CO₂ emission policies, the computation of optimal policies, the assessment of commitment, responsibility and safe operating spaces (Heitzig et al., 2016), e.g. in terms of planetary boundaries (Rockström et al., 2009).

1.1 What this paper is about

This paper presents SURFER, a tool for estimating the Sea level Uprise Response under Forcings of Emissions and solar Radiation modification, that has been designed to meet all the above requirements. SURFER is a simple carbon–climate–sea level rise model based on a novel combination of conservation laws for the masses of atmospheric and oceanic carbon and for the oceanic temperature anomalies, and of ad hoc parameterisations for the different sea level rise contributors.

It consists of nine loosely coupled ordinary differential equations (ODEs) that are easy to understand and modify. Model simulations over 10 000 years (obtained with standard numerical approximations for stiff ODEs) typically run in less than 0.005 s on a standard laptop, see Sect. 2.5.

The model is easy to calibrate and reproduces the results of high-dimensional earth system models of intermediate complexity (EMICs) and Earth system models (ESMs), capturing well the responses of global mean surface temperature anomaly, atmospheric carbon concentration, sea level rise and ocean acidification to anthropogenic forcing.

The main focus of this paper is twofold: on the one hand, we discuss and motivate the model equations and parameters, and contrast these equations and parameters to those found in similar reduced models and climate emulators. This is done in detail in Sects. 2 and 4. The aim of Sect. 2 is to establish full model transparency and, hopefully, understandability. Developing a tool for estimating the impacts of CO₂ emissions and solar radiation modification measures on sea level rise and on the acidification of the oceans over millennia necessarily means making a number of choices and compromises. These choices and compromises are motivated by the intended applications and in Sect. 4 we provide recommendations for extending the model of Sect. 2 to applications that go beyond those targeted by SURFER.

On the other hand, we provide numerical evidence that, for the intended applications, SURFER can indeed reproduce the results of more sophisticated, high-dimensional earth system models. This is done in Sect. 3.

1.2 What the paper is not about

Before turning to the specification of the model equations of SURFER, let us shortly discuss what this paper is not about.

We have argued that one intended application of SURFER is that of assessing the value of CO₂ emission policies for the next decades based on sea level rise and ocean acidification under different kinds of uncertainty. We do not provide an example of such an application here. Besides the specification of an emission policy, this would also require specifying

• the uncertainties that affect the specific problem at stake,

• a function for measuring the value of possible trajectories (sometimes called a metric) and, most importantly,

• a measure for aggregating probability distributions (often the expected value measure),

for the specific problem at stake. This would go well beyond the scope of this paper but we refer the reader to (Botta et al., 2018, 2021; Helwegen et al., 2019; Carlino et al., 2020) for a discussion of the steps involved in setting up a stylised problem. Similarly, we do not apply SURFER to the computation of optimal policies (Botta et al., 2018; Helwegen et al., 2019; Carlino et al., 2020) or to the study of commitment and responsibility (Martínez Montero et al., 2022a; Botta et al., 2021) under uncertainty here. Instead, we will do so in a dedicated companion paper.

A final caveat is at place: SURFER is not meant to be a “better” model than the reference ESMs it is intended to emulate in terms of climate metrics, and has not been calibrated to minimise a well-defined distance from a set of “trusted” models (see Sect. 3). Rather, SURFER is a tool that trades model complexity for speed and understandability. Both properties are crucial for the intended applications: speed is needed to tackle the computational complexity associated with uncertainty; and, under uncertainty, nothing is worse than applying a tool that one does not fully understand.

2 Model specification

SURFER consists of a system of nine ordinary differential equations for the masses of carbon in four different reservoirs (atmosphere, upper ocean layer, deeper ocean layer and land), the temperature anomalies of two different reservoirs (upper ocean layer and deeper ocean layer), the volume of Greenland and Antarctic ice sheets and sea level rise related to glaciers. Sea level rise due to the ocean thermal expansion is also taken into account through a parameterisation in terms of the ocean temperature anomalies.

Two external forcings drive the system: anthropogenic CO₂ emissions and solar radiation modification in the form of SO₂ injections.

The CO₂ emission rate, through the accumulation of CO₂ in the atmosphere and due to the greenhouse effect, leads to long-lasting temperature increases, while the SO₂ injections, through reflecting some of the incoming solar radiation, are responsible for fast but short-lasting temperature decreases. We refer the readers to (Moreno-Cruz and Keith, 2013; Helwegen et al., 2019; Robock, 2020) for more information on this form of geoengineering and to (Visioni et al., 2021) for the latest results from the Geoengineering Model Inter-comparison Project. The CO₂ emissions are the source in the carbon cycle sub-model, which evolves the amount of carbon in the considered carbon reservoirs. The carbon concentration in the atmosphere and the SO₂ injections are the driving forces of the climate sub-model, which evolves the temperature anomalies in two thermal reservoirs. Finally, the temperature anomaly of the upper ocean layer (assumed in equilibrium with the atmosphere and land) forces the melting of the ice sheets and glaciers, causing sea level to rise, while both ocean layers contribute to sea level rise due to the ocean's thermal expansion. Figure 1 provides a conceptual diagram for the model. We have not included other interactions between the sub-models but we address the most relevant feedbacks in Sect. 4.

Figure 1 Conceptual diagram of SURFER. The state variables are indicated by the boxes, interactions and sources are depicted by black and dark orange arrows respectively.

2.1 Carbon cycle model

The carbon cycle model is based on BEAM, a simple carbon cycle model developed by Glotter et al. (2014) for economic and policy analyses. BEAM stands for “Bolin and Eriksson Adjusted Model” in acknowledgement of Bolin and Eriksson (1959). Since we modified the model by Glotter et al. (2014) by including an extra carbon reservoir, a land reservoir¹, along with minor modifications to the original equations, we proceed with a full presentation of the carbon cycle model.

The equations for the atmosphere (A), upper ocean layer (U), deeper ocean layer (D) and land (L) carbon reservoirs read as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{A}}}{\mathrm{d}t}}& {\displaystyle }=E-F_{\mathrm{A}\to \mathrm{U}}-F_{\mathrm{A}\to \mathrm{L}},\text{(2)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{U}}}{\mathrm{d}t}}& {\displaystyle }=F_{\mathrm{A}\to \mathrm{U}}-F_{\mathrm{U}\to \mathrm{D}},\text{(3)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{D}}}{\mathrm{d}t}}& {\displaystyle }=F_{\mathrm{U}\to \mathrm{D}},\text{(4)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{L}}}{\mathrm{d}t}}& {\displaystyle }=F_{\mathrm{A}\to \mathrm{L}},\end{array}$$

where M_i is the mass of carbon in reservoir i, F_i→j the net carbon flux from reservoir i to reservoir j and E is the anthropogenic carbon emission rate. As part of the land reservoir, we consider only soils and vegetation, ignoring carbon in permafrost and fossil fuel reserves. Sinks and sources associated with carbon outgassing, weathering and sediment burial are ignored because they are of secondary importance at the timescales considered here (10 to 5000 years).

2.1.1  F_A→U

Modelling the carbon flux between the atmosphere and the ocean relies on fundamental ocean carbonate chemistry which we now summarise (see chap. 8 of textbook by Sarmiento and Gruber, 2006 for a deeper treatment).

When CO₂ in the atmosphere goes into the ocean it undergoes a series of chemical reactions:

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{gas}\right)}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*},\text{(R2)}& {\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\rightleftharpoons {\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{H}}^{+},\text{(R3)}& {\displaystyle }{\mathrm{HCO}}_{\mathrm{3}}^{-}\rightleftharpoons {\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}+\mathrm{2}{\mathrm{H}}^{+},\end{array}$$

where H₂CO${}_{\mathrm{3}}^{*}$ represents a mix of aqueous carbon dioxide, CO_2(aqueous), and carbonic acid, H₂CO₃². The distribution between the three carbon species, H₂CO${}_{\mathrm{3}}^{*}$, HCO${}_{\mathrm{3}}^{-}$ (bicarbonate), and CO${}_{\mathrm{3}}^{\mathrm{2}-}$ (carbonate), is fast with respect to the ocean's circulation timescale, and hence equilibrium is assumed. The equilibrium distribution relations,

$$\begin{array}{}\text{(5)}& K_{\mathrm{1}}={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right]}},K_{\mathrm{2}}={\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}},\end{array}$$

are dictated by the ocean's acidity, quantified by the proton concentration [H⁺]. K₁ and K₂ are dissociation constants and [H⁺], measured in moles per kilogram, relates to the ocean pH as

$$\begin{array}{}\text{(6)}& \mathrm{pH}=-{\log }_{\mathrm{10}}\left[{\mathrm{H}}^{+}\right].\end{array}$$

The total dissolved inorganic carbon (DIC) can then be written as

$$\begin{array}{ll}\text{(7)}& {\displaystyle }\mathrm{DIC}& {\displaystyle }=\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]{\displaystyle }& {\displaystyle }=\left(\mathrm{1}+{\displaystyle \frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]}}+{\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\right)\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right].\end{array}$$

Figure 2 Carbon species fractional contribution to DIC for varying pH.

In Fig. 2, we show the fractional contributions of the different carbon species to DIC for varying pH. For the present pH value of around 8, we can see in Fig. 2 that bicarbonate is the dominant carbon species in the ocean. From the chemical reactions (R1), (R2) and (R3), we see that when CO₂ dissolves in the ocean, hydrogen ions are released and ocean acidifies. This in turn means that the proportion of carbonate decreases and that of H₂CO${}_{\mathrm{3}}^{*}$ increases. The alkalinity, however, defined as the excess of bases over acids,

$$\begin{array}{ll}\text{(8)}& {\displaystyle }Q& {\displaystyle }=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]+\left[{\mathrm{OH}}^{-}\right]-\left[{\mathrm{H}}^{+}\right]+\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}\right]{\displaystyle }& {\displaystyle }+\text{minor bases},\end{array}$$

does not change with those reactions. Since no other reactions are accounted for in this carbon cycle model the alkalinity is constant. This assumption will lead to stronger than expected acidification on long timescales (∼4000 years) in which calcium carbonate production, dissolution and burial (not accounted for here) are significant. As it is usual practice, we will approximate the alkalinity³ by its dominant terms, that is by the carbonate alkalinity,

$$\begin{array}{}\text{(9)}& Q\approx \left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]=\left({\displaystyle \frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]}}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\right)\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right].\end{array}$$

The flow of CO₂ between the atmosphere and upper ocean layer is proportional to the difference in CO₂ partial pressures,

$$\begin{array}{}\text{(10)}& F_{\mathrm{A}\to \mathrm{U}}\propto \left(p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}-p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{U}}\right),\end{array}$$

and by writing the partial pressures as proportional to the corresponding carbon masses,

$$\begin{array}{}\text{(11)}& F_{\mathrm{A}\to \mathrm{U}}=-k_{\mathrm{A}\to \mathrm{U}}{M}_{\mathrm{A}}+k_{\mathrm{U}\to \mathrm{A}}{M}_{\mathrm{U}}^{\prime },\end{array}$$

where $p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{U}}$ refers to the partial pressure of carbon in the form H₂CO${}_{\mathrm{3}}^{*}$, and $M_{\mathrm{U}}^{\prime }$ to the corresponding carbon mass. The parameters k_i→j are the transport coefficients from reservoir i to reservoir j.

The equilibrium concentration of H₂CO${}_{\mathrm{3}}^{*}$ in the ocean, corresponding to an atmospheric CO₂ partial pressure $p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}$ can be determined through the CO₂ solubility constant,

$$\begin{array}{}\text{(12)}& K_{\mathrm{0}}={\displaystyle \frac{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right]}{p_{{\mathrm{CO}}_{\mathrm{2}}}}},\end{array}$$

where $p_{{\mathrm{CO}}_{\mathrm{2}}}$ is written without a reservoir index because, when in equilibrium, atmospheric and upper ocean have the same CO₂ partial pressure.

2.1.2  F_U→D

The exchange of carbon between the two ocean layers is ruled by oceanic currents and therefore depends on the total dissolved inorganic carbon in each layer,

$$\begin{array}{}\text{(13)}& F_{\mathrm{U}\to \mathrm{D}}=k_{\mathrm{U}\to \mathrm{D}}{M}_{\mathrm{U}}-k_{\mathrm{D}\to \mathrm{U}}{M}_{\mathrm{D}}.\end{array}$$

This is in contrast with the carbon exchange between the upper ocean to atmosphere which depends on the upper ocean carbon concentration in the form of H₂CO${}_{\mathrm{3}}^{*}$. Now we can write the carbon cycle equations as

$$\begin{array}{ll}\text{(14a)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{A}}}{\mathrm{d}t}}& {\displaystyle }=E-k_{\mathrm{A}\to \mathrm{U}}{M}_{\mathrm{A}}+k_{\mathrm{U}\to \mathrm{A}}{M}_{\mathrm{U}}^{\prime }-F_{\mathrm{A}\to \mathrm{L}},\text{(14b)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{U}}}{\mathrm{d}t}}& {\displaystyle }=k_{\mathrm{A}\to \mathrm{U}}{M}_{\mathrm{A}}-k_{\mathrm{U}\to \mathrm{A}}{M}_{\mathrm{U}}^{\prime }-k_{\mathrm{U}\to \mathrm{D}}{M}_{\mathrm{U}}{\displaystyle }& {\displaystyle }+k_{\mathrm{D}\to \mathrm{U}}{M}_{\mathrm{D}},\\ \text{(14c)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{D}}}{\mathrm{d}t}}& {\displaystyle }=k_{\mathrm{U}\to \mathrm{D}}{M}_{\mathrm{U}}-k_{\mathrm{D}\to \mathrm{U}}{M}_{\mathrm{D}},\text{(14d)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{L}}}{\mathrm{d}t}}& {\displaystyle }=F_{\mathrm{A}\to \mathrm{L}}.\end{array}$$

Following Bolin and Eriksson (1959) we assume that the four reservoirs were in equilibrium at pre-industrial times (with E(t_PI)=0) and we examine the equilibrium equations,

$$\begin{array}{ll}\text{(15a)}& {\displaystyle }\mathrm{0}=& {\displaystyle }-k_{\mathrm{A}\to \mathrm{U}}{M}_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)+k_{\mathrm{U}\to \mathrm{A}}{M}_{\mathrm{U}}^{\prime }\left(t_{\mathrm{PI}}\right)-F_{\mathrm{A}\to \mathrm{L}}\left(t_{\mathrm{PI}}\right),\text{(15b)}& {\displaystyle }\mathrm{0}=& {\displaystyle }{k}_{\mathrm{A}\to \mathrm{U}}{M}_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)-k_{\mathrm{U}\to \mathrm{A}}{M}_{\mathrm{U}}^{\prime }\left(t_{\mathrm{PI}}\right)-k_{\mathrm{U}\to \mathrm{D}}{M}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right){\displaystyle }& {\displaystyle }+k_{\mathrm{D}\to \mathrm{U}}{M}_{\mathrm{D}}\left(t_{\mathrm{PI}}\right),\\ \text{(15c)}& {\displaystyle }\mathrm{0}=& {\displaystyle }{k}_{\mathrm{U}\to \mathrm{D}}{M}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)-k_{\mathrm{D}\to \mathrm{U}}{M}_{\mathrm{D}}\left(t_{\mathrm{PI}}\right),\text{(15d)}& {\displaystyle }\mathrm{0}=& {\displaystyle }{F}_{\mathrm{A}\to \mathrm{L}}\left(t_{\mathrm{PI}}\right).\end{array}$$

These allow to relate transport coefficients,

$$\begin{array}{}\text{(16)}& {\displaystyle }{k}_{\mathrm{A}\to \mathrm{U}}& {\displaystyle }=k_{\mathrm{U}\to \mathrm{A}}{\displaystyle \frac{M_{\mathrm{U}}^{\prime }\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}},\text{(17)}& {\displaystyle }{k}_{\mathrm{U}\to \mathrm{D}}& {\displaystyle }=k_{\mathrm{D}\to \mathrm{U}}{\displaystyle \frac{M_{\mathrm{D}}\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)}},\end{array}$$

where we can further write

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{M}_{\mathrm{U}}^{\prime }=\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}\right]W_{\mathrm{U}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }{M}_{\mathrm{A}}=\text{moles of}{\mathrm{CO}}_{\mathrm{2}}\text{in atmosphere}\times {\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}},\\ {\displaystyle }& {\displaystyle }{M}_{(\mathrm{U},\mathrm{D})}={\mathrm{DIC}}_{(\mathrm{U},\mathrm{D})}{W}_{(\mathrm{U},\mathrm{D})}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}},\end{array}$$

where W_U and W_D stand for the whole mass of the upper and deeper ocean layers and can be approximated by

$$\begin{array}{}\text{(18)}& W_{\mathrm{U}}\approx m_{\mathrm{O}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{W}}{\displaystyle \frac{h_{\mathrm{U}}}{h_{\mathrm{U}}+h_{\mathrm{D}}}},W_{\mathrm{D}}\approx m_{\mathrm{O}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{W}}{\displaystyle \frac{h_{\mathrm{D}}}{h_{\mathrm{U}}+h_{\mathrm{D}}}},\end{array}$$

with m_O the moles of water in the ocean, ${\stackrel{\mathrm{‾}}{m}}_{\mathrm{W}}$ the molar mass of H₂O, ${\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}$ the molar mass of carbon and h_U and h_D the thicknesses of the ocean layers. Considering the equilibrium solubility relation (Eq. 12) and

$$p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}=\mathrm{1}\left(\mathrm{atm}\right){\displaystyle \frac{\text{moles of}{\mathrm{CO}}_{\mathrm{2}}\text{in atmosphere}}{\text{moles in atmosphere}}},$$

we arrive at

$$\begin{array}{}\text{(19)}& {\displaystyle }{k}_{\mathrm{A}\to \mathrm{U}}& {\displaystyle }=k_{\mathrm{U}\to \mathrm{A}}{\displaystyle \frac{M_{\mathrm{U}}^{\prime }\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}}=k_{\mathrm{U}\to \mathrm{A}}{\displaystyle \frac{W_{\mathrm{U}}{K}_{\mathrm{0}}}{m_{\mathrm{A}}}},\text{(20)}& {\displaystyle }{k}_{\mathrm{U}\to \mathrm{D}}& {\displaystyle }=k_{\mathrm{D}\to \mathrm{U}}{\displaystyle \frac{M_{\mathrm{D}}\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)}}=k_{\mathrm{D}\to \mathrm{U}}{\mathit{\delta }}_{\mathrm{DIC}}{\displaystyle \frac{h_{\mathrm{D}}}{h_{\mathrm{U}}}},\end{array}$$

where m_A are the moles of air in the atmosphere and

$$\begin{array}{}\text{(21)}& {\mathit{\delta }}_{\mathrm{DIC}}={\displaystyle \frac{{\mathrm{DIC}}_{\mathrm{D}}\left(t_{\mathrm{PI}}\right)}{{\mathrm{DIC}}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)}},\end{array}$$

where DIC_U(t_PI) and DIC_D(t_PI) are the pre-industrial DIC concentration in upper and lower ocean layers. The parameter δ_DIC specifies the DIC gradient between the two ocean layers and effectively accounts for the biological and carbonate pumps.

The next step is to express the carbon mass in H₂CO${}_{\mathrm{3}}^{*}$ form, $M_{\mathrm{U}}^{\prime }$, in Eq. (14a) and (14b) as a function of the state variables. We begin by

$${\displaystyle \frac{M_{\mathrm{U}}^{\prime }}{M_{\mathrm{U}}}}={\displaystyle \frac{[{\mathrm{H}}_{\mathrm{2}}{\mathrm{CO}}_{\mathrm{3}}^{*}{]}_{\mathrm{U}}}{{\mathrm{DIC}}_{\mathrm{U}}}}={\left(\mathrm{1}+{\displaystyle \frac{K_{\mathrm{1}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{U}}}}+{\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}_{\mathrm{U}}^{\mathrm{2}}}}\right)}^{-\mathrm{1}}.$$

Using the definitions of DIC and carbonate alkalinity (Q in Eq. 9), and the relation of DIC_i with carbon mass M_i, [H⁺]_i can be solved for in terms of M_i:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}[{\mathrm{H}}^{+}{]}_i& ={\displaystyle \frac{K_{\mathrm{1}}}{\mathrm{2}{\stackrel{\mathrm{̃}}{Q}}_i}}(\sqrt{(M_i-{\stackrel{\mathrm{̃}}{Q}}_i{)}^{\mathrm{2}}-\mathrm{4}{\displaystyle \frac{K_{\mathrm{2}}}{K_{\mathrm{1}}}}{\stackrel{\mathrm{̃}}{Q}}_i({\stackrel{\mathrm{̃}}{Q}}_i-\mathrm{2}{M}_i)}\\ & +(M_i-{\stackrel{\mathrm{̃}}{Q}}_i)),\end{array}\end{array}$$

with i=U or D and where

$${\stackrel{\mathrm{̃}}{Q}}_i=Q_i{W}_i{\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}},$$

is the carbon mass corresponding to the carbonate alkalinity. The factor tracking the ocean's buffer capacity becomes

$$\begin{array}{}\text{(23)}& \begin{array}{rl}B& \left(M_{\mathrm{U}}\right)\equiv {\displaystyle \frac{M_{\mathrm{U}}^{\prime }}{M_{\mathrm{U}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{U}}}{\mathrm{2}{M}_{\mathrm{U}}}}\\ & +{\displaystyle \frac{K_{\mathrm{1}}\sqrt{(M_{\mathrm{U}}-{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{U}}{)}^{\mathrm{2}}-\mathrm{4}\frac{K_{\mathrm{2}}}{K_{\mathrm{1}}}{\stackrel{\mathrm{̃}}{Q}}_{\mathrm{U}}({\stackrel{\mathrm{̃}}{Q}}_{\mathrm{U}}-\mathrm{2}{M}_{\mathrm{U}})}-\mathrm{4}{M}_{\mathrm{U}}{K}_{\mathrm{2}}}{\mathrm{2}{M}_{\mathrm{U}}(K_{\mathrm{1}}-\mathrm{4}{K}_{\mathrm{2}})}}.\end{array}\end{array}$$

2.1.3  F_A→L

Now we turn to the equation for the land reservoir. The proposed equation is not process based, it is instead based on the output of the Zero Emissions Commitment Model Inter-comparison Project (ZECMIP) (Jones et al., 2019; MacDougall et al., 2020) experiments by different EMICs and ESMs.

We analysed the output of the ZECMIP experiments B1 and B3 in Jones et al. (2019); MacDougall et al. (2020), and observed that the land carbon anomaly relaxes to a value proportional to the atmospheric carbon anomaly after typically 4 to 6 decades. This behaviour can be captured by the following relation:

$$\begin{array}{}\text{(24)}& \mathit{\delta }{M}_{\mathrm{L}}\left(t\right)={\mathit{\alpha }}_{\mathrm{L}}\mathit{\delta }{M}_{\mathrm{A}}\left(t\right),\end{array}$$

where $\mathit{\delta }{M}_i\left(t\right)=M_i\left(t\right)-M_i\left(t_{\mathrm{PI}}\right)$. The values of α_L approached by the different models in B1 experiments (lower cumulated emissions) tend to be higher than the ones approached by the higher cumulated emissions experiment B3. We noticed that the quantity,

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathit{\delta }{M}_{\mathrm{L}}\left(t\right)}{\mathit{\delta }{M}_{\mathrm{A}}\left(t\right)}}{\displaystyle \frac{M_{\mathrm{A}}\left(t\right)}{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}},\end{array}$$

tends to approach a model-dependent constant value, say β_L, which is independent of the cumulated emissions after 4 to 6 decades. Based on these observations we propose the following equation for the land carbon anomaly δM_L,

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{d\mathit{\delta }{M}_{\mathrm{L}}}{\mathrm{d}t}}=k_{\mathrm{A}\to \mathrm{L}}\left({\mathit{\beta }}_{\mathrm{L}}{\displaystyle \frac{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{A}}}}\mathit{\delta }{M}_{\mathrm{A}}-\mathit{\delta }{M}_{\mathrm{L}}\right).\end{array}$$

With this equation δM_L relaxes to an equilibrium value proportional to the ratio $\frac{\mathit{\delta }{M}_{\mathrm{A}}}{M_{\mathrm{A}}}$. This dependency can be interpreted as the result of two competing processes: CO₂ fertilisation (∝δM_A) and an enhanced bacterial respiration due to climate change (∝M_A).

The final equations for the carbon cycle can now be written as

$$\begin{array}{}\text{(27a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{M}_{\mathrm{A}}}{\mathrm{d}t}}& =E\left(t\right)-k_{\mathrm{A}\to \mathrm{U}}\left(M_{\mathrm{A}}-{\displaystyle \frac{m_{\mathrm{A}}}{W_{\mathrm{U}}{K}_{\mathrm{0}}}}B\left(M_{\mathrm{U}}\right)M_{\mathrm{U}}\right)\\ & -k_{\mathrm{A}\to \mathrm{L}}\left({\mathit{\beta }}_{\mathrm{L}}{M}_{\mathrm{A}}\right(t_{\mathrm{PI}})\left(\mathrm{1}-{\displaystyle \frac{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{A}}}}\right)\\ & -(M_{\mathrm{L}}-M_{\mathrm{L}}(t_{\mathrm{PI}}\left)\right)),\end{array}\text{(27b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{M}_{\mathrm{U}}}{\mathrm{d}t}}& =k_{\mathrm{A}\to \mathrm{U}}\left(M_{\mathrm{A}}-{\displaystyle \frac{m_{\mathrm{A}}}{W_{\mathrm{U}}{K}_{\mathrm{0}}}}B\left(M_{\mathrm{U}}\right)M_{\mathrm{U}}\right)\\ & -k_{\mathrm{U}\to \mathrm{D}}\left(M_{\mathrm{U}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{DIC}}}}{\displaystyle \frac{h_{\mathrm{U}}}{h_{\mathrm{D}}}}{M}_{\mathrm{D}}\right),\end{array}\text{(27c)}& {\displaystyle \frac{\mathrm{d}{M}_{\mathrm{D}}}{\mathrm{d}t}}=k_{\mathrm{U}\to \mathrm{D}}\left(M_{\mathrm{U}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{DIC}}}}{\displaystyle \frac{h_{\mathrm{U}}}{h_{\mathrm{D}}}}{M}_{\mathrm{D}}\right),\text{(27d)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{M}_{\mathrm{L}}}{\mathrm{d}t}}& =k_{\mathrm{A}\to \mathrm{L}}\left({\mathit{\beta }}_{\mathrm{L}}{M}_{\mathrm{A}}\right(t_{\mathrm{PI}})\left(\mathrm{1}-{\displaystyle \frac{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}{M_{\mathrm{A}}}}\right)\\ & -(M_{\mathrm{L}}-M_{\mathrm{L}}(t_{\mathrm{PI}}\left)\right)).\end{array}\end{array}$$

Equation (27a–d) is very similar to the ones presented by Glotter et al. (2014) but with the following important differences:

1. There is a land carbon reservoir. This update to the Glotter et al. (2014) model improves the agreement with most recent results from EMICs and ESMs.

2. The ocean buffer factor is explicitly written in terms of M_U to highlight the non-linear nature of the model.

3. The relation between the transport coefficients between the two ocean layers depends not only on the ratio of the thickness of the layers (δ) but also on the ratio of their pre-industrial concentration of dissolved inorganic carbon (δ_DIC). This allows for an equilibrium solution in which the dissolved inorganic carbon concentration is different in the upper and lower layers, which is known to be the case due to the soft tissue and carbonate pumps.

The presented carbon cycle equations for atmosphere and ocean (and also the ones in Glotter et al., 2014) are similar to the ones considered in DICE, the Dynamic Integrated model of Climate and the Economy (Nordhaus, 1992). The big difference between the two is that the upper ocean buffer factor B is considered to be constant in DICE while in SURFER it evolves with the ocean acidification. By including the non-linearities due to ocean carbonate chemistry, SURFER's carbon cycle, as the one by Glotter et al. (2014), captures the fact that as the ocean takes in CO₂ from the atmosphere it becomes a worse sink for future CO₂ intake. This enables the correct tracking of carbon concentrations up to timescales of several thousands of years, which is impossible with a linear model like DICE. One of the main objectives of the present contribution is to provide a model of sea level rise caused by ice sheet melting which is a slow process lasting several thousands of years. As explained by Archer et al. (2009), a linear carbon cycle is inadequate for such long-term purposes.

Another benefit of SURFER's carbon cycle is the tracking of the pH and the concentrations of the different carbon species in the ocean; this can be done a posteriori, using the obtained M_U(t) together with Eqs. (22), (6), (9) and (5). SURFER can thus be of use for policy analyses that deal with ocean acidification.

Ocean acidification destabilises marine ecosystems making it more difficult for shellfish and corals to grow. As such, ocean acidification is one of the nine identified planetary boundaries (Rockström et al., 2009; Steffen et al., 2015). This planetary boundary has however not been quantified in terms of pH. It is instead given in terms of the aragonite saturation state which can be approximated as

$$\begin{array}{}\text{(28)}& {\mathrm{\Omega }}_{\mathrm{ar}}\approx {\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}{]}_{\text{saturation ar}}}},\end{array}$$

where [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] is the carbonate concentration in the ocean and [CO${}_{\mathrm{3}}^{\mathrm{2}-}$]_{saturation ar} is the carbonate concentration at which aragonite is saturated⁴. Aragonite is the most vulnerable form of calcium carbonate and its saturation state relates in a more straightforward way to some forms of marine life than pH. A value of Ω_ar<1 means that aragonite dissolves but organisms struggle to grow and thrive already at bigger values of Ω_ar. The planetary boundary is set at 80 % of the pre-industrial value which was Ω_ar(t_PI)=3.44.

2.2 Climate model with solar radiation modification

SURFER's climate sub-model is a linear 2-box model, similar to those in Gregory (2000); Held et al. (2010), for the evolution of the upper and deeper ocean temperature anomalies δT_U and δT_D, respectively (measured in ^∘C),

$$\begin{array}{}\text{(29)}& {\displaystyle }{c}_{\mathrm{vol}}{h}_{\mathrm{U}}{\displaystyle \frac{d\mathit{\delta }{T}_{\mathrm{U}}}{\mathrm{d}t}}=\mathit{\beta }\left({\displaystyle \frac{F(M_{\mathrm{A}},I)}{\mathit{\beta }}}-\mathit{\delta }{T}_{\mathrm{U}}\right)-\mathit{\gamma }(\mathit{\delta }{T}_{\mathrm{U}}-\mathit{\delta }{T}_{\mathrm{D}}),\text{(30)}& {\displaystyle }{c}_{\mathrm{vol}}{h}_{\mathrm{D}}{\displaystyle \frac{d\mathit{\delta }{T}_{\mathrm{D}}}{\mathrm{d}t}}=\mathit{\gamma }(\mathit{\delta }{T}_{\mathrm{U}}-\mathit{\delta }{T}_{\mathrm{D}}).\end{array}$$

The atmosphere is assumed to be in thermal equilibrium with the upper ocean layer. Due to this, and that the upper ocean heat capacity is much bigger than that of the atmosphere, the atmosphere is not explicitly part of the climate sub-model and the radiative forcings are applied to the upper ocean layer. The constant c_vol is the seawater's volumetric heat capacity, obtained by multiplying the specific heat capacity of seawater and the density of seawater. The thicknesses h_U and h_D of the two ocean layers are the same as for the carbon cycle. γ is thermal conductivity between the ocean layers and β is the climate feedback parameter related to the equilibrium climate sensitivity (in ^∘C),

$$\begin{array}{}\text{(31)}& \mathrm{ECS}={\displaystyle \frac{F_{\mathrm{2}X}}{\mathit{\beta }}},\end{array}$$

with F_2X the radiative forcing corresponding to a doubling of CO₂. The anthropogenic radiative forcing (measured in W m⁻²) in Eq. (29) responsible for the temperature anomalies consists of two terms,

$$\begin{array}{}\text{(32)}& \begin{array}{rl}F(M_{\mathrm{A}},I)& =F_{\mathrm{2}X}{\log }_{\mathrm{2}}\left({\displaystyle \frac{M_{\mathrm{A}}}{M_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)}}\right)\\ & -{\mathit{\alpha }}_{{\mathrm{SO}}_{\mathrm{2}}}\exp \left(-({\mathit{\beta }}_{{\mathrm{SO}}_{\mathrm{2}}}/I{)}^{{\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}}\right),\end{array}\end{array}$$

a first one corresponding to the standard greenhouse effect and a second one corresponding to solar radiation modification in the form of SO₂ injections. The solar radiation modification term comes from Niemeier and Timmreck (2015). The variable I corresponds to the sulfur injection rate and is in general time-dependent. ${\mathit{\alpha }}_{{\mathrm{SO}}_{\mathrm{2}}}$, ${\mathit{\beta }}_{{\mathrm{SO}}_{\mathrm{2}}}$ and ${\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}$ are the fitting parameters considered in Niemeier and Timmreck (2015).

2.3 Sea level rise

Four different components contribute to SURFER's estimation of sea level rise: ocean thermal expansion, glaciers and two ice sheets,

$$\begin{array}{}\text{(33)}& S_{\mathrm{tot}}=S_{\mathrm{th}}+S_{\mathrm{gl}}+S_{\mathrm{GIS}}+S_{\mathrm{AIS}},\end{array}$$

where S_tot is the total sea level rise and S_th, S_gl, S_GIS and S_AIS are the ocean thermal expansion, glacier, Greenland ice sheet and Antarctic ice sheet contributions, respectively. While ocean thermal expansion and glaciers are the first contributors to sea level rise on the short timescales of decades, on the longer timescales of centuries and millennia, the biggest contributions come from the ice sheets.

2.3.1 Ocean thermal expansion

The ocean thermal expansion parameterisation relies on the ocean layer sizes and temperature anomalies that are part of SURFER's climate sub-model. This contribution to sea level rise is then computed as

$$\begin{array}{}\text{(34)}& S_{\mathrm{th}}={\mathit{\alpha }}_{\mathrm{U}}{h}_{\mathrm{U}}\mathit{\delta }{T}_{\mathrm{U}}+{\mathit{\alpha }}_{\mathrm{D}}{h}_{\mathrm{D}}\mathit{\delta }{T}_{\mathrm{D}},\end{array}$$

where α_i is the thermal expansion coefficient corresponding to the i layer (in ^∘C⁻¹), h_i is the size of the i layer (in metres) and δT_i is the temperature anomaly (in ^∘C) with respect to pre-industrial times of the i layer. We consider that the expansion coefficients of the two layers have different values to capture the fact that surface waters have bigger thermal expansion coefficients than deeper denser waters as shown, for example, in Fig. 1c of Williams et al. (2012). As a simplification, we neglect the size change of the ocean layers h_(U,D) due to sea level rise and we also assume that the expansion coefficients are constant in time.

The sea level rise contribution from ocean thermal expansion comes from both ocean layers. In the timescales of decades and a couple of centuries, the deeper ocean layer does not contribute much due to its thermal inertia. As the deeper ocean warms up, in the timescale of thousands of years, it can become the main contributor to S_th. Figure 3a shows the sea level rise commitment from ocean expansion once thermal equilibrium has been achieved between the two ocean layers, that is, with δT_U=δT_D.

2.3.2 Glaciers

The sea level rise contribution from glaciers is modelled with an ordinary differential equation that relaxes the current sea level rise value due to glaciers, S_gl, to its expected equilibrium value for the current temperature S_gl eq(δT_U),

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{gl}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{gl}}}}\left(S_{\mathrm{gl}\mathrm{eq}}\right(\mathit{\delta }{T}_{\mathrm{U}})-S_{\mathrm{gl}}),\end{array}$$

where τ_gl is the relaxation timescale. The same form of equation was used by Mengel et al. (2016), although with a differently parameterised S_gl eq. Levermann et al. (2013) analysed the sea level commitments of different sea level rise components depending on the forcing temperature. They estimate the shape of such S_gl eq(δT_U) for all land glaciers excluding ice sheets (see Fig. 1b in Levermann et al., 2013) which we will approximate as

$$\begin{array}{}\text{(36)}& S_{\mathrm{gl}\mathrm{eq}}\left(\mathit{\delta }{T}_{\mathrm{U}}\right)=S_{\mathrm{gl}\mathrm{pot}}\tanh \left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{U}}}{\mathit{\zeta }}}\right),\end{array}$$

where S_gl pot is the potential sea level rise due to glaciers, which corresponds to all the ice volume in glaciers in units of sea level rise equivalent, and ζ is a sensitivity coefficient. Figure 3b shows the shape of S_gl eq for the suggested values of S_gl pot and ζ in Table 6.

Figure 3 Equilibrium sea level rise contribution from ocean thermal expansion and glaciers.

This way of modelling the glaciers has a couple of advantages to similar methods used in other simple models. First of all, the formulation is fully transparent and nothing more than Eqs. (35) and (36) are needed. As for the behaviour, it captures some expected physics:

1. No forcing corresponds to no sea level rise from glaciers.

2. For small enough forcings, different levels of forcing lead to different levels of sea level rise.

3. Different levels of forcing on the same initial state generate different rates of sea level rise.

4. There is a cap on the maximum amount of sea level rise that can come from glaciers.

2.3.3 Ice sheets

Multi-stability regions and tipping points have been identified both for the Greenland and Antarctic ice sheets (e.g. Lenton et al., 2008; Letreguilly et al., 1991; Pattyn, 2006; Ridley et al., 2010; Robinson et al., 2012; Gregory et al., 2020; Garbe et al., 2020). The proposed ice sheet model highlights those tipping points and is easy to adapt to both ice sheets such that it captures their dynamics. The state variable is the volume fraction of ice with respect to a reference state, which we set to be the ice sheet's pre-industrial state. The ice sheets' contributions to sea level rise with respect to pre-industrial times is computed as a function of the ice sheets' melted fractions of ice and their total sea level rise potential,

$$\begin{array}{}\text{(37)}& S_{\mathrm{GIS}}=S_{\mathrm{GIS}\mathrm{pot}}(\mathrm{1}-V_{\mathrm{GIS}}(t\left)\right),S_{\mathrm{AIS}}=S_{\mathrm{AIS}\mathrm{pot}}(\mathrm{1}-V_{\mathrm{AIS}}(t\left)\right),\end{array}$$

where S_GIS pot and S_AIS pot are the sea level rise potentials of Greenland's and Antarctic ice sheet, respectively, and V_GIS and V_AIS their volume fraction with respect to their pre-industrial volume.

To capture the dynamics of an ice sheet featuring multi-stability and tipping points we propose a non linear ordinary differential equation for the ice volume fraction,

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}=\mathit{\mu }(V,\mathit{\delta }{T}_{\mathrm{U}})\left(-V^{\mathrm{3}}+a_{\mathrm{2}}{V}^{\mathrm{2}}+a_{\mathrm{1}}V+c_{\mathrm{1}}\mathit{\delta }{T}_{\mathrm{U}}+c_{\mathrm{0}}\right),\end{array}$$

where the third-order polynomial of the volume fraction and the term proportional to a forcing temperature imply a double-fold bifurcation diagram for the steady states in terms of a constant forcing temperature, see Fig. 4. In contrast to the ocean part of the carbon cycle model and climate model presented before, the ice sheet model is not explicitly derived from physical processes; it is a generic dynamical system model based on the concept of a double-fold bifurcation to be calibrated on state-of-the-art ice sheet models' output. In that sense, it is an emulator. The different terms of the polynomial are not there because they represent specific physical processes but because as a whole they produce the desired steady-state structure.

Figure 4 Ice sheet steady states for model defined in Eqs. (38) and (39).

The constant parameters $(a_{\mathrm{2}},a_{\mathrm{1}},c_{\mathrm{1}},c_{\mathrm{0}})$ are given in terms of the bifurcation points $\left(\right(T_{+},V_{+}),(T_{-},V_{-}\left)\right)$ as

$$\begin{array}{}\text{(39a)}& {\displaystyle }{a}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{3}(V_{-}+V_{+})}{\mathrm{2}}},\text{(39b)}& {\displaystyle }{a}_{\mathrm{1}}& {\displaystyle }=-\mathrm{3}{V}_{-}{V}_{+},\text{(39c)}& {\displaystyle }{c}_{\mathrm{1}}& {\displaystyle }=-{\displaystyle \frac{(V_{+}-V_{-}{)}^{\mathrm{3}}}{\mathrm{2}(T_{+}-T_{-})}},\text{(39d)}& {\displaystyle }{c}_{\mathrm{0}}& {\displaystyle }=+{\displaystyle \frac{T_{+}{V}_{-}^{\mathrm{2}}(V_{-}-\mathrm{3}{V}_{+})-T_{-}{V}_{+}^{\mathrm{2}}(V_{+}-\mathrm{3}{V}_{-})}{\mathrm{2}(T_{-}-T_{+})}},\end{array}$$

which are the quantities which determine the steady-state structure of the system, see Fig. 4. Since we want to impose the additional constraint of there existing a steady-state ice volume fraction of 1 at temperature anomaly equal to 0, the number of free parameters is reduced by one by setting,

$$\begin{array}{}\text{(40)}& V_{-}={\displaystyle \frac{-\mathrm{2}+V_{+}\left(\mathrm{1}+G^{\mathrm{1}/\mathrm{3}}+G^{-\mathrm{1}/\mathrm{3}}\right)}{-\mathrm{1}+G^{\mathrm{1}/\mathrm{3}}+G^{-\mathrm{1}/\mathrm{3}}}},\end{array}$$

with

$$\begin{array}{}\text{(41)}& G=\left({\displaystyle \frac{T_{+}+T_{-}+\mathrm{2}\sqrt{T_{-}{T}_{+}}}{T_{+}-T_{-}}}\right).\end{array}$$

Instead of fixing the 4 parameters $\left(\right(T_{+},V_{+}),(T_{-},V_{-}\left)\right)$ independently, V₋ is given in terms of the other three such that the pre-industrial reference state condition $V_{\mathrm{eq}}(\mathit{\delta }{T}_{\mathrm{U}}=\mathrm{0})=\mathrm{1}$ is satisfied.

The evolution is further affected by the forcing temperature δT_U and inverse timescale,

$$\begin{array}{}\text{(42)}& \mathit{\mu }(V,\mathit{\delta }{T}_{\mathrm{U}})=\left\{\begin{array}{ll}\mathrm{1}/{\mathit{\tau }}_{+}& \text{if}H>\mathrm{0},\\ \mathrm{1}/{\mathit{\tau }}_{-}& \text{if}H<\mathrm{0}\text{and}V>\mathrm{0},\\ \mathrm{0}& \text{if}H<\mathrm{0}\text{and}V=\mathrm{0},\end{array}\right.\end{array}$$

with

$$\begin{array}{}\text{(43)}& H=\left(-V^{\mathrm{3}}+a_{\mathrm{2}}{V}^{\mathrm{2}}+a_{\mathrm{1}}V+c_{\mathrm{1}}\mathit{\delta }{T}_{\mathrm{U}}+c_{\mathrm{0}}\right).\end{array}$$

We write μ(V,δT_U) as a function of the state variables V and δT_U in such a way that it can take three different constant values. This is done to reflect the timescale asymmetry between the processes of melting and freezing ice and to ensure that the ice fraction remains bigger or equal to zero.

2.4 Calibration and initial conditions

In this section we give the values and units of the parameters used to run the model. We also provide an explanation of how we have fixed some of them and obtained the pre-industrial initial conditions that we use in Sect. 3.

Parameters and initial conditions for the carbon cycle model can be found in Tables 1 and 2, the parameters for the climate model are in Table 3 and the ones for sea level rise can be found in Tables 4, 5 and 6 for the ocean thermal expansion, glaciers and ice sheet (adapted to Greenland and Antarctica) models, respectively. All examples in this paper start from pre-industrial conditions, which for the climate and sea level rise models are trivial, $\mathit{\delta }{T}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)=\mathit{\delta }{T}_{\mathrm{D}}\left(t_{\mathrm{PI}}\right)=\mathrm{0}$, S_gl(t_PI)=0, $V_{\mathrm{GIS}}\left(t_{\mathrm{PI}}\right)=V_{\mathrm{AIS}}\left(t_{\mathrm{PI}}\right)=\mathrm{1}$.

2.4.1 Carbon cycle

The dissociation parameters K₁ and K₂, the CO₂ solubility K₀, and the alkalinity Q are assumed to be constant. A temperature (T)- and salinity (S)-dependent alternative is provided for K₁, K₂ and K₀ in Sarmiento and Gruber (2006, p. 325)⁵. The specific values that we will be using for the constants K₁, K₂ and K₀ have been obtained from the expressions in Sarmiento and Gruber (2006). We fix the parameter δ_DIC to 1.15 in accordance with data provided in Sarmiento and Gruber (2006, p. 320).

We fix the parameters K₁, K₂ and K₀ in the carbon cycle model by ensuring that initial conditions are consistent with pre-industrial data. The parameter $\mathit{\delta }=h_{\mathrm{D}}/h_{\mathrm{U}}$, also playing a role in the initial conditions, has been set to δ=20, which is a middle ground between the thermal box sizes of Gregory (2000); Held et al. (2010) and the ratio of sizes suggested in Bolin and Eriksson (1959); Glotter et al. (2014), the reason for the chosen value is provided below. Other parameters, like k_A→U, k_U→D and k_A→L, have been adjusted to match the dynamics of more complex carbon cycle models. We now proceed with a more in-depth explanation.

In pre-industrial times the atmospheric CO₂ concentration was 280ppm which corresponds to an atmospheric carbon mass,

$$\begin{array}{ll}{\displaystyle }{M}_{\mathrm{A}}\left(t_{\mathrm{PI}}\right)& {\displaystyle }=\text{moles of}{\mathrm{CO}}_{\mathrm{2}}\text{in atmosphere}\left(t_{\mathrm{PI}}\right)\times {\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}\\ {\displaystyle }& {\displaystyle }=\mathrm{280}\times {\mathrm{10}}^{-\mathrm{6}}{m}_{\mathrm{A}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}=\mathrm{580.3}\mathrm{PgC},\end{array}$$

where ${\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}$ is the carbon molar mass. We assume that in pre-industrial conditions the carbon cycle was in equilibrium. Additionally we impose that the total mass of carbon in the ocean was 38 000 PgC. Using the relations between ocean carbon masses, the corresponding DICs and the equilibrium equations we can write

$$\begin{array}{ll}{\displaystyle }{M}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)+M_{\mathrm{D}}\left(t_{\mathrm{PI}}\right)& {\displaystyle }={\mathrm{DIC}}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)m_{\mathrm{O}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{C}}{\stackrel{\mathrm{‾}}{m}}_{\mathrm{W}}{\displaystyle \frac{\mathrm{1}+{\mathit{\delta }}_{\mathrm{DIC}}\mathit{\delta }}{\mathrm{1}+\mathit{\delta }}}\\ \text{(44)}& {\displaystyle }& {\displaystyle }=\mathrm{38}\mathrm{000}\mathrm{PgC},\end{array}$$

where we have written the carbon and water molar masses as ${\stackrel{\mathrm{‾}}{m}}_{(\mathrm{C},\mathrm{W})}$. We see that δ=20 implies a reasonable value of DIC_U(t_PI)=1973.53 µmol kg⁻¹ by comparing to global average data in Fig. 8.1.2 of Sarmiento and Gruber (2006) and to global averages from ZECMIP, see pre-industrial range of values in Fig. 9c. This value of DIC is one of the ingredients used to fix the dissociation and solubility constants.

The pH of the upper ocean in pre-industrial times was around 8.2. Here we will fix it to 8.17 in accordance with historical CMIP6 results in Gutiérrez et al. (2021) which implies a hydrogen ion concentration of $\left[{\mathrm{H}}^{+}\right]\left(t_{\mathrm{PI}}\right)={\mathrm{10}}^{-\mathrm{8.17}}$ mol kg⁻¹. Q_U and DIC_U at pre-industrial time can be written as

$$\begin{array}{rl}& Q_{\mathrm{U}}=\left({\displaystyle \frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]\left(t_{\mathrm{PI}}\right)}}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{{\left(\left[{\mathrm{H}}^{+}\right]\left(t_{\mathrm{PI}}\right)\right)}^{\mathrm{2}}}}\right)K_{\mathrm{0}}{p}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}\left(t_{\mathrm{PI}}\right),\\ & {\mathrm{DIC}}_{\mathrm{U}}\left(t_{\mathrm{PI}}\right)=\left(\mathrm{1}+{\displaystyle \frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]\left(t_{\mathrm{PI}}\right)}}+{\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}}{{\left(\left[{\mathrm{H}}^{+}\right]\left(t_{\mathrm{PI}}\right)\right)}^{\mathrm{2}}}}\right)K_{\mathrm{0}}{p}_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}\left(t_{\mathrm{PI}}\right),\end{array}$$

where $p_{{\mathrm{CO}}_{\mathrm{2}}}^{\mathrm{A}}\left(t_{\mathrm{PI}}\right)=\mathrm{280}\times {\mathrm{10}}^{-\mathrm{6}}$ atm. We use these relations to fix the dissociation and solubility constants. First, we impose an alkalinity value compatible with observations Q_U=2200 µmol kg⁻¹. Second, we impose the already fixed value of DIC_U(t_PI)=1973.53 µmol kg⁻¹. Third, we use the temperature (T)- and salinity (S)-dependent expressions for the dissociation and solubility given in Sarmiento and Gruber (2006). Last, we solve the system of two equations for T and S numerically. Such a procedure yields “effective” T=294.7 K and S=32.49 ‰ a warmer and slightly less salty ocean than the global averages, but they determine dissociation constants which, in the end, yield realistic carbon masses, concentrations and alkalinity in the pre-industrial ocean.

We ignore the temperature dependence of the carbonate concentration corresponding to aragonite saturation and we fix it to

$$\begin{array}{}\text{(45)}& [{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}{]}_{\text{saturation ar}}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left(t_{\mathrm{PI}}\right)}{{\mathrm{\Omega }}_{\mathrm{ar}}\left(t_{\mathrm{PI}}\right)}},\end{array}$$

where [CO${}_{\mathrm{3}}^{-\mathrm{2}}$](t_PI) is obtained through Eqs. (5) and (12) and Ω_ar(t_PI)=3.44 from Rockström et al. (2009); Steffen et al. (2015).

The ZECMIP B1 and B3 experiments' outputs suggest a land carbon parameter β_L between 0.5 and 2.3. We settle for 1.7 but the choice is not critical. This value is closer to those of ESMs than to those of the EMICs as can be seen in Figs. 9f and 10f. The pre-industrial mass of land carbon is set to 2200 PgC for plotting purposes but this quantity does not affect the land carbon uptake.

Finally, for the inverse timescale k_A→U we take the value recommended by Glotter et al. (2014). k_U→D is fixed to obtain a timescale for the deep ocean dynamics of 1000 years $k_{\mathrm{U}\to \mathrm{D}}=\mathit{\delta }{\mathit{\delta }}_{\mathrm{DIC}}/\mathrm{1000}$. The inverse timescale k_A→L is fixed to match the output of ZECMIP B1 and B3 experiments.

Users of the model are invited to explore other possibilities of fixing parameters but in this paper we will restrict ourselves to using SURFER only with the supplied parameters and initial conditions, and we show that despite its simplicity, its predictions are in excellent agreement with more complex models.

Table 1 Parameter values for the carbon cycle model.

Table 2 Initial equilibrium pre-industrial conditions obtained for the parameters specified in Table 1.

2.4.2 Climate

In the climate model we need to fix the parameters h_U (or h_D, since the ratio δ has already been fixed), c_vol, F_2X, β and γ, and the parameters from the aerosol forcing ${\mathit{\alpha }}_{{\mathrm{SO}}_{\mathrm{2}}}$, ${\mathit{\beta }}_{{\mathrm{SO}}_{\mathrm{2}}}$ and ${\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}$. The values adopted in the following are listed on Table 3.

The sea water volumetric heat capacity c_vol is obtained by multiplying the specific heat capacity of seawater, which is taken to be 3850 J kg^−1 ∘C⁻¹ and the density of seawater, taken to be 1027 kg m⁻³. The extra radiative forcing due to a doubling concentration of CO₂, F_2X, has been chosen according to the IPCC 6th Assessment Report (Arias et al., 2021) and parameter h_U according to Gregory (2000). β and γ have been fixed to yield equilibrium climate sensitivity and transient climate response,

$$\begin{array}{}\text{(46)}& \mathrm{ECS}={\displaystyle \frac{F_{\mathrm{2}X}}{\mathit{\beta }}}=\mathrm{3.5}{}^{\circ }\mathrm{C},\mathrm{TCR}={\displaystyle \frac{F_{\mathrm{2}X}}{\mathit{\beta }+\mathit{\gamma }}}=\mathrm{2.0}{}^{\circ }\mathrm{C},\end{array}$$

compatible with results in the 6th IPCC Assessment Report (Arias et al., 2021), i.e. ECS very likely 2 to 5 ^∘C and TCR very likely 1.4 to 2.2 ^∘C. Aerosol forcing parameters come from the work of Niemeier and Timmreck (2015).

Table 3 Parameter values for the temperature module.

2.4.3 Ocean thermal expansion

The thermal expansion coefficients α_U and α_D in Eq. (34) have been first estimated by looking at the thermal expansion coefficient profile along the water column presented in Fig. 1c of Williams et al. (2012) and taking into account the sizes of h_U and h_D in SURFER. Then they have been slightly corrected to better match the long-term trends presented by Van Breedam et al. (2020). Figures 5 and 16b show the performance of SURFER's thermal expansion parameterisation against other models both on short and long timescales.

Figure 5 Sea level rise contribution from ocean thermal expansion for historical period and projections along RCP scenarios. SURFER's estimates correspond to the solid coloured lines and observations for the historical period from Frederikse et al. (2020a) are shown with solid black lines with the corresponding uncertainty range in grey shade. BRICK's mean and uncertainty range are shown with dotted lines and coloured shades. The dashed lines correspond to results of Van Breedam et al. (2020). Sea level rise is with respect to year 2000.

Table 4 Thermal expansion coefficients.

2.4.4 Glaciers

We have fixed the total sea level rise potential of glaciers (since pre-industrial times) to 0.5 m. This is an intermediate value between those reported by Levermann et al. (2013) and Farinotti et al. (2019). The sensitivity ζ has been fixed to 2 ^∘C to mimic the glaciers' commitment curve in Levermann et al. (2013). Finally, the timescale has been fixed to 200 years which corresponds to an intermediate value for the range found by Mengel et al. (2016). Figures 6 and 16c show the performance of SURFER's glaciers model against other models both on short and long timescales.

Figure 6 Sea level rise contribution from glaciers for historical period and projections along RCP scenarios. SURFER's estimates correspond to the solid coloured lines, observations for the historical period from Frederikse et al. (2020a) are shown with solid black lines with the corresponding uncertainty range in grey shade. BRICK's mean and uncertainty range are shown with dotted lines and coloured shades. The dashed lines correspond to results of Van Breedam et al. (2020). Sea level rise is with respect to year 2000.

Table 5 Glacier model parameters.

2.4.5 Ice sheets

We now proceed to fix the values of the parameters $(T_{+},T_{-},V_{+},{\mathit{\tau }}_{+},{\mathit{\tau }}_{-})$⁶ for adapting the ice sheet model to Greenland and Antarctica.

For Greenland, we have calibrated the bifurcation points $(T_{\pm },V_{+})$ by requiring the steady states of Eq. (38) to reproduce part of the steady-state structure found in Robinson et al. (2012). In Fig. 7a, we show the upper and lower steady-state branches found by Robinson et al. (2012) together with SURFER's double-fold steady-state structure. Then, the melting timescale τ₋ has been fixed to match the constant forcing transient results of Robinson et al. (2012)⁷, see Fig. 7b. Finally, not many references present accumulation (ice sheet growth) experiments which are needed to fix τ₊. For this reason we used different references to fix τ₊ (Letreguilly et al., 1991; Pattyn, 2006). The calibration of τ₊ was done in the way as that of τ₋, i.e. seeking to reproduce the dynamic accumulation experiments as in Letreguilly et al. (1991); Pattyn (2006).

Figure 7 SURFER's Greenland ice sheet against the work of Robinson et al. (2012).

For the case of Antarctica, fitting was more complicated than for Greenland. Contrary to the previous case, we are not aware of results from the same Antarctic ice sheet model setup for both steady-state and transient experiments. For example, for the PISM model, Garbe et al. (2020) provided the steady-state structure and Winkelmann et al. (2015) performed some transient experiments but since the climate models used for the forcing were not the same, the results are not completely compatible⁸. With this less than ideal situation, we have relied on the two most recent references, Garbe et al. (2020); Van Breedam et al. (2020), to fix SURFER's Antarctica parameters $(T_{\pm },V_{+},{\mathit{\tau }}_{-})$. The steady-state experiments of Garbe et al. (2020) have oriented us towards the range of acceptable parameters values while we have finally fixed those values by attempting to fit the dynamical experiments of Van Breedam et al. (2020). Again, in the absence of time-series results from accumulation experiments in the same references, τ₊ has been fixed by fitting to the results of Pattyn (2006) and Huybrechts (1993). Figure 16e contrasts SURFER's prediction for S_AIS when forced with extended RCP scenarios to the results of Van Breedam et al. (2020).

Finally, Greenland's and Antarctica's ice sheet sea level rise potential, linking the ice volume fraction to an eustatic sea level rise, come from Van Breedam et al. (2020) and Garbe et al. (2020), respectively.

As more models provide both steady-state structure and transient experiments, fitting can be improved by following the same strategy as used for Greenland. Additionally, it would be best if results from bigger models were presented separately for East and West Antarctica since the two components have different tipping points and evolve with different characteristic timescales. Then SURFER could treat East and West Antarctica as two separate ice sheets.

Again, SURFER is intended to be tuned on state-of-the art models. Its parameters can, therefore, be revised as new simulations become available. The best scenario would be to have a model inter-comparison project on the timescales of millennia for both ice sheets.

Table 6 Parameter values used here for Greenland and Antarctic ice sheets.

2.5 Numerics

The model has been implemented in The Julia programming language (Bezanson et al., 2017) in the Jupyter lab environment. The equations have been integrated using the package DifferentialEquations.jl with the integration method Rosenbrock23() with abstol=1e-12 and reltol=1e-3. The model runs extremely fast, each run of up to 10 000 years taking ≈0.003 s on a laptop with processor Intel^® Core™ i7-9850H CPU @ 2.60 GHz × 12.

3 Numerical results and comparisons

In this section we show the model's behaviour when forced by different CO₂ emission scenarios and SO₂ injections. Only the example in Sect. 3.5 considers SO₂ injections; all other examples are forced with CO₂ emissions alone. The different examples are meant to showcase the different parts of the model. Section 3.1 illustrates the long-term behaviour of the carbon cycle components, Sect. 3.2 the different ocean acidification metrics that are derivable with SURFER, Sect. 3.3 and Sect. 3.4 focus on the short- and long-term sea level rise behaviour respectively, and finally Sect. 3.5 deals with solar radiation modification. Whenever it has been possible to retrieve outputs of other models or historical data we have done so for making the comparisons easier. All examples start from pre-industrial conditions, more details on initial conditions can be found in Sect. 2.4.

3.1 ZECMIP B1 and B3 experiments

The Zero Emissions Commitment Model Intercomparison Project (ZECMIP) (Jones et al., 2019; MacDougall et al., 2020) was proposed to quantify the amount of unrealised temperature change that occurs after CO₂ emissions cease and to investigate the geophysical drivers behind this climate response. The ZECMIP B1 and B3 experiments consist in starting with pre-industrial conditions and forcing the system with the bell-shaped emission curves in Fig. 8, corresponding to 1000 PgC cumulated emissions for the B1 experiment and 2000 PgC for the B3 experiment.

Figure 8 CO₂ emission curves for ZECMIP experiments B1 and B3.

We consider the output of 6 EMICs and 2 ESMs that participated in the ZECMIP experiment. The EMICs are Bern3D-LPX-ECS3K, DCESS1.0, MESM, MIROC-lite-LCM, PLASIM-GENIE and UVicESCM2.10, and the ESMs are GFDL-ESM2M and NorESM2-LM. Figures 9 and 10 show the good agreement of SURFER's output to that of the ZECMIP experiments. SURFER's land model is closer to the ESMs than the EMICs although this can be changed by recalibrating β_L on the carbon cycle equations.

Figure 9 ZECMIP B1 experiment outputs of 6 EMICs, 2 ESMs and SURFER.

Figure 10 ZECMIP B3 experiment outputs of 6 EMICs, 2 ESMs and SURFER.

3.2 Ocean acidification under RCP or SSP forcing scenarios

In this section, we run SURFER forced with historic CO₂ emissions followed by CO₂ emissions associated to the representative concentration pathways (RCPs) and shared socio-economic pathways (SSPs) together with their extensions (Meinshausen et al., 2011, 2020)⁹, see Fig. 11.

Figure 11 CO₂ emission rate for the RCP and SSP scenarios with their extensions.

For the RCPs emission scenarios we have considered, as Van Breedam et al. (2020), that CO₂ emissions are zero from 2300 onwards. For the SSPs emission scenarios, as in Meinshausen et al. (2020), we consider that CO₂ emissions are zero from 2250.

Figure 12 shows the evolution of different variables when forced by these scenarios. As CO₂ is emitted into the atmosphere, the CO₂ concentration in the atmosphere rises and so does the temperature and the DIC in the ocean's upper layer. Together with the increase in DIC, the pH, carbonate concentration and aragonite saturation state decrease. For the aragonite saturation state (Ω_ar)_U, in Fig. 12e, we have shadowed the region beyond the planetary boundary for ocean acidification. SURFER predicts that all considered emission scenarios cross that planetary boundary. Scenarios RCP8.5, SSP5-8.5 and SSP3-7 reach values of aragonite saturation state smaller or close to 1; in those scenarios aragonite would dissolve in the upper ocean.

Figure 12 SURFER's output when forced by the extended RCPs and SSPs emission scenarios.

Figure 13 compares SURFER's pH prediction to that of CMIP5 and CMIP6 shown in Kwiatkowski et al. (2020). SURFER exhibits more acidification than CMIP5 and CMIP6 but it is in better agreement with CMIP6 runs.

Figure 13 Change in pH in the upper ocean layer under RCPs and SSPs emission scenarios.

3.3 Historical sea level rise and centennial projections

In this section, we again run SURFER forced with historic CO₂ emissions followed by CO₂ emissions associated to the shared socio-economic pathways (SSPs) together with their extensions (Meinshausen et al., 2020) up to year 2150.

We then contrast, in Fig. 14, SURFER's sea level rise predictions to the historic observations reported in Frederikse et al. (2020a) and the projections reported in Table 9.9 of IPCC AR6 WGI, see Fox-Kemper et al. (2021). The agreement of SURFER's predictions with both data sets is remarkable given the simplicity of the model.

Figure 14 Historic and shorter term sea level rise projections. Solid coloured lines show SURFER's total sea level rise estimates. Solid black lines with shaded envelope correspond to historic observations reported by Frederikse et al. (2020a), dashed coloured lines correspond to IPCC AR6 estimates and shaded regions to IPCC AR6 likely (medium confidence) ranges. The dash–dot line in SSP5-8.5 corresponds to the low confidence IPCC AR6 estimate. SURFER's and historical results are with respect to year 2004. IPCC results are relative to a baseline of 1995–2014.

3.4 Long-term (millennial) sea level rise projections

Here we show the model's long-term behaviour when forced by historic anthropogenic CO₂ emissions followed by the extended RCPs emission scenarios, see Moss et al. (2010); Meinshausen et al. (2011), which we continue to extend from the year 2300 onwards with zero emissions, as in Van Breedam et al. (2020), see Fig. 11.

Figure 15 shows the atmospheric CO₂ concentration and the global mean temperature anomaly for a period of 10 000 years and Fig. 16 contrasts SURFER's sea level rise predictions for each of the considered contributors with the corresponding results of Van Breedam et al. (2020).

Due the absence of calcium carbonate reactions in the carbon sub-model and the fact that SURFER does not include any feedback from the ice sheet melting on the temperature, quantities from the carbon cycle and climate sub-models reach equilibrium around year 6000, see Fig. 15. After that time, only the ice sheet components of the model continue to evolve. If carbonate compensation was included, a slow lowering of the atmospheric CO₂ concentration and temperature would be observed, together with its corresponding effect on the ice sheets¹⁰. Albedo feedbacks of the ice sheets on the climate sub-model would act in the opposite direction increasing the temperature. These effects, as mentioned in Sect. 4, are expected to be relatively small.

The agreement of SURFER's total sea level rise with the results of Van Breedam et al. (2020) is good for the higher emission scenarios and worse for the lower ones, see Fig. 16. The main source for the discrepancies comes, unsurprisingly, from the Greenland ice sheet contribution. The reason for this is that we calibrated Greenland's parameters by fitting to the results of Robinson et al. (2012) and that the tipping points do not coincide: Greenland's ice sheet in Robinson et al. (2012), and hence SURFER's, seems to have a higher temperature tipping point than the one of in Van Breedam et al. (2020). The agreement for Greenland's ice sheet with (Van Breedam et al., 2020), see Fig. 16d, is therefore not good for the lower emission scenarios because the corresponding forcing temperatures lie close to this threshold in SURFER but are clearly beyond it in Van Breedam et al. (2020).

SURFER's ocean thermal expansion parameterisation yields a S_th that agrees well with the results of Van Breedam et al. (2020), see Fig 16b, similarly for the Antarctic ice sheet, see Fig 16e. For the glaciers, we have different behaviour although the same order of magnitude which, as expected, is very sub-leading with respect to the other contributors at these timescales. While for SURFER the different RCP scenarios lead to different contributions from the glaciers (higher emission scenarios leading to more melting), in Van Breedam et al. (2020) all the considered ice in glaciers eventually melts at a rate that is scenario-dependent.

Figure 15 SURFER's long-term prediction when forced by extended RCP scenarios.

Figure 16 Long-term sea level rise from SURFER and Van Breedam et al. (2020) when forced by extended RCP scenarios. Solid lines correspond to SURFER and dashed to the results in Van Breedam et al. (2020). Colours indicate different RCP scenarios as in Fig. 15.

3.5 Solar radiation modification

As a last example, we perform an experiment done by the Geoengineering Model Inter-comparison Project (GeoMIP) and compare SURFER's output to the results of Visioni et al. (2021). We focus on the G6sulfur experiment introduced in Kravitz et al. (2015). In the G6sulfur experiment, stratospheric aerosols are injected into the model with the goal of reducing the magnitude of the net anthropogenic radiative forcing from a high forcing scenario (SSP5-8.5) to match that of a medium forcing scenario (SSP2-4.5).

We run SURFER forced by the CO₂ emissions corresponding to the SSP scenarios SSP5-8.5 and SSP2-4.5 in the absence of solar radiation modification. We compute the radiative forcing difference between the two scenarios, and we perform a third run in which SURFER is forced by the CO₂ emissions corresponding to scenario SSP5-8.5 and with sulfur injections exactly compensating the extra radiative forcing in SSP5-8.5 with respect to SSP2-4.5,

$$\begin{array}{}\text{(47)}& {\displaystyle }I\left(t\right)={\mathit{\beta }}_{{\mathrm{SO}}_{\mathrm{2}}}{\left(-\log \left({\displaystyle \frac{F^{\left(\text{SSP5-8.5}\right)}\left(t\right)-F^{\left(\text{SSP2-4.5}\right)}\left(t\right)}{{\mathit{\alpha }}_{{\mathrm{SO}}_{\mathrm{2}}}}}\right)\right)}^{\mathrm{1}/{\mathit{\gamma }}_{{\mathrm{SO}}_{\mathrm{2}}}},\end{array}$$

with F^(SSP5-8.5) and F^(SSP2-4.5) the time-varying radiative forcing obtained under the corresponding scenarios.

Figure 17 shows SURFER predictions for CO₂ concentration, radiative forcing, temperature and pH of upper ocean layer. This shows that while the radiative forcing, and hence the temperature, is lowered by the injection of aerosols, CO₂ concentration is high and therefore the ocean is dangerously acidic. Notice that this behaviour appears in SURFER by construction: in the version of the model presented here, there is no feedback from climate (temperature) into the carbon cycle. In Fig. 18, we compare the rate of sulfur injections needed by SURFER and by the ESMs that participated in the GeoMIP, see Visioni et al. (2021), to accomplish the G6sulfur experiment.

Figure 17 SURFER's output for the G6sulfur experiment in the black dotted line. Orange and blue lines are the SSP5-8.5 and SSP2-4.5 emission scenarios for comparison.

Figure 18 Comparison of SO₂ injection rate needed to do the G6sulfur experiment by SURFER and by ESMs participating in GeoMIP reported in Visioni et al. (2021).

4 Discussion

In the literature, other reduced complexity models for sea level rise can be found (e.g. Wong et al., 2017; Nauels et al., 2017; Palmer et al., 2020). Wong et al. (2017) introduce BRICK, a framework for modelling sea level rise. They argue that models for risk analysis should be accessible, transparent, flexible and efficient. In this paper, we have shown that SURFER complies with all of these criteria. Additionally, compared to BRICK, SURFER models ocean acidification and incorporates tipping points in the ice sheet components. Such phenomena are important when assessing policies (Lenton and Ciscar, 2013). Nauels et al. (2017) and Palmer et al. (2020) also provide an efficient reduced complexity model and a statistical emulator, respectively, that allow for sea level rise projections. These models, however, are already less transparent than SURFER and again, they do not incorporate tipping point dynamics. Finally, SURFER has been shown to reproduce results from EMICs, ESMs and 3D ice sheet models on timescales from decades to millennia. Is is unclear to us whether the BRICK, MAGICC and the model by Palmer et al. (2020) are applicable on such long timescales. Thus, we are convinced that SURFER is a valuable addition to the literature.

SURFER is simple, easy to understand and fast to run, but it misses some processes, carbon reservoirs and feedbacks. This gives room for extensions, and here we provide relevant information for users who would like to take some of the relevant processes into account, while explaining why we have not included them in the present model.

Permafrost holds approximately 1400 PgC (Canadell et al., 2021): that is about twice as much carbon as currently contained in the atmosphere. The thawing and subsequent release of part of that carbon into the atmosphere may therefore constitute a substantial positive feedback on anthropogenic emissions. Several studies, simulations and observations have been made to quantify the strength of these effects, (MacDougall and Knutti, 2016; Chadburn et al., 2017; Burke et al., 2017; Turetsky et al., 2020; Burke et al., 2020) but the spread in the estimates is considerable. MacDougall and Knutti (2016) made projections of the release of carbon from permafrost soils using a perturbed parameter ensemble with the UVic-ESCM EMIC. Among other things, they computed the sensitivities of cumulative emitted carbon per extra degree of warming. They showed the (transient) sensitivities changed with time, obtaining 24, 39 and 47 PgC ^∘C⁻¹ for the years 2100, 2200 and 2300 under all RCP scenarios. A similar computation was done by Burke et al. (2017) with the IMOGEN model coupled to two different land models. In that case, although lower values were obtained, the sensitivities remained time-dependent (∼10, 20, 30 PgC ^∘C⁻¹ for 2100, 2200 and 2300 under all RCP scenarios). In an observation-based study, Chadburn et al. (2017) estimated an equilibrium sensitivity of permafrost area loss per degree of future global warming of 4.0${}_{-\mathrm{1.1}}^{+\mathrm{1.0}}$ million km^2 ∘C⁻¹. The climatic consequences for CO₂ emissions from thawed permafrost were not explored. In the ZECMIP experiments (Jones et al., 2019; MacDougall et al., 2020), only two of the models (NorESM2-LM and UVicESCM2.10) had a permafrost module. The presence of the module in these models is visible on the absolute size of their land carbon reservoir compared to the other models (see Figs. 9 and 10), but not on the evolution of the land reservoir anomalies, suggesting that, in these models at least, the permafrost feedback does not dominate other carbon fluxes. Burke et al. (2020) acknowledges that modelling future permafrost thaw and its resulting CO₂ emissions remains a challenge to ESMs and EMICs, and points to deeper soils and the inclusion of a representation of abrupt thawing events as the main points to be improved. For this reason we left the addition of this reservoir and its possible feedbacks for the future.

The value of the equilibrium climate sensitivity used in SURFER is compatible with that obtained by ESMs that include the effect of ice–albedo feedback from the melting of sea ice. Consequently, the ice–albedo feedback coming from sea ice can be thought of as already included in SURFER even if sea ice is not explicitly part of the model. A representation of the ice–albedo feedback due to the melting of the ice sheets, however, is missing in SURFER. Wunderling et al. (2020) studied the impact of the melting of several cryosphere elements on the global mean temperature. They found that the full melting of Greenland's ice sheet would lead to a temperature increase of 0.13 ^∘C. Approximately 50 % of which corresponds to albedo, and the rest to changes in lapse rate, water vapour and clouds. For the West Antarctic ice sheet, an increase in temperature of 0.05 ^∘C was found, again with 50 % of it coming from albedo effects. The East Antarctic ice sheet was not included in the study. Adding these effects has been left for future work in which more data from geographically explicit models become available.

SURFER's carbon cycle does not take into account the soft tissue or the carbonate pump which act to create a dissolved inorganic carbon gradient in the water column with the deeper layer containing more DIC than the upper layer. We have introduced a constant parameter, δ_DIC, which allows for this difference in DIC between layers to be captured. The soft tissue and carbonate pumps might be affected by temperature and CO₂ concentration changes, leading to a time evolving δ_DIC, while in SURFER this parameter is constant.

The feedback of temperature on the physical component of the carbon cycle can be included by substituting the solubility constant K₀ and the dissociation constants K₁ and K₂ by their temperature-dependent expressions as done in Glotter et al. (2014). This feedback, which acts in the direction of increasing atmospheric CO₂ by reducing the ocean carbon uptake, is however known to be small (Glotter et al., 2014) due to two competing processes; while solubility decreases with temperature, the dissociation constants increase.

Calcium carbonate compensation, a slower process in which the ocean's pH is neutralised by the dissolution of CaCO₃ from sediments is not present in SURFER. Calcium carbonate compensation acts on timescales of 3–7 kyr (Archer et al., 2009) and provides an extra buffer for atmospheric CO₂. Since the impact from ocean acidification on marine ecosystems takes place at much shorter timescales, and since the highest and most harmful rates of sea level rise occur before the year 4000 for the RCP scenarios, it seems unnecessary for our purposes to add this effect (together with a sediments reservoir) into SURFER. Due to this choice, its predictions will start to deviate from expected on timescales of several thousands of years and ocean acidification will be slightly intensified with respect to models which do include such processes. We have also ignored even longer timescale effects related to the CO₂ chemical reaction with silicate rocks.

In SURFER, solar radiation modification acts by construction only on the globally averaged temperature. Cao and Jiang (2017) studied the feedbacks between solar radiation modification and the carbon cycle. They found that, under solar radiation modification, the land carbon reservoir becomes more efficient in absorbing CO₂. This reduces the total atmospheric CO₂ concentration and thus the total amount of SO₂ injections needed to reach a certain temperature goal. It also reduces very slightly ocean acidification. On the other hand, Tjiputra et al. (2016) reported an increased ocean carbon uptake in the presence of solar radiation modification and a negligible global mean change in land carbon uptake. While the results by Tjiputra et al. (2016) and Cao and Jiang (2017) show the importance of investigating the subject further, we decided to wait for a more coherent picture before including such feedbacks in SURFER. Such feedbacks could potentially be included by making the land carbon equation, and the solubility and dissociation constants, temperature-dependent. The fact that in SURFER solar radiation modification only affects the temperature may implicitly cause an under-appreciation of the risks associated with this technology, especially its direct impact on atmospheric chemistry, circulation, precipitation and its indirect impact on health, food security and ecological systems (Zarnetske et al., 2021). As already said in the introduction, this inclusion is motivated by the need to put in a single coherent framework the contrasting timescales of the (short) residence time of stratospheric aerosols in the atmosphere, the (long) residence time of carbon in the ocean–atmosphere system and the different timescales involved in the sea level response. We nevertheless urge potential users of SURFER interested in evaluating the cost and benefits and related dilemmas associated with solar radiation modification to adequately incorporate its potentially severe adverse side effects in their studies and conclusions. Solar radiation modification impacts cannot be reduced to temperature and any study considering this technology as an option should account for the multifaceted risks associated with it (Robock, 2016, 2020; Zarnetske et al., 2021). Not doing so will likely lead to bad and dangerous advice.

5 Conclusions

In this paper we have presented SURFER, a new model that is easy to understand and modify, fast to apply, and hence well suited to be used for policy assessment. SURFER emulates the results of EMICs and ESMs regarding CO₂ concentration, temperature anomalies, ocean acidification and sea level rise. Since it emulates well a range of aggregated quantities, it can be used in policy assessments that value policies according to a wider criteria than just global mean temperature. Furthermore, due to its lightness and ability to correctly represent millennial timescales, it is also well suited for long-sighted decision problems and for commitment and responsibility analyses in presence of uncertainty, among other kinds of policy assessments.

We have shown that SURFER's sea level rise module performs exceptionally well on the short timescales. Moreover, this module includes ice sheet tipping points, which are particularly important for good performance on the long timescales. Being parameterised on the ice sheet tipping points, the module can easily be updated to match latest research in this fast growing area. With such a flexible ice sheet module, one may envision representing uncertainty about tipping behaviour as parameter uncertainty and investigating the effect of this non-deterministic setting on policy assessments.

Finally, we have shown that SURFER works well under a variety of forcing scenarios. These scenarios do not only include different rates of positive CO₂ emissions and a range of total cumulative emissions, but also future technologies such as solar radiation modification and, for SSP1-2.6 and RCP2.6, carbon dioxide removal (in the form of atmospheric negative CO₂ emissions). As a consequence, SURFER is well suited for policy assessments that require considering a variety of forcing scenarios. This is the case for sequential decision problems and also for commitment assessments that capture uncertainty about available future options regarding earth management. This last application will be the main focus of a companion paper that relies on SURFER as it's main computational engine.