Pathfinder v1.0.1: a Bayesian-inferred simple carbon–climate model to explore climate change scenarios

Abstract

The Pathfinder model was developed to fill a perceived gap within the range of existing simple climate models. Pathfinder is a compilation of existing formulations describing the climate and carbon cycle systems, chosen for their balance between mathematical simplicity and physical accuracy. The resulting model is simple enough to be used with Bayesian inference algorithms for calibration, which enables assimilation of the latest data from complex Earth system models and the IPCC sixth assessment report, as well as a yearly update based on observations of global temperature and atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. The model's simplicity also enables coupling with integrated assessment models and their optimization algorithms or running the model in a backward temperature-driven fashion. In spite of this simplicity, the model accurately reproduces behaviours and results from complex models – including several uncertainty ranges – when run following standardized diagnostic experiments. Pathfinder is an open-source model, and this is its first comprehensive description.

1 Introduction

Simple climate models (SCMs) typically simulate global mean temperature change caused by either atmospheric concentration changes or anthropogenic emissions of ${\mathrm{CO}}_{\mathrm{2}}$ and other climatically active species. They are most often composed of ad hoc parametric laws that emulate the behaviour of more complex Earth system models (ESMs). The emulation allows for simulating large ensembles of experiments that would be too costly to compute with ESMs. However, the SCM denomination refers to a fairly broad range of models whose complexity can go from a couple of boxes that only emulate one part of the climate system (e.g. a global temperature impulse response function; ) to hundreds of state variables representing the different cycles of greenhouse gases and their effect on climate change (e.g. the compact Earth system model OSCAR; ). Simpler models are easier and faster to solve, but they may not be adequate for all usages. Therefore, finding the “simplest but not simpler” model depends on a study's precise goals.

In our recent research, we have perceived a deficiency within the existing offer of SCMs, in spite of their large and growing number . We have therefore developed the Pathfinder model to fill this gap: it is a parsimonious ${\mathrm{CO}}_{\mathrm{2}}$-only model that carefully balances simplicity and accuracy of representation of physical processes. Pathfinder was designed to fulfil three key requirements: (1) the capacity to be calibrated using Bayesian inference, (2) the capacity to be coupled with integrated assessment models (IAMs), and (3) the capacity to explore a very large number of climate scenarios to narrow down those compatible with limiting climate impacts. The latter motivated the model's name.

While these three requirements clearly call for the simplest model possible, as they all need a fast solving model, they also imply a certain degree of complexity. The Bayesian calibration requires an explicit representation of the processes (i.e. the variables) that are used to constrain the model. Coupling with IAMs requires accurately embedding the latest advances of climate sciences to be policy relevant . Exploring future climate impacts requires the flexibility to link additional (and potentially regional) impact variables to the core carbon–climate equations.

The Pathfinder model is essentially an integration of existing formulations, adapted to our modelling framework and goals. It is calibrated on Earth system models that contributed to the Coupled Model Intercomparison Project phase 6 (CMIP6), on additional data from the sixth assessment report of the IPCC (AR6), and on observations of global Earth properties up to the year 2021. The calibration philosophy of Pathfinder is to use complex models as prior information and only real-world observations and assessments combining many lines of evidence as constraints.

Compared to other SCMs , Pathfinder is much simpler than models like MAGICC , OSCAR , or even HECTOR . It is comparable in complexity to FaIR or BernSCM , although it is closer to the latter as it trades off an explicit representation of non-${\mathrm{CO}}_{\mathrm{2}}$ species for one of the carbon cycle's main components. This choice was made to help calibration, keep the model invertible, and make the model compatible with IAMs such as DICE . While most SCMs are calibrated using procedures that resemble Bayesian inference , Pathfinder relies on an established algorithm whose implementation is fully tractable and that allows for an annual update as observations of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and global temperature become available.

Here, we present the first public release of Pathfinder and its source code. We first provide a detailed description of the model's equations. We then describe the Bayesian setup used for calibration, the sources of prior information for it, and the resulting posterior configuration. We end with a validation of the model using standard diagnostic simulations and quantitative metrics for the climate system and carbon cycle.

2 Equations

An overview of Pathfinder is presented in Fig. . The model is composed of a climate module, of three separate modules for the carbon cycle (ocean, land without land use and land permafrost), and of two additional modules describing global impacts: sea level rise (SLR) and surface ocean acidification. We do not emulate cycles of other non-CO${}_{\mathrm{2}}$ gases. Mathematically, the model is driven by prescribing time series of any combination of two of four variables: global mean surface temperature (GMST) anomaly ($T$), global atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration ($C$), global non-${\mathrm{CO}}_{\mathrm{2}}$ effective radiative forcing ($R_x$), and global anthropogenic emissions of ${\mathrm{CO}}_{\mathrm{2}}$ ($E_{{\mathrm{CO}}_{\mathrm{2}}}$). The model can therefore be run in the traditional emission-driven and concentration-driven modes but also in a temperature-driven mode (in terms of code, implemented as separate versions of the model). This is notably important for the calibration, during which it is driven by observations of GMST and atmospheric ${\mathrm{CO}}_{\mathrm{2}}$.

Figure 1

Pathfinder in a nutshell. Green blocks represent the carbon cycle, and red blocks represent the climate response. Blue blocks with dotted arrows are impacts that can be derived with the model. Grey blocks are variables that are directly related to anthropogenic activity. Possible inputs of the model are distinguishable through the bold contours of the blocks. In this scheme, arrows correspond to a forward mode where inputs would be $E_{{\mathrm{CO}}_{\mathrm{2}}}$ and $R_x$.

[Figure omitted. See PDF]

The following presents all equations of the models. Variables are noted using Roman letters and compiled in Tables  and . With a few exceptions, parameters are noted using Greek letters and are summarized in Tables  and . The model has 21 state variables that follow first-order differential equations in time. The time variable is denoted as $t$ and kept implicit unless required.

2.1 Climate

The GMST change ($T$) induced by effective radiative forcing (ERF; $R$) is represented using a widely used two-box energy balance model with deep-ocean heat uptake efficacy . The first box represents the Earth surface's temperature (including atmosphere, land, and surface ocean), and the other one is the deep ocean's temperature ($T_{\mathrm{d}}$). Their time-differential equations are

1$${\mathrm{\Theta }}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=R-{\displaystyle \frac{\mathit{\varphi }\ln \left(\mathrm{2}\right)}{T_{\mathrm{2}\times }}}T-{\mathit{ϵ}}_{\mathrm{heat}}\mathit{\theta }(T-T_{\mathrm{d}}),$$ and 2$${\mathrm{\Theta }}_{\mathrm{d}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}}=\mathit{\theta }(T-T_{\mathrm{d}}),$$ where $\mathit{\varphi }$ is the radiative parameter of ${\mathrm{CO}}_{\mathrm{2}}$, $T_{\mathrm{2}\times }$ is the equilibrium climate sensitivity (ECS) at ${\mathrm{CO}}_{\mathrm{2}}$ doubling, ${\mathrm{\Theta }}_{\mathrm{s}}$ is the heat capacity of the surface, ${\mathrm{\Theta }}_{\mathrm{d}}$ is the heat capacity of the deep ocean, $\mathit{\theta }$ is the heat exchange coefficient, and ${\mathit{ϵ}}_{\mathrm{heat}}$ is the deep-ocean heat uptake efficacy.

The global ERF is simply the sum of the ${\mathrm{CO}}_{\mathrm{2}}$ contribution ($R_{{\mathrm{CO}}_{\mathrm{2}}}$) from its change in atmospheric concentration ($C$), expressed using the IPCC AR5 formula , and that of non-${\mathrm{CO}}_{\mathrm{2}}$ climate forcers ($R_x$): 3$$R=R_{{\mathrm{CO}}_{\mathrm{2}}}+R_x,$$ with 4$$R_{{\mathrm{CO}}_{\mathrm{2}}}=\mathit{\varphi }\ln \left({\displaystyle \frac{C}{C_{\mathrm{pi}}}}\right),$$ where $C_{\mathrm{pi}}$ is the preindustrial atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration.

The above energy balance model naturally provides the ocean heat content (OHC; $U_{\mathrm{ohc}}$) as 5$$U_{\mathrm{ohc}}={\mathit{\alpha }}_{\mathrm{ohc}}({\mathrm{\Theta }}_{\mathrm{s}}T+{\mathrm{\Theta }}_{\mathrm{d}}{T}_{\mathrm{d}}),$$ and the ocean heat uptake (OHU) as 6$${\displaystyle \frac{\mathrm{d}{U}_{\mathrm{ohc}}}{\mathrm{d}t}}={\mathit{\alpha }}_{\mathrm{ohc}}\left({\mathrm{\Theta }}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}+{\mathrm{\Theta }}_{\mathrm{d}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{d}}}{\mathrm{d}t}}\right),$$ where ${\mathit{\alpha }}_{\mathrm{ohc}}$ is the fraction of energy used to warm the ocean (i.e. excluding the energy needed to heat up the atmosphere and land and to melt ice).

2.2 Sea level rise

Global SLR has been implemented in Pathfinder as a variable of interest to model climate change impacts. In this version, it is firstly a proof of concept, modelled in a simple yet sensible manner. The total sea level rise ($H_{\mathrm{tot}}$) is the sum of contributions from thermal expansion ($H_{\mathrm{thx}}$), the Greenland ice sheet (GIS; $H_{\mathrm{gis}}$), the Antarctic ice sheet (AIS; $H_{\mathrm{ais}}$), and glaciers ($H_{\mathrm{gla}}$):

7$$H_{\mathrm{tot}}=H_{\mathrm{thx}}+H_{\mathrm{gis}}+H_{\mathrm{ais}}+H_{\mathrm{gla}}.$$

The thermal expansion contribution scales linearly with the OHC : 8$$H_{\mathrm{thx}}={\mathrm{\Lambda }}_{\mathrm{thx}}{U}_{\mathrm{ohc}},$$ where ${\mathrm{\Lambda }}_{\mathrm{thx}}$ is the scaling factor of the thermosteric contribution to SLR. Note, however, that the thermal capacity of the climate module does not match that of the real-world ocean , and so this equation cannot describe equilibrium SLR over millennial timescales.

To model contributions from ice sheets and glaciers, we followed the general approach of . The SLR caused by GIS follows a first-order differential equation with its specific timescale, and the equilibrium SLR from GIS is assumed to be a cubic function of GMST: 9$${\displaystyle \frac{\mathrm{d}{H}_{\mathrm{gis}}}{\mathrm{d}t}}={\mathit{\lambda }}_{\mathrm{gis}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{gis}}}}\left({\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}T+{\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{3}}{T}^{\mathrm{3}}-H_{\mathrm{gis}}\right),$$ where ${\mathit{\lambda }}_{\mathrm{gis}}$ is an offset parameter introduced because GIS was not in a steady state at the end of the preindustrial era, ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}$ is the linear term of equilibrium of GIS SLR, ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{3}}$ is the cubic term of equilibrium of GIS SLR, and ${\mathit{\tau }}_{\mathrm{gis}}$ is the timescale of the GIS contribution. The motivation for replacing the quadratic term of with a cubic one is the oddness of the cubic function that leads to negative (and not positive) SLR for negative $T$ (which happens during the earlier years of the calibration run).

The contribution from glaciers is also a first-order differential equation with an equilibrium inspired by . We expanded it with a cubic term to account for the fact that we aggregate all glaciers together and allow more skewness in the curve describing the equilibrium SLR as a function of $T$. In addition, we added an exponential sensitivity to speed up the convergence to equilibrium under a warmer climate: 10$$\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{gla}}}{\mathrm{d}t}}& ={\mathit{\lambda }}_{\mathrm{gla}}+{\displaystyle \frac{\exp \left({\mathit{\gamma }}_{\mathrm{gla}}T\right)}{{\mathit{\tau }}_{\mathrm{gla}}}}\\ & \left({\mathrm{\Lambda }}_{\mathrm{gla}}\left(\mathrm{1}-\exp \left(-{\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{1}}T-{\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{3}}{T}^{\mathrm{3}}\right)\right)-H_{\mathrm{gla}}\right),\end{array}$$ where ${\mathit{\lambda }}_{\mathrm{gla}}$ is an offset parameter accounting for the lack of initial steady state, ${\mathrm{\Lambda }}_{\mathrm{gla}}$ is the SLR potential if all glaciers melted, ${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{1}}$ is the linear sensitivity of glaciers' equilibrium to climate change, ${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{3}}$ is the cubic sensitivity of glaciers' equilibrium to climate change, ${\mathit{\tau }}_{\mathrm{gla}}$ is the timescale of the glaciers contribution, and ${\mathit{\gamma }}_{\mathrm{gla}}$ is the sensitivity of glaciers' timescale to climate change.

Following , the contribution from AIS is further divided into two terms, one for surface mass balance (SMB; $H_{\mathrm{ais},\mathrm{smb}}$) and one for solid ice discharge (SID; $H_{\mathrm{ais},\mathrm{sid}}$), so that $H_{\mathrm{ais}}=H_{\mathrm{ais},\mathrm{smb}}+H_{\mathrm{ais},\mathrm{sid}}$. It is expected that precipitation will increase over Antarctica under higher GMST, leading to an increase in SMB and to a negative sea level rise contribution modelled as 11$${\displaystyle \frac{\mathrm{d}{H}_{\mathrm{ais},\mathrm{smb}}}{\mathrm{d}t}}=-{\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}T,$$ where ${\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}$ is the AIS SMB sensitivity to climate change (expressed in sea level equivalent). At the same time, increasing surface ocean temperatures will cause more SID through basal melting, which we model using a first-order differential equation assumed to be independent of the SMB effect, and with a term that speeds up the effect the more SID happened: 12$${\displaystyle \frac{\mathrm{d}{H}_{\mathrm{ais},\mathrm{sid}}}{\mathrm{d}t}}={\mathit{\lambda }}_{\mathrm{ais}}+{\displaystyle \frac{\mathrm{1}+{\mathit{\alpha }}_{\mathrm{ais}}{H}_{\mathrm{ais},\mathrm{sid}}}{{\mathit{\tau }}_{\mathrm{ais}}}}\left({\mathrm{\Lambda }}_{\mathrm{ais}}T-H_{\mathrm{ais},\mathrm{sid}}\right),$$ where ${\mathit{\lambda }}_{\mathrm{ais}}$ is an offset parameter accounting for the lack of initial steady state, ${\mathrm{\Lambda }}_{\mathrm{ais}}$ is the SLR equilibrium of AIS SID, ${\mathit{\tau }}_{\mathrm{ais}}$ is the timescale of the AIS SID contribution, and ${\mathit{\alpha }}_{\mathrm{ais}}$ is the sensitivity of the timescale to past SID. In the model's code, however, we directly solve for the total AIS contribution as 13$$\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{ais}}}{\mathrm{d}t}}& =-{\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}T+{\mathit{\lambda }}_{\mathrm{ais}}+{\displaystyle \frac{\mathrm{1}+{\mathit{\alpha }}_{\mathrm{ais}}(H_{\mathrm{ais}}-H_{\mathrm{ais},\mathrm{smb}})}{{\mathit{\tau }}_{\mathrm{ais}}}}\\ & \left({\mathrm{\Lambda }}_{\mathrm{ais}}T-(H_{\mathrm{ais}}-H_{\mathrm{ais},\mathrm{smb}})\right).\end{array}$$

