Role of groundwater discharge in sustaining baseflow balance in the middle reaches of the Yarlung Zangbo River

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 × 10⁹ m³/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.

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⁻⁶–10⁻³ m/s, the specific storage is (4.9–10) × 10⁻⁴ m⁻¹, 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⁻⁸–10⁻¹¹ m/s, the specific storage is <3.3 × 10⁻⁶ m⁻¹, 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.

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). ²²²Rn 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).

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.

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⁻¹⁰ m/s. The main component of the permafrost unit was envisioned as a consistently weakly permeable layer with an AHC of 5.2 × 10⁻¹¹ 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⁻¹¹ m/s (Figure S4). In summer (May–October), the seasonally frozen layer melts, forming an aquifer with an AHC of 3.0 × 10⁻⁴ m/s, while the thawed areas outside the river valleys exhibit a lower AHC of 7.7 × 10⁻⁶ m/s. The ‘Discrete Feature’ module was used to generalize the fracture zone with an AHC of 10⁻³ 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.

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.

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 × 10⁶ m³/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 × 10⁶ m³/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 × 10⁶ m³/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 × 10⁶ m³/km/a. The river reach is approximately 200 km long, meaning that approximately 30–40 × 10⁸ m³ 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 × 10⁸ m³). Further discussion regarding the water source and the circulation mechanism of such a large discharge is still required.

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 × 10⁸ m³ 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 × 10⁸ m³, 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.