Rock weathering can counteract river CO2 emissions induced by permafrost thaw

Study area and field sampling

Spanning about 30% of the entire QTP, this 780,000 km² study area is situated between the boundary of the Loess Plateau (Gansu Province, 103° E) and the border of China (Tibet Autonomous Region, 78° E), and includes the headwater of eight major river basins: the Yellow, Yangtze, Lancang–Mekong, Nu–Salween, Derung (Irrawaddy), Yarlung Tsangpo–Brahmaputra, Marja Tsangpo (Ganges) and Indus Rivers (Fig. 1a and Supplementary Table 1). This region spans gradients of climate, topography, vegetation cover, permafrost extent and lithology, with elevations between 1,600 and 8,000 m above sea level. The rivers drain the major tectonic terranes of the QTP that contain a wide range of sedimentary, igneous and metamorphic rocks^51,52 (Extended Data Fig. 1b). In the northeast of the study area, the Yellow and Yangtze Rivers drain marine carbonate and siliciclastic sediments of the eastern Kunlun–Qaidam Terrane and the Triassic deep-marine sediments of the Songpan–Ganzi flysch complex⁵¹. The Lancang, Nu and the southwestern Yangtze catchments drain mostly terrestrial siliciclastic and shallow marine carbonate rocks of the eastern Qiangtang and northern Lhasa Terrane⁵¹. In particular, the sediments of the Qiangtang Terrane contain abundant evaporites⁴⁵. By contrast, the catchments of the Yarlung Tsangpo and Indus Rivers are underlain by a series of plutonic and volcanic rocks of the southern Lhasa Terrane apart from the marine Tethyan sedimentary series of the Himalayas (Extended Data Fig. 1b). Permafrost extent across the study area varies from continuous to isolated with the most extensive permafrost extent concentrated in the Yangtze and eastern Yellow Rivers, and high-elevation areas of all other river basins⁵³. The climate is characterized by a cold and dry winter season and an ice-free season between April and October with heavy monsoon rains during the study period⁵⁴. Potential glacial influence is minimized by sampling streams at least 20 km downstream from glaciers. Flow distance and tributaries collectively reduce potential biogeochemical signatures of glacial meltwaters.

We collected 175 individual samples from 50 rivers across this area during the daytime in spring (May–June), summer (July–August) and autumn (September–October) between 2016 and 2018, and in spring ice-out (early April) in 2023. The Yellow River basin was visited eight times; the Yangtze River basin five times; and the Lancang, Nu and Yarlung River basins were all sampled on four dates, the Derung and Indus River basins twice, and the Marja River basin once. We measured, compiled, interpolated and/or calculated all variables first for each of the 175 samples (Supplementary Table 2) and then calculated mean values and their standard deviations for sites with multiple samples to obtain one estimate for each of the 50 sampling sites (Supplementary Table 5).

Topography, hydrology, lithology and climate parameters

All topographic analyses were performed on a 90-m resolution digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM)⁵⁵ using version pre2.5 of the TopoToolbox⁵⁶. We extracted a river network and estimated, for each sample location, drainage-basin-averaged slopes, S_bsn, and normalized steepness index, k_sn, with a reference concavity of 0.45, a segment length of 1,000 m and a drainage-area threshold of 100 km². For all rivers in the network, we also extracted the stream order (with drainage-area threshold of 5 km²) and calculated the total length of rivers, \({L}_{{\mathrm{riv}}_{i}}\) with a given order, that is, upstream of each sampling point. All these analyses were performed on a DEM projected in ArcGIS-Pro using the inbuilt Asia North Lambert Conformal Conic projection. To provide an estimate of the total surface area of rivers upstream of a sampling point, A_riv, we multiplied the river length in each stream order, \({L}_{{\mathrm{riv}}_{i}}\), with an average width for rivers in that stream order, \({\bar{W}}_{i}\):