2.3 Ocean carbon

To calculate the ocean carbon sink, we use the classic mixed-layer impulse response function model from , updated to the equivalent box model formulation of and extended in places to introduce parameter adjustments for calibration. In the model, the mixed layer is split into five boxes (subscript $j$), as represented in Fig. , so that the total carbon in the mixed-layer pool ($C_{\mathrm{o}}$) is

14$$C_{\mathrm{o}}=\sum _j{C}_{\mathrm{o},j}.$$ This total carbon mass is converted into a molar concentration of dissolved inorganic carbon (DIC; $c_{\mathrm{dic}}$) as follows: 15$$c_{\mathrm{dic}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{dic}}}{{\mathit{\beta }}_{\mathrm{dic}}}}{C}_{\mathrm{o}},$$ where ${\mathit{\alpha }}_{\mathrm{dic}}$ is a fixed conversion factor and ${\mathit{\beta }}_{\mathrm{dic}}$ is a scaling factor for the conversion. The latter can be seen as a factor multiplying the mixed-layer depth: it is 1 if the depth is unchanged from the original model.

Figure 2

The ocean sink model in Pathfinder follows the structure of the mixed-layer pulse response function introduced by . The mixed layer is represented through five sub-pools that each have a different timescale for transport to the deep-ocean carbon pool.

[Figure omitted. See PDF]

The non-linear carbonate chemistry in the mixed layer is emulated in two steps. First, the model's original polynomial function is used to determine the partial pressure of ${\mathrm{CO}}_{\mathrm{2}}$ affected by changes in DIC only ($p_{\mathrm{dic}}$): 16$${\displaystyle }\begin{array}{rl}{p}_{\mathrm{dic}}=& (\mathrm{1.5568}-\mathrm{0.013993}{T}_{\mathrm{o}})c_{\mathrm{dic}}\\ +& (\mathrm{7.4706}-\mathrm{0.20207}{T}_{\mathrm{o}}){\mathrm{10}}^{-\mathrm{3}}{c_{\mathrm{dic}}}^{\mathrm{2}}\\ -& (\mathrm{1.2748}-\mathrm{0.12015}{T}_{\mathrm{o}}){\mathrm{10}}^{-\mathrm{5}}{c_{\mathrm{dic}}}^{\mathrm{3}}\\ +& (\mathrm{2.4491}-\mathrm{0.12639}{T}_{\mathrm{o}}){\mathrm{10}}^{-\mathrm{7}}{c_{\mathrm{dic}}}^{\mathrm{4}}\\ -& (\mathrm{1.5768}-\mathrm{0.15326}{T}_{\mathrm{o}}){\mathrm{10}}^{-\mathrm{10}}{c_{\mathrm{dic}}}^{\mathrm{5}},\end{array}$$ where $T_{\mathrm{o}}$ is the preindustrial surface ocean temperature. Second, the actual partial pressure of ${\mathrm{CO}}_{\mathrm{2}}$ ($p_{{\mathrm{CO}}_{\mathrm{2}}}$) is calculated using an exponential climate sensitivity : 17$$p_{{\mathrm{CO}}_{\mathrm{2}}}=(p_{\mathrm{dic}}+C_{\mathrm{pi}})\exp \left({\mathit{\gamma }}_{\mathrm{dic}}T\right),$$ where ${\mathit{\gamma }}_{\mathrm{dic}}$ is the sensitivity of $p_{{\mathrm{CO}}_{\mathrm{2}}}$ to climate change.

The flux of carbon between the atmosphere and the ocean ($F_{\mathrm{ocean}}$, defined positively if it is a carbon sink) is caused by the difference in partial pressure of ${\mathrm{CO}}_{\mathrm{2}}$ in the atmosphere and at the oceanic surface, following an exchange rate that varies linearly with GMST, that is here used as a proxy for wind changes: 18$$F_{\mathrm{ocean}}={\mathit{\nu }}_{\mathrm{gx}}(\mathrm{1}+{\mathit{\gamma }}_{\mathrm{gx}}T)(C-p_{{\mathrm{CO}}_{\mathrm{2}}}),$$ where ${\mathit{\nu }}_{\mathrm{gx}}$ is the preindustrial gas exchange rate and ${\mathit{\gamma }}_{\mathrm{gx}}$ is its sensitivity to climate change.

This flux of carbon entering the ocean is split between the mixed-layer carbon sub-pools, and this added carbon is subsequently transported towards the deep ocean at a rate specific to each sub-pool. This leads to the following differential equations: 19$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{o},j}}{\mathrm{d}t}}=-{\displaystyle \frac{C_{\mathrm{o},j}}{{\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{o}}}{\mathit{\tau }}_{\mathrm{o},j}}}+{\mathit{\alpha }}_{\mathrm{o},j}{F}_{\mathrm{ocean}},\forall j,$$ where ${\mathit{\alpha }}_{\mathrm{o},j}$ are the sub-pools' splitting shares (with ${\sum }_j{\mathit{\alpha }}_{\mathrm{o},j}=\mathrm{1}$), ${\mathit{\tau }}_{\mathrm{o},j}$ are the sub-pools' timescales for transport to the deep ocean, and ${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{o}}}$ is a scaling factor applied to all sub-pools. Finally, the deep-ocean carbon pool ($C_{\mathrm{d}}$) is obtained through mass balance: 20$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{d}}}{\mathrm{d}t}}=\sum _j{\displaystyle \frac{C_{\mathrm{o},j}}{{\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{o}}}{\mathit{\tau }}_{\mathrm{o},j}}}.$$

2.4 Ocean acidification

While in the real world, ocean acidification is directly related to the carbonate chemistry and the ocean uptake of anthropogenic carbon, we do not have a simple formulation at our disposal that could link it to our ocean carbon cycle module. We therefore use a readily available emulation of the surface ocean acidification ($\mathrm{pH}$) that links it directly to the atmospheric concentration of ${\mathrm{CO}}_{\mathrm{2}}$ with the following polynomial approximation:

21$$\begin{array}{rl}\mathrm{pH}& ={\mathit{\kappa }}_{\mathrm{pH}}\left(\mathrm{8.5541}-\mathrm{0.00173}C+\mathrm{1.3264}{\mathrm{10}}^{-\mathrm{6}}{C}^{\mathrm{2}}\right.\\ & \left.-\mathrm{4.4943}{\mathrm{10}}^{-\mathrm{10}}{C}^{\mathrm{3}}\right),\end{array}$$ where ${\mathit{\kappa }}_{\mathrm{pH}}$ is a scaling factor (that defaults to 1). We note that this approach is reasonable for the surface ocean, as it quickly equilibrates with the atmosphere (but it would not work for the deep ocean).

2.5 Land carbon

The land carbon module of Pathfinder is a simplified version of the one in OSCAR . It is shrunk down to four global carbon pools: vegetation, litter, and active and passive soil (see Fig. ). All terrestrial biomes are lumped together, and there is therefore no accounting for the impact of land use change on the land carbon cycle in this version of Pathfinder. This is an extreme assumption – although very common in SCMs – motivated by simplicity, and it implies that ${\mathrm{CO}}_{\mathrm{2}}$ emissions from fossil fuel burning and land use change are assumed to behave in the exact same way, in spite of their not doing so in reality .

Figure 3

The land sink model in Pathfinder is derived from OSCAR and represents the biosphere as a set of four carbon pools: vegetation, litter, and active and passive soil. These pools exchange carbon through fluxes whose direction is given by the arrows.

[Figure omitted. See PDF]

The vegetation carbon pool ($C_{\mathrm{v}}$) results from the balance between net primary productivity (NPP; $F_{\mathrm{npp}}$), emission from wildfires ($E_{\mathrm{fire}}$), emission from harvest and grazing ($E_{\mathrm{harv}}$), and loss of carbon from biomass mortality ($F_{\mathrm{mort}}$): 22$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{v}}}{\mathrm{d}t}}=F_{\mathrm{npp}}-E_{\mathrm{fire}}-E_{\mathrm{harv}}-F_{\mathrm{mort}}.$$

NPP is expressed as its own preindustrial value multiplied by a function of ${\mathrm{CO}}_{\mathrm{2}}$ and of GMST ($r_{\mathrm{npp}}$). This function thus embeds the so-called ${\mathrm{CO}}_{\mathrm{2}}$ fertilization effect, whereby NPP increases with atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, described using a generalized logarithmic functional form: 23$$F_{\mathrm{npp}}=F_{\mathrm{npp}\mathrm{0}}{r}_{\mathrm{npp}},$$ with 24$$r_{\mathrm{npp}}=\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{npp}}}{{\mathit{\alpha }}_{\mathrm{npp}}}}\left(\mathrm{1}-{\left({\displaystyle \frac{C}{C_{\mathrm{pi}}}}\right)}^{-{\mathit{\alpha }}_{\mathrm{npp}}}\right)\right)(\mathrm{1}+{\mathit{\gamma }}_{\mathrm{npp}}T),$$ where $F_{\mathrm{npp}\mathrm{0}}$ is the preindustrial NPP, ${\mathit{\beta }}_{\mathrm{npp}}$ is the ${\mathrm{CO}}_{\mathrm{2}}$ fertilization sensitivity, ${\mathit{\alpha }}_{\mathrm{npp}}$ is the ${\mathrm{CO}}_{\mathrm{2}}$ fertilization shape parameter for saturation, and ${\mathit{\gamma }}_{\mathrm{npp}}$ is the sensitivity of NPP to climate change (that can be positive or negative). The generalized logarithmic functional form implies that $r_{\mathrm{npp}}\to (\mathrm{1}+{\mathit{\beta }}_{\mathrm{npp}}\ln (C/C_{\mathrm{pi}}\left)\right)(\mathrm{1}+{\mathit{\gamma }}_{\mathrm{npp}}T)$ as ${\mathit{\alpha }}_{\mathrm{npp}}\to {\mathrm{0}}^{+}$.

Harvesting and mortality fluxes are taken proportional to the carbon pool itself even though in reality the mortality fluxes are climate dependent. For simplicity we assume a constant mortality following the equations in OSCAR : 25$$E_{\mathrm{harv}}={\mathit{\nu }}_{\mathrm{harv}}{C}_{\mathrm{v}},$$ and 26$$F_{\mathrm{mort}}={\mathit{\nu }}_{\mathrm{mort}}{C}_{\mathrm{v}},$$ where ${\mathit{\nu }}_{\mathrm{harv}}$ is the harvesting or grazing rate and ${\mathit{\nu }}_{\mathrm{mort}}$ is the mortality rate.

Wildfires emissions are also assumed to be proportional to the vegetation carbon pool, but with an additional linear dependency of the emission rate on ${\mathrm{CO}}_{\mathrm{2}}$ (as a proxy of changes in leaf area index and evapotranspiration) and GMST ($r_{\mathrm{fire}}$): 27$$E_{\mathrm{fire}}={\mathit{\nu }}_{\mathrm{fire}}{r}_{\mathrm{fire}}{C}_{\mathrm{v}},$$ with 28$$r_{\mathrm{fire}}=\left(\mathrm{1}+{\mathit{\beta }}_{\mathrm{fire}}\left({\displaystyle \frac{C}{C_{\mathrm{pi}}}}-\mathrm{1}\right)\right)(\mathrm{1}+{\mathit{\gamma }}_{\mathrm{fire}}T),$$ where ${\mathit{\nu }}_{\mathrm{fire}}$ is the wildfire rate, ${\mathit{\beta }}_{\mathrm{fire}}$ is the sensitivity of wildfires to ${\mathrm{CO}}_{\mathrm{2}}$, and ${\mathit{\gamma }}_{\mathrm{fire}}$ is their sensitivity to climate change.

Soil carbon is divided into three pools. The litter carbon pool ($C_{\mathrm{s}\mathrm{1}}$) receives the mortality flux as sole input, it emits part of its carbon through heterotrophic respiration ($E_{\mathrm{rh}\mathrm{1}}$), and it transfers another part to the next pool through stabilization ($F_{\mathrm{stab}}$): 29$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{s}\mathrm{1}}}{\mathrm{d}t}}=F_{\mathrm{mort}}-F_{\mathrm{stab}}-E_{\mathrm{rh}\mathrm{1}}.$$ Similarly, the active soil carbon pool ($C_{\mathrm{s}\mathrm{2}}$) receives the stabilization flux, is respired ($E_{\mathrm{rh}\mathrm{2}}$), and transfers carbon to the last pool through passivization ($F_{\mathrm{pass}}$): 30$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{s}\mathrm{2}}}{\mathrm{d}t}}=F_{\mathrm{stab}}-F_{\mathrm{pass}}-E_{\mathrm{rh}\mathrm{2}}.$$ The passive carbon pool ($C_{\mathrm{s}\mathrm{3}}$) receives this final input flux and is respired ($E_{\mathrm{rh}\mathrm{3}}$): 31$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{s}\mathrm{3}}}{\mathrm{d}t}}=F_{\mathrm{pass}}-E_{\mathrm{rh}\mathrm{3}}.$$ Although information pertaining to this fourth pool is not commonly provided by ESMs, it was introduced in Pathfinder to adjust the complex models' turnover time of soil carbon to better match isotopic data . For completeness, we note that the total heterotrophic respiration is $E_{\mathrm{rh}}=E_{\mathrm{rh}\mathrm{1}}+E_{\mathrm{rh}\mathrm{2}}+E_{\mathrm{rh}\mathrm{3}}$, and the total soil carbon pool is $C_{\mathrm{s}}=C_{\mathrm{s}\mathrm{1}}+C_{\mathrm{s}\mathrm{2}}+C_{\mathrm{s}\mathrm{3}}$.

All soil-originating fluxes are taken proportional to their pool of origin and multiplied by a function ($r_{\mathrm{rh}}$) explained hereafter. For the litter pool, this gives 32$$E_{\mathrm{rh}\mathrm{1}}={\mathit{\nu }}_{\mathrm{rh}\mathrm{1}}{r}_{\mathrm{rh}}{C}_{\mathrm{s}\mathrm{1}},$$ and 33$$F_{\mathrm{stab}}={\mathit{\nu }}_{\mathrm{stab}}{r}_{\mathrm{rh}}{C}_{\mathrm{s}\mathrm{1}},$$ where ${\mathit{\nu }}_{\mathrm{rh}\mathrm{1}}$ is the litter respiration rate and ${\mathit{\nu }}_{\mathrm{stab}}$ is the stabilization rate. For the active soil pool, we have 34$$E_{\mathrm{rh}\mathrm{2}}={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}-{\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}{\mathit{\alpha }}_{\mathrm{pass}}}{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{pass}}}}{r}_{\mathrm{rh}}{C}_{\mathrm{s}\mathrm{2}},$$ and 35$$F_{\mathrm{pass}}={\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{pass}}}{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{pass}}}}{r}_{\mathrm{rh}}{C}_{\mathrm{s}\mathrm{2}},$$ and for the passive soil pool: 36$$E_{\mathrm{rh}\mathrm{3}}={\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}{r}_{\mathrm{rh}}{C}_{\mathrm{s}\mathrm{3}},$$ where ${\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}$ is the soil respiration rate (averaged over active and passive pools), ${\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}$ is the passive soil respiration rate, and ${\mathit{\alpha }}_{\mathrm{pass}}$ is the fraction of passive carbon (over active plus passive soil carbon). This slightly convoluted formulation is motivated by the lack of information regarding the active/passive split in ESMs, which we alleviate using additional data during calibration.

