The Tibetan Plateau (TP) plays a vital role in supplying water to approximately two billion individuals residing downstream, as it serves as the source of numerous significant rivers across Asia. Nevertheless, the significance of groundwater in sustaining the plateau rivers remains uncertain, hence jeopardizing the long-term water resource sustainability of Asia. Here, we assess the groundwater contribution to river flow in the middle reaches of the Yarlungzangbo River (YR), a typical TP location, and we pinpoint the vital role that plateau groundwater plays in sustaining the water balance of the river. The findings indicate that the YR mainstream has a runoff excess of approximately 36–42 × 108 m3/a, which accounts for 22.8–26.5% of the mainstream flow and is strongly correlated with groundwater discharge. The research findings also indicate that the movement of groundwater on the plateau is primarily influenced by the variation in land elevation and replenishment at higher elevations. The contribution of groundwater is inversely proportional to the topographic steepness and precipitation, and its recharge mechanism is associated with the large-scale active fault zones in the vicinity of the YR.

  • Groundwater supports the base flow of the Yarlung Zangbo River (YR), providing 22.8–26.5% annually to the middle reaches of the Yarlungzangbo mainstream.

  • The groundwater flow in the plateau is driven and sustained by the topographic gradient and recharge at high elevations.

  • A series of tensile fault structures in the YR basin are potential channels for groundwater circulation and discharge.

The Tibetan Plateau (TP) serves as the origin of numerous major rivers across Asia and plays a crucial role in providing freshwater to the extensive low-lying regions of China, Bhutan, Nepal, and India (Bookhagen 2012). The elevated terrain of the TP functions as a crucial source of freshwater for a significant portion of the global population. However, climate change is threatening this massive reservoir (Huss & Hock 2018). Due to an increase in rainfall and a decrease in frozen precipitation, the water levels of various lakes, rivers, and underground water sources on the TP and in other mountainous regions have increased significantly since the late 1990s (Langman et al. 2022; Yao et al. 2022; Lorenzi et al. 2024). The groundwater reserves grew at a rate of 5.0 × 109 m3/a from 2003 to 2009, a growth rate that is similar to that of lake water storage (Jiao et al. 2015). Hence, it is of utmost importance to understand the mechanisms through which groundwater originates and its influence on rivers across the TP.

The groundwater system has a crucial function in sustaining river flow during the dry season and storing precipitation throughout the monsoon season (from late May to September), thereby mitigating the peak river flow (Cuthbert et al. 2019). The proportion of groundwater in alpine rivers surpasses that of snow and glacial meltwater (Bookhagen 2012). For instance, groundwater plays a crucial role in sustaining the flow of rivers downstream in the Pamir Mountains, particularly during the winter season when groundwater serves as the primary water source. Groundwater accounts for approximately 40% of the total annual river flow (Pohl et al. 2015). However, the water budget of TP rivers is usually studied in terms of surface water composition, and the consensus is that precipitation and the thawing of ice and snow are the primary contributors to maintaining water equilibrium. The influence of groundwater on rivers is seldom considered, greatly limiting the precision of water resource evaluations within the basin (Li et al. 2020). While the contribution of groundwater discharge to rivers in terms of available water resources is acknowledged, there remains limited understanding regarding the extent to which rivers on the TP rely on groundwater. Recent research has investigated the significance of groundwater systems in mountainous regions on the TP, specifically focusing on faulted areas. Additionally, deep circulating groundwater originating from thermal springs can also contribute to river flow upon emerging at the surface (Andermann et al. 2012). This particular subsurface water movement is typically propelled and sustained by variations in land elevation and replenishment in elevated areas through the infiltration of precipitation and melted snow (Ge et al. 2008; Schmidt et al. 2020). However, there is still a lack of comprehensive knowledge regarding the mechanisms through which groundwater flows into rivers. The quantification of the recharge discharge process and the contribution of groundwater to river discharge remain insufficiently determined. These essential yet absent details are crucial for predicting the future availability of water and alterations in the ecohydrological system of the TP.

The Yarlung Zangbo River (YR) is a representative example of a major river on the TP, and it flows along the collisional boundary between the Indian Plate and the Eurasian Plate (Shen et al. 2022). The most extensive transverse rift zone in the basin is located in the middle reaches of the YR. Significant tectonic activity has the potential to modify the flow pattern of groundwater controlled by topography, leading to an increase in the amount of groundwater contributing to runoff. Consequently, this alteration may impact the composition of rivers. However, previous relevant studies focused mainly on the surface water cycle level, and the abnormal variation in river isotopic compositions was thought to be caused by the spatial zoning of precipitation-forming water vapor and that this was only marginally influenced by the considerable contribution of groundwater to runoff (Yao et al. 2022). Alternatively, other studies have focused on the vertical circulation level of groundwater and have found that groundwater aquifers are recharged by atmospheric precipitation infiltration and discharged by evaporation, concluding that the discharge of shallow groundwater causes changes in the isotopic composition of rivers (Ge et al. 2008). The YR mainstream exhibits an atypical composition during winter despite the dissimilarity in precipitation and water vapor, which can be attributed to more than 90% of the basin's rainfall occurring during the summer season. In alpine basins, the primary process responsible for recharging and discharging groundwater is the vertical and horizontal circulation of water within a specific structure. Therefore, the patterns of groundwater flow under tectonic constraints, as well as its impact on rivers, remain poorly constrained. Our previous research revealed anomalous physical and chemical characteristics in the isotopes, hydrochemistry, and water temperature of the middle reaches of the YR. It was concluded that high-elevation surface water enters the groundwater system through fault zones and is discharged into the YR valley. However, to date, the contribution of groundwater to the runoff of the YR remains undocumented, and the dynamic mechanisms governing groundwater recharge are not yet fully understood.

