Historical and future contributions of inland waters to the Congo Basin carbon balance

As the second largest area of contiguous tropical rainforest and second largest river basin in the world, the Congo Basin has a significant role to play in the global carbon (C) cycle. For the present day, it has been shown that a significant proportion of global terrestrial net primary productivity (NPP) is transferred laterally to the land–ocean aquatic continuum (LOAC) as dissolved CO₂, dissolved organic carbon (DOC), and particulate organic carbon (POC). Whilst the importance of LOAC fluxes in the Congo Basin has been demonstrated for the present day, it is not known to what extent these fluxes have been perturbed historically, how they are likely to change under future climate change and land use scenarios, and in turn what impact these changes might have on the overall C cycle of the basin. Here we apply the ORCHILEAK model to the Congo Basin and estimate that 4 % of terrestrial NPP (NPP = 5800±166 Tg C yr⁻¹) is currently exported from soils and vegetation to inland waters. Further, our results suggest that aquatic C fluxes may have undergone considerable perturbation since 1861 to the present day, with aquatic CO₂ evasion and C export to the coast increasing by 26 % (186±41 to 235±54 Tg C yr⁻¹) and 25 % (12±3 to 15±4 Tg C yr⁻¹), respectively, largely because of rising atmospheric CO₂ concentrations. Moreover, under climate scenario RCP6.0 we predict that this perturbation could continue; over the full simulation period (1861–2099), we estimate that aquatic CO₂ evasion and C export to the coast could increase by 79 % and 67 %, respectively. Finally, we show that the proportion of terrestrial NPP lost to the LOAC could increase from approximately 3 % to 5 % from 1861–2099 as a result of increasing atmospheric CO₂ concentrations and climate change. However, our future projections of the Congo Basin C fluxes in particular need to be interpreted with some caution due to model limitations. We discuss these limitations, including the wider challenges associated with applying the current generation of land surface models which ignore nutrient dynamics to make future projections of the tropical C cycle, along with potential next steps.

1 Introduction

As the world's second largest area of contiguous tropical rainforest and second largest river, the Congo Basin has a significant role to play in the global carbon (C) cycle. Current estimates of its C stocks and fluxes are limited by a paucity of field data and therefore have substantial uncertainties, both quantified and unquantified (Williams et al., 2007; Lewis et al., 2009; Dargie et al., 2017). Nevertheless, it has been estimated that there is approximately 50 Pg C stored in its above-ground biomass (Verhegghen et al., 2012) and up to 100 Pg C contained within its soils (Williams et al., 2007). Moreover, a recent study estimated that around 30 (6.3–46.8) Pg C is stored in the peats of the Congo alone (Dargie at al., 2017). Field data suggest that storage in tree biomass increased by 0.34 (0.15–0.43) Pg C yr⁻¹ in intact African tropical forests between 1968 and 2007 (Lewis et al., 2009) due in large part to a combination of increasing atmospheric CO₂ concentrations and climate change (Ciais et al., 2009; Pan et al., 2015), while satellite data indicate that terrestrial net primary productivity (NPP) has increased by an average of 10 g C m⁻² yr⁻¹ per year between 2001 and 2013 in tropical Africa (Yin et al., 2017).

At the same time, forest degradation, clearing for rotational agriculture, and logging are causing C losses to the atmosphere (Zhuravleva et al., 2013; Tyukavina et al., 2018), while droughts have reduced vegetation greenness and water storage over the last decade (Zhou et al., 2014). A recent estimate of above-ground C stocks in tropical African forests, mainly in the Congo, indicates a minor net C loss from 2010 to 2017 (Fan et al., 2019). Moreover, recent field data suggest that the above-ground C sink in tropical Africa was relatively stable from 1985 to 2015 (Hubau et al., 2020).

There are large uncertainties associated with projecting future trends in the Congo Basin terrestrial C cycle, firstly related to predicting which trajectories of future CO₂ levels and land use changes will occur and secondly to our ability to fully understand and simulate these changes and in turn their impacts. Future model projections for the 21st century agree that temperature will significantly increase under both low and high emission scenarios (Haensler et al., 2013), while precipitation is only projected to substantially increase under high emission scenarios, with the basin mean remaining more or less unchanged under low emission scenarios (Haensler et al., 2013). Uncertainties in future land use change projections for Africa are among the highest for any continent (Hurtt et al., 2011).

For the present day at the global scale, it has been estimated that between 1 and 5 Pg C yr⁻¹ is transferred laterally to the land–ocean aquatic continuum (LOAC) as dissolved CO₂, dissolved organic carbon (DOC), and particulate organic carbon (POC) (Cole at al., 2007; Tranvik et al., 2009; Regnier et al., 2013; Drake et al., 2018; Ciais et al., 2020). This C can subsequently be evaded back to the atmosphere as CO₂, undergo sedimentation in wetlands and inland waters, or be transported to estuaries or the coast. The tropical region is a hotspot area for inland water C cycling (Richey et al., 2002; Melack et al., 2004; Abril et al., 2014; Borges et al., 2015a; Lauerwald et al., 2015) due to high terrestrial NPP and precipitation, and a recent study used an upscaling approach based on observations to estimate present-day CO₂ evasion from the rivers of the Congo Basin at 251±46 Tg C yr⁻¹ and the lateral C (TOC + DIC) export to the coast at 15.5 (13–18) Tg C yr⁻¹ (Borges at al., 2015a, 2019). To put this into context, the estimate of aquatic CO₂ evasion represents 39 % of the global value estimated by Lauerwald et al. (2015, 650 Tg C yr⁻¹) or 14 % of the global estimate of Raymond et al. (2013, 1800 Tg C yr⁻¹). Note that while Lauerwald et al. (2015) and Raymond et al. (2013) relied largely on the same database of partial pressure of CO₂ (pCO₂) measurements (GloRiCh; Hartmann et al., 2014) as the basis for their estimates, they took different, albeit both empirically led, approaches. Moreover, both approaches were limited by a relative paucity of data from the tropics, which also explains the high degree of uncertainty associated with our understanding of global riverine CO₂ evasion.

Whilst the importance of LOAC fluxes in the Congo Basin has been demonstrated for the present day, it is not known to what extent these fluxes have been perturbed historically, how they are likely to change under future climate change and land use scenarios, and in turn what impact these changes might have on the overall C balance of the Congo. In light of these knowledge gaps, we address the following research questions.

• What is the relative contribution of LOAC fluxes (CO₂ evasion and C export to the coast) to the present-day C balance of the basin?

• To what extent have LOAC fluxes changed from 1861 to the present day, and what are the primary drivers of this change?

• How will these fluxes change under future climate and land use change scenarios (RCP6.0, which represents the “no mitigation” scenario), and what are the limitations associated with these future projections?

Understanding and quantifying these long-term changes requires a complex and integrated mass-conservation modelling approach. The ORCHILEAK model (Lauerwald et al., 2017), a new version of the land surface model ORCHIDEE (Krinner et al., 2005), is capable of simulating observed terrestrial and aquatic C fluxes in a consistent manner for the present day in the Amazon (Lauerwald et al., 2017) and Lena (Bowring et al., 2019, 2020) basins, albeit with limitations including a lack of explicit representation of POC fluxes and in-stream autotrophic production (see Lauerwald et al., 2017; Bowring et al., 2019, 2020; and Hastie et al., 2019, for further discussion). Moreover, it was recently demonstrated that this model could recreate observed seasonal and interannual variation in Amazon Basin aquatic and terrestrial C fluxes (Hastie et al., 2019).

In order to accurately simulate aquatic C fluxes, it is crucial to provide a realistic representation of the hydrological dynamics of the Congo River, including its wetlands. Here, we develop new wetland forcing files for the ORCHILEAK model from the high-resolution dataset of Gumbricht et al. (2017) and apply the model to the Congo Basin. After validating the model against observations of discharge, flooded area, DOC concentrations, and pCO₂ for the present day, we then use the model to understand and quantify the long-term (1861–2099) temporal trends in both the terrestrial and aquatic C fluxes of the Congo Basin.

2 Methods

