Articles | Volume 15, issue 21
Model description paper
09 Nov 2022
Model description paper |  | 09 Nov 2022

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

Marina Martínez Montero, Michel Crucifix, Victor Couplet, Nuria Brede, and Nicola Botta

We present SURFER, a novel reduced model for estimating the impact of CO2 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 CO2 emission and solar radiation modification policies, for supporting the computation of optimal (CO2 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 CO2 in the atmosphere–ocean system (Archer2005; Archer et al.2009). While reducing CO2 emissions is necessary, different geoengineering methods have been proposed to reduce some of the negative effects of climate change and complement the efforts of CO2 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 CO2 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 CO2 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 CO2 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; Shepherd2019) 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 CO2 concentrations in the atmosphere are not free from risks (see Gardiner2010; Moreno-Cruz and Keith2013; Robock2016, 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 CO2 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 CO2 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 CO2 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 CO2 emissions and solar radiation modification in the form of SO2 injections.

The CO2 emission rate, through the accumulation of CO2 in the atmosphere and due to the greenhouse effect, leads to long-lasting temperature increases, while the SO2 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 Keith2013; Helwegen et al.2019; Robock2020) 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 CO2 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 SO2 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 1Conceptual 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 reservoir1, 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:


where Mi is the mass of carbon in reservoir i, Fij 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).


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 Gruber2006 for a deeper treatment).

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


where H2CO3* represents a mix of aqueous carbon dioxide, CO2(aqueous), and carbonic acid, H2CO32. The distribution between the three carbon species, H2CO3*, HCO3- (bicarbonate), and CO32- (carbonate), is fast with respect to the ocean's circulation timescale, and hence equilibrium is assumed. The equilibrium distribution relations,

(5) K 1 = [ H + ] [ HCO 3 - ] [ H 2 CO 3 * ] , K 2 = [ H + ] [ CO 3 2 - ] [ HCO 3 - ] ,

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

(6) pH = - log 10 [ H + ] .

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


Figure 2Carbon 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 CO2 dissolves in the ocean, hydrogen ions are released and ocean acidifies. This in turn means that the proportion of carbonate decreases and that of H2CO3* increases. The alkalinity, however, defined as the excess of bases over acids,

(8)Q=[HCO3-]+2[CO32-]+[OH-]-[H+]+[B(OH)4-]+minor bases,

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 alkalinity3 by its dominant terms, that is by the carbonate alkalinity,

(9) Q [ HCO 3 - ] + 2 [ CO 3 2 - ] = K 1 [ H + ] + 2 K 1 K 2 [ H + ] 2 [ H 2 CO 3 * ] .

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

(10) F A U p CO 2 A - p CO 2 U ,

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

(11) F A U = - k A U M A + k U A M U ,

where pCO2U refers to the partial pressure of carbon in the form H2CO3*, and MU to the corresponding carbon mass. The parameters kij are the transport coefficients from reservoir i to reservoir j.

The equilibrium concentration of H2CO3* in the ocean, corresponding to an atmospheric CO2 partial pressure pCO2A can be determined through the CO2 solubility constant,

(12) K 0 = [ H 2 CO 3 * ] p CO 2 ,

where pCO2 is written without a reservoir index because, when in equilibrium, atmospheric and upper ocean have the same CO2 partial pressure.


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,

(13) F U D = k U D M U - k D U M D .

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 H2CO3*. Now we can write the carbon cycle equations as


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


These allow to relate transport coefficients,


where we can further write

MU=[H2CO3*]WUmCMA=moles ofCO2in atmosphere×mC,M(U,D)=DIC(U,D)W(U,D)mC,

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

(18) W U m O m W h U h U + h D , W D m O m W h D h U + h D ,

with mO the moles of water in the ocean, mW the molar mass of H2O, mC the molar mass of carbon and hU and hD the thicknesses of the ocean layers. Considering the equilibrium solubility relation (Eq. 12) and

pCO2A=1(atm)moles ofCO2in atmospheremoles in atmosphere,

we arrive at


where mA are the moles of air in the atmosphere and

(21) δ DIC = DIC D ( t PI ) DIC U ( t PI ) ,

where DICU(tPI) and DICD(tPI) 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 H2CO3* form, MU, in Eq. (14a) and (14b) as a function of the state variables. We begin by


Using the definitions of DIC and carbonate alkalinity (Q in Eq. 9), and the relation of DICi with carbon mass Mi, [H+]i can be solved for in terms of Mi:

(22) [ H + ] i = K 1 2 Q ̃ i ( ( M i - Q ̃ i ) 2 - 4 K 2 K 1 Q ̃ i ( Q ̃ i - 2 M i ) + ( M i - Q ̃ i ) ) ,

with i=U or D and where


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

(23) B ( M U ) M U M U = 1 2 - Q ̃ U 2 M U + K 1 ( M U - Q ̃ U ) 2 - 4 K 2 K 1 Q ̃ U ( Q ̃ U - 2 M U ) - 4 M U K 2 2 M U ( K 1 - 4 K 2 ) .


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:

(24) δ M L ( t ) = α L δ M A ( t ) ,

where δMi(t)=Mi(t)-Mi(tPI). 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,