Therefore, we aimed to examine the hydrological mechanisms of plateau groundwater by focusing on the middle reaches of the YR as a case study. The systematic use of physical and chemical factors, Rn isotopes, and numerical simulation methods reveals the coevolution mechanism of geomorphic structure and groundwater recharge–runoff–discharge processes in alpine valleys, with the goals of (i) identifying the groundwater flow patterns and spatial distribution characteristics of typical alpine basins and (ii) reevaluating the groundwater discharge and its contribution to the runoff of the middle YR, (iii) uncovering the mechanism behind groundwater circulation in typical alpine basins. The originality of this paper lies in the coupling of isotopic and numerical modeling methodologies in order to quantify the crucial role that plateau groundwater plays in sustaining river base flow.

YR, which is long from east to west and narrow from north to south, lies in the middle of the Qinghai-Tibet Plateau, with the world's highest and youngest Himalayas to the south and the Gangdise Mountains and Nianqing Tanggula Mountains to the north (Figure 1(a)). The basin is dominated by an alpine valley landform, and the tectonic unit belongs to the lateral part of the Himalayas. With the ongoing N-S convergence between India and the Asia plate, a series of large-scale thrust and strike-slip faults developed at the boundaries of the terranes (Zhongming et al. 2005). The YR flows from west to east along the fault, cutting a deep canyon in the Himalayas. This research focuses on the section from Lazi to Nugesha in the middle reaches of the YR, spanning approximately 199.48 km, which is the key development area of the ‘one river two tributaries’ project in Tibet. Alpine permafrost and predominantly continuous permafrost are distributed along the southern and northern mountain ridge areas, with seasonally frozen ground occurring along the river valleys in between the mountain ridges (Figure 1(c)).
Figure 1

(a) Geographic location of the TP and YR basin. (b) The YR basin's terrain is depicted with a color scheme that corresponds to elevation using a 1 km digital elevation model. The labels indicate the positions of the sample sites employed in this research. (c) Estimated distribution of various categories of permafrost and seasonally frozen ground in the study area.

Figure 1

(a) Geographic location of the TP and YR basin. (b) The YR basin's terrain is depicted with a color scheme that corresponds to elevation using a 1 km digital elevation model. The labels indicate the positions of the sample sites employed in this research. (c) Estimated distribution of various categories of permafrost and seasonally frozen ground in the study area.

Close modal
Under the action of the tectonic stress field of subduction and pushing to the north of the ancient Indian land, a series of associated and derived N-S-trending structures were formed in addition to the large-scale Yarlung suture zone and other E-W-trending structures. Nearly N-S-trending grabens and rifts developed in en-echelon arrays terminating at strike-slip faults in the center of the plateau and dying out on the southern flank of the Himalayas. Among these fault zones, the Dingjie-Xietongmen-Shenzha (DXS) and Yadong-Dangxiong-Gulu (YDG) fault zones cut the middle valley of the YR at 89°E–92°E, which is the largest transverse rift zone in the basin (Han 1986) (Figure 2(a)). Our study area is the middle reach from Lazi to Nugesha in the YR basin, which includes the southern section of the DXS fault and the YDG fault. Based on the geological map of the basin and the inferred geological sequence permeability, the hydrogeological units of the study area are grouped and simplified as five hydrogeological units (Figure 2(b), Yao et al. 2021): Intrusive rocks unit, Paleogene volcanic unit, Pliocene-Quaternary unit, Upper Cretaceous terrigenous unit, Paleozoic to Mesozoic terrigenous unit (see the Supplementary Material for more details).
Figure 2

Simplified hydrologic map depicting the YR basin.

Figure 2

Simplified hydrologic map depicting the YR basin.

Close modal

According to the formation lithology and aquifer media, the groundwater in the study area is primarily quaternary loose rocks pore water and bedrock fissure water (Liu 2014). The quaternary loose rocks pore water is primarily distributed in the valley and fore-mountain zone in the intermountain region, with the lithology consisting primarily of gravel and gravel layers with loose structures. The depth of the groundwater level is generally <30 m, the aquifer thickness ranges from several meters to hundreds of meters, the hydraulic conductivity is 10−6–10−3 m/s, the specific storage is (4.9–10) × 10−4 m−1, and the specific yield is 0.025–0.35. The bedrock fissure water is primarily found in the bedrock mountain ranges south and north of YR. The aquifer is made up of granite, biotite monzonitic granite, and diorite granite, which are medium acid and ultrabasic intrusion rocks that range from the Yanshanian to the alpine stage. The rock is relatively broken, the hydraulic conductivity is 10−8–10−11 m/s, the specific storage is <3.3 × 10−6 m−1, and the specific yield is 0.001. The rock mass topography is steep, with a thin weathering crust. Furthermore, ground vegetation is sparse, topographic relief is extensive, precipitation is excreted smoothly, and residence duration is brief. Therefore, groundwater recharging takes a short period.

