A fresh look at the pre-industrial air–sea carbon flux using the alkalinity budget

Disparities between observational and model-based estimates of the ocean carbon sink persist, highlighting the need for improved understanding and methodologies to reconcile differences in both magnitude and trends over recent decades. A potential key source of uncertainty lies in the pre-industrial air–sea carbon flux, which is essential for isolating the anthropogenic component from observations. This flux, thought to result globally from an imbalance between riverine discharge and sediment burial of carbon, remains highly uncertain, limiting the confidence in impactful applications such as the Global Carbon Budget (GCB). In this study, we present a new theoretical framework that enables direct estimation of the riverine/burial-driven pre-industrial carbon outgassing using both carbon and alkalinity budgets. This approach is validated with a series of ocean biogeochemical simulations, which also highlight the main factors influencing its regional distribution. We then demonstrate the utility of the framework through two proof-of-concept applications. The first revisits the pre-industrial riverine/burial-driven air–sea carbon flux using existing carbon and alkalinity budgets, offering a simple method for reassessment as these budgets are updated. The second application leverages sensitivity simulations to construct a composite simulated estimate that aligns with both carbon and alkalinity budgets to assess the regional distribution of the pre-industrial riverine/burial-driven air–sea carbon flux. This approach is well suited for model intercomparisons, enabling an efficient reassessment of regional flux patterns and helping to reduce biases related to ocean model physics or biogeochemical parameterizations.

1 Introduction

Accurately estimating the anthropogenic carbon sink in the ocean is crucial for gaining a deeper understanding of the underlying mechanisms and is a prerequisite for projecting its future evolution and the climate response to future emissions scenarios (Canadell et al., 2021). This anthropogenic carbon flux is currently assessed with yearly updates by the Global Carbon Budget (GCB; Friedlingstein et al., 2025), using both observational products and model simulations employing Global Ocean Biogeochemical Models (GOBMs). However, reconstructions derived from surface ocean pCO₂ data – currently the main observation-based approach – tend to yield higher estimates than models, both globally and regionally. This mismatch has grown steadily since the early 2000s, now reaching a 10 %–20 % difference over the past decade (e.g. Hauck et al., 2020; DeVries et al., 2023; Gruber et al., 2023; Friedlingstein et al., 2025).

Recent studies have begun to investigate the origins of discrepancies in both the magnitude and trend of observational versus model-based estimates of the ocean anthropogenic carbon sink. Analyses of GOBM-derived estimates (Terhaar et al., 2024) and pCO₂-based products (Ford et al., 2024) point to multiple sources of uncertainty, including methodological differences, model biases, and data limitations. On the modelling side, GOBMs have been shown to underestimate the global sink magnitude (Terhaar et al., 2022), as well as decadal variability, especially in the Southern Ocean (Mayot et al., 2023, 2024). On the observational side, the sparse and uneven spatial and temporal coverage of surface ocean pCO₂ measurements remains a major limitation (Hauck et al., 2023; Landschützer et al., 2023; Dong et al., 2024). While the causes of these mismatches are likely multifaceted, one less-discussed contributor is the uncertainty surrounding the pre-industrial air–sea carbon flux and its influence on pCO₂-based estimates.

The net air–sea carbon flux derived from pCO₂-based data encompasses both anthropogenic and natural components. The natural component originates, at the global scale, from the balance between riverine discharge and the burial of organic matter (OM) and calcium carbonate (CaCO₃). As these external fluxes together represent a net source of carbon for the ocean, they result in a net carbon outgassing at steady state during the pre-industrial era. Consequently, assessing the anthropogenic carbon flux through observations requires determining the pre-industrial riverine/burial-driven air–sea carbon flux and its spatial distribution (e.g. Hauck et al., 2020; Friedlingstein et al., 2025).

Assessing this outgassing carbon flux remains highly uncertain, with estimates ranging from 0.23 to 0.78 PgC yr⁻¹ (Aumont et al., 2001; Jacobson et al., 2007; Resplandy et al., 2018; Lacroix et al., 2020; Regnier et al., 2022), depending on the modelling approach used to derive them (forward or inverse) and estimates of riverine and burial fluxes (see Table E1). Specifically, the most recent estimate of 0.65±0.30 PgC yr⁻¹ is that used in the latest GCB release (Regnier et al., 2022; Friedlingstein et al., 2025, Table 1).

The spatial distribution of this riverine/burial-driven air–sea carbon flux is also highly uncertain. It strongly depends on the assumptions and methods used to assess it, including how sediment burial processes are represented, and both the magnitude and characteristics of riverine carbon inputs – particularly the balance between organic and inorganic forms, as well as the lability of terrestrial organic matter – (see Table E1). The earlier estimate, derived from a modelling analysis (Aumont et al., 2001), distributed this flux as follows: 49 % in the southern region, 25 % in the inter-tropical region, and 26 % in the northern region. In contrast, the most recent estimate, currently used in the GCB and also based on a modelling study (Lacroix et al., 2020, Table 1), suggests a very different partitioning, 14 %, 64 %, and 22 %, respectively, reshaping our understanding of the regional distribution of this flux.

Uncertainties in estimating the riverine/burial-driven pre-industrial outgassing may contribute to the persistent discrepancies between observation-based and model-derived estimates of the anthropogenic ocean carbon sink, both globally and regionally (Friedlingstein et al., 2025, their Figs. 11 and 14). These disparities have fluctuated over time, largely in response to stepwise adjustments made by the GCB team following new reassessments of the magnitude and spatial distribution of this flux in the literature (Fig. 1, Table 1). For instance, a substantial decrease in the global estimate of the pre-industrial outgassing from 2019 to 2020 contributed to a notable narrowing of the global observation-model gap. More recently, from 2022 to 2023, a redistribution of the flux between regions, from the southern region to the tropics, led to a reduced Southern Hemisphere bias and a compensating increase in the inter-tropical discrepancy.

Table 1 Review of the air–sea carbon outgassing from riverine/burial fluxes of carbon and Alk used in the GCBs. Both the global values and their regional distribution are presented, along with the associated references.

Figure 1 Schematic illustration of how a bias in evaluating the pre-industrial ocean carbon outgassing affects the assessment of the anthropogenic carbon flux based on pCO₂-products. A downward revision of the pre-industrial outgassing would decrease anthropogenic carbon flux estimates based on pCO₂-products, while an upward revision would increase it. This effect applies both globally and regionally.

Enhancing our understanding of the riverine/burial-driven air–sea carbon flux is critical to achieving more precise estimates of the anthropogenic carbon flux and its distribution from data-driven assessments. Numerical models hold great promise in addressing this challenge, particularly for estimating the spatial distribution of the flux. However, at present, the representation of the pre-industrial air–sea carbon flux remains uncertain in inter-model comparison exercises like CMIP6 (the 6th phase of the Coupled Model Intercomparison Project,  Eyring et al., 2016) and the 2024 GCB (Friedlingstein et al., 2025), likely due to differences in model setups and various/incomplete representations of sediment burial and riverine discharge (Terhaar et al., 2024). The magnitude of this global net flux ranges from −0.73 to 0.38 PgC yr⁻¹, while its inter-hemispheric gradient, defined as the difference between its values in the Northern Hemisphere and Southern Hemisphere, ranges from −0.09 to 0.82 PgC yr⁻¹ (Fig. 2a).

The methods employed thus far to estimate the riverine/burial-driven air–sea carbon flux at the global scale mostly rely on closing the ocean carbon budget. However, they often exhibit limitations in addressing the ocean alkalinity budget. Alkalinity (Alk), defined as the excess of proton acceptors over proton donors, or of positive conservative charges over negative ones, plays a pivotal role in driving air–sea carbon exchanges, which are strongly dependent on the relative balance between Alk and dissolved inorganic carbon (DIC; e.g. Humphreys et al., 2018). Similar to carbon, the Alk budget is controlled by both sources and sinks at the boundaries of the oceanic domain (Middelburg et al., 2020). Conventionally, it is hypothesized that the global Alk inventory has been conserved during the pre-industrial era, with the burial of CaCO₃ balancing the Alk riverine discharge (e.g. Revelle and Suess, 1957; Aumont et al., 2001; Planchat et al., 2023).

Nonetheless, the hypothesis of an imbalanced Alk budget during the pre-industrial era remains plausible, based on estimates of global Alk sources and sinks (e.g. Milliman, 1993; Middelburg et al., 2020). Paleoclimatic studies suggest that such an imbalance could arise from additional CaCO₃ burial (e.g. Cartapanis et al., 2018) or from a carbonate compensation mechanism involving biological processes alongside riverine inputs (Boudreau et al., 2018). Unlike carbon, Alk is not exchanged with the atmosphere, so balancing its budget depends on processes acting over longer timescales, particularly through interactions with the continents (e.g. erosion) and marine biogeochemistry (e.g. sediment dynamics). As a result, Alk budget balancing is slow, yet interactive with the carbon cycle through changes in atmospheric CO₂ and ocean acidity (Hain et al., 2014). An imbalance in the Alk budget would induce an additional air–sea carbon flux beyond that directly inferred from the ocean carbon budget, ultimately resulting in a non-conserved global ocean carbon inventory.

Here, we take a fresh look at the pre-industrial air–sea carbon flux by introducing a new theoretical framework that enables direct estimation of the riverine/burial-driven pre-industrial carbon outgassing, based on both carbon and Alk budgets. We validate this approach using a suite of ocean biogeochemical simulations, which also help identify the key drivers of its regional distribution. We then demonstrate the utility of this framework through two proof-of-concept applications: (i) revisiting the global magnitude of the pre-industrial riverine/burial-driven air–sea carbon flux using existing carbon and Alk budgets and (ii) reassessing its regional distribution using sensitivity simulations to construct a composite simulated estimate consistent with both budgets.

Figure 2 Pre-industrial air–sea carbon flux from models. The globally integrated pre-industrial flux (x axis) and its inter-hemispheric gradient (y axis) are indicated for (a) CMIP6 Earth system models (ESMs) and GCB Global Ocean Biogeochemical Models (GOBMs), as well as for (b) the NEMO-PISCES sensitivity simulations. (a) The 15 CMIP6 ESMs (10 GCB GOBMs; see Appendix B1) are plotted with red squares (orange circles). The black square and circle refer to the CMIP6 and GCB ensemble means. The CMIP6 and GCB ensemble ranges (line), mean (major tick), and quartiles (minor ticks) are respectively displayed to the top and right in red and orange. The star refers to the reference value used in the GCB 2024 (Table 1; Friedlingstein et al., 2025). The inter-hemispheric air–sea carbon flux gradient is defined as the difference between its values in the Northern Hemisphere and Southern Hemisphere (Sect. 2.1.4).

2 Methods

In this study, we use “steady state” to refer to the temporal stability of the globally integrated air–sea carbon flux. We describe the carbon and Alk budgets as “balanced” or “imbalanced” according to whether fluxes into and out of the ocean are quantitatively balanced. An imbalanced budget drives a deviation in the global inventory: a positive (negative) imbalance leads to an increase (decrease) in the global inventory. We refer to a “deviation” in a given variable when the system is in a steady state, but a persistent trend is identified for that variable (e.g. the global carbon inventory). In contrast, we use the term “drift” when a trend in a variable disrupts the steady state. Unless otherwise stated, all simulations and results refer to pre-industrial conditions.

Figure 3 Schematics of the theoretical framework introduced in this paper. (a) Conceptual representation of a process X (e.g. CaCO₃ burial), which affects carbon (F^C, X) and Alk (F^Alk, X), thereby induces, at steady state, a compensating carbon flux (F^{C, air–sea}(X)) and the resulting deviations in carbon and Alk inventories ($D^{\mathrm{C}/\mathrm{Alk}}\left(X\right)$). The equilibrium of the Alk:DIC pair with atmospheric CO₂ is represented by a solid grey line. Carbon (Alk) fluxes are represented through solid (dashed) arrows. (b) Schematic diagram of carbon and Alk budgets by ocean region. $F^{\mathrm{C}/\mathrm{Alk}}$ and $T^{\mathrm{C}/\mathrm{Alk}}$ respectively refer to the total external fluxes (directed into the ocean) and to the northward transport of carbon and Alk. $D^{\mathrm{C}/\mathrm{Alk}}$ corresponds to the regional carbon and Alk inventory deviations in each basin. S, I, and N refer to the different ocean regions, i.e. the Southern Hemisphere, the inter-tropical zone, and the Northern Hemisphere.

2.1 Theoretical framework

2.1.1 Governing equation of the pre-industrial riverine/burial-driven air–sea carbon flux

The collection of surface ocean pCO₂ data, and associated statistical methods, only allows for the direct reconstruction of the contemporary air–sea carbon flux ($F_{\mathrm{cont}.}^{\text{C, air–sea}}$), which encompasses both an anthropogenic ($F_{\mathrm{ant}.}^{\text{C, air–sea}}$) and a natural ($F_{\mathrm{nat}.}^{\text{C, air–sea}}$) component (e.g. Hauck et al., 2020), as follows:

