Recent intensified erosion and massive sediment deposition in Tibetan Plateau rivers

Domain: earth_science
Published: 2024-01-24
Editorial outcome: minor_revision

Abstract
Recent climate change has caused an increase in warming-driven erosion and sediment transport processes on the Tibetan Plateau (TP). Yet a lack of measurements hinders our understanding of basin-scale sediment dynamics and associated spatiotemporal changes. Here, using satellite-based estimates of suspended sediment, we reconstruct the quantitative history and patterns of erosion and sediment transport in major headwater basins from 1986 to 2021. Out of 13 warming-affected headwater regions, 63% of the rivers have experienced significant increases in sediment flux. Despite such intensified erosion, we find that 30% of the total suspended sediment flux has been temporarily deposited within rivers. Our findings reveal a pronounced spatiotemporal heterogeneity within and across basins. The recurrent fluctuations in erosion-deposition patterns within river channels not only result in the underestimation of erosion magnitude but also drive continuous transformations in valley morphology, thereby endangering local ecosystems, landscape stability, and infrastructure project safety.

Introduction
Rivers originating from the Tibetan Plateau (TP) not only serve as crucial water resources but also play essential ecological roles1,2. They transport indispensable terrigenous materials—including freshwater, sediment, carbon, and nutrients—fostering the vitality of local ecosystems and supporting the livelihood of nearly two billion people in neighboring downstream regions. Yet, these river systems, uniquely positioned at the world’s highest elevation within the mid-latitude zone, are now confronted with intensified threats from amplified climate change, rapid degradation of the cryosphere, and emerging human activities2,3. With the TP standing as Earth’s paramount yet fragile water tower4, it experiences a warming trend of 0.32°C per decade3,5, which is twice the global average. Such a pronounced warming rate, when juxtaposed with the plateau’s notable altitudinal range, accentuates the susceptibility of these rivers to intricate physical and ecological processes, such as increased warming-driven erosion6–8, deteriorating water quality9,10, and escalating risks of natural hazards3,11.

Existing studies highlight a substantial increase in water and sediment fluxes in TP rivers, with consequential effects on regional ecosystems and biogeochemical cycles1,12. For example, elevated sediment transport can curtail water clarity, potentially suppressing aquatic photosynthesis and regional primary production13. These sediments may also convey bioavailable micronutrients, influencing phytoplankton proliferation, and habitat suitability14. Notably, shifting sediment dynamics can reshape terrestrial landscapes15, driving riverine reconfigurations16, and destabilizing landscapes throughout the region17,18. Whereas riverine suspended sediment assessments are pivotal, reconstructions specific to the TP and its surrounding regions are notably absent due to insufficient high-resolution monitoring18, even as indications underscore their markedly elevated sediment flux relative to other cryospheric regions1,17. For example, between 1995-2015, TP rivers discharged sediment at rates roughly 1.8-fold13,19 and 4.5-fold20,21 than Greenland and pan-Arctic rivers, respectively. Furthermore, sediment transport model development is limited due to a lack of comprehensive and publicly high-quality available data22,23. For instance, there are limited long-term hydrological stations, with less than ~30% of TP rivers being consistently monitored (Supplementary Table 1)1,12. These in situ records primarily exist at discrete locations of the headwater outlets, suggesting that the warming-driven sediment dynamics within the basins are still unknown. In addition, records from hydrological stations cannot fully reflect the actual basin-scale sediment source-to-sink processes24. Therefore, there is an urgent need for a comprehensive spatiotemporal view of changes in sediment transport processes on the TP.

Recent advancements in remote sensing techniques13,25,26, data availability27, and computational resources28 have facilitated the use of satellite imagery to enhance our understanding of sediment dynamics in ungauged headwaters. Here, we combined remote sensing and in situ river monitoring to quantitatively reconstruct the history of suspended sediment source-to-sink processes in major headwater regions on the TP. We focused on mainstreams and major tributaries of both the large rivers of Asia (e.g., the Yangtze, Yellow, Ganges, Indus, Salween, and Mekong Rivers) and inland rivers (e.g., the Amu Darya and Tarim Rivers), with an average drainage area of >1.9×10^4 km^2. We utilized ~76,000 cloud-free Landsat images and billions of purely water river pixels to assess 36-year changes to fluvial suspended sediment concentration (SSC) and flux (Methods). We employed a suite of algorithms that are trained by the most accurate and comprehensive in situ datasets available1,17. We revealed a strong west-east spatiotemporal heterogeneity in basin-scale sediment yield and transport patterns, which is primarily driven by the combined effects of climate change and specific geomorphic processes (e.g., glacial erosion or collapses, permafrost thaw and associated landslides, rock-ice avalanches). All our findings are based on state-of-the-art advances in techniques and satellite remote sensing approaches, which can provide valuable insights into sediment mobilization, transport, and delivery in cold regions.

Results
Spatial heterogeneity in SSC and suspended sediment flux
We found strong spatial heterogeneity in SSC and sediment yield patterns (Fig. 1), which is gradually enlarged with the differences in climate change, local soil texture, and hydrological conditions. Nearly four-decade satellite records suggested that there was a marked decline in SSCs from northwest to southeast, with an area-weighted-averaged difference of 1.78±0.27 kg/m^3 (Fig. 1a; mean±SE). However, the spatial variations of SSCs were not consistent with runoff (Fig. 1b), in contrast to the global discharge pattern of suspended sediment in river outlets (where rivers flow into the oceans and major inland seas)18. Satellite-estimated sediment yields also revealed significant spatial variations across the headwater basins of the TP, as well as within individual river channels. For instance, in the TP-west basins (Westerlies area), the basin-weighted sediment yield showed an annual average ranging from 229.83±56.01 t/km^2/yr in the Amu Darya basin to 891.69±125.86 t/km^2/yr in the Indus River basin. In contrast, the eastern basins of the TP (East-Asia monsoon area) displayed a smaller magnitude, with the annual average sediment yield varying from 30.67±12.98 t/km^2/yr in the Yellow River basin to 137.55±29.14 t/km^2/yr in the Salween River basin (Supplementary Tables 2 and 3). The pronounced spatial variability in sediment yields stems from the varying response levels of distinct river systems to climate change, particularly evident in the unique warming and weathering rates28, erosion dynamics29, and the magnitude of changes in runoff1,2.