Where N is the highest stream order in each catchment. For stream orders 3–8, \({\bar{W}}_{i}\) and its uncertainty were based on the mean and standard deviation of 100 individual width measurements (only 30 measurements for stream order 8). Measurements from hydrologic stations at our sampling sites were combined with widths of randomly selected rivers in the study area mapped on Google Earth (Supplementary Table 6). Average width for the first- and second-order streams were estimated by extrapolating an exponential fit to rivers of orders 3–6, with an uncertainty based on the prediction band at these stream orders (Supplementary Fig. 1). We used only stream orders 3–6 for the extrapolation, because many higher-order rivers in the study area are inset into narrow bedrock gorges and break the exponential increase in stream-widths¹⁸ (Supplementary Fig. 1). Our interpolation predicted the first-order stream widths of 4.3 (1.3–13.8) m, which covers the range reported in ref. ⁵⁷ (1.9 ± 1.1 m), but is wider than the value reported more recently in ref. ⁵⁸ (0.32 ± 0.07 m), leading to a potential overestimate in CO₂ emissions.

The same DEM was analysed in TopoToolbox after projecting it with the Asia North Albers Equal Area Conic projection in ArcGIS-Pro to extract drainage areas, A_bsn, and to estimate an areal fraction of permafrost extent in the landscape upstream of each sampling point. We based the latter analysis on a published map of permafrost extent⁵³ that classifies surfaces as being underlain by continuous (≥90%), discontinuous (≥50% and <90%), sporadic (≥10% and <90%) or isolated (<10%) permafrost. Because this classification is based on ranges, we assigned values of 95%, 70%, 30% and 5%, respectively, to each permafrost class (0% for no permafrost), and then calculated the average pixel value upstream of each sampling point in TopoToolbox. Finally, we re-classified all catchment-averaged permafrost-extent values into the same permafrost categories. We generated a TIFF image from the global lithologic map⁵⁹ to extract the fraction of carbonate lithologies upstream of each sampling point in TopoToolbox. In the fraction of carbonate rocks, we included the lithologic categories ‘sc—carbonate sedimentary rocks’ and ‘sm—mixed sedimentary rocks’. We did not consider ‘mt—metamorphics’ here, but we found that including these would not make a substantial difference to our analysis. We estimated the fraction of wetlands, the mean annual air temperature and the mean annual precipitation values for the watersheds draining each of our sites by matching the coordinates of each site with the corresponding watershed in HydroATLAS⁶⁰, in which we obtained the relevant attributes.

For sampling sites within the cross-section of the hydrometric station, we obtained a measurement of discharge, Q_spl (m³ s⁻¹), and suspended sediment concentrations (Supplementary Table 2) at the time of sampling. We were also able to obtain average annual discharge, Q_annual (m³ yr⁻¹), for the most downstream sampling sites of all sampled catchments (basin IDs 9, 16, 18, 19, 22, 27, 31, 37, 45, 46, 47, 48, 49 and 50) and annual suspended sediment for the four largest outlet basins (basin IDs 7, 16, 22 and 27; 74% of the total study area). Note that suspended sediment concentrations were not available for the sampling site with ID9, but for the more upstream sampling site ID7.

Sample collection and analyses

CO₂ concentrations and fluxes

Triplicate river-water samples for \({p}_{{\mathrm{CO}}_{2}}\) were measured in situ from different locations at each site using a Qubit Dissolved CO₂ System equipped with an internal pump and a \({p}_{{\mathrm{CO}}_{2}}\) probe (S157-P, Qubit Biology). The \({p}_{{\mathrm{CO}}_{2}}\) probe was coated with a polytetrafluoroethylene membrane sleeve that is permeable to CO₂ gas, but not to water. The system was calibrated with standard CO₂ gas ranging from 0 μatm to 10,000 μatm. Recalibration of the system was done before each round of fieldwork. Water–air CO₂ emission rates were measured with \({p}_{{\mathrm{CO}}_{2}}\) collection simultaneously using four floating chambers deployed for 1 h. At each transect, chambers were held in place above a range of water depths from shallow water within 2 m of the river bank to deep water in the mid-channel. After mixing the contents in the chamber three times, 50 ml of gas was retrieved from the chambers after 0 min, 2 min, 5 min, 10 min, 30 min and 60 min of measurement and transferred to air-tight gas sampling bags for concentration determination in the laboratory by gas chromatograph equipped with a flame ionization detector (Agilent 7890B GC-FID). Characteristics and construction of the chambers followed previously tested setups that showed minimally biased results⁶¹. Local ambient air samples were also taken and used to back-calculate CO₂ concentrations in water in equilibrium with the atmosphere.

