# Rock organic carbon oxidation CO2 release offsets silicate weathering sink

Jesse Zondervan, Robert Hilton, Mathieu Dellinger, Fiona Clubb, Tobias

Roylands, Mateja Ogrič

► To cite this version:

Jesse Zondervan, Robert Hilton, Mathieu Dellinger, Fiona Clubb, Tobias Roylands, et al.. Rock organic carbon oxidation CO2 release offsets silicate weathering sink. Nature, 2023, 623, pp.329-333. 10.1038/s41586-023-06581-9 . hal-04242697

HAL Id: hal-04242697

https://hal.science/hal-04242697v1

Submitted on 22 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Title: Rock organic carbon oxidation CO₂ release offsets silicate weathering sink

**Authors:** Jesse R. Zondervan¹†*, Robert G. Hilton¹*, Mathieu Dellinger², Fiona Clubb³, Tobias Roylands³, Mateja Ogrič³

## Affiliations:

¹Department of Earth Sciences, University of Oxford; Oxford, OX1 3AN, United Kingdom

²EDYTEM-CNRS-University Savoie Mont Blanc (USMB), Chambéry, 73000, France

³Department of Geography, Durham University; Durham, DH1 3LD, United Kingdom

*Corresponding authors. Email: j.zondervan@ucl.ac.uk, robert.hilton@earth.ox.ac.uk

†Now at: Department of Earth Sciences, University College London; London, WC1E 6BT, United Kingdom

## Summary paragraph

Mountain uplift and erosion have regulated the balance of carbon between Earth's interior and atmosphere, where prior focus has been placed on the role of silicate mineral weathering in CO₂ drawdown and keeping Earth's climate in a habitable state ¹⁻⁵. However, weathering can release CO₂ as rock-organic carbon is oxidised at the near surface⁶,⁷ and this important geological CO₂ flux has remained poorly constrained³,⁸. Here we use the trace element rhenium to quantify this flux across global river catchments³,⁹, and combine this with a spatial extrapolation model. We find a CO₂ release of 68⁺¹⁸⁻⁶ megatons of carbon from weathering of organic carbon in near-surface rocks annually, rivalling or even exceeding the CO₂ drawdown by silicate weathering at the global scale¹⁰. Hotspots of CO₂ release are found in mountain ranges with high uplift rates exposing fine-grained sedimentary rock, such as the eastern Himalayas, the Rocky Mountains, and the Andes. Our results demonstrate that rock organic carbon is far from inert and causes weathering in regions to be either net sources or sinks of CO₂. We propose that this calls into question how erosion and weathering drive the long-term carbon cycle and contribute to the fine balance of carbon fluxes between the atmosphere, biosphere and lithosphere²,¹¹.
31

## Main

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

The tectonic activity that builds mountains can result in the uplift and exposure of organic carbon in rocks (OC$_{petro}$) alongside silicate mineral phases. The OC$_{petro}$ represents carbon stored in rocks that has accumulated over millions of years, previously sequestered from the atmosphere by photosynthesis and buried in sedimentary basins¹². Indeed, sedimentary and metasedimentary lithologies presently dominate bedrock geology of Earth surface, occupying ~64% of Earth surface¹³ and have OC concentrations of ~0.25 to > 1.0 weight percent by mass (%), whereas igneous rocks have much lower values (~0% or in the case of some marine basalts < 0.1%)¹⁵. Denudation supplies OC$_{petro}$ to the surface through physical plus chemical weathering³,¹⁵, and varies with rock type, relief, tectonic uplift, climate, and vegetation¹⁶, ¹⁷. Previous work has revealed the presence of OC$_{petro}$ in soils and rivers⁶,¹⁸⁻²⁰ and quantified the erosion of unweathered OC$_{petro}$¹⁴ and its global flux using data from the solid load of rivers at $43_{-25}^{+61}$ MtC yr$^{-1}$ ¹⁴,¹⁹. In contrast, estimated global rates of OC$_{petro}$ oxidation and CO$_2$ release during rock weathering currently derive from carbon cycle mass balance arguments, or ballpark upscaling of global river trace element fluxes⁸, with a range of estimates from 38 MtC yr$^{-1}$ ²¹ to 100 MtC yr$^{-1}$ ²². The uncertainty on OC$_{petro}$ oxidation fluxes is highlighted by recent work that cites a potential overall range for CO$_2$ release of 0 – 300 MtC yr$^{-1}$ ²³.

49

50

51

52

53

54

55

56

Consequently, to determine the role of rock weathering in the carbon cycle we require a robust, global quantification of OC$_{petro}$ oxidation over Earth's surface. Here we combine a compilation of OC$_{petro}$ oxidation proxy data from dissolved rhenium in well-studied catchments around the world; new probabilistic models of global OC$_{petro}$ stock and denudation; and a spatially explicit OC$_{petro}$ oxidation model with quantified uncertainty. This approach derives a global flux by extrapolating proxy derived OC$_{petro}$ oxidation data, whilst accounting for sampling bias across variables such as denudation rate and underlying geology.

57

## Rhenium as an OC<sub>petro</sub> oxidation proxy

58

59

60

61

62

63

64

65

66

The development of the trace element rhenium (Re) as a proxy to trace the oxidation of OC$_{petro}$ across landscapes²⁴,²⁵ has been underpinned by (i) the link between OC accumulation in marine sediments and organic matter being a host of Re ²⁶,²⁷; (ii) the paired loss of Re and OC$_{petro}$ during weathering of sedimentary rocks⁷,²⁵,²⁸; and (iii) the geochemical behaviour of Re that means it is exported as a dissolved oxyanion²⁹, flushed from a near-surface, oxidative weathering zone²⁵. In addition, studies tracking the fate of carbon released from the lithosphere during rock organic carbon weathering have found it can: directly enter the atmosphere as CO$_2$³⁰,³¹ or first dissolve as inorganic carbon in water³², and some can be incorporated into microbial biomass⁸.

67

68

69

70

71

In this study we compile published estimates of OC$_{petro}$ oxidation using the dissolved Re proxy, supplemented with new estimates derived from published dissolved Re concentrations⁷,⁹,²⁴,³³ (see Methods). A forward mixing model is used to quantify the proportion of dissolved Re from OC$_{petro}$ oxidation using ion ratios³⁴,²⁴, while constraints on organic carbon to rhenium contents of rocks being weathered ([$OC_{petro}$]/Re) come from new
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Spatial patterns of global OC<sub>petro</sub> stock and oxidation. (a) Locations of Re proxy samples and their upstream catchments. (b) Spatially explicit estimates of OC<sub>petro</sub> stocks in the upper 1 m of bedrock. (c) Spatial model of rock denudation derived from <sup>10</sup>Be data and a global raster of topographic slope. (d) Rock organic carbon oxidation fluxes extrapolated by our calibrated spatial model over the global surface.</figcaption></figure>

and published measurements (Methods, Supplementary Table S3). Our compilation includes 59 river basins, covering a range of drainage areas (50 - 5,900,000 km²), denudation rates and
climate regimes (Fig. 1), excluding river basins with high Re pollution levels such as the Danube, Yangtze and Mississippi⁹ (see Methods). The total OC$_{petro}$ weathering flux constrained from the Re-proxy across the drainage area of river basins in the dataset is 18 MtC yr$^{-1}$ (17 – 23 MtC yr$^{-1}$ within one standard deviation). The river basins in this study cover 18% of the Earth’s continental surface, and this flux would thus scale to $98^{+28}_{-9}$ MtC yr$^{-1}$ globally. However, OC$_{petro}$ stocks are spatially heterogeneous. To obtain a robust representative total OC$_{petro}$ weathering flux, in the next section, we use a spatial extrapolation model that considers patterns in OC$_{petro}$ stock and denudation rates.

## Distribution of OC$_{petro}$ availability

We spatially quantify the carbon stock and weathering flux of OC$_{petro}$ at Earth’s surface using a data-driven modelling approach. Our model incorporates topographic and lithological factors to estimate OC$_{petro}$ stocks, denudation rates and oxidative weathering rates, and is calibrated using our Re-proxy compilation (Supplementary Table S1). Unlike silicate weathering, which quickly becomes kinetically limited with increasing mineral supply by denudation³⁵, OC$_{petro}$ weathering appears to be predominately a supply-limited process⁸. This is reflected in oxidation rates which scale with erosion up to some of the highest erosion rates found on Earth, such as Taiwan and the European Alps⁷,²⁵. Recent work at the rock outcrop scale has shown that temperature and hydrology can control OC$_{petro}$ oxidation and CO$_2$ release over time in locations with very high rates of denudation³⁰,³¹. However, whilst the spatial control of denudation rates is well demonstrated on inter-catchment OC$_{petro}$ oxidation rates⁷,²⁵ our spatial catchment-scale Re-proxy compilation does not express other environmental controls (Methods). While temperature and hydrology controls likely operate, their spatial predictive power is small based on the available data. Here, oxidative weathering is modelled at a 1-km² grid scale, resolving at the scale of catchments constrained by the Re proxy (Supplementary Table S1).