Temporal trends in SSCs across the TP rivers have exhibited higher increasing rates (6.02±2.85%/10yr; mean±SE) compared to runoff (5.53±0.51%/10yr) from 1986 to 2021 and have again revealed distinct spatial variations within each river. In highly glacierized basins in the western Himalayas, the increase rate of upstream SSCs was significantly higher (4.53±3.26 %/10yr) than downstream (1.53±0.48%/10yr), despite the increasing rate of upstream runoff (4.93±3.10 %/10yr) being lower than downstream (6.05±3.25 %/10yr). In contrast, in the East-Asia monsoon basins, most rivers have retained a semblance of their natural variability in SSC in response to runoff. In line with changes in runoff, upstream SSCs demonstrate higher increasing rates (1.50±0.67 %/10yr) than downstream SSCs (0.63±0.23%/10yr), and the timing of peak SSC remains regionally consistent with that of peak runoff. The observed spatiotemporal discrepancies between SSCs and their association with runoff suggest that hydraulic factors may not be the dominant driving force of erosion in the western region. Instead, other factors that are closely linked to the underlying surface characteristics responses to climate change, such as glacial erosion, freeze-thaw erosion of the active layer, and the pronounced influence of westerly winds during the period of low vegetation cover, may exert a more significant influence (Fig. 1). Conversely, in most areas of the eastern and southern TP, runoff continues to play an external driving role in shaping the sediment yield pattern.

High-concentration sediment export hotspots were concentrated in headwater basins dominated by glacier and permafrost (>0.5 kg/m^3, Fig. 1). The glacierization (percentage of glacier cover per headwater basin) and snowmelt runoff ratios of six exorheic rivers were evaluated, revealing a positive correlation with the average SSCs (P<0.05, Supplementary Fig. 2). At a global scale, increasing sediment yield has been the key signal of cryosphere degradation since the 1960s12,27, primarily because of climate change-driven landscape changes, such as glacier retreat, permafrost thaw, and snowpack reduction. Such landscape disturbances have altered the magnitude and frequency of glacial erosion and thermokarst dynamics17, resulting in increased heterogeneity in sediment yield patterns. For example, the headwater basins with the highest SSCs were the Indus basin (27.61% glacier cover) and Tarim basin (20.61% glacier cover), which are highly glaciated, with the annual average SSCs of 6.80±0.35 kg/m^3 and 4.06±1.24 kg/m^3, respectively (Fig. 1). These values were an order of magnitude greater than those of other ice-free erodible landscapes30.

Moreover, permafrost thaw has led to the expansion of thermokarst landscapes (e.g., thaw slump, thermokarst landslide, and thaw subsidence) and escalating sediment discharge7,12, which resulted in relatively higher SSCs in satellite imagery across permafrost regions17,31 (permafrost cover >50%, Fig. 1 and Supplementary Figs. 7 and 8). At the hillslope scale, the heterogeneity of SSCs was further amplified by the interplay between permafrost soil properties, vegetation cover, and land cover characteristics32. Specifically, basins with higher soil bulk density and sand content are prone to significant water erosion during peak rainfall events or peak meltwater events, leading to pronounced spatial gradients in SSCs (Supplementary Fig. 3). In addition, vegetation development in proglacial and permafrost thawing areas can stabilize slopes. However, this process is complicated by deglaciation and thawing-triggered slope instability, as well as the intensification of erosive rainfall29,33. At the basin scale, the interplay between sediment yield and annual precipitation followed an inverted parabolic curve, with the minimum sediment yield observed at intermediate levels of annual precipitation (400–600 mm), and an increase on both sides of this minimum (Supplementary Fig. 4). With the exception of the arid Qaidam and Inner basins, areas with sparse vegetation and low ground coverage are strongly associated with high-intensity glacier and permafrost regions, where basin-scale sediment yields appear to the highest. Notably, as influences of glacier and permafrost diminish and vegetation coverage increases, sediment yield decreases accordingly. In contrast, in the eastern areas with shrub and forest vegetation, water erosion driven by rainfall becomes the primary factor of sediment mobilization, with erosion rates increasing in tandem with higher rainfall intensity.

Temporal changes in suspended sediment flux
Over the past six decades, climate change has indisputably caused significant increases in water and sediment fluxes on the TP1,12. At the basin scale, the satellite record spanning four decades revealed a significant increase in suspended sediment flux (SSF) in 63% of the river sections (P<0.05, Fig. 2). However, this overall increasing trend masks substantial spatiotemporal variations within individual headwater basins due to limited river sampling and inadequate in situ observations17,18. Compared to previous studies1,12,34, we offer several orders of magnitude more data to elucidate the pronounced spatiotemporal heterogeneity in sediment fluxes, predominantly evidenced by alterations in sediment transport dynamics.

The sediment fluxes in the western Himalayas have increased substantially in recent years, particularly in the Indus (17.87±7.41 %/10yr; mean±SE), Amu Darya (24.57±19.32 %/10yr), and Ganges (11.27±6.42 %/10yr) basins, due to the rapid retreat of glaciers35 (Supplementary Table 3). Glacial erosion can mobilize and transport large amounts of sediment downstream (Fig. 2a), but the transport efficacy hinges on the connectivity between glaciers and river channels, as well as the glacier melt volume17,36. Satellite observations revealed that 20.98±9.59% of the glacial sediment (P<0.05) remains in the proglacial zone as moraines, debris cones, and alluvial fans, rather than being rapidly transported downstream (Supplementary Fig. 5 and Supplementary Table 3), resulting in a transport-limited environment. For instance, in the upper Indus basin, 41.67±22.87% of annual sediment supply (reported by Darband station34) is temporarily stored in the proglacial zone, particularly in glacial lakes, floodplain, alluvial fans, and wider valleys that have a relatively lower terrain ruggedness indices than downstream areas (Supplementary Figs. 6 and 7). Moreover, sediment deposited near the glacier terminus and subglacial channel may be transported downstream during intermittent ice and snow melt floods or extreme rainstorm events. Such potential spatiotemporal variations further increase the uncertainty of sediment transport dynamics, resulting in a non-significant decadal trend in the midstream areas of glacierized basins, in contrast to the significant increase observed upstream and downstream (Fig. 2).