Riverine CO₂ emission rates, \({F}_{{{\rm{CO}}}_{2}}\) (mmol m⁻² day⁻¹) were derived as follows:

where k is the gas transfer velocity (m day⁻¹), C_water is the water gas concentration (mol m⁻³) and C_eq is the gas concentration (mol m⁻³) in water in equilibrium with the local atmosphere corrected for temperature-induced changes in solubility according to Henry’s law. We assumed that \({F}_{{{\rm{CO}}}_{2}}\) is solely diffusive (that is, CO₂ ebullition is negligible). We obtained k from a linear regression of \({p}_{{\mathrm{CO}}_{2}}\) against time to eliminate possible bias due to gas accumulation in the chamber headspace. Fluxes from all chambers were averaged to yield one measurement per sample.

Physicochemical analyses

Water samples were collected simultaneously with gas samples for physicochemical analyses at each site. Water temperature, pH and dissolved oxygen were measured in situ with portable field probes (Hach HQ40d). Air temperature, air pressure and wind speed were measured in situ with a portable anemometer (Testo 480).

Water samples were filtered through 0.22 µm polyethersulfone syringe filters to measure cation and anion concentrations. We rinsed each polyethylene bottle with filtered water before filling, and acidified each bottle for cation analysis with ultrapure HNO₃. We measured major cations (K⁺, Na⁺, Mg²⁺ and Ca²⁺) and anions (Cl⁻ and SO₄²⁻) on an ion chromatograph (Dionex ICS-1100). NH₄⁺ and NO₃⁻ were determined colorimetrically with an Autoanalyzer-3 (Bran+Luebbe) following the methods in ref. ⁶². We estimated HCO₃⁻ by charge balance. Dissolved silica (Si) was analysed on a Varian 720 ICP-OES. DOC concentrations were measured on duplicate filtered water samples acidified with H₂SO₄ following ref. ¹⁸. DIC concentrations were calculated from the sum of HCO₃⁻ and \({p}_{{\mathrm{CO}}_{2}}\) concentrations. To determine the POC content of the suspended sediments, water samples were filtered in duplicates through pre-weighed combusted 0.45 µm filters. The filters were then dried at 105 °C for at least 12 h. Before POC analysis, the samples were acidified with 2 M HCl to remove inorganic carbon. No rinsing step was applied after acidification to avoid washing fine particulate material or acid-labile organic matter off the filters. The filters were subsequently dried again at 105 °C for over 5 days to ensure complete removal of residual acid and moisture. Afterwards, the filters were reweighed, and POC mass on the filters was quantified using a Vario Micro Cube elemental analyser (Elementar) with the furnace temperature set at 980 °C.

Sulfur and oxygen isotope analyses

Stable isotopic composition of sulfur, oxygen and carbon is reported in the text using delta notations:

where R_sample and R_standard are ratios of heavy to light isotopes (³⁴S:³²S, ¹⁸O:¹⁶O or ¹³C:¹²C) in samples (R_sample) and in standards (R_standard) of standard mean ocean water or Vienna PeeDee Belemnite, respectively.