In addition, the function $r_{\mathrm{rh}}$, describing the dependency of respiration (and related fluxes) on temperature and on the availability of fresh organic matter to be decomposed, is defined as follows: 37$$\begin{array}{rl}& r_{\mathrm{rh}}=\\ & \underset{\left(\mathrm{1}+{\mathit{\beta }}_{\mathrm{rh}}\left(\frac{C_{\mathrm{s}\mathrm{1}}}{C_{\mathrm{s}}}\frac{C_{\mathrm{s}}\left(t_{\mathrm{0}}\right)}{C_{\mathrm{s}\mathrm{1}}\left(t_{\mathrm{0}}\right)}-\mathrm{1}\right)\right)}{\underbrace{\left(\mathrm{1}+{\mathit{\beta }}_{\mathrm{rh}}\left({\displaystyle \frac{C_{\mathrm{s}\mathrm{1}}}{C_{\mathrm{s}\mathrm{1}}+C_{\mathrm{s}\mathrm{2}}+C_{\mathrm{s}\mathrm{3}}}}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{stab}}}{{\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}}}\right)-\mathrm{1}\right)\right)}}\exp \left({\mathit{\gamma }}_{\mathrm{rh}}T\right),\end{array}$$ where ${\mathit{\beta }}_{\mathrm{rh}}$ is the sensitivity of the respiration to fresh organic matter availability (expressed here as the relative change in the $C_{\mathrm{s}\mathrm{1}}/C_{\mathrm{s}}$ ratio with regard to preindustrial times), and ${\mathit{\gamma }}_{\mathrm{rh}}$ is its sensitivity to climate change (equivalent to a “$Q_{\mathrm{10}}$” formulation with $Q_{\mathrm{10}}=\exp \left(\mathrm{10}{\mathit{\gamma }}_{\mathrm{rh}}\right)$).

Finally, the net carbon flux from the atmosphere to the land ($F_{\mathrm{land}}$, defined positively if it is a carbon sink) is obtained as the net budget of all pools combined: 38$$F_{\mathrm{land}}=F_{\mathrm{npp}}-E_{\mathrm{fire}}-E_{\mathrm{harv}}-E_{\mathrm{rh}},$$ and this system of equations leads to the following preindustrial steady state: 39$$\left\{\begin{array}{l}{C}_{\mathrm{v}}\left(t_{\mathrm{0}}\right)=\frac{F_{\mathrm{npp}\mathrm{0}}}{{\mathit{\nu }}_{\mathrm{fire}}+{\mathit{\nu }}_{\mathrm{harv}}+{\mathit{\nu }}_{\mathrm{mort}}}\\ C_{\mathrm{s}\mathrm{1}}\left(t_{\mathrm{0}}\right)=C_{\mathrm{v}}\left(t_{\mathrm{0}}\right)\frac{{\mathit{\nu }}_{\mathrm{mort}}}{{\mathit{\nu }}_{\mathrm{rh}\mathrm{1}}+{\mathit{\nu }}_{\mathrm{stab}}}\\ C_{\mathrm{s}\mathrm{2}}\left(t_{\mathrm{0}}\right)=C_{\mathrm{s}\mathrm{1}}\left(t_{\mathrm{0}}\right)\frac{{\mathit{\nu }}_{\mathrm{stab}}}{{\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}}(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{pass}})\\ C_{\mathrm{s}\mathrm{3}}\left(t_{\mathrm{0}}\right)=C_{\mathrm{s}\mathrm{1}}\left(t_{\mathrm{0}}\right)\frac{{\mathit{\nu }}_{\mathrm{stab}}}{{\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}}{\mathit{\alpha }}_{\mathrm{pass}}.\end{array}\right.$$

2.6 Permafrost carbon

As the land carbon cycle described in the previous section does not account for permafrost carbon, we implemented this feedback using the emulator developed by but aggregated into a unique global region. Figure  gives a representation of the permafrost module as described in the following. The emulation starts with a theoretical thawed fraction ($\overline{a}$) that represents the fraction of thawed carbon under steady-state for a certain level of local warming. It is formulated with a sigmoid function (that equals 0 at preindustrial and 1 under very high GMST):

40$$\begin{array}{rl}& \overline{a}=\\ & -a_{\min }+{\displaystyle \frac{(\mathrm{1}+a_{\min })}{{\left(\mathrm{1}+\left({\left(\mathrm{1}+\frac{\mathrm{1}}{a_{\min }}\right)}^{{\mathit{\kappa }}_a}-\mathrm{1}\right)\exp (-{\mathit{\gamma }}_a{\mathit{\kappa }}_a{\mathit{\alpha }}_{\mathrm{lst}}T)\right)}^{\frac{\mathrm{1}}{{\mathit{\kappa }}_a}}}},\end{array}$$ where $-a_{\min }$ is the minimum thawed fraction (corresponding to 100 % frozen soil carbon), ${\mathit{\kappa }}_a$ is a shape parameter determining the asymmetry of the function, ${\mathit{\gamma }}_a$ is the sensitivity of the theoretical thawed fraction to local climate change, and ${\mathit{\alpha }}_{\mathrm{lst}}$ is the proportionality factor between local and global climate change.

Figure 4

The permafrost carbon model in Pathfinder is taken from . The frozen pool dynamic lags behind a theoretical value that is determined by the temperature anomaly. Thawed carbon is then split between three pools that are emitted to the atmosphere at different rates.

[Figure omitted. See PDF]

The actual thawed fraction ($a$) then moves towards its theoretical value at a speed that depends on whether it is thawing (i.e. $a<\overline{a}$) or freezing (i.e. $a>\overline{a}$). This is written as a non-linear differential equation: 41$$\begin{array}{rl}{\displaystyle \frac{\mathrm{d}a}{\mathrm{d}t}}& =\mathrm{0.5}({\mathit{\nu }}_{\mathrm{thaw}}+{\mathit{\nu }}_{\mathrm{froz}})(\overline{a}-a)\\ & +\mathrm{0.5}\left|\right({\mathit{\nu }}_{\mathrm{thaw}}-{\mathit{\nu }}_{\mathrm{froz}}\left)\right(\overline{a}-a\left)\right|,\end{array}$$ where ${\mathit{\nu }}_{\mathrm{thaw}}$ is the rate of thawing and ${\mathit{\nu }}_{\mathrm{froz}}$ is the rate of freezing. Because ${\mathit{\nu }}_{\mathrm{thaw}}>{\mathit{\nu }}_{\mathrm{froz}}$, the absolute value in the equation leads to the right-hand side being ${\mathit{\nu }}_{\mathrm{thaw}}(\overline{a}-a)$ if $a<\overline{a}$, or ${\mathit{\nu }}_{\mathrm{froz}}(\overline{a}-a)$ if $a>\overline{a}$. The change in the pool of frozen carbon ($C_{\mathrm{fr}}$) naturally follows: 42$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{fr}}}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{d}a}{\mathrm{d}t}}{C}_{\mathrm{fr}\mathrm{0}},$$ where $C_{\mathrm{fr}\mathrm{0}}$ is the amount of frozen carbon at preindustrial times.

Thawed carbon is not directly emitted to the atmosphere: it is split into three thawed carbon sub-pools ($C_{\mathrm{th},j}$) that have their own decay time but are all affected by an additional function ($r_{\mathrm{rt}}$). This leads to the following budget equations: 43$${\displaystyle \frac{\mathrm{d}{C}_{\mathrm{th},j}}{\mathrm{d}t}}=-{\mathit{\alpha }}_{\mathrm{th},j}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{fr}}}{\mathrm{d}t}}-{\displaystyle \frac{C_{\mathrm{th},j}}{{\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{th}}}{\mathit{\tau }}_{\mathrm{th},j}}}{r}_{\mathrm{rt}},\forall j,$$ where ${\mathit{\alpha }}_{\mathrm{th},j}$ is the sub-pools' splitting shares (with ${\sum }_j{\mathit{\alpha }}_{\mathrm{th},j}=\mathrm{1}$), ${\mathit{\tau }}_{\mathrm{th},j}$ is the sub-pools' decay times, and ${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{th}}}$ is a scaling factor applied to all sub-pools. The additional $r_{\mathrm{rt}}$ function describes the sensitivity of heterotrophic respiration to climate change in boreal regions using a Gaussian formula: 44$$r_{\mathrm{rt}}=\exp \left({\mathit{\kappa }}_{\mathrm{rt}}{\mathit{\gamma }}_{\mathrm{rt}\mathrm{1}}{\mathit{\alpha }}_{\mathrm{lst}}T-{\mathit{\kappa }}_{\mathrm{rt}}{\mathit{\gamma }}_{\mathrm{rt}\mathrm{2}}({\mathit{\alpha }}_{\mathrm{lst}}T{)}^{\mathrm{2}}\right),$$ where ${\mathit{\kappa }}_{\mathrm{rt}}$ is a factor scaling the sensitivity of thawed carbon against that of regular soil carbon, ${\mathit{\gamma }}_{\mathrm{rt}\mathrm{1}}$ is the sensitivity to local temperature change (i.e. a Q10), and ${\mathit{\gamma }}_{\mathrm{rt}\mathrm{2}}$ is the quadratic term in the latter sensitivity that represents a saturation effect. Noting that all the emitted carbon is assumed to be ${\mathrm{CO}}_{\mathrm{2}}$, the global emission from permafrost ($E_{\mathrm{pf}}$) is thus 45$$E_{\mathrm{pf}}=\sum _j{\displaystyle \frac{C_{\mathrm{th},j}}{{\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{th}}}{\mathit{\tau }}_{\mathrm{th},j}}}{r}_{\mathrm{rt}}.$$

2.7

Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$

The change in atmospheric concentration of ${\mathrm{CO}}_{\mathrm{2}}$ is the budget of all carbon cycle fluxes to which we add the exogenous anthropogenic emissions ($E_{{\mathrm{CO}}_{\mathrm{2}}}$): 46$${\mathit{\alpha }}_{\mathrm{C}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}t}}=E_{{\mathrm{CO}}_{\mathrm{2}}}+E_{\mathrm{pf}}-F_{\mathrm{land}}-F_{\mathrm{ocean}},$$ where ${\mathit{\alpha }}_{\mathrm{C}}$ is the conversion factor from volume fraction to mass for ${\mathrm{CO}}_{\mathrm{2}}$.

3 Bayesian calibration

3.1 Principle

Bayesian inference is a powerful tool for assimilating observational data into reduced-complexity models such as Pathfinder . The approach consists in deducing probability distributions of parameters from a priori knowledge on those distributions and on distributions of observations of some of the model's variables, using Bayes' theorem . Summarily, the Bayesian calibration updates the joint distribution of parameters to make it as compatible with the constraints as possible given their prior estimates, which increases the internal coherence of Pathfinder by excluding combination of parameters that are unlikely.

Such a Bayesian calibration is vulnerable to the possibility that the priors draw on the same information as the constraints. However, given that Pathfinder is a patchwork of emulators whose parameters are obtained independently from one another and following differing experimental setups, we expect that the coherence of information contained within the priors and the constraints is very low. Our choice of using only complex models as prior information and only observations and assessments as constraints also aims at limiting this vulnerability.

Concretely, the posterior probability ${\mathcal{P}}_{\mathrm{post}}$ of a sample $k$ from the joint parameters distribution ${\mathit{\xi }}_k$, conditional to a set of observations $\mathit{x}$, is proportional (symbol $\propto$) to its own prior probability ${\mathcal{P}}_{\mathrm{pre}}$ and to the likelihood $\mathcal{L}$ of the model simulating $\mathit{x}$ given ${\mathit{\xi }}_k$:

47$${\mathcal{P}}_{\mathrm{post}}\left({\mathit{\xi }}_k\mathrm{|}\mathit{x}\right)\propto \mathcal{L}\left(\mathit{x}\mathrm{|}{\mathit{\xi }}_k\right){\mathcal{P}}_{\mathrm{pre}}\left({\mathit{\xi }}_k\right).$$ Here, we assume all observations are independently and identically distributed following a normal distribution (with mean values ${\mathit{\mu }}_x$, and standard deviations ${\mathit{\sigma }}_x$ expressed in real physical units), which leads to the following likelihood: 48$$\mathcal{L}\left(\mathit{x}\mathrm{|}{\mathit{\xi }}_k\right)=\prod _{i=\mathrm{1}}^{n_x}{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{x,i}\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{\left({\mathcal{F}}_i\right({\mathit{\xi }}_k)-{\mathit{\mu }}_{x,i}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{x,i}^{\mathrm{2}}}}\right),$$ where ${\mathcal{F}}_i\left({\mathit{\xi }}_k\right)$ is the model's output for the $i$th observable (out of $n_x$) with input parameters ${\mathit{\xi }}_k$.

3.2 Implementation

The Pathfinder model is a set of differential equations with a number of input parameters, of which $n_{\mathit{\xi }}$ are calibrated through Bayesian inference, and an additional two input variables provided as time series (i.e. one value per time step required). While the two input time series can be any combination of two out of four variables (anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions, non-${\mathrm{CO}}_{\mathrm{2}}$ ERF, atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration, or GMST), for calibration we use the two most well-constrained variables that are direct physical observations of the global Earth system: atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and GMST. These input time series cover the historical period from 1751 to 2020. Therefore, the ${\mathit{\xi }}_k$ vector is

49$${\mathit{\xi }}_k={\left\{\mathit{\{}{\mathit{\xi }}_j{\mathit{\}}}_{j=\mathrm{1}}^{n_{\mathit{\xi }}},\mathit{\{}C\left(t\right){\mathit{\}}}_{t=\mathrm{1751}}^{\mathrm{2021}},\mathit{\left\{}T\right(t){\mathit{\}}}_{t=\mathrm{1751}}^{\mathrm{2021}}\right\}}_k.$$ However, to ease the computation by reducing the dimension of the system, we do not use annual time series of observations as inputs, but we assume that each input time series (for variable $X$ being $C$ or $T$) follows 50$$X\left(t\right)=X_{\mathit{\mu }}\left(t\right)+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_X{X}_{\mathit{\sigma }}\left(t\right)+{\mathit{ϵ}}_X\mathrm{AR}\mathrm{1}(t;{\mathit{\rho }}_X),$$ where $X_{\mathit{\mu }}$ and $X_{\mathit{\sigma }}$ are fixed exogenous annual time series (i.e. structural parameters), ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_X$ is the relative standard deviation of the time series (without noise), ${\mathit{ϵ}}_X$ is the noise intensity, and $\mathrm{AR}\mathrm{1}$ is an autoregressive process of order 1 and autocorrelation parameter ${\mathit{\rho }}_X$. This assumption leads to the final expression of the ${\mathit{\xi }}_k$ vector: 51$${\mathit{\xi }}_k={\left\{\mathit{\{}{\mathit{\xi }}_j{\mathit{\}}}_{j=\mathrm{1}}^{n_{\mathit{\xi }}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{C}},{\mathit{ϵ}}_{\mathrm{C}},{\mathit{\rho }}_{\mathrm{C}},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_T,{\mathit{ϵ}}_T,{\mathit{\rho }}_T,\right\}}_k.$$