(25) δ M L ( t ) δ M A ( t ) M A ( t ) M A ( t PI ) ,

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 δML,

(26) d δ M L d t = k A L β L M A ( t PI ) M A δ M A - δ M L .

With this equation δML relaxes to an equilibrium value proportional to the ratio δMAMA. This dependency can be interpreted as the result of two competing processes: CO2 fertilisation (δMA) and an enhanced bacterial respiration due to climate change (MA).

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


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 MU 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 (Nordhaus1992). 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 CO2 from the atmosphere it becomes a worse sink for future CO2 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 MU(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

(28) Ω ar [ CO 3 2 - ] [ CO 3 2 - ] saturation ar ,

where [CO32-] is the carbonate concentration in the ocean and [CO32-]saturation ar is the carbonate concentration at which aragonite is saturated4. 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(tPI)=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 δTU and δTD, respectively (measured in C),


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 cvol is the seawater's volumetric heat capacity, obtained by multiplying the specific heat capacity of seawater and the density of seawater. The thicknesses hU and hD 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),

(31) ECS = F 2 X β ,

with F2X the radiative forcing corresponding to a doubling of CO2. The anthropogenic radiative forcing (measured in W m−2) in Eq. (29) responsible for the temperature anomalies consists of two terms,

(32) F ( M A , I ) = F 2 X log 2 M A M A ( t PI ) - α SO 2 exp - ( β SO 2 / I ) γ SO 2 ,

a first one corresponding to the standard greenhouse effect and a second one corresponding to solar radiation modification in the form of SO2 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. αSO2, βSO2 and γSO2 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,

(33) S tot = S th + S gl + S GIS + S AIS ,

where Stot is the total sea level rise and Sth, Sgl, SGIS and SAIS 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

(34) S th = α U h U δ T U + α D h D δ T D ,

where αi is the thermal expansion coefficient corresponding to the i layer (in C−1), hi is the size of the i layer (in metres) and δTi 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 Sth. 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 δTU=δTD.

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, Sgl, to its expected equilibrium value for the current temperature Sgl eq(δTU),

(35) d S gl d t = 1 τ gl ( S gl eq ( δ T U ) - S gl ) ,

where τgl is the relaxation timescale. The same form of equation was used by Mengel et al. (2016), although with a differently parameterised Sgl 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 Sgl eq(δTU) for all land glaciers excluding ice sheets (see Fig. 1b in Levermann et al.2013) which we will approximate as

(36) S gl eq ( δ T U ) = S gl pot tanh δ T U ζ ,

where Sgl 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 Sgl eq for the suggested values of Sgl pot and ζ in Table 6.

Figure 3Equilibrium 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; Pattyn2006; 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,

(37) S GIS = S GIS pot ( 1 - V GIS ( t ) ) , S AIS = S AIS pot ( 1 - V AIS ( t ) ) ,

where SGIS pot and SAIS pot are the sea level rise potentials of Greenland's and Antarctic ice sheet, respectively, and VGIS and VAIS 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,

(38) d V d t = μ ( V , δ T U ) - V 3 + a 2 V 2 + a 1 V + c 1 δ T U + c 0 ,

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 4Ice sheet steady states for model defined in Eqs. (38) and (39).


The constant parameters (a2,a1,c1,c0) are given in terms of the bifurcation points ((T+,V+),(T-,V-)) as


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,

(40) V - = - 2 + V + 1 + G 1 / 3 + G - 1 / 3 - 1 + G 1 / 3 + G - 1 / 3 ,


(41) G = T + + T - + 2 T - T + T + - T - .

Instead of fixing the 4 parameters ((T+,V+),(T-,V-)) independently, V is given in terms of the other three such that the pre-industrial reference state condition Veq(δTU=0)=1 is satisfied.

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

(42) μ ( V , δ T U ) = 1 / τ + if H > 0 , 1 / τ - if H < 0 and V > 0 , 0 if H < 0 and V = 0 ,


(43) H = - V 3 + a 2 V 2 + a 1 V + c 1 δ T U + c 0 .

We write μ(V,δTU) as a function of the state variables V and δTU 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, δTU(tPI)=δTD(tPI)=0, Sgl(tPI)=0, VGIS(tPI)=VAIS(tPI)=1.

2.4.1 Carbon cycle

The dissociation parameters K1 and K2, the CO2 solubility K0, and the alkalinity Q are assumed to be constant. A temperature (T)- and salinity (S)-dependent alternative is provided for K1, K2 and K0 in Sarmiento and Gruber (2006, p. 325)5. The specific values that we will be using for the constants K1, K2 and K0 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 K1, K2 and K0 in the carbon cycle model by ensuring that initial conditions are consistent with pre-industrial data. The parameter δ=hD/hU, 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 kA→U, kU→D and kA→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 CO2 concentration was 280ppm which corresponds to an atmospheric carbon mass,

MA(tPI)=moles ofCO2in atmosphere(tPI)×mC=280×10-6mAmC=580.3PgC,

where mC 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


where we have written the carbon and water molar masses as m(C,W). We see that δ=20 implies a reasonable value of DICU(tPI)=1973.53µmol kg−1 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 [H+](tPI)=10-8.17 mol kg−1. QU and DICU at pre-industrial time can be written as


where pCO2A(tPI)=280×10-6 atm. We use these relations to fix the dissociation and solubility constants. First, we impose an alkalinity value compatible with observations QU=2200µmol kg−1. Second, we impose the already fixed value of DICU(tPI)=1973.53µmol kg−1. 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

(45) [ CO 3 2 - ] saturation ar = [ CO 3 2 - ] ( t PI ) Ω ar ( t PI ) ,

where [CO3-2](tPI) is obtained through Eqs. (5) and (12) and Ωar(tPI)=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 kA→U we take the value recommended by Glotter et al. (2014). kU→D is fixed to obtain a timescale for the deep ocean dynamics of 1000 years kUD=δδDIC/1000. The inverse timescale kA→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 1Parameter values for the carbon cycle model.

Download Print Version | Download XLSX

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

Download Print Version | Download XLSX

2.4.2 Climate

In the climate model we need to fix the parameters hU (or hD, since the ratio δ has already been fixed), cvol, F2X, β and γ, and the parameters from the aerosol forcing αSO2, βSO2 and γSO2. The values adopted in the following are listed on Table 3.

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

(46) ECS = F 2 X β = 3.5 C , TCR = F 2 X β + γ = 2.0 C ,

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 3Parameter values for the temperature module.

Download Print Version | Download XLSX

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 hU and hD 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 5Sea 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 4Thermal expansion coefficients.

Download Print Version | Download XLSX

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 6Sea 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 5Glacier model parameters.

Download Print Version | Download XLSX

2.4.5 Ice sheets

We now proceed to fix the values of the parameters (T+,T-,V+,τ+,τ-)6 for adapting the ice sheet model to Greenland and Antarctica.

For Greenland, we have calibrated the bifurcation points (T±,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)7, 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; Pattyn2006). 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 7SURFER'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 compatible8. 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±,V+,τ-). 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 SAIS 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 6Parameter values used here for Greenland and Antarctic ice sheets.