Stable isotopes of δ³⁴S-SO₄²⁻ and δ¹⁸O-SO₄²⁻ of dissolved SO₄²⁻ were measured following ref. ⁶³ for 34 river samples collected during the fall 2018 sampling campaign. In brief, sulfate was eluted off a Dowex resin column with 20 ml 0.8 M HCl and mixed with BaCl₂ to precipitate barite (BaSO₄). The precipitate was subsequently cleaned with 6 M HCl and rinsed three times with ultrapure water. After removing impurities with 10 ml 0.05 M DTPA, the barite was then cleaned three times with ultrapure water and dried at 70 °C. For δ³⁴S-SO₄²⁻ analysis, 400 μg of barite was combusted in excess oxygen with vanadium pentoxide in a Flash EA coupled by continuous flow and analysed by an isotope ratio mass spectrometry (Thermo Scientific 253 Plus). For δ¹⁸O-SO₄²⁻ analysis, 180 μg of barite was pyrolysed in a thermal conversion element analyser and passed through continuous He flow into an isotope ratio mass spectrometer (Thermo Scientific 253 Plus) using a ConFlo 3. We determined δ¹⁸O-H₂O using a wave scanning cavity ring-down spectrometer (Picarro L2130i) as well. Additional triplet δ³⁴S-SO₄²⁻, δ¹⁸O-SO₄²⁻ and δ¹⁸O-H₂O data for rivers within the Yellow, Yangtze, Lancang, Nu, Yarlung Tsangpo and Indus catchments were obtained from the literature (Supplementary Table 7). We had measurements from only 34 of the 50 sampling sites. For the remaining 16 sites, we substituted data from the same or closest sites reported in the literature. We then calculated the average for sites with multiple measurements to obtain a single estimate for each of the 50 sampling sites.

Carbon isotope analysis

We collected samples for δ¹³C–Δ¹⁴C analyses of DOC, DIC and CO₂ both near the river bank and in the mid-channel from the sampled rivers from 2017 to 2023 (except 2019 because of COVID restrictions). For Δ¹⁴C and δ¹³C measurements of DOC and DIC, the river-water samples were filtered using 0.45 μm polysulfone filters, then the near-bank and mid-channel samples were mixed and collected in acid-washed (10% HCl v/v, 24 h) 1 l HDPE Nalgene bottles that were rinsed three times with filtered river water before collection. Four submerged bubble traps were deployed from the near bank to the mid-channel to collect gas samples for Δ¹⁴C-CO₂ and δ¹³C-CO₂ analyses. The traps remained in place for 2–4 h before sampling. Accumulated gas was transferred into 10-l air-tight Tedlar PVF gas sampling bags to prevent contamination by contemporary atmospheric CO₂. δ¹³C was analysed on an isotope ratio mass spectrometer (Thermo Scientific 253 Plus), and Δ¹⁴C was measured at the accelerator mass spectrometry facility at Peking University, Guangzhou Institute of Geochemistry and Beta Analytic Laboratory.

We also compiled 92 concurrent δ¹³C–Δ¹⁴C measurements each for DOC and DIC within the QTP from the literature. Finally, we were able to obtain 378 total observations of fluvial DOC (n = 163), DIC (n = 147) and CO₂ (n = 68) Δ¹⁴C data, including our own measurements outlined above.

Endmembers of carbon isotopes

For Fig. 2a, we used the following δ¹³C ranges for endmember carbon sources: terrestrial OC derived from globally dominant C3 plant photosynthesis (−32‰ to −24‰)⁶⁴; autochthonous, in-stream OC production (−23‰ to −17‰)¹³; and rock sources, including kerogen (−32‰ to −17‰), carbonate rocks (−2.5‰ to 0‰) and solid Earth degassed CO₂ (−6.5‰ to 10‰ dependent on source rocks)⁶⁴.

The Stable Isotope Mixing Models in R (simmr) package was used to estimate possible contributions of three endmember sources for all Δ¹⁴C samples as detailed described in ref. ⁴⁶. Endmember sources for the QTP were further constrained by using measured Δ¹⁴C values obtained from QTP soils (Supplementary Table 8). A mass balance fitting the potential contributions of the endmember sources was solved by using a Markov Chain Monte Carlo method:

where C14_sample represents the ¹⁴C content of the different waterborne carbon components of each sample (in pmC), which is the sum of the mass balance among the proportional contributions of three endmember sources. To solve this mass balance, the relative contributions of the three endmember sources were assumed to sum to 1. We divided the soil column into top (<50 cm) and deep (≥50 cm) pools to where the topsoil source (C14_topsoil) is a simplified representation of the uppermost soil carbon accumulation in the active layer, reflecting the seasonally thawed, often younger soils. The deep soil source (C14_{deep soil}) represents older carbon from the deep active layer to the permafrost, including the Holocene and late Pleistocene aged sources from decaying over 5,000 years bp carbon, and decaying more than 10,000 years bp carbon, respectively. The rock source (C14_rock) represents Δ¹⁴C-dead carbon sources from kerogen and carbonates.

Unmixing major elements

For all 175 samples, we used the inverse model mixing elements and dissolved isotopes in rivers (MEANDIR)⁶⁵ to unmix contributions of solute sources to our measured riverine chemistry.

Considered solutes and endmembers

In the inversion, we considered measured concentrations of sodium [Na⁺], calcium [Ca²⁺], magnesium [Mg²⁺], chloride [Cl⁻], sulfate [SO₄²⁻] and silicate [Si], as well as the isotopic composition of sulfur in sulfate (δ³⁴S-SO₄²⁻). All elemental concentrations were normalized to the sum of the major ions [Na⁺], [Ca²⁺], [Mg²⁺] and [SO₄²⁻] (refs. ^45,65). We considered contributions from atmospheric sources and from the weathering of evaporites, sulfides, carbonates and silicates. Anthropogenic inputs were considered negligible because the area is far from major industrial activity or human settlements, and we did not consider the inputs from hot springs (Supplementary Text 3). The size of the study area and lithologic heterogeneity is large, so we used endmember compositions of silicate and carbonate weathering that reflect global averages^23,65,66 (Supplementary Text 4 and Supplementary Table 9). The range of precipitation endmembers was derived from the 5th and 95th percentiles of average precipitation compositions in each of the drainage basins (Supplementary Table 10). A range of different evaporite minerals contribute to the solute budget of the QTP rivers, which complicates the distinction between evaporite and non-evaporite solute sources⁴⁵. To address this challenge, we defined a very broad evaporite endmember that allowed for a wide range of compositions (Supplementary Table 9). Moreover, to distinguish evaporite and silicate weathering sources, we included silica in the inversion and defined a range of allowed silica concentrations for the silicate endmember using the 5th and 95th percentiles of a compilation of unmixed silicate weathering concentrations from small mountain streams for which evaporite contributions are negligible³⁵. Finally, we used published δ³⁴S compositions of sulfides and sulfur-bearing evaporites from our study area, to unmix sulfur sources (Supplementary Table 11).

MEANDIR model parameters

MEANDIR uses a Monte Carlo approach to estimate the range of fractional endmember contributions that fit the measured data. At each iteration, the model picks a set of endmember compositions randomly from user-defined endmember ranges and then fits the fractional contributions of endmembers that best fit the observed data. The final result is commonly obtained as the median of a given fraction or number of model runs with the lowest misfit between the model and the observations. Whereas this approach produces a solution to all samples, it also allows very poor fits to the data. Hence, we iterated over each sample until we had 100 solutions that fit all of the ion data within ±15% and the δ³⁴S-SO₄²⁻ data within ±2‰. We allowed a maximum of 100,000 iterations to run before the inversion was stopped. Sixteen out of the 175 samples did not have all of the major element analyses data and 10 samples did not produce the necessary 100 solutions within the defined error bounds during the inversion. Hence, only 149 samples were considered in the endmember unmixing and carbon balance. All MEANDIR input parameters are reported in the Supplementary Code.

Treatment of potassium