$$\begin{array}{}\text{(1)}& F_{\mathrm{cont}.}^{\text{C, air–sea}}=F_{\mathrm{ant}.}^{\text{C, air–sea}}+F_{\mathrm{nat}.}^{\text{C, air–sea}},\end{array}$$

where positive fluxes are directed into the ocean (consistent throughout this paper). Therefore, it is crucial to determine the natural component to extract the anthropogenic carbon flux from pCO₂-based products. Within the anthropogenic carbon flux, we incorporate the perturbation of the natural carbon flux in response to climate change, ensuring that $F_{\mathrm{ant}.}^{\text{C, air–sea}}$ fully reflects the carbon sink resulting from all human-induced disturbances (e.g. Hauck et al., 2020). Accordingly, $F_{\mathrm{nat}.}^{\text{C, air–sea}}$ is directly defined as the riverine/burial-driven pre-industrial air–sea carbon flux ($F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}$), i.e.

$$\begin{array}{}\text{(2)}& F_{\mathrm{nat}.}^{\text{C, air–sea}}=F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}.\end{array}$$

The anthropogenic carbon flux can then be derived from pCO₂-based data as follows:

$$\begin{array}{}\text{(3)}& F_{\mathrm{ant}.}^{\text{C, air–sea}}=F_{\mathrm{cont}.}^{\text{C, air–sea}}-F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}.\end{array}$$

Assuming a steady-state pre-industrial air–sea carbon flux and a balanced Alk budget, the global riverine/burial-driven air–sea carbon flux can be directly inferred by closing the ocean carbon budget (e.g. Regnier et al., 2022):

$$\begin{array}{}\text{(4)}& F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}+F^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}=\mathrm{0}\end{array}$$

with

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{F}^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}& =F^{\mathrm{C},\mathrm{riv}.}+F^{\mathrm{C},\mathrm{bur}.\mathrm{org}.}+F^{\mathrm{C},\mathrm{bur}.\mathrm{inorg}.}\\ & +F^{\mathrm{C},\mathrm{minor}\mathrm{components}},\end{array}\end{array}$$

where “riv.” stands for “riverine discharge”, “bur. org.” for “OM burial”, and “bur. inorg.” for “CaCO₃ burial”, and “minor components” encompass other minor external fluxes, such as carbon release by mid-ocean ridges and groundwater discharge. Thus, assuming a conserved pre-industrial global Alk inventory, the riverine/burial-driven air–sea carbon flux is the opposite of the riverine and burial fluxes of carbon:

$$\begin{array}{}\text{(6)}& F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}=-F^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}.\end{array}$$

However, when considering a pre-industrial imbalanced Alk budget (i.e. a non-conserved global ocean Alk inventory), it becomes necessary to account for the Alk budget explicitly to infer the riverine/burial-driven air–sea carbon flux:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{F}^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}& =F^{\mathrm{Alk},\mathrm{riv}.}+F^{\mathrm{Alk},\mathrm{bur}.\mathrm{org}.}\\ & +F^{\mathrm{Alk},\mathrm{bur}.\mathrm{inorg}.}+F^{\mathrm{Alk},\mathrm{minor}\mathrm{components}},\end{array}\end{array}$$

where “minor components” this time encompass other minor external fluxes such as anaerobic processes, silicate weathering, and groundwater discharge. Importantly, under the assumption of a steady-state system – that is, with a stable air–sea carbon flux – any imbalance in the Alk budget induces a compensating carbon flux. To estimate this flux, we extend the conceptual framework introduced by Humphreys et al. (2018) by utilizing a phase diagram (Alk, DIC) in the form of an Alk and DIC flux diagram, while operating on a global scale (Fig. 3a). For any flux affecting carbon and/or Alk, it is possible to derive an air–sea carbon flux and the associated global carbon and Alk inventory imbalances. This approach relies on the equilibrium relationship between the Alk:DIC pair and atmospheric CO₂. At global scale, at steady state, any deviation in Alk is directly proportional to a DIC anomaly, and this proportionality coefficient can be estimated with high precision, as follows:

$$\begin{array}{}\text{(8)}& Q_{\mathrm{inv}}\simeq {\displaystyle \frac{\mathrm{Alk}}{\mathrm{3}\cdot \mathrm{Alk}-\mathrm{2}\cdot \mathrm{DIC}}},\end{array}$$

where Q_inv, as previously defined in Planchat et al. (2023), represents the inverse of the “isocapnic quotient” approximation introduced by Humphreys et al. (2018) (see Appendix A). In this study, Q_inv is defined based on the mean surface values of Alk and DIC. In the case of a steady-state air–sea carbon flux (see Fig. B2), every external process X (e.g. riverine discharge) that exerts an impact on carbon (F^C, X) and/or Alk (F^Alk, X) results in a global imbalance, shifting the surface ocean away from equilibrium with the atmosphere. Specifically, this requires an air–sea carbon flux (F^{C, air–sea}(X); Fig. 3a) to maintain global equilibrium with respect to the atmospheric CO₂. This also leads to deviations in global carbon and Alk inventories (D^C and D^Alk, respectively). In summary, for any given process X, we can define

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{F}^{\text{C, air–sea}}\left(X\right)=F^{\mathrm{Alk},X}\cdot Q_{\mathrm{inv}}-F^{\mathrm{C},X}\\ D^{\mathrm{C}}\left(X\right)=F^{\mathrm{Alk},X}\cdot Q_{\mathrm{inv}}\\ D^{\mathrm{Alk}}\left(X\right)=F^{\mathrm{Alk},X}.\end{array}\right.\end{array}$$

Applying this theoretical framework to the total external carbon and Alk fluxes ($F^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}$ and $F^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}$, respectively), we can deduce the riverine/burial-driven air–sea carbon flux ($F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}$) and the respective deviations in global carbon and Alk inventories as follows:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}{F}_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}=F^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}\cdot Q_{\mathrm{inv}}-F^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}\\ D^{\mathrm{C}}=F^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}\cdot Q_{\mathrm{inv}}\\ D^{\mathrm{Alk}}=F^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}.\end{array}\right.\end{array}$$

It is worth noting that this general expression also applies to the specific case where the global Alk inventory is conserved ($F^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}=\mathrm{0}$).

2.1.2 Regional proxies for the spatial distribution of the pre-industrial riverine/burial-driven air–sea carbon flux

A direct relationship between carbon and Alk fluxes and the net global air–sea carbon flux can be established under the assumption of a steady-state air–sea carbon flux. However, this approach does not apply directly at the regional scale, where ocean circulation transports both Alk and DIC, and biogeochemical processes also generate regional sources and sinks.

The concept of the regional carbon:Alk budget imbalance

To gain a deeper understanding of the factors shaping the spatial distribution of the riverine/burial-driven air–sea carbon flux, we expand upon the theoretical framework previously outlined for the global scale (refer to Sect. 2.1.1) and adapt it as a proxy for application at regional scales. This is essential to understanding the extent to which specific regional carbon:Alk budget imbalances can drive the global air–sea carbon flux as well as deviations in carbon and Alk inventories.

The air–sea carbon flux calculated by applying Eq. (10) to the riverine and burial fluxes of a specific region can only be considered a potential air–sea carbon flux, i.e. a capacity to generate such a flux at the global scale, without any guarantee that it fully occurs within the same region. Due to ocean circulation and the associated transport of carbon and Alk, regional carbon:Alk budget imbalances in riverine and burial fluxes explain the regional distribution of the drivers of the global air–sea carbon flux. However, they only partially explain the regional distribution of the flux itself.

The inter-hemispheric flux gradient and transport

Understanding the spatial distribution of the pre-industrial riverine/burial-driven air–sea carbon flux is crucial for understanding the biases between observational and model-based estimates of the anthropogenic carbon sink. Yet, ocean circulation and carbon pumps within the ocean induce an asymmetry in the ocean on either side of the Intertropical Convergence Zone (ITCZ), which serves as an inter-hemispheric transport barrier (e.g. Murnane et al., 1999; Aumont et al., 2001; Resplandy et al., 2018). To assess the significance of this asymmetry on the air–sea carbon flux, particularly its components associated with riverine and burial fluxes, we provide two metrics for large-scale inter-hemispheric fluxes: (i) the inter-hemispheric air–sea carbon flux gradient (G), which is defined as the integrated net flux north of 20° N ($F_N^{\text{C, air–sea}}$) minus that south of 20° S ($F_S^{\text{C, air–sea}}$),

$$\begin{array}{}\text{(11)}& G=F_N^{\text{C, air–sea}}-F_S^{\text{C, air–sea}},\end{array}$$

and (ii) the inter-hemispheric ocean transport of carbon (T^C) and Alk (T^Alk), both directed northward, defined as the mean transport between 20° N ($T_N^{\mathrm{C}/\mathrm{Alk}}$) and 20° S ($T_S^{\mathrm{C}/\mathrm{Alk}}$):