ORCHILEAK (Lauerwald et al., 2017) is a branch of the ORCHIDEE land surface model (LSM), building on past model developments such as ORCHIDEE-SOM (Camino-Serrano et al., 2018), and represents one of the first LSM-based approaches which fully integrates the aquatic C cycle within the terrestrial domain. ORCHILEAK simulates DOC production in the canopy and soils, the leaching of dissolved CO₂ and DOC to the river from the soil, the mineralization of DOC, and in turn the evasion of CO₂ to the atmosphere from the water surface. Moreover, it represents the transfer of C between litter, soils, and water within floodplains and swamps (see Sect. 2.2). Once within the river-routing scheme, ORCHILEAK assumes that the lateral transfer of CO₂ and DOC is proportional to the volume of water. DOC is divided into a refractory and labile pool within the river, with half-lives of 80 and 2 d, respectively. The refractory pool corresponds to the combined slow and passive DOC pools of the soil C scheme, and the labile pool corresponds to the active soil pool (see Sect. 2.4.1). The concentration of dissolved CO₂ and the temperature-dependent solubility of CO₂ are used to calculate pCO₂ in the water column. In turn, CO₂ evasion is calculated based on pCO₂, along with a diurnally variable water surface area and a gas exchange velocity. Fixed gas exchange velocities of 3.5 and 0.65 m d⁻¹, respectively, are used for rivers (including open floodplains) and forested floodplains.

In this study, as in previous studies (Lauerwald et al., 2017; Hastie et al., 2019; Bowring et al., 2019, 2020), we run the model at a spatial resolution of 1^∘ and use the default time step of 30 min for all vertical transfers of water, energy, and C between vegetation, soil, and the atmosphere, as well as the daily time step for the lateral routing of water. Until now, in the tropics, ORCHILEAK has been parameterized and calibrated only for the Amazon River Basin (Lauerwald et al., 2017; Hastie et al., 2019). To adapt and apply ORCHILEAK to the specific characteristics of the Congo River Basin (Sect. 2.1), we had to establish new forcing files representing the maximal fraction of floodplains (MFF) and the maximal fraction of swamps (MFS) (Sect. 2.2) and to recalibrate the river-routing module of ORCHILEAK (Sect. 2.3). All of the processes represented in ORCHILEAK remain identical to those previously represented for the Amazon ORCHILEAK (Lauerwald et al., 2017; Hastie et al., 2019). In the following methodology sections, we describe the following: Sect. 2.1 – Congo Basin description; Sect. 2.2 – development of floodplain and swamp forcing files; Sect. 2.3 – calibration of hydrology; Sect. 2.4 – simulation set-up; Sect. 2.5 – evaluation and analysis of simulated fluvial C fluxes; and Sect. 2.6 – calculating the net carbon balance of the Congo Basin. For a full description of the ORCHILEAK model, please see Lauerwald et al. (2017).

2.1 Congo Basin description

The Congo Basin is the world's second largest area of contiguous tropical rainforest and second largest river basin (Fig. 1), covering an area of 3.7×10⁶ km², with a mean discharge of around 42 000 m⁻³ s⁻¹ (O'Loughlin et al., 2013) and a variation between 24 700 and 75 500 m⁻³ s⁻¹ across months (Coynel et al., 2005).

Figure 1 Extent of the Congo Basin, central quadrant of the Cuvette Centrale, and locations of sampling stations used for validation (of DOC, discharge, and partial pressure of CO₂) along the Congo and Oubangui River (in italics).

Figure 2 Present-day (1981–2010) spatial distribution of the principal climate and land use drivers used in ORCHILEAK across the Congo Basin; (a) mean annual temperature (^∘C), (b) mean annual rainfall (mm yr⁻¹), (c–h) mean annual maximum vegetated fraction for PFTs 2, 3, 10, 11, 12, and 13, respectively; (i) river area; and (j) poor soils. All are at a resolution of 1^∘ except for river area (0.5^∘).

The major climate (ISMSIP2b; Frieler et al., 2017; Lang et al., 2017) and land cover (LUH-CMIP5) characteristics of the Congo Basin for the present day (1981–2010) are shown in Fig. 2. The mean annual temperature is 25.2 ^∘C but with considerable spatial variation from a low of 18.4 ^∘C to a high of 27.2 ^∘C (Fig. 2a), while mean annual rainfall is 1520 mm, varying from 733 to 4087 mm (Fig. 2b). ORCHILEAK prescribes 13 different plant functional types (PFTs). Land use is mixed with tropical broadleaved evergreen (PFT2, Fig. 1c), tropical broadleaved rain green (PFT3, Fig. 1d), C₃ grass (PFT10, Fig. 2e), and C₄ grass (PFT11, Fig. 2f) covering a maximum of 26 %, 35 %, 8 %, and 25 % of the basin area, respectively. Most published estimates for land cover follow national boundaries, so we can make broad comparisons with published estimates for the Democratic Republic of Congo (DRC). For example, our value for total forest cover for the DRC (65 %) is close to the 67 % and 68 % values estimated by the Congo Basin Forest Partnership (CBFP, 2009) and Potapov et al. (2012), respectively. Agriculture covers only a small proportion of the basin according to the LUH dataset that is based on FAO cropland area statistics, with C₃ (PFT12, Fig. 2g) and C₄ (PFT13, Fig. 2h) agriculture making up a maximum basin area of 0.5 % and 2 %, respectively. In reality, a larger fraction of the basin is composed of small-scale and rotational agriculture (Tyukavina et al., 2018). The ORCHILEAK model also has a “poor soils” forcing file (Fig. 2j) which prescribes reduced decomposition rates in soils with low-nutrient and low-pH soils such as Podzols and Arenosols (Lauerwald et al., 2017). This file is developed from the Harmonized World Soil Database (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2009).

2.2 Development of floodplain and swamp forcing files

In ORCHILEAK, water in the river network can be diverted to two types of wetlands: floodplains and swamps. In each grid wherein a floodplain exists, a temporary water body can be formed adjacent to the river and is fed by the river once bankfull discharge (see Sect. 2.3) is exceeded. In grids wherein swamps exist, a constant proportion of river discharge is fed into the base of the soil column; ORCHILEAK does not explicitly represent a groundwater reservoir, so this imitates the hydrological coupling of swamps and rivers through the groundwater table. The maximal proportions of each grid which can be covered by floodplains and swamps are prescribed by the maximal fraction of floodplain (MFF) and the maximal fraction of swamp (MFS) forcing files, respectively (Guimberteau et al., 2012). See also Lauerwald et al. (2017) and Hastie et al. (2019) for further details. We created an MFF forcing file for the Congo Basin derived from the Global Wetlands^v3 database, which includes the 232 m resolution tropical wetland map of Gumbricht et al. (2017) (Fig. 3a and b). We firstly amalgamated all the categories of wetland (which include floodplains and swamps) before aggregating them to a resolution of 0.5^∘ (the resolution at which the floodplain and swamp forcing files are read by ORCHILEAK), assuming that this represents the maximum extent of inundation in the basin. This results in a mean MFF of 10 %; i.e. a maximum of 10 % of the surface area of the Congo Basin can be inundated with water. This is identical to the mean MFF value of 10 % produced with the Global Lakes and Wetlands Database, GLWD (Lehner and Döll, 2004; Borges et al., 2015b). We also created an MFS forcing file from the same dataset (Fig. 3c and d), merging the “swamps” and “fens” wetland categories (although note that there are virtually no fens in the Congo Basin) from the Global Wetlands^v3 database (Gumbricht et al., 2017) and again aggregating them to a 0.5^∘ resolution. Please see Table 1 of Gumbricht et al. (2017) for further details.

Figure 3 (a) Wetland extent (from Gumbricht et al., 2017). (b) The new maximal fraction of floodplain (MFF) forcing file developed from (a). (c) The swamp (including fens) category within the Congo Basin from Gumbricht et al. (2017). (d) The new maximal fraction of swamp (MFS) forcing file developed from (c). Panels (a) and (b) are at the same resolution as the Gumbricht dataset (232 m), while panels (b) and (d) are at a resolution of 0.5^∘. Note that 0.5^∘ is the resolution of the subunit basins in ORCHILEAK (Lauerwald et al., 2015), with each 1^∘ grid containing four sub-basins.

2.3 Calibration of hydrology

