GoKawiilStudy area and field sampling
Spanning about 30% of the entire QTP, this 780,000 km2 study area is situated between the boundary of the Loess Plateau (Gansu Province, 103° E) and the border of China (Tibet Autonomous Region, 78° E), and includes the headwater of eight major river basins: the Yellow, Yangtze, Lancang–Mekong, Nu–Salween, Derung (Irrawaddy), Yarlung Tsangpo–Brahmaputra, Marja Tsangpo (Ganges) and Indus Rivers (Fig. 1a and Supplementary Table 1). This region spans gradients of climate, topography, vegetation cover, permafrost extent and lithology, with elevations between 1,600 and 8,000 m above sea level. The rivers drain the major tectonic terranes of the QTP that contain a wide range of sedimentary, igneous and metamorphic rocks51,52 (Extended Data Fig. 1b). In the northeast of the study area, the Yellow and Yangtze Rivers drain marine carbonate and siliciclastic sediments of the eastern Kunlun–Qaidam Terrane and the Triassic deep-marine sediments of the Songpan–Ganzi flysch complex51. The Lancang, Nu and the southwestern Yangtze catchments drain mostly terrestrial siliciclastic and shallow marine carbonate rocks of the eastern Qiangtang and northern Lhasa Terrane51. In particular, the sediments of the Qiangtang Terrane contain abundant evaporites45. By contrast, the catchments of the Yarlung Tsangpo and Indus Rivers are underlain by a series of plutonic and volcanic rocks of the southern Lhasa Terrane apart from the marine Tethyan sedimentary series of the Himalayas (Extended Data Fig. 1b). Permafrost extent across the study area varies from continuous to isolated with the most extensive permafrost extent concentrated in the Yangtze and eastern Yellow Rivers, and high-elevation areas of all other river basins53. The climate is characterized by a cold and dry winter season and an ice-free season between April and October with heavy monsoon rains during the study period54. Potential glacial influence is minimized by sampling streams at least 20 km downstream from glaciers. Flow distance and tributaries collectively reduce potential biogeochemical signatures of glacial meltwaters.
We collected 175 individual samples from 50 rivers across this area during the daytime in spring (May–June), summer (July–August) and autumn (September–October) between 2016 and 2018, and in spring ice-out (early April) in 2023. The Yellow River basin was visited eight times; the Yangtze River basin five times; and the Lancang, Nu and Yarlung River basins were all sampled on four dates, the Derung and Indus River basins twice, and the Marja River basin once. We measured, compiled, interpolated and/or calculated all variables first for each of the 175 samples (Supplementary Table 2) and then calculated mean values and their standard deviations for sites with multiple samples to obtain one estimate for each of the 50 sampling sites (Supplementary Table 5).
Topography, hydrology, lithology and climate parameters
All topographic analyses were performed on a 90-m resolution digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM)55 using version pre2.5 of the TopoToolbox56. We extracted a river network and estimated, for each sample location, drainage-basin-averaged slopes, S bsn , and normalized steepness index, k sn , with a reference concavity of 0.45, a segment length of 1,000 m and a drainage-area threshold of 100 km2. For all rivers in the network, we also extracted the stream order (with drainage-area threshold of 5 km2) and calculated the total length of rivers, \({L}_{{\mathrm{riv}}_{i}}\) with a given order, that is, upstream of each sampling point. All these analyses were performed on a DEM projected in ArcGIS-Pro using the inbuilt Asia North Lambert Conformal Conic projection. To provide an estimate of the total surface area of rivers upstream of a sampling point, A riv , we multiplied the river length in each stream order, \({L}_{{\mathrm{riv}}_{i}}\), with an average width for rivers in that stream order, \({\bar{W}}_{i}\):
$${A}_{\mathrm{riv}}=\mathop{\sum }\limits_{1}^{N}{L}_{{\mathrm{riv}}_{i}}{\bar{W}}_{i}$$ (1)
Where N is the highest stream order in each catchment. For stream orders 3–8, \({\bar{W}}_{i}\) and its uncertainty were based on the mean and standard deviation of 100 individual width measurements (only 30 measurements for stream order 8). Measurements from hydrologic stations at our sampling sites were combined with widths of randomly selected rivers in the study area mapped on Google Earth (Supplementary Table 6). Average width for the first- and second-order streams were estimated by extrapolating an exponential fit to rivers of orders 3–6, with an uncertainty based on the prediction band at these stream orders (Supplementary Fig. 1). We used only stream orders 3–6 for the extrapolation, because many higher-order rivers in the study area are inset into narrow bedrock gorges and break the exponential increase in stream-widths18 (Supplementary Fig. 1). Our interpolation predicted the first-order stream widths of 4.3 (1.3–13.8) m, which covers the range reported in ref. 57 (1.9 ± 1.1 m), but is wider than the value reported more recently in ref. 58 (0.32 ± 0.07 m), leading to a potential overestimate in CO 2 emissions.
The same DEM was analysed in TopoToolbox after projecting it with the Asia North Albers Equal Area Conic projection in ArcGIS-Pro to extract drainage areas, A bsn , and to estimate an areal fraction of permafrost extent in the landscape upstream of each sampling point. We based the latter analysis on a published map of permafrost extent53 that classifies surfaces as being underlain by continuous (≥90%), discontinuous (≥50% and <90%), sporadic (≥10% and <90%) or isolated (<10%) permafrost. Because this classification is based on ranges, we assigned values of 95%, 70%, 30% and 5%, respectively, to each permafrost class (0% for no permafrost), and then calculated the average pixel value upstream of each sampling point in TopoToolbox. Finally, we re-classified all catchment-averaged permafrost-extent values into the same permafrost categories. We generated a TIFF image from the global lithologic map59 to extract the fraction of carbonate lithologies upstream of each sampling point in TopoToolbox. In the fraction of carbonate rocks, we included the lithologic categories ‘sc—carbonate sedimentary rocks’ and ‘sm—mixed sedimentary rocks’. We did not consider ‘mt—metamorphics’ here, but we found that including these would not make a substantial difference to our analysis. We estimated the fraction of wetlands, the mean annual air temperature and the mean annual precipitation values for the watersheds draining each of our sites by matching the coordinates of each site with the corresponding watershed in HydroATLAS60, in which we obtained the relevant attributes.
For sampling sites within the cross-section of the hydrometric station, we obtained a measurement of discharge, Q spl (m3 s−1), and suspended sediment concentrations (Supplementary Table 2) at the time of sampling. We were also able to obtain average annual discharge, Q annual (m3 yr−1), for the most downstream sampling sites of all sampled catchments (basin IDs 9, 16, 18, 19, 22, 27, 31, 37, 45, 46, 47, 48, 49 and 50) and annual suspended sediment for the four largest outlet basins (basin IDs 7, 16, 22 and 27; 74% of the total study area). Note that suspended sediment concentrations were not available for the sampling site with ID9, but for the more upstream sampling site ID7.
Sample collection and analyses
... continue reading