$$\begin{array}{}\text{(12)}& T^{\mathrm{C}/\mathrm{Alk}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left(T_N^{\mathrm{C}/\mathrm{Alk}}+T_S^{\mathrm{C}/\mathrm{Alk}}\right).\end{array}$$

These two metrics rely on the subdivision of the ocean into two poleward basins, one south of 20° S and the other north of 20° N, separated by an intertropical basin (Fig. 3b, and see Appendix C). A decomposition of the inter-hemispheric air–sea carbon flux gradient (G) into components associated with carbon transport and with riverine and burial processes is provided in Appendix D2.

Table 2 Summary of the NEMO-PISCES sensitivity simulations with a short description (see Appendix B2.1 and Table B1 for more details).

2.2 Model and simulations

2.2.1 Model and configuration

As part of the NEMO (Nucleus for European Modelling of the Ocean) suite of models, here we used the marine biogeochemical model PISCES (Pelagic Interactions Scheme for Carbon and Ecosystem Studies) to take a fresh look at the pre-industrial air–sea carbon flux. This involved a comprehensive consideration of both the carbon and Alk budgets, with a specific focus on external fluxes, notably CaCO₃ burial. While globally resembling PISCES-v2, as detailed in Aumont et al. (2015) and utilized in IPSL-CM6A-LR (Boucher et al., 2020), we introduced two key modifications in PISCES: (i) an adjustment to the N-fixation parameterization, following Bopp et al. (2022), and (ii) an adaptation of the burial fraction of CaCO₃ to balance the Alk budget and conserve the global Alk inventory without necessitating an Alk restoring scheme – periodically restoring the global Alk inventory to a reference value by adding/removing the required amount, either uniformly or in a weighted manner (see Planchat et al., 2023, their Appendix A2, for details). Our simulations were conducted offline using a tripolar ORCA (orthogonal curvilinear ocean mesh) grid with a nominal resolution of 2° and included 30 vertical levels. The ocean physics were derived from pre-industrial simulations of IPSL-CM5A-LR (Dufresne et al., 2013, based on NEMOv3.2), with a repeated 500-year period, and a fixed and homogeneous atmospheric CO₂ concentration of 284 ppm at the ocean surface.

To ensure model stability and attainment of a steady state (i.e. stable air–sea carbon flux; e.g. Orr et al., 2017; see Fig. B2), all simulations presented below used the same initial conditions and have been run 2550 years after an initial 500-year spin-up using the standard configuration (Sect. 2.2.2). We calculated the carbon and Alk budgets related to their associated external sources/sinks using data from the last 50 years of the simulations. The carbon and Alk inventory deviations were estimated through linear regression over the same period.

2.2.2 Standard simulation (std) and its riverine/burial component

The standard simulation (referred to as “std”), based on the standard configuration described above, involves carbon and Alk riverine supply as well as OM and CaCO₃ burial. Riverine supply of carbon and Alk is based on output from the Global Erosion Model (GEM) of Ludwig et al. (1996) and considers both inorganic and organic carbon riverine discharge (0.37 and 0.14 PgC yr⁻¹, respectively). Carbon and Alk are added at river mouths using a monthly climatology that is applied recursively. The inorganic fraction is added as bicarbonate ions, thus affecting both DIC and Alk in a similar manner. The organic fraction is assumed to be fully labile and remineralizes instantaneously at the river mouth, thus impacting only DIC. This simulation also includes the burial of OM and CaCO₃ produced by pelagic organisms, which are exported to the ocean interior and only partially remineralized or dissolved in the water column and at the seafloor (e.g. Planchat et al., 2023). These combined fluxes constitute the riverine and burial fluxes (Eqs. 5 and 7), which, as introduced in Sect. 2.1.1, lead to the riverine/burial-driven air–sea carbon flux.

To isolate the riverine/burial-driven component of the air–sea carbon flux, a simulation without riverine and burial fluxes was conducted (referred to as “norivbur”), simulating only the component of the flux associated with the ocean carbon pumps. Indeed, while at the global scale, the net air–sea carbon flux directly corresponds to the riverine/burial-driven air–sea carbon flux (Eq. 2), at the regional scale (N, S, or I, Fig. 3b), the air–sea carbon flux ($F_{\mathrm{nat}.}^{\text{C, air–sea}}$) can be decomposed at first approximation into two components: one internal component linked to the functioning of the ocean carbon pumps ($F_{\mathrm{pumps}}^{\text{C, air–sea}}$) and a boundary component associated with the riverine and burial fluxes ($F_{\mathrm{riv}./\mathrm{bur}.}^{\text{C, air–sea}}$) – our primary focus:

$$\begin{array}{}\text{(13)}& F_{\mathrm{nat}.,N/S/I}^{\text{C, air–sea}}=F_{\mathrm{pumps},N/S/I}^{\text{C, air–sea}}+F_{\mathrm{riv}./\mathrm{bur}.,N/S/I}^{\text{C, air–sea}}.\end{array}$$

Subsequently, by taking the difference between the std and norivbur simulations, this allows us to determine the riverine/burial-driven air–sea carbon flux of our standard configuration:

$$\begin{array}{}\text{(14)}& F_{\mathrm{riv}./\mathrm{bur}.,N/S/I}^{\text{C, air–sea}}=F_{\mathrm{nat}.,N/S/I}^{\text{C, air–sea}}\left(\mathrm{std}\right)-F_{\mathrm{nat}.,N/S/I}^{\text{C, air–sea}}\left(\mathrm{norivbur}\right),\end{array}$$

where the “nat.” label was omitted since the simulations were conducted under pre-industrial conditions, and therefore, no anthropogenic component was included.

2.2.3 Sensitivity simulations

The set of sensitivity simulations considered covers a broad range of perturbations to the carbon and Alk riverine and burial fluxes. These simulations aim to assess the effects of different assumptions regarding these external fluxes on the riverine/burial-driven air–sea carbon flux (Table 2, Fig. 2b; see Appendix B2.1). Importantly, within the context of our study, the absolute values of the fluxes – whether they align with literature estimates or not – are not of primary concern. What matters are the relative differences between these values across simulations, which reflect the assumptions being tested and briefly outlined below (see Appendix B2.1 for further details).

First, we introduced variations in riverine discharge to account for uncertainties in its magnitude and partly unresolved characteristics (e.g. labile/refractory, organic/inorganic partitioning). By closing the Alk budget, these variations influenced CaCO₃ burial. In “rivref”, the OM riverine discharge was considered fully refractory (i.e. persisting on a timescale longer than that of ocean circulation), in contrast to the labile assumption in the standard simulation. We explored fully organic and inorganic riverine discharges in “rivorg” and “rivinorg”, respectively, and also increased riverine discharge by a factor of 1.5 in “riv1p5”, while maintaining the same partitioning as std. Second, to assess the effect of a non-conserved Alk inventory or an Alk restoration scheme, we disabled OM and CaCO₃ burial, artificially restoring Alk in “nobur-resto” or assuming a non-conserved Alk inventory in “nobur-diseq”. Third, we varied CaCO₃ burial to address uncertainties in its pre-industrial magnitude and spatial distribution. We added CaCO₃ burial/dissolution between the Atlantic and Pacific, balancing the Alk budget in “atlpac”, or not in “atlpac-diseq”, and we also added CaCO₃ burial in the tropics, resulting in an imbalanced Alk budget in “tropics-diseq”.

We classify our simulations as either “equilibrated” (suffix “-eq”) or “disequilibrated” (suffix “-diseq”), based on the conservation of global carbon and Alk inventories. In both cases, the air–sea carbon flux has reached a steady state. Equilibrated simulations are characterized by balanced global carbon and Alk budgets, resulting in conserved global ocean inventories over time. In contrast, disequilibrated simulations exhibit imbalanced budgets, leading to evolving global ocean inventories, which are therefore not conserved. It is important to note that the variations applied in our set of sensitivity simulations directly affected only carbon and Alk fluxes, while nutrient fluxes were held constant in order to avoid perturbing OM and CaCO₃ production. Finally, we report that at the global scale, for the standard simulation, Q_inv≃0.797 (Eq. 8), and this coefficient shows minimal variation across all sensitivity simulations considered (< 0.002 in absolute terms).

Figure 4 Description of the standard NEMO-PISCES simulation (std; see Fig. B3 for additional elements). (a) Map of the pre-industrial air–sea carbon flux, where positive values indicate ocean ingassing. (b) The zonally integrated air–sea carbon flux (dark blue) and the riverine/burial-driven air–sea carbon flux (aquamarine). When the riverine/burial-driven flux exceeds (is less than) the simulated one, the area in between is shaded in rose (cyan). (c) Partitioning of the riverine (orange) and burial (dark gold) fluxes by ocean region (southern, inter-tropical, and northern). The fluxes, in petagrams of carbon per year (PgC yr⁻¹) for carbon (in bold) and Alk (in normal font), are directed by arrows, with orientation indicating the sign and size reflecting the absolute magnitude of the flux. The regional partitioning of the riverine/burial-driven air–sea carbon flux (aquamarine) and of the potential air–sea carbon flux from regional carbon:Alk budget imbalances (light blue) is also provided above. (d) Partitioning of the integrated external sources and sinks of carbon (shaded) and Alk (hatched). The negative impact of OM burial on Alk is attributed to the release of ammonium when OM is remineralized at the seafloor rather than buried. Detailed descriptions of (c) and (d) can be found in supplementary figures (Figs. S1 and S2 on Zenodo).

3 Results

3.1 Pre-industrial air–sea carbon flux and its riverine/burial-driven component

In the pre-industrial era, the simulated air–sea carbon flux, derived from the standard simulation (std), remains stable at $-\mathrm{0.27}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$ (see Fig. B4a), indicating a net global ocean outgassing. However, the distribution of the surface air–sea carbon flux exhibits considerable heterogeneity among different regions, primarily driven by ocean circulation patterns (Fig. 4a). Regions characterized by carbon-rich deep-water upwellings, such as equatorial and Southern Ocean upwelling zones, tend to show carbon outgassing. Conversely, poleward heat transport within the ocean, exemplified by western boundary currents like the Gulf Stream, promotes carbon absorption as surface waters cool. Consequently, the outcome is a pronounced meridional air–sea carbon flux gradient, with ingassing in the Northern Hemisphere ($+\mathrm{0.57}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$), outgassing in the inter-tropical zone ($-\mathrm{0.91}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$), and minimal outgassing in the Southern Hemisphere ($-\mathrm{0.06}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$), primarily due to the impact of Southern Ocean upwelling between 45–65° S (Fig. 4b, c). The inter-hemispheric air–sea carbon flux gradient (G; Sect. 2.1.2) is $+\mathrm{0.51}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$, while the inter-hemispheric transport of carbon and Alk (T^C and T^Alk; Sect. 2.1.2) amounts to $-\mathrm{0.35}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$ and $-\mathrm{0.07}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$, respectively (see Fig. B4a).

The air–sea carbon flux can be subsequently decomposed into a component associated with the functioning of the ocean carbon pumps and a component associated with riverine and burial fluxes (Sect. 2.2.2). More specifically, the simulation without riverine and burial fluxes (norivbur) is used to isolate the distribution of air–sea carbon fluxes associated with the internal functioning of the ocean carbon pumps. Indeed, regional air–sea carbon fluxes are primarily influenced by these pumps, which establish and sustain vertical and horizontal carbon gradients within the ocean (e.g. Sarmiento and Gruber, 2006; Murnane et al., 1999; Aumont et al., 2001; Resplandy et al., 2018). Thus, both the physical pump (involving ocean circulation and air–sea carbon exchange) and the biological pump (comprising processes like production, export, and the remineralization/dissolution of OM and CaCO₃) play pivotal roles in elucidating the overall distribution of the air–sea carbon flux. These air–sea carbon fluxes exhibit significant ingassing in the Northern Hemisphere (+0.67 PgC yr⁻¹) and outgassing in the inter-tropical zone (−0.79 PgC yr⁻¹), with minimal ingassing in the Southern Hemisphere (+0.16 PgC yr⁻¹; see Table B2). Overall, the air–sea carbon flux associated with the oceanic carbon pumps is expected to be net zero when integrated at the global scale, although norivbur shows a small residual component (+0.05 PgC yr⁻¹; see Fig. B4a). This residual component is attributed to a residual carbon budget imbalance due to internal ocean processes (see Appendix B2.4).

Finally, by taking the difference between our standard simulation and the simulation without riverine and burial fluxes (std minus norivbur), we isolate the component of interest, i.e. that induced by riverine and burial fluxes (Sect. 2.2.2). This riverine/burial-driven air–sea carbon flux results in a net global outgassing of 0.31 PgC yr⁻¹, distributed among the northern, inter-tropical, and southern regions as follows: 0.10, 0.12, and 0.10 PgC yr⁻¹ (Fig. 4c).

Figure 5 The role of riverine and burial fluxes of carbon and Alk in determining the pre-industrial air–sea carbon flux. (a, b) Comparison between the net global air–sea carbon flux^a and (a) the integrated riverine fluxes of carbon or (b) the integrated riverine and burial fluxes of carbon. When the net air–sea carbon flux balances the considered external fluxes (on the 1:1 line), simulation names are indicated in black. This applies to (a) simulations that do not account for burial and conserve the global Alk inventory (norivbur and nobur-resto) and (b) all simulations conserving the global Alk inventory (excluding nobur-diseq, atlpac-diseq, and tropics-diseq). (c) Theoretical framework that accounts for Alk and carbon budgets to reconstruct the net air–sea carbon flux. The net air–sea carbon flux (filled contours) is determined by multiplying the integrated riverine and burial fluxes of Alk (x axis) by Q_inv and then subtracting the integrated riverine and burial fluxes of carbon (y axis). The deviation of the net air–sea carbon flux from this relationship in the NEMO-PISCES sensitivity simulations is small (less than 0.01 PgC yr⁻¹ for all, except nobur-resto: less than 0.03 PgC yr⁻¹). Simulations with a conserved global Alk inventory align with the zero x-axis line. The most recent carbon and Alk budgets (Table 3) provide estimates of riverine and burial fluxes of carbon^b and Alk^c, as shown at the top and on the right in grey. The net air–sea carbon flux reconstructed from these flux estimates is indicated as grey rectangles, with confidence intervals at 75 %, 50 %, and 25 %, and projected on the colour bar. ^a The net air–sea carbon flux of the NEMO-PISCES sensitivity simulations was adjusted for their respective residual carbon budget imbalances (see Appendix B2.4). ^b This distribution also includes fluxes from groundwater discharge. ^c This distribution also includes fluxes from anaerobic processes, groundwater discharge, and reverse weathering.

3.2 The global riverine/burial-driven air–sea carbon flux

3.2.1 Role of sediment burial fluxes

Accounting for the riverine carbon input alone in the standard simulation ($+\mathrm{0.52}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$) is insufficient to explain the simulated air–sea carbon outgassing of 0.27 PgC yr⁻¹. Burial-associated carbon fluxes, from both OM ($-\mathrm{0.17}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$) and CaCO₃ ($-\mathrm{0.04}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$; Figs. 4d and B4a), act to partially offset this input, thereby reducing the net outgassing to 0.31 PgC yr⁻¹.

The importance of burial fluxes in driving the riverine/burial-driven air–sea carbon flux is furthermore exemplified by our set of sensitivity simulations (Fig. 5a, b). Increasing the riverine input by a factor of 1.5, while maintaining its partitioning (riv1p5), results in an increase in carbon outgassing of 0.17 PgC yr⁻¹. This is less than the increase in riverine carbon discharge (+0.26 PgC yr⁻¹); as part of this, additional carbon is buried as CaCO₃ (−0.09 PgC yr⁻¹) to maintain a balanced Alk budget (see Fig. B4b). Similarly, a change in the partitioning of the riverine input between organic and inorganic forms (rivorg and rivinorg) does not alter the total magnitude of the riverine carbon input compared to the standard configuration, but it does affect the air–sea carbon outgassing. It reaches 0.47 and 0.20 PgC yr⁻¹, respectively (see Table B2), as the associated decrease (−0.38 PgC yr⁻¹) or increase (+0.14 PgC yr⁻¹) in the riverine Alk discharge relative to std (see Fig. B4b) leads to corresponding changes in CaCO₃ burial (+0.19 and −0.07 PgC yr⁻¹, respectively) to maintain a balanced Alk budget. This highlights the pivotal role of CaCO₃ burial in shaping the air–sea carbon flux under the assumption of a balanced Alk budget, where riverine Alk inputs are offset by CaCO₃ burial (Fig. 4d).

However, such a carbon budget – which deduces the pre-industrial air–sea carbon flux from riverine and burial fluxes of carbon – is only valid under the condition of a balanced Alk budget (Fig. 5b). When this assumption does not hold, it becomes necessary to account for both the carbon and Alk budgets to correctly assess the riverine/burial-driven air–sea carbon flux (Sect. 2.1.1).

3.2.2 Impact of an imbalanced alkalinity budget

The possibility of an imbalanced Alk budget during the pre-industrial era has been hypothesized multiple times over the past 3 decades (e.g. Milliman, 1993; Middelburg et al., 2020; Cartapanis et al., 2018; Boudreau et al., 2018). The simulations atlpac-diseq and tropics-diseq allow us to assess the implications of such a deviation in the global Alk inventory (Table 2; see also Table B1), by controlling both the magnitude and spatial distribution of CaCO₃ burial in a way that better reflects current paleoceanographic reconstructions (e.g. Cartapanis et al., 2018). Both simulations implement an imbalanced Alk budget (−0.10 PgC yr⁻¹) via additional CaCO₃ burial, either in the deep Atlantic (atlpac-diseq) or in the shallow tropics to represent coral reef processes (tropics-diseq). They lead to the same increase in steady-state air–sea carbon outgassing relative to std (+0.07 PgC yr⁻¹). This outcome may seem counterintuitive when applying a simple carbon budget, since both simulations prescribe extra carbon removal from the ocean (to the sediments) yet result in enhanced carbon loss to the atmosphere (Fig. 5b). In fact, the associated outgassing leads to a net decrease in the global ocean carbon inventory (−0.16 PgC yr⁻¹), which exceeds, in absolute terms, the additional CaCO₃ burial (−0.10 PgC yr⁻¹; see Fig. B4b).

3.2.3 Validating the governing equation of the pre-industrial riverine/burial-driven air–sea carbon flux

Importantly, even with an imbalanced Alk budget that drives deviations in the global carbon and Alk inventories, the ocean can maintain a steady-state air–sea carbon flux (Fig. B2). Overall, the theoretical framework introduced in Sect. 2.1.1 is fully validated by our set of sensitivity simulations. At the pre-industrial steady state, the net air–sea carbon flux (F^{C, air–sea}) can be expressed as the product of the integrated riverine and burial Alk flux (F^Alk) multiplied by Q_inv, minus the integrated riverine and burial carbon flux ($F^{\mathrm{C},\mathrm{bur}./\mathrm{riv}.}$):

$$\begin{array}{}\text{(15)}& F^{\text{C, air–sea}}=\underset{D^{\mathrm{C}}}{\underbrace{Q_{\mathrm{inv}}\cdot \stackrel{D^{\mathrm{Alk}}}{\overbrace{F^{\mathrm{Alk}}}}}}-F^{\mathrm{C},\mathrm{bur}./\mathrm{riv}.},\end{array}$$

where D^C and D^Alk represent the global carbon and Alk inventory deviations, respectively (Fig. 5c). This formulation highlights the critical role of pre-industrial Alk budget assumptions – as well as the persistent uncertainties and unknowns – in estimating the pre-industrial air–sea carbon flux.

Figure 6 Spatial distribution of the riverine/burial-driven air–sea carbon flux. Comparison between the riverine/burial-driven air–sea carbon flux (y axis, PgC yr⁻¹), the fraction of this flux occurring south of 20° S (x axis, %), and its inter-hemispheric gradient (colour dots, PgC yr⁻¹). The fraction of this flux occurring south of 20° S is also shown for Aumont et al. (2001) and Lacroix et al. (2020) (black stars), assuming the same global riverine/burial-driven air–sea carbon flux as our standard simulation (std).

3.3 The spatial distribution of the riverine/burial-driven air–sea carbon flux

3.3.1 Contrasting regional fluxes

The inter-hemispheric gradient of the pre-industrial air–sea carbon flux is primarily controlled by ocean interior processes and the functioning of the ocean carbon pumps. Specifically, in an ocean without any riverine and burial carbon fluxes (norivbur), the inter-hemispheric gradient amounts to +0.51 PgC yr⁻¹ (see Fig. B4a). The biological pump contributes to carbon uptake in the Northern Hemisphere through surface biological activity and leads to carbon release in the Southern Hemisphere due to the upwelling of carbon-rich deep waters, as documented in previous studies (e.g. Murnane et al., 1999; Aumont et al., 2001; Resplandy et al., 2018). When subtracting the gradient estimated from simulation norivbur to all other sensitivity simulations, we find that only a fraction of the inter-hemispheric air–sea carbon flux gradient is accounted for by riverine and burial fluxes (ranging from −0.18 to +0.11 PgC yr⁻¹).

Our set of sensitivity simulations, which explore various assumptions about riverine and burial fluxes (Sect. 2.2.3), encompass the uncertainty range associated with the inter-hemispheric gradient of the riverine/burial-driven air–sea carbon flux (Fig. 6). The main point of contention regarding this gradient lies in the fraction of the flux occurring in the Southern Hemisphere, where the largest discrepancies in estimates of the anthropogenic carbon sink between pCO₂-based methods and model simulations were located in GCBs (from 2018 to 2022; e.g. Hauck et al., 2020; Friedlingstein et al., 2022b), before being mostly shifted to the inter-tropical region (since 2023; e.g. Friedlingstein et al., 2023, 2025). In our simulations, the fraction of the flux occurring in the Southern Hemisphere ranges from less than 5 % (nobur-diseq) to more than 50 % (rivref). By comparison, it was estimated at 49 % (Aumont et al., 2001) and then revised to 14 % (Lacroix et al., 2020) in the GCBs (Table 1 and Fig. 6) and even as low as 4 % in the literature (Jacobson et al., 2007). This is particularly intriguing, as one might expect this distribution to be primarily governed by the strength of the meridional overturning circulation – and its role in transporting riverine/burial-derived carbon southward – yet our sensitivity simulations, despite identical ocean dynamics, reveal highly contrasting distributions.

3.3.2 Influencing factors

There is no direct correlation between the magnitude of the riverine/burial-driven air–sea carbon flux and the proportion of this flux occurring south of 20° S (Fig. 6). Notably, the substantial uncertainty on the refractory nature of organic riverine discharge (e.g. Aumont et al., 2001; Gruber et al., 2009) is demonstrated to result in a significant shift in the proportion of the riverine/burial-driven air–sea carbon flux occurring in the Southern Ocean (54 % in rivref vs. 31 % in std; Fig. 6), even though the total flux remains the same. Conversely, when the riverine discharge is increased by 50 % (riv1p5), the distribution of the riverine/burial-driven air–sea carbon flux remains unchanged compared to std, while the total outgassing increases from 0.32 PgC yr⁻¹ (std) to 0.49 PgC yr⁻¹ (riv1p5; see Table B2).

The decoupling between the magnitude of the net riverine/burial-driven air–sea carbon flux and its inter-hemispheric gradient is primarily linked to the distribution, both horizontally and vertically, of the carbon:Alk budget imbalance resulting from riverine and burial fluxes (Sect. 2.1.2). When an excess of CaCO₃ burial is considered at the bottom of the Atlantic (primarily in the Northern Hemisphere; atlpac-diseq), the resulting impact of the carbon:Alk budget imbalance on the riverine/burial-driven air–sea carbon flux occurs remotely, in the Southern Hemisphere, due to the meridional overturning circulation. This results in a relative outgassing compared to std (−0.07 PgC yr⁻¹) and an increase in the inter-hemispheric riverine/burial-driven air–sea carbon flux gradient (+0.07 PgC yr⁻¹; see Fig. B4b). Conversely, when the surplus of CaCO₃ burial is in the shallow tropics (tropics-diseq), the riverine/burial-driven air–sea carbon flux anomaly compared to std is equivalent to the one reported for atlpac-diseq, but the inter-hemispheric gradient is this time nearly not impacted relative to std (+0.01 PgC yr⁻¹) since the flux anomaly is concentrated in the shallow tropics, primarily affecting the regional air–sea carbon flux. Similarly, flux anomalies resulting from carbon:Alk budget imbalances with respect to the riverine fluxes tend to manifest regionally (i.e. primarily in the Northern Hemisphere): (i) a fully organic riverine discharge (rivorg) leads to a relative outgassing compared to std (−0.19 PgC yr⁻¹), aligned with a decrease in the inter-hemispheric gradient (−0.18 PgC yr⁻¹), and (ii) a fully inorganic riverine discharge (rivinorg) leads to a relative ingassing compared to std (+0.07 PgC yr⁻¹), aligned with an increase in the inter-hemispheric gradient (+0.05 PgC yr⁻¹; see Fig. B4b).

4 Proof-of-concept applications and discussions

4.1 The global flux

4.1.1 Approach

The theoretical framework introduced in this study (Sect. 2.1.1) has been validated by our set of sensitivity simulations (Sect. 3.2). It is therefore now possible to estimate the global magnitude of the pre-industrial riverine/burial-driven air–sea carbon flux and to investigate the associated global carbon and Alk inventory deviations (Eq. 10) based on existing carbon and Alk budgets, which encompass all external oceanic sources and sinks of carbon and Alk. For consistency with the literature, we rely on the most recent carbon (Regnier et al., 2022) and Alk (Middelburg et al., 2020) budgets, even though they were derived independently and are partly inconsistent (Table 3). We carefully accounted for the uncertainties and extreme values associated with the various external sources/sinks of carbon and Alk (Table 3).

Table 3 Literature-based estimates of riverine/burial fluxes of carbon and Alk, from the most recent carbon and Alk budgets (Sect. 4.1.1 and 4.2.1, and see Fig. E3), including the calculation of the corresponding air–sea carbon flux, as well as carbon and Alk content deviations. The values are presented in petagrams of carbon per year. Values in brackets represent the uncertainty or extreme range, while the bold value indicates the best estimate or average. The intervals are arranged with the smallest absolute value first, except when both positive and negative values are present in the range.

^a Including fluxes from groundwater discharge and anaerobic processes. ^b Including fluxes from reverse weathering.

4.1.2 Findings

Using the theoretical framework introduced in this paper and literature-based estimates of riverine/burial fluxes of carbon and Alk, based on the most recent carbon and Alk budgets, we derive, from Eq. (15), a pre-industrial riverine/burial-driven air–sea carbon flux estimate of −0.49 $\left[-\mathrm{0.34};-\mathrm{0.70}\right]$ PgC yr⁻¹ (Table 3 and Fig. 5c). This pre-industrial riverine/burial-driven air–sea carbon flux is associated with global carbon and Alk inventory deviations of 0.06 $\left[-\mathrm{0.05};\mathrm{0.11}\right]$ PgC yr⁻¹ and 0.07 $\left[-\mathrm{0.06};\mathrm{0.14}\right]$ PgC yr⁻¹, respectively (see Fig. E1). This estimate is based on an integrated external flux of 0.55 [0.45;0.65] PgC yr⁻¹ for carbon and 0.07 $\left[-\mathrm{0.06};\mathrm{0.14}\right]$ PgC yr⁻¹ for Alk.

4.1.3 Discussion

This new estimation of the pre-industrial riverine/burial-driven air–sea carbon flux represents a downward revision of the latest value of $-\mathrm{0.65}\pm \mathrm{0.30}$ PgC yr⁻¹ currently adopted in the GCB (Friedlingstein et al., 2025), which was derived from a comprehensive assessment of the global land-to-ocean carbon continuum (Regnier et al., 2022). Applying our revised estimate in the calculation of the anthropogenic carbon uptake based on pCO₂-based methods would reduce the overall discrepancy between observation-based and model-derived oceanic carbon uptake estimates by 0.16 PgC yr⁻¹ over the historical period, thus alleviating a portion of the present offset (Fig. 1; Friedlingstein et al., 2025).

The discrepancy between our reassessment of the riverine/burial-driven outgassing and the value currently used in the GCB underscores the crucial importance of clearly defining ocean boundary conditions and the pressing need to develop a combined and consistent carbon and Alk budget for the ocean to achieve a robust estimate. Part of this discrepancy arises because atmospheric carbon uptake by continental shelves (0.10 PgC yr⁻¹; Regnier et al., 2022) is fully integrated into our net pre-industrial riverine/burial-driven air–sea carbon flux as we also consider OM and CaCO₃ burial in these regions, reducing this flux by 0.10 PgC yr⁻¹.

The current inconsistencies between the independently developed carbon and Alk budgets make our estimate less robust and highlight the need for a combined revision of both. Beyond the 0.10 PgC yr⁻¹ reduction in outgassing due to differing ocean boundary definitions relative to GCB, the remaining 0.06 PgC yr⁻¹ decrease in our new estimate is linked to a slight imbalance in the Alk budget (+0.07 PgC yr⁻¹ Middelburg et al., 2020). However, the discrepancy in CaCO₃ burial estimates between the most recent carbon and Alk budgets (Regnier et al., 2022; Middelburg et al., 2020) would translate into a 0.22 PgC yr⁻¹ difference in the Alk budget (Table 3). If the carbon flux associated with CaCO₃ burial were aligned with the Alk budget from Middelburg et al. (2020), the outgassing would decrease by an additional 0.18 PgC yr⁻¹. Conversely, aligning the Alk flux associated with CaCO₃ burial with the carbon value from Regnier et al. (2022) would reduce the outgassing by 0.11 PgC yr⁻¹. Thus, reconciling CaCO₃ burial fluxes in both carbon and Alk budgets is expected to further reduce the current outgassing offset (Friedlingstein et al., 2025). Establishing a combined and internally consistent carbon and Alk budget is therefore essential to confidently reassess the pre-industrial outgassing within the theoretical framework presented here.

4.2 The flux distribution

4.2.1 Approach

The set of sensitivity simulations conducted to validate our theoretical framework spans a wide range of assumptions regarding riverine and burial fluxes of carbon and Alk. This provides all the necessary tools to reassess the spatial distribution of the riverine/burial-driven air–sea carbon flux. As in our global estimate (Sect. 4.1), this reassessment strategy is grounded in the most recent global budgets of carbon and Alk. By logically selecting and weighting some of our sensitivity simulations, we construct a composite simulation whose riverine and burial fluxes match those reported in these global budgets (Table 3). In this way, the composite simulation also combines the associated air–sea carbon fluxes, both at the global scale and regionally. It is this regional aspect that enables a revised estimate of the spatial distribution of the pre-industrial riverine/burial-driven air–sea carbon flux.

First, for each literature-based estimate of the external sources and sinks of carbon and Alk, we constructed a skewed Gaussian probability density function (PDF) that captures the median/mean value and the reported uncertainty range. This was achieved in two steps for each literature estimate of the various external sources/sinks of carbon and Alk. A triangular distribution was first generated using the estimated central value and minimum/maximum bounds via the “random.triangular” function from the Python library numpy. This triangular distribution was then fitted with a skewed normal PDF using the “stats.skewnorm.fit” function from the scipy library. This approach allowed us to preserve the essential characteristics of the literature values (median/mean and extremes) while working with continuous distributions.

Second, we constructed a composite simulation that isolates the effect of riverine and burial fluxes on the air–sea carbon flux (i.e. excluding the influence of internal carbon pumps). This was achieved by linearly combining a subset of our sensitivity simulations. Throughout the remainder of the paper, we refer to the composite simulated estimate as the pre-industrial riverine/burial-driven air–sea carbon flux derived from this combined simulation. The following four-step workflow is designed to ensure that the riverine and burial fluxes in the composite simulation are consistent with the latest literature estimates (Middelburg et al., 2020; Regnier et al., 2022; see Fig. E2 and Table 3).

• Step 1: We initialized our composite simulation by isolating the effect of riverine and burial fluxes on the air–sea carbon flux, removing the contribution of internal carbon pumps. This was done by subtracting the norivbur simulation from the standard simulation (“std”–“norivbur”).

• Step 2: Next, we adjusted the carbon fluxes associated with riverine discharge and OM burial to match literature estimates. This was achieved by weighting the simulation where riverine inputs were increased by a factor of 1.5 (“riv1p5”–“std”), as the overall riverine flux amplitude was the first variable that needed to be tuned. At the end of this stage, the composite simulated estimate was a linear combination of “std”–“norivbur” and “riv1p5”–“std”.

• Step 3: We then adjust the Alk fluxes associated with riverine discharge and OM burial to match literature estimates. This was achieved by weighting the simulation where all riverine discharge was considered inorganic (“rivinorg”–“std”), which did not alter the carbon values already matched in Step 2. The composite simulated estimate became a linear combination of the result from Step 2 and “rivinorg”–“std”.

• Step 4: Finally, we ensured that Alk fluxes associated with CaCO₃ burial also matched the literature estimate. This was done using the simulation with enhanced CaCO₃ burial/dissolution and a global Alk imbalance (“atlpac-diseq”–“std”), without affecting the fluxes adjusted in previous steps. Given the 2:1 stoichiometric ratio between Alk and DIC in CaCO₃ processes, this step simultaneously ensured consistency for both the carbon and Alk components of CaCO₃ burial. The final composite estimate was a linear combination of the result from Step 3 and “atlpac-diseq”–“std”.

• Extra step (correction): Due to inconsistencies between the most recent carbon and Alk budgets – specifically in the CaCO₃ burial flux (Table 3) – an additional correction step was required. This correction, applied similarly to Step 3, again uses “rivinorg”–“std” to consider an increased carbon sink via CaCO₃ burial, while maintaining Alk balance. This adjustment targets only the excess CaCO₃ burial of carbon needed to reconcile our composite simulation with the carbon budget from Regnier et al. (2022). Note that this step would not be necessary if the carbon and Alk budgets were internally consistent.

In summary, this composite simulated estimate, built as a weighted linear combination of targeted sensitivity simulations and constrained by the latest literature estimates of riverine and burial fluxes, provides a spatially explicit representation of the pre-industrial riverine/burial-driven air–sea carbon flux. By design, the integrated value of this flux in the composite simulation is consistent with that obtained by applying the theoretical framework to the most recent carbon and Alk budgets from the literature (Sect. 4.1).

Figure 7 Description of the composite simulated estimate resulting from a linear combination of the NEMO-PISCES sensitivity simulations and literature-based estimates of riverine/burial fluxes of carbon and Alk (Sect. 4.2.1, and see Fig. E3). (a) PDF illustrating the total riverine and burial fluxes of carbon (shaded) and Alk (hatched) in the composite simulated estimate, along with the associated PDF for the resulting deviation in carbon (solid) and Alk (dashed) content. (b, c) PDFs of the net air–sea carbon flux and the inter-hemispheric air–sea carbon flux gradient. Within each of these sub-panels, the PDF associated with no riverine and burial fluxes of carbon and Alk (norivbur; cyan line) is juxtaposed with the one corresponding to only riverine and burial fluxes of carbon and Alk (composite simulated estimate minus norivbur; aquamarine) to obtain the total value (composite simulated estimate; dark blue). Further details on the residual component where no riverine and burial fluxes are considered are explained in Appendix B2.4. (d, e, f) The associated spatial distribution for the southern, inter-tropical, and northern regions: (d) without riverine and burial fluxes of carbon, (e) exclusively related to riverine and burial fluxes of carbon and Alk, and (f) the overall distribution. In (e), the percentage of each component is provided in brackets.

4.2.2 Findings

The construction of a composite simulated estimate resulting from a linear combination of our sensitivity simulations to align with the literature-based estimates for carbon and Alk budgets (Fig. 7a; Sect. 4.2.1, and see Fig. E3) enables a reassessment of the distribution of this riverine/burial-driven air–sea carbon outgassing (0.15±0.13, 0.20±0.10, and 0.16±0.08 PgC yr⁻¹ for the southern, inter-tropical, and northern regions, respectively; Fig. 7e). The uncertainty associated with these values is primarily linked to the uncertainties/extremes in literature-based estimates (see Fig. E1). Such a distribution implies that 29 % of the outgassing occurs in the southern region, 40 % in the inter-tropical region, and 31 % in the northern region.

4.2.3 Discussion

The distribution we found corresponds to an intermediate distribution compared to those adopted in the GCB over time, falling between the most recent estimate of 14 %, 64 %, and 22 % (Lacroix et al., 2020) and the earlier estimate of 49 %, 25 %, and 26 % (Aumont et al., 2001, Table 1). This would partially confirm the reduction in the discrepancy between pCO₂-based and model estimates in the southern region, while avoiding the introduction of a bias in the inter-tropical region, as noted in GCB 2023 (Friedlingstein et al., 2023) compared to previous GCBs (e.g. Friedlingstein et al., 2022b).

By summing the fluxes from the composite simulated estimate and the simulation without riverine and burial fluxes (norivbur), the total inter-hemispheric air–sea carbon flux gradient can be obtained. Notably, this amounts to 0.50±0.15 PgC yr⁻¹ (Fig. 7c), which aligns with the inter-hemispheric CO₂ concentration gradient in the atmosphere between the South Pole and Mauna Loa during the pre-industrial era. It was historically assessed at +0.82 ppm (Keeling et al., 1989) and more recently reevaluated at $+\mathrm{0.55}\pm \mathrm{0.15}$ ppm (Resplandy et al., 2018) through interpolation of atmospheric CO₂ concentration measurements.

A more comprehensive characterization of riverine and burial fluxes of carbon and Alk remains a critical challenge for accurately constraining the spatial distribution of the riverine/burial-driven air–sea carbon flux. This is particularly true for the fate of terrestrial organic carbon and its associated lability, which remains highly uncertain (Aumont et al., 2001; Jacobson et al., 2007; Gruber et al., 2009). Nevertheless, the approach proposed in this study is flexible and can accommodate future revisions of these external fluxes. Fundamentally, the selection of sensitivity simulations used to construct the composite simulation (Sect. 4.2.1; see also Fig. E2) can be revisited as scientific understanding progresses or as model representations evolve. For instance, in NEMO-PISCES, burial tends to occur predominantly near coastal margins. To counterbalance this biased feature in the composite simulation, we selected the sensitivity simulation with extra CaCO₃ burial in the deep Atlantic basin (atlpac-diseq), rather than the one with increased burial in the shallow tropics (tropics-diseq). A limitation of our approach is that a substantial revision of the spatial distribution of a flux – such as riverine inputs – would require rerunning a simulation, as it cannot be addressed through our current framework alone.

Nonetheless, the use of sensitivity simulations to build a composite simulation underscores the method's potential for reassessing the distribution of the pre-industrial air–sea carbon flux. By drawing from a set of pre-existing simulations and grounding the reassessment in the theoretical framework developed in this study, the spatial pattern of the flux can be revised in a consistent and coherent manner, without the need for additional model runs. This approach is particularly well suited for model intercomparison exercises, as it allows for efficient re-evaluation of regional fluxes and contributes to reducing biases linked to differences in ocean circulation or biogeochemical parameterizations across models.

5 Conclusion and perspectives

We have offered a fresh perspective on the pre-industrial air–sea carbon flux by leveraging the ocean alkalinity budget. The theoretical framework we introduced, validated through sensitivity simulations conducted with NEMO-PISCES, demonstrates both its robustness and its practical relevance for assessing the riverine/burial-driven pre-industrial air–sea carbon flux in the context of the GCB.

Through two proof-of-concept applications, we demonstrate the potential of this theoretical framework to identify biases between observation-based and model-derived estimates of the oceanic carbon sink at both global and regional scales and to partially correct persistent offsets. In the first application, we revisit the global magnitude of the pre-industrial riverine/burial-driven air–sea carbon flux using existing carbon and alkalinity budgets. This yields a simple and rapid method for reassessment whenever these budgets are revised. In the second application, we propose a method to reassess the spatial distribution of the pre-industrial riverine/burial-driven air–sea carbon flux. This is achieved by constructing a composite simulation, based on a linear combination of sensitivity simulations, that aligns with both carbon and alkalinity budgets. This approach is particularly well suited for model intercomparison exercises, as it enables efficient reassessment of regional fluxes while helping to mitigate biases related to ocean physics or biogeochemical parameterizations.

These flexible applications now call for four key efforts from the community regarding the pre-industrial riverine/burial-driven air–sea carbon flux:

• i.

  To reduce uncertainty in its global magnitude:

  • Clarify the definition of ocean domain boundaries at the coastal interface within the land-to-ocean continuum, where multiple fluxes intersect (riverine discharge and part of organic matter and CaCO₃ burial).

  • Establish a combined and internally consistent carbon and alkalinity budget, as current independently developed estimates remain inconsistent (e.g. CaCO₃ burial).

• ii.

  To reduce uncertainty in its regional distribution:

  • Improve our understanding of the intrinsic properties of riverine and burial fluxes (e.g. the fate of terrestrial organic matter).

  • Promote intermodel comparison efforts to identify systematic biases and improve robustness across modelling approaches.

Appendix A: Theoretical framework

As a complement to the theoretical framework introduced in Sect. 2.1.1, here we outline how to derive Q_inv, the inverse of the “isocapnic quotient” approximation introduced by Humphreys et al. (2018). Specifically, we develop the method proposed by Planchat et al. (2023) and subsequently demonstrate its full consistency with the approach employed by Humphreys et al. (2018).

For a fixed salinity (S) and temperature (T), pCO₂ – the partial pressure of CO₂ in seawater – can be differentiated as follows:

$$\begin{array}{}\text{(A1)}& \mathrm{d}p{\mathrm{CO}}_{\mathrm{2}}={\left.{\displaystyle \frac{\partial p{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{Alk}}}\right|}_{\mathrm{DIC},S,T}\cdot \mathrm{dAlk}+{\left.{\displaystyle \frac{\partial p{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{DIC}}}\right|}_{\mathrm{Alk},S,T}\cdot \mathrm{dDIC}.\end{array}$$

Assuming pCO₂ is fixed – for instance, at equilibrium with atmospheric CO₂ – leads to

$$\begin{array}{}\text{(A2)}& {\left.{\displaystyle \frac{\mathrm{dDIC}}{\mathrm{dAlk}}}\right|}_{p{\mathrm{CO}}_{\mathrm{2}},S,T}=-{\left.{\displaystyle \frac{\partial p{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{Alk}}}\right|}_{\mathrm{DIC},S,T}\cdot {\left.{\displaystyle \frac{\partial \mathrm{DIC}}{\partial p{\mathrm{CO}}_{\mathrm{2}}}}\right|}_{\mathrm{Alk},S,T}.\end{array}$$

Yet, pCO₂ is defined by

$$\begin{array}{}\text{(A3)}& p{\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{K_{\mathrm{2}}}{K_{\mathrm{0}}\cdot K_{\mathrm{1}}}}\cdot {\displaystyle \frac{{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}^{\mathrm{2}}}{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}},\end{array}$$

where K₀, K₁, and K₂ are the stoichiometric equilibrium/dissociation constants of the CO₂ system (e.g. Sarmiento and Gruber, 2006).

We then introduce the simplifying assumption:

$$\begin{array}{}\text{(A4)}& \left\{\begin{array}{l}\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]\simeq \mathrm{2}\mathrm{DIC}-\mathrm{Alk}\\ \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\simeq \mathrm{Alk}-\mathrm{DIC}.\end{array}\right.\end{array}$$

This assumption is reasonable given that $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ together typically account for over 99 % of DIC and over 97 % of Alk (Sarmiento and Gruber, 2006; Humphreys et al., 2018). Under this assumption, Eq. (A3) can be approximated as

$$\begin{array}{}\text{(A5)}& p{\mathrm{CO}}_{\mathrm{2}}\simeq {\displaystyle \frac{K_{\mathrm{2}}}{K_{\mathrm{0}}\cdot K_{\mathrm{1}}}}\cdot {\displaystyle \frac{{\left(\mathrm{2}\mathrm{DIC}-\mathrm{Alk}\right)}^{\mathrm{2}}}{\mathrm{Alk}-\mathrm{DIC}}}.\end{array}$$

Accordingly, the partial derivatives of pCO₂ with respect to Alk and DIC at constant pCO₂, S, and T become

$$\begin{array}{}\text{(A6)}& \left\{\begin{array}{l}{\left.\frac{\partial p{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{Alk}}\right|}_{\mathrm{DIC},S,T}\simeq \frac{K_{\mathrm{2}}}{K_{\mathrm{0}}\cdot K_{\mathrm{1}}}\cdot \frac{-\mathrm{Alk}\cdot \left(\mathrm{2}\mathrm{DIC}-\mathrm{Alk}\right)}{{\left(\mathrm{Alk}-\mathrm{DIC}\right)}^{\mathrm{2}}}\\ {\left.\frac{\partial p{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{DIC}}\right|}_{\mathrm{Alk},S,T}\simeq \frac{K_{\mathrm{2}}}{K_{\mathrm{0}}\cdot K_{\mathrm{1}}}\cdot \frac{\left(\mathrm{3}\mathrm{Alk}-\mathrm{2}\mathrm{DIC}\right)\cdot \left(\mathrm{2}\mathrm{DIC}-\mathrm{Alk}\right)}{{\left(\mathrm{Alk}-\mathrm{DIC}\right)}^{\mathrm{2}}}\end{array}\right..\end{array}$$

Substituting these expressions into Eq. (A2) gives

$$\begin{array}{}\text{(A7)}& {\left.{\displaystyle \frac{\mathrm{dDIC}}{\mathrm{dAlk}}}\right|}_{p{\mathrm{CO}}_{\mathrm{2}},S,T}={\displaystyle \frac{\mathrm{1}}{Q}}=Q_{\mathrm{inv}}\simeq {\displaystyle \frac{\mathrm{Alk}}{\mathrm{3}\mathrm{Alk}-\mathrm{2}\mathrm{DIC}}}.\end{array}$$

It is worth noting that the same expression can also be derived following the method presented in Humphreys et al. (2018, their Appendix C). Using the same approximation as in Eq. (A4), they arrive at the following form for Alk (see their Eq. C.6):

$$\begin{array}{}\text{(A8)}& \mathrm{Alk}\simeq \mathrm{2}\mathrm{DIC}+{\displaystyle \frac{\mathit{\beta }}{\mathrm{2}}}-\sqrt{{\displaystyle \frac{{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{4}}}+\mathit{\beta }\mathrm{DIC}}\end{array}$$

with

$$\begin{array}{}\text{(A9)}& \mathit{\beta }={\displaystyle \frac{{\left(\mathrm{2}\mathrm{DIC}-\mathrm{Alk}\right)}^{\mathrm{2}}}{\left(\mathrm{Alk}-\mathrm{DIC}\right)}}=p{\mathrm{CO}}_{\mathrm{2}}\cdot {\displaystyle \frac{K_{\mathrm{0}}\cdot K_{\mathrm{1}}}{K_{\mathrm{2}}}}.\end{array}$$

Differentiating Alk with respect to DIC at constant pCO₂, S, and T then yields

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}& {\left.{\displaystyle \frac{\mathrm{dAlk}}{\mathrm{dDIC}}}\right|}_{p{\mathrm{CO}}_{\mathrm{2}},S,T}\simeq \mathrm{2}-{\displaystyle \frac{\mathit{\beta }}{\mathrm{2}\cdot \sqrt{\frac{{\mathit{\beta }}^{\mathrm{2}}}{\mathrm{4}}+\mathit{\beta }\mathrm{DIC}}}}\\ & =\mathrm{2}-{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{1}-\frac{\mathrm{4}\mathrm{DIC}}{\mathit{\beta }}}}}.\end{array}\end{array}$$

After rearrangement, this leads to the same expression for Q_inv:

$$\begin{array}{}\text{(A11)}& {\left.{\displaystyle \frac{\mathrm{dAlk}}{\mathrm{dDIC}}}\right|}_{p{\mathrm{CO}}_{\mathrm{2}},S,T}=Q={\displaystyle \frac{\mathrm{1}}{Q_{\mathrm{inv}}}}\simeq {\displaystyle \frac{\mathrm{3}\mathrm{Alk}-\mathrm{2}\mathrm{DIC}}{\mathrm{Alk}}}.\end{array}$$

An exact formulation of Q is also provided by Humphreys et al. (2018, their Appendix D).

Appendix B: Model and simulations

B1 CMIP6 ESMs and GCB GOBMs

We present an evaluation of the representation of the pre-industrial air–sea carbon flux in ESMs and GOBMs that participated in the CMIP6 exercise (Eyring et al., 2016) and the 2024 GCB exercise (Friedlingstein et al., 2025). This assessment offers valuable insights into the current state of the art regarding the modelling of this flux in the models utilized for intercomparison studies. To ensure comparability, we regridded the CMIP6 data to a regular 1°×1° grid using the distance-weighted average remapping method “remapdis” provided by the Climate Data Operator (CDO). This step was taken because the data available from the 2024 GCB (Hauck et al., 2022) were already on a regular 1°×1° grid. However, it is important to note that this regridding process introduced a minor error in the integrated air–sea carbon flux values.

Figure B1 Towards a controlled adjustment of extra CaCO₃ burial/dissolution. (a) Map depicting CaCO₃ burial in the standard simulation (std). (b, c) Masks employed to drive (b) a balanced (atlpac) or an imbalanced (atlpac-diseq) additional CaCO₃ burial/dissolution between the deep Atlantic and Pacific, as well as (c) an extra CaCO₃ burial in the tropics. Red (blue) shading represents an addition (removal) of DIC and Alk in the grid cell at a 1:2 ratio. The grid cells considered for this addition/removal are located at 4750 m for the deep Atlantic and Pacific masks and between 0 and 100 m for the tropics mask.

We assessed 15 CMIP6 ESMs from 12 different climate modelling centres (Eyring et al., 2016): CanESM5 (r1i1p2f1) and CanESM5-CanOE (r1i1p2f1) from CCCma, with two distinct marine biogeochemical models; CMCC-ESM2 (r1i1p1f1) from CMCC; CNRM-ESM2-1 (r1i1p1f2) from CNRM-CERFACS; ACCESS-ESM1-5 (r1i1p1f1) from CSIRO; IPSL-CM6A-LR (r1i1p1f1) from IPSL; MIROC-ES2L (r1i1p1f2) from MIROC; UKESM1-0-LL (r1i1p1f2) from MOHC; MPI-ESM1-2-LR (r1i1p1f1) and MPI-ESM1-2-HR (r1i1p1f1) from MPI-M, with two different resolutions; MRI-ESM2-0 (r1i2p1f1) from MRI, CESM2-WACCM (r1i1p1f1) from NCAR; NorESM2-LM (r1i1p1f1) from NCC; and GFDL-CM4 (r1i1p1f1) and GFDL-ESM4 (r1i1p1f1) from NOAA-GFDL, with two distinct marine biogeochemical models. Only the air–sea CO₂ flux (positive downward, “fgCO₂” in kgC m⁻² s⁻¹) of the pre-industrial control simulations was considered, from 1850 to 2100, and yearly averaged. Each ESM was weighted in the calculation of the CMIP6 mean such that each modelling group has the same total contribution. We also assessed the 10 GOBMs used in the 2024 GCB exercise (Friedlingstein et al., 2025): NEMO3.6-PISCESv2-gas (CNRM), NEMO4.2-PISCES (IPSL), MPIOM-HAMOCC6, MRI-ESM2-3, ACCESS, MICOM-HAMOCC (NorESM-OC), MOM6-COBALT (Princeton), FESOM-2.1-REcoM3, NEMO-PlanckTOM12, and CESM-ETHZ). Once again, only the air–sea CO₂ flux (positive downward, “fgCO₂” in mol m⁻² s⁻¹) of the control simulations (i.e. constant atmospheric CO₂, no climate change and variability) was considered, from 1959 to 2023, and yearly averaged. We found that the drift in the ESMs and GOBMs in the net air–sea carbon flux was consistently less than 0.10 PgC (100 years)⁻¹, and as such, it had a negligible impact on the related results (see Figs. 2a and D1a).

B2 NEMO-PISCES sensitivity simulations

B2.1 Configurations

Here we provide additional details regarding the various configurations of the sensitivity simulations conducted using NEMO-PISCES (Table B1). In the standard configuration, the slight deviation (−0.02 PgC yr⁻¹) between Alk riverine discharge (+0.35 PgC yr⁻¹) and inorganic carbon riverine discharge (+0.37 PgC yr⁻¹) arises from the supply of inorganic nitrogen by rivers, presumed to be in the form of nitrate, which has a negative impact on Alk (Fig. B4). It is worth noting that the global values and latitudinal distribution of riverine inputs are based on Ludwig et al. (1996) and have recently been revised (Li et al., 2017; Liu et al., 2024), although the human imprint on these fluxes cannot be removed. Lastly, we emphasize that we did not evaluate the implications of partitioning riverine inputs between inorganic and organic components on biological production and, consequently, its effects on the air–sea carbon flux, as we only altered DIC and Alk in the various configurations. Finally, we accounted for atmospheric deposition in our sensitivity simulations, since atmospheric nitrogen deposition is considered a nitrate source, which impacts Alk. This has however a negligible effect, as does the dilution effect (see Fig. S2).

The paper has been crafted to be accessible and comprehensible for both observationalists and modellers. However, the deviations mentioned for the carbon and Alk inventories manifest themselves in model outputs in the form of drifts. Furthermore, all the sensitivity simulations conducted also address modelling issues. In particular, a case that can be encountered in marine biogeochemistry models, both historically and even today, is the consideration, or lack thereof, of the OM and CaCO₃ burial. This can have consequences on the carbon flux, depending on whether the global Alk inventory is conserved through a global-scale Alk restoration scheme or left deviating (nosed-resto, nosed-diseq; Planchat et al., 2023). Finally, the choice of the different configurations and their resulting impact on the air–sea carbon flux also serve as a reminder of the importance of carefully considering the global Alk inventory in models and controlling its potential deviation/drift according to desired hypotheses (e.g. a balanced Alk budget or not).

Figure B2 Spin-up of the NEMO-PISCES sensitivity simulations. (a) Time series of the net air–sea carbon flux with a 50-year rolling mean throughout the 2550 years of the simulations. (b) Same time series in relative to std and without smoothing. The thin black lines refer to the combined exponential fits ($y=\mathit{\alpha }\cdot e^{-\frac{t}{\mathit{\tau }}}+\mathit{\beta }$, where α is the net air–sea carbon flux offset, τ is the time constant, and β is the baseline; using the curve_fit function from the scipy python library): (i) one for the short-term considering the first 50 years and (ii) one for the long-term, considering the remaining 2500 years. (c, d) For the short-term (c) and long-term (d) exponential fits, the net air–sea carbon flux offset (α) is displayed as a function of the time constant (τ) with their associated uncertainties.

From a practical standpoint, in NEMO-PISCES, CaCO₃ burial predominantly occurs in coastal areas (Fig. B1a), with limited differentiation in burial at depth between the Atlantic (less acidic) and Pacific (more acidic) regions (Sarmiento and Gruber, 2006; Cartapanis et al., 2018; Ridgwell and Zeebe, 2005, Fig. B1a). To address this limitation, we introduced the configuration “atlpac” wherein we constrain extra CaCO₃ burial in the deep Atlantic while simulating extra CaCO₃ dissolution in the deep Pacific. This adjustment aims to enhance the representation of CaCO₃ burial while maintaining a balanced Alk budget (i.e. conserving the global Alk inventory; Fig. B1b). Additionally, considering the possibility of an imbalanced Alk budget during the pre-industrial era due to extra CaCO₃ burial at depth (Cartapanis et al., 2018), we created two configurations to account for this extra carbon burial (0.10 PgC yr⁻¹): (i) in the deep Atlantic (atlpac-diseq) and (ii) in the shallow tropical regions (tropics-diseq), simulating the accumulation of carbon by coral reefs (Fig. B1b).

Table B1 Summary of the sensitivity simulations led with NEMO-PISCES with their full description and characteristics (supplement to Table 2).

B2.2 Spin-up

We track here the evolution of the net air–sea carbon flux during the spin-up for all the NEMO-PISCES sensitivity simulations (Table B1 and Fig. B2a), which are initially branched to a simulation at a quasi-steady state equivalent to our standard simulation (std). Two characteristic timescales emerge (Fig. B2b): (i) a short-term stabilization over the first 50 years and (ii) a long-term stabilization beyond 50 years. The short-term (long-term) stabilization primarily corresponds to the response of the surface (deep) ocean to the modifications associated with the configuration regarding the DIC and Alk external fluxes (Fig. B2c, d). Thus, for the simulation where we constrain extra CaCO₃ burial in the shallow tropics, only a stabilization of the surface ocean is generally needed, resulting in only a short-term stabilization. Conversely, in the case where this extra CaCO₃ burial is constrained in the deep Atlantic, only a stabilization of the deep ocean is generally needed, resulting in only a long-term stabilization. Finally, in the case where riverine organic matter input is considered to be entirely refractory (rivref), a significant anomaly in external fluxes is induced at the surface compared to the standard simulation (std), as well as in the deep ocean because this organic carbon input is spread all over the ocean. This results in both short-term and long-term responses.

Figure B3 Description of the standard NEMO-PISCES simulation (std; continued from Fig. 4). Zonally integrated (a) carbon and (b) Alk fluxes in supplement to Fig. 4b. (c) Latitudinal distribution of the northward transport of carbon (solid) and Alk (dashed). When the regional imbalance exceeds (falls behind) the simulated air–sea carbon flux, the area in between is shaded in rose (cyan).

B2.3 Standard simulation (std)

We provide additional details here regarding the standard simulation (std, Fig. B3) to offer points of comparison with historical modelling studies that have initiated research efforts on this pre-industrial carbon flux (Aumont et al., 2001; Murnane et al., 1999).

B2.4 Residual carbon budget imbalance

A minor imbalance in the carbon budget from external sources and sinks persists in our sensitivity simulations, yet without any associated deviation in the ocean carbon content. This residual is particularly evident in the standard simulation (std; see “Total” and “Drift” in Figs. S2 and B4a) but is also visible in several others (e.g. rivref, rivorg, rivinorg, riv1p5, atlpac, atlpac-diseq, and tropics-diseq). Despite extensive checks, the origin of this small residual in the total carbon flux could not be determined with confidence. Possible explanations that remain include (i) a missing component in the OM burial, (ii) a bug in the model code producing an artificial carbon sink within the internal oceanic fluxes, or (iii) a transient enhancement of the oceanic carbon sink, associated with a rearrangement of the organic carbon pools and related to internal variability over the 50-year averaging period. The latter appears most consistent with the model behaviour, as air–sea carbon fluxes fluctuate around their mean value by about ±0.05 PgC yr⁻¹ even under steady-state conditions. Further investigations are therefore required to fully characterize this residual imbalance in the carbon budget, which likely explains the slight offset between the simulated air–sea carbon flux and the flux inferred from the global carbon–Alk imbalance (e.g. +0.04 PgC yr⁻¹ for std; Fig. B4a).

Table B2 Comprehensive description of the net air–sea carbon flux in the NEMO-PISCES sensitivity simulations. The values provided in parentheses are expressed relative to the simulation without riverine and burial fluxes, representing the riverine/burial-driven air–sea carbon flux or carbon transport.

Another form of imbalance in the carbon budget appears in the simulation without an external ocean source or sink (norivbur), which exhibits non-conserved ocean carbon and Alk inventories (see “Total” and “Drift” in Figs. S2 and Fig. B4a). This imbalance most likely originates from the representation of nitrogen reactions in NEMO-PISCES, which includes the restoration of nitrate content in the ocean. An imbalance between nitrification (which decreases Alk) and denitrification (which increases Alk) would generate an internal Alk imbalance. A similar issue related to nitrogen-cycle reactions has also been reported for COBALTv2 (Stock et al., 2020). Since CaCO₃ burial is not considered in this simulation (see Sect. 2.2), this imbalance cannot be compensated by the mechanism that conserves the global Alk inventory. Consequently, at steady state, the resulting positive global Alk inventory deviation (D^Alk) induces an air–sea carbon flux (F^{C, air–sea}) and an ocean carbon content deviation (D^C) of equal magnitude, such that $F^{\text{C,\hspace{0.17em}air–sea}}=D^{\mathrm{C}}=Q_{\mathrm{inv}}\cdot D^{\mathrm{Alk}}$ (see Fig. 3a and Sect. 3.2.3). Thus, the carbon budget imbalance in norivbur arises from an air–sea carbon flux driven by an internal Alk imbalance, which also results in a non-conserved ocean carbon content. As expected, this imbalance is nearly identical in the simulation without burial and with a non-conserved global Alk inventory (nobur–diseq).

Since we were not able to fully characterize or access the spatial distribution of these residual imbalances in the carbon and Alk budgets, we did not attempt to correct them – except in Fig. 5, where only the total air–sea carbon flux value is considered, without spatial distribution. Other very minor residual inconsistencies may persist, such as small deviations in the budgets over the 50-year averaging period. In addition, the modelling scheme of the physical component of NEMO–PISCES induces a slight Alk deviation and a somewhat larger carbon content deviation (respectively −0.002 PgC yr⁻¹ and +0.01 PgC yr⁻¹ in std).

B2.5 NEMO-PISCES sensitivity simulation ensemble

We provide a comprehensive overview of the global-scale carbon and Alk budgets for all NEMO-PISCES sensitivity simulations (Fig. B4). Even more detailed information can be found in Figs. S1 and S2 (https://doi.org/10.5281/zenodo.8421898, Planchat, 2023b). Finally, we also provide a comprehensive characterization of the air–sea carbon flux in the NEMO-PISCES sensitivity simulations, including both the total flux and the riverine/burial-driven component (Table B2).

Figure B4 Global-scale carbon and Alk budgets for all NEMO-PISCES sensitivity simulations. Carbon and Alk budgets (a) in absolute values for the standard simulation (std) and the simulation without riverine and burial fluxes of carbon and Alk (norivbur) or (b) relative to std for the other NEMO-PISCES sensitivity simulations. The type of representation is close to the one shared in Fig. 4c, but integrated over the whole ocean. All fluxes, in petagrams of carbon per year (PgC yr⁻¹) for carbon (in bold) and Alk (in normal font), are directed by arrows, with orientation indicating the sign and size reflecting the absolute magnitude of the flux. In (b), only the fluxes (riverine discharge, as well as OM and CaCO₃ burial) with a significant anomaly are displayed, along with their associated changes relative to the standard simulation (std) in brackets, for both carbon (bold) and Alk (normal font). Additionally, values for carbon and Alk deviations (for simulations with a “-diseq” suffix), net addition flux (for atlpac, atlpac-diseq, and tropics-diseq), or the term of Alk restoration (for nosed-resto) is/are also shown when applicable (Table B1). In (a), for the standard simulation (std), a first approximation of the impact of OM and CaCO₃ production in the surface waters is also inferred from particulate organic carbon (POC) and particulate inorganic carbon (PIC) export at 100 m (in brackets with a star). Finally, in addition to the air–sea carbon flux (dark blue), the air–sea carbon flux stemming from global imbalance (light blue; Sect. 2.1.1 and Appendix B2.4 for an explanation of the residual imbalance) is also shared, as well as the associated inter-hemispheric air–sea carbon flux gradient (dark and light cyan). In (b), as the values are shown relative to the standard simulation (std), the simulated air–sea carbon flux anomalies are equivalent to the ones of the riverine/burial-driven air–sea carbon flux. A detailed description of the NEMO-PISCES sensitivity simulations can be found in Figs. S1 and S2.

Appendix C: Ocean regions, and boundary conditions

The boundaries chosen to demarcate the southern, inter-tropical, and northern regions at 20° S and 20° N (see Figs. 3b, 4b, c, as well as Tables B2, E1, and E2) have indeed been previously employed in the literature (e.g. Aumont et al., 2001; Resplandy et al., 2018). These boundaries primarily align with physical features of the ocean, especially concerning air–sea carbon fluxes. It is in, or very close to, these latitudes that the air–sea carbon flux resulting from regional carbon:Alk budget imbalance reconciles with the simulated one (see Fig. 4c). By employing these boundaries, the air–sea carbon flux from regional carbon:Alk budget imbalances (see Sect. 2.1.2) closely matches the simulated values for each oceanic region (see Fig. S1). This alignment deteriorates when the boundaries are shifted away from 20° S and 20° N. Consequently, we have opted for a consistent approach, maintaining the 20° S/N boundary to delineate distinct oceanic regions, despite the shift to 30° S/N boundaries in the GCB, primarily to correspond with terrestrial biomes (Friedlingstein et al., 2025). However, for potential use in the GCB, we share values of the spatial distribution with boundaries at 30° S/N in Table E2.

Appendix D: Inter-hemispheric air–sea carbon flux gradient

D1 Partitioning between the northern and southern components

We share additional insights regarding the inter-hemispheric air–sea carbon flux gradient, which is crucial for the global carbon cycle in its connection with the atmosphere and land (e.g. Keeling et al., 1989; Resplandy et al., 2018). It is thus valuable to distinguish in this inter-hemispheric gradient the component associated with the net air–sea carbon flux in both the southern and northern regions (Fig. D1).

Figure D1 Decomposition of the inter-hemispheric air–sea carbon flux gradient (supplement to Fig. 2). Decomposition of the inter-hemispheric air–sea carbon flux gradient into the net southern and northern air–sea carbon fluxes for (a) CMIP6 and GCB and (b) the NEMO-PISCES sensitivity simulations. Filled contours correspond to the inter-hemispheric air–sea carbon flux gradient. (a) The 15 CMIP6 ESMs (10 GCB GOBMs) are plotted with red squares (orange circles). The black square and circle refer to the CMIP6 and GCB ensemble means.

Figure D2 Drivers of the spatial distribution of the riverine/burial-driven air–sea carbon flux. Inter-hemispheric gradient of the riverine/burial-driven air–sea carbon flux (filled contours) and its two components, from carbon:Alk budget imbalances (see Sect. 2.1.2). One component (x axis) is associated with the inter-hemispheric gradient of air–sea carbon flux driven by northern and southern carbon:Alk budget imbalances (and inventory deviations), while the other component (y axis) corresponds to the inter-hemispheric gradient of air–sea carbon flux associated with the inter-hemispheric transport of the carbon:Alk budget imbalance (Eq. D7 and Fig. D1). The deviation in the simulated inter-hemispheric gradient in NEMO-PISCES sensitivity simulations as compared to the reconstructed ones using the two components is minimal ($<\mathrm{0.01}\mathrm{PgC}{\mathrm{yr}}^{-\mathrm{1}}$, not shown).

D2 Partitioning between carbon transport and riverine and burial processes

D2.1 Expression

In addition to defining the inter-hemispheric air–sea carbon flux gradient (G) and the inter-hemispheric oceanic transports of carbon (T^C) and Alk (T^Alk), here we propose a decomposition of G into contributions associated with carbon transport and with riverine and burial processes.

To this end, we recall that the total regional fluxes of carbon and Alk can be expressed as follows:

$$\begin{array}{}\text{(D1)}& \left\{\begin{array}{l}{F}_{N/S/I}^{\mathrm{C}}=F_{N/S/I}^{\text{C, air–sea}}+F_{N/S/I}^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}\\ F_{N/S/I}^{\mathrm{Alk}}=F_{N/S/I}^{\mathrm{Alk},\mathrm{riv}./\mathrm{bur}.}\end{array}\right.\end{array}$$

Specifically, by considering these fluxes ($F_{N,S,I}^{\mathrm{C}/\mathrm{Alk}}$) along with regional carbon and Alk deviations ($D_{N,S,I}^{\mathrm{C}/\mathrm{Alk}}$) and assuming a steady-state ocean (see Fig. B2), we derive two expressions for the ocean transport of carbon and Alk through their respective budget closure equations ($T_{N/S}^{\mathrm{C}/\mathrm{Alk}}$):

$$\begin{array}{}\text{(D2)}& \left\{\begin{array}{l}{T}_{N/S}^{\mathrm{C}/\mathrm{Alk}}+F_{\ge \mathrm{20}\mathit{°}\mathrm{N}/\mathrm{20}\mathit{°}\mathrm{S}}^{\mathrm{C}/\mathrm{Alk}}+D_{\ge \mathrm{20}\mathit{°}\mathrm{N}/\mathrm{20}\mathit{°}\mathrm{S}}^{\mathrm{C}/\mathrm{Alk}}=\mathrm{0}\\ -T_{N/S}^{\mathrm{C}/\mathrm{Alk}}+F_{<\mathrm{20}\mathit{°}\mathrm{N}/\mathrm{20}\mathit{°}\mathrm{S}}^{\mathrm{C}/\mathrm{Alk}}+D_{<\mathrm{20}\mathit{°}\mathrm{N}/\mathrm{20}\mathit{°}\mathrm{S}}^{\mathrm{C}/\mathrm{Alk}}=\mathrm{0}\end{array}\right..\end{array}$$

Hence, we define the ocean transport of carbon and Alk as the average of its two expressions (Eq. D2; Fig. 3b):

$$\begin{array}{}\text{(D3)}& \begin{array}{rl}{T}_{N/S}^{\mathrm{C}/\mathrm{Alk}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left[\left(F_{<\mathrm{20}\mathit{°}N/\mathrm{20}\mathit{°}S}^{\mathrm{C}/\mathrm{Alk}}+D_{<\mathrm{20}\mathit{°}N/\mathrm{20}\mathit{°}S}^{\mathrm{C}/\mathrm{Alk}}\right)\right.\\ & \left.-\left(F_{\ge \mathrm{20}\mathit{°}N/\mathrm{20}\mathit{°}S}^{\mathrm{C}/\mathrm{Alk}}+D_{\ge \mathrm{20}\mathit{°}N/\mathrm{20}\mathit{°}S}^{\mathrm{C}/\mathrm{Alk}}\right)\right],\end{array}\end{array}$$

in particular

$$\begin{array}{}\text{(D4)}& \left\{\begin{array}{l}\begin{array}{rl}{T}_N^{\mathrm{C}/\mathrm{Alk}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left[\left(F_S^{\mathrm{C}/\mathrm{Alk}}+D_S^{\mathrm{C}/\mathrm{Alk}}+F_I^{\mathrm{C}/\mathrm{Alk}}+D_I^{\mathrm{C}/\mathrm{Alk}}\right)\right.\\ & \left.-\left(F_N^{\mathrm{C}/\mathrm{Alk}}+D_N^{\mathrm{C}/\mathrm{Alk}}\right)\right]\end{array}\\ \begin{array}{rl}{T}_S^{\mathrm{C}/\mathrm{Alk}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left[\left(F_S^{\mathrm{C}/\mathrm{Alk}}+D_S^{\mathrm{C}/\mathrm{Alk}}\right)\right.\\ & \left.-\left(F_I^{\mathrm{C}/\mathrm{Alk}}+D_I^{\mathrm{C}/\mathrm{Alk}}+F_N^{\mathrm{C}/\mathrm{Alk}}+D_N^{\mathrm{C}/\mathrm{Alk}}\right)\right],\end{array}\end{array}\right.\end{array}$$

from which an expression of the inter-hemispheric transport of carbon and Alk (Eq. 12) can be derived:

$$\begin{array}{}\text{(D5)}& T^{\mathrm{C}/\mathrm{Alk}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left[\left(F_S^{\mathrm{C}/\mathrm{Alk}}+D_S^{\mathrm{C}/\mathrm{Alk}}\right)-\left(F_N^{\mathrm{C}/\mathrm{Alk}}+D_N^{\mathrm{C}/\mathrm{Alk}}\right)\right].\end{array}$$

Specifically, using Eq. (D1), the inter-hemispheric transport of carbon can be rewritten as follows:

$$\begin{array}{}\text{(D6)}& \begin{array}{rl}{T}^{\mathrm{C}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left[\left(F_S^{\text{C, air–sea}}+F_S^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}+D_S^{\mathrm{C}}\right)\right.\\ & \left.-\left(F_N^{\text{C, air–sea}}+F_N^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}+D_N^{\mathrm{C}}\right)\right].\end{array}\end{array}$$

After rearrangement, the inter-hemispheric air–sea carbon flux gradient can be expressed as

$$\begin{array}{}\text{(D7)}& G=\underset{\mathrm{Transport}\mathrm{component}}{\underbrace{-\mathrm{2}\cdot T^{\mathrm{C}}}}\underset{\mathrm{Riverine}\mathrm{and}\mathrm{burial}\mathrm{component}}{\underbrace{-\left(F_N^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}-F_S^{\mathrm{C},\mathrm{riv}./\mathrm{bur}.}\right)-\left(D_N^{\mathrm{C}}-D_S^{\mathrm{C}}\right)}}.\end{array}$$

At first glance, it may appear that this expression is exclusively formulated in terms of carbon, seemingly without any consideration of Alk. However, Alk plays a subtle yet integral role in this equation: firstly, because T^C depends on both the southern and northern air–sea carbon fluxes (Eqs. D5 and D1), and these regional fluxes are chemically driven by the relative imbalance between Alk and DIC. Secondly, the deviations in the carbon content of the northern and southern oceans ($D_N^{\mathrm{C}}$ and $D_S^{\mathrm{C}}$, respectively) are directly linked to the deviations in Alk content (Eq. 10). Thus, the role of Alk is intricately interwoven within the formulation of G (Eq. D7).

D2.2 Results

It is possible to decompose the inter-hemispheric gradient of the riverine/burial-driven air–sea flux into a component associated with the inter-hemispheric transport of a carbon:Alk budget imbalance and a component associated with a carbon:Alk budget imbalance stemming from riverine and burial fluxes (including inventory deviations; Eq. (D7) and Figs. D2 and D3. Focusing solely on the effect of riverine and burial fluxes (i.e. relative to norivbur), the component associated with these external fluxes (−0.07 PgC yr⁻¹) is offset by the transport-related component (+0.07 PgC yr⁻¹) in the standard configuration, resulting in a null inter-hemispheric riverine/burial-driven air–sea carbon flux gradient for std (Figs. B4a and D2). In the case of a surplus of CaCO₃ burial in the deep Atlantic (atlpac-diseq), the increase in the inter-hemispheric riverine/burial-driven air–sea carbon flux gradient relative to std (+0.07 PgC yr⁻¹; see Fig. B4b) is primarily attributed to the transport of a carbon:Alk budget imbalance (Fig. D2). Conversely, when the riverine discharge is entirely organic (rivorg), it is mostly the external flux component that causes the decrease in the inter-hemispheric air–sea carbon flux gradient relative to std (−0.18 PgC yr⁻¹; see Figs. B4b and D2), and the same outcome occurs when the riverine discharge is entirely inorganic (rivinorg). This emphasizes that the spatial distribution of the carbon:Alk budget imbalance stemming from external fluxes, in conjunction with oceanic transport, plays a significant role in shaping the pre-industrial inter-hemispheric riverine/burial-driven air–sea carbon flux gradient.

Figure D3 Partitioning of the inter-hemispheric air–sea carbon flux gradient from inter-hemispheric carbon:Alk budget imbalances (supplement to D2). (a, b) Visual construction of the inter-hemispheric air–sea carbon flux gradient resulting from regional carbon:Alk budget imbalances due to (a) riverine and burial fluxes (including inventory deviations) and (b) inter-hemispheric transport, defining the values used in Fig. D2. The reference was set on the simulation without riverine and burial fluxes (norivbur), so that the combination of the arrows of (a) and (b) results in the inter-hemispheric riverine/burial-driven air–sea carbon flux gradient. (c) Synthetic characterization for the whole set of NEMO-PISCES sensitivity simulations of the inter-hemispheric air–sea carbon flux gradient: with a southern/northern decomposition (as in Fig. D1b) and the partitioning resulting from regional carbon:Alk budget imbalances due to riverine and burial fluxes (including inventory deviations) constructed in (a) and inter-hemispheric transport constructed in (b).

Appendix E: Applications

E1 Literature review

Here, we provide a literature review on (i) the evolution of the assessment and characterization of the air–sea carbon flux since the late 1990s (Table E1) and (ii) the evolution of the estimation and characterization of the riverine/burial-driven air–sea carbon flux in comparison with our composite simulated estimate (Table E2). We also provide the PDFs of the literature-based estimates for the ocean's external sources/sinks of carbon and Alk, derived from the most recent carbon and Alk budgets (Regnier et al., 2022; Middelburg et al., 2020), which were used to construct the composite simulated estimate (Fig. E1; see Table 3 as well as Sect. 4.2.1 and 4.2.2).

Figure E1 Literature-based estimates of the riverine discharge^a, OM burial^b, and CaCO₃ burial, with their associated uncertainties/extremes through normalized PDFs. ^a Including groundwater discharge for both carbon and Alk and anaerobic processes for Alk. ^b Including reverse weathering for Alk.

(Sarmiento et al., 1992)Tans et al. (1990)Jacobson et al. (2007)Jacobson et al. (2007)(Khatiwala et al., 2009)Jacobson et al. (2007)Takahashi et al. (2009)Wanninkhof et al. (2013)Landschützer et al. (2014)Rödenbeck et al. (2013)(Landschützer et al., 2014)(DeVries et al., 2019)Jacobson et al. (2007)

Table E1 Literature review of the net air–sea carbon flux and its characterization. Bold lines are those accounting for external fluxes and boundaries at 20° S/N for the spatial distribution, and the line shown in italics corresponds to the present study.

^b Boundaries at 15° S/N. ^d Boundaries at 13° S/N. ^e The values were extracted from the supplement figure of Mikaloff Fletcher et al. (2007). ^f Values from Resplandy et al. (2018). ^g Accounting for all the pCO₂-based estimates, except the one by Wanninkhof et al. (2013), which is incomplete.

Regnier et al. (2022)Lacroix et al. (2020)Friedlingstein et al. (2023, 2025)Aumont et al. (2001)Friedlingstein et al. (2022b)Murnane et al. (1999)Resplandy et al. (2018)Resplandy et al. (2018)Aumont et al. (2001)Jacobson et al. (2007)Jacobson et al. (2007)Jacobson et al. (2007)Friedlingstein et al. (2020, 2022a)Lacroix et al. (2020)Resplandy et al. (2018)Jacobson et al. (2007)Le Quéré et al. (2018b)Friedlingstein et al. (2019)Jacobson et al. (2007)Resplandy et al. (2018)Jacobson et al. (2007)Resplandy et al. (2018)Aumont et al. (2001)

Table E2 Literature review of the riverine/burial-driven air–sea carbon flux estimates, sharing also values and a description of the riverine discharge, OM burial, and CaCO₃ burial. The values obtained from the composite simulated estimate are also shared in italics. Note that its riverine discharge, OM burial and CaCO₃ fluxes (see Table 3) differ from those considered in the GCB and derived from Regnier et al. (2022), as we also account for flux on continental shelves.

^a Rounded at 0.05 PgC yr⁻¹ in their paper. ^b Imbalance due to the consideration of a disturbance in the biological pump (production minus burial: +0.45 PgC yr⁻¹) and the consideration of a free global Alk inventory. ^c Assumption that carbon from riverine discharge is not transported and that part of it is buried in the same region as well. ^d e.g. Le Quéré et al. (2016, 2018a).

E2 Reassessing the regional distribution of the riverine/burial-driven air–sea carbon flux

We provide a schematic of the practical framework introduced in this paper (Sect. 4.2.1), which follows a four-step workflow to construct a composite simulation aligning with the most recent carbon and Alk budgets (Fig. E2).

In addition, we also share (Fig. E3) the various components of the composite simulated estimate creation process as described in Sect. 4.2.1, the results of which are presented in Sect. 4.2.2 (see Fig. 7 and Table E2).

Figure E2 Schematic of the four-step workflow to construct the composite simulation. The composite simulated estimate of the riverine/burial-driven air–sea carbon flux is built from a linear combination of our sensitivity simulations to match the most recent carbon and Alk budgets (riverine discharge, OM burial, and CaCO₃ burial). An additional step, equivalent to Step 3 but not shown in this general workflow, was necessary to adjust the CaCO₃ burial of carbon in excess (Regnier et al., 2022), compared to the values accounted for Alk (Middelburg et al., 2020). In the schematic, the exponents indicate the source of each variable/parameter: E denotes the composite simulated estimate, S refers to the various sensitivity simulations, and L stands for literature estimates.

Figure E3 Components of the composite simulated estimate (supplement to Fig. 7). Each of the components represents the elements added at the different stages of the composite simulated estimate construction process (Steps 1, 2, 3, 4, and the extra step; Fig. E2). (a, b) Decomposition of the composite simulated estimate PDF associated with the total (a) carbon and (b) Alk external fluxes. (c, d, e, f) Characterization of the riverine/burial air–sea carbon flux in the composite simulated estimate, showing the various components for (c) the total value, as well as the (d) southern, (e) inter-tropical, and (f) northern regions. The black solid lines represent the total values for the composite simulated estimate, while the black dotted lines (a, b) correspond to the total carbon and Alk external fluxes from literature estimates (Fig. E1).