As the main driver of the export of C from the terrestrial to the aquatic system, it is crucial that the model can represent present-day hydrological dynamics, at the very least on the main stem of the Congo. As this study is primarily concerned with decadal to centennial timescales our priority was to ensure that the model can accurately recreate observed mean annual discharge at the most downstream gauging station at Brazzaville. We also tested the model's ability to simulate observed discharge seasonality and flood dynamics. Moreover, no data are available with which to directly evaluate the simulation of DOC and CO₂ leaching from the soil to the river network, and thus we tested the model's ability to recreate the spatial variation of observed riverine DOC concentrations and pCO₂ at specific stations where measurements are available (Borges at al., 2015b; Bouillon et al., 2012, 2014; locations shown in Fig. 1), with the river DOC and CO₂ concentration being regarded as an integrator of the C transport at the terrestrial–aquatic interface.

We first ran the model for the present-day period, defined as 1990 to 2005–2010 depending on which climate forcing data were applied, using four climate forcing datasets, namely ISIMIP2b (Frieler et al., 2017), Princeton GPCC (Sheffield et al., 2006), GSWP3 (Kim, 2017), and CRUNCEP (Viovy, 2018). We used ISIMIP2b for the historical and future simulations as it is the only climate forcing dataset to cover the full period (1861–2099). However, we compared it to other climate forcing datasets for the present day in order to gauge its ability to simulate observed discharge on the Congo River at Brazzaville (Table A1). Without calibration, the majority of the different climate forcing model runs performed poorly and were unable to accurately represent the seasonality and mean monthly discharge at Brazzaville (Table A1). The best-performing climate forcing dataset was ISIMIP2b, followed by Princeton GPCC, with root mean square errors (RMSEs) of 29 % and 40 % and a Nash–Sutcliffe efficiency (NSE) of 0.20 and −0.25, respectively. NSE is a statistical coefficient specifically used to test the predictive skill of hydrological models (Nash and Sutcliffe, 1970).

For ISIMIP2b we further calibrated key hydrological model parameters, namely the constants (tau, τ) which help to control the water residence time of the groundwater (slow reservoir), headwaters (fast reservoir), and floodplain reservoirs in order to improve the simulation of observed discharge at Brazzaville and Oubangui (Table 2). To do so, we tested different combinations of τ values for the three reservoirs, eventually settling on 1, 0.5, and 0.5 (d) for the slow, fast, and floodplain reservoirs, respectively, with all three being reduced compared to those values used in the original ORCHILEAK calibration for the Amazon (Lauerwald at al., 2017). The actual residence time of each reservoir is calculated at each time step. The residence time of the flooded reservoir, for example, is a product of τ_flood, a topographical index, and the flooded fraction of the grid cell.

In order to calibrate the simulated discharge against observations, we first modified the flood dynamics of ORCHILEAK in the Congo Basin for the present day by adjusting bankfull discharge (streamr_50th; Lauerwald et al., 2017) and the 95th percentile of water level heights (floodh_95th). As in previous studies on the Amazon Basin (Lauerwald et al., 2017; Hastie et al., 2019) we defined bankfull discharge, i.e. the threshold discharge at which floodplain inundation starts (i.e. overtopping of banks), as the median discharge (50th percentile i.e. streamr_50th) of the present-day climate forcing period (1990 to 2005). After re-running each model parametrization (different τ values) to obtain those bankfull discharge values, we calculated floodh_95th over the simulation period for each grid cell (Table 1). This value is assumed to represent the water level over the riverbanks at which the maximum horizontal extent of floodplain inundation is reached. We then ran the model for a final time and validated the outputs against discharge data at Brazzaville (Cochonneau et al., 2006; Fig. 1). This procedure was repeated iteratively with the ISIMIP2b climate forcing with modifying the τ value of each reservoir in order to find the best-performing parametrization.

Table 1 Main forcing files used for simulations.

Table 2 Performance statistics for the modelled versus observed seasonality of discharge at Brazzaville and Bangui, as well as flooded area in Cuvette Centrale. Observed flooded area is from GIEMS (Papa et al., 2010; Becker et al., 2018).

We firstly compared simulated versus observed discharge at Brazzaville (NSE, RMSE; Table 2), before using the data of Bouillon et al. (2014) to further validate discharge at Bangui (Fig. 1) on the main tributary Oubangui. In addition, we compared the simulated seasonality of flooded area against the satellite-derived dataset GIEMS (Prigent et al., 2007; Becker et al., 2018) within the Cuvette Centrale wetlands (Fig. 1).

2.4 Simulation set-up

A list of the main forcing files used, along with data sources, is presented in Table 1. The derivation of the floodplains and swamps (MFF and MFS) is described in Sect. 2.2, while the calculation of bankfull discharge (streamr_50th) and the 95th percentile of water table height over floodplain (floodh_95th) (Table 1) is described in Sect. 2.3.

2.4.1 Soil carbon spin-up

ORCHILEAK includes a soil module, primarily derived from ORCHIDEE-SOM (Camino-Serrano et al., 2018). The soil module has three different pools of soil DOC: the passive, slow, and active pool, and these are defined by their source material and residence times (τ_carbon). ORCHILEAK also differentiates between flooded and non-flooded soils as well as the decomposition rates of DOC, soil organic carbon (SOC), and litter being reduced (3 times lower) in flooded soils. In order for the soil C pools to reach steady state, we spun up the model for around 9000 years with fixed land use representative of 1861 and looping over the first 30 years of the ISMSIP2b climate forcing data (1861–1890). During the first 2000 years of spin-up, we ran the model with an atmospheric CO₂ concentration of 350 µ atm and default soil C residence times (τ_carbon) halved, which allowed it to approach steady state more rapidly. Following this, we ran the model for a further 7000 years, reverting to the default τ_carbon values. At the end of this process, the soil C pools had reached approximately steady state, with <0.02 % change in each pool over the final century of the spin-up.

2.4.2 Transient simulations

After the spin-up, we ran a historical simulation from 1861 until the present day, which is 2005 in the case of the ISIMIP2b climate forcing data. We then ran a future simulation until 2099 using the final year of the historical simulation as a restart file. In both of these simulations, climate, atmospheric CO₂, and land cover change were prescribed as fully transient forcings according to the RCP6.0 scenario. For climate variables, we used the IPSL-CM5A-LR model outputs for RCP6.0 bias-corrected by the ISIMIP2b procedure (Frieler et al., 2017; Lange, 2017), while land use change was taken from the 5th Coupled Model Intercomparison Project (CMIP5). As our aim is to investigate long-term trends, we calculated 30-year running means of simulated C flux outputs in order to smooth interannual variations. RCP6.0 is an emissions pathway that leads to a stabilization of radiative forcing at 6.0 watts per square metre (W m⁻²) in the year 2100 without exceeding that value in prior years (Masui et al., 2011). It is characterized by intermediate energy intensity, substantial population growth, middle to high C emissions, increasing cropland area to 2100, and decreasing natural grassland area (van Vuuren et al., 2011). In the paper which describes the development of the future land use change scenarios under RCP6.0 (Hurtt et al., 2011), it is shown that land use change is highly sensitive to land use model assumptions, such as whether or not shifting cultivation is included. The LUH1 reconstruction, for instance, indicates shifting cultivation affecting all of the tropics with a residence time of agriculture of 15 years, whereas the review from Heinimann et al. (2017) revised the area of this type of agriculture downwards, with generally low values in the Congo except in the northeast and southeast, but suggested a shorter turnover of agriculture of only 2 years. In view of such uncertainties, we did not include shifting agriculture in the model. Moreover, there is considerable uncertainty associated with the effect of future land use change in Africa (Hurtt et al., 2011). We chose RCP6.0 as it represents a no mitigation (middle to high emissions) scenario. Moreover, the ISIMIP2b data only provided two Representative Concentration Pathways (RCPs) at the time we performed the simulations: RCP2.6 (low emission) and RCP6.0.

With the purpose of separately evaluating the effects of land use change, climate change, and rising atmospheric CO₂, we ran a series of factorial simulations. In each simulation, one of these factors was fixed at its 1861 level (the first year of the simulation), or in the case of fixed climate change, we looped over the years 1861–1890. The outputs of these simulations (also 30-year running means) were then subtracted from the outputs of the main simulation (original run with all factors varied) so that we could determine the contribution of each driver (Fig. 11, Table 1).