Samples and analyses

In May 2018, field observations and sampling activities were carried out. The electrical conductivity (± 1.0%), pH (± 0.1), and temperature (± 0.1°C) were measured in the field at the time of sampling using calibrated hand-held field meters (HACH: HQ40D). 222Rn was analyzed in the field using a DURRIDGE USA (RAD 7) radon detector. The reported data were averaged from five measurements taken within approximately 1 hour. The detection limit was 0.4000 cpm/pCi/L, and the analytical precision was within ±5% (1σ). The sampling points are shown in Figure 1(b).

Model data

The 90 m Shuttle Radar Topography Mission digital elevation model (DEM) was employed to delineate the terrain features and derive the catchment boundaries for the YR basin. The study area shares the same north‒south boundary as the YR basin, while its western boundary is defined by a river, and its eastern boundary is demarcated by another subbasin (Yangcun Basin). The high topography between these two subbasins serves as a natural hydrological divide (Figure S1). The monthly precipitation and temperature records for Lazi, Rikaze, and Nugesha from 1982 to 2022 were acquired from the China Meteorological Information Center (http://data.cma.cn/). The monthly gridded precipitation product of the TP (0.1° × 0.1°, http://data.tpdc.ac.cn) served as a means of validating the interpolated precipitation of the YR Basin, as the meteorological stations are all located below 4,500 m above sea level (Figure S2). The monthly runoff data from the Rikaze and Nugesha stations for the period of 2013–2018 were acquired from the Hydrology and Water Resources Survey Bureau of Rikaze. Daily groundwater level data for 12 boreholes from 2018 to 2022 were obtained from the China Institute of Geological Environment Monitoring. Data on the distribution of glaciers and permafrost in China were obtained from the second dataset of China's Glacier Inventory.

Groundwater model

We identified five primary hydrogeology–permafrost units in the study area by integrating lithology and permafrost characteristics (Figure 3, see the Supplementary Material for more details). During the period of winter freezing, the ground layer containing seasonal ice is considered to have limited permeability. Nevertheless, this seasonally frozen ground transitions into an aquifer when the summer thaw period occurs (Evans et al. 2015). As the YR basin lies between the Himalayas and the TP, the north‒south boundary can be considered to act as an impermeable boundary. The west is bounded by rivers, which can be regarded as a boundary with fixed water levels. Most of the eastern boundary, except the YR outlet, is the watershed line between adjacent basins, which can be generalized as a confining boundary. The flow boundary in the outlet region is determined by the proportional relationship between the difference between the groundwater level and surface elevation beyond the model area.
Figure 3

Schematic of the conceptual model (Refer to Yao et al. 2021).

Figure 3

Schematic of the conceptual model (Refer to Yao et al. 2021).

Close modal

The numerical groundwater model is constructed based on the Finite Element subsurface FLOW system (FeFlow) software developed by WASY, a German company. The finite element method is employed to discretize the structure, which transforms a complex object into a finite number of easily analyzable units. Additionally, the unknown function within the study area is approximated by employing fragment functions on each unit. The unit discretization process was applied to the model, resulting in the division of the study area into 2,510,236 triangular units. These units were further characterized by a total of 1,272,392 nodes (Figure S3). The vertical dispersion is five layers with a total thickness of approximately 3.5 km. The initial model was developed in a state of equilibrium to establish the fundamental pattern of subsurface water movement within the research region. Recurring models were subsequently created using daily time intervals, building upon the established steady-state model to evaluate the yearly and seasonal proportions of groundwater discharge concerning overall river currents (detailed information on model development is provided in the Supplementary Material).

The model enables the assignment of time series parameters to accurately simulate the changeable climate and hydrological conditions in reality. Therefore, time series parameters are used when parameter assignments are made. In winter (November–April), the seasonally frozen aquifer can be regarded as a weakly permeable layer (Figure 1(c)), with an average hydraulic conductivity (AHC) of 8.2 × 10−10 m/s. The main component of the permafrost unit was envisioned as a consistently weakly permeable layer with an AHC of 5.2 × 10−11 m/s (Yao et al. 2021). The alpine permafrost layer was envisioned as a somewhat less cohesive, weakly permeable layer than the northern uninterrupted permafrost layer, possessing an AHC of 6.5 × 10−11 m/s (Figure S4). In summer (May–October), the seasonally frozen layer melts, forming an aquifer with an AHC of 3.0 × 10−4 m/s, while the thawed areas outside the river valleys exhibit a lower AHC of 7.7 × 10−6 m/s. The ‘Discrete Feature’ module was used to generalize the fracture zone with an AHC of 10−3 m/s. Darcy's law was applied to define the fluid movement within the fault. Groundwater recharge was determined as the product of precipitation and the rainfall infiltration coefficient (α) for different months. The precipitation infiltration coefficient is based on the empirical values of different land distribution types. The α values are 0.4–0.8 in the river valley area, 0.1–0.2 in most areas of the mountain area, and 0.2–0.4 in the alpine source area, which is characterized by extensive wind compounds and moraine-covered regions with significant permeability. The initial values of recharge and hydraulic conductivity in the study area were obtained from Ge et al. (2008), Chen et al. (2021), and Yao et al. (2021) and adjusted by manual calibration and the PEST inversion algorithm (Doherty & Hunt 2010). Relevant data from the study area, including winter base flow records (from November to April) and borehole groundwater level monitoring data, were utilized as calibration targets for the transient model.

Model calibration

River flow primarily originates from precipitation, glacial ice melting, and the discharge of groundwater (Andermann et al. 2012). However, the majority of river flow is attributed to glacial meltwater during the summer (May–October), while its contribution decreases significantly in the dry winter period (November–April). As glaciers are situated in both the southern and northern parts of the basin and are located in the upper source region, the significance of glacial meltwater is limited to this specific area. Furthermore, the amount of winter precipitation primarily consists of snowfall and is less than 10 mm per month. Therefore, the recharge of river flow by surface runoff derived from precipitation and meltwater is assumed to be negligible in the dry winter months. During the winter season, the main contributor to river flow in the middle reaches of the YR basin is groundwater (Hayashi 2020). Stream baseflow, derived from groundwater discharge, can be most accurately estimated by analyzing the observed streamflow data during this period (Walvoord et al. 2012). Therefore, the model parameters are adjusted and calibrated based on the base flow in the Nugesha and Rikaze subbasins in winter and the groundwater level monitoring data at each borehole point.

To assess the effectiveness of the groundwater flow models, evaluation criteria such as relative error (RE) and Fit-NRMSE were employed (Equations (3) and (4) in the Supplementary Material). Figure 4 shows the simulated and observed groundwater levels of the model (2018–2023 daily data and groundwater level data from the China Institute of Geological Environment Monitoring), with RE values ranging from −0.4% to 0.19%, with an average of approximately −0.13%. The mean value of the Fit-NRMSE coefficient was 0.85 (Table S1). Figure 5(a) is a fitting diagram of the observed groundwater level and the simulated value in a civil well during the field survey in 2018. The RE values range from −1.28% to −0.46%, with an average of −0.78%. In general, the groundwater level is well-fitted, and the simulated groundwater level is slightly lower than the actual water level. Figure 5(b) shows a comparison of the simulated base flow and the observed base flow (2013–2018) in both subbasins within the model. The fitted RE value of the Nugesha subbasin is 3.08%, and the Fit-NRMSE coefficient is 0.79; the fitted RE value of the Rikaze subbasin is 3.96%, and the Fit-NRMSE coefficient is 0.78. The simulated baseflow exhibits a high level of concurrence with the recorded data. Hence, the model offers plausible and dependable outcomes, as its approximated baseflows serve as a reliable indicator for real groundwater discharge into streams.
Figure 4

Simulated and observed (during 2018–2022) borehole water levels for the steady-state model (the location of drilling points is shown in Figure S5).

Figure 4

Simulated and observed (during 2018–2022) borehole water levels for the steady-state model (the location of drilling points is shown in Figure S5).

Close modal
Figure 5

Contrast between simulated and observed data. (a) Water levels in pumping wells within the study area were compared between computed and observed values under steady-state conditions (the position of the observation point is shown in Figure S5). (b) Baseflow values were determined and compared between the two hydrological stations in relation to the steady-state model. The blue triangles represent the baseflow calculated by the dynamic model, while the boxplot displays the observed streamflow during periods of low winter precipitation.

Figure 5

Contrast between simulated and observed data. (a) Water levels in pumping wells within the study area were compared between computed and observed values under steady-state conditions (the position of the observation point is shown in Figure S5). (b) Baseflow values were determined and compared between the two hydrological stations in relation to the steady-state model. The blue triangles represent the baseflow calculated by the dynamic model, while the boxplot displays the observed streamflow during periods of low winter precipitation.

Close modal

Groundwater flow patterns

The representation of groundwater head patterns and flow routes in the simulation reveals interconnected groundwater flow systems, including a smaller-scale system characterized by shorter distances, a larger catchment-scale system, and an even broader regional-scale system with longer flow paths. The distribution of the water head is greatly influenced by the topography and hydrogeological context (Figure 6(a)). The simulated water pressure gradually decreases as it moves from the southern Himalayas and northern Tanggula Mountains toward the central YR valley area, aligning perfectly with the surrounding topographical features. The groundwater head in the northern and southern mountains exceeds 4,000 m, covering an area amounting to approximately 90.7% of the YR middle reaches. The eastern exit of the middle reaches of the YR features the minimum water head, with a value of 3,769 m. The Pathline module of FeFlow identifies 106,732 groundwater paths terminating by discharge into rivers and provides particle tracking distances. Short-distance (<10 km) flow paths (over 70% of the total path lines) are mainly distributed near the valleys of the YR and its tributaries, representing the topography of general underground flow patterns and indicating that the groundwater in these areas has a shorter renewal cycle and younger groundwater age. The permafrost regions in the north and south also show thinner and shorter paths. The moderate-length (10–20 km) flow paths account for approximately 23% of the total path lines and are mainly distributed in the southern mountainous area with elevations of 4,000–5,000 m. The distribution of flow paths exceeding 20 km in length is relatively limited, comprising merely 7% of the overall path lines. These longer paths are primarily concentrated in the northern region at elevations surpassing 5,000 m. Revealing the presence of a theoretical nested groundwater cycle in an actual watershed, simulated flow paths are depicted across local and regional scales (Gleeson & Manning 2008). This depiction further confirms that the flows within the catchment area have a great influence on the regional groundwater systems and interactions between surface water and groundwater (Jasechko et al. 2016; Yao et al. 2021).
Figure 6

(a) Geographical distribution of groundwater flow. (b) The simulated average annual baseflow discharge (Qg) per kilometer along streams and groundwater recharge (R) per unit basin area.

Figure 6

(a) Geographical distribution of groundwater flow. (b) The simulated average annual baseflow discharge (Qg) per kilometer along streams and groundwater recharge (R) per unit basin area.

Close modal

Groundwater recharge is the process by which an underground water-bearing system receives water from the surrounding environment. This term refers to the recharging of underground water via precipitation, meltwater, irrigation, underground runoff, or river seepage. Groundwater discharge refers to the release of water from subterranean aquifers into the surrounding ecosystem, and springs serve as natural locations where this discharge occurs. Hot springs are widespread in the middle reaches of the YR (Figure S1), and groundwater discharge varies greatly from less than 0.1 to more than 30 × 106 m3/km/a (expressed in unit river length). The maximum discharge is mainly concentrated in the river valley where the north–south-trending fault zones meet the YR, reflecting the water-transport effect of the north–south-trending tensile fault zones. The study area experienced an increase in replenishment, with the highest rates observed in the central YR valley and lower rates observed toward the south and north, from 0 to 450 mm/a (expressed in unit basin area). In the middle of the YR, in 82.3% of the area, the groundwater discharge is less than 10 × 106 m3/km/a, and in 78.2% of the area, the groundwater system recharge is less than 200 mm/a (Figure 6((b)). Areas with high groundwater discharge to rivers (>10 × 106 m3/km/a) are concentrated in regions characterized by elevated recharge rates (300–500 mm/a). These areas are primarily situated within the alluvial sector at the confluence of the Rikaze and Nugesha subbasins. The significance of this location stems from its role as a central hub for population and agriculture in the middle reaches of the YR. Consequently, irrigation practices and reservoir management contribute significantly to augmenting both recharge and discharge rates. Furthermore, based on the simulation results, the YR mainstream area from the Lazi to Nugesha stations has a groundwater discharge of 15–20 × 106 m3/km/a. The river reach is approximately 200 km long, meaning that approximately 30–40 × 108 m3 of groundwater is discharged annually to sustain the base flow of the YR. This result confirms the significant contribution of groundwater to the flow of the middle YR, as it is similar to the groundwater discharge calculated in our earlier study (36 × 108 m3). Further discussion regarding the water source and the circulation mechanism of such a large discharge is still required.

Implication of physicochemical parameters and FBM and 222Rn MBM

The chemical and physical parameters and radon isotopic compositions of water bodies are useful for revealing groundwater migration patterns and quantifying groundwater discharge. Water temperature is a direct indicator of the hydraulic interaction between surface water, shallow groundwater, and geothermal water. The temperature fluctuations in the river water and groundwater are significant in the study area. The temperature of the surface water (7–19.5°C) was close to that of the shallow groundwater (8.4–18.5°C) but significantly lower than that of the geothermal water (10.9–80+°C) (Figure 7(a) and Table S2). The YR temperature fluctuated and increased in the runoff direction, with typical regional characteristics. The temperature of the YR mainstream increases significantly after passing through the Gucuo fault zone, and then the low-temperature Nianchu River near 89°E enters the mainstream, causing an apparent decrease (10.1 °C) in the mainstream temperature. The river then passes through the YDG fault zone, where the mainstream temperature reaches an obvious peak (19.5 °C). The B-spline curve clearly displays the obvious temperature increase in the YR mainstream after passing through the fault zone, indicating that the river receives high-temperature groundwater discharge near the fault zone. The pH of the YR river water and groundwater ranges from 7.6 to 9.2, indicating neutral to weakly alkaline conditions. The pH value of the mainstream displays distinct spatial characteristics and is much lower in the fault zone (Figure 7(b)). In general, a decrease in pH can result from deep degassing of groundwater in a fault zone. The concentration of volatile dissolved components decreases, and CO2, H2S, and other gases separate from the groundwater due to temperature and pressure changes during the upwelling process, the differentiation and evolution of silicate minerals in the water, and the processes of cooling and precipitation. Furthermore, the confluences with the Nianchu River (89°E) and Lhasa River (91°E) are associated with changes in pH, and the pH values of the two rivers as they join the YR are greater than those of the YR mainstream. The lower pH of the river water at this point is not affected by the input of the higher-pH tributaries, indicating that the tributaries have minimal effects on the YR mainstream and are most likely the result of groundwater input into the Nianchu and Lhasa Rivers.
Figure 7

(a) Spatial variation in water temperature in the YR middle reaches; (b) spatial variation in pH in the middle reaches of the YR.

Figure 7

(a) Spatial variation in water temperature in the YR middle reaches; (b) spatial variation in pH in the middle reaches of the YR.

Close modal
The YR exhibits clear regional patterns in terms of drainage and the presence of 222Rn isotopes. From the Lazi station to the Nugesha station, the runoff increased significantly. The average annual YR runoff at the Nugesha station in the eastern Rikaze subbasin (500.27 m3/s) was approximately three times that at the Lazi station in the west (151.46 m3/s) (Figure 8(a)). Although there are tributaries in between, none of them, except the Nianchu River, are large enough to influence the flow of the YR mainstream. Furthermore, the annual average precipitation in Lazi (370.7 mm) and Nugesha (435.6 mm) is similar, and the temperature in the east of Rikaze is higher, with significantly greater evaporation than that in the west in Lazi County. It is conservatively estimated that if no other water source exists, the total runoff at the Nugesha station cannot increase significantly by relying only on the inherited inflow of the Lazi and Nianchu River tributaries and the precipitation in this section. As a result, additional sources (such as groundwater) must be introduced to balance runoff at the Nugesha station. According to the flow balance model (FBM, Equation (7)), the exchange flux and recharge mode between groundwater and river water in the middle YR (Lazi–Nugesha) can be further examined to determine the amount of groundwater discharged to the YR valley. The discharge, precipitation, and water collection areas used in this model are derived from the hydrographic station records in Lazi and Nugesha. The data and results obtained from the calculations are presented in Table 1. The results show that an average of 36.0 × 108 m3 of groundwater is discharged into the river every year between the Lazi station and the Nugesha station, accounting for 22.8% of the total annual runoff at the Nugesha station (157.76 × 108 m3). In terms of the radon isotopes in the river water, the 222Rn concentration at the Nugesha station is approximately 150 Bq/m3 higher than that at the Lazi station (Figure 8(b)). The considerable increase in 222Rn concentrations in river water from the Lazi station to the Nugesha station is also due to the discharge of groundwater with high 222Rn concentrations into the river. The 222Rn isotope mass balance model (MBM, Equation (12)) is calculated using a field survey and observed data on runoff velocity, river depth, and 222Rn content. Because no 222Rn concentration data for groundwater were obtained in the field, the average concentration (Cg) of groundwater 222Rn in Equation (12) is estimated using the average 222Rn concentration of geothermal water recorded near the YR fault zone (DXS) (Table S2). The groundwater discharge from Lazi to Nugesha is 41.9 × 108 m3, amounting to 26.5% of the yearly runoff at Nugesha station. In this model, geothermal water is used to represent the 222Rn concentration in groundwater. As the water temperature increases, radioactive elements such as U and Th in rocks migrate into the water body, resulting in 222Rn concentrations that may be higher than those observed in low-temperature groundwater. If relatively low concentrations of low-temperature groundwater are considered, the annual groundwater recharge may be higher. Overall, the amount of groundwater discharge in the YR mainstream, calculated using the FBM and 222Rn MBM, is approximately 36.0–41.9 × 108 m3. Because the influences of evaporation, the groundwater runoff coefficient, and various secondary tributaries are not considered in the calculation, there may be errors in the results. Nevertheless, this result closely approximates the groundwater discharge (30–40 × 108 m3) previously calculated using a numerical model. This finding further underscores the significance of these data as a reference and highlights the crucial role played by groundwater in sustaining base flow in the middle reaches of the YR.
Table 1

The calculation parameters and results of the FBM and 222Rn MBM

FBMLongitude/(°E)Latitude/(°N)Runoff/Q (m3/s)Precipitation/P0 (mm)catchment area/S (×104 km2)Groundwater discharge/Vg (×108 m3)
LZ 87.68 29.16 151.46 370.7 6.720 35.982 
NGS 89.78 29.38 500.27 435.6 13.52  
MBM222Rn concentration/C (Bq/m3)Flow velocity/v (m/s)Depth/h (m)Length of reach/l (m)groundwater discharge per unit length/q (×10−4 m3/s/m)Groundwater discharge/Vg (×108 m3)
LZ 485.07 4.5 340.19 3.902 41.865 
NGS 708.83      
FBMLongitude/(°E)Latitude/(°N)Runoff/Q (m3/s)Precipitation/P0 (mm)catchment area/S (×104 km2)Groundwater discharge/Vg (×108 m3)
LZ 87.68 29.16 151.46 370.7 6.720 35.982 
NGS 89.78 29.38 500.27 435.6 13.52  
MBM222Rn concentration/C (Bq/m3)Flow velocity/v (m/s)Depth/h (m)Length of reach/l (m)groundwater discharge per unit length/q (×10−4 m3/s/m)Groundwater discharge/Vg (×108 m3)
LZ 485.07 4.5 340.19 3.902 41.865 
NGS 708.83      

Note. LN refers to Lazi station; NGS refers to Nugesha station.

Figure 8

(a) Spatial variation in runoff in the YR; (b) spatial variation in Rn isotopes in the middle reaches of the YR.

Figure 8

(a) Spatial variation in runoff in the YR; (b) spatial variation in Rn isotopes in the middle reaches of the YR.

Close modal

Control of groundwater flow pattern and discharge

We conducted additional analyses to investigate how topography and recharge gradients influence the groundwater flow distribution and discharge in the central YR using the simulated results. Topographic roughness is a metric that is employed to evaluate the degree of steepness of terrain. As the topographic roughness decreases, the terrain becomes flatter (Gleeson & Manning 2008; Yao et al. 2021). The correlation between topography and the length of groundwater flow paths is significant. In the case of a similar elevation difference, the length of the flow path increases with flatter terrain characterized by lower topographic roughness, while it decreases with steeper terrain characterized by higher topographic roughness (Figure 9). Our findings directly support the topographic control of groundwater flow paths. In addition, the groundwater flow pattern depicted in Figure 10(a) indicates that the groundwater level exceeds the river stage, providing further support for the discharge of groundwater into the middle reaches of the YR. The relationship between the distribution of the groundwater flow field and DEM elevation (Figure S6) provides evidence that the topographic gradient and high-elevation recharge sources play significant roles in generating and sustaining groundwater flow on the plateau.
Figure 9

Relationships between the maximum variations in elevation along flow paths and their corresponding distances.

Figure 9

Relationships between the maximum variations in elevation along flow paths and their corresponding distances.

Close modal
Figure 10

A 3D model of the YR where the color palette represents elevation. (a) Groundwater level in the middle reaches of the research site; (b) conceptual model of groundwater circulation in the study area.

Figure 10

A 3D model of the YR where the color palette represents elevation. (a) Groundwater level in the middle reaches of the research site; (b) conceptual model of groundwater circulation in the study area.

Close modal

In addition to topographic control, basin fault structure has a considerable impact on groundwater flow. The groundwater discharge in the middle reaches of the YR, particularly in the valley area, is greater than that of large rivers on plains. Over a river length of less than 200 km, approximately 36–42 × 108 m3 of groundwater is discharged annually to maintain the base flow of the YR. If such a large flow is realized only by the slow infiltration of precipitation along the unsaturated zone in traditional hydrology, the whole cycle would take a long time, which contradicts the groundwater age and the weak water‒rock interaction of the groundwater observed in a previous study (Tan et al. 2014). Therefore, the groundwater replenishment of the YR mainstream may include not only water moving slowly through the unsaturated zone or surface fractures but also concentrated groundwater drainage with rapid circulation of alpine meltwater and precipitation along deep and large fractures in the middle reaches of the YR. Recent research suggests that deep groundwater may represent the world's largest ‘reservoir’. Because the majority of this groundwater is deep and is stored in rocks with low permeability, it has difficulty circulating or flowing to the surface. If the high-elevation recharge areas are proximal to the low-elevation discharge areas and a significant fault is present, substantial pressure differentials will occur. Consequently, this pressure disparity will facilitate the movement of shallow water toward deeper regions, where it can actively participate in deep circulation before eventually resurfacing through water conduction within the fault zone (McIntosh & Ferguson 2021; Shi et al. 2021). Extensive deep faults with nearly north‒south orientations were observed in the study region, particularly in the southern section of the basin (Figures 2 and S1). In terms of mechanical characteristics, these fault structures are mostly defined by tension or tensional shear, which considerably aids the circulation of precipitation and meltwater from the southern Himalayas to the YR valley via fault zones (Figure 10(b)). The distance of lateral migration of high-elevation recharge water along a fault may exceed 150 km (Su et al. 2020), and the water flow in each extensional fault will be driven by the large difference in the water head generated by the topography and by the density difference between hot and cold water. The recharge water will be discharged as spring water at the fault's lowest point where the fault meets the YR.

An FeFlow-based groundwater flow model was developed to investigate the dynamics of groundwater flow in the Lazi–Nugesha reaches of the YR basin. The model considers a variety of elements, including precipitation, fault structure, and permafrost, and combines river water physicochemical properties and Rn isotopes to identify groundwater discharge mechanisms. The movement of groundwater in the examined region is related to variations in topography, in which a greater disparity in elevation leads to a reduced distance traveled by the flow. In the middle reaches of the YR, short (<10 km) flow paths account for more than 70% of the spatial flow pattern of the groundwater system, indicating rapid groundwater renewal. The recharge of groundwater occurs mainly through precipitation or meltwater at high elevations, followed by discharge in valleys and faults. Furthermore, data on water temperature, pH, and Rn isotopes indicate the driving effect of water conduction in fault zones on groundwater recharge and discharge. The annual discharge of the groundwater that reaches the mainstream of the YR along fractures is approximately 36–42 × 108 m3, accounting for 22.8–26.5% of the discharge at the Nugesha station. The significance of groundwater in sustaining the continuous flow of alpine rivers is underscored, facilitating investigations into regional hydrological phenomena. However, due to the complexity of actual samples and the impact of simulation parameters, the estimation of groundwater's contribution to YR flow may be biased. Therefore, additional research into groundwater migration patterns is required in the future to provide more accurate groundwater resource assessments and planning.

D. S. conceptualized the process, designed the work, contributed in the data collection, and wrote the original draft, developed the methodology, and provided the software; S. H. wrote the review and edited the article. All authors have read and agreed to the published version of the manuscript.

This study was financially supported by the School Scientific Research Development Fund project (2022LFR091).

We appreciate Yu Zhang, and Peixin Cong for their assistance with the field work. We also thank the editors for their valuable comments.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Andermann
C.
,
Longuevergne
L.
,
Bonnet
S.
,
Crave
A.
,
Davy
P.
&
Gloaguen
R.
(
2012
)
Impact of transient groundwater storage on the discharge of Himalayan rivers
,
Nature Geoscience
,
5
,
127
132
.
Bookhagen
B.
(
2012
)
Himalayan groundwater
,
Nature Geoscience
,
5
,
97
98
.
Chen
J.
,
Kuang
X.
,
Lancia
M.
,
Yao
Y.
&
Zheng
C.
(
2021
)
Analysis of the groundwater flow system in a high-altitude headwater region under rapid climate warming: Lhasa River Basin, Tibetan Plateau. Journal of Hydrology: Regional Studies, 36, 100871
.
Cuthbert
M. O.
,
Gleeson
T.
,
Moosdorf
N.
,
Befus
K. M.
,
Schneider
A.
,
Hartmann
J.
&
Lehner
B.
(
2019
)
Global patterns and dynamics of climate–groundwater interactions
,
Nature Climate Change
,
9
,
137
141
.
Doherty
J. E.
&
Hunt
R. J.
(
2010
)
Approaches to Highly Parameterized Inversion: A Guide to Using PEST for Groundwater-Model Calibration
.
Wisconsin: US Geological Survey
, pp.
2328
0328
.
Evans
S. G.
,
Ge
S.
&
Liang
S.
(
2015
)
Analysis of groundwater flow in mountainous, headwater catchments with permafrost
,
Water Resources Research
,
51
,
9564
9576
.
Ge
S.
,
Wu
Q. B.
,
Lu
N.
,
Jiang
G. L.
&
Ball
L.
(
2008
)
Groundwater in the Tibet Plateau, western China
,
Geophysical Research
,
35
,
L18403
.
Han
T.
(
1986
)
Active tectonic belts in Tibet and their control on earthquakes
,
Acta Geologica Sinica
,
60
(
4
),
23
29
.
Huss
M.
&
Hock
R.
(
2018
)
Global-scale hydrological response to future glacier mass loss
,
Nature Climate Change
,
8
,
135
140
.
Jasechko
S.
,
Kirchner
J. W.
,
Welker
J. M.
&
McDonnell
J. J.
(
2016
)
A substantial proportion of global streamflow is less than three months old
,
Nature Geoscience
,
9
(
2
),
126
129
.
Langman
J. B.
,
Martin
J.
,
Gaddy
E.
,
Boll
J.
&
Behrens
D.
(
2022
)
Snowpack aging, water isotope evolution, and runoff isotope signals, Palouse Range, Idaho, USA
,
Hydrology
,
9
(
6
),
94
.
Li
Z.
,
Li
Z.
,
Feng
Q.
,
Zhang
B.
,
Gui
J.
,
Xue
J.
&
Gao
W.
(
2020
)
Runoff dominated by supra-permafrost water in the source region of the Yangtze River using environmental isotopes
,
Journal of Hydrology
,
582
,
124506
.
Liu
Z.
(
2014
)
Study on the formation mechanism of the typical tropical high-temperature geothermal system in Nimu-Nagqu, Tibet
.
Doctoral dissertation. Beijing: Chinese Academy of Geological Sciences
.
McIntosh
J. C.
&
Ferguson
G.
(
2021
)
Deep meteoric water circulation in earth's crust
,
Geophysical Research Letters
,
48
(
5
),
e2020GL090461
.
Schmidt
A. H.
,
Lüdtke
S.
&
Andermann
C.
(
2020
)
Multiple measures of monsoon-controlled water storage in Asia
,
Earth and Planetary Science Letters
,
546
,
116415
.
Shen
T.
,
Wang
G.
,
van der Beek
P.
,
Bernet
M.
,
Chen
Y.
,
Zhang
P.
,
Wang
A.
&
Cao
K.
(
2022
)
Impacts of late Miocene normal faulting on Yarlung Tsangpo River evolution, southeastern Tibet
,
Bulletin
,
134
,
3142
3154
.
Shi
D.
,
Tan
H.
,
Chen
X.
,
Rao
W.
&
Basang
R.
(
2021
)
Uncovering the mechanisms of seasonal river–groundwater circulation using isotopes and water chemistry in the middle reaches of the Yarlungzangbo River, Tibet
,
Journal of Hydrology
,
603
,
127010
.
Tan
H.
,
Zhang
Y.
,
Zhang
W.
,
Kong
N.
,
Zhang
Q.
&
Huang
J.
(
2014
)
Understanding the circulation of geothermal waters in the Tibetan Plateau using oxygen and hydrogen stable isotopes
,
Applied Geochemistry
,
51
,
23
32
.
Yao
T.
,
Bolch
T.
,
Chen
D.
,
Gao
J.
,
Immerzeel
W.
,
Piao
S.
&
Zhao
P.
(
2022
)
The imbalance of the Asian water tower
,
Nature Reviews Earth & Environment
,
3
(
10
),
618
632
.
Yao
Y.
,
Zheng
C.
,
Andrews
C. B.
,
Scanlon
B. R.
,
Kuang
X.
,
Zeng
Z.
,
Jeong
S. J.
,
Lancia
M.
,
Wu
Y.
&
Li
G.
(
2021
)
Role of groundwater in sustaining northern Himalayan rivers
,
Geophysical Research Letters
,
48
(
10
),
e2020GL092354
.
Zhongming
Z.
,
Linong
L.
,
Xiaona
Y.
,
Wangqiang
Z.
&
Wei
L.
(
2005
)
Nyainqentanglha Shan: A window into the tectonic, thermal, and geochemical evolution of the Lhasa block, southern Tibet, Journal of Geophysical Research: Solid Earth, 110 (B8), B08413
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Supplementary data