Global carbon sequestration through continental chemical weathering in a climatic change context

This study simulates carbon dioxide (CO2) sequestration in 300 major world river basins (about 70% of global surface area) through carbonates dissolution and silicate hydrolysis. For each river basin, the daily timescale impacts under the RCP 2.6 and RCP 8.5 climate scenarios were assessed relative to a historical baseline (1969–1999) using a cascade of models accounting for the hydrological evolution under climate change scenarios. Here we show that the global temporal evolution of the CO2 uptake presents a general increase in the annual amount of CO2 consumed from 0.247 ± 0.045 Pg C year−1 to 0.261 and 0.273 ± 0.054 Pg C year−1, respectively for RCP 2.6 and RCP 8.5. Despite showing a general increase in the global daily carbon sequestration, both climate scenarios show a decrease between June and August. Such projected changes have been mapped and evaluated against changes in hydrology, identifying hot spots and moments for the annual and seasonal periods.


2−
), which is later exported by rivers to other water bodies [2][3][4][5][6][7][8][9][10][11] . Riverine dissolved loadings may be used as a proxy for chemical weathering assessment 4,11 , and it has been recognised that hydrology impacts these flux-discharge relationships [5][6][7] . Nevertheless, it is unclear how the potential impacts of climate change will affect these dissolved mass balances. Forecasting how potential shortterm (daily or monthly) shifts in hydrology under a changing climate may alter these fluxes of riverine matter is needed to assess the potential evolution of the global carbon cycle under a wide range of future scenarios. Here we hypothesized that a more significant seasonality of the hydrological cycle would cause a heterogeneous spatial and temporal pattern in CO 2 sequestration at the global scale.
Climate change scenarios show increases in air temperatures and significant changes in the hydrological cycle 12 . Still, its effect on chemical weathering and related CO 2 uptake dynamics and dissolved solids exportation by rivers is poorly understood 13 . For instance, on one side, an increase in the soils microbial activity is expected, implying an increase in CO 2 production through respiration 14,15 . However, higher air temperatures may decrease CO 2 dissolution in water, as carbonate weathering follows a boomerang-shaped evolution with increasing temperature 16 . Hydrology is the dominant driver of matter transport at large scale 17 . Nevertheless, the potential implication of respiration and temperature in the global balance remains poorly understood 13 .
How will the annual soil CO 2 consumption evolve under these hydrological shifts? Where and when will these soil CO 2 consumption changes be more relevant? The answer to these scientific questions needs a largescale comprehensive field study. However, such a study is challenging regarding the in situ resources needed to measure the variables involved in these biogeochemical cycles.
In this complex framework, modelling arises as an alternative approach to overcome this challenge by yielding insights into where and when hot spots and hot moments (places and times with disproportionally high chemical weathering rates) may be located. Identifying potential hot spots and moments of change is important to target areas needing the deployment of resources to understand in situ processes, which could yield relevant  19,20 , link thermodynamics, kinetics and transport processes to hydrology and vegetation. They have yielded insights, for instance, that 40% of the total increase in CO 2 consumption in the Mackenzie River basin is related to the direct effect of climate change on hydrology alteration 15 . Moreover, they highlighted that a "chemostatic" behaviour of riverine dissolved load suggests that hydrological changes impact the dissolved element concentration in rivers 19,20 . Empirical models are commonly developed at large spatial scales 7,21 , and their results are used in relevant assessments such as the ESCOBA program 22 and the IPCC Assessment Reports 12 . The objective of this study was to assess the impact of future modifications of the hydrological cycle on the global CO 2 sequestration by weathering. Based on the largest basins in the world, we aimed at identifying the watersheds or areas with higher CO 2 consumption compared to the rest of the globe (hot spots 23 ) and the periods showing larger sequestration (hot moments 23 ) on which to focus research efforts.
In this study, a cascade of models was set up to cover the 300 largest basins in the world, obtaining the spatial and temporal potential evolution of CO 2 sequestration under two climate change scenarios (RCP 2.6 and RCP 8.5). Historical and future dynamics of cations and anions derived from chemical weathering of rocks and associated carbon sequestration were modelled by coupling the hydrological model VIC with the geochemical model ICWR. Then, the results from this set-up were used as an input in the MEGA model to estimate CO 2 consumption. Observed data have been used to evaluate the results at each step of the model cascade. The reader is referred to the methods section and the Extended Data Fig. 2.