Download Print Version | Download XLSX

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 CO2 emission scenarios and SO2 injections. Only the example in Sect. 3.5 considers SO2 injections; all other examples are forced with CO2 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 CO2 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 8CO2 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 9ZECMIP B1 experiment outputs of 6 EMICs, 2 ESMs and SURFER.


Figure 10ZECMIP 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 CO2 emissions followed by CO2 emissions associated to the representative concentration pathways (RCPs) and shared socio-economic pathways (SSPs) together with their extensions (Meinshausen et al.2011, 2020)9, see Fig. 11.

Figure 11CO2 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 CO2 emissions are zero from 2300 onwards. For the SSPs emission scenarios, as in Meinshausen et al. (2020), we consider that CO2 emissions are zero from 2250.

Figure 12 shows the evolution of different variables when forced by these scenarios. As CO2 is emitted into the atmosphere, the CO2 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 12SURFER'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 13Change 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 CO2 emissions followed by CO2 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 14Historic 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 CO2 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 CO2 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 CO2 concentration and temperature would be observed, together with its corresponding effect on the ice sheets10. 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 Sth 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 15SURFER's long-term prediction when forced by extended RCP scenarios.


Figure 16Long-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 CO2 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 CO2 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,

(47) I ( t ) = β SO 2 - log F ( SSP5-8.5 ) ( t ) - F ( SSP2-4.5 ) ( t ) α SO 2 1 / γ SO 2 ,

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 CO2 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, CO2 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 17SURFER'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 18Comparison of SO2 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 Ciscar2013). 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 Knutti2016; 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−1 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−1 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-1.1+1.0 million km2C−1. The climatic consequences for CO2 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 CO2 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 CO2 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 K0 and the dissociation constants K1 and K2 by their temperature-dependent expressions as done in Glotter et al. (2014). This feedback, which acts in the direction of increasing atmospheric CO2 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 CaCO3 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 CO2. 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 CO2 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 CO2. This reduces the total atmospheric CO2 concentration and thus the total amount of SO2 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 (Robock2016, 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 CO2 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 CO2 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 CO2 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.

Code availability

The exact version of SURFER used to produce the results used in this paper is archived on Zenodo (, Martínez Montero et al.2022b) under MIT license, as is the input data to run the model and produce the plots for all the simulations presented in this paper.

Data availability

Data from other references used in this paper for comparison with SURFER's output are as follows:

Author contributions

MMM, MC, NiB, and NuB conceptualised the project. MMM, MC, and VC contributed to the methodology by developing the model. MC managed and supervised the project. MMM wrote the model code (software) and carried out model validation. MMM and NuB did the visualisations. MMM wrote the original draft. MC, NiB, NuB, and VC substantially reviewed and edited the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank the editors and reviewers, whose comments have lead to significant improvements of the original paper. The authors also thank Jonas van Breedam, Alex Robinson and Daniele Visioni for sharing their data. This project is TiPES contribution #146: this project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 820970. Michel Crucifix is funded as Research Director by the Belgian National Fund of Scientific Research. Victor Couplet is funded as Research Fellow by the Belgian National Fund of Scientific Research (F.S.R.-FNRS).

Financial support

This project has received funding from the European Union’s Horizon 2020 research and innovation programme (grant no. 820970). Michel Crucifix is funded as Research Director by the Belgian National Fund of Scientific Research. Victor Couplet is funded as Research Fellow by the Belgian National Fund of Scientific Research (F.S.R. – FNRS).

Review statement

This paper was edited by Ludovic Räss and reviewed by Daniele Visioni and two anonymous referees.


Archer, D.: Fate of fossil fuel CO2 in geologic time, J. Geophys. Res., 110, C09S05,, 2005. a

Archer, D., Eby, M., Brovkin, V., Ridgwell, A., Cao, L., Mikolajewicz, U., Caldeira, K., Matsumoto, K., Munhoven, G., Montenegro, A., and Tokos, K.: Atmospheric Lifetime of Fossil Fuel Carbon Dioxide, Annu. Rev. Earth Pl. Sc., 37, 117–134,, 2009. a, b, c

Arias, P.A., Bellouin, N., Coppola, E., Jones, R., Krinner, G., Marotzke, J., Naik, V., Palmer, M., Plattner, G.-K., Rogelj, J., Rojas, M., Sillmann, J., Storelvmo, T., Thorne, P., Trewin, B., Rao, K. A., Adhikary, B., Allan, R., Armour, K., Bala, G., Barimalala, R., Berger, S., Canadell, J., Cassou, C., Cherchi, A., Collins, W., Collins, W., Connors, S., Corti, S., Cruz, F., Dentener, F., Dereczynski, C., Luca, A. D., Niang, A. D., Doblas-Reyes, F., Dosio, A., Douville, H., Engelbrecht, F., Eyring, V., Fischer, E., Forster, P., Fox-Kemper, B., Fuglestvedt, J., Fyfe, J., Gillett, N., Goldfarb, L., Gorodetskaya, I., Gutierrez, J., Hamdi, R., Hawkins, E., Hewitt, H., Hope, P., Islam, A., Jones, C., Kaufman, D., Kopp, R., Kosaka, Y., Kossin, J., Krakovska, S., Lee, J.-Y., Li, J., Mauritsen, T., Maycock, T., Meinshausen, M., Min, S.-K., Monteiro, P., Ngo-Duc, T., Otto, F., Pinto, I., Pirani, A., Raghavan, K., Ranasinghe, R., Ruane, A., Ruiz, L., Sallée, J.-B., Samset, B., Sathyendranath, S., Seneviratne, S., Sörensson, A., Szopa, S., Takayabu, I., Tréguier, A.-M., van den Hurk, B., Vautard, R., von Schuckmann, K., Zaehle, S., Zhang, X., and Zickfeld, K.: 2021: Technical Summary. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, in press, 2021. a, b

Bakker, A. M. R., Wong, T. E., Ruckert, K. L., and Keller, K.: Sea-level projections representing the deeply uncertain contribution of the West Antarctic ice sheet, Scientific Reports, 7, 3880,, 2017. a

Bakker, A. M. R., Wong, T. E., Ruckert, K. L., and Keller, K.: BRICK, The Pennsylvania State University [data set],, last access: 27 October 2022. a

Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.: Julia: A Fresh Approach to Numerical Computing, SIAM Rev., 59, 65–98,, 2017. a

Bolin, B. and Eriksson, E.: Distribution of matter in the sea and atmosphere: Changes in the Carbon Dioxide Content of the Atmosphere and Sea due to Fossil Fuel Combustion, in: The atmosphere and the sea in motion, The Rockefeller Institute Press, 130–142, ISBN 978-0-87470-033-6, (last access: 27 October 2022), 1959. a, b, c, d

Botta, N., Jansson, P., and Ionescu, C.: The impact of uncertainty on optimal emission policies, Earth Syst. Dynam., 9, 525–542,, 2018. a, b

Botta, N., Brede, N., Crucifix, M., Ionescu, C., Jansson, P., Li, Z., Martínez-Montero, M., and Richter, T.: Responsibility Under Uncertainty: Which Climate Decisions Matter Most?, Environ. Model. Assess., accepted,, 2021. a, b, c

Burke, E. J., Ekici, A., Huang, Y., Chadburn, S. E., Huntingford, C., Ciais, P., Friedlingstein, P., Peng, S., and Krinner, G.: Quantifying uncertainties of permafrost carbon–climate feedbacks, Biogeosciences, 14, 3051–3066,, 2017. a, b

Burke, E. J., Zhang, Y., and Krinner, G.: Evaluating permafrost physics in the Coupled Model Intercomparison Project 6 (CMIP6) models and their sensitivity to climate change, The Cryosphere, 14, 3155–3174,, 2020. a, b

Canadell, J.G., Monteiro, P. M. S., Costa, M. H., Cotrim da Cunha, L., Cox, P. M., Eliseev, A. V., Henson, S., Ishii, M., Jaccard, S., Koven, C., Lohila, A., Patra, P. K., Piao, S., Rogelj, J., Syampungani, S., Zaehle, S., and Zickfeld, K.: 2021: Global Carbon and other Biogeochemical Cycles and Feedbacks. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, in press, 2021. a

Cao, L. and Jiang, J.: Simulated Effect of Carbon Cycle Feedback on Climate Response to Solar Geoengineering, Geophys. Res. Lett., 44, 12484–12491,, 2017. a, b

Carlino, A., Giuliani, M., Tavoni, M., and Castelletti, A.: Multi-objective optimal control of a simple stochastic climate-economy model, IFAC-PapersOnLine, 53, 16593–16598,, 2020. a, b

Chadburn, S. E., Burke, E. J., Cox, P. M., Friedlingstein, P., Hugelius, G., and Westermann, S.: An observation-based constraint on permafrost loss as a function of global warming, Nat. Clim. Change, 7, 340–344,, 2017. a, b

Clark, P. U., Shakun, J. D., Marcott, S. A., Mix, A. C., Eby, M., Kulp, S., Levermann, A., Milne, G. A., Pfister, P. L., Santer, B. D., Schrag, D. P., Solomon, S., Stocker, T. F., Strauss, B. H., Weaver, A. J., Winkelmann, R., Archer, D., Bard, E., Goldner, A., Lambeck, K., Pierrehumbert, R. T., and Plattner, G. K.: Consequences of twenty-first-century policy for multi-millennial climate and sea-level change, Nat. Clim. Change, 6, 360–369,, 2016. a, b

Farinotti, D., Huss, M., Fürst, J. J., Landmann, J., Machguth, H., Maussion, F., and Pandit, A.: A consensus estimate for the ice thickness distribution of all glaciers on Earth, Nat. Geosci., 12, 168–173,, 2019. a

Fox-Kemper, B., Hewitt, H., Xiao, C., Aðalgeirsdóttir, G., Drijfhout, S., Edwards, T., Golledge, N., Hemer, M., Kopp, R., Krinner, G., Mix, A., Notz, D., Nowicki, S., Nurhati, I., Ruiz, L., Sallée, J.-B., Slangen, A., and Yu, Y.: Ocean, Cryosphere and Sea Level Change. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1211–1362, (last access: 27 October 2022), 2021. a, b

Frederikse, T., Landerer, F., Caron, L., Adhikari, S., Parkes, D., Humphrey, V. W., Dangendorf, S., Hogarth, P., Zanna, L., Cheng, L., and Wu, Y. H.: The causes of sea-level rise since 1900, Nature, 584, 393–397,, 2020a. a, b, c, d, e

Frederikse, T., Landerer, F., Caron, L., Adhikari, S., Parkes, D., Humphrey, V. W., Dangendorf, S., Hogarth, P., Zanna, L., Cheng, L., and Wu, Y.-H.: data supplement of “The causes of sea-level rise since 1900”, Zenodo [data set],, 2020b. a

Garbe, J., Albrecht, T., Levermann, A., Donges, J. F., and Winkelmann, R.: The hysteresis of the Antarctic Ice Sheet, Nature, 585, 538–544,, 2020. a, b, c, d, e, f

Gardiner, S. M.: Is “Arming the Future” with Geoengineering Really the Lesser Evil?: Some Doubts about the Ethics of Intentionally Manipulating the Climate System, in: Climate Ethics: Essential Readings, Oxford University Press, ISBN 978-0-19-539962-2, 978-0-19-756284-0,, 2010. a

Glotter, M. J., Pierrehumbert, R. T., Elliott, J. W., Matteson, N. J., and Moyer, E. J.: A simple carbon cycle representation for economic and policy analyses, Climatic Change, 126, 319–335,, 2014. a, b, c, d, e, f, g, h, i, j, k

Gregory, J. M.: Vertical heat transports in the ocean and their effect on time-dependent climate change, Clim. Dynam., 16, 501–515,, 2000. a, b, c

Gregory, J. M., George, S. E., and Smith, R. S.: Large and irreversible future decline of the Greenland ice sheet, The Cryosphere, 14, 4299–4322,, 2020. a

Gutiérrez, J. M., Jones, R. G., Narisma, G. T., Alves, L. M., Amjad, M., Gorodetskaya, I. V., Grose, M., Klutse, N. A. B., Krakovska, S., Li, J., Martínez-Castro, D., Mearns, L. O., Mernild, S. H., Ngo-Duc, T., van den Hurk, B., , and Yoon, J.-H.: 2021: Atlas. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, in press, (last access: 27 October 2022), 2021. a

Heitzig, J., Kittel, T., Donges, J. F., and Molkenthin, N.: Topology of sustainable management of dynamical systems with desirable states: from defining planetary boundaries to safe operating spaces in the Earth system, Earth Syst. Dynam., 7, 21–50,, 2016. a

Held, I. M., Winton, M., Takahashi, K., Delworth, T., Zeng, F., and Vallis, G. K.: Probing the fast and slow components of global warming by returning abruptly to preindustrial forcing, J. Climate, 23, 2418–2427,, 2010. a, b

Helwegen, K. G., Wieners, C. E., Frank, J. E., and Dijkstra, H. A.: Complementing CO2 emission reduction by solar radiation management might strongly enhance future welfare, Earth Syst. Dynam., 10, 453–472,, 2019. a, b, c, d