The flux of CO$_2$ release by OC$_{petro}$ oxidation, $J_{ox}$ (mass × length$^{-2}$ × time$^{-1}$), can be expressed by a mass-balance of the form:

$$J_{ox} = \varepsilon \cdot \rho \cdot [OC_{petro}] \cdot \chi \qquad (\text{Eq. 1})$$

where $\varepsilon$ (length × time$^{-1}$) is the denudation rate, $\rho$ is rock density (mass × length$^{-3}$), [OC$_{petro}$] is the OC concentration in rock (mass × mass$^{-1}$) and $\chi$ is the weathering intensity as the fraction of OC weathered from rock. Weathering intensity $\chi$ has been shown to vary between low values of 0.2 in highly erosive settings⁷ and very high values of 0.98 in slow denudation settings⁸ with most falling in a range of 0.6 – 0.9 ⁷,³³,³⁴. Thus, $\chi$ presents a significantly smaller variance across environments in contrast to denudation rates and [OC$_{petro}$] which vary spatially by more than four orders of magnitude.

To constrain the stock of OC$_{petro}$ in the near surface we use [OC$_{petro}$] from the US Geological Survey rock geochemical database, combined with global lithological maps¹³ and spatial chemical lithology classifications³⁶. Our geospatial model simulates a large global near-surface OC$_{petro}$ stock, with the estimate and its interquartile range at $1490^{+2580}_{-980}$ Gt OC$_{petro}$ in the first meter of bedrock. This estimate is consistent with a global estimate of 1100 Gt OC$_{petro}$ within the first meter of sedimentary rocks¹⁴, a reassessment of deep soil radiocarbon data which provides evidence for OC$_{petro}$ inputs²⁰, and is of comparable magnitude to that of global soils (2060 ± 215 Gt OC)³⁷ and marine sediments (2322 ± 75 Gt OC)³⁸. As opposed to soil OC stocks, the distribution of OC$_{petro}$ is primarily controlled by the geological history of continents. Whilst the highest [OC$_{petro}$] is found in black shales (Extended Data Fig. 4), such rocks compose a tiny fraction of the Earth’s surface¹³, and instead, most OC$_{petro}$ is found in fine-grained sedimentary deposits such as shales (Fig. 1b). Geospatial patterns reveal low
OC$_{petro}$ stocks on the African continent (Fig. 1b), owing to a low occurrence of fine-grained sedimentary rocks. In contrast, substantial portions of Eurasia, South America, and the middle of North America east of the Rocky Mountains contain shales. The overlap of OC$_{petro}$ stocks and patterns of denudation driven mostly by rock uplift in mountains determines the exposure of this OC stock to oxidative weathering. We estimated denudation using a probabilistic spatial model that incorporates catchment-scale cosmogenic radionuclide (CRN) denudation rates³⁹, digital topography⁴⁰, and lithological maps¹³. The resultant modelled global denudation rate and its interquartile range is 11$^{+13}_{-6}$ Gt yr$^{-1}$, within range of recent estimates of global denudation at 28$^{+64}_{-20}$ Gt yr$^{-1}$ and 15 ± 2.8 Gt yr$^{-1}$.

## Spatial model of OC<sub>petro</sub> oxidation

Rather than the more classical mean-field parameterization schemes previously employed to model OC<sub>petro</sub> supply rates¹⁴, we use a probabilistic approach⁴¹ that accounts for the uncertainties in both variables in a spatial model. In each cell, empirical probability distributions of [OC<sub>petro</sub>] based on rock type (Extended Data Fig. 4) and probability distributions of denudation based on rock type and topographic slope (Extended Data Fig. 6) are sampled in 10,000 Monte Carlo simulations. We calibrated the geospatial model by minimizing the residuals between the modelled cell values of OC<sub>petro</sub> oxidation rates (J<sub>ox</sub>) and our compilation of Re-proxy data at the river catchment scale (Methods). Thus, our approach describes the spatial patterns of oxidative weathering rate as a function of topographic slope and rock type, which leads to simulations that are consistent with an assessment of global rock nitrogen weathering patterns which are dominated by denudation of fine-grained sedimentary rocks⁴¹.

Using our spatial model, we estimate that oxidative weathering of OC<sub>petro</sub> releases 68$^{+18}_{-6}$ MtC yr$^{-1}$ as CO₂ from the land-surface environment. The flux is lower than our spatially uncorrected extrapolation of Re-proxy measurements (98$^{+28}_{-9}$ MtC yr$^{-1}$), consistent with the slight bias towards high denudation rate settings in the river basin dataset. The best estimate of the oxidative weathering flux is higher than an independent estimate of OC<sub>petro</sub> erosion and river transport (i.e., the export of OC<sub>petro</sub> that has not been weathered) of 43$^{+61}_{-25}$ MtC yr$^{-1}$ based on river solid load composition and flux¹⁹, even though the uncertainties overlap. While a direct comparison of these estimates is difficult based on their quantification from dissolved vs particulate river chemistry and flux, they suggest an average weathering intensity of OC<sub>petro</sub> of ~60%, which is consistent with studies from large river basins¹⁹ and intensities measured in soils⁸,⁴².

The OC<sub>petro</sub> oxidation model can be used to estimate the turnover time of rock organic carbon at the surface. When combined with OC<sub>petro</sub> stocks, the model suggests that 0.05$^{+0.12}_{-0.03}$ % yr$^{-1}$ of the global OC stock in the first 10 cm of bedrock may be oxidized during denudation and weathering. A global OC<sub>petro</sub> loss rate of ~0.05 % yr$^{-1}$ equates to a carbon turnover time (the ratio of total OC<sub>petro</sub> to carbon outputs by oxidation) of ~2,000 years. This is about double the value of global soils⁴³, but shorter than turnover times in tundras of ~3,900 years⁴⁴. Given the large stock of OC<sub>petro</sub> that we report (~150 Pg C in the upper 10 cm) and its turnover time, it appears that OC<sub>petro</sub> cannot be assumed to be inert and passive in the shallow subsurface. The input of OC<sub>petro</sub> into soils can also impact soil residence time estimations and can lead to an underestimation of soil carbon exchange fluxes with the atmosphere²⁰.

Across the land surface, OC<sub>petro</sub> weathering is relatively focused (Fig. 1d), with variations in rock type and relief, which drive OC<sub>petro</sub> content and denudation respectively, determining the magnitude of OC<sub>petro</sub> oxidation and CO₂ release. Large regions of the African continent have lower OC stocks in bedrock and have lower relief, which together limit OC weathering. In
contrast, higher OC$_{petro}$ oxidation rates are estimated for northern latitudes, where OC-rich rock and high-relief landscapes are more prevalent. Overall, 10% of the Earth surface with the highest OC$_{petro}$ oxidation rates account for 60% of the global flux in our model. The world average rate is 0.5 tC km$^{-2}$ yr$^{-1}$, hotspots (10 times world average) and hyperactive areas (all areas >5 times world average) are responsible for 32% and 44% of CO$_2$ emissions, while representing only 1.2% and 3% of ice-free terrestrial land area, respectively. OC$_{petro}$ weathering rates in our model are more spatially concentrated than a 1-km resolution spatial model of silicate weathering$^{47}$, where hotspots (0.51% by area) and hyperactive areas (2.9% by area) accounted for 8.6% and 19.6% of total CO$_2$-consumption. This outcome appears reasonable because OC$_{petro}$ is less common spatially than silicate minerals and react faster$^{3,25}$.

## Weathering CO$_2$ sources versus sinks

The OC$_{petro}$ weathering flux and release of CO$_2$ to the atmosphere of $68_{-6}^{+18}$ MtC yr$^{-1}$ is similar to global terrestrial CO$_2$ uptake by silicate weathering (94 - 143 MtC yr$^{-1}$)$^{10}$. Silicate weathering involves dissolved and gaseous CO$_2$ uptake through bicarbonate production and the release of dissolved ions, some of which then precipitate as marine carbonate rocks$^4$. The resultant total transfer of carbon from the atmosphere to the lithosphere by silicate weathering is 47 - 72 MtC year$^{-1}$. In addition to their opposing impacts on the transfer of carbon between the atmosphere and lithosphere, fluxes of silicate weathering vs OC$_{petro}$ oxidation may have contrasting responses to climate. Silicate weathering is invoked as a negative feedback to climate warming through increased rates of silicate weathering from increased temperature and a more vigorous hydrological cycle, drawing down more CO$_2$$^{35,46}$. In contrast, in high denudation rate settings the CO$_2$ release from OC$_{petro}$ oxidation may increase with temperature$^{30,31}$, while links to glacial erosion processes complicate the feedback between oxidative weathering and climate$^{33}$.

Silicate and OC$_{petro}$ weathering processes may overlap, as sedimentary rocks contain silicate minerals as well as OC. However, the relative magnitude of these fluxes will vary spatially with climate, rock type and denudation$^{35,46}$. We assess the net balance of rock weathering within major river basins (Fig. 2), using our OC$_{petro}$ oxidation model and silicate weathering estimates$^{10}$.