During this Bayesian assimilation, the Pathfinder model is run solely over the historical period (from 1750 to 2021), as the constraints concern only preindustrial or historical years. For the computation, the time-differential system of Pathfinder is solved using an implicit–explicit numerical scheme (also called IMEX), with a time step of a quarter of a year. This solving scheme relies on writing the differential equations of all state variables $X$ as 52$${\displaystyle \frac{\mathrm{d}X}{\mathrm{d}t}}=-\mathit{\nu }X+\mathcal{R},$$ where $\mathit{\nu }$ is the constant speed of the linear part of the differential equation and $\mathcal{R}$ is its non-linear part, discretizing these equations as 53$${\displaystyle \frac{X(t+\mathit{\delta }t)-X\left(t\right)}{\mathit{\delta }t}}=-\mathit{\nu }X(t+\mathit{\delta }t)+\mathcal{R}\left(t\right),$$ where $\mathit{\delta }t$ is the solving time step (which is $\frac{\mathrm{1}}{n_t}$ times the annual time step of the model's inputs and outputs, here $n_t=\mathrm{4}$), and finally explicitly solving for all $X(t+\mathit{\delta }t)$. We note this is also the default solving scheme for regular simulations with the model, although the value of $n_t$ can be altered and alternative schemes are available.

The Bayesian procedure itself is implemented using the Python computer language and specifically the PyMC3 package . The solving of Eq. () and its normalization are done using the package's full-rank automatic differentiation variational inference (ADVI) algorithm , with 100 000 iterations (and default algorithm options). The choice of variational inference instead of Markov chain Monte Carlo is motivated by the significant size of our model and the speed of ADVI. An additional strength of the full-rank version of the ADVI algorithm is its ability to generate correlated posterior distributions even if the prior ones are uncorrelated. Convergence of the algorithm was controlled through convergence of the ELBO metric . All results presented hereafter are obtained through drawing 2000 sets of parameters – which we call configurations – from the posterior or prior distributions.

3.3 Constraints

We use a set of 19 constraints related to all aspects of the model that correspond to the set of observations $\mathit{x}$ in the Bayesian calibration. Many of the constraints are observations, but some are ranges assessed by expert panels such as the Global Carbon Project or the IPCC. They cover either a recent point in time or an assumed preindustrial equilibrium, and they are typically taken over a period of at least a few years to reduce the effect of natural variability.

Table  summarizes these constraints, the periods over which they are considered, and their distributions. The following subsections provide further details on the constraints, and the constraint distributions are shown in Fig. .

Table 1

Constrained variables in Pathfinder, with values before and after calibration. Variables are noted under their text notation, and Tables  and provide the corresponding notation in code. The uncertainty correspond to the 1$\mathit{\sigma }$ uncertainty range.

Variable	Period	Method	Prior	Posterior	Constraints	Unit
$E_{{\mathrm{CO}}_{\mathrm{2}}}$	2011–2020	Mean	9.1 $\pm$ 1.3	10.0 $\pm$ 0.7	11.0 $\pm$ 0.9	PgC yr${}^{-\mathrm{1}}$
$\frac{\mathrm{d}C}{ct}$	2011–2020	Mean	2.41 $\pm$ 0.06	2.40 $\pm$ 0.01	2.40 $\pm$ 0.01	ppm yr${}^{-\mathrm{1}}$
$F_{\mathrm{land}}$	1960–2020	Sum	95 $\pm$ 52	123 $\pm$ 26	145 $\pm$ 35	PgC
$F_{\mathrm{ocean}}$	1960–2020	Sum	89$\pm$ 12	97 $\pm$ 13	105 $\pm$ 20	PgC
$C_{\mathrm{v}}$	1750	Mean	407 $\pm$ 54	407 $\pm$ 37	450 $\pm$ 50	PgC
$C_{\mathrm{s}}$	1750	Mean	1181 $\pm$ 735	1086 $\pm$ 284	1088 $\pm$ 249	PgC
$F_{\mathrm{NPP}}$	1998–2002	Mean	60.0 $\pm$ 7.9	59.5 $\pm$ 3.9	55.0 $\pm$ 5.0	PgC yr${}^{-\mathrm{1}}$
$C$	2012–2021	Mean	403.6 $\pm$ 0.3	403.6 $\pm$ 0.1	403.6 $\pm$ 0.1	ppm
$R_x$	2010–2019	Mean	0.01 $\pm$ 0.47	0.33 $\pm$ 0.37	0.56$\pm$ 0.53	W m${}^{-\mathrm{2}}$
$T$	2001–2020	Mean	0.96 $\pm$ 0.08	0.97 $\pm$ 0.06	1.00 $\pm$ 0.07	K
$\frac{\mathrm{d}T}{\mathrm{d}t}$	2000–2019	Mean	0.028 $\pm$ 0.003	0.028 $\pm$ 0.002	0.029 $\pm$ 0.002	K yr${}^{-\mathrm{1}}$
$\frac{\mathrm{d}{U}_{\mathrm{ohc}}}{\mathrm{d}t}$	2006–2018	Mean	16 $\pm$ 0.10	0.62 $\pm$ 0.09	12 $\pm$ 0.17	W m${}^{-\mathrm{2}}$
$\frac{\mathrm{d}{H}_{\mathrm{thx}}}{\mathrm{d}t}$	2006–2018	Mean	1.02 $\pm$ 0.22	1.14 $\pm$ 0.21	1.39 $\pm$ 0.40	mm yr${}^{-\mathrm{1}}$
$\frac{\mathrm{d}{H}_{\mathrm{gla}}}{\mathrm{d}t}$	2006–2018	Mean	0.63 $\pm$ 0.24	0.62 $\pm$ 0.04	0.62 $\pm$ 0.03	mm yr${}^{-\mathrm{1}}$
$\frac{\mathrm{d}{H}_{\mathrm{ais}}}{\mathrm{d}t}$	2006–2018	Mean	$-$0.02 $\pm$ 0.23	0.30 $\pm$ 0.10	0.37 $\pm$ 0.08	mm yr${}^{-\mathrm{1}}$
$\frac{\mathrm{d}{H}_{\mathrm{gis}}}{\mathrm{d}t}$	2006–2018	Mean	0.36 $\pm$ 0.12	17 $\pm$ 0.10	0.63 $\pm$ 0.07	mm yr${}^{-\mathrm{1}}$
$H_{\mathrm{tot}}$	1901–1990	Difference	72 $\pm$ 17	83 $\pm$ 10	89 $\pm$ 32	mm
$H_{\mathrm{lia}}$	1750	Mean	45 $\pm$ 17	45 $\pm$ 11	30 $\pm$ 13	mm
$\mathrm{logit}\left(\mathrm{ff}\right)$	1750	Mean	1.69 $\pm$ 0.38	1.47 $\pm$ 0.28	1.38 $\pm$ 0.37	1

3.3.1 Climate system

To constrain the temperature response, we use the same five data sets of observed GMST as in Sect.  to derive average and standard deviation of two constraints: the average GMST change and the average GMST yearly trend obtained through second-order accuracy gradient , both over the latest 20 years of data (2002–2021). Because these data sets are already used as input to the Bayesian setup, albeit in a different way, they do not provide much of a constraint and are used mostly to ensure that the ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_T$ and ${\mathit{ϵ}}_T$ parameters remain within sensible range.

To further constrain the climate system, we use the mean OHU assessed by the IPCC AR6 over 2006–2018 Table 2.7 and the non-${\mathrm{CO}}_{\mathrm{2}}$ ERF (averaged over 2010-2019) also estimated for the AR6. The central value of the latter is taken from Table AIII.3, and corresponding GitHub repository, and its uncertainty is constructed using data from Table 7.8 and assuming the ERF of all species are normally distributed and uncorrelated, but fully correlated in time for each separate species (which likely overestimates the uncertainty).

To better align with the IPCC AR6, we also constrain the ECS of our model (i.e. the $T_{\mathrm{2}\times }$ parameters). To do so, because the distribution of ECS cannot be assumed to be normal, we follow the framework of , who define the climate feedback factor $\mathrm{ff}$ so that $T_{\mathrm{2}\times }=T_{\mathrm{2}\times }^{*}/(\mathrm{1}-\mathrm{ff})$, where $T_{\mathrm{2}\times }^{*}$ is the minimal ECS value (roughly corresponding to the Planck feedback). We assume this feedback factor follows a logit-normal distribution, which implies $\mathrm{logit}\left(\mathrm{ff}\right)$ $=$ $\ln (\mathrm{ff}/(\mathrm{1}-\mathrm{ff}\left)\right)$ $=$ $\ln (T_{\mathrm{2}\times }/T_{\mathrm{2}\times }^{*}-\mathrm{1})$ follows a normal distribution. We therefore constrain $\mathrm{logit}\left(\mathrm{ff}\right)$, using distribution parameters and a value of $T_{\mathrm{2}\times }^{*}$ calibrated to fit the probabilistic ranges of ECS provided by the AR6. This fit of the ECS distribution is illustrated in Fig. .

3.3.2 Carbon cycle

Similarly to what is done with GMST, we constrain the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ level over the latest 10 years of data (2012–2021) using the NOAA/ESRL data . The rest of the global ${\mathrm{CO}}_{\mathrm{2}}$ budget is constrained using the 2021 Global Carbon Budget (GCB; ). We use the net atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ growth and total anthropogenic emissions (fossil and land use) over the last 10 years and the ocean and land carbon sinks accumulated since the beginning of the instrumental measurement period (1960–2020). Note that our definition of the land carbon sink ignoring land use change is consistent with that of the GCB.

Given its number of parameters and their inconsistent sources, we further constrain the land carbon module by considering present-day (mean over 1998–2002) NPP and preindustrial vegetation and soil carbon pools. These preindustrial pools are taken from the AR6 for the central value Fig. 5.12, but their relative uncertainty is taken from the AR5 Fig. 6.1 since it is lacking in the AR6. In addition, the soil carbon pool constraint is corrected downward by estimates of peatland carbon Table 1, since it is an ecosystem missing in TRENDY models (and in ours) but not in the IPCC assessments.

3.3.3 Sea level rise

To constrain the separate SLR contributions from thermal expansion, GIS, AIS, and glaciers, we use the model-based SLR speed estimates over the recent past (averaged over 2006–2018) reported in the AR6 Table 9.5. To constrain the total contribution, we also use the historical (1901–1990) sea level rise inferred from tide gauges from the same source, although the value is corrected upward for the missed impact of uncharted glaciers .

Contrarily to all other modules, the SLR module is not assumed to start at steady state in 1750, which is represented through the ${\mathit{\lambda }}_{\mathrm{ice}}$ ($\mathrm{ice}\in [\mathrm{gla},\mathrm{gis},\mathrm{ais}]$) parameters. We assume this is entirely due to the so-called Little Ice Age (LIA) relaxation, which we assume can be simply modelled in Pathfinder through exponential decay of our three ice-related contributions since $t_{\mathrm{0}}=\mathrm{1750}$. This gives a net LIA contribution of $H_{\mathrm{lia}}={\sum }_{\mathrm{ice}}{\mathit{\lambda }}_{\mathrm{ice}}{\mathit{\tau }}_{\mathrm{ice}}\exp \left(-\frac{t-t_{\mathrm{0}}}{{\mathit{\tau }}_{\mathrm{ice}}}\right)$. We constrain this diagnostic variable using the global SLR reported by over 1900–2005 for their control experiment.

3.4 Parameters (prior distributions)

Out of the model's 77 parameters, 33 are assumed to be fixed (i.e. they are structural parameters), and the remaining $n_{\mathit{\xi }}=\mathrm{44}$ parameters are estimated through Bayesian inference. Prior distributions of the ${\mathit{\xi }}_j$ parameters are assumed log-normal if the physical parameter must be defined positive, logit-normal if it must be between 0 and 1, and normal otherwise. To avoid extreme parameter values that could make the model diverge during calibration, the posterior distributions are bound to ${\mathit{\mu }}_{\mathit{\xi },j}\pm \mathrm{5}{\mathit{\sigma }}_{\mathit{\xi },j}$, where ${\mathit{\mu }}_{\mathit{\xi },j}$ and ${\mathit{\sigma }}_{\mathit{\xi },j}$ are the mean and standard deviation of the $j$th parameter's prior distribution. These two values are taken from the literature, deduced from multi-model ensembles, or in a few instances arbitrarily set, as described in the following subsections. Note that when parameters are deduced from multi-model ensembles, there are effectively two rounds of calibration: first, a calibration on individual models using ordinary least-squares regressions to obtain prior distributions and, second, the Bayesian calibration itself that leads to the posterior distributions. In addition, the prior distributions of ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_X$, ${\mathit{ϵ}}_X$ and ${\mathit{\rho }}_X$ are assumed normal, half-normal, and uniform, respectively. All prior distributions are assumed independent, meaning that the prior joint distribution $\mathit{\xi }$ does not exhibit any covariance.

All parameters are summarized in Tables  and along with their properties and values. The following subsections further explain how the prior distributions of the parameters are established, and these distributions are shown in Fig. .

Figure 5

Parameter distributions before (black lines) and after (blue lines) the Bayesian calibration. Parameters are noted under their text notation, and Tables  and provide the corresponding notation in code.

[Figure omitted. See PDF]

3.4.1 Climate

All the parameters of the climate module are calibrated. The prior distribution of the radiative parameter $\mathit{\varphi }$ is taken from the AR5 Table 8.SM.1. All other prior distributions of the parameters of the climate module (i.e. $T_{\mathrm{2}\times }$, ${\mathrm{\Theta }}_{\mathrm{s}}$, ${\mathrm{\Theta }}_{\mathrm{d}}$, $\mathit{\theta }$, and ${\mathit{ϵ}}_{\mathrm{heat}}$) are taken from 35 CMIP6 models whose climate responses were derived for the AR6 using the “abrupt-4xCO2” experiment Sect. 7.SM.2.1, and corresponding GitHub repository. Here, $T_{\mathrm{2}\times }$ is simply assumed to be half the reported equilibrium temperature at quadrupled ${\mathrm{CO}}_{\mathrm{2}}$. In addition, the prior distribution of the ocean warming fraction ${\mathit{\alpha }}_{\mathrm{ohc}}$ is taken from the AR6 Table 7.1.

3.4.2 Sea level rise

Some parameters from the SLR module are structural: the maximum SLR contribution from glaciers (${\mathrm{\Lambda }}_{\mathrm{gla}}$) is taken from Sect. 9.6.3.2; the equilibrium AIS SLR (${\mathrm{\Lambda }}_{\mathrm{ais}}$) is from Fig. 13.14; and the ${\mathit{\tau }}_{\mathrm{gis}}$, ${\mathit{\tau }}_{\mathrm{gla}}$, and ${\mathit{\tau }}_{\mathrm{ais}}$ timescales are the mean values from Table S1 (assuming they provide the 90 % range of a log-normal distribution). All other parameters are calibrated. The prior distribution of the thermosteric parameter ${\mathrm{\Lambda }}_{\mathrm{thx}}$ is taken from the AR6 Sect. 9.2.4.1, as are the prior distributions of the preindustrial offset parameters ${\mathit{\lambda }}_{\mathrm{gis}}$, ${\mathit{\lambda }}_{\mathrm{gla}}$ and ${\mathit{\lambda }}_{\mathrm{ais}}$ earliest period of Table 9.5. For the remaining parameters, we derive prior distributions using SLR projections compiled by for a number of ice sheets and glaciers models over various Representative Concentration Pathway (RCP) and Shared Socioeconomic Pathway (SSP) scenarios. Using the models' outputs, we apply Eq. () to estimate the ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}$ and ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{3}}$ parameters; Eq. () for the ${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{1}}$, ${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{3}}$, and ${\mathit{\gamma }}_{\mathrm{gla}}$ parameters; and Eq. () for the ${\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}$ and ${\mathit{\alpha }}_{\mathrm{ais}}$ parameters. During these fits, all other parameters are assumed to take their default value if structural, and their best-guess value otherwise. Results of this calibration on the individual models compiled by are shown for each SLR contribution in Figs. , and .

3.4.3 Ocean carbon

The ocean carbon cycle module has a number of structural parameters: ${\mathit{\alpha }}_{\mathrm{dic}}$, all ${\mathit{\alpha }}_{\mathrm{o},j}$, and all ${\mathit{\tau }}_{\mathrm{o},j}$ are taken from Tables A2 and A3, based on the Princeton model. The prior distribution of the adjustment factor ${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{o}}}$ is arbitrarily taken to apply a 20 % uncertainty on the oceanic transport timescales. All other prior distributions for this module's parameters are derived from 12 CMIP6 models with interactive carbon cycle that contributed to C4MIP . $T_{\mathrm{o}}$ is taken on average over the piControl simulation. ${\mathit{\nu }}_{\mathrm{gx}}$ and ${\mathit{\gamma }}_{\mathrm{gx}}$ are calibrated by applying Eq. () to the models' outputs for the 1pctCO2, 1pctCO2-rad, and 1pctCO2-bgc experiments, while ${\mathit{\beta }}_{\mathrm{dic}}$ and ${\mathit{\gamma }}_{\mathrm{dic}}$ are calibrated by applying Eqs. ()–() and (). Results of this calibration on the individual CMIP6 models are shown in Figs.  and .