2.5 Evaluation and analysis of simulated fluvial C fluxes

We first evaluated DOC concentrations and pCO₂ at several locations along the Congo main stem (Fig. 1) and on the Oubangui River against the data of Borges at al. (2015b) and Bouillon et al. (2012, 2014). We also compared the various simulated components of the net C balance (e.g. NPP) of the Congo against values described in the literature (Williams et al., 2007; Lewis et al., 2009; Verhegghen et al., 2012; Valentini et al., 2014; Yin et al., 2017). In addition, we assessed the relationship between the interannual variation in present-day (1981–2010) C fluxes of the Congo Basin and variation in temperature and rainfall. This was done through linear regression using STATISTICA^TM. We found trends in several of the fluxes over the 30-year period (1981–2010) and thus detrended the time series with the “Detrend” function, part of the “SpecsVerification” package in R (R Core Team, 2013), before undertaking the statistical analysis focused on the climate drivers of interannual variability.

2.6 Calculating the net carbon balance of the Congo Basin

We calculated net ecosystem production (NEP) by summing the terrestrial and aquatic C fluxes of the Congo Basin (Eq. 1), while we incorporated disturbance fluxes (land use change flux and harvest flux) to calculate net biome production (NBP) (Eq. 2). Positive values of NBP and NEP equate to a net terrestrial C sink.

NEP is defined as follows:

$$\begin{array}{}\text{(1)}& \mathrm{NEP}=\mathrm{NPP}+\mathrm{TF}-\mathrm{SHR}-{\mathrm{FCO}}_{\mathrm{2}}-L{E}_{\mathrm{Aquatic}},\end{array}$$

where NPP is terrestrial net primary production, TF is the throughfall flux of DOC from the canopy to the ground, SHR is soil heterotrophic respiration (only that evading from the terra firme soil surface), FCO₂ is CO₂ evasion from the water surface, and LE_Aquatic is the lateral export flux of C (DOC + dissolved CO₂) to the coast. NBP is equal to NEP except with the inclusion of the C lost (or possibly gained) via land use change (LUC) and crop harvest (HAR). Wood harvest is not included for logging and forestry practices, but during deforestation LUC, a fraction of the forest biomass is harvested and channelled to wood product pools with different decay constants. LUC includes land conversion fluxes and the lateral export of wood product biomass, assuming that wood products from deforestation are not consumed and released as CO₂ over the Congo but in other regions:

$$\begin{array}{}\text{(2)}& \mathrm{NBP}=\mathrm{NEP}-(\mathrm{LUC}+\mathrm{HAR}).\end{array}$$

3 Results

3.1 Simulation of hydrology and aquatic carbon fluxes

The final model configuration is able to closely reproduce the mean monthly discharge at Brazzaville (Fig. 4a, Table 2) and captures the seasonality moderately well (Fig. 4a, Table 2; RMSE = 23 %, R²=0.84 versus RMSE = 29 %, R²=0.23 without calibration; Table A1). At Bangui on the Oubangui River (Fig. 1), the model is able to closely recreate observed seasonality (Fig. 4b; RMSE = 59 %, R²=0.88) but substantially underestimates the mean monthly discharge, with our value being only 50 % of the observed. We produce reasonable NSE values of 0.66 and 0.31 for Brazzaville and Bangui, respectively, indicating that the model is moderately accurate in its simulation of seasonality.

Figure 4 Seasonality of simulated versus observed discharge at (a) Brazzaville on the Congo (for the 1990–2005 monthly mean; Cochonneau et al., 2006), (b) Bangui on the Oubangui (2010–2012; Bouillon et al., 2014), and (c) flooded area in the central (approx. 200 000 km²) part of the Cuvette Centrale wetlands versus GIEMS (1993–2007; Becker et al., 2018). The observed flooded area data represent the maximum minus minimum (permanent water bodies such as rivers) GIEMS inundation. See Fig. 1 for locations.

We also evaluated the simulated seasonal change in flooded area in the central (approx. 200 000 km²; Fig. 1) part of the Cuvette Centrale wetlands against the GIEMS inundation dataset (1993–2007, maximum inundation minus minimum or permanent water bodies; Prigent et al., 2007; Becker et al., 2018). While our model is able to represent the seasonality in flooded area relatively well (R²=0.75; Fig. 4c), it considerably overestimates the magnitude of flooded area relative to GIEMS (Fig. 4c, Table 2). However, the dataset that we used to define the MFF and MFS forcing files (Gumbricht et al., 2017) is produced at a higher resolution than GIEMS and will capture smaller wetlands than the GIEMS dataset; thus, the greater flooded area is to be expected. GIEMS is also known to underestimate inundation under vegetated areas (Prigent et al., 2007; Papa et al., 2010) and has difficulties capturing small inundated areas (Prigent et al., 2007; Lauerwald et al., 2017). Indeed, with the GIEMS data we produce an overall flooded area for the Congo Basin of just 3 %, less than one-third of that produced with the Gumbricht dataset (Gumbricht et al., 2017) or the GLWD (Lehner and Döll, 2004). As such, it is to be expected that there is a large RMSE (272 %; Table 2) between simulated flooded area and GIEMS; more importantly, the seasonality of the two is highly correlated (R²=0.67; Table 2).

In Fig. 5, we compare simulated DOC concentrations at six locations (Fig. 1) along the Congo River and Oubangui tributary against the observations of Borges at al. (2015b). We show that we can recreate the spatial variation in the DOC concentration within the Congo Basin relatively closely with an R² of 0.74 and an RMSE of 24 % (Fig. 5). We are also able to simulate the broad spatial pattern of pCO₂ measured in the main-stem Congo reported by Borges et al. (2019). During the high-flow season (2009–2019 mean of 6 consecutive months of highest flow to account for interannual variation) we simulate a mean pCO₂ of 3373 and 5095 ppm at Kisangani and Kinshasa (Brazzaville), respectively, compared to the observed values of 2424 and 5343 ppm during high water (measured in December 2013; Borges et al., 2019) (Table 3). Similarly, during the low-flow season (mean of 6 consecutive months of lowest flow, 2009–2019) we simulate a mean pCO₂ of 1563 and 2782 ppm at Kisangani and Kinshasa, respectively, compared to the observed values of 1670 and 2896 ppm during falling water (June 2014; Borges et al., 2019) (Table 3).

Figure 5 Observed (Borges et al., 2015a) versus simulated DOC concentrations at several sites along the Congo and Oubangui rivers. See Fig. 1 for locations. The simulated and observed DOC concentrations represent the median values across the particular sampling period at each location detailed in Borges et al. (2015a).

Table 3 Observed (Borges et al., 2019) and modelled pCO₂ (in ppm) at Kinshasa (Brazzaville) and Kisangani on the Congo River at various water levels.

While we are able to recreate observed spatial differences in DOC and pCO₂, as well as broad seasonal variations, we are not able to correctly predict the exact timing of the simulated highs and lows, a reflection of not fully capturing the hydrological seasonality. For example, our mean June pCO₂ at Kinshasa (Brazzaville) is 4470 ppm, while Borges et al. measured a mean of 2896 ppm (Table 3). However, our value for July of 2621 ppm is much closer, and moreover our mean value for December of 5154 ppm is relatively close to the observed value of 5343 ppm. Similarly, we fail to predict the timing of the June falling water at Kisangani (Table 3).

In Fig. 6, we compare simulated pCO₂ against the observed monthly time series at Bangui on the Oubangui River (Bouillon et al., 2012, 2014), which as far as we are aware is the longest time series of pCO₂ published (and accessible) from the Congo Basin, spanning March 2010 to March 2012 (with only the single month of June 2010 missing). Again, while the model fails to correctly predict the precise timing of the peak as with the Kinshasa and Kisangani datasets, the broad seasonal variation in pCO₂ is captured, with the observed and modelled times series ranging from 227–4040 and 415–2928 ppm, respectively (Fig. 6).

Figure 6 Time series of observed versus simulated pCO₂ at Bangui on the River Oubangui. Observed data are from Bouillon et al. (2012, 2014).