Within uncertainties, rock weathering in about a third of the major river basins is a net source of CO$_2$ after OC$_{petro}$ oxidation is considered, even while using the values of initial atmospheric CO$_2$ consumption of silicate weathering rather than the smaller quantity of CO$_2$ eventually locked up in the lithosphere through carbonate precipitation of the associated released dissolved ions (Fig. 2a). The Yangtze (Changjiang) and Mekong draining the eastern flanks of the Himalayas, and the Mackenzie River draining shales west of the Rockies in Canada are major sources of CO$_2$ from rock weathering. These high-emitting basins have in common some of the highest
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Earth's major river basins, their carbon sinks, and sources from silicate weathering¹⁰ and OCpetro weathering (this study), and their overall rock weathering budget based on these fluxes. (a) OCpetro vs silicate weathering fluxes and their net values. Basins that produce a net source of CO₂ are shown in the shaded half of the plot, with the net magnitude of the weathering CO₂ flux illustrated by the symbol color (in MtC yr⁻¹). (b) Net weathering balance vs basin average denudation (red arrow: cross-over at ~30 mm kyr⁻¹) and (c) vs basin average OCpetro stock. Error bars represent uncertainty of OCpetro oxidation model outputs based on the uncertainty of the training data (see section 1.2.3). Variable rates of uplift and erosion, climate and OCpetro stocks across Earth's surface impact OCpetro and silicate weathering rates differently, leading to regions where rock weathering is (d) a source or (e) sink of CO₂.</figcaption></figure>

basin-average denudation rates and OCpetro stocks (Fig. 2b-c), which is consistent with OCpetro oxidation being driven by OCpetro stocks and denudation (Eq. 1; 7,8,25).

Hotspots of CO₂ release during rock weathering appear to lie at the edges of major active mountain ranges where relatively young, marine sedimentary deposits are uplifted and
supplied to oxygen through denudation. Examples include the shales of the Himalayan collision zones and east of the Rocky Mountains (Fig. 1b, 1d). On the other hand, basins where rock weathering is the biggest net sink of CO₂ do not necessarily lie at the extremes of low denudation and low OCpetro stocks. Whilst the tropical Congo River and volcanics-dominated Godavari River basins have low basin-average denudation rates and low OCpetro stocks, the biggest weathering sink, the Amazon River, lies in the global middle range of denudation rates and OCpetro stocks (Fig. 2b, 2c). There, the kinetically limited silicate weathering reaction benefits from long sediment residence times and a warm, humid climate. While the Andes is a hotspot for OCpetro oxidation fluxes (Figure 1d), the exceptionally large lowland drainage area of the Amazon means that OCpetro oxidation may be supply limited. In a third of river basins weathering remains carbon neutral within uncertainty such as the volcanic-rich Columbia River catchment.

To avoid large swings in atmospheric CO₂ over millions of years and maintain an apparent close balance of CO₂ sources and sinks²,¹¹, any potential imbalances in weathering-derived carbon fluxes must be addressed by accounting for other components in the long-term carbon cycle. Solid-Earth degassing associated with volcanoes and diffuse release from metamorphism in subduction zones is responsible for 79 ± 9 MtC yr⁻¹ released into the atmosphere (Fig. 3a; ⁴⁷), while any additional (non-subduction) global CO₂ release during orogenic metamorphism and sulfide oxidation and inorganic C uptake during seafloor weathering are more poorly constrained³. As our results show that the weathering of OCpetro offsets silicate weathering in the long-term carbon cycle, a large additional sink of CO₂ is needed. This may be provided by burial of organic matter in ocean sediments, which could contribute as much as 170 MtC yr⁻¹ (Figure 3b; ⁴⁸). In addition, as OCpetro fluxes can overtake silicate weathering during periods of more intense uplift and erosion (Fig. 2b,d; Extended Data Fig. 8), the question whether orogenic periods in Earth history are sources or sinks of atmospheric CO₂ is now a reopened question³,³¹,⁴⁹,⁵⁰. The net balance will depend on factors such as transport of terrestrial biospheric carbon to oceans (157±74 MtC yr⁻¹) and its burial¹⁹. A global comparison of catchment-scale OCpetro oxidation yields and estimated terrestrial biospheric OC burial (Extended Data Figure 8) suggest the OC burial can apparently offset or even overcompensate CO₂ release from OCpetro oxidation. This understanding persists when the additional marine OC burial sink in sediment is factored into global flux estimates (Fig. 3b; ⁴⁸). The dynamics of Earth's thermostat thus need to be revisited to account for variation in all these fluxes and consider how their relative importance may have changed as life evolved and the OC stocks of sedimentary rocks have increased³,²².
<figure><img src="image_4.png" /><figcaption>**Figure 3.** A shift in understanding the geological carbon cycle. (a) The inorganic geological carbon cycle relies on a global balance between solid-Earth CO₂ degassing and silicate weathering. (b) The emerging understanding of the role of organic carbon in the global geological carbon cycle, supported by the high flux of OCpetro oxidation reported in this study. Hence the biological, chemical, and physical processes of biospheric organic carbon production, burial and release control long-term climate variability and stability.</figcaption></figure>

**Data availability**