Inclusion of potassium [K⁺] in the inversion doubles the number of samples that cannot be well-modelled. Because [K⁺] and [Na⁺] are well-correlated (Supplementary Table 2), we made the simplifying assumption that fractional contributions to [K⁺] scale linearly with the fractional contributions to [Na⁺]. Because the concentrations of potassium are generally low, the uncertainty introduced by this approximation is small.

Areal CO₂ fluxes from weathering

From the cation charge contributed by carbonate and silicate weathering and the anion charges contributed by sulfate from sulfide oxidation, we estimated the balance of CO₂ sinks and sources from weathering in the landscapes drained by each sampled river. Any carbonate that was already precipitated in streams or soils upstream of the sampling point did not contribute to our CO₂ budget and is not considered here. We followed previous approaches^28,67 to estimate the following fluxes of CO₂, \({F}_{{{\rm{CO}}}_{2}}\), in mass of carbon per area per time (tC km⁻² yr⁻¹) with positive numbers being CO₂ production and negative numbers being CO₂ sequestration:

1. 1.

   CO₂ consumed by silicate weathering, \({F}_{{(\text{C}{\text{O}}_{2})}_{\text{sil}}}\):

   $${F}_{{({\mathrm{CO}}_{2})}_{\mathrm{sil}}}=-(1-{f}_{\mathrm{sulf}})\times {M}_{{\rm{C}}}\times {[\mathrm{Cat}]}_{\mathrm{sil}}^{\mathrm{eq}}\times \frac{{Q}_{{\rm{w}}}}{{A}_{\mathrm{bsn}}}$$

   (5)

2. 2.

   CO₂ consumed by carbonate weathering, \({F}_{{(\text{C}{\text{O}}_{2})}_{\text{carb}}}\):

   $${F}_{{({\mathrm{CO}}_{2})}_{\mathrm{carb}}}=-0.5(1-{f}_{\mathrm{sulf}})\times {M}_{{\rm{C}}}\times {[\mathrm{Cat}]}_{\mathrm{carb}}^{\mathrm{eq}}\times \frac{{Q}_{{\rm{w}}}}{{A}_{\mathrm{bsn}}}$$

   (6)

3. 3.

   CO₂ produced by sulfuric acid weathering of carbonate, \({F}_{{(\text{C}{\text{O}}_{2})}_{\text{sulf}}}\):

where f_sulf () is the fraction of weathering by sulfuric acid, \({f}_{\text{sulf}}=\frac{{[{\mathrm{SO}}_{4}^{2-}]}_{\text{sulf}}^{\text{eq}}}{{[\mathrm{Cat}]}_{\text{sil}}^{\text{eq}}+{[\mathrm{Cat}]}_{\text{carb}}^{\text{eq}}}\), \({[\mathrm{Cat}]}_{\text{sil}}^{\text{eq}}\) and \({[\mathrm{Cat}]}_{\text{carb}}^{\text{eq}}\) (charge equivalents in mol m⁻³) are the cation charges contributed by silicate and carbonate weathering, respectively, M_c (g mol⁻¹) is the molar mass of carbon, Q_w (m³ yr⁻¹) is the annual average discharge, and A_bsn (m²) is the drainage basin area. Equations (5)–(7) assume that the sulfuric acid partitions between silicate and carbonate phases in the proportion of the fraction of cation charges contributed by each reaction. Depending on the conditions in the soil, this simplification may not apply⁶⁸. However, changing this assumption would not affect the net CO₂ balance from weathering reactions, that is, the sum of equations (5)–(7), but would instead trade off CO₂ drawdown fluxes between silicates and carbonates. We note that in the calculation of the CO₂ balance, we excluded one (out of 175) sample with f_sulf > 1, which is impossible. This estimate must have arisen because of the high evaporite load that was not correctly estimated by the inversion model. We also excluded one sample that had an anomalously high sulfate concentration compared with all other 175 samples (almost 10 times higher than the 75th percentile of all samples and twice as high as the next highest concentration), and disproportionately dominated the upscaling described in the next sections. Finally, we note that we did not have annual discharge estimates for all sampling sites (see section ‘Topography, hydrology, lithology and climate parameters’), and instead upscaled Q_w (m³ yr⁻¹) from the discharge estimate at the time of sampling, Q_spl (m³ s⁻¹), for each of the 175 samples by simply multiplying the value by the number of seconds in a year. Therefore, Q_spl for the 50 sampling sites is effectively a weighted average flux. Note that these fluxes were used only to compare relative changes across the permafrost gradient. To obtain the total fluxes from the study area, we used actual measured annual averages for the most downstream sites of the catchments (see sections ‘Topography, hydrology, lithology and climate parameters’ and ‘Total carbon fluxes from the study area’).