3.2 Carbon fluxes along the Congo Basin for the present day

For the present day (1981–2010) we estimate a mean annual terrestrial net primary production (NPP) of 5800±166 Tg C yr⁻¹ (standard deviation, SD) (Fig. 7), corresponding to a mean areal C fixation rate of approximately 1500 g C m⁻² yr⁻¹ (Fig. 8a). We find a significant positive correlation between the interannual variation of NPP and rainfall (detrended R²=0.41, p<0.001; Table A2) and a negative correlation between annual NPP and temperature (detrended R²=0.32, p<0.01; Table A2). We also see considerable spatial variation in NPP across the Congo Basin (Fig. 8a).

Figure 7 Annual C budget (NBP, net biome production) for the Congo Basin for the present day (1981–2010) simulated with ORCHILEAK; NPP is terrestrial net primary productivity, TF is throughfall, SHR is soil heterotrophic respiration, FCO₂ is aquatic CO₂ evasion, C leaching is C leakage to the land–ocean aquatic continuum (FCO₂+LE_Aquatic), LUC is the flux from land use change, and LE_Aquatic is the export C flux to the coast. Range represents the standard deviation (SD) from 1981–2010.

Figure 8 Present-day (1981–2010) spatial distribution of (a) terrestrial net primary productivity (NPP), (b) dissolved organic carbon export from soils and floodplain vegetation into the aquatic system (DOC_inp), (c) CO₂ leaching from soils and floodplain vegetation into the aquatic system (CO_2inp), and (d) aquatic CO₂ evasion (FCO₂). Main rivers are in blue. All are at a resolution of 1^∘.

We simulate a mean soil heterotrophic respiration (SHR) of 5300±99 Tg C yr⁻¹ across the Congo Basin (Fig. 7). Contrary to NPP, interannual variation in annual SHR is positively correlated with temperature (detrended R²=0.57, p<0.0001; Table A2) and inversely correlated with rainfall (detrended R²=0.10), though the latter relationship is not significant (p>0.05). We estimate a mean annual aquatic CO₂ evasion rate of 1363±83 g C m⁻² yr⁻¹, amounting to a total of 235±54 Tg C yr⁻¹ across the total water surfaces of the Congo Basin (Fig. 7), and attribute 85 % of this flux to flooded areas, meaning that only 32 Tg C yr⁻¹ is evaded directly from the river surface. Interannual variation in aquatic CO₂ evasion (1981–2010) shows a strong positive correlation with rainfall (detrended R²=0.75, p<0.0001; Table A2) and a weak negative correlation with temperature (detrended R²=0.09, not significant, p>0.05). Aquatic CO₂ evasion also exhibits substantial spatial variation (Fig. 8d), displaying a pattern similar to both terrestrial DOC leaching (DOC_inp) (R²=0.81, p<0.0001; Fig. 8b) and terrestrial CO₂ leaching (CO_2inp) (R²=0.96, p<0.0001, Fig. 8c) into the aquatic system, but not terrestrial NPP (R²=0.01, p<0.05; Fig. 8a). We simulate a mean annual flux of DOC throughfall from the canopy of 27±1 Tg C yr⁻¹ and C (DOC + dissolved CO₂) export flux to the coast of 15±4 Tg C yr⁻¹ (Fig. 7).

Figure 9 Simulation results for various C fluxes and stocks from 1861–2099 using IPSL-CM5A-LR model outputs for RCP6.0 (Frieler et al., 2017). All panels except for atmospheric CO₂, biomass, and soil C correspond to 30-year running means of simulation outputs. This was done in order to suppress interannual variation, as we are interested in longer-term trends.

For the present day (1981–2010) we estimate a mean annual net ecosystem production (NEP) of 277±137 Tg C yr⁻¹ and a net biome production (NBP) of 107±133 Tg C yr⁻¹ (Fig. 7). Interannually, both NEP and NBP exhibit a strong inverse correlation with temperature (detrended NEP R²=0.55, p<0.0001, detrended NBP R²=0.54, p<0.0001) and weak positive relationship with rainfall (detrended NEP R²=0.16, p<0.05, detrended NBP R²=0.14, p<0.05). Furthermore, we simulate a present-day (1981–2010) living biomass of 41±1 Pg C and a total soil C stock of 109±1 Pg C.

3.3 Long-term temporal trends in carbon fluxes

We find an increasing trend in aquatic CO₂ evasion (Fig. 9a) throughout the simulation period, rising slowly at first until the 1960s when the rate of increase accelerates. In total CO₂ evasion rose by 79 % from 186 Tg C yr⁻¹ at the start of the simulation (1861–1890 mean) (Fig. 10) to 333 Tg C yr⁻¹ at the end of this century (2070–2099 mean; Fig. 10), while the increase until the present day (1981–2010 mean) is +26 % (to 235 Tg C yr⁻¹), though these trends are not uniform across the basin (Fig. A1). The lateral export flux of C to the coast (LE_Aquatic) follows a similar relative change (Fig. 9b), rising by 67 % in total from 12 Tg C yr⁻¹ (Fig. 10) to 15 Tg C yr⁻¹ for the present day and finally to 20 Tg C yr⁻¹ (2070–2099 mean; Fig. 10). This is greater than the equivalent increase in DOC concentration (24 %; Fig. 9b) due to the concurrent rise in rainfall (by 14 %; Fig. 9h) and in turn discharge (by 29 %; Fig. 9h).

Figure 10 Annual C budget (NBP, net biome production) for the Congo Basin for (a) 1861–1890 and (b) 2070–2099, simulated with ORCHILEAK; NPP is terrestrial net primary productivity, TF is throughfall, SHR is soil heterotrophic respiration, FCO₂ is aquatic CO₂ evasion, C leaching is C leakage to the land–ocean aquatic continuum (FCO₂+LE_Aquatic), LUC is the flux from land use change, and LE_Aquatic is the export C flux to the coast. Range represents the standard deviation (SD).

Terrestrial NPP and SHR also exhibit substantial increases of 35 % and 26 %, respectively, across the simulation period and similarly rise rapidly after 1960 (Fig. 9c). NEP, NBP (Fig. 9d), and living biomass (Fig. 9e) follow roughly the same trend as NPP, but NEP and NBP begin to slow down or even level off around 2030, and in the case of NBP, we actually simulate a decreasing trend over approximately the final 50 years. Interestingly, the proportion of NPP lost to the LOAC also increases from approximately 3 % to 5 % (Fig. 9c). We also find that living biomass stock increases by a total of 53 % from 1861 to 2099. Total soil C also increases over the simulation but only by 3 % from 107 to 110 Pg C yr⁻¹ (Fig. 9e). Emissions from land use change (LUC) show considerable decadal fluctuation, increasing rapidly in the second half of the 20th century and decreasing in the mid-21st century before rising again towards the end of the simulation (Fig. 9f). The harvest flux (Fig. 9f) rises throughout the simulation with the exception of a period in the mid-21st century during which it stalls for several decades. This is reflected in the change in land use areas from 1861–2099 (Fig. A2) during which the natural forest and grassland PFTs marginally decrease, while both C₃ and C₄ agricultural grassland PFTs increase.

3.4 Drivers of simulated trends in carbon fluxes

The dramatic increase in the concentration of atmospheric CO₂ (Fig. 9g) and subsequent fertilization effect on terrestrial NPP have the greatest overall impact on all of the fluxes across the simulation period (Fig. 11). They are responsible for the vast majority of the growth in NPP, SHR, aquatic CO₂ evasion, and the flux of C to the coast (Fig. 11a–d). The effect of LUC on these four fluxes is more or less neutral, while the impact of climate change is more varied. The aquatic fluxes (Fig. 11c and d) respond positively to an acceleration in the increase of both rainfall (and in turn discharge; Fig. 9h) and temperature (Fig. 9g) starting around 1970. From around 2020, the impact of climate change on the lateral flux of C to the coast (Fig. 11d) reverts to being effectively neutral, likely a response to a slowdown in the rise of rainfall and indeed a decrease in discharge (Fig. 9h), as well as perhaps the effect of temperature crossing a threshold. The response of the overall loss of terrestrial C to the LOAC (i.e. the ratio of LOAC to NPP; Fig. 11e) is relatively similar to the response of the individual aquatic fluxes, but crucially, climate change exerts a much greater impact, contributing substantially to an increase in the loss of terrestrial NPP to the LOAC in the 1960s and again in the second half of the 21st century. These changes closely coincide with the pattern of rainfall and in particular with changes in discharge (Fig. 9h).