All dissolved rhenium sample data are available in the Supplementary Materials Tables S1-S5, which in addition to geospatial data including those for Figure 1 are deposited in a Zenodo repository (https://doi.org/10.5281/10.5281/zenodo.8144244). Source data for Figure 2 and Extended Data Figures 4–6 are also available in separate Excel and CSV files respectively (see Figure legends).

**Code availability**

The spatial OCpetro oxidation model contains two spatial probabilistic subroutines of (i) OCpetro stocks in surface rocks; and (ii) spatially defined denudation rates. These are combined in a Monte Carlo framework alongside the Re-proxy river catchment data to
optimize the model and then extrapolate OC$_{petro}$ oxidation rates (Extended Data Fig. 1).

Model simulations were implemented at 1-km grid scale (Mollweide projection, WGS84 datum) in the Python programming language. R and python code, their environments and the necessary data files to run these are all deposited in a Zenodo depository
(https://doi.org/10.5281/zenodo.8144244) and additionally code is available from Github:
https://github.com/jessezondervan/Global OCpetro Oxidation.

References

1. Chamberlin, T. C. An Attempt to Frame a Working Hypothesis of the Cause of Glacial Periods on an Atmospheric Basis. *The Journal of Geology* **7**, 545-584 (1899).
https://doi.org/10.1086/608449

2. Berner, R. A., Caldeira, K. The need for mass balance and feedback in the geochemical carbon cycle. *Geology* **25**, 955-956 (1997). https://doi.org/10.1130/0091-7613(1997)025<0955:Tnfmba>2.3.Co;2

3. Hilton, R. G., West, A. J. Mountains, erosion and the carbon cycle. *Nature Reviews Earth & Environment* **1**, 284-299 (2020). https://doi.org/10.1038/s43017-020-0058-6

4. Gaillardet, J., Dupré, B., Louvat, P., Allègre C. J. Global silicate weathering and CO2 consumption rates deduced from the chemistry of large rivers. *Chemical Geology* **159**, 3-30 (1999). https://doi.org/10.1016/S0009-2541(99)00031-5

5. Raymo, M. E., Ruddiman, W. F. Tectonic forcing of late Cenozoic climate. *Nature* **359**, 117-122 (1992). https://doi.org/10.1038/359117a0

6. Galy, V., Beyssac, O., France-Lanord, C., Eglinton, T. Recycling of graphite during Himalayan erosion: A geological stabilization of carbon in the crust. *Science* **322**, 943-945 (2008). https://doi.org/10.1126/science.1161408

7. Hilton, R. G., Gaillardet, J., Calmels, D., Birck, J.-L. Geological respiration of a mountain belt revealed by the trace element rhenium. *Earth and Planetary Science Letters* **403**, 27-36 (2014). https://doi.org/10.1016/j.epsl.2014.06.021

8. Petsch, S. T. in *Treatise on Geochemistry* (eds Holland, H. D., Turekian, K. K.) 217-238 (Elsevier, 2014).

9. Miller, C.A., Peucker-Ehrenbrink, B., Walker, B.D, Marcantonio, F. Re-assessing the surface cycling of molybdenum and rhenium. *Geochimica et Cosmochimica Acta* **75**, 7146-7179 (2011). https://doi.org/10.1016/j.gca.2011.09.005

10. Moon, S., Chamberlain, C. P., Hilley, G. E. New estimates of silicate weathering rates and their uncertainties in global rivers. *Geochimica et Cosmochimica Acta* **134**, 257-274 (2014). https://doi.org/10.1016/j.gca.2014.02.033

11. Derry, L.A. Carbonate weathering, CO₂ redistribution, and Neogene CCD and pCO₂ evolution. *Earth and Planetary Science Letters* **597**, 117801 (2022).
https://doi.org/10.1016/j.epsl.2022.117801

12. Derry, L.A., France-Lanord, C. Neogene growth of the sedimentary organic carbon reservoir. *Paleoceanography* **11**, 267-275 (1996). https://doi.org/10.1029/95PA03839

13. Hartmann, J., Moosdorf, N. The new global lithological map database GLiM: A representation of rock properties at the Earth surface. *Geochemistry, Geophysics, Geosystems* **13**, Q12004 (2012). https://doi.org/10.1029/2012GC004370
327

328

329

14. Copard, Y., Amiotte-Suchet, P., Di-Giovanni, C. Storage and release of fossil organic carbon related to weathering of sedimentary rocks. *Earth and Planetary Science Letters* **258**, 345-357 (2007). https://doi.org/10.1016/j.epsl.2007.03.048

330

331

332

15. Bolton, E.W., Berner, R.A., Petsch, S.T. The weathering of sedimentary organic matter as a control on atmospheric O₂: II. Theoretical modeling. *American Journal of Science* **306**, 575-615 (2006). https://doi.org/10.2475/08.2006.01

333

334

16. Larsen, I. J., Montgomery, D. R., Greenberg, H. M. The contribution of mountains to global denudation. *Geology* **42**, 527-530 (2014). https://doi.org/10.1130/g35136.1

335

336

337

17. Wittmann, H., Oelze, M., Gaillardet, J., Garzanti, E., von Blanckenburg, F. A global rate of denudation from cosmogenic nuclides in the Earth's largest rivers. *Earth-Science Reviews* **204**, 103147 (2020). https://doi.org/10.1016/j.earscirev.2020.103147

338

339

340

341

18. Blair, N. E., Leithold, E. L., Ford, S. T., Peeler, K. A., Holmes, J. C., Perkey, D. W. The persistence of memory: The fate of ancient sedimentary organic carbon in a modern sedimentary system. *Geochimica et Cosmochimica Acta* **67**, 63-73 (2003). https://doi.org/10.1016/S0016-7037(02)01043-8

342

343

344

19. Galy, V., Peucker-Ehrenbrink, B., Eglinton, T. Global carbon export from the terrestrial biosphere controlled by erosion. *Nature* **521**, 204-207 (2015). https://doi.org/10.1038/nature14400

345

346

347

20. Grant, K.E., Hilton, R.G., Galy, V.V. Global patterns of radiocarbon depletion in subsoil linked to rock-derived organic carbon. *Geochem. Persp. Let.* **25**, 36–40 (2023). https://doi.org/10.7185/geochemlet.2312

348

349

350

21. Garrels, R.M., Lerman, A. Coupling of the sedimentary sulfur and carbon cycles; an improved model. *American Journal of Science* **284**, 989-1007 (1984). https://doi.org/10.2475/ajs.284.9.989

351

352

353

22. Kump, L.R., Arthur, M.A. Interpreting carbon-isotope excursions: carbonates and organic matter. *Chemical Geology* **161**, 181-198 (1999). https://doi.org/10.1016/S0009-2541(99)00086-8

354

355

356

23. Regnier, P., Resplandy, L., Najjar, R. G., Ciais, P. The land-to-ocean loops of the global carbon cycle. *Nature* **603**, 401-410 (2022). https://doi.org/10.1038/s41586-021-04339-9

357

358

359

360

361

24. Dalai, T. K., Singh, S. K., Trivedi, J. R., Krishnaswami, S. Dissolved rhenium in the Yamuna river system and the Ganga in the Himalaya: role of black shale weathering on the budgets of Re, Os, and U in rivers and CO₂ in the atmosphere. *Geochimica et Cosmochimica Acta* **66**, 29-43 (2002). https://doi.org/10.1016/S0016-7037(01)00747-5

362

363

364

365

366

25. Hilton, R. G., Turowski, J. M., Winnick, M., Dellinger, M., Schleppi, P., Williams, K. H., Lawrence, C. R., Maher, K., West, M., Hayton, A. Concentration-Discharge Relationships of Dissolved Rhenium in Alpine Catchments Reveal Its Use as a Tracer of Oxidative Weathering. *Water Resources Research* **57**, e2021WR029844 (2021). https://doi.org/10.1029/2021WR029844

367

368

369

26. Selby, D., Creaser, R. A. Re-Os geochronology of organic rich sediments: an evaluation of organic matter analysis methods. *Chemical Geology* **200**, 225-240 (2003). https://doi.org/10.1016/S0009-2541(03)00199-2
370

371

372

27. Morford, J.L., Martin, W.R., Carney, C.M. Rhenium geochemical cycling: Insights from continental margins. *Chemical Geology* **324**, 73-86 (2012).
https://doi.org/10.1016/j.chemgeo.2011.12.014

373

374

375

376

28. Jaffe, L. A., Peucker-Ehrenbrink, B., Petsch, S. T. Mobility of rhenium, platinum group elements and organic carbon during black shale weathering. *Earth and Planetary Science Letters* **198**, 339-353 (2002). https://doi.org/10.1016/S0012-821X(02)00526-5

377

378

379

29. Colodner, D., Sachs, J., Ravizza, G., Turekian, K., Edmond, J., Boyle, E. The geochemical cycle of rhenium: a reconnaissance. *Earth and Planetary Science Letters* **117**, 205-221 (1993). https://doi.org/10.1016/0012-821X(93)90127-U

380

381

382

383

30. Soulet, G., Hilton, R. G., Garnett, M. H., Roylands, T., Klotz, S., Croissant, T., Dellinger, M., Le Bouteiller, C. Temperature control on CO₂ emissions from the weathering of sedimentary rocks. *Nature Geoscience* **14**, 665-671 (2021).
https://doi.org/10.1038/s41561-021-00805-1

384

385

386

387

388

31. Roylands, T., Hilton, R. G., Garnett, M. H., Soulet, G., Newton, J. A., Peterkin, J. L., Hancock P. Capturing the short-term variability of carbon dioxide emissions from sedimentary rock weathering in a remote mountainous catchment, New Zealand. *Chemical Geology* **608**, 121-204 (2022).
https://doi.org/10.1016/j.chemgeo.2022.121024

389

390

391

392

32. Tune, A.K., Druhan, J.L., Wang, J., Bennett, P.C. and Rempe, D.M. Carbon dioxide production in bedrock beneath soils substantially contributes to forest carbon cycling. *Journal of Geophysical Research: Biogeosciences*, **125**, e2020JG005795 (2020).
https://doi.org/10.1029/2020JG005795

393

394

395

33. Horan, K., Hilton, R.G., Selby, D., Ottley, C.J., Gröcke, D.R., Hicks, M., Burton, K.W. Mountain glaciation drives rapid oxidation of rock-bound organic carbon. *Science Advances* **3**, e1701107 (2017). https://doi.org/10.1126/sciadv.1701107

396

397

398

399

400

34. Horan, K., Hilton, R.G., Dellinger, M., Tipper, E., Galy, V., Calmels, D., Selby, D., Gaillardet, J., Ottley, C.J., Parsons, D.R., Burton, K.W. Carbon dioxide emissions by rock organic carbon oxidation and the net geochemical carbon budget of the Mackenzie River Basin. *American Journal of Science* **319**, 473-499 (2019).
https://doi.org/10.2475/06.2019.02

401

402

403

35. West, A. J. Thickness of the chemical weathering zone and implications for erosional and climatic drivers of weathering and for carbon-cycle feedbacks. *Geology* **40**, 811-814 (2012). https://doi.org/10.1130/g33041.1

404

405

406

407

36. Amiotte Suchet, P., Probst, J.-L., Ludwig, W. Worldwide distribution of continental rock lithology: Implications for the atmospheric/soil CO₂ uptake by continental weathering and alkalinity river transport to the oceans. *Global Biogeochemical Cycles* **17**, 1038 (2003). https://doi.org/10.1029/2002GB001891

408

409

410

37. Batjes, N. H. Harmonized soil property values for broad-scale modelling (WISE30sec) with estimates of global soil carbon stocks. *Geoderma* **269**, 61-68 (2016). https://doi.org/10.1016/j.geoderma.2016.01.034

411

412

413

38. Atwood, T.B., Witt, A., Mayorga, J., Hammill, E., Sala, E. Global patterns in marine sediment carbon stocks. *Frontiers in Marine Science* **7**, 165 (2020).
https://doi.org/10.3389/fmars.2020.00165
414

415

416

39. Codilean, A.T., Munack, H., Cohen, T.J., Saktura, W.M., Gray, A., Mudd, S.M.
OCTOPUS: An open cosmogenic isotope and luminescence database. *Earth System
Science Data* **10**, 2123-2139 (2018). https://doi.org/10.5194/essd-10-2123-2018

417

418

419

420

40. Amatulli, G., McInerney, D., Sethi, T., Strobl, P., Domisch, S. Geomorpho90m,
empirical evaluation and accuracy assessment of global high-resolution
geomorphometric layers. *Scientific Data* **7**, 162 (2020).
https://doi.org/10.1038/s41597-020-0479-6

421

422

423

41. Houlton, B. Z., Morford, S. L., Dahlgren, R. A. Convergent evidence for widespread rock nitrogen sources in Earth's surface environment. *Science* **360**, 58-62 (2018).
https://doi.org/doi:10.1126/science.aan4399

424

425

426

427

42. Hemingway, J.D., Rothman, D.H., Grant, K.E., Rosengard, S.Z., Eglinton, T.I., Derry, L.A., Galy, V.V. Mineral protection regulates long-term global preservation of natural organic carbon. *Nature* **570**, 228-231 (2019). https://doi.org/10.1038/s41586-019-1280-6

428

429

430

431

43. Shi, Z., Allison, S.D., He, Y., Levine, P.A., Hoyt, A.M., Beem-Miller, J., Zhu, Q., Wieder, W.R., Trumbore, S., Randerson, J.T. The age distribution of global soil carbon inferred from radiocarbon measurements. *Nature Geoscience* **13**, 555-559 (2020). https://doi.org/10.1038/s41561-020-0596-z

432

433

434

44. Luo, Z., Wang, G., Wang, E. Global subsoil organic carbon turnover times
dominantly controlled by soil properties rather than climate. *Nature Communications*
**10**, 3688 (2019). https://doi.org/10.1038/s41467-019-11597-9

435

436

437

438

45. Hartmann, J., Jansen, N., Dürr, H.H., Kempe, S., Köhler, P. Global CO₂-consumption by chemical weathering: What is the contribution of highly active weathering regions? *Global and Planetary Change* **69**, 185-194 (2009).
https://doi.org/10.1016/j.gloplacha.2009.07.007

439

440

441

46. Walker, J. C. G., Hays, P. B., Kasting, J. F. A negative feedback mechanism for the long-term stabilization of Earth's surface temperature. *Journal of Geophysical Research: Oceans* **86**, 9776-9782 (1981). https://doi.org/10.1029/JC086iC10p09776

442

443

47. Plank, T., Manning, C. E. Subducting carbon. *Nature* **574**, 343-352 (2019).
https://doi.org/10.1038/s41586-019-1643-z

444

445

446

48. Smith, R. W., Bianchi, T. S., Allison, M. C. Savage, V. Galy, High rates of organic carbon burial in fjord sediments globally. *Nature Geoscience* **8**, 450-453 (2015).
https://doi.org/10.1038/ngeo2421

447

448

449

49. Evans, M.J., Derry, L.A. and France-Lanord, C. Degassing of metamorphic carbon dioxide from the Nepal Himalaya. *Geochemistry, Geophysics, Geosystems*, **9**, 4 (2008). https://doi.org/10.1029/2007GC001796

450

451

452

453

454

50. Märki, L., Lupker, M., France-Lanord, C., Lavé, J., Gallen, S., Gajurel, A.P., Haghipour, N., Leuenberger-West, F. and Eglinton, T. An unshakable carbon budget for the Himalaya. *Nature Geoscience*, **14**, 745-750 (2021).
https://doi.org/10.1038/s41561-021-00815-z
## Methods

The workflow of the materials and methods (Extended Data Fig. 1) starts with the compilation and derivation of rhenium-derived OC<sub>petro</sub> oxidation flux estimates (Section 1), before detailing the global spatial OC<sub>petro</sub> oxidation model (Section 2), its denudation and OC<sub>petro</sub> stock sub-models (Sections 2.1, 2.2) and the Monte Carlo routines used in the model approach (Section 2.4).

### Section 1. Rhenium-based river catchment estimates of OC<sub>petro</sub> oxidation

From a series of dissolved rhenium measurements (typically completed by ICP-MS), the dissolved Re flux J<sub>Re</sub> (t yr<sup>-1</sup>) can be used to estimate OC<sub>petro</sub> oxidation flux, J<sub>OCPetro-ox</sub> (tC yr<sup>-1</sup>) using:

$$J_{OCPetro-ox} = J_{Re} \times \left( \frac{[OC_{petro}]}{[Re]} \right)_i \times f_c \quad (\text{Eqn. S1})$$

where f<sub>c</sub> is the fraction of Re derived from OC<sub>petro</sub> oxidation<sup>34</sup> and ([OC<sub>petro</sub>]/[Re])<sub>i</sub> is the organic carbon to rhenium concentration ratio (g g<sup>-1</sup>) in rocks undergoing weathering. In some of the catchments where it may be important, an additional term, not shown in Equation S1, has been applied to correct for the presence of graphite which may not undergo alteration during weathering<sup>33</sup>.

### 1.1 Compiled published measurements

In this study we compile estimates of OC<sub>petro</sub> oxidation using the dissolved Re proxy from published literature (Supplementary Table S1). These include: the Yamuna River, India<sup>24</sup>; 10 Taiwanese Rivers<sup>7</sup>; 4 rivers from the western Southern Alps<sup>33</sup>; 4 rivers from the Mackenzie Basin, Canada<sup>3,34</sup>; and 6 rivers draining the Peruvian Andes<sup>51</sup>. Two Swiss catchments<sup>25</sup> are not included because of their very small catchment area compared to the geospatial scales over which we complete the upscaling.

For some of these case studies dissolved rhenium flux has been estimated from repeated sampling and discharge records<sup>34</sup>, while earlier studies all include single snapshot samples<sup>7,24</sup>, and all include measurements of the local sedimentary rock composition. Most of these compiled studies have used dissolved ion ratios to estimate the source of dissolved Re, i.e. f<sub>c</sub> (Eqn. S1), apart from the Taiwan dataset<sup>7</sup>. While uncertainties on the OC<sub>petro</sub> oxidation yields appear relatively large (Supplementary Table S1), it is important to note that the measured range in yields is much larger than the uncertainties.

### 1.2 New estimates of OC<sub>petro</sub> oxidation

To expand the 25 estimates of OC<sub>petro</sub> oxidation from river catchments described previously, we build on a previous study of dissolved Re fluxes in large rivers that reports dissolved Re concentrations and fluxes for major basins around the world<sup>9</sup>. We use these measurements and combine them with estimates of f<sub>c</sub> and ([OC<sub>petro</sub>]/[Re])<sub>i</sub>, discussed in the following sections, to calculate OC<sub>petro</sub> oxidation yields with associated uncertainties using published approaches<sup>34,25</sup>

In locations with significant local sources of fossil fuel combustion (e.g., coal-fired power plants or steel works), rainwater can contain concentrations of Re that approach those of river water<sup>8,52</sup>, whereas locations that have minimal impacts from local pollution sources have Re concentrations in rainwater that are below detection<sup>25,33</sup>. In the large river dataset<sup>9</sup> some large rivers are noted to have markedly increased Re concentrations and fluxes, with the conclusion that this was due to anthropogenic Re inputs. In a first order catchment of the Mississippi Basin, this has been confirmed by a detailed Re mass balance<sup>52</sup>. A study of Re
across Indian catchments suggests that whilst Re in Himalayan catchments and the mainstem Ganges and Brahmaputra behave conservatively, peninsular lower relief catchments with denser populations and industrial activity suggest anthropogenic inputs⁵³. For this study's purposes to quantify weathering reactions, we only use Himalayan rivers and the mainstem Ganges and Brahmaputra in India and we have further excluded Re data from the Danube, Mississippi, and Yangtze rivers from our analysis. Our addition of catchment Re data to the Miller dataset includes a large contribution of small upland catchments with higher average erosion rates, where the authors of these studies selected sites with minimal human disturbance (Supplementary Table S1). We further consider the role of anthropogenic Re in our model results in section 3.

### 1.2.1 Estimation of Re source and $f_c$

To estimate the fraction of dissolved Re sourced from OC$_{petro}$, $f_c$, for the rivers in the Re flux dataset⁹, we follow a previously used forward model mixing approach³⁴,²⁵:

$$Re_{OC} = Re_{tot} \times f_c = Re_{tot} - Re_{sulf} - Re_{sil} \quad (\text{Eqn. S2})$$

Where Re$_{OC}$ is the rhenium concentration of OC$_{petro}$-derived Re in the dissolved load, Re$_{tot}$ is the measured concentration, Re$_{sulf}$ and Re$_{sil}$ are the concentrations derived from weathering of sulfide and silicate minerals respectively. These unknowns can be quantified as:

$$Re_{sulf} = SO_4 \times \left(\frac{Re}{SO_4}\right)_{sulf} \quad (\text{Eqn. S3})$$

$$Re_{sil} = Na \times \left(\frac{Re}{Na}\right)_{sil} \quad (\text{Eqn. S4})$$

where the element ratios of the end members for silicate, $(Re/Na)_{sil}$, and sulfide $(Re/SO_4)_{sulf}$ are defined and with the assumption that the dissolved sulfate (SO$_4$) and sodium (Na) respectively derive only from sulfide oxidation and silicate weathering and are conservative. This returns an upper bound on the Re$_{sulf}$ and Re$_{sil}$ components (Supplementary Table S2). Following recent work³⁴, we use a range of values for each, where $(Re/SO_4)_{sulf}$ ranges from 0.2 x 10⁻³ pmol/µmol to 4 x 10⁻³ pmol/µmol (23), and $(Re/Na)_{sil}$ ranges from 0.4 x 10⁻³ pmol/µmol to 2 x 10⁻³ pmol/µmol. Here, we correct Na+ concentrations for atmospheric-derived Na, $[Na+]^*$, where $[Na+]^* = [Na+] - [Cl^-] \times 0.8$, assuming all Cl⁻ derives from precipitation, which has a molar $[Na+]/[Cl^-]$ ratio of 0.8. We similarly correct for atmospheric SO$_4$ inputs.

### 1.2.2 Constraints on ($[OC_{petro}]/[Re]_i$)

A recent compilation⁵⁴ provides measurements of $[OC_{petro}]/[Re]_i$ from rock samples of different ages around the world. However, most of these measurements were made on black shales, with OC$_{petro}$ contents > 1% which occur on only 0.3% of the Earth surface (Fig. S5). Riverbed material sediments from erosive catchments provide an alternative way to capture landscape-scale average rock composition, albeit with some potential for weathering to alter the primary signal. Here we compile measurements of $[OC_{petro}]/[Re]$ on bed materials from rivers around the world (Supplementary Table S3; Extended Data Fig.2) and supplement this dataset with additional samples from mudrocks of the Eastern Cape, New Zealand, and the Peruvian Andes measured using methods described in detail previously²⁵. We find that regions with lower OC$_{petro}$ concentrations that are more typical of sedimentary rocks at the
continental surface - units including fine-grained sedimentary rocks make up > 35 % of Earth surface (Extended Data Fig.5) - have lower and more consistent ratios of OC and Re in their rocks. The samples from the Peel River in the Mackenzie River basin³⁴, overlap the lower end of the published black shale values. Since this is the catchment with the highest proportion of black shales in our Re dataset, these samples allow us to capture the imprint of this important marginal lithology at the landscape-scale.

The bedrock composition in the catchments of rivers studied in the Re flux dataset⁹ are not reported. However, we note the good geographic coverage and number of samples that we have from riverbed materials from erosive settings around the world. These provide constraint on the initial organic carbon to rhenium ratio in the rocks. To conservatively quantify uncertainty in the range of OC<sub>petro</sub> oxidation rates from dissolved Re data we perform a Monte Carlo simulation in which we uniformly sample the entire range of measured [OC<sub>petro</sub>]/[Re] values, from low values indicative of carbon-poor and/or metamorphic rocks 2.5 x 10<sup>-8</sup> g g<sup>-1</sup>³³ towards 1.26 x 10<sup>-6</sup> ³⁴ in catchments with higher OC in rocks (Supplementary Table S3; Extended Data Fig.3).

### 1.2.3 OC<sub>petro</sub> oxidation yields and uncertainties

Equation 1 is used for each basin in the Re flux dataset⁹. Uncertainties in f<sub>c</sub> derive from the range of values used in the sulfide and silicate end member compositions (Eqns. S3, S4). For ([OC<sub>petro</sub>]/[Re]) we use the range of values discussed in Section 1.2.2 (Extended Data Fig.3).

A Monte Carlo uncertainty propagation is used on these variables, with 10,000 randomly selected combinations of input values (with uniform sampling) are used to estimate J<sub>OCpetro-ox</sub> for each basin. The median value of the Monte Carlo simulation and the interquartile range are reported (Supplementary Table S1).

## 1.3 Geospatial catchment boundaries

To derive the catchment outlines and areas of Re-proxy samples in our compiled dataset, we used the HYDROSHEDS flow direction grid at 3 arc-second resolution⁵⁵, and ArcGIS Pro⁵⁶.

For catchments in Iceland and outside the latitudinal cover of HYDROSHEDS, catchments were derived from ALOS AW3D using TauDEM functionality in OpenTopography⁵⁷ and from the HYDRO1K flow direction grid respectively⁵⁸. Whilst most published sample coordinates Supplementary Table S1) give the correct location on the cited drainage systems, in a handful of cases coordinates had to be amended by up to a few km, which may reflect errors in transcribing (e.g., Kikori and Purari⁹). Final quality control included a comparison of the extracted drainage basin areas and those published, with good agreement overall (< 2% residual). However, some drainage areas cited in the Re flux dataset⁹ refer to the river mouth, rather than the river catchment upstream of the Re sample location. In these cases, we use the Re sample location and its upstream catchment. Finalised coordinates of Re samples snapped to the drainage system are given in Supplementary Table S1, together with the Re locations' upstream drainage area. Spatial files of upstream drainage boundaries and Re sample locations are available on Zenodo (Available from: https://doi.org/10.5281/zenodo.8144244).

To convert dissolved Re concentrations into Re fluxes, average annual water discharge was calculated using published numbers at gauges (Supplementary Table S1) and scaled to the upstream drainage area of Re sample locations.

In addition to spatial catchment boundaries for the Re proxy dataset, we compare our spatial model output to published estimates of silicate weathering¹⁰ that use the GRDC dataset of Major River Basins of the World⁵⁹. Drainage areas used by¹⁰ have slight discrepancies with those found in the GRDC dataset. We account for these in our analysis of major river basin net weathering flux (Supplementary Table S4).
## Section 2. Global spatial OC<sub>petro</sub> oxidation model

The following section provides additional rationale and details of the modelling approaches. In brief, they contain two spatial probabilistic subroutines of (i) OC<sub>petro</sub> stocks in surface rocks; and (ii) spatially defined denudation rates. These are combined in a Monte Carlo framework alongside the Re-proxy river catchment data to optimize the model and then extrapolate OC<sub>petro</sub> oxidation rates (Extended Data Fig.1). Model simulations were implemented at 1-km grid scale (Mollweide projection, WGS84 datum) in the Python programming language<sup>60</sup>.

### 2.1 OC<sub>petro</sub> stocks

Rock samples from the USGS Rock Geochemical Database sorted into lithological categories (Supplementary Table S5) were mapped onto units of the highest-resolution global lithological maps currently available<sup>13</sup>. Extended Data Fig.4 (below) shows the OC<sub>petro</sub> concentration of key lithologies in the USGS Rock Geochemical database. Weight percentage values from the USGS Rock Geochemical Database were converted to OC<sub>petro</sub> stock, using rock densities (Supplementary Table S5). In our Monte Carlo framework, OC<sub>petro</sub> stocks at each grid cell were sampled independently using the empirical distributions of rock OC<sub>petro</sub> content derived from the USGS database (Extended Data Fig.4) and our unit classification (Supplementary Table S5). In our lithology model, complex mapped units present in GLiM consist of a combination of carbonates, silicates and various grainsizes (Extended Data Fig.5, Supplementary Table S5). To calculate the OC<sub>petro</sub> reservoir among these units, we derive the fractional abundance of lithology types (F<sub>n</sub>) from continental-scale literature estimates<sup>36</sup>:

$$ [OC_{petro}]_{rock} = F_1([OC_{petro}]_{lithology,1}) + F_2([OC_{petro}]_{lithology,2}) + \cdots + F_n([OC_{petro}]_{lithology,n}) \qquad (\text{Eqn. S5}) $$

$$ F_1 + F_2 + \cdots + F_n = 1 \qquad (\text{Eqn. S6}) $$

### 2.2 Denudation model

The denudation model is parameterised using a regression approach, similar to techniques applied elsewhere<sup>16,41</sup>. We regressed a compilation of long-term catchment scale <sup>10</sup>Be denudation estimates<sup>39</sup> against mean local slope generated from the Geomorpho90m dataset<sup>40</sup>. Mean local slope was calculated using the focal statistics tool in ArcGIS Pro<sup>56</sup> and the Geomorpho90m slope dataset with a 5-km moving radius. Slope values were then matched to <sup>10</sup>Be denudation estimates at a single cell based on the reported longitude and latitude. A quantile regression approach<sup>41,61,62</sup>, allows us to avoid over- and underestimations inherent in using a mean model fit to the global land surface<sup>16</sup> (Extended Data Fig.6). For each unique slope value in the global raster, denudation quantiles were used to construct a cumulative distribution function which could be sampled in each Monte Carlo run (c.f. <sup>41</sup>). We account for differential erodibilities of sedimentary, crystalline metamorphic, and igneous rock types by running regressions between slope values and <sup>10</sup>Be values for each rock type (Extended Data Fig.6). Thus, only <sup>10</sup>Be values from catchments dominated (>80%) by one rock type are used in this regression. This accounting of erodibilities is important as OC-rich shales are weaker and more erodible than OC-poor strong igneous rocks. Residuals between the CRN denudation dataset and the modelled denudation do not change when differential rock erodibility is considered. However, when combined with our OC<sub>petro</sub> stock model, the rock erodibility-corrected OC supply rate model results in 20% higher rates. We also consider the grid-scale bias considered by previous workers<sup>16,41</sup>: as DEM resolution decreases, slope -
as the spatial derivative of elevation - decreases, resulting in an artificial flattening effect¹⁶.

As our Monte Carlo framework is computationally intensive, using a 90-m resolution global raster input would not be feasible. However, we use a 90-m resolution slope dataset to run regression curves shown in Extended Data Fig.6, after which we output a 90-m resolution raster dataset of estimated denudation rates using the median regression curve. By resampling the raster dataset of estimated denudation rates to 1-km resolution after conversion from slope values, we avoid the bias that can lead to an underestimation of denudation by the flattening effect. In our Monte Carlo framework, the quantile regression curves for each raster value can then be sampled to draw a representative denudation value out of the empirical distribution of denudation rates.

## 2.4 Model calibration

The global model is calibrated by minimising the residual with the Re-proxy based estimates of $OC_{petro}$ oxidation (tC yr⁻¹) from 59 globally distributed river basins (Supplementary Table S1). Model selection was performed by running a Monte Carlo simulation (10,000 runs), using the variable $OC_{petro}$ stock and denudation models described above, to find the output which minimises total residuals across all 59 calibration basins simultaneously, such that the sum of all basin residuals was < 1%. These simulations were run on the University of Oxford's Advanced Research Computing (ARC) facility, taking about 24 core hours per simulation. The residuals of individual basins can be quite large for the biggest catchments (e.g., the Amazon basin), but the relative residual especially for the larger basins falls within the uncertainty of model outputs, while accurately predicting the total $OC_{petro}$ oxidation flux globally (Extended Data Fig.7). We note that overall, in basins with moderate $OC_{petro}$ oxidation fluxes, the model may return conservative estimates. However, because this model has the advantage of being globally and spatially explicit, regional over- and underestimation of rock OC oxidation, especially at a local scale (< 10,000 km²), trade off, while we are able to capture regional differences due to tectonics and geological setting (Fig.1c).

## Section 3. Limitations and uncertainties

There is a temporal mismatch between CRN denudation data which informs our probabilistic denudation model, and our Re-proxy calibration data. The Re proxy based $OC_{petro}$ oxidation fluxes used to calibrate our spatial extrapolation model capture fluxes from global rivers within the past decade or less. The CRN technique integrates denudation fluxes that span a millennium or more. Given that anthropogenic land-use change has doubled erosion and weathering since the early 1900s⁶³, our global scale estimates of $OC_{petro}$ oxidation rates therefore reflect the combined influence of natural and anthropogenic activities on global weathering rates which cannot be deconvolved in this present study.

Results of model vs Re-predicted $OC_{petro}$ oxidation fluxes help us to assess the potential for anthropogenic Re input to impact our estimates (Extended Data Fig.7). We have considered anthropogenic Re inputs by removing three large river basins from a previous compilation⁹ and by adding carefully selected river catchment sites to our Re dataset (Section 1.2). In addition, our conversion of Re fluxes to $OC_{petro}$ oxidation is conservative because we uniformly sample the range of Re/OC ratios starting at the lowest measured Re/OC ratio (Section 1.2.2). This leads to error bars within our estimates that are conservatively large.

Most notably, the model outputs of $OC_{petro}$ oxidation vs the Re-estimated fluxes for each basin (Extended Data Figure 7) show a tendency for the model to underpredict smaller catchments more than larger catchments. Our confidence in the weathering signal from Re in the small upland catchments is highest, and the upland, high erosion rate regions that these catchments sample contribute a dominant proportion of the global $OC_{petro}$ flux in our model. Consequently, whilst we cannot completely deconvolve the effect of anthropogenic Re in our
696

constraints, we have confidence that the effect is unlikely to result in a significant overprediction of global OC$_{petro}$ flux estimates.

697

698

The global extrapolation of OC<sub>petro</sub> oxidation proxy data attempts to account for the dataset's underlying heterogeneities in denudation and OC<sub>petro</sub> stocks. However, it does not consider variability in temperature or precipitation which may control weathering – as seen in small scale field measurements of OC<sub>petro</sub> oxidation at sites of high denudation<sup>30,31</sup>. This is primarily due to the size of the Re-proxy catchment database, its spatial coverage and uncertainties inherent in any proxy approach. While the Re-proxy dataset is latitudinally variable (Fig.1a) the model misfit minimization procedure shows the first order controls on flux by OC<sub>petro</sub> stock and denudation (Extended Data Fig.7) meaning that any climatic controls on weathering could be not resolved at the global scale. We note that any bias introduced by extrapolating the global Re-proxy data without including climatic spatial controls on weathering intensity is likely to be minimal, because the underlying dataset spans from the tropics to Arctic locations.

699

700

701

702

703

704

705

706

707

708

709

710

Previous work has suggested that OC<sub>petro</sub> oxidation and cycling of other OC pools can take place during floodplain transport in large fluvial systems<sup>64,65</sup>. While the Re-proxy dataset includes large basins with extensive floodplain areas (Fig.1a), our model-data misfit approach may attribute the downstream fluxes to any higher denudation parts of the catchment. When the model is then upscaled, in lowland floodplain areas where fluvial processing and recycling of sediment<sup>65</sup> can allow biogeochemical reactions to continue<sup>65</sup>, we may predict conservative OC<sub>petro</sub> oxidation rates. At present, we lack empirical data to partition weathering between mountain and floodplain sections<sup>3</sup>. An alternative way to view this is that the removal of alluvial domains from contributing to denudation (since these are depositional) holds minimal control on the overall estimate (&lt; 1%). This result comes from a comparison of model outputs under two parameterization schemes: one where denudation occurs over all ice-free lands vs one where denudation only occurs in ice-free non-alluvial landscapes. Overall, the model's largest contributor to uncertainty is in the conversion of dissolved Re fluxes to OC<sub>petro</sub> oxidation estimates, which are extrapolated in our spatial model. This conversion depends on [OC<sub>petro</sub>]/[Re] ratios which introduce most of the uncertainty in the resulting OC<sub>petro</sub> oxidation rates (Section 1.2.3), and therefore in the model's global output. More constraints on the relative "grey shale" vs "black shale" contribution to catchment Re fluxes could help tighten uncertainty in [OC<sub>petro</sub>]/[Re] ratios (see <sup>66</sup> for a discussion).

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

Finally, our model includes implicit assumptions and features of the datasets which must be acknowledged. First, the model assumes a steady state, which might not accurately describe OC<sub>petro</sub> oxidation in regions responding to changes in uplift, deglaciation, or human activities, which may yet to have reached steady-state conditions.

730

731

732

733

Second, the catchment-scale CRN denudation data used in our model generally derive from lithologies that are quartz-rich and coarse-grained. These typically have lower erodibilities, potentially leading to an underestimation of denudation rates of softer shales which contain the majority of OC<sub>petro</sub> stocks.

734

735

736

737

738

## Methods references:

739

51. Dellinger, M., Hilton, R. G., Torres, M. A., Burt, E. I., Baronas, J. J., Clark, K. E., Galy, V., Ccahuana Quispe, A. J., West, A. J. High rates of rock organic carbon oxidation sustained as Andean sediment transits the Amazon foreland-floodplain. In review in *Proceedings of the National Academy of Sciences*. Preprint at https://eartharxiv.org/repository/list/ (2023).

740

741

742

743
744

745

746

747

52. Ogrič, M., Dellinger, M., Grant, K.E., Galy, V., Gu, X., Brantley, S.L., Hilton, R.G. Low rates of rock organic carbon oxidation and anthropogenic cycling of rhenium in a slowly denuding landscape. *Earth Surface Processes and Landforms* **48**, 1202-1218 (2023). https://doi.org/10.1002/esp.5543

748

749

750

53. Rahaman, W., Singh, S.K. and Shukla, A.D., 2012. Rhenium in Indian rivers: Sources, fluxes, and contribution to oceanic budget. *Geochemistry, Geophysics, Geosystems*, 13(8). https://doi.org/10.1029/2012GC004083

751

752

753

754

54. Sheen, A.I., Kendall, B., Reinhard, C.T., Creaser, R.A., Lyons, T.W., Bekker, A., Poulton, S.W., Anbar, A.D. A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia. *Geochimica et Cosmochimica Acta* **227**, 75-95 (2018). https://doi.org/10.1016/j.gca.2018.01.036

755

756

757

55. Lehner, B., Verdin, K., Jarvis, A. New global hydrography derived from spaceborne elevation data. *Eos, Transactions*, **89**, 93-94 (2008). Data available at https://www.hydrosheds.org.

758

759

56. ESRI, ArcGIS Pro, Release 2.8. Redlands, CA: Environmental Systems Research Institute, 2021.

760

761

57. ALOS World 3D 30 meter DEM, V3.2, Japan Aerospace Exploration Agency, (2021). Distributed by OpenTopography. https://doi.org/10.5069/G94M92HB.

762

763

764

58. HYDRO1k Elevation Derivative Database, Cent. for Earth Resour. Obs. and Sci., Sioux Falls, S. D.: U.S. Geological Survey, (2000). Available from: https://doi.org/10.5066/F77P8WN0

765

766

767

768

769

59. GRDC, Major River Basins of the World / Global Runoff Data Centre, GRDC, 2nd, rev. ext. ed. Koblenz, Germany: Federal Institute of Hydrology (BfG), (2020). Available from: https://www.bafg.de/GRDC/EN/02_srvcs/22_gslrs/221_MRB/riverbasins_node.html# doc2731742bodyText8

770

771

60. Python Software Foundation, Python Language Reference, version 3.9, (2020). Available from: http://www.python.org

772

773

61. Koenker, R. quantreg: Quantile Regression, R package version 5.94, (2022). Available from: https://CRAN.R-project.org/package=quantreg

774

775

62. R Core Team, R: A Language and Environment for Statistical Computing, Vienna, Austria, (2021). Available from: https://www.R-project.org/

776

777

778

63. Syvitski, J., Ángel, J. R., Saito, Y., Overeem, I., Vörösmarty, C. J., Wang, H., Olago, D. Earth's sediment cycle during the Anthropocene. *Nature Reviews Earth & Environment* **3**, 179-196 (2022). https://doi.org/10.1038/s43017-021-00253-w

779

780

781

782

64. Bouchez, J., Beyssac, O., Galy, V., Gaillardet, J., France-Lanord, C., Maurice, L., Moreira-Turcq, P. Oxidation of petrogenic organic carbon in the Amazon floodplain as a source of atmospheric CO2. *Geology* **38**, 255-258 (2010). https://doi.org/10.1130/G30608.1

783

784

785

786

65. Repasch, M., Scheingross, J.S., Hovius, N., Lupker, M., Wittmann, H., Haghipour, N., Gröcke, D.R., Orfeo, O., Eglinton, T.I., Sachse, D. Fluvial organic carbon cycling regulated by sediment transit time and mineral protection. *Nat. Geosci.* **14**, 842-848 (2021). https://doi.org/10.1038/s41561-021-00845-7
787

788

789

66. Charbonnier, Q., Bouchez, J., Gaillardet, J., Calmels, D., Dellinger, M. The influence of black shale weathering on riverine barium isotopes. *Chemical Geology* **594**, 120741 (2022). https://doi.org/10.1016/j.chemgeo.2022.120741

790

791

792

67. Crawford, I., Layton-Matthews, D., Peter, J.M., Gadd, M.G., Voinot, A., Rogers, N. in *Targeted Geoscience Initiative: 2018 Report of Activities* (ed. Rogers, N.) 139-161 (Geological Survey of Canada, 2019). https://doi.org/10.4095/313647

793

794

795

68. Peucker-Ehrenbrink, B. Land2Sea database of river drainage basin sizes, annual water discharges, and suspended sediment fluxes. *Geochemistry, Geophysics, Geosystems* **10**, (2009). https://doi.org/10.1029/2008GC002356

796

Acknowledgments

797

798

The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work.

799

800

http://dx.doi.org/10.5281/zenodo.22558. In addition, this work is based on digital elevation products and high-performance topographic analysis services provided by the

801

OpenTopography Facility with support from the National Science Foundation under NSF

802

Award Numbers 1948997, 1948994 & 1948857. This study was supported by the European

803

Research Council Starting Grant ROC-CO2 678779 (R.G.H.)

804

Author Information

805

Authors and Affiliations

806

University of Oxford, Oxford, UK

807

Jesse R. Zondervan

808

Robert G. Hilton

809

EDYTEM-CNRS-University Savoie Mont Blanc (USMB), Chambéry, France

810

Matthieu Dellinger

811

Durham University, Durham, UK

812

Fiona Clubb

813

Tobias Roylands

814

Mateja Ogrič

815

Contributions

816

J.R.Z. and R.G.H. conceived the study (conceptualisation). J.R.Z., R.G.H., and F.C.

817

developed the methodologies. J.R.Z. and R.G.H. conducted the investigation. Resources were

818

acquired by R.G.H., M.D., T.R., and M.O. Data visualization was done by J.R.Z. and R.G.H.,

819

while data curation was handled by J.R.Z., R.G.H., and M.D. The funding was acquired by

820

R.G.H. J.R.Z. and R.G.H. wrote the original draft of the paper. The paper was reviewed and

821

edited by J.R.Z., R.G.H., M.D., F.C., T.R., and M.O.

822

Corresponding authors

823

Correspondence to Jesse Zondervan or Robert Hilton

824

Ethics Declarations

825

Competing interests

826

Authors declare that they have no competing interests.

827

Extended Data

828

Extended Data Figures 1 to 8
<figure><img src="image_5.png" /></figure>

Extended Data Figure 1.

**Flow chart of data, subroutines and outcomes and the relevant sections of the Methods section.** Probabilistic geospatial models of OC<sub>petro</sub> stocks and denudation are used to extrapolate a global compilation of river catchment OC<sub>petro</sub> oxidation rates across the Earth's surface.

830

831

832

833

834

835

836

837
<figure><img src="image_6.png" /><figcaption>Extended Data Figure 2.</figcaption></figure>

Rhenium (Re) to organic carbon (OC) concentration ratio (10⁻⁷ g g⁻¹) in published black shales⁵⁴,⁶⁷ alongside those measured in riverbed materials (Supplementary Table S3).

Peel riverbed materials capture the input of OC from black shale lithologies, which in other basins are much less common.

838

839

840

841

842

843

844

845
<figure><img src="image_7.png" /></figure>

Extended Data Figure 3.

Rhenium (Re) to organic carbon (OC) concentration ratio (10⁻⁷ g g⁻¹) of riverbed materials, organized by sample location (Supplementary Table S3). Conservatively, the Monte Carlo analysis of OCpetro oxidation rates from the dissolved rhenium proxy uniformly samples the complete measured range of Re/OC values presented here.
<figure><img src="image_8.png" /></figure>

Extended Data Figure 4.

Organic carbon concentrations (OC in weight %) of lithological types in the USGS Rock Geochemical Database (a) and the same figure in log scale (b). Black lines show median values, boxes show the interquartile range (IQR), whiskers show IQR + 1.5 × IQR, and outliers are shown as open symbols (although still included in the analysis).
<figure><img src="image_9.png" /><figcaption>Extended Data Figure 5.</figcaption></figure>

Total continental area (m²) of each mapped units from the GLIM model¹³, which include mixed units of different grain-size siliciclastic units and mixes with carbonates and siliciclastics.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

Extended Data Figure 6.

Cosmogenic isotope-derived denudation rates (¹⁰Be erosion rate mm/kyr) from the OCTOPUS dataset⁴¹ analysed here. Sites are grouped by dominant lithology and the denudation rate as a function of basin average slope (m/km) returned by quantile regression for (a) igneous, (b) metamorphic, and (c) sedimentary dominated catchments. The quantiles corresponding to each grid cell slope value constitute cumulative density functions (CDFs) which can be sampled in each Monte Carlo simulation that quantifies global oxidative weathering rates (Extended Figure 1)
<figure><img src="image_13.png" /></figure>

Extended Data Figure 7.

Rock organic carbon (OC<sub>petro</sub>) oxidation rate (tC yr<sup>-1</sup>) from the Re-proxy versus the model upscaled output for the global set of catchments used in this study

(Supplementary Table S1). Catchments are coloured by the latitude of the sampling point (0° to 80° N and S). Error bars represent uncertainty in Re-proxy values for OC<sub>petro</sub> oxidation (see section 1 & 1.2.3) and the resultant uncertainty in extrapolation model outputs. The dashed line and grey shading shows fit to the data (95% confidence intervals) and the solid black line shows 1:1 line.

877

878

879

880

881

882

883

884

885

886

887

888
<figure><img src="image_14.png" /></figure>

Extended Data Figure 8.

CO₂ release through OCpetro oxidation (black dots, this study; YOCpetro, black line) is more sensitive to sediment yield than CO₂ sequestration through silicate weathering (Ycsil, grey dotted-dashed line)⁴, whilst terrestrial biospheric OC burial (YCorg, green dashed line)¹⁹ could be even more sensitive. CO₂ sequestration trendlines are from Galy et al. ¹⁹. The black dots show CO₂ release through OCpetro oxidation for all the Re sample catchments which have suspended sediment yield data available from the Land2Sea database⁶⁸ (Supplementary Table S1) to allow comparison with the previous assessment of silicate weathering vs biospheric OC burial¹⁹. The exponent for YOCpetro is 0.38 ± 0.08 with r² = 0.43, reflecting the scatter in data owing to the spatial dependence of OCpetro oxidation on rock OC content and denudation rates (Fig. 2). Therefore, comparisons based on the global trends in this plot, without considering their spatial context (see Fig. 2 for spatially explicit comparisons), should be done with care considering the low certainty.