Introduction

Carbon dioxide (CO2) emissions from rivers and streams account for about 85% of global inland water CO2 emissions1, making it a crucial component of the global carbon cycle. Permafrost covers up to 21 million km2 of the Northern Hemisphere (NH) land area2. Due to the cold and wet environmental conditions, these permafrost areas store large amounts of soil organic carbon, accounting for approximately 50% of global carbon storage3. Climate warming accelerates the widespread thawing of permafrost, causing ground temperature increase, ground ice melting, and permafrost loss4,5. Dramatic permafrost degradation has increased the lateral carbon exports to adjacent freshwater ecosystems6, where it is partly mineralized and evaded as greenhouse gases into the atmosphere, further enhancing the strength of permafrost carbon feedback to climate change. For now, although terrestrial ecosystem in the NH permafrost regions serves as a carbon sink, the CO2 emissions from inland waters largely offset the sink7. However, terrestrial carbon budgets rarely incorporate rivers8,9, and temporally and spatially explicit estimate of riverine CO2 emissions and their impacts on land carbon sink in the NH permafrost regions is still missing.

Several preceding studies have been dedicated to assessing the magnitude of CO2 emissions from rivers and streams on regional, national and global scales10,11,12,13. Meanwhile, riverine CO2 emissions in permafrost regions11,14,15,16,17 are gaining recognitions due to their important contributions to global carbon emissions. Based on the global river carbon flux dataset, the riverine CO2 emissions in the pan-Arctic permafrost regions were estimated at 107–223 Tg CO2–C yr−1 7. Nevertheless, for the whole NH permafrost regions, current estimate of riverine CO2 release focuses on the multiyear-average emission assessment18, but the annual spatio-temporal changes of riverine CO2 emissions remain unconstrained. This greatly hinders the understanding of the role of riverine CO2 emission in offsetting the terrestrial carbon sink under future climate scenarios.

It is likely that impacts of climate warming on the carbon cycle will vary across the permafrost zone. Notably, the extensive permafrost regions in both the Arctic and Tibetan Plateau (the largest mountain permafrost area globally) experience amplified warming compared to the global average19,20. Compared to the Arctic, permafrost in the Tibetan Plateau is characterized by higher ground temperatures, deeper active layer thicknesses and lower organic carbon storage21,22. Meanwhile, permafrost degradation in the Tibetan Plateau has been more severe and extensive than in the Arctic5,23. However, previous studies have often treated these permafrost regions separately, overlooking the opportunity to compare their riverine CO2 emissions in a unified framework. Thus, the driving mechanisms for riverine CO2 emissions between the Arctic and Tibetan Plateau permafrost regions are poorly understood, implying uncertainty in understanding and projecting the future coupled land-water carbon cycle.

Riverine CO2 emissions mainly depend on the supply of carbon from soil and shallow subsurface input6. This process is influenced by factors such as climate24, permafrost14 and land use25. However, for the NH permafrost zones, knowledge of the influence of these factors on the distribution and variations of riverine CO2 emissions remain limited. Furthermore, although aquatic CO2 emissions may offset a considerable portion of the terrestrial carbon sink, terrestrial carbon budgets for permafrost regions in Earth System Models (ESMs) typically do not account for the contribution of riverine carbon emissions26. These limitations hinder accurate estimations of riverine CO2 emissions from permafrost regions and future projections of permafrost carbon–climate feedback. Therefore, given the accelerating permafrost degradation under future climate scenarios4, there is a pressing need to assess how climatic and environmental factors affect riverine CO2 emissions, and what changes have been seen in riverine CO2 emissions across the NH permafrost region.

Here, we conducted field observations on the Tibetan Plateau and synthesized previous measurements in the NH permafrost regions, presenting a compilation of 5685 in situ riverine CO2 partial pressure (pCO2) data between May and October from 2000 to 2020 (Fig. 1, Table S2). We quantified the relative importance of environmental variables influencing riverine pCO2 and emission across the Arctic and Tibetan Plateau. We then upscaled the riverine pCO2 and efflux during 2000 to 2020 in NH permafrost regions using a random forest model. To estimate the impacts of riverine CO2 emissions on land carbon sink, we synthesized the model outputs of land carbon uptake from the Coupled Model Intercomparison Project Phase 6 (CMIP6). Finally, we examined the spatio-temporal changes of riverine CO2 emissions from 2000 to 2020, and quantified the contributions of environmental factors to the inter-annual variations using a statistical model.

Fig. 1: Observed riverine pCO2 and their influencing factors across the Northern Hemisphere permafrost region.
figure 1