Figure 11 Contribution of the anthropogenic drivers atmospheric CO₂ concentration (CO_2 atm), climate change (CC), and land use change (LUC) to changes in the various carbon fluxes along the Congo Basin under IPSL-CM5A-LR model outputs for RCP6.0 (Frieler et al., 2017).

Overall temperature and rainfall increase by 18 % and 14 % from 24 to 28 ^∘C and 1457 to 1654 mm, respectively, but in Fig. A2 one can see that this increase is non-uniform across the basin. Generally speaking, the greatest increase in temperature occurs in the south of the basin, while it is the east that sees the largest rise in rainfall (Fig. A2). Land use changes are similarly non-uniform (Fig. A2).

The response of NBP and NEP (Fig. 11f and g) to anthropogenic drivers is more complex. The simulated decrease in NBP towards the end of the run is influenced by a variety of factors; LUC and climate begin to have a negative effect on NBP (contributing to a decrease in NBP) at a similar time, while the positive impact (contributing to an increase in NBP) of atmospheric CO₂ begins to slow down and eventually level off (Fig. 11g). LUC continues to have a positive effect on NEP (Fig. 11f) due to the fact that the expanding C₄ crops have a higher NPP than forests, while it has an overall negative effect on NBP at the end of the simulation due to the inclusion of emissions from crop harvest.

4 Discussion

4.1 Congo Basin carbon balance

We simulate a mean present-day terrestrial NPP of approximately 1500 g C m⁻² yr⁻¹ (Fig. 6), substantially larger than the MODIS-derived value of around 1000 g C m⁻² yr⁻¹ from Yin et al. (2017) across central Africa, though it is important to note that satellite-derived estimates of NPP can underestimate the impact of CO₂ fertilization, namely its positive effect on photosynthesis (De Kauwe et al., 2016; Smith et al., 2020). Our stock of the present-day living biomass of 41.1 Pg C is relatively close to the total Congo vegetation biomass of 49.3 Pg C estimated by Verhegghen et al. (2012) based on the analysis of MERIS satellite data. Moreover, our simulated Congo Basin soil C stock of 109±1.1 Pg C is consistent with the approximately 120–130 Pg C across Africa between the latitudes 10^∘ S and 10^∘ N in the review of Williams et al. (2007), within which the Congo represents roughly 70 % of the land area. Therefore, their estimate of soil C stocks across the Congo only would likely be marginally smaller than ours. It is also important to note that neither estimate of soil C stocks explicitly takes into account the newly discovered peat store of 30 Pg C (Dargie et al., 2017), and therefore both are likely to represent conservative values. In addition, Williams et al. (2007) estimate the combined fluxes from conversion to agriculture and cultivation to be around 100 Tg C yr⁻¹ in tropical Africa (largely synonymous with the Congo Basin), which is relatively close to our present-day estimate of harvesting + land use change flux of 170 Tg C yr⁻¹.

Our results suggest that CO₂ evasion from the water surfaces of the Congo is sustained by the transfer of dissolved CO₂ and DOC with 226 and 73 Tg C, respectively, from wetland soils and vegetation to the aquatic system each year (1980–2010; Fig. 8). Moreover, we find that a disproportionate amount of this transfer occurs within the Cuvette Centrale wetland (Figs. 1 and 8) in the centre of the basin, in agreement with a recent study by Borges et al. (2019). In our study, this is due to the large areal proportion of inundated land, facilitating the exchange between soils and aquatic systems. Borges et al. (2019) conducted measurements of DOC and pCO₂, amongst other chemical variables, along the Congo main stem and its tributaries from Kinshasa in the west of the basin (beside Brazzaville; Fig. 1) through the Cuvette Centrale to Kisangani in the east (close to station d in Fig. 1). They found that both DOC and pCO₂ approximately doubled from Kisangani downstream to Kinshasa (Table 3) and demonstrated that this variation is overwhelmingly driven by fluvial–wetland connectivity, highlighting the importance of the vast Cuvette Centrale wetland in the aquatic C budget of the Congo Basin.

Our estimate of the integrated present-day aquatic CO₂ evasion from the river surface of the Congo Basin (32 Tg C yr⁻¹) is the same as that estimated by Raymond et al. (2013) (also 32 Tg C yr⁻¹), downscaled over the same basin area, but smaller than the 59.7 Tg C yr⁻¹ calculated by Lauerwald et al. (2015) and far smaller than that of Borges et al. (2015a) at 133–177 Tg C yr⁻¹ or Borges et al. (2019) at 251±46 Tg C yr⁻¹. The recent study by Borges et al. (2019) is based on the most extensive dataset of Congo Basin pCO₂ measurements to date and thus suggests that we substantially underestimate total riverine CO₂ evasion. As previously discussed, we simulate the broad spatial and temporal variation in observed DOC and pCO₂ (Borges et al., 2015a, b; Fig. 5, Table 3) relatively well. It is therefore somewhat surprising that our basin-wide estimate of riverine CO₂ evasion is so different. Below we discuss some possible explanations for this discrepancy related to methodological differences and limitations.

One potential cause for the differences could be the river gas exchange velocity k. We applied a mean riverine gas exchange velocity k₆₀₀ of 3.5 m d⁻¹, which is similar to the 2.9 m d⁻¹ used by Borges et al. (2015a) but substantially smaller than the mean of approximately 8 m d⁻¹ estimated across Strahler orders 1–10 in Borges et al. (2019) (taking the contributing water surface area of each Strahler order into account). A sensitivity analysis was performed in Lauerwald et al. (2017), which showed that in the physical approach of ORCHILEAK, CO₂ evasion is not very sensitive to the k value, unlike data-driven models. Namely, Lauerwald et al. (2017) showed that an increase or decrease in k₆₀₀ for rivers and swamps (flooded forests) of 50 % only led to a 1 % and −4 % change in total CO₂ evasion, respectively. In ORCHILEAK, k does have an important impact on pCO₂; i.e. a lower k value will increase pCO₂, but this will also lead to a steeper water–air CO₂ gradient and so ultimately to approximately the same FCO₂ over time. In other words, over the scales covered in this research (the large catchment area and water residence times of the Congo), FCO₂ is mainly controlled by the allochthonous inputs of carbon to the river network because by far the largest fraction of these C inputs is leaving the system via CO₂ emission to the atmosphere (as opposed to being laterally transferred downstream). Therefore, we do not consider k to be a major source of the discrepancy. Additionally, our k₆₀₀ value of 0.65 m d⁻¹ for forested floodplains (based on Richey et al., 2002) compares well to recent a study which directly measured k₆₀₀ on two different flooded forest sites in the Amazon Basin, observing a range of 0.24 to 1.2 m d⁻¹ (MacIntyre et al., 2019).

Another potential reason for our smaller riverine CO₂ evasion could be river surface area. We simulate a mean present-day (1980–2010) total river surface area of 25 900 km² compared to the value of 23 670 km² used in Borges et al. (2019, Supplement), so similarly we think that this can be discounted as a major source of discrepancy. However, it should be noted that both estimates are high compared to the recent estimate of 17 903 km² based on an analysis of Landsat images (Allen and Pavelsky, 2018).