Areal CO₂ fluxes from river emissions

We estimated areal CO₂ emission fluxes from each of the sampled catchments (in mass of carbon per area per time) based on CO₂ emission flux rates (\({F}_{{\mathrm{CO}}_{2}}\)) measured with the floating chambers as

where \({F}_{{({{\rm{CO}}}_{2})}_{{\rm{rem}}}}\) (mol m⁻² day⁻¹) is the upscaled riverine emission (rem) flux.

Total carbon fluxes from the study area

We estimated the total CO₂ fluxes from weathering and river emissions as well as downstream fluxes of DOC and DIC across the study area in TgC yr⁻¹. To this end, we summed the contributions of CO₂ fluxes from the 14 most downstream sampling sites of each sampled basin. We calculated the summed weathering, DOC, and DIC fluxes based on the annual discharge estimates, Q_annual multiplied by the average concentrations in the basin, and we used the standard deviations of the concentrations to estimate the uncertainty of the summed value (Supplementary Table 5). Here, we were able to use actual measured annual discharge measurements and did not use the Q_w values upscaled from individual discharge measurements at the time of sampling. Note that we used the solute and carbon concentrations of basin 7 instead of the further downstream basin 9 because we did not have any major element chemistry for basin 9. To account for uncertainty in the stream-area calculations, we estimated the summed CO₂ emission fluxes, \({F}_{{(\text{C}{\text{O}}_{2})}_{\text{rem}}}\), with 100,000 Monte Carlo runs. First, we randomly picked values of river widths from a normal distribution for stream orders 3–8 (Supplementary Table 5) and from a log-normal distribution for the interpolated stream orders 1–2. Multiplied by stream lengths (equation (1)), these values yield 100,000 estimates of stream areas upstream of each sampling site. We then randomly picked 100,000 \({F}_{{({{\rm{CO}}}_{2})}_{{\rm{rem}}}}\) values from a normal distribution based on the average and standard deviation for each sample site (Supplementary Table 5) and solved equation (8) for each of the 100,000 Monte Carlo runs. The final reported value is the median and standard deviation from all 100,000 runs. We also multiplied annual suspended sediment fluxes by POC content of the largest four of the outlet basins to obtain POC fluxes. Because our measurements were limited to the ice-free season between April and October, and because our first- and second-order stream widths are higher than expected, we probably overestimate CO₂ emissions over the year.

Statistical analysis

Using R (v.4.4.2), we performed a PCA to explore the main environmental drivers and their effects on our variables of interest. The drivers considered were catchment size, average k_sn, permafrost cover, average annual precipitation, average annual air temperature, wetland cover and fraction of carbonate rocks within the catchment. These variables were scaled and centred before the PCA. We then used the two main axes of the PCA (PC1 and PC2) to explore how they were related to several variables using linear regression models. Apart from the PCA, we used an MLR model, in which each response variable was regressed against all key environmental variables (Extended Data Table 1 and Supplementary Table 4). All predictor and response variables were standardized before analysis to facilitate the comparison of effect sizes using standardized regression coefficients. We used a bidirectional stepwise variable selection procedure based on the Akaike information criterion to identify the most parsimonious set of drivers for each response variable. Model performance and the significance of individual predictors were assessed using adjusted R² and P values.