a Locations of the observed riverine pCO2 across permafrost zones. The bar chart shows the comparison (Mann–Whitney U test) of pCO2 between the Arctic and Tibetan Plateau (TP). The background of Northern Hemisphere permafrost distribution is from ref. 2. b–c Distribution patterns of riverine pCO2 and CO2 efflux with different permafrost zone types in the Arctic and TP. The scatter plot illustrates the magnitude and distribution density of riverine pCO2 and CO2 efflux of individual samples. The upper and lower ends of boxes denote the 0.25 and 0.75 percentiles, respectively. The upper and lower whisker caps indicate the 1.5 interquartile range of upper and lower quartile, respectively. Dots outside whiskers indicate outliers. The black line indicates the median value, and the black square represents the mean value. d Schematic illustration of factors influencing spatial distribution riverine pCO2 across the Northern Hemisphere permafrost zone.

Full size image

Results and discussion

Riverine pCO2 and efflux

Riverine pCO2 and CO2 efflux exhibit distinct patterns across different permafrost types in the Arctic and the Tibetan Plateau (TP). River pCO2 in the Arctic (1,532 μatm) was significantly higher than in the TP (470 μatm) (p < 0.05) (Fig. 1a). The average river pCO2 in the Arctic was higher among the continuous, discontinuous and sporadic permafrost zones compared to the isolated permafrost zone (Fig. 1b). In the TP, river pCO2 was similar across the permafrost zones. The continuous permafrost zone in the Arctic and TP had high riverine CO2 efflux (1,167 g C m⁻2 yr⁻1, 187 g C m⁻2 yr⁻1) compared to the discontinuous (375 g C m⁻2 yr⁻1, 53 g C m⁻2 yr⁻1), sporadic (285 g C m⁻2 yr⁻1, 12 g C m⁻2 yr⁻1) and isolated (119 g C m⁻2 yr⁻1, 58 g C m⁻2 yr⁻1) permafrost zones (Fig. 1c). These observations are in line with the regional observations of Arctic riverine CO2 effluxes11,27,28. The difference in spatial distribution of CO2 efflux across different permafrost types may be attributed to variations in soil organic carbon content29 and respiration rates30, highlighting the critical influence of permafrost distribution on riverine CO2 emissions across the whole NH permafrost region.

To upscale the riverine CO2 emissions on a hemisphere scale, we examine the effects of environmental factors on spatial distribution of riverine pCO2 using 15 environmental variables including six categories of climate, terrain, permafrost, vegetation, human and soil properties (Methods). Results show that soil properties have the greatest influence (44.0%) in the NH permafrost regions, followed by permafrost (21.4%) and climate (21.3%) (Fig. 1d). In the Arctic, soil properties (37.5%) have the highest contribution to the riverine pCO2 (Figure S1a), whereas the effect of climate (28.6%) in the TP is comparable to the terrain (23.8%) and permafrost properties (22.9%) (Figure S1b). Our results suggest different controls on distribution of riverine CO2 emission between high-latitude and high-altitude NH regions, and enable more detailed assessments of riverine CO2 emissions across the region.

Spatial distribution of riverine CO2 emissions

Based on simulated riverine pCO2 from the random forest model (Figure S2), and gas transfer velocity and global monthly river surface area calculated by Liu et al.12,31 using the HydroBASIN level 4 framework dataset32, we calculated riverine CO2 emissions and their spatial distributions. Riverine CO2 emission rates show a range of between 14 and 41 Tg C per month in the NH permafrost regions (Table S1). Unlike previous studies18, we estimated the riverine CO2 emissions by separately upscaling pCO2 and CO2 efflux in the Arctic and TP permafrost region considering the differences in environmental controls between the two regions (Figure S1). Moreover, compared with the global-scale study13, significantly more monitoring sites with 5685 monthly data were synthesized for accurate estimates of CO2 emission for the Arctic and TP. Taken together, our results suggest that average riverine CO2 emissions in May–October from 2000 to 2020 amount to 200 ± 15 Tg C yr⁻1, with 188 ± 14 Tg C yr⁻1 in the Arctic and 11.7 ± 0.8 Tg C yr⁻1 in the TP (Fig. 2). Our result is similar to the Arctic riverine CO2 emissions during the same period calculated from an earlier global study (202 Tg C yr⁻1)12 and Northern Hemisphere (226 Tg C yr⁻1)18 studies, and slightly lower than the riverine CO2 emissions in the TP during 2010s (18.6 Tg C yr⁻1)10. This may be attributed to the differences in study area between this and previous works. Specifically, our study focuses on the rivers within the TP permafrost region, whereas earlier studies encompassed the whole TP, which may cause the slightly discrepancies in the estimates. Notably, compared with the total estimation using bootstrapping approach18, our result shows the annual spatiotemporal distribution of riverine CO2 emission and regional disparities in the NH permafrost region. The estimated riverine CO2 emissions across the NH permafrost regions account for 10–11% of global riverine CO2 emissions12,13, which is slightly lower than the area proportion of northern cryosphere rivers (20%, Table S1) to the area of global rivers12,33. The results imply that riverine CO2 release play a critical role in the permafrost carbon feedback to climate warming. Our findings emphasize the importance of riverine carbon cycling in the permafrost regions through comparing and quantifying the river area and carbon emission on a global scale. Moreover, our study shows a significant advancement in the simulation by supplementing more observed CO2 data on the TP and considering the impacts of permafrost conditions and other environmental variables. Most significantly, given the regional discrepancy between the Arctic and TP cold regions23, our regionalized simulation in the Arctic and TP can better reveal the spatial distribution of riverine CO2 emissions.