The difference in our simulated riverine CO₂ evasion compared to the empirically derived estimate of Borges et al. (2019) could be caused by the lack of representation of aquatic plants in the ORCHILEAK model. Borges et al. (2019) used the stable isotope composition of δ¹³C–DIC to determine the origin of dissolved CO₂ in the Congo River system and found that the values were consistent with a DIC input from the degradation of organic matter, in particular from C₄ plants. Crucially, they further found that the δ¹³C–DIC values were unrelated to the contribution of terra firme C₄ plants; rather, they were more consistent with the degradation of aquatic C₄ plants, namely macrophytes. ORCHILEAK does not represent aquatic plants, and the wider LSM ORCHIDEE does not have an aquatic macrophyte PFT (though root respiration of floodplain plants for the PFTs represented is accounted for as a C source). This could at the very least partly explain our conservative estimate of river CO₂ evasion given that tropical macrophytes have relatively high NPP. Rates as high as 3500 g C m⁻² yr⁻¹ have been measured on floodplains in the Amazon (Silva et al., 2009). While this value is higher than the values simulated in the Cuvette Centrale by ORCHILEAK (Fig. 8), they are of the same order of magnitude, so this alone cannot fully explain the discrepancy compared to the results of Borges et al. (2019). In the Amazon Basin it has been shown that wetlands export approximately half of their gross primary production (GPP) to the river network compared to upland (terra firme) ecosystems, which only export a few percent (Abril et al., 2014). More importantly, Abril et al. (2014) found that tropical aquatic macrophytes export 80 % of their GPP compared to just 36 % for flooded forest. Therefore, the lack of a bespoke macrophyte PFT is indeed likely to be one reason for the discrepancy between our results and those of Borges but largely due to their particularly high export efficiency to the river–floodplain network as opposed to differences in NPP. While being a significant limitation, creating and incorporating a macrophyte PFT would be a substantial undertaking given that the authors are unaware of any published dataset which has systematically mapped their distribution and abundance. It is important to note that while ORCHILEAK does not include the export of C from aquatic macrophytes, it also neglects their NPP. Moreover, most aquatic macrophytes described in the literature have short (<1 year) life cycles (Mitchel and Rogers, 1985). As such, while this model limitation is likely one of the causes for our relatively low estimate of riverine CO₂ evasion, it will only have a limited net effect on our estimate of the overall annual C balance (NBP, NEP) of the Congo Basin.

Finally, another cause for the difference in riverine CO₂ evasion could be that the resolution of ORCHILEAK (0.5^∘ river network and 1^∘ for C fluxes) is not sufficient to fully capture the dynamics of the smallest streams of the Congo Basin, which have been shown to have the highest DOC and CO₂ concentrations (Borges et al., 2019). Indeed, ORCHILEAK typically does not simulate the highest observed pCO₂ measurements of the smallest tributaries (i.e. >16 000 ppm). This is partly because for the fast reservoir (headwaters) in ORCHILEAK we assume full pCO₂ equilibrium with the atmosphere over one full day, which prevents very high pCO₂ values from building in the water column.

Despite these limitations, it is important to note that in our simulations, the evasion flux from rivers only contributes 15 % of total aquatic CO₂ evasion, and including the flux from wetlands and floodplains, we produce a total of 235 Tg C yr⁻¹. Moreover, the majority of this evasion occurs in the Cuvette Centrale (Fig. 8), which suggests that while ORCHILEAK fails to attribute a large portion of this flux to small rivers (owing to the coarse resolution of the river network) we nonetheless do capture the source of carbon. In other words, in ORCHILEAK the majority of this carbon evades directly from the floodplain and wetlands of the Cuvette Centrale as opposed to the small rivers.

Our simulated export of C to the coast of 15 (15.3) Tg C yr⁻¹ is virtually identical to the TOC + DIC export estimated by Borges et al. (2015a) of 15.5 Tg C yr⁻¹, which is consistent with the fact that we simulate a similar spatial variation of DOC concentrations (see Figs. 8 and 1 for locations). It is also relatively similar to the 19 Tg C yr⁻¹ (DOC + DIC) estimated by Valentini et al. (2014) in their synthesis of the African carbon budget. Valentini et al. (2014) used the largely empirically based global nutrient export from the WaterSheds (NEWS) model framework, and they point out that Africa was underrepresented in the training data used to develop the regression relationships which underpin the model, and this could thus explain the small disagreement.

Of the total 15 Tg C yr⁻¹ exported to the coast, we simulate a 2.4 Tg C yr⁻¹ component of dissolved CO₂, which is relatively similar to the empirically derived estimate of the total DIC export of 3.3 Tg C yr⁻¹ calculated in Wang et al. (2013). According to Wang et al., dissolved CO₂ accounts for the majority (1.9 Tg C yr⁻¹), with the rest being the weathering derived flux of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$. Thus, the discrepancy between the two estimates is likely to be largely caused by our lack of accounting for the weathering-derived flux (${\mathrm{HCO}}_{\mathrm{3}}^{-}$), which they estimate at 1.4 Tg C yr⁻¹. In summary, despite this model limitation the results of Wang et al. (2013) suggest that we still capture the majority of the DIC flux.

4.2 Trends in terrestrial and aquatic carbon fluxes

There are relatively sparse observed data available on the long-term trends of terrestrial C fluxes in the Congo. Yin et al. (2017) used MODIS data to estimate NPP between 2001 and 2013 across central Africa. They found that NPP increased on average by 10 g C m⁻² per year, while we simulate an average annual increase of 4 g C m⁻² yr⁻¹ over the same period across the Congo Basin. The two values are not directly comparable as they do not cover precisely the same geographic area, but it is encouraging that our simulations exhibit a similar trend as remote sensing data. As previously noted, MODIS-derived estimates of NPP do not fully include the effect of CO₂ fertilization (de Kauwe et al., 2016), whereas ORCHILEAK does. Thus, the MODIS NPP product may underestimate the increasing trend in NPP, which would bring our modelled trend further away from this dataset. On the other hand, forest degradation effects and recent droughts have been associated with a decrease in greenness (Zhou et al., 2014) and above-ground biomass loss (Qie et al., 2019) in tropical forests.

Up to a point, our results also concur with estimates based on the upscaling of biomass observations (Lewis et al., 2009; Hubau et al., 2020). Lewis et al. (2009) upscaled forest plot measurements to determine that intact tropical African forests represented a net uptake of approximately 300 Tg C yr⁻¹ between 1968 and 2007, and this is consistent with our NEP estimate of 275 Tg C yr⁻¹ over the same period. However, more recently an analysis based on an extension of the same dataset found that the above-ground C sink in tropical Africa was relatively stable from 1985 to 2015 (Hubau et al., 2020).

A major source of the uncertainty associated with future projections of NPP and NEP comes from our limited understanding and representation of the CO₂ fertilization effect. Recent analysis of data from some of the longest-running Free-Air CO₂ Enrichment (FACE) sites, consisting of early-successional temperate ecosystems, found a 29.1±11.7 % stimulation of biomass over a decade (Walker et al., 2019). A meta-analysis (Liu et al., 2019) of seven temperate FACE experiments combined with process-based modelling also found substantial sensitivity (0.64±0.28 Pg C yr⁻¹ per 100 ppm) of biomass accumulation to atmospheric CO₂ increase, and the same study showed that ORCHIDEE model simulations were largely consistent with the experiments. However, other FACE experiments on mature temperate forests (Körner et al., 2005) and eucalyptus forests bring into question whether the fertilization effects observed in temperate FACE experiments can be extrapolated to other ecosystems. For example, the Swiss FACE study, a deciduous mature forest, found no significant biomass increase with enhanced CO₂ (Körner et al., 2005), while a FACE experiment on a mature eucalyptus forest in Australia found that while CO₂ stimulated an increase in C uptake through GPP, this did not carry to the ecosystem level, largely as a result of a concurrent increase in soil respiration (Jiang et al., 2020). Unfortunately, no results are yet available from any tropical FACE experiments, though the Amazon FACE experiment is underway and the eventual results will be crucial in developing our understanding of the CO₂ fertilization effect beyond the temperate zone.