3.4.4 Ocean acidification

In this version of Pathfinder, ${\mathit{\kappa }}_{\mathrm{pH}}$ is a structural parameter set to 1.

3.4.5 Land carbon

Parameters related to the passive soil carbon pool are taken from Table S5: ${\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}$ is structural, while ${\mathit{\alpha }}_{\mathrm{pass}}$ is not. All the prior distribution of the parameters related to the preindustrial steady state of the land carbon (i.e. $F_{\mathrm{npp}\mathrm{0}}$, ${\mathit{\nu }}_{\mathrm{fire}}$, ${\mathit{\nu }}_{\mathrm{harv}}$, ${\mathit{\nu }}_{\mathrm{mort}}$, ${\mathit{\nu }}_{\mathrm{stab}}$, ${\mathit{\nu }}_{\mathrm{rh}\mathrm{1}}$, and ${\mathit{\nu }}_{\mathrm{cs}}$) are derived from 11 TRENDYv7 models , exactly as for OSCAR v3.1 except that all biomes and regions are lumped together. The prior distribution of the remaining parameters are derived from 12 CMIP6 models that contributed to C4MIP . Using the models' outputs for the 1pctCO2, 1pctCO2-rad, and 1pctCO2-bgc experiments, we calibrated ${\mathit{\beta }}_{\mathrm{npp}}$, ${\mathit{\alpha }}_{\mathrm{npp}}$, and ${\mathit{\gamma }}_{\mathrm{npp}}$ through Eq. (); ${\mathit{\beta }}_{\mathrm{fire}}$ and ${\mathit{\gamma }}_{\mathrm{fire}}$ through Eq. (); and ${\mathit{\beta }}_{\mathrm{rh}}$ and ${\mathit{\gamma }}_{\mathrm{rh}}$ through Eq. (). Results of this calibration on the individual CMIP6 models is shown in Figs. , , and .

3.4.6 Permafrost carbon

The permafrost module's parameters are recalibrated using the same algorithm as used by but adapted to the global formulation of Pathfinder. First, the algorithm is run once to obtain a set of parameters reproducing the behaviour of the global average of five permafrost models (with data from UVic added to the four original models). This gives the values of the structural parameters (i.e. ${\mathit{\alpha }}_{\mathrm{lst}}$, ${\mathit{\gamma }}_{\mathrm{rt}\mathrm{1}}$, ${\mathit{\gamma }}_{\mathrm{rt}\mathrm{2}}$, ${\mathit{\kappa }}_{\mathrm{rt}}$, $a_{\min }$, all ${\mathit{\alpha }}_{\mathrm{th},j}$, all ${\mathit{\tau }}_{\mathrm{th},j}$, ${\mathit{\nu }}_{\mathrm{thaw}}$, and ${\mathit{\nu }}_{\mathrm{froz}}$). Second, the algorithm is run five additional times for each of the five permafrost models separately, with the structural parameters established in the first step, to obtain prior distributions of the remaining parameters (i.e. $C_{\mathrm{fr}\mathrm{0}}$, ${\mathit{\kappa }}_a$, ${\mathit{\gamma }}_a$, and ${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{th}}}$).

3.4.7

Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$

The conversion factor ${\mathit{\alpha }}_{\mathrm{C}}$ is a structural parameter whose value is taken from the latest GCBs e.g.. The prior distribution of preindustrial ${\mathrm{CO}}_{\mathrm{2}}$ concentration ($C_{\mathrm{pi}}$) is taken from the AR6 Sect.  2.2.3.2.1, assuming the difference between minimum and maximum over the 0–1850 period is representative of the 90 % uncertainty range.

3.4.8

Historical ${\mathrm{CO}}_{\mathrm{2}}$ and GMST

The structural $X_{\mathit{\mu }}$ and $X_{\mathit{\sigma }}$ time series are taken from the latest observations, as follows. $T_{\mathit{\mu }}$ and $T_{\mathit{\sigma }}$ are taken as the average and standard deviation of five observational GMST data sets: HadCRUT5 , Berkeley Earth , GISTEMP , NOAAGlobalTemp , and . We use the 1850–1900 period to define our preindustrial baseline, and GMST change is assumed to be zero before the earliest date available in each data set. Regarding atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, $C_{\mathit{\mu }}$ is taken as the global value reported by NOAA/ESRL and $C_{\mathit{\sigma }}$ as a constant $\pm$1 $\mathrm{ppm}$ uncertainty for 1980 onward (this uncertainty is arbitrarily taken higher than the actual uncertainty estimated through instrumental measures to increase freedom in the calibration). Before that period, $C_{\mathit{\mu }}$ comes from the IPCC AR6 Table AIII.1a, and $C_{\mathit{\sigma }}$ is linearly interpolated backwards from the instrumental uncertainty in 1980 to the preindustrial one in 1750. Finally, the prior distribution of ${\mathit{\rho }}_X$ is set to uniform over $[\mathrm{0},\mathrm{1}]$, that of ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_X$ is a unit normal distribution, and that of ${\mathit{ϵ}}_X$ is set arbitrarily to a half-normal of parameter 0.05 $\mathrm{K}$ for GMST and 0.5 $\mathrm{ppm}$ for ${\mathrm{CO}}_{\mathrm{2}}$.

3.5 Results (posterior distributions)

The following subsections discuss the adjustments between the prior and posterior parameters that are the results of the Bayesian calibration, as well as the matching of the constraints. These sections constantly refer to Fig.  that shows the prior and posterior distributions of the model's parameters, Fig.  that shows those of the constrained variables, and Fig.  that displays the correlation matrix of the posterior parameters (there is no correlation among the prior parameters). Prior and posterior values of the parameters can also be retrieved from Table .

Figure 6

Distributions of the constrained variables. Dashed lines give the distributions used to constrain. Black lines give the distribution before calibration, while blue lines give the distribution posterior to calibration. Under a variable's name, we give the period over which the constraint is estimated and the data processing method: mean over the period, difference between last and first year, or sum of all the years over the period. (1750 is the preindustrial period). Variables are noted under their text notation, and Tables  and provide the corresponding notation in code.

[Figure omitted. See PDF]

3.5.1 Climate system

Our climate-related constraints lead to adjusting all the parameters of the climate module. As explained in Sect. , the constraints for present-day GMST change and its derivative are met by construction.

The ECS ($T_{\mathrm{2}\times }$) is the parameter with the strongest adjustment, since it is directly constrained. Its precise value is discussed hereafter in Sect. , but we note that it is unsurprisingly decreased, as the CMIP6 model ensemble tends to overestimate the ECS compared to the IPCC-assessed value. Consequently, our posterior $\mathrm{logit}\left(\mathrm{ff}\right)$ matches the constraint well. The adjustment of the ECS significantly reduces the gap between our posterior distribution of the non-${\mathrm{CO}}_{\mathrm{2}}$ ERF and its constraint, although the posterior central value remains 41 % lower (but well within uncertainty range).

Among the dynamic parameters that are adjusted, we note that the deep-ocean heat capacity ${\mathrm{\Theta }}_{\mathrm{d}}$ is somewhat increased compared to the prior and that the heat exchange coefficient $\mathit{\theta }$ is also increased. These dynamic parameters are likely adjusted through our OHU constraint that is corrected in the posterior so the difference in the central values is lowered from 22 % to 14 %, which remains well within the uncertainty range.

Table 2

Diagnostics of climate and carbon cycle responses in Pathfinder before and after Bayesian calibration. Comparison with AR6 and CMIP6 data is shown. For AR6 data we give the median and 90 % confidence interval, while for every other value we give the mean $\pm$1$\mathit{\sigma }$.

2$\times$ CO${}_{\mathrm{2}}$	Pathfinder unconstrained	Pathfinder constrained	CMIP6	AR6
ECS (K)	4.1 $\pm$ 1.3	3.3 $\pm$ 0.7	3.7 $\pm$ 1.1	3.0 (2.0, 4.5)
TCR (K)	2.2 $\pm$ 0.5	1.9 $\pm$ 0.3	2.0 $\pm$ 0.4	1.8 (1.4, 2.2)
TCRE (K EgC${}^{-\mathrm{1}}$)	2.20 $\pm$ 0.63	1.65 $\pm$ 0.32	1.77 $\pm$ 0.37	1.65 (1.0, 2.3)
${\mathit{\beta }}_{\mathrm{ocean}}$ (PgC ppm${}^{-\mathrm{1}}$)	0.81 $\pm$ 0.10	0.87 $\pm$ 0.11	0.91 $\pm$ 0.09
${\mathit{\gamma }}_{\mathrm{ocean}}$ (PgC K${}^{-\mathrm{1}}$)	$-$12.9 $\pm$ 5.4	$-$12.5 $\pm$ 6.0	$-$8.6 $\pm$ 2.9
${\mathit{\beta }}_{\mathrm{land}}$ (PgC ppm${}^{-\mathrm{1}}$)	1.05 $\pm$ 0.5	1.26 $\pm$ 0.30	1.22 $\pm$ 0.40
${\mathit{\gamma }}_{\mathrm{land}}$ (PgC K${}^{-\mathrm{1}}$)	$-$33.2 $\pm$ 26.6	$-$25.3 $\pm$ 24.2	$-$34.1 $\pm$ 38.4

In addition, a number of weak but physically meaningful correlations across the climate module's parameters are found, such as a positive correlation between $T_{\mathrm{2}\times }$ and ${\mathit{ϵ}}_{\mathrm{heat}}$ see, e.g., a positive correlation between $T_{\mathrm{2}\times }$ and ${\mathrm{\Theta }}_{\mathrm{d}}$ (that tends to exclude configuration that would warm fast and high), and a negative correlation between $T_{\mathrm{2}\times }$ and $\mathit{\varphi }$ (to match the GMST and ERF constraints together).

3.5.2 Carbon cycle

Similarly to GMST, the posterior distribution of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration matches the constraint by construction. Its derivative, however, is (slightly overly) corrected to match the GCB estimate. Global anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions are significantly increased to get closer to the GCB constraint, but their central value remains 9 % too low. Since these emissions are determined through mass balance and the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ matches observations, this implies that the total carbon sinks (i.e. $F_{\mathrm{land}}+F_{\mathrm{ocean}}$) must be weaker.

This is confirmed for the ocean sink, as the posterior central value of $F_{\mathrm{ocean}}$ is 8 % lower than the constraint but still noticeably corrected if compared to the prior. This correction is explained by small adjustments in some parameters of the ocean carbon module. The mixed-layer depth is slightly increased through ${\mathit{\beta }}_{\mathrm{dic}}$. All other parameters remain mostly unaffected by the calibration, and only minor correlations are found. These results, along with the fact that our prior distribution spans only about half of the constraint's distribution, suggest that there is a structural limitation in our ocean carbon module that warrants further investigation.

It is also confirmed that the posterior land sink is weaker than the constraint, by 15 % for the central value, which is nevertheless a significant reduction of the prior gap of 34 %. To explain this adjustment, we observe that the ${\mathrm{CO}}_{\mathrm{2}}$ fertilization sensitivities (${\mathit{\beta }}_{\mathrm{npp}}$ and ${\mathit{\gamma }}_{\mathrm{npp}}$) are adjusted upwards. However, our constraint on present-day NPP prevents these adjustments from being too important, as the posterior distribution of this variable is similar to the prior and its central value remains 8 %–9 % higher than its constraint. An increased preindustrial NPP mechanically leads to an increase in preindustrial carbon pools, but these require further adjustments of the land carbon turnover rates, and most notably the mortality rate ${\mathit{\nu }}_{\mathrm{mort}}$ and the passive carbon fraction ${\mathit{\alpha }}_{\mathrm{pass}}$, to better match their constraints (of which the one on total soil carbon is perfectly met).

Table 3

Comparison of SSP scenarios for GMST change projections (w.r.t. 1850–1900) to AR6 Table 4.5, and for ocean and land carbon storage projections to CMIP6 . Land carbon storage projections were corrected using the land use change emissions data from SSPs . For GMST data we give the median and the 90 % confidence interval while for every other values we give the mean $\pm$1$\mathit{\sigma }$.

Experiment	Model	GMST	GMST	GMST	Ocean carbon	Land carbon
		(K)	(K)	(K)	storage (PgC)	storage (PgC)
		2021–2040	2041–2060	2081–2100	2015–2100	2015–2100
ssp119	Pathfinder	1.5 (1.3, 1.8)	1.6 (1.4, 1.9)	1.5 (1.2, 1.7)	132 $\pm$ 21	49 $\pm$ 33
ssp119	AR6 or CMIP6	1.5 (1.2, 1.7)	1.6 (1.2, 2.0)	1.4 (1.0, 1.8)	111 $\pm$ 11	73 $\pm$ 33
ssp126	Pathfinder	1.5 (1.3, 1.8)	1.8 (1.6, 2.1)	1.9 (1.6, 2.2)	179 $\pm$ 28	109 $\pm$ 45
ssp126	AR6 or CMIP6	1.5 (1.2, 1.8)	1.7 (1.3, 2.2)	1.8 (1.3, 2.4)	162 $\pm$ 8	120 $\pm$ 50
ssp245	Pathfinder	1.6 (1.4, 1.8)	2.1 (1.8, 2.4)	2.8 (2.4, 3.3)	265 $\pm$ 41	225 $\pm$ 76
ssp245	AR6 or CMIP6	1.5 (1.2, 1.8)	2.0 (1.6, 2.5)	2.7 (2.1, 3.5)	252 $\pm$ 11	178 $\pm$ 76
ssp370	Pathfinder	1.6 (1.4, 1.8)	2.2 (1.9, 2.6)	3.7 (3.2, 4.3)	354 $\pm$ 53	330 $\pm$ 112
ssp370	AR6 or CMIP6	1.5 (1.2, 1.8)	2.1 (1.7, 2.6)	3.6 (2.8, 4.6)	338 $\pm$ 15	269 $\pm$ 124
ssp585	Pathfinder	1.7 (1.5, 1.9)	2.5 (2.2, 2.9)	4.4 (3.8, 5.2)	420 $\pm$ 63	409 $\pm$ 148
ssp585	AR6 or CMIP6	1.6 (1.3, 1.9)	2.4 (1.9, 3.0)	4.4 (3.3, 5.7)	398 $\pm$ 17	311 $\pm$ 162