Results and discussion
Interannual and seasonal fluctuations. The annual mean inorganic C sequestration through chemical weathering during the Historical period  amounts to 0.247 Pg C year −1 (1 Pg = 10 15 g), in the lower range of previous research 4, [8][9][10][11]21,[24][25][26][27] , spanning from 0.22 to 0.30 Pg C year −1 . The climate change scenarios present a potential increase in the amount of C sequestered, reaching mean values of 0.261 and 0.273 Pg C year −1 for RCP 2.6 and RCP 8.5, respectively. The annual series (Fig. 1a) shows a similar temporal evolution for the middle part of the century, but such a difference enlarges within the "End of Century Projection" (ECP) period (2069-2099). An annual average increase in C sequestration of 6% (RCP 2.6) and 10% (RCP 8.5) is forecasted if the Historical period is considered the baseline. The annual increases represent < 0.1% of the anthropogenic emissions projected for 2100 in the RCP 8.5 scenario and ~ 0.7% of the "negative emissions" projected for the RCP2.6 scenario 28 . These C sequestration increases are low at the global scale; however, these become more relevant when evaluating the change percentages at a lower time-step (daily).
When comparing the mean daily temporal evolution in the ECP and the Historical periods (Fig. 1b), more considerable differences are found for the RCP 8.5 scenario. In both scenarios, a positive difference is found for most of the year, agreeing with a more considerable amount of C sequestered during the ECP period compared to the Historical one. Such increases may reach a positive 20% value in May for the RCP 8.5 scenario, followed by a shift towards a decrease between June and August. It suggests a lower CO 2 consumption through continental chemical weathering during these months. Previous studies highlighted that hydrology is a dominant driver in the temporal evolution of CO 2 consumption 29 . Therefore, a decrease in CO 2 consumption suggests that a lower amount of water is being exported in this period. This decrease coincides with the low discharge period in the www.nature.com/scientificreports/ northern latitudes, suggesting that these regions significantly influence the global amount of soil CO 2 consumed through chemical weathering. The Northern Hemisphere has a more considerable amount of continental land and the hydrological cycle is expected to intensify in these areas. This suggests that the annual shift found between June and August regarding the global carbon weathering CO 2 uptake could be a consequence of hydrological changes at these latitudes.
Hot spots and moments for CO 2 uptake. On a global scale, the spatial distribution of carbon sequestration during the Historical period ( Fig. 2a) is comparable to what is described in previous literature 7 . Similarly, values are analogous to regional cases, such as the Amazon 11,29,30 , Congo 5,21 , Niger 31 , Garonne 21 , the Alps region 32 or the United States continuum 33 , even if the lateritic soils have lower CO 2 uptake, as shown by Boeglin & Probst 31 . Alterations in this spatial pattern (Fig. 2b) are considered primary consequences of hydrological changes since the geochemical processes take place over a much larger time scale than biogeochemical processes on the Earth surface. More significant changes are found in the northern latitudes, where an increase in the annual discharges in these river basins is expected due to a larger precipitation amount, an earlier onset and a more intensive snowmelt in spring 34 . Even though a higher seasonality will change the hydrological cycle within the year.
A general increasing trend in CO 2 sequestration is found during the January-February-March (JFM) period (Extended Data Fig. 1), especially in basins such as the Yenisei and the Ob, which present changes > 25%, particularly in the RCP 8.5 scenario. Increasing temperatures are expected to cause an accelerated snowmelt 34 , which could lead to the increase of dissolved inorganic carbon river fluxes. In contrast, heterogeneous projections are www.nature.com/scientificreports/ found for southern latitudes, as in the Murray-Darling River basin in Australia and the Orange River in South Africa. These basins are located under mixed climatic zones 35 , which could explain the heterogeneous results found in the scenarios. The potential evolution in the river discharge seasonality in these areas may cause different geochemical fluxes since a discharge decrease of more than 25% in the low-flow period is expected 34 . However, further analysis using mechanistic approaches or in situ measurements is needed to understand and assess potential changes in biogeochemical cycles in these regions during the January to March (JFM) period. A different pattern is found for the July-August-September (JAS) season, when most of the northern latitudes present negative changes, suggesting a decrease in the CO 2 uptake due to chemical weathering. This shift in the northern latitudes is attributed to a projected decline in discharge, notably on the RCP 8.5 scenario 36 .
As discharge intensity is a primary critical factor of CO 2 flux changes in the short-term (daily, monthly or seasonal) assessment, we compared the relative changes in CO 2 uptake in a subset of river basins for both scenarios (Fig. 3). Such a comparison illustrates how watersheds located in tropical and cold climates present significant differences for CO 2 uptake, especially basins in the Siberian region (e.g. Kolyma, Khatanga or Lena) or active arc-islands in New Guinea (e.g. Fly or Sepik). These changes show that, even with their smaller draining areas and discharges, the role of arc-islands should be taken into consideration and should be a focus of future research as their relative role in carbon sequestration through chemical weathering is relevant 7 . Besides, the potential impacts of climate change in these tropical areas are expected to be more relevant in these smaller river basins than in larger ones like the Amazon or Congo 7 .
Similarly, the relative change in CO 2 sequestration on watersheds with mixed climates (e.g. the Nile, Krishna, Niger or Murray Darling) is large. It could be related to higher amounts of carbonate rock outcrops in these areas, releasing significant amounts of CO 2 during weathering 37 . Regarding cold and polar basins, permafrost processes were not included in the present analysis. They are expected to significantly impact the carbon sequestration by weathering and the loads to rivers 38 .
Even though there are some relative changes noted at the annual scale, these differences are not maintained throughout the year. The short-term temporal evolution should be a focus of interest for future research, especially in those areas located in tropical and cold climates since they seem to be more sensitive to climate change.
Modelling approach limitations. This study uses a model output as input to another model, which implies a propagation of the error in each step. We assessed the uncertainties of our approach with previously published data. The results show a very significant correlation to observed data 27 for the ICWR model (ρ SPEARMAN = 0.83, p < 0.01, n = 173, see Extended Data Fig. 5), as well as for the MEGA model results in the Historical period 6 (ρ SPEARMAN = 0.71, p < 0.01, n = 48, see Extended Data Fig. 6). In summary, even though a cascade of models usually accumulates errors, the present modelling approach fits well with current observed data.
Continental weathering occurs at the critical zone, where slow geological processes interact with faster biogeochemical ones. The present modelling approach includes a macroscale hydrological model, which resolves the energy and water balance. The main uncertainties are discussed elsewhere 36 . Concerning the ICWR model, the parameters have been fitted to observed data for a given situation. These laws may change under a climate change scenario. The MEGA model balance is based on a set of hypotheses described in the literature and summarised in the Extended Data Fig. 4. The main limitation of this modelling set-up is that it does not simulate transient states, as the chemical weathering rates are described through linear equations. Regolith thickness is www.nature.com/scientificreports/ kept constant during simulation even though it is coupled to physical erosion, which could change chemical weathering rates 39 . Temporal variation of below-ground CO 2 concentration is not included in the ICWR model, even though it has been noted as relevant in silicate weathering and affected by vegetation changes 40 . Further uncertainties are related to the resolution of the lithological map (see "Methods" section). It has already been shown that using a map with a finer resolution helped improve the modelling results at local scale 41 .