In response to a warming climate, permafrost degradation and thermokarst activities have also increased sediment availability and transport capacity in the eastern permafrost headwaters. The annual increase rate in basin-weighted SSF for the Brahmaputra, Salween, Yangtze, and Mekong basins are 9.46±5.42 %/10yr, 18.36±8.23 %/10yr, 15.21±6.83 %/10yr, and 6.89±3.71 %/10yr (Supplementary Table 3), respectively, which can largely attribute to the expansion of new flow paths and active erodible landscapes22,36,37. The creation of new gullies and the expansion of flow paths during extreme rainstorm events have increased sediment conveyance by enhancing hillslope-channel coupling, reworking previously-stored sediment, and linking disconnected sediment sources, leading to a steady increase in sediment upstream for decades (Fig. 2 and Supplementary Figs. 5, 17). However, the signal generated by upstream permafrost thaw can be disconnected from downstream areas and the basin outlets by local sediment sinks, including expanded areas of thaw subsidence, thaw gullies, and thermokarst lakes. In arid and cold regions that comprise the upstream Yangtze, Inner (0.41±0.19 %/10yr), Yellow (4.35±3.51 %/10yr), and Hexi Corridor (1.24±0.67 %/10yr) basins, the smaller thermokarst lakes (less than 1000 m^2) account for 24% of the total TP, while larger ones (greater than 150,000 m^2) account for 50%, which potentially act as important sediment sinks, effectively offsetting the downstream transport of sediment flux (Supplementary Figs. 5 and 7). Additionally, the heterogeneity of sediment transport in permafrost basins is further enhanced by changing soil permeability (Supplementary Fig. 3)38, surface/subsurface flow paths39, and hydrogeomorphic connectivity37.

Changes in riverine erosion and deposition pattern
As demonstrated in the previous studies29,40, the changes observed in fluvial sediment in cold regions result from a combination of multiple climate change and cryosphere degradation-driven geomorphic processes. Drawing on extensive satellite data, we observed continual annual changes in riverine erosion and deposition patterns within the TP basins. Such shifts intensify the spatiotemporal heterogeneity in sediment dynamics. To gain a comprehensive understanding of the changes and to quantify the regime shift, we calculated the sediment-flux budget along the major rivers.

Compared to existing reports based on sediment monitoring stations at basin outlets1,8,12,34,41,42, we found that the overall sediment flux of each major river has been significantly underestimated (P<0.05). We conducted stringent screening to ensure data integrity, with most of our data originating from robust in situ measurements employing depth-integrated field measurement methods (Supplementary Table 2). While these data may not entirely encompass the sediment discharge trends in the large rivers, their extensive spatial coverage and temporal range enable meaningful comparisons with satellite-calibrated datasets. Specifically, the underestimation of SSC ranged from 7.81% to 43.52%, while the suspended sediment flux was underestimated between 9.01% and 46.85% (Supplementary Table 2). This is particularly true for large basins that rely solely on a single station at their outlets to represent the entire basin’s erosion amount (e.g., Xiangda station in Mekong basin and Daojieba station in Salween basin, Supplementary Fig. 10). The underestimation can primarily be attributed to the dynamic redistribution of sediment within each headwater basin, due to the occurrences of sediment deposition or a lack of transport capability at specific locations, which is not adequately captured by the available station data.

We found that sediment deposition occurred in most basins (nearly >90%, Fig. 2), and we categorized these sedimentary landscapes into three types and provided a notable catalog of locations that have experienced or have the potential to experience sediment deposition (Supplementary Fig. 8). The first type of sedimentary landscapes occurred in areas where the runoff driving force was insufficient to exceed topographic and transport thresholds, primarily in the upper and middle reaches of glacierized basins around the Himalayas, including the Indus, Tarim, Ganges, and upper Brahmaputra basins. These basins are characterized by high rock uplift, aridity, high relief, and large variability in river runoff and sediment transport43, which have led to the development of numerous large-scale structural sedimentary basins along the TP margin, such as the Thakkhola–Mustang Graben in Nepal and the Zhada basin in the upper Sutlej region44. Satellite-estimated sediment budget showed that ~51.86 Mt/yr sediment has been deposited in these basins, accounting for ~18.21% of the headwater outlets export (Supplementary Table 3). More importantly, once extreme rainstorms or meltwater floods occur, these temporarily stored sediment could be remobilized and flushed downstream, causing rapid changes in regional erosion-deposition modes and significant cryosphere hazards (e.g., glacial/landslide lake outburst floods, rock–ice avalanches, rockfalls, and debris flows). For example, in 2016, a glacial lake outburst floods (GLOF) event in the Bhotekoshi/Sunkoshi River (Ganges basin, Nepal) triggered a 30-fold increase in local sediment flux through channel-connected landslides45.

Furthermore, the second type is primarily controlled by geomorphic systems and associated hydrological connectivity, resulting in the formation of macro-distributed braided river systems and their floodplains. These sedimentary responses tend to concentrate in valley basins with steep profiles, such as the Indus, upper and middle Brahmaputra, lower Salween, and midstream and downstream areas of the Yangtze and Yellow basins. These basins are characterized by narrow gorges that alternate with wide alluvial reaches, extensive river terraces, and remnants of large lakes that are either perched or located close to the rivers (Supplementary Figs. 5 and 7). As sediment is transported to these specific locations, the river channel slope decreases, resulting in reduced sediment carrying capacity and likely deposition. This process leads to the widening and shallowing of the river channel and the accumulation of sediment onto sand bars (Fig. 3b and Supplementary Fig. 7). For instance, in the past few decades, ~21.35 Mt of total sediment (representing ~21.34% of sediment discharged downstream) has been deposited in the braided rivers and wide-narrow valley nodes of the middle reaches of the upper Yangtze River (Fig. 3c). Meanwhile, in the wide valleys of the midstream Brahmaputra River, ~45.29% of the annual sediment supply has been deposited (Fig. 3d). During the summer flood seasons, increased sediment deposition worsens the rise in the riverbed, further constraining the water cross-section and ultimately hindering sediment transport and flood discharge capacity (Supplementary Fig. 9). Our estimates indicate that over the past few decades, ~26.34 Mt/yr of sediment has been deposited at these specific geomorphic units, representing ~34.33% of the total sediment transport downstream (Supplementary Table 3).

Vegetation development can also play a role in hindering sediment transport and the expansion of erodible landscapes. In the warming TP, the overall greening and vegetation restoration can stabilize soil, reduce erosion, intercept and slow down water flow, and contribute to the formation, maintenance, and stabilization of deposition, thereby slowing erosion rates and maintaining slope stability and ecological balance12 (Supplementary Fig. 4). However, emerging human activities such as over-exploitation, overgrazing, damming, and deforestation in some local areas can cause shifts in erosion-deposition patterns in alpine mountainous areas (Supplementary Fig. 10), exacerbating soil erosion and deposition in rivers26,40. Despite the significance of changes in vegetation cover, we are unable to directly detect their contribution to changes in the sediment budget. This is particularly true given that they are highly coupled with other factors; for instance, climate change, water-sediment connectivity, and related rapid adjustments in underlying surfaces, which exert more important controls on basin-scale sediment delivery processes. Overall, a series of changes in the cryospheric landscape and associated hydrogeomorphic processes can explain the observed increase and higher spatiotemporal heterogeneity in sediment transport patterns on the TP. However, it is important to note that effective sediment export still depends on regional erosion and deposition patterns.

The operation of Earth observation satellites facilitates the high-resolution mapping of suspended sediment concentrations and fluxes in remote ungauged headwater basins. Here, we used long-term Landsat data to estimate sediment transport on the TP, to quantify and capture the complex, interconnected sediment erosion/deposition patterns induced by geomorphic changes associated with global warming. We found strong spatiotemporal heterogeneity in SSC and sediment flux patterns, which is gradually enlarged with the difference in magnitude and frequency of regional landscape changes in response to climate change. Due to enhanced glacier retreat, snowpack melting, and permafrost disturbances, high amounts of sediments were produced and transported into river channels, while differences in geomorphic patterns for each headwater basin determined sediment transport vary spatially to a large extent. ~30% of the total suspended sediment flux has been temporarily deposited in sedimentary basins, proglacial areas, and river systems with limited sediment-water connectivity or transport capacity. Once reactivated and remobilized, these previously deposited sediments can be remobilized and constitute a significant threat to the local ecology and the development of hydropower projects.

The fluvial sediment erosion and deposition we focused on were by no means confined to the TP. In fact, sedimentary and geomorphic climate change responses are expected to accelerate worldwide, especially in high mountain areas and polar regions1,17,18,40. Permafrost and glaciers are degrading rapidly in the high latitude polar and sub-polar areas such as western Norway, northern Swedish, British Columbia, and northern Canada13,32,46,47. Fluvial sediment in the Arctic induced by glacier recession and permafrost thaw has increased by many orders of magnitude48–51. If sediment fluxes increased as a result of glacier recession, permafrost degradation, and extreme weather events, the consequences downstream would be serious and manifold, including increased carbon loading, reduced water quality and storage, bed aggradation, and rising flood stage. Thus, the ability of remote sensing technology to monitor and reconstruct fluvial sediment responses to progressive climate change has broader implications for human health, hydropower development, and ecosystem security.

Methods
Data and quality control
In situ measurements. To obtain reliable suspended sediment concentration (SSC) estimates for large rivers on the TP, we collected in situ water and sediment data from pristine headwater regions (Supplementary Table 1). These data were carefully selected to avoid anthropogenic impacts (e.g., dams, reservoirs, and hydropower stations)1 and environmental noise52, and to ensure they were consistent with satellite image acquisition frequencies. We utilized long-term, depth-integrated sampling approaches conducted by local hydrological stations, which are considered the most comprehensive method. Other sample data obtained from surface water, a single depth, and automatic pump samples were not considered due to the potential for sampling technique changes to introduce uncertainty in SSC calibrations53. In situ runoff and sediment data for the Chinese headwaters are sourced from the Hydrological Data Yearbook, the Ministry of Water Resources, China (http://www.mwr.gov.cn/sj/). Data for other international rivers come from the Water and Power Development Authority, Pakistan (http://www.wapda.gov.pk/) and refs. 34,41,42,54,55. The in situ observational data have undergone stringent quality control at the local monitoring stations. Further details can be found in the refs. 1,5.

Satellite image collections. The Landsat surface reflectance tier 1 image collections were obtained from the Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land imager (OLI) sensors28. The TM and ETM+ collections consist of four visible and near-infrared (VNIR) bands and two short-wave infrared (SWIR) bands that have been processed for orthorectified surface reflectance with a spatial resolution of 30m/pixel. The OLI collections consist of five VNIR bands and two SWIR bands. To address the scan-line corrector failure in Landsat 7, we adopted the local linear histogram matching technique as recommended by the USGS, ensuring the integrity and completeness of Landsat imagery56. In light of variances in band configurations and responses among OLI, TM, and ETM+, we used designated band adjustment factors and algorithms to harmonize surface reflectance across them57. Using overlapping image dates and areas, we cross-calibrated Landsat 5, 7, and 8 reflectance values prior to preprocessing, ensuring consistency and eliminating systematic satellite biases.

The Landsat SR product provides essential information, including basic land cover classification, and quality assurance codes for each pixel in the ‘pixel_qa’ band. All image collections were atmospherically corrected and included a cloud, shadow, water, and snow mask, as well as a per-pixel saturation mask. A total of 145,908 Landsat images were applied to the Google Earth Engine (GEE, https://earthengine.google.com/) for further cloud removal, calibration, and mask extraction (Supplementary Fig. 1). Due to the limitations caused by Landsat’s 16-day repeat period and challenges with mountainous shadows and clouds, we leveraged in situ observations from recent years to calibrate and interpolate the data for each region of interest (ROI) on the TP58.

Geographic and hydrography datasets. We used the current best available global geographic and hydrological datasets to spatially couple and analyze sediment transport processes, including Multi-Error-Removed-Improved-Terrain (MERIT) DEM59, MERIT Hydro60, Global River Widths from Landsat (GRWL)27, Global Land Ice Measurements from Space (GLIMS)61, The GlObal geOreferenced Database of Dams (GOODD) dataset62. The MERIT DEM and Hydro, with a resolution of 3 arc-seconds, provided information on river flow direction, elevations, upstream/downstream links, and drainage area. We calibrated the extracted river masks and calculated the mean river width using the GRWL vector product, which has a width greater than 30 meters and was used to prevent seasonal or geomorphic disruption and evolution of each retrieved river. The GLIMS provides the most detailed global inventory of glaciers, we clipped it to obtain glacier area changes in the ROIs. We excluded the impacts of human activities by using GOODD, the largest georeferenced database of global dams, and by obtaining hydropower development data on the TP obtained from the International Hydropower Association (https://www.hydropower.org/status-report).

Supplementary datasets. Due to the limited in situ observations on the TP, we supplemented our analysis by collecting additional data to investigate the effects of climate change and geomorphic processes on the sediment yield and transport processes. These data included soil63,64, land use65, and vegetation changes66, which are related to the sediment characteristics, as well as land surface temperature67,68, precipitation69, runoff70,71, permafrost72–74, glacier35,75 and snow76, which are related to thermal and hydraulic erosion. Note that the spatiotemporal resolution of these datasets was not entirely consistent with the sediment data retrieved from satellite imagery. Consequently, we endeavored to augment our dataset by incorporating on-site meteorological observations to the greatest extent possible. We also leveraged the fifth-generation European Center for Medium-Range Weather Forecasts (ECMWF) Reanalysis (ERA5) dataset as a supplementary source of meteorological intelligence in areas where on-site station data were not available.

Preprocess of satellite imagery
Alpine environments are known for their complexity in optical images, with pixels containing various sources of environmental noise such as clouds, snow, ice, shadows, forests, and image artifacts77. Previous studies have attempted to extract surface water using Landsat images through techniques such as single-band thresholding and dual-band indexing, but these methods suffer from low accuracy and have limited applicability in alpine regions52,78. For example, in the western TP, these scenes exhibited diverse features due to elevation differences over short horizontal distances. To address these limitations, we employed multi-index methods to generate a reliable water mask for each image that overlapped one or more ROIs. Our approach included the use of several indices, such as the Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Modified NDWI (MNDWI), and Automated Water Extraction Index (AWEI)79, to improve accuracy and ensure adequate coverage of the targeted area.

To ensure high-quality water mask extraction from the satellite images, we began by eliminating images with excessive cloud cover (>20%). Next, we adopted the method outlined by ref. 80 to partition the TP region within the glacier-rich regions along the Himalayan mountain area (e.g., the Indus, Ganges, and upper Brahmaputra basins) and other cold regions. Regarding the main basins around the Himalayas, we segmented each ROI into areas above and below 665 m elevation and determined the optimal thresholds for each water index using the Otsu method81 (Supplementary Fig. 11). Furthermore, for other cold regions, we employed the MNDWI and selected the Landsat band thresholds that distinguish between water and land based on ref. 53. To address the issue of sub-pixel level inundation changes or overestimating SSC during flood events, we utilized an automated method82 to determine the water fraction in the extracted water masks and ensure the root mean square error (RMSE) of less than 7%.

To minimize the uncertainty caused by the varying geochemical properties and spectral signatures of water and sediment53, we applied an unsupervised K-Means clustering approach to classify the optical properties of each river mask. We generated our ROIs by buffering the regions of each water mask by 90m. Sub-pixel rivers (<30m wide) were not preserved, although this may have caused spatial discontinuity. We focused only on the mainstreams and large tributaries of the large headwater basins on the TP, where water flow is continuous enough to obtain a long-term sediment sequence. For each image of the ROI, we calculated monthly spectral signatures using the subsets of all the image collections and generated the K-Means groupings by calling the Cluster function (ee.Clusterer.wekaKMeans{}) in GEE (Supplementary Fig. 12).

Calculating suspended sediment concentration and flux
We compared the daily median pixel values of individual Landsat bands and band ratios. By setting a threshold for the minimum spectral feature distance from the K-Means cluster centroids, we assigned samples from each river’s monthly subset to their respective K-Means groups. This dynamic assignment allowed for flexible, seasonally sensitive estimates of suspended sediment concentration (SSC), which strongly reduced errors and bias relative to using a single model18. Furthermore, we mapped the spatial distribution of the five optimum clusters, which can capture monthly observable river sections across the entire TP (Supplementary Fig. 12a). This approach additionally guarantees the classification of neighboring river sections into distinct categories based on variations in their spectral characteristics. Consequently, it ensures that the spatial variability of SSCs does not result from differences in the assessment of a single model. After normalizing values from distinct satellite sensor bands, marked variations in band values or their rations were evident among these clusters (P<0.001, Supplementary Fig. 12b). We also plotted the true-color visuals of selected rivers based on their average RGB reflectance (Supplementary Fig. 14a). Subtle variations between each cluster will be showcased in terms of average true color juxtaposed with their respective SSC ranges.

Following Dethier et al.18, we employed least-squares multiple regression to construct the SSC calibration model. After identifying the optimal five K-Means clusters, we used variables automatically selected from a list of spectral bands and band ratios using cross-validated Lasso Regression83. We aligned 75% of in situ SSC measurements with their corresponding band values or ratios, based on the consistent sampling coordinates and timings. These associations provided the foundation for further training the SSC calibration approach according to ref. 53, wherein we selected the model with the least number of explanatory variables that fell within one standard error of the most successful iteration of the model. This form of multiple regression eliminates low-power explanatory variables to avoid over-fitting and multi-collinearity. Furthermore, we reallocated river pixel clusters to location-matched satellite imagery for each model, enabling maximal extrapolation of local calibration models due to their consistent spectral characteristics in rivers84. This suite of algorithms aims to minimize the uncertainties inherent in satellite products, and the integration of machine learning further enhances the robustness of the models (Supplementary Fig. 11). The governing equation for the algorithm is:

\beta = \text{argmin} \sum_{i=1}^{N} (SSC_i - \beta^T x_i)^2 + \lambda ||\beta||_1 (1)

where the subscript i refers to the sequence number corresponding to each cluster, and SSC is the response variable. x is the matrix of independent variables, representing the ratios or values between K-Means grouping, for example, Band 4 and 5; band ratios B2/B1, B4/B3/B2. \beta is the vector of regression coefficients, and \lambda is the regularization parameter. The first term of the objective function is the square error loss, and the second term is the L1 regularization term, which penalizes the sum of the absolute values of the regression coefficients. When \lambda is larger, the model punishes the regression coefficients more heavily, and their values become smaller, or even zero. Therefore, Lasso regression can achieve feature selection and dimensionality reduction effects.

Using this algorithm, we deduced SSC values for specific moments by extracting surface reflectance band values and ratios from regions yet to be sampled. The resulting calibration models enabled us to generate maps of the median, mean, and summary of SSCs for each ROI. We amalgamated satellite-estimated SSCs with monthly average discharge data (in m^3/s), predominantly derived from in situ station measurements. Yet, many rivers on the TP feature incomplete records, although the dataset we collected represents the most comprehensive long-term discharge dataset available for the region (Supplementary Table 1), as acknowledged in previous studies1,12. We also supplemented the discharge data by incorporating existing products specific to the entire TP70,71, further enriched with ERA5 datasets. Temporally, we generated annual discharge time series (Q_A) based on monthly discharge records to align the temporal resolution with the satellite-estimated monthly SSC dataset. In cases of data gaps for river discharge, we filled the gaps using 3 or 5-year running annual averages.

In our analysis, we focused on the mainstreams and significant tributaries, characterized by a minimum river width of 90 meters, while smaller streams without available discharge data were excluded from consideration. To improve the precision of our annual sediment flux calculations, we determined the discharge-averaged SSC for each extractable river cross-section, aggregating these values into a one-year statistical average map (referred to as SSC_A, Eq. 2). Sediment flux was computed by multiplying the annual SSC (SSC_A) by the annual runoff (Q_A) for the river segments scaled to a 5 km river network. This calculation was performed for each river cross-section, following the downstream flow from the source to the outlet. This approach allowed us to present an extensive 36-year record for the major rivers across the TP. Many of these rivers lack contemporary published sediment observations, and some are even devoid of any historical records—remarkably, historical measurements align well with our reconstructed observations (Supplementary Fig. 16).

SSC_A = \frac{1}{Q_A} \sum_{i=1}^{12} SSC_i \times Q_i (2)

where subscript i refers to the monthly average value of the SSC and discharge, the bar denotes the whole-record average of that parameter, and subscript A refers to the annual average value of runoff.

Validation and uncertainty
We aimed to minimize the uncertainty of our calibration models by utilizing unsupervised K-Means clustering, implementing river-specific corrections to the calibration models, and applying additional related parameters. We applied the procedures53 to minimize both relative error and bias during their application (Eqs. 3–4; Supplementary Fig. 13). We also calculated the complexity of root mean square error on log-transformed values. Using the median, it resists distortions from outliers—a frequent occurrence in this extensive dataset owing to factors such as sampling errors, timing discrepancies between sampling and Landsat image capture, image anomalies, and pixel misclassification.

rel:error = 10^{M(|lg SSC_{estimated i} - lg SSC_{measured i}|)} - 1 (3)

model residual_i = lg SSC_{estimated i} - lg SSC_{measured i} (4)

where SSC_{measured} and SSC_{estimated} are taken from the holdout dataset that was not used for calibration, and where M denotes the median of SSC.

We validated models for each cluster using a total of 19,690 verified pairs of satellite-estimated SSCs and in situ measurements, achieving relative errors between 0.19 and 0.32 (Supplementary Fig. 13). We also calculated the median relative error for each cluster and the relative bias with respect to in situ data for each station (Supplementary Fig. 14). Site-specific validations between monitoring stations and satellite imagery further yielded a relative error of 0.24 and an average RMSE (root mean square error) of 0.50 (Supplementary Fig. 15). We rigorously analyzed the model uncertainties, rooted in data quality, environment noise, long-term water-body changes, flood events, as well as the transferability of local calibration. We employed a suite of data quality controls and uncertainty assessment techniques to minimize uncertainties to acceptable thresholds. Further details can be found in the Supplementary Information.

Sediment budget model
Tian et al.86 developed a basin-wide Sediment Budget Method to systematically quantify sediment contributions from various tributaries and sections of the Yangtze River. We utilized this method’s universal concept to quantify the sediment budget of each river we focused on, which can be monitored by Landsat satellites for almost four decades. We improved the accuracy of point calculations (i.e., hydrological stations) by incorporating satellite-estimated sediment data, allowing us to capture the amount of deposition with greater precision. From the source to sink, we computed sediment flux variations for every 5 km segment along the river network, assessing whether sedimentation occurs or not based on Eq. 5. We centered on the mainstream and major tributaries, as signals from minor tributaries are susceptible to seasonal variations and shifts in erosion and deposition. Adhering to our flux computation guidelines, we cataloged the annual erosion and deposition patterns across the large TP rivers. We also quantified and mapped locations with either extant or potential for massive sediment deposition, combining these with the spatiotemporal SSC dataset we provided (Supplementary Fig. 8). In this method, proportional deposition (P) is a parameter used to simplify sediment deposition in river channels, which is defined by the following equations:

P = D / (SSF_a + SSF_b) (5)

where SSF_a and SSF_b are suspended sediment fluxes (Mt/yr) in river channels a and b, respectively. D represents the total amount of sediment deposition in the downstream river section after the convergence of channels a and b, and is calculated as the difference between the upstream and downstream sediment fluxes. We employed this coefficient to compute the annual erosion (or deposition) pattern in river systems, enabling us to assess the annual deposition events in the river channels. Stratifying the rivers into upper, middle, and lower sections based on their length, we systematically aggregated all measurable sediment fluxes for each respective year. This summation provided us with a comprehensive annual sediment erosion (or deposition) view of the entire river. Additional information on the Sediment Budget Method can be found in ref. 86.

Analyzing trends in suspended sediment dynamics
We utilized the Mann–Kendall non-parametric trend test to investigate changes in annual SSC and SSF for each river, and we used a P-value to assess statistical significance. To facilitate statistical analysis, we established a sampling time interval from January 1, 1986, to December 31, 2021, as Landsat 5 began imaging the TP at the end of 1985 (https://developers.google.com/earth-engine/datasets/catalog/landsat). We chose an annual sampling interval due to the large-scale spatiotemporal resolution and the need for high precision and operating speed. We obtained single-day SSC data at ROIs. To account for glaciers and permafrost areas (e.g., Amu Darya, Indus basins), we restricted the retrieval time to the thaw periods (land surface temperature >0°C). For areas with seasonally frozen ground or unfrozen areas, we did not impose time limitations. Therefore, we were able to reconstruct the sediment dynamics of each river for 36 years and analyze their trends.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Li et al. under review describes a regional-scale study in which sediment loads and fluxes were estimated for major rivers starting in the Tibetan Plateau. The study relied primarily on remote sensing imagery (Landsat) to investigate spatiotemporal trends in the headwaters of major Asian rivers, finding a general tendency towards increasing loads since the late 1980s, for which the authors (and others prior to them) make a strong argument that this is happening in association with global warming. Overall, the results of this study are of great importance to the understanding of some of the most important rivers on Earth, which support a disproportionately large fraction of the world’s population.\nThis work is certainly of significance to the field, as it largely corroborates the message that Li et al (2021; well acknowledged in this manuscript) postulated in their estimates of sediment load increases using in situ data only. That said, the study does not present a completely new message, but instead, it presents an alternative -and more comprehensive- set of evidence for sediment load increases in the Tibetan Plateau based on remote sensing imagery.\nThe work presented supports the conclusions and claims stated. It was good to see that the focus of the study was on spatial and temporal trends, which is certainly the major gap for the Tibetan Plateau, rather than on drivers, which the authors skillfully controlled for by limiting the analysis to areas of little human development.\nI could not find any evident flaws in the data analysis or interpretation of this study. However, I would encourage the authors to make their data and scripts open as soon as possible to allow reviewers and other scientists to reproduce and ratify their analysis.\nThe methodology of the study is sound, in general, although there are a few aspects that deserve some more clarity in the main text. With regards to the Landsat imagery, I was glad to read the SI subsection on Calibration and Uncertainty, but I wonder if some of these correct assertions about the uncertainty associated with remote sensing imagery analysis can be postulated in the discussion. Moreover, I did not see an explanation for how the authors dealt with the 2003 Landsat 7’s Scan Line Corrector failure, and how this might have affected the trend analysis.\nAnother important methodological aspect I was not clear about is the use of the calibration curves for the SSC calculation. The methods section broadly explains the process used, but it was not clear if all the curves shown in Fig S14 were used, or if only a selected group of them were selected based on some criteria. At least five of the curves have an error over 0.5.\nAn aspect of the work that I found particularly interesting and novel was the deposition and erosion patterns, which were not only quantified numerically, but also displayed geographically at a remarkably fine scale. These patterns, however, are purely based on a proportional deposition parameter (P), which is empirically derived from the SSF fluxes estimates. The evidence for such patterns would be much stronger if this empirical calculation was corroborated with a more physically-based representation of erosion and deposition, which are certainly well acknowledged in the geomorphology literature.\nMauricio E. Arias", "reviewer_2": "General comments\nThis paper estimated SSC and sediment flux across large rivers in the Tibetan Plateau using Landsat and then looked at trends in SSC/flux over time, over space from upstream to downstream, calculated a sediment budget to infer erosion/deposition, and compared flux estimates to factors such as glacial cover and land cover. Overall, this is a well-conceived study with excellent figures and interesting analyses.\nI think this manuscript has great potential but requires more focus and refining. My two major recommendations are 1) using more precise and specific language. Often the language used was too general, or repeating the same concept with different words, that I could not understand the authors interpretations or framing. 2) Providing more details on methods when needed and providing quality checks to give us more confidence in the data and big interpretations made from these data. See more comments below.\nSpecific comments\nTitle: If the title is “rapid changes to erosion and sediment transport” I think it could be good to show that the rates of change in flux are faster than other studies from other parts of the world. OR consider changing the title to reflect the focus of the manuscript, which appeared to be more about spatial patterns in sediment flux and changes in sediment flux in cold regions than it was about rapid changes in erosion/transport.\nValidation\nThe model validation process needs to be clearer. How many matchups were used to train models? Was 1 model built and applied to all points shown on the map? Was the same model used when applying to river pixels? Was there one model for each “cluster”? Is figure S14 the hold-out validation data that was used to evaluate the models? I find some of the relative error values strange considering the scatter shown in some of the plots (for example 48-Tuandukou. I wouldn’t use a model with a validation that looked like this). Please elaborate on model validation and evaluation process and data and show the validation plots for each model used (if more than 1 model used).\nIf using more than 1 model, how can we be sure spatial variability in SSC or flux isn’t due to differences in model predictions in adjacent river reaches/points/pixels?\nHow were differences in surface reflectance between landsat 5, 7, and 8 dealt with or tell us why this doesn’t matter for your methods/study?\nFluxes\nThe flux calculation could benefit from more description. For example, what is the temporal resolution of the discharge data (daily?). The equation showed monthly SSC and Q averages were used which sounds good, but there will surely be months (particularly during monsoon) where there are no Landsat SSC observations. How was annual flux computed from sparse data and these biases in data quantity addressed?\nThe sediment budget model using the flux data could also use more description. Were the budget calculations made along a network of “river reaches”? IF so, what is the typically reach length used for the budget model? How can smaller tributaries without SSC or Q measurements be accounted for? What is the temporal scale of flux estimates used in the sediment budget model (long-term mean flux, annual flux, instantaneous flux from each image?)\n-I was often confused when point based SSC/flux data, versus reach-based (if used at all?), versus pixel- based SSC was used in different analysis. Please be clear how data was aggregated over space and time for different analyses.\n-In all analyses, how was temporal bias in Landsat observations dealt with (e.g. there will only be SSC in non-rainy seasons for portion of TP with seasonal cloud cover, ice cover). Was data filtered to sites with representative data first? Was only summer season data used?\n-The phrase “source to sink processes” is used many times and is a good phrase to frame the paper at the beginning, but I recommend getting more specific later on in the paper as this phrase is too general and can mean many things (e.g. are you talking about bed/bank erosion, deposition, flux, burial etc.)\nLines 25-27: Are the % increases in SSC and flux in comparison to “previously reported” studies or over the time period within this study? Please clarify. Also, if referring to previous studies, this is a difficult comparison to make in an abstract considering differences in methods/time periods etc. in previous studies. This result reported in the abstract didn’t seem to be a focus and I couldn’t find it in the manuscript. I would recommend reporting only major results from your study in abstract. What are the main one or two results you really want people to know?\nLine 29-31: “Rapid shifts in the erosion-deposition patterns have regulated the sediment source to sink processes”… This sentence is an example of the type of wording that is frequently used but is too vague to have meaning. While I understand the general idea here, and its good to be general in the abstract, I recommend getting a little more specific as this sentence is difficult to interpret.\nLines 34-36: Consider revising this opening sentence and wording. 1) I assume “water energy” is referring to hydropower? Please clarify. 2) Hydropower is the only item in this list that is provided to people, while ecological function and sediment/carbon/nutrients are provided to ecosystems not people. Please revise this phrasing. 3) Is 2 billion people the population in the TP? If so please clarify and note TP rivers do not supply these things “worldwide” but only in the TP basin, so remove the word worldwide.\nLine 41: What is a sedimentary effect?\nLine 43: Is sediment availability refer to sediment supply? And availability in the landscape or within the river channel? Please be specific when possible.\nLine 46-49: changes in sediment delivery in all environments is the balance of mobilization, transport, and deposition complicated by climate change. Consider getting more specific what this means in the context of “alpine cold environments”.\nLine 50: Note, sediment mobilization and transport are examples of “source-to-sink” processes, consider removing phrase source to sink here and/or be more specific.\nLine 54: What is mean by “geomorphic feedback on these changes”? Consider being more specific or removing this phrase.\nLine 85: This first paragraph seems to be more about spatial heterogeneity or do some of these statistics and sediment yields also refer to temporal trends? Please clarify and perhaps focus this paragraph on spatial patterns if that is your goal. So remove the word “spatiotemporal” if this is only spatial.\nLine 137: “Furthermore, the thawing of permafrost has led to the expansion of thermokarst landscapes and hillslope mass wasting, which resulted in a relative higher SSCs where permafrost is developed on the TP”. Is this an interpretation from the SSC data or relationships you found between thawing permafrost and SSC? Please clarify.\nFigure 1. Are these data long-term average or annual averages as noted in caption?\nLine 167-170: Consider rephrasing this sentence and be more specific. This is an example of a general statement found throughout the manuscript that is difficult to interpret because it’s too general.\nLine 222 and Table S2: Does the underestimation of sediment flux refer to comparison between yields in this study and other studies? If so, where are the sediment yields in Table 2 from other studies? Also if so, I assume other studies were not satellite based but hopefully they are based on surface SSC samples because satellite based SSC and flux will almost always be underestimated compared to a depth-integrated field based measurement. Just making sure your interpretation of underestimated fluxes is not caused by differences in data and satellite biases instead of the deposition processes mentioned.\nSome papers that may be useful background or citations:\nBalasubramanian S V et al 2020 Robust algorithm for estimating total suspended solids (TSS) in inland and nearshore coastal waters Remote Sens. Environ. 246 111768\nGardner et al 2023.Human activities change suspended sediment concentration along rivers. Environ. Res. Lett. 18 064032 https://iopscience.iop.org/article/10.1088/1748-9326/acd8d8\nClaverie M, Vermote E F, Franch B and Masek J G 2015 Evaluation of the Landsat-5 TM and Landsat-7 ETM+ surface reflectance products Remote Sens. Environ. 169 390–403"}, "author_response": "We have thoroughly addressed the referee’s comments in the revised manuscript. A detailed, point-by-point response to each comment can be found within this document.\n\n[... Author responses to Reviewer 1 and Reviewer 2 ...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed this manuscript during the first round. I have now read through the authors response to comments, and I was glad to see their detail response to my earlier comments. I have no further major issues with this study; the manuscript reads well and it will be an important contribution to the geophysical sciences.", "reviewer_2": "The authors have done an great job revising this manuscript. I have no further comments.\nSpecific comments:\nLine 349: I think there is a small typo: \"deep-integrated\" should be \"depth-integrated\", correct?"}, "author_response": ""}]

Decision letter