Figure 7

Correlation matrix of Pathfinder's parameters after the Bayesian calibration. Parameters are classified according to the equations they are related to, i.e. climate system, sea level, ocean carbon, land carbon, and permafrost carbon. Parameters are noted under their text notation, and Tables  and provide the corresponding notation in code.

[Figure omitted. See PDF]

Figure 8

Historical time series of key variables from Pathfinder. Red lines are observations, black lines are the model's outputs before calibration, and blue lines are the same after calibration. Shaded areas and vertical bars correspond to the $\mathrm{1}\mathit{\sigma }$ uncertainty range. Temperature observations are taken from HadCRUT5 , , Berkeley Earth , GISTEMP , and NOAA/MLOST . Other sources are NOAA/ESRL , GCB 2021 , and AR6 .

[Figure omitted. See PDF]

Figure 9

Time series of GMST change, integrated land carbon uptake and integrated ocean carbon uptake for idealized experiments (abrupt-2xCO2, 1pctCO2, 1pctCO2-bgc and 1pctCO2-rad) and projections according to SSP scenarios in Pathfinder. Shaded areas give the $\mathrm{1}\mathit{\sigma }$ uncertainty range.

[Figure omitted. See PDF]

The land carbon module exhibits significant correlations among posterior parameters. This is likely a consequence of all the constraints combined as they dictate both the preindustrial steady state of the module and its transient response over the historical period. Eliminated configurations are those, for instance, that would show high initial carbon pools that are very sensitive to climate change (as these would lead to a very weak land sink) or that would exhibit a weak ${\mathrm{CO}}_{\mathrm{2}}$ fertilization effect associated with a fast turnover time (that would also lead to a weak sink).

3.5.3 Sea level rise

The prior parameters of the SLR module are the least informed of our Bayesian setup. The model initially underestimates the thermal expansion, as well as the GIS and AIS SLR rates. The calibration brings the posterior distributions closer to their respective constraints but it always remains in the lower end of the uncertainty range. The correction is done by adjusting many of the module's parameters (most notably ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}$, ${\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}$, ${\mathit{\lambda }}_{\mathrm{ais}}$ or ${\mathit{\lambda }}_{\mathrm{gla}}$) and by finding strong correlations among them (thus excluding physically unrealistic combinations).

The historical SLR is markedly corrected by the constraint: from a 19 % gap between the central values of the constraint and the prior estimate, to only 7 % after calibration. Here, we also note that the sum of individual contributions to historical SLR reported in AR6 do not match that total SLR Table 9.5, which likely has some impact on the consistency between our constraints. Finally, although the LIA relaxation contribution is not altered by the calibration, as its central value remains 50 % too high, it is the likely source of the strong correlations found among the parameters of this module because it straightforwardly links the individual SLR contributions together.

4 Model diagnosis

4.1 Historical period

Because in the Bayesian setup we do not use annual time series of observations as constraints, the posterior distributions given in Fig.  do not inform on the whole dynamic of the model over the historical period. To further diagnose the model's behaviour, Fig.  gives the time series from 1900 to 2021 of six key variables. GMST and atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ match the historical observations very well via construction of these input time series. The non-CO${}_{\mathrm{2}}$ ERF exhibits a very high variability, owing to our temperature-driven setup and the natural variability in the GMST input. Beyond that, the ERF time series is consistent with the AR6 estimates , albeit somewhat lower on average in the recent past, as seen with the posterior distribution. Consistently with the interpretation of carbon cycle variables in the calibration results, anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions, and the ocean and land carbon sinks are slightly underestimated compared to the GCB estimates . Several reasons could explain this discrepancy, from the lack of land use change in Pathfinder to the inconsistency of the GCB figures (that do not close the budget, while ours do). Nevertheless, the interest of the calibration is clearly illustrated, as the posterior uncertainty range covers observations much better than the prior one.

4.2 Idealized simulations

To complete the diagnosis of our model with common metrics used with climate and carbon models, we ran a set of standard idealized experiments, corresponding to the CMIP6 abrupt-2xCO2, 1pctCO2, 1pctCO2-bgc, and 1pctCO2-rad. A summary of these metrics' values is given in Table , and the resulting time series are shown in Fig. .

The abrupt-2xCO2 experiment sees an abrupt doubling of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, and it is used to diagnose the model's ECS that is defined as the equilibrium temperature for a doubling of the preindustrial atmospheric concentration of ${\mathrm{CO}}_{\mathrm{2}}$ (we acknowledge that it is superfluous with this version of Pathfinder since it is also a parameter). Using the GMST anomaly at the end of 1500 years of this experiment leads to an unconstrained estimate of ECS of 4.1 $\pm$ 1.3 $\mathrm{K}$ and a constrained estimate of 3.3 $\pm$ 0.7 $\mathrm{K}$. Consistently, the latter value is between the ECS value extracted from CMIP6 models that is higher (3.7 $\pm$ 1.1 $\mathrm{K}$) and the final value assessed in the AR6 that is lower (3.0 $\mathrm{K}$, with a 67 % confidence interval between 2.5 and 4.0 $\mathrm{K}$).

Using the 1pctCO2 experiment that sees a 1 % yearly increase in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, we can estimate the model's transient climate response (TCR), which is defined as the GMST change after 70 years, when atmospheric concentration ${\mathrm{CO}}_{\mathrm{2}}$ has just doubled. The CMIP6 models have a TCR of 2.0 $\pm$ 0.4 $\mathrm{K}$ . Pathfinder's unconstrained value is higher, at 2.2 $\pm$ 0.5 $\mathrm{K}$, while the constrained one goes down to 1.9 $\pm$ 0.3 $\mathrm{K}$. If we divide the TCR by the cumulative anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions compatible with the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ increase in this experiment, we obtain an estimate of the transient climate response to emissions (TCRE). Similarly to the TCR, it is higher in the unconstrained ensemble and lower in the constrained one, when compared to CMIP6 models . Both downward adjustments of the TCR and TCRE are consistent with that of ECS, with the posterior TCRE matching the AR6 assessed range very well .

To look more closely at the carbon cycle, we perform two variants of the latter experiment: in 1pctCO2-rad, atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ only has a radiative effect, as it is kept at preindustrial level for the carbon cycle, whereas in 1pctCO2-bgc, atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ only has a biogeochemical effect, as the climate system sees only the preindustrial level. These three experiments are used to calculate the carbon–concentration ($\mathit{\beta }$) and carbon–climate ($\mathit{\gamma }$) feedback metrics that measure the carbon sinks' sensitivities to changes in atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ and GMST, respectively. We apply the same method as to calculate these, which leads to metrics at the time of ${\mathrm{CO}}_{\mathrm{2}}$ doubling that are in line with CMIP6 models . As both carbon sinks were adjusted upwards by the Bayesian calibration, the constraints logically increased both ${\mathit{\beta }}_{\mathrm{ocean}}$ and ${\mathit{\beta }}_{\mathrm{land}}$, to values fairly close to those of the complex models. The ${\mathit{\gamma }}_{\mathrm{ocean}}$ is not affected by the calibration and remains 45 % too low, which again suggests a structural limitation in our formulation of the ocean sink. This is, however, compensated for during calibration by the ${\mathit{\gamma }}_{\mathrm{land}}$ being 26 % higher than in complex models.

4.3 Scenarios

To further validate Pathfinder, we run the five SSP scenarios for which climate and carbon cycle projections were reported by a large enough number of models in the AR6 (namely, ssp119, ssp126, ssp245, ssp370, and ssp585). These simulations are run with prescribed ${\mathrm{CO}}_{\mathrm{2}}$ concentration and non-${\mathrm{CO}}_{\mathrm{2}}$ ERF (the latter is taken from ). Time series of GMST and cumulative land and ocean sinks are shown in Fig. . Table  shows a comparison of the projected changes in GMST to the CMIP6 estimates Table 4.2 and of carbon pools to (since this was not directly reported in the AR6).

Be it on short-, mid- or long-term scales, Pathfinder's projections of GMST are very much in line with the one assessed by the IPCC in the AR6 based on multiple lines of evidence Table 4.5. The only significant difference is a smaller uncertainty range in our projections for the longer-term periods. Although this is the result of the efficiency of the Bayesian calibration, one might wonder whether the climate module is over-constrained (or equivalently, too limited in its number of parameters).

The ocean carbon storage appears to be overestimated by 5 % to 20 % by Pathfinder across SSP scenarios. This is consistent with the upward adjustment of the ocean carbon sink stemming from our Bayesian calibration. To compare the land carbon storage with CMIP6 models, because our land carbon module does not include land use change processes, we correct the value reported by complex models by the cumulative land use change emissions of each scenario . While the land carbon storage of Pathfinder is well in line under ssp126 (a scenario consistent with the 2 ${}^{\circ }\mathrm{C}$ target), it is underestimated in ssp119 (consistent with the 1.5 ${}^{\circ }\mathrm{C}$ target) and increasingly overestimated in higher warming scenarios. A likely explanation is that the climate–carbon feedback on land is underestimated in Pathfinder, as suggested by the $\mathit{\gamma }$ metric seen in Sect. . Alternatively (or concurrently), the absence of loss of sink capacity caused by land cover change can explain the overestimation of the land sink under high ${\mathrm{CO}}_{\mathrm{2}}$. The Pathfinder model's estimates of both sinks nonetheless remain well within the CMIP6 models' uncertainty ranges.

Our SLR emulator gives estimates (Table ) that are always on the lower end of the range reported in the AR6 Table 9.9. This can be explained by the fact that our individual SLR rate estimates are on the lower end of their respective constraints (see Sect. ). This discrepancy also highlights potential structural limitations in the SLR module (e.g. too few separate contributions) and the difficulty of calibrating the module given the short time period of data available, both from complex models (that cover the 21st century only) and observations, compared to the long timescale of the SLR processes. Nevertheless, our estimates remain within the uncertainty range of the IPCC assessment, especially as we do not account for contribution from land water storage that causes an additional 0.03 [0.01, 0.04] $\mathrm{m}$ of SLR in all scenarios in 2100 .

Table 4

Comparison of SSP scenarios between Pathfinder and AR6 for SLR (with respect to 1995–2014) and SLR speed projections Table 9.9. We give the median value and the 90 % confidence interval in parentheses.

Experiment	Model	SLR (m)	SLR (m)	SLR rate (mm yr${}^{-\mathrm{1}}$)	SLR rate (mm yr${}^{-\mathrm{1}}$)
		2050	2100	2040–2060	2080–2100
ssp119	Pathfinder	0.15 (0.13, 0.18)	0.30 (0.25, 0.36)	3.5 (2.9, 4.2)	2.7 (2.2, 3.4)
ssp119	AR6	0.18 (0.15, 0.23)	0.38 (0.28, 0.55)	4.1 (2.8, 6.0)	4.2 (2.4, 6.6)
ssp126	Pathfinder	0.16 (0.14, 0.19)	0.35 (0.30, 0.43)	4.1 (3.5, 5.0)	3.6 (2.9, 4.5)
ssp126	AR6	0.19 (0.16, 0.25)	0.44 (0.32, 0.62)	4.8 (3.5, 6.8)	5.2 (3.2, 8.0)
ssp245	Pathfinder	0.17 (0.15, 0.20)	0.46 (0.39, 0.56)	5.0 (4.2, 6.0)	6.2 (5.1, 8.0)
ssp245	AR6	0.20 (0.17, 0.26)	0.56 (0.44, 0.76)	5.8 (4.4, 8.0)	7.7 (5.2, 11.6)
ssp370	Pathfinder	0.18 (0.15, 0.21)	0.56 (0.47, 0.69)	5.5 (4.7, 6.7)	9.1 (7.4, 11.7)
ssp370	AR6	0.22 (0.18, 0.27)	0.68 (0.55, 0.90)	6.4 (5.0, 8.7)	10.4 (7.4, 14.8)
ssp585	Pathfinder	0.19 (0.17, 0.23)	0.67 (0.56, 0.82)	6.4 (5.4, 7.8)	11.4 (9.1, 15.0)
ssp585	AR6	0.23 (0.20, 0.29)	0.77 (0.63, 1.01)	7.2 (5.6, 9.7)	12.1 (8.6, 17.6)

5 Concluding remarks

In this paper, we have presented the Pathfinder model: a simple global carbon–climate model with selected impact variables, carefully designed to balance accuracy of representation and simplicity of formulation and calibrated through Bayesian inference on the latest data from Earth system models and observations. Pathfinder has been shown to perform very well in comparison to complex models, although there remains room for further improvement of the model and its calibration setup. We identify four main avenues to improve the model.

First, some parts of the model may well lean too much on the complexity side of the simplicity–accuracy balance we aimed to strike, owing to the creation process of Pathfinder that mostly compiled existing formulations. Future development should therefore strive to reduce complexity wherever possible. The ocean carbon sub-pools and perhaps the land carbon pools are potential leads in this respect.

Second, the ocean carbon module alone appears to be limited by its structure, which has been inherited from a 25-year-old (yet seminal) article . Although it is undoubtedly a significant undertaking, developing an alternative formulation of the ocean carbon dynamic, calibrated on state-of-the-art ocean models and properly connected to ocean pH and the ocean of the climate module, would benefit more than just the SCM community.

Third, integration of land use and land cover change in such a model appears warranted, despite the difficulty of doing so in a physically sensible yet simple manner. Given our expertise with the OSCAR model and its bookkeeping module , we are confident that this can be done, although it will demand extra care to keep the model compatible with the IAMs it is also meant to be linked to.

Fourth, the Bayesian setup can be extended, notably by including more time periods for the existing constraints but also by introducing and constraining entirely new variables such as isotopic ratios or inter-hemispheric gradients ; however, a balance must be struck with respect to the calibration's computation time. Here, we caution against including complex models' results as constraints in the Bayesian calibration, as was done for the IPCC AR6 , as it goes against the philosophy of Pathfinder to use complex models' results as prior information only.

Fifth, although our model is restricted to ${\mathrm{CO}}_{\mathrm{2}}$ by design because of how IAMs like DICE are also limited to ${\mathrm{CO}}_{\mathrm{2}}$ emissions, we can imagine many reasons why one would want to add non-${\mathrm{CO}}_{\mathrm{2}}$ climate forcers into Pathfinder. We would suggest doing so by following the model's philosophy: that is, by taking existing reduced-complexity formulations such as something between FaIR and OSCAR and adding the new parameters into the Bayesian setup with the relevant observational constraints.