Huybrechts, P.: Glaciological Modelling of the Late Cenozoic East Antarctic Ice Sheet: Stability or Dynamism?, Geogr. Ann. A, 75, 221–238,, 1993. a

Jones, C. D., Frölicher, T. L., Koven, C., MacDougall, A. H., Matthews, H. D., Zickfeld, K., Rogelj, J., Tokarska, K. B., Gillett, N. P., Ilyina, T., Meinshausen, M., Mengis, N., Séférian, R., Eby, M., and Burger, F. A.: The Zero Emissions Commitment Model Intercomparison Project (ZECMIP) contribution to C4MIP: quantifying committed climate changes following zero carbon emissions, Geosci. Model Dev., 12, 4375–4385,, 2019. a, b, c, d

Kravitz, B., Robock, A., Tilmes, S., Boucher, O., English, J. M., Irvine, P. J., Jones, A., Lawrence, M. G., MacCracken, M., Muri, H., Moore, J. C., Niemeier, U., Phipps, S. J., Sillmann, J., Storelvmo, T., Wang, H., and Watanabe, S.: The Geoengineering Model Intercomparison Project Phase 6 (GeoMIP6): simulation design and preliminary results, Geosci. Model Dev., 8, 3379–3392,, 2015. a

Kwiatkowski, L., Torres, O., Bopp, L., Aumont, O., Chamberlain, M., Christian, J. R., Dunne, J. P., Gehlen, M., Ilyina, T., John, J. G., Lenton, A., Li, H., Lovenduski, N. S., Orr, J. C., Palmieri, J., Santana-Falcón, Y., Schwinger, J., Séférian, R., Stock, C. A., Tagliabue, A., Takano, Y., Tjiputra, J., Toyama, K., Tsujino, H., Watanabe, M., Yamamoto, A., Yool, A., and Ziehn, T.: Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections, Biogeosciences, 17, 3439–3470,, 2020. a

Lawrence, M. G., Schäfer, S., Muri, H., Scott, V., Oschlies, A., Vaughan, N. E., Boucher, O., Schmidt, H., Haywood, J., and Scheffran, J.: Evaluating climate geoengineering proposals in the context of the Paris Agreement temperature goals, Nat. Commun., 9, 2041–1723,, 2018. a

Lenton, T. M. and Ciscar, J. C.: Integrating tipping points into climate impact assessments, Climatic Change, 117, 585–597,, 2013. a

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., and Schellnhuber, H. J.: Tipping elements in the Earth's climate system, P. Natl. Acad. Sci. USA, 105, 1786–1793,, 2008. a

Letreguilly, A., Huybrechts, P., and Reeh, N.: Steady-state characteristics of the Greenland ice sheet under different climates, J. Glaciol., 37, 149–157,, 1991. a, b, c

Levermann, A., Clark, P. U., Marzeion, B., Milne, G. A., Pollard, D., Radic, V., and Robinson, A.: The multimillennial sea-level commitment of global warming, P. Natl. Acad. Sci. USA, 110, 13745–13750,, 2013. a, b, c, d

MacDougall, A. and Eby, M.: ZECMIP, School of Earth and Ocean Sciences University of Victoria [data set],, last access: 27 October 2022. a

MacDougall, A. H. and Knutti, R.: Projecting the release of carbon from permafrost soils using a perturbed parameter ensemble modelling approach, Biogeosciences, 13, 2123–2136,, 2016. a, b

MacDougall, A. H., Frölicher, T. L., Jones, C. D., Rogelj, J., Matthews, H. D., Zickfeld, K., Arora, V. K., Barrett, N. J., Brovkin, V., Burger, F. A., Eby, M., Eliseev, A. V., Hajima, T., Holden, P. B., Jeltsch-Thömmes, A., Koven, C., Mengis, N., Menviel, L., Michou, M., Mokhov, I. I., Oka, A., Schwinger, J., Séférian, R., Shaffer, G., Sokolov, A., Tachiiri, K., Tjiputra , J., Wiltshire, A., and Ziehn, T.: Is there warming in the pipeline? A multi-model analysis of the Zero Emissions Commitment from CO2, Biogeosciences, 17, 2987–3016,, 2020. a, b, c, d

Martinez Montero, M., Crucifix, M., Botta, N., and Brede, N.: Commitment as Lost Opportunities, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5928,, 2022a. a

Martínez Montero, M., Crucifix, M., Couplet, V., Brede, N., and Botta, N.: SURFER v2.0.0, Zenodo [code],, 2022b. a

Meinshausen, M., Smith, S. J., Calvin, K., Daniel, J. S., Kainuma, M. L., Lamarque, J., Matsumoto, K., Montzka, S. A., Raper, S. C., Riahi, K., Thomson, A., Velders, G. J., and van Vuuren, D. P.: The RCP greenhouse gas concentrations and their extensions from 1765 to 2300, Climatic Change, 109, 213–241,, 2011. a, b