Conclusion
This study presents a forecast of the potential evolution of CO 2 consumption during the chemical weathering of continental rocks at the global scale under two climate change scenarios. The research provides two main new insights on the temporal evolution of the carbon cycle in the short-term, derived from shifts in the water cycle: (a) Even though there is a trend towards increasing the amount of carbon sequestered at the annual and global scale, this pattern is not homogeneously distributed neither in space nor time, and (b) the river basins located in cold climates (such as the Kolyma, Lena, Yukon, Kuskokwim or Yenisei) and in tropical climates with a lower draining area (like the Fly, Sepik or Narmada) may be more sensitive to the effects of climate change.
The insights presented suggest that research on the terrestrial part of the carbon cycle should focus on these areas to constrain the main cycle-shaping drivers better, targeting an understanding of the carbon fluxes' temporal dynamics. This study provides the first snapshot of carbon sequestration's potential evolution through the chemical weathering of rocks on a global scale, which may be used for comparison in future studies in this field. Even if the carbon fluxes from chemical weathering are small compared to the total CO2 emissions, they are in the same order of magnitude with the natural CO2 emissions to the atmosphere. Local increases in chemical weathering could lead to acidification of local surface waters, which could affect the equilibrium of freshwater and coastal ecosystems.

Methods
Four Representative Concentration Pathways (RCP) were selected by the Intergovernmental Panel on Climate Change (IPCC) in its fifth Assessment Report 12 regarding climate change impacts on the carbon cycle. Two of those four RCPs are included in the present analysis (RCP 2.6 and RCP 8.5). RCP 2.6 considers that the radiative forcing (i.e. the strength of change drivers) will peak at about 2.6 W m −2 in 2100 and then decline, whereas RCP 8.5 considers that the radiative forcing will reach > 8.5 W m −2 by 2100. The highest and lowest RCPs were selected to capture the broadest range of future expected atmospheric greenhouse gas concentrations. Biascorrected climate model outputs of five general circulation models (GCM) were included in the VIC hydrological model to simulate impacts on daily river discharge on 0.5° × 0.5° at the global scale 36 48 . Validation results showed that the observed hydrological conditions were realistically represented by the VIC hydrological model for most river basins. The average daily discharges from the five GCM and its respective uncertainty were taken as an input for the geochemical model.
The input data handling and modelling in the present study follow the workflow summarised in Extended Data Fig. 2. First, a subset of the 300 largest river basins covering ~ 70% of the global land area was selected from the HydroBASINS dataset 49 , including only exoreic river basins (Extended Data Fig. 3). The catchment area selected was derived from the HydroBASINS dataset. Second, the outlet points of river basins from the simulated river discharge VIC database were selected manually by exploring the time series for each sampling location and the surrounding points (± 0.5°) to choose the correct outlet. When a catchment presented several outlets (especially for endorheic catchments draining to a lake or inner sea), all of them were considered, and their discharge were added to account for all water flowing out of the drained area. Deconvolution of the daily discharge signal for each outlet is applied to distinguish between the surface runoff from the soil and groundwater flow. Third, draining area characteristics are summarised for each river basin by considering the relative abundance of lithological classes contained in the Global Lithological Map 50 (GLIM), the soil classes as defined by the Harmonized World Soil Database 51 (HWSD), and the climatic zones present according to the Köppen-Geiger classification presented by Beck et al. 35 . The modelling approach consists of two models in cascade; the flux of ions derived from the chemical weathering of rocks is computed using the ICWR 17 model at a daily time step for each river basin under each discharge time series. Then, the Major Element Geochemical Approach 37 (MEGA) model uses these loadings to estimate the CO 2 consumed through chemical weathering.
From this set-up, a daily time series of major ion river fluxes and CO 2 consumption through chemical weathering from the 300 largest river basins were obtained from 1965 to 2100. Two periods are established to estimate the potential effect of hydrological shifts on weathering CO 2 consumption: the Historical (October 1969-September 1999) and the ECP (October 2069-September 2099). A comparison between the two periods is also accomplished. Daily discharge time series have been taken from all these scenario-model combinations while temporal aggregation is performed after applying the ICWR and MEGA models.
The stream discharge time series have been deconvoluted to separate the influence of the surface runoff from the baseflow and interflow discharges. Such deconvolution is performed following the digital filter shown in Eq. (1), as proposed by Eckhardt 52 . In the present study, only interflow and baseflow are considered significantly transporting ions, while surface runoff is the principal agent for dilution. A digital filter based on the slopes of the increasing and decreasing parts of the hydrograph was selected. It is sensitive to changes in the amount of water and seasonality and flood events. Nonetheless, this is a source of uncertainty as the infiltration process depends on catchment-to-catchment properties, while this equation assumes linearity between the groundwater outflow (baseflow) and its storage 52 . www.nature.com/scientificreports/ where b(t) is the baseflow time series, BFI mx represents the long-term ratio of baseflow to total streamflow, a is the filter parameter, and Q(t) is the total discharge. When b(t) > Q(t), then b(t) is replaced by Q(t). The parameters that were selected for the separation were BFI mx = 0.8 and a = 0.95 (as recommended by Xie et al. 53 ). These parameters were kept constant in all catchments, which also became a source of uncertainty. The ICWR model is an empirical model developed to estimate major ion riverine fluxes released by chemical weathering of rocks at the global scale, based on Eq. (2). It has been validated at the global scale under static conditions 17 and at the local scale under dynamic evolution 41 . Computing it requires a description of the draining area of each river basin in terms of soil cover and lithological distribution. The parameters of the equation take into account the concentration for each ion x drained from a lithological class i. The ICWR model parameters were calibrated after an atmospheric input correction. Thus, the loadings computed from this equation only relate to the chemical weathering process, which is needed as an input for the MEGA model.
where F x represents the specific flux of ion x in mol km −2 year −1 , b represents the baseflow and interflow obtained from the deconvolution of the total discharge signal, f sx is a factor considering the soil shielding effect on CW (adim), L i is the relative abundance of a lithological class i in terms of the total area of the drainage basin (expressed in the 0-1 range), and C ix is the calibrated parameter that represents the concentration of the ion x draining from the lithological class (expressed in mol L −1 ). This equation is applied to the time series of the deconvoluted discharge, including the baseflow and interflow.
The MEGA model is based on the mass balance shown in Extended Data Fig. 4 37 . The input data are the riverine molar loadings of each major ion, after the atmospheric input correction. We assume that anthropogenic influence is not captured in the measurements of major ion concentrations. R pyr and R sil molar ratios for each catchment account for the Ca 2+ + Mg 2+ loads coming from the weathering of silicate rocks and the SO 4 2− originating from pyrite oxidation, respectively. The results from this model relate to the CO 2 consumed by the CW of rocks.
Considering the river loadings of all major ions, the CO 2 consumed is computed as follows: first, the atmospheric deposition (wet and dry) must be removed from these loads, remaining the molar fluxes derived from CW. All Cl − is associated with the dissolution of halite (NaCl), though the remaining Na + is associated with Nasilicates. If the Cl − > Na + , the remaining Cl − is linked to the dissolution of sylvite (KCl). Evaporites do not uptake CO 2 when dissolving, while Na + and K − silicate rocks (e.g. albite NaAl 3 O 8 and orthoses KAlSi 3 O 8 ) require 1 mol CO 2 for each ion mol released to water. Then, using R pyr , it is possible to discriminate the SO 4 2− load released respectively by gypsum (CaSO 4 ) dissolution and pyrite (FeS 2 ) oxidation. Gypsum does not consume CO 2 during dissolution, while pyrite releases 2 mol CO 2 for each SO 4 2− ion released. Later, the C b (molar flux of Ca 2+ and Mg 2+ released by carbonates dissolution) is computed using R sil to account for the quantity of Ca 2+ and Mg 2+ derived from the weathering of the carbonates; a 60% of Ca 2+ is linked to calcite (CaCO 3 ) dissolution while the remaining 40% is released by dolomite (CaMgCO 3 ) dissolution. The dissolution of calcite and dolomite also takes up to 1 mol of CO 2 . The remaining Ca 2+ and Mg 2+ fluxes are linked to silicates that consume 2 mol CO 2 . A further description of the model calculation is found in Amiotte-Suchet 54 , Amiotte Suchet and Probst 37 and Donnini et al. 32 . The R pyr and R sil values used in the present study are compiled in the Extended Data Table 1.
Two 30-year periods are considered in the present study: a Historical (1969-1999) and a ECP (2069-2099), both computed at the daily time scale. The potential impacts derived from climate change are computed by quantifying the relative difference in the ECP period to the Historical period, following Eq. (3). Where ΔCO 2 denotes the relative change for each basin (expressed in %), and L CO2 represents the mean amount of CO 2 consumed through CW of rocks during the ECP and Historical period, expressed in Mg C year −1 . This result analysis is divided into two steps: at the annual scale, where L CO2 is the mean annual amount of CO 2 consumed in each period; and at the season scale, where L CO2 symbolises the mean amount of CO 2 consumed during all the springs, summers, autumns, and winters of each period.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.