In spite of these few shortcomings and potential development leads, Pathfinder v1.0.1 is a powerful tool that fits the niche it has been created for perfectly. We will further demonstrate the strengths and flexibility of Pathfinder in other publications. Meanwhile, we invite the community to seize this open-source model and use it in any study that could benefit from a simple, fast, and accurate carbon–climate model aligned with the latest climate science.

Appendix A Additional information about the model

A1 Technical requirements

Pathfinder has been developed and run in Python (v3.7.6) (https://docs.python.org/3.7/, last access: 5 December 2022), preferentially using IPython (v7.19.0) . Currently, packages required to run it are NumPy (v1.19.2) , SciPy (v1.5.2) , and Xarray (v0.16.0) , and it has hard-coded dependencies on PyMC3 (v3.8) and Theano (v1.0.4) that are in fact used only for calibration. Other versions of Python or these packages were not tested.

The calibration procedure takes about 9 h to run on a desktop computer (with a base speed of 3.4 $\mathrm{GHz}$). Simple use of the model is much faster: the idealized experiments and SSP scenarios for this description paper, which represent 2984 simulated years, were run in about 20 min for all 2000 configurations and on a single core. A single simulated year takes a few tenths of a second, although a number of options in the model can drastically alter this performance. Note also that this scales sub-linearly with the amount of configurations or scenarios because of the internal workings of the Xarray package, albeit at the cost of increased demand in random-access memory.

A2 Known issues

Two relatively benign issues that have been identified during development remain unsolved. First, the model requires a high number of sub-time steps (i.e. high $n_t$) to remain stable under high ${\mathrm{CO}}_{\mathrm{2}}$ because of the ocean carbon cycle. Second, the version of the model that is driven by $T$ and $R_x$ time series is extremely sensitive to its inputs because it mathematically requires the first two derivatives of $T$ and the first derivative of $R_x$.

A3 Changelog

Brief description of the successive versions of Pathfinder. v1.0.1.

Exact same physical equations and numerical values as v1.0. Added best-guess parameters calculated as the average of the posterior distribution, and corresponding historical outputs, for single-configuration runs. Improved README and MANUAL files.

v1.0.

First release. Exact model described in the preprint version of this paper .

Appendix B Additional figures and tables

Figure B1

Calibration to estimate prior ${\mathit{\nu }}_{\mathrm{gx}}$ and ${\mathit{\gamma }}_{\mathrm{gx}}$ from CMIP6 time series of $F_{\mathrm{ocean}}$. We fit our equation on the results of the $+$1 % CO${}_{\mathrm{2}}$ (1pctCO2) experiment (in blue) and its variants 1pctCO2-rad (in green) and 1pctCO2-bgc (in orange). Coloured markers are CMIP6 model data, while the solid black lines show the fit from Pathfinder. Panels without a black line indicate that at least one of the required variables was not reported by the complex model.

[Figure omitted. See PDF]

Figure B2

Calibration to estimate prior ${\mathit{\beta }}_{\mathrm{dic}}$ and ${\mathit{\gamma }}_{\mathrm{dic}}$ from CMIP6 time series of $p_{{\mathrm{CO}}_{\mathrm{2}}}$. We fit our equation on the results of the $+$1 % CO${}_{\mathrm{2}}$ (1pctCO2) experiment (in blue) and its variants 1pctCO2-rad (in green) and 1pctCO2-bgc (in orange). Coloured markers are CMIP6 model data, while the solid black lines show the fit from Pathfinder. Panels without a black line indicate that at least one of the required variables was not reported by the complex model.

[Figure omitted. See PDF]

Figure B3

Calibration to estimate prior ${\mathit{\beta }}_{\mathrm{ef}}$ and ${\mathit{\gamma }}_{\mathrm{ef}}$ from CMIP6 time series of $r_{\mathrm{fire}}$, shown as a ratio over its preindustrial value. We fit our equation on the results of the $+$1 % CO${}_{\mathrm{2}}$ (1pctCO2) experiment (in blue) and its variants 1pctCO2-rad (in green) and 1pctCO2-bgc (in orange). Coloured markers are CMIP6 model data, while the solid black lines show the fit from Pathfinder. Panels without a black line indicate that at least one of the required variables was not reported by the complex model.

[Figure omitted. See PDF]

Figure B4

Calibration to estimate prior ${\mathit{\beta }}_{\mathrm{rh}}$ and ${\mathit{\gamma }}_{\mathrm{rh}}$ from CMIP6 time series of $r_{\mathrm{rh}}$, shown as a ratio over its preindustrial value. We fit our equation on the results of the $+$1 % CO${}_{\mathrm{2}}$ (1pctCO2) experiment (in blue) and its variants 1pctCO2-rad (in green) and 1pctCO2-bgc (in orange). Coloured markers are CMIP6 model data, while the solid black lines show the fit from Pathfinder. Panels without a black line indicate that at least one of the required variables was not reported by the complex model.

[Figure omitted. See PDF]

Figure B5

Calibration to estimate prior ${\mathit{\beta }}_{\mathrm{npp}}$, ${\mathit{\alpha }}_{\mathrm{npp}}$, and ${\mathit{\gamma }}_{\mathrm{npp}}$ from CMIP6 time series of $r_{\mathrm{npp}}$, shown as a ratio over its preindustrial value. We fit our equation on the results of the $+$1 % CO${}_{\mathrm{2}}$ (1pctCO2) experiment (in blue) and its variants 1pctCO2-rad (in green) and 1pctCO2-bgc (in orange). Coloured markers are CMIP6 model data, while the solid black lines show the fit from Pathfinder. Panels without a black line indicate that at least one of the required variables was not reported by the complex model.

[Figure omitted. See PDF]

Figure B6

Calibration to estimate the prior of GIS SLR module parameters (${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}$, ${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{3}}$, and ${\mathit{\lambda }}_{\mathrm{gis}\mathrm{0}}$). We fit our equation on the compiled outputs from for which there is more than one scenario available. Each panel's title displays the name of the institute, model, and configuration used. Coloured markers are the model data, while the solid black lines show the fit from Pathfinder.

[Figure omitted. See PDF]

Figure B7

Calibration to estimate the prior of glacier SLR module parameters (${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{1}}$, ${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{3}}$, ${\mathit{\gamma }}_{\mathrm{gla}}$, and ${\mathit{\lambda }}_{\mathrm{gla}\mathrm{0}}$). We fit our equation on the compiled outputs from for which there is more than one scenario available. Each panel's title displays the name of the model or study used. Coloured markers are the model data, while the solid black lines show the fit from Pathfinder.

[Figure omitted. See PDF]

Figure B8

Calibration to estimate the prior of AIS SLR module parameters (${\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}$, ${\mathit{\alpha }}_{\mathrm{ais}}$, and ${\mathit{\lambda }}_{\mathrm{ais}\mathrm{0}}$). We fit our equation on the compiled outputs from for which there is more than one scenario available. Each panel's title displays the name of the institute, model and configuration used. Coloured markers are the models data while the solid black lines show the fit from Pathfinder.

[Figure omitted. See PDF]

Figure B9

Distribution of the logit of the ECS ($T_{\mathrm{2}\times }$) inferred from AR6. Blue points are the assessments from AR6, the solid line is the cumulative distributive function (CDF) fitted on those assessments, and the dashed line is the PDF associated with the CDF (arbitrary scale). The value of the fitted parameters is given above the plot.

[Figure omitted. See PDF]

Table B1

Summary of Pathfinder's equation variables in climate, sea level, and ocean carbon modules.

Manual	Code	Description	Units
$R_{\mathrm{c}}$	RFco2	${\mathrm{CO}}_{\mathrm{2}}$ (effective) radiative forcing	W m${}^{-\mathrm{2}}$
$R_x$	ERFx	Non-${\mathrm{CO}}_{\mathrm{2}}$ effective radiative forcing	W m${}^{-\mathrm{2}}$
$R$	ERF	Effective radiative forcing	W m${}^{-\mathrm{2}}$
$T$	T	Global surface temperature anomaly	K
$T_{\mathrm{d}}$	Td	Deep-ocean temperature anomaly	K
$\mathrm{logit}\left(\mathrm{ff}\right)$	logit_ff	Logit of the climate feedback factor (for calib.)	1
$U_{\mathrm{ohc}}$	OHC	Ocean heat content (anomaly)	W yr m${}^{-\mathrm{2}}$
$H_{\mathrm{thx}}$	Hthx	Thermosteric sea level rise	mm
$H_{\mathrm{gla}}$	Hgla	Glaciers' contribution to sea level rise	mm
$H_{\mathrm{gis}}$	Hgis	Greenland ice sheet's contribution to sea level rise	mm
$H_{\mathrm{ais},\mathrm{smb}}$	Hais_smb	Surface mass balance component of Hais	mm
$H_{\mathrm{ais}}$	Hais	Antarctic ice sheet's contribution to sea level rise	mm
$H_{\mathrm{tot}}$	Htot	Total sea level rise	mm
$H_{\mathrm{lia}}$	Hlia	Sea level rise from relaxation after LIA between 1900 and 2005 (for calib.)	mm
$C_{\mathrm{o},j}$	Co_j	Change in surface ocean carbon sub-pools	PgC	$j\in \left[\right[\mathrm{1},\mathrm{5}\left]\right]$
$C_{\mathrm{o}}$	Co	Change in surface ocean carbon pool	PgC
$C_{\mathrm{d}}$	Cd	Change in deep-ocean carbon pool
$c_{\mathrm{dic}}$	dic	Change in surface DIC	$\mathrm{\mu }\mathrm{mol}$ kg${}^{-\mathrm{1}}$
$p_{\mathrm{dic}}$	pdic	Subcomponent of pCO2	ppm
$p_{{\mathrm{CO}}_{\mathrm{2}}}$	pCO2	${\mathrm{CO}}_{\mathrm{2}}$ partial pressure at the ocean surface	ppm
$F_{\mathrm{ocean}}$	Focean	Ocean carbon sink	PgC yr${}^{-\mathrm{1}}$

Table B2

Summary of Pathfinder's equation variables for land carbon, permafrost, and atmospheric modules.

Manual	Code	Description	Units
$r_{\mathrm{npp}}$	r_npp	Relative change in NPP	1
$r_{\mathrm{fire}}$	r_fire	Relative change in wildfire intensity	1
$r_{\mathrm{rh}}$	r_rh	Relative change in heterotrophic respiration rate	1
$F_{\mathrm{npp}}$	NPP	Net primary productivity	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{fire}}$	Efire	Emissions from wildfire	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{harv}}$	Eharv	Emissions from harvest and grazing	PgC yr${}^{-\mathrm{1}}$
$F_{\mathrm{mort}}$	Fmort	Mortality flux	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{rh}\mathrm{1}}$	RH1	Litter heterotrophic respiration	PgC yr${}^{-\mathrm{1}}$
$F_{\mathrm{stab}}$	Fstab	Stabilization flux	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{rh}\mathrm{2}}$	RH2	Active soil heterotrophic respiration	PgC yr${}^{-\mathrm{1}}$
$F_{\mathrm{pass}}$	Fpass	Passivization flux	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{rh}\mathrm{3}}$	RH3	Passive soil heterotrophic respiration	PgC yr${}^{-\mathrm{1}}$
$F_{\mathrm{land}}$	Fland	Land carbon sink	PgC yr${}^{-\mathrm{1}}$
$E_{\mathrm{rh}}$	RH	Heterotrophic respiration	PgC yr${}^{-\mathrm{1}}$
$C_{\mathrm{v}}$	Cv	Vegetation carbon pool	PgC
$C_{\mathrm{s}\mathrm{1}}$	Cs1	Litter carbon pool	PgC
$C_{\mathrm{s}\mathrm{2}}$	Cs2	Active soil carbon pool	PgC
$C_{\mathrm{s}\mathrm{3}}$	Cs3	Passive soil carbon pool	PgC
$C_{\mathrm{s}}$	Cs	Total soil carbon pool	PgC
$r_{\mathrm{rt}}$	r_rt	Relative change in permafrost respiration rate	1
$\overline{a}$	abar	Theoretical thawed fraction	1
$a$	a	Actual thawed fraction	1
$E_{\mathrm{pf}}$	Epf	Emissions from permafrost	PgC yr${}^{-\mathrm{1}}$
$C_{\mathrm{th},j}$	Cth_j	Thawed permafrost carbon sub-pools	PgC	$j\in \left[\right[\mathrm{1},\mathrm{3}\left]\right]$
$C_{\mathrm{fr}}$	Cfr	Frozen permafrost carbon pool	PgC
$E_{{\mathrm{CO}}_{\mathrm{2}}}$	Eco2	Anthropogenic ${\mathrm{CO}}_{\mathrm{2}}$ emissions	PgC yr${}^{-\mathrm{1}}$
$C$	${\mathrm{CO}}_{\mathrm{2}}$	Atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration	ppm
$\mathrm{pH}$	pH	Surface ocean pH	1

Table B3

Parameters used in the climate, ocean carbon, and sea level modules.