Meinshausen, M., Nicholls, Z. R. J., Lewis, J., Gidden, M. J., Vogel, E., Freund, M., Beyerle, U., Gessner, C., Nauels, A., Bauer, N., Canadell, J. G., Daniel, J. S., John, A., Krummel, P. B., Luderer, G., Meinshausen, N., Montzka, S. A., Rayner, P. J., Reimann, S., Smith, S. J., van den Berg, M., Velders, G. J. M., Vollmer, M. K., and Wang, R. H. J.: The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500, Geosci. Model Dev., 13, 3571–3605,, 2020. a, b, c, d

Mengel, M., Levermann, A., Frieler, K., Robinson, A., Marzeion, B., and Winkelmann, R.: Future sea level rise constrained by observations and long-term commitment, P. Natl. Acad. Sci. USA, 113, 2597–2602,, 2016. a, b

Mengel, M., Nauels, A., Rogelj, J., and Schleussner, C. F.: Committed sea-level rise under the Paris Agreement and the legacy of delayed mitigation action, Nat. Commun., 9, 601,, 2018. a

Moreno-Cruz, J. B. and Keith, D. W.: Climate policy under uncertainty: A case for solar geoengineering, Climatic Change, 121, 431–444,, 2013. a, b

Moss, R. H., Edmonds, J. A., Hibbard, K. A., Manning, M. R., Rose, S. K., Van Vuuren, D. P., Carter, T. R., Emori, S., Kainuma, M., Kram, T., Meehl, G. A., Mitchell, J. F., Nakicenovic, N., Riahi, K., Smith, S. J., Stouffer, R. J., Thomson, A. M., Weyant, J. P., and Wilbanks, T. J.: The next generation of scenarios for climate change research and assessment, Nature, 463, 747–756,, 2010. a

Nauels, A., Meinshausen, M., Mengel, M., Lorbacher, K., and Wigley, T. M. L.: Synthesizing long-term sea level rise projections – the MAGICC sea level model v2.0, Geosci. Model Dev., 10, 2495–2524,, 2017. a, b

Nauels, A., Gütschow, J., Mengel, M., Meinshausen, M., Clark, P. U., and Schleussner, C. F.: Attributing long-term sea-level rise to Paris Agreement emission pledges, P. Natl. Acad. Sci. USA, 116, 23487–23492,, 2019. a

Niemeier, U. and Timmreck, C.: What is the limit of climate engineering by stratospheric injection of SO2?, Atmos. Chem. Phys., 15, 9129–9141,, 2015. a, b, c

Nordhaus, W. D.: Cowles Foundation Discussion Paper 1009: The “DICE”' Model: Background and Structure of a Dynamic Integrated Climate-Economy Model of the Economics of Global Warming, (last access: 26 October 2022), 1992. a

Palmer, M. D., Gregory, J. M., Bagge, M., Calvert, D., Hagedoorn, J. M., Howard, T., Klemann, V., Lowe, J. A., Roberts, C. D., Slangen, A. B., and Spada, G.: Exploring the Drivers of Global and Local Sea‐Level Change Over the 21st Century and Beyond, Earth's Future, 8, 2328–4277,, 2020. a, b, c

Pattyn, F.: GRANTISM: An Excel™ model for Greenland and Antarctic ice-sheet response to climate changes, Computers and Geosciences, 32, 316–325,, 2006. a, b, c, d

Ridley, J., Gregory, J. M., Huybrechts, P., and Lowe, J.: Thresholds for irreversible decline of the Greenland ice sheet, Clim. Dynam., 35, 1065–1073,, 2010. a

Robinson, A., Calov, R., and Ganopolski, A.: Multistability and critical thresholds of the Greenland ice sheet, Nat. Clim. Change, 2, 429–432,, 2012. a, b, c, d, e, f, g, h

Robock, A.: Albedo enhancement by stratospheric sulfur injections: More research needed, Earth's Future, 4, 644–648,, 2016. a, b

Robock, A.: Benefits and risks of stratospheric solar radiation management for climate intervention (Geoengineering), The Bridge by the National Academy of Engineering, 50, 59–67, (last access: 27 October 2022), 2020. a, b, c

Rockström, J., Steffen, W., Noone, K., Persson, Å., Chapin, F. S., Lambin, E. F., Lenton, T. M., Scheffer, M., Folke, C., Schellnhuber, H. J., Nykvist, B., de Wit, C. A., Hughes, T., van der Leeuw, S., Rodhe, H., Sörlin, S., Snyder, P. K., Costanza, R., Svedin, U., Falkenmark, M., Karlberg, L., Corell, R. W., Fabry, V. J., Hansen, J., Walker, B., Liverman, D., Richardson, K., Crutzen, P., and Foley, J. A.: A safe operating space for humanity, Nature, 461, 472–475,, 2009. a, b, c, d

Sarmiento, J. and Gruber, N.: Ocean Biogeochemical Dynamics, Princeton University Press, (last access: 26 October 2022), 2006. a, b, c, d, e, f, g, h

Shepherd, T. G.: Storyline approach to the construction of regional climate change information, P. Roy. Soc. A-Math., Phy., 475, 20190013,, 2019. a