With these limitations in our understanding of tropical forest ecosystems in mind, over the entire simulation period (1861–2099) we estimate that aquatic CO₂ evasion will increase by 79 % and the export of C to the coast by 67 %. While there are no long-term observations of aquatic CO₂ evasion in the Congo, a recent paper examined trends in observed DOC fluxes in the Congo at Brazzaville–Kinshasa over the last 30 years (N'kaya et al., 2020). The authors found a 45 % increase in the annual flux of DOC from 11.1 Tg C yr⁻¹ (mean from 1987–1993) to 16.1 Tg C yr⁻¹ (mean from 2006–2017). Comparing the same two periods, we find a smaller increase of 15 % from 12.3 to 14.2 Tg C yr⁻¹. While our increase is substantially smaller, these observations are still over relatively short timescales, and thus interannual variations could have a considerable influence over the means of the two periods. Irrespectively, it is encouraging that observations concur with the overall simulated increasing trend. Perhaps most interesting is that N'kaya et al. (2020) attribute this increase to hydrological changes, specifically an increase in flood events in the central basin (including the Cuvette Centrale). Over this period, we too attribute the increase in carbon fluxes to the coast in part to climate change (Fig. 11d), and over the full simulation period, the largest increase in DOC and CO₂ leaching into the aquatic system occurs within the Cuvette Centrale (Fig. A1).

Comparing our results to models of other basins, our simulated increases in outgassing (79 %) and the export of C to the coast (67 %) are considerably greater than the 23 % and 27 % rises predicted for the Amazon Basin (Lauerwald et al., 2020) over the same period and under the same scenario. This is largely due to the fact climate change is predicted to have a substantial negative impact on the aquatic C fluxes in the Amazon, something that we do not find for the Congo where rainfall is projected to increase over the 21st century (RCP6.0). In the Amazon, Lauerwald et al. (2020) show that while there are decadal fluctuations in precipitation and discharge, total values across the basin remain unchanged in 2099 compared to 1861. However, changes in the spatial distribution of precipitation mean that the total water surface area actually decreases in the Amazon. Indeed, while we find an increase in the ratio of C exports to LOAC ∕ NPP from 3 % to 5 %, Lauerwald et al. (2020) find a comparative decrease.

Our simulated increase in DOC export to the coast up to the present day is smaller than findings recently published for the Mississippi River using the Dynamic Land Ecosystem Model (DLEM; Ren at al., 2016). In addition, the Mississippi study identified LUC including land management practices (e.g. irrigation and fertilization), followed by change in atmospheric CO₂, as the biggest factors in the 40 % increase in DOC export to the Gulf of Mexico (Ren et al., 2016). Another recent study (Tian et al., 2015) found an increase in DIC export from eastern North America to the Atlantic Ocean from 1901–2008 but no significant trend in DOC. The authors demonstrated that climate change and increasing atmospheric CO₂ had a significant positive effect on long-term C export, while land use change had a substantial negative impact.

4.3 Limitations and further model developments

It is important to note that we can have greater confidence in the historic trend (until the present day), as future changes are more reliant on the skill of Earth system model predictions and of course on the accuracy of the RCP6.0 scenario. As discussed above, our understanding and representation of CO₂ fertilization, especially in the tropics, is a major limitation. Moreover, the majority of land surface models, ORCHILEAK included in its current iteration, do not represent the effect of nutrient limitation on plant growth, meaning that estimates of land C uptake may be too large (Goll et al., 2017). There are also considerable uncertainties associated with future climate projections in the Congo Basin (Haensler et al., 2013). Nutrient limitation on growth and a better representation of the effects of enhanced CO₂, particularly with regards to soil respiration (Jiang et al., 2020) and tree mortality (Hubau et al., 2020), are two crucial aspects which need to be further developed.

Additionally, we do not account for methane fluxes from Congo wetlands, estimated at 1.6 to 3.2 Tg (CH₄) per year (Tathy et al., 1992), and instead assume that all C is evaded in the form of CO₂. Another limitation is the lack of accounting for bespoke peatland dynamics in the ORCHILEAK model. ORCHILEAK is able to represent the general reduction in C decomposition in water-logged soils, and indeed Hastie et al. (2019) demonstrated that increasing the maximum floodplain extent in the Amazon Basin led to an increase in NEP despite fuelling aquatic CO₂ evasion because of the effect of reducing soil heterotrophic respiration. Furthermore, ORCHILEAK uses a “poor soils” mask forcing file (Fig. 2j) based on the Harmonized World Soil Database (FAO/IIASA/ISRIC/ISS-CAS/JRC, 2009), which prescribes reduced decomposition rates in low-nutrient and low-pH soils (e.g. Podzols and Arenosols). The effect of the “poor soils” forcing can clearly be seen in the spatial distribution of the soil C stock in Fig. A3, where the highest C storage coincides with the highest proportion of poor soils. Interestingly, this does not include the Cuvette Centrale wetlands (Fig. 1), an area which was recently identified as containing the world's largest intact tropical peatland and a stock of around 30 Pg C (Dargie at al., 2017). One potential improvement that could be made to ORCHILEAK would be the development of a new tailored “poor soils” forcing file for the Congo Basin which explicitly includes Histosols, perhaps informed by the Soil Grids database (Hengl et al., 2014), to better represent the Cuvette Centrale. This could, in turn, be validated and/or calibrated against the observations of Dargie et al. (2017). A more long-term aim could be the integration or coupling of the ORCHIDEE-PEAT module with ORCHILEAK. ORCHIDEE-PEAT (Qiu et al., 2019) represents peat as an independent sub-grid hydrological soil unit in which peatland soils are characterized by peat-specific hydrological properties and multi-layered transport of C and water. Thus far, it has only been applied to northern peatlands, and calibrating it to tropical peatlands, along with integrating it within ORCHILEAK, would require considerable further model development but would certainly be a valuable longer-term aspiration. This could also be applied across the tropical region and would allow us to comprehensively explore the implications of climate change and land use change for tropical peatlands. In addition, ORCHILEAK does not simulate the erosion and subsequent burial of POC within river and floodplain sediments. Although it does not represent the lateral transfer of POC, it does incorporate the decomposition of inundated litter as an important source of DOC and dissolved CO₂ to the aquatic system; i.e. it is assumed that POC from submerged litter decomposes locally in ORCHILEAK. Moreover, previous studies have found that DOC as opposed to POC (Spencer et al., 2016; Bouillon et al., 2012) overwhelmingly dominates the total load of C in the Congo.

The representation of the rapid C loop of aquatic macrophytes should also be made a priority in terms of improving models such as ORCHILEAK, particularly in the tropics. As previously discussed, ORCHILEAK also fails to account for the weathering-derived flux (${\mathrm{HCO}}_{\mathrm{3}}^{-}$). Finally, the issue of shifting cultivation demands further attention; at least for the present day a shifting cultivation forcing file could be developed based on remote sensing data (Tyukavina et al., 2018). For an additional discussion of the limitations of ORCHILEAK, please also see Lauerwald et al. (2017) and Hastie et al. (2019).

5 Conclusions

For the present day, we show that aquatic C fluxes, in particular CO₂ evasion, are important components of the Congo Basin C balance that are larger than, for example, the combined fluxes from LUC and harvesting, with around 4 % of terrestrial NPP being exported to the aquatic system each year. Our simulations show that these fluxes may have undergone considerable perturbation since 1861 to the present day and that under RCP6.0 this perturbation could continue; over the entire simulation period (1861–2099), we estimate that aquatic CO₂ evasion will increase by 79 % and the export of C to the coast by 67 %. We further find that the ratio of C exports to the LOAC ∕ NPP could increase from 3 % to 5 %, driven by both rising atmospheric CO₂ concentrations and climate change. This calls for long-term monitoring of C levels and fluxes in the rivers of the Congo Basin and further investigation of the potential impacts of such change. Our results also highlight the limitations of the current generation of land surface models and call for investment in further model development.

Appendix A

Table A1 Performance statistics for the modelled versus observed seasonality of discharge on the Congo at Brazzaville.

Table A2 Pearson correlation coefficient (r) between detrended carbon fluxes and detrended climate variables.

Figure A1 Change (Δ, 2099 minus 1861) in the spatial distribution of (a) terrestrial NPP, (b) DOC leaching into the aquatic system, (c) CO₂ leaching into the aquatic system, and (d) aquatic CO₂ evasion. All are at a resolution of 1^∘.

Figure A2 Change (Δ, 2099 minus 1861) in the spatial distribution of the principal climate and land use drivers across the Congo Basin; (a) mean annual temperature (^∘C), (b) mean annual rainfall (mm yr⁻¹), and (c–h) mean annual maximum vegetated fraction for PFTs 2, 3, 10, 11, 12, and 13. All are at a resolution of 1^∘.

Figure A3 Spatial distribution of simulated total carbon stored in soils for the present day (1981–2020).