/ Manual	Code	Description	Units
$\mathit{\varphi }$	phi	Radiative parameter of ${\mathrm{CO}}_{\mathrm{2}}$	W m${}^{-\mathrm{2}}$
$T_{\mathrm{2}\times }$	T2x	Equilibrium climate sensitivity	K
${\mathrm{\Theta }}_{\mathrm{s}}$	THs	Heat capacity of the surface	W yr m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$
${\mathrm{\Theta }}_{\mathrm{d}}$	THd	Heat capacity of the deep ocean	W yr m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$
$\mathit{\theta }$	th	Heat exchange coefficient	W m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$
${\mathit{ϵ}}_{\mathrm{heat}}$	eheat	Deep-ocean heat uptake efficacy	1
$T_{\mathrm{2}\times }^{*}$	T2x0	Minimal value of the ECS distribution (for calib.)	K
${\mathit{\alpha }}_{\mathrm{ohc}}$	aOHC	Fraction of energy warming the ocean	1
${\mathrm{\Lambda }}_{\mathrm{thx}}$	Lthx	Proportionality factor of thermosteric SLR	mm m${}^{\mathrm{2}}$ W${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$
${\mathit{\lambda }}_{\mathrm{gla}}$	lgla0	Initial imbalance in SLR from glaciers	mm yr${}^{-\mathrm{1}}$
${\mathrm{\Lambda }}_{\mathrm{gla}}$	Lgla	Maximum contribution to SLR from glaciers	mm
${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{1}}$	Ggla1	Linear sensitivity of steady-state glaciers SLR to climate	K${}^{-\mathrm{1}}$
${\mathrm{\Gamma }}_{\mathrm{gla}\mathrm{3}}$	Ggla3	Cubic sensitivity of steady-state glaciers SLR to climate	K${}^{-\mathrm{3}}$
${\mathit{\tau }}_{\mathrm{gla}}$	tgla	Timescale of glaciers' contribution to SLR	year
${\mathit{\gamma }}_{\mathrm{gla}}$	ggla	Sensitivity of glaciers' timescale to climate	K${}^{-\mathrm{1}}$
${\mathit{\lambda }}_{\mathrm{gis}}$	lgis0	Initial imbalance in SLR from GIS	mm yr${}^{-\mathrm{1}}$
${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{1}}$	Lgis1	Linear sensitivity of steady-state GIS SLR to climate	mm K${}^{-\mathrm{1}}$
${\mathrm{\Lambda }}_{\mathrm{gis}\mathrm{3}}$	Lgis3	Cubic sensitivity of steady-state GIS SLR to climate	mm K${}^{-\mathrm{3}}$
${\mathit{\tau }}_{\mathrm{gis}}$	tgis	Timescale of GIS contribution to SLR	year
${\mathrm{\Lambda }}_{\mathrm{ais},\mathrm{smb}}$	Lais_smb	Sensitivity of AIS SMB increase due to climate	mm yr${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$
${\mathit{\lambda }}_{\mathrm{ais}}$	lais	Initial imbalance in SLR from AIS	mm yr${}^{-\mathrm{1}}$
${\mathrm{\Lambda }}_{\mathrm{ais}}$	Lais	Sensitivity of steady-state AIS SLR to climate	mm K${}^{-\mathrm{1}}$
${\mathit{\tau }}_{\mathrm{ais}}$	tais	Timescale of AIS contribution to SLR	year
${\mathit{\alpha }}_{\mathrm{ais}}$	aais	Sensitivity of AIS timescale to AIS SLR	mm${}^{-\mathrm{1}}$
${\mathit{\alpha }}_{\mathrm{dic}}$	adic	Conversion factor for DIC	µmol kg${}^{-\mathrm{1}}$ PgC${}^{-\mathrm{1}}$
${\mathit{\beta }}_{\mathrm{dic}}$	bdic	Inverse scaling factor for DIC	1
${\mathit{\gamma }}_{\mathrm{dic}}$	gdic	Sensitivity of pCO2 to climate	K${}^{-\mathrm{1}}$
$T_{\mathrm{o}}$	To	Preindustrial surface ocean temperature	°C
${\mathit{\nu }}_{\mathrm{gx}}$	vgx	Surface ocean gas exchange rate	yr${}^{-\mathrm{1}}$
${\mathit{\gamma }}_{\mathrm{gx}}$	ggx	Sensitivity of gas exchange to climate	K${}^{-\mathrm{1}}$
${\mathit{\alpha }}_{\mathrm{o},j}$	aoc_j	Surface ocean sub-pool fractions	1 $j\in \left[\right[\mathrm{1},\mathrm{5}\left]\right]$
${\mathit{\tau }}_{\mathrm{o},j}$	toc_j	Timescales of surface ocean sub-pools	year $j\in \left[\right[\mathrm{1},\mathrm{5}\left]\right]$
${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{o}}}$	k_toc	Scaling factor for timescales of surface ocean sub-pools	1

Table B4

Parameters used in the permafrost, land carbon modules, and for calibration.

Manual	Code	Description	Units
${\mathit{\beta }}_{\mathrm{npp}}$	bnpp	Sensitivity of NPP to ${\mathrm{CO}}_{\mathrm{2}}$ (fertilization effect)	1
${\mathit{\alpha }}_{\mathrm{npp}}$	anpp	Shape parameter for fertilization effect	1
${\mathit{\gamma }}_{\mathrm{npp}}$	gnpp	Sensitivity of NPP to climate	K${}^{-\mathrm{1}}$
${\mathit{\beta }}_{\mathrm{fire}}$	bfire	Sensitivity of wildfire intensity to ${\mathrm{CO}}_{\mathrm{2}}$	1
${\mathit{\gamma }}_{\mathrm{fire}}$	gfire	Sensitivity of wildfire intensity to climate	K${}^{-\mathrm{1}}$
${\mathit{\beta }}_{\mathrm{rh}}$	brh	Sensitivity of heterotrophic respiration to fresh organic matter	1
${\mathit{\gamma }}_{\mathrm{rh}}$	grh	Sensitivity of heterotrophic respiration to climate	K${}^{-\mathrm{1}}$
$F_{\mathrm{npp},\mathrm{0}}$	npp0	Preindustrial NPP	PgC yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{fire}}$	vfire	Wildfire intensity	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{harv}}$	vharv	Harvest and grazing rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{mort}}$	vmort	Mortality rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{stab}}$	vstab	Stabilization rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{rh}\mathrm{1}}$	vrh1	Litter heterotrophic respiration rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{rh}\mathrm{23}}$	vrh23	Soil (active and passive) respiration rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{rh}\mathrm{3}}$	vrh3	Passive soil respiration rate	yr${}^{-\mathrm{1}}$
${\mathit{\alpha }}_{\mathrm{pass}}$	apass	Fraction of passive soil	1
${\mathit{\alpha }}_{\mathrm{lst}}$	aLST	Climate scaling factor over permafrost regions	1
${\mathit{\gamma }}_{\mathrm{rt}\mathrm{1}}$	grt1	Sensitivity of (boreal) heterotrophic respiration to climate	K${}^{-\mathrm{1}}$
${\mathit{\gamma }}_{\mathrm{rt}\mathrm{2}}$	grt2	Sensitivity of (boreal) heterotrophic respiration to climate (quadratic)	K${}^{-\mathrm{2}}$
${\mathit{\kappa }}_{\mathrm{rt}}$	krt	Scaling factor for sensitivity of permafrost respiration to climate	1
$a_{\min }$	amin	Minimal thawed fraction	1
${\mathit{\kappa }}_a$	ka	Shape parameter for theoretical thawed fraction	1
${\mathit{\gamma }}_a$	ga	Sensitivity of theoretical thawed fraction to climate	K${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{thaw}}$	vthaw	Thawing rate	yr${}^{-\mathrm{1}}$
${\mathit{\nu }}_{\mathrm{froz}}$	vfroz	Freezing rate	yr${}^{-\mathrm{1}}$
${\mathit{\alpha }}_{\mathrm{th},j}$	ath_j	Thawed permafrost carbon sub-pool fractions	1 $j\in \left[\right[\mathrm{1},\mathrm{3}\left]\right]$
${\mathit{\tau }}_{\mathrm{th},j}$	tth_j	Timescales of thawed permafrost carbon sub-pools	year $j\in \left[\right[\mathrm{1},\mathrm{3}\left]\right]$
${\mathit{\kappa }}_{{\mathit{\tau }}_{\mathrm{th}}}$	k_tth	Scaling factor for timescales of surface ocean sub-pools	1
$C_{\mathrm{fr},\mathrm{0}}$	Cfr0	Preindustrial frozen permafrost carbon pool	PgC
${\mathit{\alpha }}_{\mathrm{C}}$	aCO2	Conversion factor for atmospheric ${\mathrm{CO}}_{\mathrm{2}}$	PgC ppm${}^{-\mathrm{1}}$
$C_{\mathrm{pi}}$	${\mathrm{CO}}_{\mathrm{2}}$pi	Preindustrial ${\mathrm{CO}}_{\mathrm{2}}$ concentration	ppm
${\mathit{\kappa }}_{\mathrm{pH}}$	k_pH	Scaling factor for surface ocean pH	1
${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{C}}$	std_${\mathrm{CO}}_{\mathrm{2}}$	Relative standard deviation of the historical ${\mathrm{CO}}_{\mathrm{2}}$ time series (for calib.)	1
${\mathit{ϵ}}_{\mathrm{C}}$	ampl_${\mathrm{CO}}_{\mathrm{2}}$	Noise amplitude of the historical ${\mathrm{CO}}_{\mathrm{2}}$ time series (for calib.)	ppm
${\mathit{\rho }}_{\mathrm{C}}$	corr_${\mathrm{CO}}_{\mathrm{2}}$	Autocorrelation of the historical ${\mathrm{CO}}_{\mathrm{2}}$ time series (for calib.)	1
${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_T$	std_T	Relative standard deviation of the historical $T$ time series (for calib.)	1
${\mathit{ϵ}}_T$	ampl_T	Noise amplitude of the historical $T$ time series (for calib.)	K
${\mathit{\rho }}_T$	corr_T	Autocorrelation of the historical $T$ time series (for calib.)	1

Table B5

Structural parameter values. Parameters are noted under their code notation, and Tables  and provide the corresponding notation in text.

Parameters	Prior	Unit
Lgla	380	mm
Lais	1200	mm
tgla	190	year
tgis	481	year
tais	2090	year
T2x0	0.61	K
adic	4.49	$\mathrm{\mu }\mathrm{mol}$ kg${}^{-\mathrm{1}}$ PgC${}^{-\mathrm{1}}$
aoc_1	0.87	1
aoc_2	0.06	1
aoc_3	0.04	1
aoc_4	0.02	1
aoc_5	0.01	1
toc_1	1.29	year
toc_2	16.7	year
toc_3	65.1	year
toc_4	348	year
toc_5	10${}^{\mathrm{9}}$	year
vrh3	$\mathrm{8.27}\times {\mathrm{10}}^{-\mathrm{5}}$	yr${}^{-\mathrm{1}}$
aLST	1.87	1
grt1	0.12	K${}^{-\mathrm{1}}$
grt2	0.0029	K${}^{-\mathrm{2}}$
krt	1.34	1
amin	0.98	1
vthaw	0.14	yr${}^{-\mathrm{1}}$
vfroz	0.011	yr${}^{-\mathrm{1}}$
ath_1	0.05	1
ath_2	0.12	1
ath_3	0.83	1
tth_1	18.2	year
tth_2	252	year
tth_3	3490	year
aCO2	2.124	PgC ppm${}^{-\mathrm{1}}$
k_pH	1	1

Table B6

Calibrated parameter values before and after Bayesian calibration. Parameters are noted under their code notation, and Tables  and provide the corresponding notation in text. The uncertainty corresponds to the 1$\mathit{\sigma }$ uncertainty range.

Parameters	Prior	Posterior	Unit	Parameters	Prior	Posterior	Unit
phi	5.35 $\pm$ 0.54	5.29 $\pm$ 0.54	W m${}^{-\mathrm{2}}$	T2x	4.13 $\pm$ 1.37	3.37 $\pm$ 0.77	K
THs	8.14 $\pm$ 0.99	8.21 $\pm$ 1.06	W yr m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$	THd	108.6 $\pm$ 61.8	123.8 $\pm$ 57.8	W yr m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$
th	0.61 $\pm$ 0.13	0.67 $\pm$ 0.12	W m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$	eheat	1.35 $\pm$ 0.40	1.41 $\pm$ 0.43	1
aOHC	0.91 $\pm$ 0.02	0.91 $\pm$ 0.02	1	Lthx	1.82 $\pm$ 0.21	1.85 $\pm$ 0.23	mm m${}^{\mathrm{2}}$ W${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$
lgla0	0.59 $\pm$ 0.24	0.40 $\pm$ 0.21	mm yr${}^{-\mathrm{1}}$	Ggla1	0.34 $\pm$ 0.18	0.34 $\pm$ 0.05	mm K${}^{-\mathrm{1}}$
Ggla3	0.022 $\pm$ 0.013	0.022 $\pm$ 0.013	mm K${}^{-\mathrm{3}}$	ggla	0.12 $\pm$ 0.09	0.11 $\pm$ 0.07	K${}^{-\mathrm{1}}$
Lgis1	82 $\pm$ 45	189 $\pm$ 55	mm K${}^{-\mathrm{1}}$	Lgis3	5.7 $\pm$ 1.4	5.8 $\pm$ 1.5	mm K${}^{-\mathrm{3}}$
Lais_smb	0.61 $\pm$ 0.19	0.40 $\pm$ 0.10	mm K${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$	lais0	0.00 $\pm$ 0.11	0.07 $\pm$ 0.08	mm yr${}^{-\mathrm{1}}$
aais	0.002 $\pm$ 0.003	0.004 $\pm$ 0.003	mm${}^{-\mathrm{1}}$	lgis0	0.33 $\pm$ 0.14	0.35 $\pm$ 0.14	mm yr${}^{-\mathrm{1}}$
k_toc	1.00 $\pm$ 0.20	0.91 $\pm$ 0.18	1	vgx	0.19 $\pm$ 0.06	0.20 $\pm$ 0.07	PgC ppm${}^{-\mathrm{1}}$ yr${}^{-\mathrm{1}}$
ggx	0.018 $\pm$ 0.029	0.019 $\pm$ 0.033	K${}^{-\mathrm{1}}$	To	18.0 $\pm$ 0.5	18.0 $\pm$ 0.5	K
bdic	0.87 $\pm$ 0.08	0.90 $\pm$ 0.09	1	gdic	0.04 $\pm$ 0.02	0.04 $\pm$ 0.02	K${}^{-\mathrm{1}}$
npp0	48.2 $\pm$ 5.1	46.5 $\pm$ 3.3	PgC yr${}^{-\mathrm{1}}$	vfire	0.006 $\pm$ 0.003	0.006 $\pm$ 0.002	yr${}^{-\mathrm{1}}$
vharv	0.003 $\pm$ 0.003	0.003 $\pm$ 0.002	yr${}^{-\mathrm{1}}$	vmort	0.11 $\pm$ 0.01	0.11 $\pm$ 0.01	yr${}^{-\mathrm{1}}$
vstab	0.32 $\pm$ 0.28	0.30 $\pm$ 0.22	yr${}^{-\mathrm{1}}$	vrh1	0.33 $\pm$ 0.29	0.27 $\pm$ 0.20	yr${}^{-\mathrm{1}}$
vrh23	0.024 $\pm$ 0.009	0.024 $\pm$ 0.008	yr${}^{-\mathrm{1}}$	bnpp	0.93 $\pm$ 0.37	1.09 $\pm$ 0.25	1
anpp	0.48 $\pm$ 0.57	0.36 $\pm$ 0.38	1	gnpp	$-$0.014 $\pm$ 0.023	$-$0.005 $\pm$ 0.024	K${}^{-\mathrm{1}}$
bfire	$-$0.05 $\pm$ 0.12	$-$0.06 $\pm$ 0.14	1	gfire	0.052 $\pm$ 0.072	0.044 $\pm$ 0.088	K${}^{-\mathrm{1}}$
brh	1.06 $\pm$ 0.43	1.01 $\pm$ 0.41	1	grh	0.056 $\pm$ 0.053	0.042 $\pm$ 0.035	K${}^{-\mathrm{1}}$
apass	0.69 $\pm$ 0.19	0.63 $\pm$ 0.20	1	ga	0.13 $\pm$ 0.04	0.13 $\pm$ 0.04	K${}^{-\mathrm{1}}$
ka	2.6 $\pm$ 2.0	2.4 $\pm$ 1.8	1	k_tth	0.96 $\pm$ 0.93	1.00 $\pm$ 0.87	1
Cfr0	546 $\pm$ 120	538 $\pm$ 122	PgC	CO2pi	278 $\pm$ 3	279 $\pm$ 3	ppm

Code and data availability

The source code of Pathfinder is openly available at https://github.com/tgasser/Pathfinder (last access: 5 December 2022). A frozen version of the code and data as developed in the paper can be found on Zenodo at 10.5281/zenodo.7003848 .

Author contributions

TG developed the model, with contributions from TB regarding sea level rise and PC regarding land carbon cycle. TG coded the model and its calibration. TB ran the diagnostic simulations and made the final figures. TB and TG wrote the manuscript with input from PC.

Competing interests

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

Disclaimer

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

Acknowledgements

We thank Andrew H. MacDougall for sharing permafrost simulations made with UVic, and Côme Cheritel for the first version of Figs. 1–4.

Financial support

Development of Pathfinder was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement no. 820829 (CONSTRAIN project) and by the Austrian Science Fund (FWF) under grant agreement no. P31796-N29 (ERM project).

Review statement

This paper was edited by Marko Scholze and reviewed by Ian Enting and one anonymous referee.