Shepherd, T. G., Boyd, E., Calel, R. A., Chapman, S. C., Dessai, S., Dima-West, I. M., Fowler, H. J., James, R., Maraun, D., Martius, O., Senior, C. A., Sobel, A. H., Stainforth, D. A., Tett, S. F., Trenberth, K. E., van den Hurk, B. J., Watkins, N. W., Wilby, R. L., and Zenghelis, D. A.: Storylines: an alternative approach to representing uncertainty in physical aspects of climate change, Climatic Change, 151, 555–571,, 2018. a

Steffen, W., Richardson, K., Rockstrom, J., Cornell, S. E., Fetzer, I., Bennett, E. M., Biggs, R., Carpenter, S. R., de Vries, W., de Wit, C. A., Folke, C., Gerten, D., Heinke, J., Mace, G. M., Persson, L. M., Ramanathan, V., Reyers, B., and Sorlin, S.: Planetary boundaries: Guiding human development on a changing planet, Science, 347, 1259855,, 2015. a, b, c

Tjiputra, J. F., Grini, A., and Lee, H.: Impact of idealized future stratospheric aerosol injection on the large-scale ocean and land carbon cycles, J. Geophys. Res.-Biogeo., 121, 2–27,, 2016. a, b

Turetsky, M. R., Abbott, B. W., Jones, M. C., Anthony, K. W., Olefeldt, D., Schuur, E. A., Grosse, G., Kuhry, P., Hugelius, G., Koven, C., Lawrence, D. M., Gibson, C., Sannel, A. B. K., and McGuire, A. D.: Carbon release through abrupt permafrost thaw, Nat. Geosci., 13, 138–143,, 2020. a

Van Breedam, J., Goelzer, H., and Huybrechts, P.: Semi-equilibrated global sea-level change projections for the next 10 000 years, Earth Syst. Dynam., 11, 953–976,, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t

Visioni, D., MacMartin, D. G., Kravitz, B., Boucher, O., Jones, A., Lurton, T., Martine, M., Mills, M. J., Nabat, P., Niemeier, U., Séférian, R., and Tilmes, S.: Identifying the sources of uncertainty in climate model simulations of solar radiation modification with the G6sulfur and G6solar Geoengineering Model Intercomparison Project (GeoMIP) simulations, Atmos. Chem. Phys., 21, 10039–10063,, 2021.  a, b, c, d, e

Williams, R. G., Goodwin, P., Ridgwell, A., and Woodworth, P. L.: How warming and steric sea level rise relate to cumulative carbon emissions, Geophys. Res. Lett., 39, 4–9,, 2012. a, b

Winkelmann, R., Levermann, A., Ridgwell, A., and Caldeira, K.: Combustion of available fossil fuel resources sufficient to eliminate the Antarctic Ice Sheet, Science Advances, 1, 1–6,, 2015. a, b

Wong, T. E., Bakker, A. M. R., Ruckert, K., Applegate, P., Slangen, A. B. A., and Keller, K.: BRICK v0.2, a simple, accessible, and transparent model framework for climate and regional sea-level projections, Geosci. Model Dev., 10, 2741–2760,, 2017. a, b

Wunderling, N., Willeit, M., Donges, J. F., and Winkelmann, R.: Global warming due to loss of large ice masses and Arctic summer sea ice, Nat. Commun., 11, 5177,, 2020. a

Zarnetske, P. L., Gurevitch, J., Franklin, J., Groffman, P. M., Harrison, C. S., Hellmann, J. J., Hoffman, F. M., Kothari, S., Robock, A., Tilmes, S., Visioni, D., Wu, J., Xia, L., and Yang, C. E.: Potential ecological impacts of climate intervention by reflecting sunlight to cool Earth, P. Natl. Acad. Sci. USA, 118, e1921854118,, 2021. a, b, c


In Bolin and Eriksson (1959), they include humus and vegetation reservoirs. We have followed a different approach for obtaining the equations and fluxes corresponding to the land reservoir.


It is common practice to consider these two species of carbon together into a single variable because they are difficult to distinguish from each other (Sarmiento and Gruber2006, chap. 8.2).


See e.g. chap. 8 Sect. 8.2 of Sarmiento and Gruber (2006) for more details on this.


[CO32-]saturation ar depends strongly on pressure and hence changes along the water column.


Similar expressions are given in Glotter et al. (2014).


The parameter V has been fixed by using Eq. (40).


Notice that δTU is the global mean temperature anomaly and that their plots are for regional summer temperature anomaly. On their supplementary information they explain how to relate the two.


More melting happens in the setup of Winkelmann et al. (2015) than it would be expected by the results in Garbe et al. (2020).


For the SSPs we were only able to find the emission data up to 2100. We modelled the extension with a linear function going to zero by 2250. This is what was done by Meinshausen et al. (2020) except for the lowest emission scenario SSP1-2.6.


We expect this effect to be small because most of the ice sheet melting occurs in the first thousands of years where carbonate compensation still plays a very subdominant role.

Short summary
We present SURFER, a lightweight model that links CO2 emissions and geoengineering to ocean acidification and sea level rise from glaciers, ocean thermal expansion and Greenland and Antarctic ice sheets. The ice sheet module adequately describes the tipping points of both Greenland and Antarctica. SURFER is understandable, fast, accurate up to several thousands of years, capable of emulating results obtained by state of the art models and well suited for policy analyses.