Fig. 2: Riverine CO2 emissions and land carbon uptake between 2000-2020.
figure 2

a The distribution of riverine CO2 emissions at a spatial resolution of 1 km between 2000 and 2020. The land boundary of the Tibetan Plateau is derived from ref. 72, while the boundaries for other regions are obtained from Natural Earth (https://www.naturalearthdata.com/). b–g Annual riverine CO2 emissions and land CO2 uptake during May to October in Alaska, East Siberia, Tibetan Plateau, West Siberia, North Europe and Canada, and (h) in the whole Northern Hemisphere permafrost regions. The error bars represent the 95% confidence interval. The percentages represent the average proportions of land CO2 uptake offset by riverine CO2 emissions during 2000-2020. The dashed lines indicate the linear fit for the annual offsetting ratios through the years and the S represents the changing trends. The shadow represents a standard deviation of riverine CO2 emissions offsetting carbon sink. The significances are given as *P < 0.05, **P < 0.01, and ***P < 0.001.

Full size image

To estimate the relative importance of riverine CO2 emission in the larger carbon cycle, we used the outputs of land carbon uptake between May and October from the Coupled Model Intercomparison Project Phase 6 (CMIP6) including CMCC–CM2–SR5 and CMCC–ESM2. We show that average land carbon uptake between May and October from 2000 to 2020 in the NH permafrost regions was –711.6 (–783.8 to –639.5) (mean with the 95% confidence interval) Tg C yr⁻1, with –661.5 (–733.4 to –589.7) Tg C yr⁻1 in the Arctic and –50.1 (–51.9 to –48.4) Tg C yr⁻1 in the TP (Figure S5). Unlike the spatially inexplicit emission estimate in a previous study18, we show that riverine CO2 emissions offset on average 28.1 ± 2.1% of annual land carbon sink across the whole NH permafrost regions, with significant regional differences in the offset ratios: 24.7 ± 2.3% in Alaska, 63.1 ± 3.8% in East Siberia, 32.3 ± 2.3% in Canada, 39.3 ± 1.9% in North Europe, 13.1 ± 0.8% in West Siberia and 23.4 ± 1.6% in TP (Fig. 2). Notably, we find that the rate of increase in the proportion of land carbon uptake offset by riverine CO2 emissions in the Arctic is about ten times higher than that in the TP (Figure S7). Although uncertainties in estimated riverine CO2 emission and land carbon uptake, and exact numbers should be treated with caution, these data suggest that riverine CO2 emissions constitute a significant component in the carbon cycle of NH permafrost regions. Similarly, field observations in Western Siberia have shown that CO₂ emissions from inland waters (including rivers, lakes, and ponds) offset 35–50% of the terrestrial carbon sink11. On a global scale, the estimated current regional ratio of annual riverine CO2 emissions to terrestrial gross primary production can reach up to 18%, with high ratios in the Arctic and TP12.

Changes in the emissions between 2000 to 2020

To reveal the changes of riverine CO2 emissions over time, we simulated the riverine pCO2 year by year from 2000 to 2020 using the time-varying watershed variables (Methods) identified to influence riverine pCO2. Results show that the riverine CO2 efflux increased with an average rate of 1.05 g C m⁻2 yr⁻1 (1248 to 1283 g C m⁻2 yr⁻1, Fig. 3) in the NH permafrost zone between 2000 to 2020. In 2020, the river CO2 emissions are highest in both the Arctic and TP permafrost regions, thus we compare anomalies relative to the multiyear mean CO2 efflux during 2000–2020 and key environmental conditions to identify drivers. We further find that increases in riverine CO2 efflux are observed in about 58% of the NH permafrost area from 2000 to 2020, with the three hotpots of Central Siberia, Northeastern East Siberia and South-central TP (Fig. 3a). Notably, we find the temporal and spatial variations of riverine CO2 emissions have similar trends with precipitation and permafrost ground temperature. Specifically, the anomaly emission hotspots in 2020 primarily overlap with regions of high precipitation in the Arctic (Fig. 3b), and of high permafrost ground temperatures in the TP (Fig. 3c). Additionally, these hotspots have also extensive peatlands with vast organic carbon stocks34, likely contributing to the pronounced responses in riverine CO2 emissions to changes in precipitation and permafrost degradation in these regions.

Fig. 3: Simulated temporal and spatial variations in riverine CO2 emissions during 2000 to 2020.
figure 3

a CO2 efflux anomalies (FCO2) in 2020 relative to the multiyear mean CO2 efflux during 2000 to 2020. b–c Precipitation and permafrost ground temperature at 1 m depth (GT) in 2020 relative to multiyear mean values from 2000 to 2020. The dashed boxes show the hotpots of their changes in Central Siberia, Northeastern East Siberia and South-central Tibetan Plateau. d–f Changes of annual riverine CO2 emissions during 2000-2020 in the Northern Hemisphere (d), Arctic (e), Tibetan Plateau (f) permafrost zone. The dashed lines indicate the linear fit for the riverine CO2 emissions and the S represents the changing trends. The shadow represents a standard deviation of riverine CO2 emissions. The significances of trend are given as *P < 0.05, **P < 0.01, and ***P < 0.001.

Full size image

In the past 20 years, riverine CO2 emissions has increased at a rate of 0.42 ± 0.16 Tg C yr⁻1 in the NH permafrost regions (Fig. 3d), with a rate of 0.33 Tg C yr⁻1 in the Arctic and 0.08 Tg C yr⁻1 in the TP (Fig. 3e, f). Notably, a robust increase of riverine CO2 emissions in the TP between 2010-2020 is consistent with the recent drastic mountain permafrost degradation5,23. Given the expected/projected changes in precipitation and permafrost ground temperatures over the next decades35,36, our result suggests a significant acceleration of riverine CO2 emissions in the future. The results of this study may also be subject to various sources of uncertainty. Although machine learning-based methods are widely used in upscaling and simulating regional carbon fluxes37,38, the representativeness of riverine pCO2 measurements and the inadequate selection of explanatory variables will cause uncertainty in the riverine CO2 emission estimations. For example, current observations do not involve the central and eastern Siberia. Moreover, although total 15 environmental factors were used in the Random Forest model, the explanatory variables applied in the random forest models may not be sufficient to adequately capture the changes in riverine CO2 emissions, such as the river floodplain parameters. Taken together, despite the various sources of uncertainty, our estimates of gridded CO2 emissions from rivers have reasonable accuracy (Figures S2–S4) and uncertainty of the annual riverine CO2 emissions are relatively low.

Drivers of inter-annual variations

To further clarify the drivers for changes in riverine CO2 efflux in the Arctic and TP over the recent 20 years, we quantified the contributions of environmental factors to the inter-variability and trends using a statistical model (Methods)39. The results show that climate factors (57.3%), in particular precipitation (Fig. 4a), influence the most of interannual variability of riverine CO2 efflux in the Arctic, followed by permafrost characteristics (34.4%). This is attributed to the fact that air temperature and precipitation play crucial roles in the key processes of terrestrial carbon budget40. More precipitation enhances the land-water connections, resulting in more lateral exports of terrestrial organic and inorganic carbon24. Meanwhile, increasing temperature can cause the thawing of permafrost, enhancing soil respiration and decomposition of organic carbon in rivers41,42. It is also supported by that soil respiration contributes most to the spatial distribution of riverine pCO2 in the Arctic (Figure S1a). Our results imply that the enhanced CO2 emissions from the Arctic rivers are mainly influenced by the direct and indirect effects of climate change.

Fig. 4: Environmental drivers of trends in riverine CO2 release.
figure 4

Upscaled CO2 efflux annual anomalies (solid lines) and linear fit (dashed lines) using all drivers and iteratively excluding the temporal dynamics for each group of drivers in the Arctic (a) and Tibetan Plateau (b). The bar and doughnut chart presents the contribution percentage of climate (precipitation, air temperature), permafrost factor as ground temperature at 1 m–depth (GT), active layer thickness (ALT), vegetation and human factors. The size of the bar chart also represents the difference in slope between each test and the upscaled CO2.

Full size image

The increasing trend of riverine CO2 efflux in the TP is mainly driven by permafrost factors (55.0%) including ground temperature at 1–m depth and active layer thickness, followed by climate (20.9%) (Fig. 4b). Given that permafrost thawing is more dramatic in the TP than in the Arctic43, it makes sense to infer that permafrost degradation is the dominant driver of carbon release. Permafrost ground ice melting promotes rainfall infiltration and increases groundwater storage, enhancing groundwater discharge and carbon export into rivers44,45. Particularly, dissolved organic carbon can be mineralized or transported after being released from soils affected by permafrost degradation46. The whole TP is warming and wetting since 2000, while the increase in precipitation can further accelerate the thawing of permafrost47,48. Although these processes suggest that riverine CO2 emissions is accelerated by permafrost degradation, the influencing mechanism still need to be explored.

Implications

Our study shows the spatio-temporal distribution of riverine CO2 emission in the NH permafrost regions and quantifies its impact on the terrestrial carbon sink at an interdecadal scale. The findings yield valuable insights with far-reaching implications for the understanding of the role of riverine CO2 emission on the permafrost carbon feedback. The results not only advance our knowledge of the impacts of river carbon dynamics on the land carbon sink, but also highlight the urgency of integrating these emissions in ESMs in order to make accurate assessments of the contemporary and future climate change49.

By integrating field observations, machine learning modeling, and mechanistic analysis, our study has laid the groundwork for accurate assessment of riverine CO2 emissions. Yet, due to the limitations in in–situ observations, our estimates of CO2 emissions are constrained to a monthly scale. However, significant CO2 fluxes can occur over a short period in the Arctic, particularly during spring freshets50. Additionally, although we consider the seasonal variations in river surface area using monthly resolved open water surface areas, the temporary expansion of river surface area such as floodplain inundation during high-flow seasons is not included in the estimates51. Thus, future refined monitoring data is essential for the accurate estimation of riverine carbon release. It is imperative to further expand land-water monitoring networks, continue data synthesis efforts, and better integrate field observations with ecosystem models. The uneven distribution of observations across sites results in additional uncertainties, thus assessing the observation’s representativeness within the heterogeneous permafrost regions and improving site-level coverage will be crucial to constrain both the contemporary and future net carbon budgets in the Northern Hemisphere permafrost regions.

Methods

Observed river CO2 data

Field observation of riverine CO2 emissions was conducted in August and September of 2017 and from June to October of 2018 in the Heihe River Basins, northeastern Tibetan Plateau (TP). CO2 concentrations were measured at 93 sites (133 observations) using a LI–8100 CO2/H2O gas analyzer (LI–COR) with a self–made floating chamber (Table S2). At each sampling site three replicate 30 second measurements were carried out. Also, water temperature (Tw, ◦C), pH, river flow velocity (FV), atmospheric pressure, dissolved oxygen saturation, turbidity and conductivity were measured at each site. Based on the atmospheric pressure at each site, we calculated the pCO252.

For CO2 observations in Northern Hemisphere (NH) permafrost rivers, our study primarily relies on the GRiMeDB dataset53, dataset from Liu et al.31 and other published data10,11,16,54. Literature data were obtained from the Web of Science using following keywords “stream,” “river,” “carbon dioxide (CO2),” “pCO2,” “CO2 flux,” “carbon emissions,” “carbon release,” “water chemistry,” “water quality,” “permafrost region,” “Tibetan Plateau,” and “Northern Hemisphere”. Given the data limitation of ice-off period in the study area, we used data from May to October. In total, we obtained in-situ pCO2 data at 4917 sites (Fig. 1a), with 9096 daily measurement data of riverine pCO2 (Supplementary Materials). For multiple observations at the same site, we averaged the data to ensure that each site had a value for each month. Consequently, a total of 5685 monthly observed riverine pCO2 data were used in this study, with 5388 data in the Arctic, and 297 data in the TP. For CO2 concentration data obtained from the GRiMeDB dataset, we converted the unit from μmol L⁻1 to μatm according to the Eqs. (1–3)53. The van’t Hoff Eq. (4) modifying Henry’s law constant to account for the effects of temperature, was used to calculate the solubility of gases at non-standard temperatures.

$${{{\rm{CO}}}_2}\,{{\rm{concetration}}}(\mu mol/L)={k_H}\times {{\rm{Pressure}}}(atm)\times {{{\rm{CO}}}_2}\,{{\rm{concentration}}}(ppm)$$
(1)
$${{{\rm{CO}}}_2}\,{{\rm{concentration}}}(\mu mol/L)={k_H}\times p{{{\rm{CO}}}_2}(\mu atm)$$
(2)
$${{\rm{Pressure}}} \, (atm)={(1-(0.0000225577\times {{\rm{Elevation}}}(m)))}^{5.25588}$$
(3)
$${k_H}={{k_H} \, _{prime}}\times {e \,}^{C\times \left(\frac{1}{{{\rm{Temperature}}}(K)}-\frac{1}{298}\right)}$$
(4)

where kH prime is the original Henry’s law constant (3.4×10⁻2) (mol L-1 atm-1) (http://www.henrys-law.org), C is the temperature correction constant 2400 used for the van’t Hoff equation, temperature is the given water temperature and 298 is the standard temperature.

Water surface area of streams and rivers

We used a global monthly river surface area data created by Liu et al.12, produced under the HydroBASIN level 4 framework32. This dataset separates land surface area and accounts for both downstream and at–a–station hydraulic geometry, providing monthly surface area estimates that include streams as small as 0.3 m in width. We then clipped this dataset to the extent of the NH permafrost regions, calculating the percentage of permafrost coverage within each basin. By multiplying this percentage by the total river surface area of that basin, we derived the monthly river surface area for the NH permafrost regions during May to October. The final monthly river surface area ranges from 69,737 to 200,304 km2 (Table S1). Additionally, we used the Global Reach-level A priori Discharge Estimates for SWOT (GRADES) river networks55 (http://hydrology.princeton.edu/data/mpan/MERIT_Basins/), a new representation of global river networks derived from the fine resolution (~90 meters) Multiple–Error Removed Improved-Terrain Digital Elevation Model (MERIT DEM), and related hydrography datasets56, as the hydrographic infrastructure for riverine CO2 estimates.

Predictor variables

For modeling riverine pCO2 (μatm) in the NH permafrost regions, we collected 25 variables from global geospatial datasets (Table S3). These datasets included climate variables (monthly air temperature at 2 m (K) and monthly precipitation (mm)), vegetation parameter (monthly gross primary production (g C m−2 d−1)), yearly permafrost variables (active layer thickness (m) and permafrost ground temperature at 1 m depth (K)), topographic variables (elevation (m) and slope (%)), and 15 soil variables from the Harmonized World Soil Database (HWSD, v1.2), including gravel content (%vol.), sand fraction (% wt.), silt fraction (% wt.), clay fraction (% wt.), textures (name), bulk density (kg dm−3), soil organic carbon (SOC) (% weight), pH (−log(H+)), cation exchange capacity (CEC) (cmol kg−1), base saturation (%), total exchangeable bases (TEB) (cmol kg−1), calcium carbonate (CaCO3) (% weight), gypsum (CaSO4) (% weight), sodicity (%), salinity (dS m−1), soil respiration rate (g C m−2 yr−1), and peatland carbon stocks (hg m−2). The yearly population density (people km−2) was from WorldPop (https://hub.worldpop.org/). Utilizing stepwise regression, we addressed multicollinearity by removing several predictor variables with high intercorrelation (Variance Inflation Factor > 10) We retained the variables deemed most representative of riverine carbon release processes. Finally, we used the soil variables including soil organic carbon (SOC) (% weight), pH (−log(H+)), cation exchange capacity (CEC) (cmol kg−1), total exchangeable bases (TEB) (cmol kg−1), and bulk density (kg dm−3). Details of the environmental variables and their sources were shown in Table S3. All geospatial data were unified to the same spatial resolution (1 km) before the upscaling of pCO2 to facilitate subsequent comparisons and modeling. Specifically, we applied bilinear interpolation57 to resample the geospatial data to a 1–km resolution. This method, commonly used for continuous data, provides smoother transitions and superior performance compared to other techniques, such as nearest neighbor interpolation. Although some error is inevitably introduced during the resampling process, this approach is deemed acceptable considering the research objectives and the need for consistency across datasets.

Spatial upscaling of riverine pCO2

The Random Forest model was used to simulate the yearly and monthly riverine pCO2 in the NH permafrost regions. Compared to other ML algorithms, RF has shown to have better accuracy and lower uncertainty58. This approach has been previously applied to upscaling CO2, CH4 and N2O fluxes in rivers59,60. Considering the limited observations and highly heterogeneous landscapes between the Arctic and the TP, we separately conducted the simulation in the two regions. Firstly, for the yearly riverine pCO2 upscaling, the observed riverine pCO2 and corresponding environmental factors in the Arctic were divided into a 7:3 ratio and used for model training and testing. The 7:3 split is a commonly used practice in machine learning applications, particularly for datasets with sufficient sample sizes, as it provides a robust test set while ensuring adequate training data for model learning61. In the TP, the 8:2 ratio was adopted to maximize the training dataset, which was essential for capturing the region’s complex and heterogeneous environmental characteristics despite the data limitations. Then, we extended the default hyperparameter ranges of random forest model, including the minimum size of data points before splitting a tree (min_n), the number of variables to use in each split (mtry), and the number of trees (n_trees). Meanwhile, random search and standard 10-fold cross-validation were used to obtain a stable performance model, with R2 and Root Mean Square Error (RMSE) used for the evaluation of model performance. Finally, the number of predictive variables to use for the min_n, mtry, and n_trees were set to be 10, 1, 1000 in the Arctic, respectively, and 12, 4, 1000 in the TP, respectively. The performances of models were shown in Figure S2. The models yielded good performance for yearly riverine pCO2 both in the Arctic (R2 = 0.76; RMSE = 419.63 µatm; P < 0.001) and in the TP (R2 = 0.47; RMSE = 203.48 µatm; P < 0.01).

In addition, to quantify the variations in riverine pCO2 between May and October, we constructed monthly models using the hyperparameter of corresponding yearly model. Unlike the yearly model, the training and testing data used for monthly model incorporated observations from the corresponding and adjacent months, which was employed to increase the observations available for each model. The modeling framework yielded reasonably good seasonal riverine pCO2 both in the Arctic (R2 = 0.19–0.88; RMSE = 287–780.73 µatm; P < 0.001, 0.01, 0.05) and the TP (R2 = 0.40–0.96; RMSE = 42.48–244.14 µatm; P < 0.001) permafrost regions (Figure S3).

To upscale riverine pCO2 across the permafrost regions, the trained yearly and monthly random forest models were coupled with the relevant gridded predicting variables. Then, the yearly riverine pCO2 of the Arctic and TP permafrost regions were estimated from 2000 to 2020, while monthly pCO2 was estimated from May to October.

Changes of riverine pCO2 from 2000 to 2020

We trained the Random Forest model using the complete data from 2000 to 2020. To ensure accurate temporal and spatial alignment, we matched each pCO2 observation with its corresponding environmental variables for the specific time and location, which improved the relationship between pCO2 and its predictors, thereby enhancing model performance. Using the trained model, we substituted the environmental variables with year-specific values to predict pCO2 on a raster–by–raster basis for each year between 2000 and 2020. Key variables that vary annually include temperature, precipitation, GPP, population density, ALT, and permafrost ground temperature at 1 m depth. Since the ALT and permafrost ground temperature at 1 m depth data only cover the period from 2000 to 2019, we used the values from 2019 as proxies for 2020. All other variables, such as slope and soil properties, remained constant over time. This method has also been successfully used to estimate carbon and/or water fluxes in different periods37,62,63.

Drivers of the riverine pCO2

To quantify the relative importance of each variable for riverine pCO2, we conducted the statistical analyzes following the three steps. Firstly, we employed the “varImp” method to obtain the feature importance ranking of yearly predictive models using the “Caret” package in R version 4.2.2. Secondly, the environmental variables were divided into six major categories, including climate (air temperature at 2 m (K), monthly precipitation (mm)), terrain (elevation (m) and slope (%)), permafrost (ground temperature at 1 m depth (K) and active layer thickness (m)), vegetation (monthly gross primary production (g C m−2 d−1)), human (population density (people km−2)) and soil (soil respiration rate (g C m−2 yr−1), peatland carbon stocks (hg m−2), soil organic carbon (% weight), soil bulk density (kg dm−3), soil pH (−log(H+)), soil cation exchange capacity (cmol kg−1), and soil total exchangeable bases (cmol kg−1). Thirdly, we calculated the percentage of each category based on the relative importance of each variable. Besides, to quantify the relative importance of each environmental variable for the whole NH permafrost region, we synthesized the records from the Arctic and TP as the input data, and constructed a random forest model suitable for upscaling the riverine pCO2 in the NH (Figure S4).

Estimates of riverine CO2 fluxes

We calculated the riverine CO2 efflux (FCO2, g C m⁻2 yr⁻1) across the water–air interface using the following equation14:

$$FC{O}_{2}=k\times ({[C{O}_{2}]}_{w}-{[C{O}_{2}]}_{a})\times 12\times {{\rm{N}}}$$
(5)

where k is the gas transfer velocity (m d−1) obtained from Liu et al.31, [CO2]w and [CO2]a represent CO2 concentration (mmol L−1) in water and in equilibrium with a constant atmospheric pCO2 of 402 μatm, representing the average pCO2 from May to June 2016, sourced from 129 worldwide stations (https://community.wmo.int/wmo-greenhouse-gas-bulletins), consistent with the 2016 record at Tiksi station64 (404 ± 0.9 ppm). The number of 12 is the molar mass of carbon (12 g mol−1) and N is the number of days during the ice-free period (summed to be 180 d from May to October). A temperature-sensitive Henry’s Law constant (kH, in atm × L/mol) was used to convert partial pressures to dissolved CO2 concentrations65. River water temperature (Tw) is calculated using the monthly air temperature (Ta).

$$[C{O}_{2}]=pC{O}_{2}\times {k}_{H}\times 0.001$$
(6)
$$-{\log }_{10}({k}_{H})=-7\times {{10}^{5}}{{T}_{w}}({\circ } {{\rm{C}}})^{2}+0.016{T}_{w}({\circ } {{\rm{C}}})+1.11$$
(7)
$${T}_{w}=0.67\times {T}_{a}+7.45$$
(8)

The total CO2 emissions (FCO2 total, Tg C yr⁻1) from all streams and rivers in the HN permafrost regions was computed as:

$${{{\rm{FCO}}}}_{2{{\rm{total}}}}=\sum \, \left({{{\rm{FCO}}}}_{2}\times {{\rm{SA}}}\times {10}^{-12}\right)$$
(9)

where the 10-12 is the convertor factor and SA is the surface area of rivers (km2).

To better show the spatial distribution of riverine pCO2 and CO2 efflux, we used the moving window method to smooth the expanded GeoTIFF data. First, a moving window mean function was defined, which applies a specified window size (17 × 17) to each pixel, calculating the mean of its neighborhood to smooth the data. During the process, the original values at river network locations were preserved to ensure that critical data remained unaffected by the smoothing.

Quantified contributions of variables to trends

We used a statistical linear regression model39 to quantify the contributions of various factors to the trends of river pCO2 in the NH permafrost regions. Firstly, we establish a model, referred to as Modelall, which is driven by all factors. To understand the impact of different factors, we iteratively held one group constant while allowing others to vary. Specifically, we analyzed the scenarios in which climate (temperature and precipitation), permafrost (ground temperature and active layer thickness), vegetation factors (GPP), and human (population density) are constant respectively, with other factors change over time. For instance, we used Modelclimate to represent the results that held temperature and precipitation constant in 2000 and allowed other factors change over time. Similarly, we derived Modelpermafrost, Modelvegetation and Modelhuman. The differences ∆CO2 (Climate), ∆CO2 (Permafrost), ∆CO2 (Vegetation), and ∆CO2 (Human) (Eqs. 10–13) were regarded as their impacts on riverine CO2 emission. The trend differences were regarded as the contributions of each group to the increasing riverine CO2 emissions.

$$\Delta {{\rm{CO}}}_2({{\rm{Climate}}})={{\rm{Model}}}_{{\rm{all}}}-{{\rm{Model}}}_{{\rm{climate}}}={\beta}_{{\rm{Climate}}}\Delta {{\rm{Climate}}}$$
(10)
$$\Delta {{\rm{CO}}}_2({{\rm{Permafrost}}})={{\rm{Model}}}_{{\rm{all}}}-{{\rm{Model}}}_{{\rm{permafrost}}}={\beta} _{{\rm{Permafrost}}}\Delta {{\rm{Permafrost}}}$$
(11)
$$\Delta {{\rm{CO}}}_2({{\rm{Vegetation}}})={{\rm{Model}}}_{{\rm{all}}}-{{\rm{Model}}}_{{\rm{vegetation}}}={\beta }_{{\rm{Vegetation}}}\Delta {{\rm{Vegetation}}}$$
(12)
$$\Delta {{\rm{CO}}}_2({{\rm{Human}}})={{\rm{Model}}}_{{\rm{all}}}-{{\rm{Model}}}_{{\rm{human}}}={\beta}_{{\rm{Human}}}\Delta {{\rm{Human}}}$$
(13)

Statistics

In this study, the Ordinary Least Squares (OLS) method was used to calculate the magnitude of the linear fit trend. We used the slope of the regression model to quantitatively assess the magnitude and direction of the trends. To assess the statistical significance of the fitted results, we employed the Mann-Kendall Test, a nonparametric trend test method that has been widely used in hydrometeorological time series analysis66. This approach was used to estimate the level of significance of the interannual trends in riverine CO2 emissions and their offset of land carbon uptake. To evaluate the significance of the fit between the model predictions and the observed values, a significance test of the Pearson correlation coefficient was conducted. This was performed using the (cor.test) function in R, which calculates the p-value based on a t-distribution to assess the statistical significance of the correlation. We considered the differences to be statistically significant when P < 0.05.

Terrestrial net ecosystem exchange

To estimate the impacts of riverine CO2 emissions on land carbon sinks, we synthesized contemporaneous model outputs of land carbon sinks from the Coupled Model Intercomparison Project Phase 6 (CMIP6)67. This is because compared to the remote sensing datasets such as SMAP L4 Global Daily 9 km EASE–Grid Carbon Net Ecosystem Exchange, version 7 product68 (SPL4CMDL) and a terrestrial net ecosystem exchange dataset inferred from the ACOS GOSAT v9 XCO2 retrievals (GCAS2021)69, CMIP6 data has a long temporal period during 2000–2020. While the machine learning-based Global 250 m Terrestrial Ecosystem Net Primary Productivity Remote Sensing Big Data Estimation Product (2001–2020)70 provides a longer temporal range, its annual resolution made it difficult to align with the monthly temporal scale of the river carbon emissions. Thus, we used the outputs of land carbon sinks from CMIP6 to estimate the percentage of riverine CO2 release offsetting land carbon sink between 2000 and 2020.

Furthermore, considering the large variability in Net Ecosystem Productivity (NEP) results across models, we chose two available models (CMCC–CM2–SR5 and CMCC–ESM2) with same experimental settings, the finest spatial resolution (100 km) and the availability for both historical and future projections. Since historical simulations end in 2014, they were extended to 2020 using SSP2–4.5 projections. SSP2–4.5 scenario represents a moderate emission pathway more aligned with current development trends. To ensure equal weight for each model, we used the average of the two models. We averaged the NEP data for the months of May to October each year to represent the annual land carbon exchange level. Following standard conventions, NEP values were converted to Net Ecosystem Exchange (NEE) by taking their negative, where positive NEE indicates carbon emissions, and negative NEE indicates carbon uptake. Then we used Python for data cropping and resampling to ensure consistency in spatial resolution and alignment with the study region. To quantify the offset by riverine CO2 emissions, we excluded river networks by setting grid cell values containing rivers to zero when calculating land carbon uptake. The 95% confidence interval of the estimated land carbon uptake was used to represent uncertainty.

Uncertainty analysis

For each year, we conducted 100 iterations of random forest (RF) model training. Each RF model was trained with the same optimized hyper–parameters and different bootstrap samples (i.e., random sampling of the original data with put–back). These 100 models were then applied to the spatial grid predictors to generate pCO2 predictions for each iteration. The mean of the 100 model predictions was used to estimate the spatially upscaled pCO2 values, while the standard deviation of these predictions quantified the associated prediction uncertainty71. Finally, the predicted pCO2 mean and uncertainty were multiplied by the area to quantify CO2 emissions and their uncertainties.

Software used for the analysis

For data analysis and organization, we used SPSS software (available at https://www.ibm.com/spss) and R version 4.2.2 (available at https://www.R-project.org/), resampling and cropping, etc. were done using python version 3.11 (available at https://www.python.org/), model construction, importance ranking and visualization were done using R statistical software. GIS processing and presentation was done using ArcMap 10.8. (available at https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview)