Strange IndiaStrange India

The workflow of materials and methods (Extended Data Fig. 1) starts with the compilation and derivation of rhenium-concentration-based OCpetro oxidation flux estimates, discussed just below. Then we detail the ‘Global spatial OCpetro oxidation model’ and its application of submodels for OCpetro stocks and denudation, and the Monte Carlo routines used in the model’s approach. Finally, we discuss the ‘Limitations and uncertainties’ involved in our methods and calculations.

Rhenium-based river catchment estimates of OCpetro oxidation

From a series of dissolved rhenium measurements (typically completed by ICP-MS), the dissolved Re flux JRe (t yr−1) can be used to estimate OCpetro oxidation flux, JOCpetro-ox (tC yr−1) using:

$${J}_{{\rm{O}}{\rm{C}}{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}-{\rm{o}}{\rm{x}}}=\,{J}_{{\rm{R}}{\rm{e}}}\times {\left(\frac{[{{\rm{O}}{\rm{C}}}_{{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}}]}{[{\rm{R}}{\rm{e}}]}\right)}_{i}\times {f}_{{\rm{c}}}$$


where fc is the fraction of dissolved Re derived from OCpetro oxidation34 and ([OCpetro]/[Re])i is the organic carbon to rhenium concentration ratio (g g−1) in rocks of type i undergoing weathering. In some catchments where it may be important, an additional term, not shown in equation (2), has been applied to correct for the presence of graphite, which may not undergo alteration during weathering33.

Compiled published measurements

In this study we compile estimates of OCpetro oxidation using the dissolved Re proxy from published literature (Supplementary Table 1). These include the Yamuna River, India24; ten Taiwanese rivers7; four rivers from the western Southern Alps33; four rivers from the Mackenzie Basin, Canada3,34; and six rivers draining the Peruvian Andes51. Two Swiss catchments25 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 records34, while earlier studies all include single snapshot samples7,24, 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, akin to fc (equation (2)), apart from the Taiwan dataset7. While uncertainties on the OCpetro oxidation yields appear relatively large (Supplementary Table 1), it is important to note that the measured range in yields is much larger than the uncertainties.

New estimates of OCpetro oxidation

To expand the 25 estimates of OCpetro 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 world9. We use these measurements and combine them with estimates of fc and ([OCpetro]/[Re])i, discussed in the following sections, to calculate OCpetro oxidation yields with associated uncertainties using published approaches25,34.

In locations with substantial local sources of fossil fuel combustion (for example, coal-fired power plants or steel works), rainwater can contain concentrations of Re that approach those of river water8,52, whereas locations that have minimal impacts from local pollution sources have Re concentrations in rainwater that are below detection25,33. In the large river dataset9, some large rivers are noted to have markedly increased Re concentrations and fluxes; the conclusion is 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 balance52. A study of Re across Indian catchments suggests that while Re in Himalayan catchments and the mainstem Ganges and Brahmaputra behave conservatively, peninsular lower relief catchments with denser populations and industrial activity suggest anthropogenic inputs53. For this study’s purpose, 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 1). We further consider the role of anthropogenic Re in our model results in the ‘Limitations and uncertainties’ section.

Estimation of Re source and f

To estimate the fraction fc of dissolved Re sourced from OCpetro for the rivers in the Re flux dataset9, we follow a previously used forward model mixing approach25,34:

$${{\rm{Re}}}_{{\rm{OC}}}={{\rm{Re}}}_{{\rm{tot}}}\times {f}_{{\rm{c}}}={{\rm{Re}}}_{{\rm{tot}}}-{{\rm{Re}}}_{{\rm{sulf}}}-{{\rm{Re}}}_{{\rm{sil}}}$$


where ReOC is the rhenium concentration of OCpetro-derived Re in the dissolved load, Retot is the measured Re concentration, Resulf and Resil are the concentrations derived from weathering of sulfide and silicate minerals, respectively. These unknowns can be quantified as:

$${{\rm{Re}}}_{{\rm{sulf}}}=\left[{{\rm{SO}}}_{4}\right]\times {\left(\frac{\left[{\rm{Re}}\right]}{\left[{{\rm{SO}}}_{4}\right]}\right)}_{{\rm{sulf}}}$$


$${{\rm{Re}}}_{{\rm{sil}}}=\left[{\rm{Na}}\right]\times {\left(\frac{\left[{\rm{Re}}\right]}{\left[{\rm{Na}}\right]}\right)}_{{\rm{sil}}}$$


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

Constraints on ([OCpetro]/[Re])i

A recent compilation54 provides measurements of ([OCpetro]/[Re])i from rock samples of different ages around the world. However, most of these measurements were made on black shales with OCpetro contents greater than 1%, which occur on only 0.3% of the Earth surface (Extended Data Fig. 5). 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 [OCpetro]/[Re] on bed materials from rivers around the world (Supplementary Table 3 and 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 previously25. We find that regions with lower OCpetro concentrations that are more typical of sedimentary rocks at the continental surface—units including fine-grained sedimentary rocks that make up more than 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 basin34 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 dataset9 is 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 OC to rhenium ratio in the rocks. To conservatively quantify uncertainty in the range of OCpetro oxidation rates from dissolved-Re data, we perform a Monte Carlo simulation in which we uniformly sample the entire range of measured [OCpetro]/[Re] values, from low values indicative of carbon-poor and/or metamorphic rocks 2.5 × 10−8 g g−1 (ref. 33) towards 1.26 × 10−6 g g−1 (ref. 34) in catchments with higher OC in rocks (Supplementary Table 3 and Extended Data Fig. 3).

OCpetro oxidation yields and uncertainties

Equation (2) is used for each basin in the Re flux dataset9. Uncertainties in fc derive from the range of values used in the sulfide and silicate end member compositions (equations (4) and (5)). For [OCpetro]/[Re], we use the range of values discussed in the just-previous section on constraints (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 JOCpetro-ox for each basin. The median value of the Monte Carlo simulation and the interquartile range are reported (Supplementary Table 1).

Geospatial catchment boundaries

To derive the catchment outlines and areas corresponding to the Re-proxy samples in our compiled dataset, we used the HYDROSHEDS flow direction grid at 3 arc-second resolution55 and ArcGIS Pro56. Catchments outside the latitudinal cover of HYDROSHEDS were derived from the HYDRO1K flow direction grid57 and catchments in Iceland were derived from ALOS AW3D using TauDEM functionality in OpenTopography58. While most published sample coordinates (Supplementary Table 1) give the correct location on the cited drainage systems, in a handful of cases, coordinates had to be amended by up to a few kilometres, which may reflect errors in transcribing (for example, Kikori and Purari9). Final quality control included a comparison of the extracted drainage basin areas and those published, with good agreement overall (less than 2% residual). However, some drainage areas cited in the Re flux dataset9 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. Finalized coordinates of Re samples determined for each drainage system, with the corresponding upstream drainage area, are given in Supplementary Table 1. Spatial files of upstream drainage boundaries and Re sample locations are available on Zenodo (available from To convert dissolved Re concentrations into Re fluxes, average annual water discharge was calculated using published numbers at gauges (Supplementary Table 1) 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 weathering10 that use the GRDC dataset in Major River Basins of the World59. Drainage areas used by ref. 10 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 4).

Global spatial OCpetro oxidation model

In the following three sections, we provide additional rationale and details of the modelling approaches. The model procedures apply two spatial probabilistic subroutines; one deals with OCpetro stocks in surface rocks and the other with 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 OCpetro oxidation rates (Extended Data Fig. 1). Model simulations were implemented at 1-km grid scale (Mollweide projection, WGS84 datum) in the Python programming language60.

OCpetro stocks

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

$${[{{\rm{O}}{\rm{C}}}_{{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}}]}_{{\rm{r}}{\rm{o}}{\rm{c}}{\rm{k}}}={F}_{1}({[{{\rm{O}}{\rm{C}}}_{{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}}]}_{{\rm{l}}{\rm{i}}{\rm{t}}{\rm{h}}{\rm{o}}{\rm{l}}{\rm{o}}{\rm{g}}{\rm{y}},1})+{F}_{2}({[{{\rm{O}}{\rm{C}}}_{{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}}]}_{{\rm{l}}{\rm{i}}{\rm{t}}{\rm{h}}{\rm{o}}{\rm{l}}{\rm{o}}{\rm{g}}{\rm{y}},2})+\,\cdots \,+{F}_{n}({[{{\rm{O}}{\rm{C}}}_{{\rm{p}}{\rm{e}}{\rm{t}}{\rm{r}}{\rm{o}}}]}_{{\rm{l}}{\rm{i}}{\rm{t}}{\rm{h}}{\rm{o}}{\rm{l}}{\rm{o}}{\rm{g}}{\rm{y}},n})$$


$${F}_{1}+{F}_{2}+{\rm{\cdots }}+{F}_{n}=1$$


Denudation model

The denudation model is parametrized using a regression approach, similar to techniques applied elsewhere16,41. We regressed a compilation of long-term catchment-scale 10Be denudation estimates39 against mean local slope generated from the Geomorpho90m dataset40. Mean local slope was calculated using the focal statistics tool in ArcGIS Pro56 and the Geomorpho90m slope dataset with a 5-km moving radius. Slope values were then matched to 10Be denudation estimates at a single cell based on the reported longitude and latitude. A quantile regression approach41,61,62, allows us to mitigate over- and underestimations inherent in using a mean model fit to the global land surface16 (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 (compare ref. 41).

We account for differential erodibilities of sedimentary, crystalline metamorphic and igneous rock types by running regressions between slope values and 10Be values for each rock type (Extended Data Fig. 6). Thus, only 10Be values from catchments dominated (more than 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 OCpetro 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 workers16,41: as DEM resolution decreases, slope—as the spatial derivative of elevation—decreases, resulting in an artificial flattening effect16. 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 as 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.

Model calibration

The global model is calibrated by minimizing the residual with the Re-proxy-based estimates of OCpetro oxidation (tC yr−1) from 59 globally distributed river basins (Supplementary Table 1). Model selection was performed by running a Monte Carlo simulation (10,000 runs), using the variable OCpetro stock and denudation models described above, to find the output which minimizes total residuals across all 59 calibration basins simultaneously, such that the sum of all basin residuals was less than 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 (for example, the Amazon basin), but the relative residual, especially for the larger basins, falls within the uncertainty of model outputs, while accurately predicting the total OCpetro oxidation flux globally (Extended Data Fig. 7). We note that, overall, in basins with moderate OCpetro 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 OCpetro oxidation found mostly at a local scale (less than 10,000 km2) tend to trade off while we are able to capture larger regional differences due to tectonics and geological setting (Fig. 1d).

Limitations and uncertainties

There is a temporal mismatch between the CRN denudation data that inform our probabilistic denudation model, and our Re-proxy calibration data. The Re-proxy-based OCpetro 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. Anthropogenic land-use change has doubled erosion and weathering since the early 1900s (ref. 63); hence, our global scale estimates of OCpetro oxidation rates reflect the combined influence of natural and anthropogenic activities on global weathering rates, which cannot be deconvolved in this present study.

Results of model versus Re-predicted OCpetro oxidation fluxes help us 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 compilation9 and by adding carefully selected river catchment sites to our Re dataset (see the Methods section headed ‘Rhenium-based river catchment estimates of OCpetro oxidation’). In addition, our conversion of Re fluxes to OCpetro oxidation is conservative because we uniformly sample the range of Re/OC ratios starting at the lowest measured Re/OC ratio (see the Methods section headed ‘Constraints on ([OCpetro]/[Re])’). This leads to error bars within our estimates that are conservatively large. Most notably, the model outputs of OCpetro oxidation versus the Re-estimated fluxes for each basin (Extended Data Fig. 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 OCpetro flux in our model. While we cannot completely deconvolve the effect of anthropogenic Re in our constraints, we have confidence that the effect is unlikely to result in a significant overprediction of global OCpetro flux estimates.

The global extrapolation of OCpetro oxidation proxy data attempts to account for the dataset’s underlying heterogeneities in denudation and OCpetro stocks. However, it does not consider variability in temperature or precipitation which may control weathering—as seen in small-scale field measurements of OCpetro oxidation at sites of high denudation30,31. 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 OCpetro 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.

Previous work has suggested that OCpetro oxidation and cycling of other OC pools can take place during floodplain transport in large fluvial systems64,65. 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 sediment65 can allow biogeochemical reactions to continue65, we may predict conservative OCpetro oxidation rates. At present, we lack empirical data to partition weathering between mountain and floodplain sections3. 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 (less than 1%). This result comes from a comparison of model outputs under two parametrization schemes: one where denudation occurs over all ice-free lands versus 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 OCpetro oxidation estimates, which are extrapolated in our spatial model. This conversion depends on [OCpetro]/[Re] ratios which introduce most of the uncertainty in the resulting OCpetro oxidation rates (see the Methods section ‘OCpetro oxidation yields and uncertainties’) and therefore in the model’s global output. More constraints on the relative ‘grey shale’ versus ‘black shale’ contribution to catchment Re fluxes could help tighten uncertainty in [OCpetro]/[Re] ratios (see ref. 66 for a discussion).

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 OCpetro oxidation in regions responding to changes in uplift, deglaciation or human activities, which may not yet have reached steady-state conditions.

Second, most catchment-scale CRN denudation data used in our model 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 OCpetro stocks.

Source link


Leave a Reply

Your email address will not be published. Required fields are marked *