
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12242-6
10.1016/j.heliyon.2024.e36211
e36211
Research Article
Sources and pathways of spring flow and climate change effects in the Dungju Ri & Yude Ri catchments, Bhutan Himalaya
Dendup Tshewang tshewang.sherubtse@rub.edu.bt
a⁎
Tshering Dendup dendup.sherubtse@rub.edu.bt
b
Tobgay Sonam tobgay88@gmail.com
c
Liu Fengjing fliu7@mtu.edu
d
a Department of Physical Science, Sherubtse College, Royal University of Bhutan, Kanglung, 42002, Trashigang, Bhutan
b Department of Environment and Life Science, Sherubtse College, Royal University of Bhutan, Kanglung, 42002, Trashigang, Bhutan
c Divisional Forest Office, Department of Forest & Park Services, Trashigang, 42001, Bhutan
d College of Forest Resources and Environmental Science, Michigan Technological University, Houghton, MI, 49931, USA
⁎ Corresponding author. tshewang.sherubtse@rub.edu.bt
14 8 2024
30 8 2024
14 8 2024
10 16 e3621123 3 2024
11 8 2024
12 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Springs and streams are vital water sources for supporting the livelihood of Himalayan residents. Escalating climate change, population growth, and economic development strain the region's freshwater resources. A national survey reveals declining spring and stream flows in Bhutan, necessitating an improved understanding of their generation. Monthly grab water samples were collected during April 2022–January 2023 from main streams, springs and other source waters at various elevations at Yude Ri and Dungju Ri catchments, Bhutan Himalayas. Samples were analyzed for pH, specific conductance, and major ions and end-member mixing analysis in combination with diagnostic tools of mixing models was used to determine sources, relative contributions, and recharge dynamics of spring flows. The results indicated that direct precipitation dominated spring flows (0.59 ± 0.21), followed by shallow groundwater (0.31 ± 0.18), and soil subsurface water (0.10 ± 0.15). The contributions of spring flow components followed an elevation gradient, with higher and lower fractions, respectively, of direct precipitation and shallow groundwater at higher elevations, e.g., 0.90 ± 0.1 to 0.13 ± 0.08 for direct precipitation and 0.03 ± 0.03 to 0.37 ± 0.19 for shallow groundwater from 3266 m to 1558 m. Spring flows primarily relied on precipitation (∼70 % from both direct precipitation and soil water), making them very sensitive to changes in precipitation. Significant contributions of shallow groundwater also indicated the vulnerability of spring flows to decreased snowfall relative to rainfall and the earlier onset of snowmelt, particularly for those located in the snow-rain transition zone (∼2500 m). Our results suggest high vulnerability of spring flows to the climate change in the Himalayas.

Graphical abstract

Image 1

Highlights

• Direct precipitation (∼60 %) & soil subsurface water (∼10 %) dominated spring flow paths.

• Shallow groundwater contributions to springs increase with a decreasing elevation.

• Spring flows are vulnerable to changes in precipitation pattern & earlier snowmelt.

Keywords

Spring flow
Water sources
Flow paths
End-member mixing analysis
Diagnostic tools of mixing models
Bhutan Himalayas
==== Body
pmc1 Introduction

Mountain springs and streams have traditionally served as the primary water sources for millions of Himalayan people for their drinking, domestic, and agricultural needs besides its cultural and religious value, and sustaining rich biodiversity [1,2]. However, there has been a reported decline in the water discharges from these resources since the 1980s across the Himalayan range, resulting in acute water shortages [[1], [2], [3], [4], [5], [6]]. Recent years have witnessed an unprecedented reduction in springs in Bhutan due to changing climate conditions [7]. It is found that about 25 % of 7399 springs assessed in Bhutan have been declining in flow, while about 1 % have completely dried up. There are numerous anecdotal evidences and newspaper articles reporting how drying up of springs and water scarcity are impacting rural population in Bhutan. While larger river basins in the Himalaya do not face immediate water-shortage problems, local communities living on mountain slopes, who rely on small spring flows and rivulets, are experiencing acute water-related issues [8].

With global warming, the 10-year mean global surface temperature was 1.1 °C higher in 2011–2020 than in 1850–1900 [9]. Mountain regions, especially the Himalayas, was much more sensitive to climate warming than lower-lying areas due to the elevation-dependent warming phenomenon [10,11]. During 1980–2018, warming of the Hindu Kush–Karakoram–Himalayan system was 0.42 °C per decade, twice the global average rate [12]. The warming has led to uncertainties regarding the availability of mountain water resources, with glaciers losing mass, snowmelt dynamics being disrupted, and precipitation and evapotranspiration patterns shifting [[12], [13], [14]]. Chandel and Ghosh [15] demonstrated that with increasing temperatures, the contribution of snowmelt to stream flow decreases substantially due to declining snowpack storage despite higher melt rates of glaciers in the Himalayas. A recent study in the Chamkhar Chhu River basin in Bhutan Himalaya showed a significant decrease in annual minimum snow and ice extent [16]. The changing precipitation patterns due to rising temperatures may further impact the hydrologic system, resulting in significant changes in water resources in the region, including small streams and mountain springs [17].

The situation in Bhutan is particularly concerning because among the total employed work force, the highest percentage is engaged in the agriculture sector (43.5 % of the population as of 2022 [18]), contributing ∼15 % of the nation's economy [19]. Despite having a high per capita water use capacity, most rivers and streams in Bhutan are inaccessible as they flow at the bottom of gorges and ravines, while settlements and farmland occupy higher ground [20]. Consequently, springs and mountain streams are the primary water sources for domestic and agricultural uses, including urban water supplies. However, climate change, population growth, and economic development are exerting pressures on these water resources [20,21]. Therefore, there is an urgent need to gain insights into spring and stream flow generation, including dominant water sources and their recharge areas. Due to the limited research capacity, budget constraints and scarcity of data, however, comprehensive hydrological studies, including research on climate change impacts on water resources in Bhutan, are lacking [22]. This problem is also exacerbated by difficult terrains and remote access with little infrastructure to support field studies. As such, only preliminary investigations on springs, such as collecting and archiving discharge data, have been conducted for a very small portion of mountain springs in Bhutan [8]. Nevertheless, a few spring revival experiments in eastern Himalayan regions including Bhutan have been initiated to adopt geohydrologic approaches by identifying the recharge areas (springshed) and taking up artificial groundwater recharge works in an attempt to solve this intractable challenge [23]. Though the springshed development approach has proven to be promising to rejuvenate springs in some Himalayan catchments, the major challenge lies in the accurate identification of the recharge areas [2]. Due to the complexity of mountainous terrain, particularly, the groundwater flow system feeding springs is poorly understood, which makes the springshed approach difficult to scale up over larger areas. Thus, information on water sources of spring flow (such as groundwater, snowmelt, rainwater and soil water) and their pathways (local or regional) is urgently needed for an in-depth understanding of inner working of mountain aquifers system [1,2,[24], [25], [26]].

Hydrologic mixing models that rely on natural isotopic and geochemical tracers have been widely and successfully applied to identify stream flow sources and catchment flow paths in mountain catchments [17,22,[27], [28], [29]], including remote headwater river basins in the Himalaya [17,22]. It serves as a unique tool in understanding how catchment hydrologic and biogeochemical properties and processes vary across scale. End-member mixing analysis (EMMA) [30], in particular, is a well-established methodology that has been widely used to identify sources of water responsible for stream flow generation in hillslopes, small catchments, and larger watersheds spanning a range of geographic, geologic, climatic, and environmental conditions [27,28,[30], [31], [32]]. Nevertheless, EMMA concept has also been extended to the interpretation of geochemical fluxes using spring flow chemistry. Frisbee et al. [29] successfully employed EMMA to illustrate that spring flow generation, like stream flow generation by Liu et al. [32] and Suecker et al. [33], integrates many different sources of water reflecting solute concentrations obtained along many different geochemical weathering pathways or sub-surface flow paths.

In this study, diagnostic tools of mixing models Hooper [31] and endmember mixing analysis (EMMA) were combined to determine the spring flow generation processes using a yearlong geochemical dataset collected in the remote and climate sensitive catchment in Bhutan, Eastern Himalaya. The specific objectives of the study were to quantify the contributions of direct precipitations (snowmelt and rainwater), groundwater and other source waters to spring flows using geochemical tracers, and understand the recharge mechanism of spring flows. Furthermore, this study highlighted the long-term impact of climate change on the fate of mountain springs to assist in developing climate-resilient intervention measures.

2 Materials and methods

2.1 Study area

The study area encompasses portions of Merak, Radhi, and Phongmey Gewogs (administrative unit or block) within the Trashigang District of Eastern Bhutan. Spanning 66.43 km2, the majority of the study area lies within Radhi Gewog (42.53 %), followed by Merak (31.07 %) and Phongmey (26.40 %). Elevations range from 3421 m above sea level (m asl) at Mindrula peak to 1006 m asl at the Gamri River. The study area forms part of the Gamri watershed, a crucial tributary of the Drangme Chhu, which ultimately flows into the Brahmaputra River in India and then into the Bay of Bengal (as depicted in Fig. 1). This study focuses on the water resources within two tributary catchments of the Gamri: Streams Yude Ri and Dungju Ri (Fig. 1). These streams originate from high-altitude Mindrula, Cheabling, and Shetymi rangelands of Merak Gewog. Springs scattered across the area contribute to the streams or dissipate before reaching them. Natural ponds and wetlands are also present within the study area.Fig. 1 Research location map showing sampling sites in the study area. See Tables S1 and S2 for all site acronym and abbreviation.

Fig. 1

Farming and livestock husbandry are the primary economic activities in the study area. Notably, Radhi is renowned for its rice production, earning the moniker "rice bowl of Eastern Bhutan" [34]. This high-water-demand crop underlines the vital role of water resources for subsistence farmers. The annual paddy production by about 932 households residing within the study area was estimated to be 2800 metric tons (1 metric ton = 1000 kg) in 2022. Meanwhile, the upper rangelands, particularly Cheabling, Mindrula and Shetymi serve as winter grazing grounds for over a hundred Brokpa communities. These semi-nomadic tribes rely on livestock rearing, particularly yaks, for their livelihood, indicating a substantial demand for fresh water in the upper study area as well [35]. Dungju Ri serves as the primary source of irrigation water for several downstream villages in Phongmey Gewog and an important drinking water source for livestock in the upstream areas. While Yude Ri itself does not directly contribute to irrigation, its tributaries play a critical role in irrigating paddy fields in Radhi Gewog [8,34]. Additionally, Yude Ri's tributary streams and springs are main source for drinking water and livestock rearing for the nomads and upstream villages. Recognizing the importance of paddy cultivation and livestock rearing in supporting local livelihoods, various water source protection and land management initiatives have been implemented in the study area. However, the study area lacks hydrologic data (e.g., spring flow discharges) and knowledge for implementing effective water resource management programs.

2.2 Hydroclimate

Hydroclimate in the study area is heavily influenced by the Indian summer monsoon, a characteristic feature of the eastern Himalayan region. The monsoon season typically commences in June, peaks in July and August, and concludes in September. This period contributes approximately 70 % of the annual precipitation, leading to exceptionally high river flows coinciding with the warmest months of a year (Fig. 2). A 28-year record (October 1, 1996 to September 30, 2023) from the Radhi meteorological station (Latitude: 27.37⁰; Longitude: 91.70⁰; and Elevation: 1576 m asl), the only available station near the study area, reveals an average annual rainfall of 1113.8 mm, ranging from 323.4 mm in 2013 to 2224.7 mm in 2017 water years. Historically, daily maximum and minimum temperatures at the station varied from 2.0 °C to 37.0 °C and from −1.0 °C to 30 °C, respectively. Additionally, higher elevations within the study area experience snowfall usually about 30 cm thick (snow depth measurements data are not available) during winter (December to February) above 2500 m asl when minimum daily temperatures fall below 0 °C. This accumulated snow begins to melt by spring (March to May) as temperatures rise. Notably, spring snowmelt coincides with the onset of monsoon rains, further intensifying the high flow period.Fig. 2 Hydroclimatic data from 1996 to 2023, along with data from water year 2022, showing annual fluctuations for (top) precipitation and (bottom) air temperature for the study area.

Fig. 2

2.3 Land use and land cover and geology

The study area is characterized by a significant altitudinal gradient, resulting in diverse microclimates, biodiversity, and agricultural practices (Fig. 3). Forest cover dominates the landscape, accounting for 63.6 % of the land area, with broadleaf forests (predominantly Quercus, Castanopsis and Acer species followed by Abies spp. and Betula spp., and sporadic chirpine forest in the lower region) being the most prevalent. Cultivated agriculture encompasses 22.6 %, with dry land ("Kamzhing") and paddy fields ("Chhuzhing") being the main types. Meadows, concentrated primarily within Merak Gewog, cover approximately 6.3 % of the area and experience intensive grazing pressure. Unsustainable rangeland management practices like tree lopping for winter fodder have contributed to forest degradation, overgrazing, and increased vulnerability to landslides and flash floods during the monsoon season [35].Fig. 3 Land use and land cover (LULC) of the study area in 2016. Data were provided by the Forest Resources Management Division, Department of Forests & Park Services, Ministry of Agriculture and Forests, Bhutan.

Fig. 3

The study area lies within the Tethyan Himalaya tectonostratigraphic zone, characterized by rocks of the Ordovician- or younger Cheka Formation exhibiting a remarkable thickness of 2.2–3.5 km and constituting a crucial lithological element in Bhutan. The dominant lithologies are tan to grey, thick-bedded, fine-to medium-grained micaceous quartzites, forming imposing cliffs. Interlayered with these quartzites are biotite-muscovite-garnet schists, as documented by Long et al. [36].

2.4 Sample collection and analysis

The specific conductivity (SC), pH, and water temperature of main streams, and potential source waters such as tributary streams, springs, and natural ponds encountered in the inaugural field expedition across the Yude Ri and Dungju Ri catchments were measured in situ in April 2022. Based on the results of these field measurements, 42 sampling sites were selected, comprising main streams (n = 12), tributary streams (n = 7), perennial springs (n = 18), natural ponds (n = 4), and wetland outflow (n = 1) along the elevation range of 1045 m to 3471 m (Fig. 1, Tables S1 and S2). In addition, five rainwater collectors (P1–P5) were installed between 1582 and 3378 m asl in the end of April and early May 2022 before the onset of the monsoon. The rainwater samples were collected monthly from the end of May to September 2022, which represent monthly, aggregated precipitation samples for chemical analysis. Also, five snow core samples were collected in the Yudi Ri catchment, four in February 2022 and one in January 2023, respectively, with the COVID-19 restriction in the area and logistical limitations. Three homemade zero-tension pan-lysimeters (Abbreviated ID with LY in Table S1) were installed in May 2022 to collect soil water samples at the depth of 0.1 m, 1 m, and 1.5 m near a stream bank or riparian zone at DS11 (2445 m) of Dungju Ri catchment. The shallowest lysimeter (LY3) was installed just under the O-horizon near the edge of a riparian area. Two piezometers (Abbreviated ID with PZ in Table S1) and two wells (Abbreviated ID with WL in Table S1) were co-installed at the depths of 1 m and 1.5 m for sampling soil water and integrated subsurface water, respectively. The soil water and subsurface water samples were collected from June 2022 through January 2023.

The grab water samples were filtered through Millipore 47-mm glass-fiber filters (pore size = 0.45 μm) and stored in 125-mL clean high-density polyethylene (HDPE) bottles following the protocols of Wilson et al. [37]. The samples were kept frozen in the freezer until analysis. All water samples were analyzed for major ions (Ca2+, Mg2+, Na+, K+, Cl−, NO3−, and SO42−) at the water chemistry laboratory in Sherubtse College of the Royal University in Bhutan. Analyses of major cations and anions were performed using a Metrohm 930 Compact Ion Chromatograph. For the quantification of anions, a Metrohm Metrosep A Supp 5–150/4.0 column (6.1006.520) was utilized with an eluent consisting of 3.2 mmolL−1 sodium carbonate and 1 mmolL−1 bicarbonate. The eluent was pumped at a flow rate of 0.7 mL/min at a column temperature of 25 °C and a conductivity detector. The system has Metrohm CO2 Suppressor (MCS) that lowers the background conductivity and improves detection sensitivity. For quantification of cations, a Metrohm Metrosep C4 - 250/4.0 column (6.1050.430) was used with an eluent composed of 1.7 mmolL−1 nitric acid and 1.0 mmolL−1 dipicolinic acid, pumped at a flow rate of 0.9 mL/min at 30 °C. The ultrapure water with a specific resistance no lesser than 18.2 MΩ cm (25 °C) was used. The detection limits for all ions were less than 1 μEqL−1. Acid neutralizing capacity (ANC) was calculated by charge balance of major ions as ANC = [ Ca2+ + Mg2+ + Na+ + K+] – [Cl− + NO3− + SO42−], all in μEqL−1. A strong correlation between calculated ANC and measured Ca2+ (R2 = 0.86, p < 0.001), without obvious outliers (figure not provided), ensured that calculated ANC values are reasonably accurate.

2.5 Diagnostic tools of mixing models and end-member mixing analysis

The present study utilized a combination of diagnostic tools of mixing models (DTMM) [31] and End-Member Mixing Analysis (EMMA) [30] to assess conservative tracers, determine the number of end-members, and quantify the contribution of end-members to spring flows. For DTMM, briefly, a principal component analysis (PCA) was conducted with all analytes to extract eigenvectors and analyte values in spring flows were projected using the eigenvectors without using any information from end-members. Relative root mean square error (RRMSE) and distributions of residuals between measured and projected analyte values were used to determine the dimension of a lower mixing space following Hooper [31]. In the meantime, the cumulative variances explained by dimensional mixing spaces were calculated by the cumulative, squared, scaled eigenvectors in conjunction with eigenvalues and used as benchmark or expected R2 values in the EMMA validation (described below) following Liu et al. [38] and Porter et al. [39]. In EMMA, a PCA was conducted again but using conservative tracers determined by DTMM. A mixing diagram was developed using principal components derived from the conservative tracers to identify eligible end-members contributing to spring flows [27]. Four criteria were utilized to screen end-members, as described in Liu et al. [32,38,40] and Porter et al. [39] (a) whether or not they form a polygon (e.g., a triangle in the case of three end-members): with all spring water samples inside the polygon; (b) end-member distances must be relatively short compared with others in proximity for all tracers used in the analysis; and (c) spring water chemistry must be well simulated using fractions of end-member contributions and their solute concentrations in end-members; and (d) the EMMA results must be hydrologically meaningful. In the case of three end-members, the first two principal components were used to solve for end-member contributions as in a traditional two-tracer for three-component mixing model following the simultaneous equations of Liu et al. [38]:(1) 1=f1+f2+f3

(2) As=f1A1+f2A2+f3A3

(3) Bs=f1B1+f2B2+f3B3

where f is the fraction of total spring flow due to an end-member; A and B represent the first two principal components; subscripts 1, 2, 3 and s represent end-members 1, 2, 3 and spring flow. Once three end-members are identified, their sequence to plug into the equations does not matter as long as fractions and concentrations of all end-members keep consistent manner. The simultaneous equations were solved by a matrix, consisting of end-member compositions and constraints, following the procedures described in Shaw et al. [41]. The EMMA results were validated by simulating spring flow chemistry following Liu et al. [27] and Liu et al. [32]. Mathematically, the spring flow chemistry was projected using end-member fractions determined above, but the principal components in the above equations were substituted by measured chemical concentrations in end-members. The R2 and slope values obtained from a correlation between measured and projected solute concentrations were then compared against the benchmark or expected R2 values from DTMM discussed above.

3 Results

3.1 Hydrochemical characteristics

3.1.1 Solute concentrations in stream water, precipitation, soil water & springs

Median solute concentrations and SC values were much higher in stream water (main streams and tributaries), soil water and springs than in precipitation (rain and snow), natural pond and wetland outflow for all solutes except K+ (in natural pond only) and NH4+ (Table S3, Fig. 4). Particularly, the median solute concentrations of weathering products such as ANC, Ca2+, Mg2+ and Na+ were much higher in springs and soil water due to an extended rock-water interaction. For example, the median ANC concentrations were 406.5 and −0.82 μEqL−1 in spring and snow, respectively. The SC values and Ca2+, Mg2+, Na+ and ANC concentrations followed a similar pattern of variation from stream and spring waters to soil water and precipitation, with highest medians in spring water. The median Cl− concentration in springs was also highest, 11.5 μEqL−1, followed by that in soil waters (7.7 μEqL−1) and main stream water (6.2 μEqL−1). The median NH4+ concentration was 18.9 μEqL−1 in natural pond samples, followed by snow (17.1 μEqL−1), but 6–10 times of those in rainfall, stream water, soil water, springs and wetland outflow. The median NO3− concentrations in springs, soil water and stream water varied minimally, ranging from 17.2 to 21.6 μEqL−1, but 2–7 times of those in all other samples except natural pond water (>25 times). Water temperature did not vary as much as solute concentrations for all types of samples. Other than water temperature and NH4+, solute concentrations did not vary much between rainfall and snow samples.Fig. 4 Box-whisker plots showing variation of water temperature and SC values and solute concentrations in stream water at main streams and tributaries, precipitation, soil water, natural ponds, wetland outflows, and spring water.

Fig. 4

Solute concentrations in spring water varied over locations and seasons. Stream temperature, SC value, and solute concentrations tended to increase with a decrease in their sampling elevations for all but NH4+ (Fig. 5). For example, the median concentrations of Ca2+ and ANC were 71.1 and 115.8 μEqL−1 from the spring YS16 (3266 m), respectively, but 428.1 and 933.8 μEqL−1 from the spring YS3 (1558 m). At lower elevations, however, some sites did not follow well the trend, with lower concentrations mostly for SC, Ca2+, Mg2+, Cl−, NO3−, and SO42− at YS2, DS18, DS19, DS1, and DS20 than those at similar elevations.Fig. 5 Spatial variation of temperature and SC values and solute concentrations in spring water. Sites are arranged from higher to lower elevations from left to right on each panel.

Fig. 5

The seasonal SC value and solute concentrations in most springs were higher in pre-monsoon and post-monsoon seasons coinciding with lower flows (Fig. 6; flows were not measured but expected to be much lower during the pre- and post-monsoon seasons than the monsoon season). However, solute concentrations in lower elevation springs tended to be higher during the monsoon season at DS2, DS3, and DS4, whose median concentrations of SC, Ca2+, Mg2+, Cl−, NO3−, and SO42− in Fig. 5 were usually higher than those at similar elevations.Fig. 6 Seasonal variation of SC values and solute concentrations in spring water (ionic charges were omitted for accommodating the limited space on x-axis).

Fig. 6

3.1.2 Hydrochemical facies

The spring water samples were analyzed and represented on a Piper trilinear diagram using AqQA software (Fig. 7). This diagram allowed for a comprehensive assessment of the hydrogeochemical facies and the hydrochemical evolution of the springs, following the methodology introduced by Piper [42]. The main hydrogeochemical facies identified were characterized by Ca–Mg–HCO3 type. Most of spring samples were clustered, indicating a similar origin and evolution of their geochemistry. Within the cationic field of the Piper diagram, most spring samples were found in zone B, with a smaller portion falling into zones A and C, indicating the absence of a single dominant cation water type. Similarly, most of the plotted spring samples fall into bicarbonate type (zone E) in the anionic field, consistent with Ravikumar and Somashekar [43].Fig. 7 Piper diagram depicting hydrogeochemical facies of spring water.

Fig. 7

3.2 End-member mixing analysis

3.2.1 Conservative solutes and the number of end-members

All eight solutes, SC, and water temperature listed in Table 1 were included in diagnostic tools of mixing models to identify conservative tracers and the number of end-members following Hooper [31]. Spring samples collected at various locations were combined to extract eigenvectors in a principal component analysis, as the origins of their hydrochemistry were similar (Fig. 7). The results indicated that the residuals between projected and measured values did not have a random distribution (p > 0.001) with measured values under a lower dimensional mixing space (usually 2-D) for all but NO3− (Table 1). This result suggests that other than NO3− all did not ideally behave conservatively upon mixing of three end-members in spring flow. However, the relative root-mean square errors (RRMSE) were lower than 5 % and 4 % for all analytes in 1-D and 2-D mixing spaces, respectively (Table 1), which were comparable or lower than those in other studies [38,44]. In addition, the first two principal components together (PC1 + PC2) explained on average 79 % of the total variance of ten tracers, while the third (PC3) added only 7 % (not listed in Table 1 but can be calculated as the sum of the differences between 3-D and 2-D values and then divided by the number of tracers). For individual tracers, PC1 and PC2 together explained at least 61 % of the variance with higher values for NO3−, SC, ANC, Ca2+, and Cl− (84–95 %). For this study, the primary purpose was to identify major end-members contributing to spring flow rather than all of them. Thus, 2-D mixing space (three end-members) was selected and all ten tracers including SC and water temperature were used for subsequent EMMA analysis. Following Porter et al. [39], the values of cumulative variances explained in 2-D were then converted to benchmark or expected R2 values between measured and projected values using EMMA results as shown in 3.2.2 below to evaluate the performance of EMMA and determine if the major end-members were indeed correctly identified (the conversion can be done simply by dividing the 2-D values by 100 for individual tracers; see the benchmark R2 listed in Table 2).Table 1 The results of diagnostic tools of mixing models using all analytes (except NH4+) in spring water (n = 162).

Table 1	Dimension	Temperature	SC	Ca2+	Mg2+	Na+	K+	Cl−	NO3−	SO42-	ANC	
Cumulative Variance Explained (%)	1-D	54	88	86	77	74	61	62	29	62	86	
2-D	63	89	86	79	76	61	84	95	64	89	
3-D	70	91	91	95	76	80	91	96	73	97	
4-D	97	91	92	95	76	88	93	96	83	97	
5-D	97	96	93	96	86	95	93	96	97	99	
R2 from Residual distributions	1-D	0.46	0.12	0.14	0.23	0.26	0.39	0.38	0.71	0.38	0.14	
2-D	0.36	0.11	0.14	0.22	0.25	0.39	0.16	0.05	0.36	0.11	
3-D	0.30	0.09	0.09	0.05	0.24	0.20	0.09	0.04	0.27	0.03	
4-D	0.03	0.09	0.08	0.05	0.24	0.12	0.07	0.04	0.17	0.03	
5-D	0.03	0.04	0.07	0.04	0.14	0.05	0.07	0.04	0.03	0.01	
RRMSE (%)	1-D	0.98	1.33	1.58	2.45	1.95	2.82	4.38	5.05	3.82	1.56	
2-D	0.89	1.28	1.56	2.36	1.89	2.81	2.84	1.35	3.72	1.39	
3-D	0.79	1.16	1.27	1.11	1.89	2.00	2.09	1.23	3.20	0.71	
4-D	0.25	1.13	1.18	1.10	1.87	1.58	1.90	1.20	2.59	0.66	
5-D	0.23	0.75	1.15	0.98	1.44	0.99	1.89	1.19	1.03	0.34	
Note: The bold numbers indicate p > 0.001.

Table 2 Summary of EMMA model evaluation, including expected R2 and actual R2, slope and intercept from regression between measured and projected water temperature and SC values and solute concentrations using EMMA results.

Table 2	Temperature	SC	Ca2+	Mg2+	Na+	K+	Cl−	NO3−	SO42-	ANC	
Expected R2, essentially the cumulative percent of variance explained by 2-D in Table 1 (but divided by 100)	
R2	0.63	0.89	0.86	0.79	0.76	0.61	0.84	0.95	0.64	0.89	
M1: Precipitation, Soil Water (LY2, January 7, 2022), and Baseflow (YS3, 5/31/2022)	
R2	0.61	0.91	0.83	0.77	0.71	0.53	0.82	0.86	0.62	0.87	
Slope	0.75	0.90	0.81	0.74	0.77	0.74	0.43	0.97	0.37	0.97	
Intercept	0.65	5.87	36.78	45.72	26.10	15.04	7.19	3.13	28.12	62.33	
M2: Precipitation, Soil Water (LY2, January 7, 2022), and Baseflow (DS0, February 7, 2022)	
R2	0.68	0.81	0.77	0.63	0.60	0.55	0.85	0.90	0.61	0.69	
Slope	1.18	0.69	0.70	0.41	0.40	0.78	0.44	1.15	0.18	0.61	
Intercept	−1.76	10.18	48.02	45.93	39.57	19.03	6.55	0.96	26.19	112.81	

3.2.2 Determination of end-members and their contributions

Per the above DTMM results, a mixing diagram was developed using PC1 and PC2 developed from all ten tracers in spring flow, along with projected principal components of precipitation (rain and snow), soil water and baseflow samples (Fig. 8). Many combinations of three potential end-members could be considered to be viable triangles from the geometrical perspective. Two triangles were identified to be the best to bound most of spring samples, with three vertices represented by precipitation (mean snow value), soil water (from mid lysimeter, LY2, January 7, 2022) and either a spring sample from YS3 (5/31/2022) (M1 model in Table 2) or a baseflow sample from DS0 (mainstream outflow, February 7, 2022) (M2 model in Table 2). It appears that the baseflow sample from DS0 plotted in a better position to serve as one of the triangle vertices (M2 model) as it bounds more spring samples inside the triangle than the one from YS3 (M1 model). Based on the end-member distances (Table 3), however, YS3 was surprisingly much better than DS0. The end-member distances of YS3 calculated for all the tracers (0 %–31 %) were much lower than those of DS0 (12 %–732 %).Fig. 8 Mixing diagram of all springs using the first two principal components (PC1 and PC2) constructed from all ten tracers. Individual samples of potential end-members were shown, except for precipitation, which was averaged for each sampling location. Triangles show examples of potential models with varying end-members. Samples of main streams, tributaries, and natural ponds were not shown for clarity.

Fig. 8

Table 3 Endmember distances of tracers used in EMMA for selected potential endmembers.

Table 3Sample Type	Location Code	Date	Temperature (%)	SC (%)	Ca2+(%)	Mg2+(%)	Na+(%)	K+(%)	Cl−(%)	NO3−(%)	SO42-(%)	ANC (%)	
Precipitation (Selected)	Snow	Mean	60	17	10	Inf	28	20	192	71	134	39	
Precipitation (Others)	P1–P5	Mean	29	112	337	1333	492	143	173	94	76	4259	
Soil Water (Selected)	LY2	7/1/22	1	11	4	16	59	45	81	16	44	14	
Soil Water (Others)	All	Mean	25	13	24	30	307	77	258	103	57	40	
Baseflow 1	YS3	5/31/22	9	0	5	3	8	9	31	18	13	10	
Baseflow 2	DS0	7/2/22	27	19	12	47	25	22	34	732	18	28	
Note: “Inf” refers to infinite value due to dividing by zero (Mg2+ was not detected in snow samples).

Precipitation (snow) and soil water sample from the mid lysimeter were well positioned as two triangle vertices to bound most of the spring samples and characterized as surface runoff directly from precipitation and soil subsurface flow, respectively. Their end-member distances were reasonably short for all tracers, particularly in comparison with other studies [32,38,39]. Not only was snow geometrically located in a better position than other rainwater samples (Fig. 8), but also its end-member distances were much shorter (Table 3). However, these results did not suggest snow is the sole end-member representing direct precipitation. The hydrochemistry of rainwater did not differ much from that of snow (Fig. 4). In fact, rainwater was located in the vicinity of snow in the mixing diagram (Fig. 8). Using nine tracers without water temperature, they became even closer. Thus, snow end-member represents direct precipitation from both snowmelt and rainwater and their contributions to spring flow cannot be separated with hydrochemical tracers used in this study.

Projection of spring flow chemistry indicated that EMMA reasonably well reproduced the measured concentrations (Fig. 9 and Table 2). However, the projected spring flow chemistry was much better for M1 than M2, as the projected and measured solute concentrations had R2 values much closer to the expected R2 values from M1 model than from M2 model for most tracers such as SC, Ca2+, Mg2+, Na+, ANC, and SO42− (Table 2). The slopes from M1 model were also closer to 1.0 than those from M2 model, particularly for SC, NO3−, and ANC. In the mixing diagram, the position of YS3 spring was relatively far from snow and soil water end-members (Fig. 8). Solute concentrations in YS3 spring were highest among all spring sites (Fig. 5). Thus, YS3 (5/31/2022) was characterized as shallow groundwater or baseflow for this study in the Dungju Ri and Yude Ri stream catchments. Conceptually, we think spring at YS3 is primarily fed by deeper groundwater over longer flowpaths and a baseflow sample from YS3 collected during a dry period, during which rainfall impact dissipates, best characterizes the chemistry of deeper groundwater. The improved representation of spring flow chemistry in M1 compared to M2 bolstered the confidence in the selected end-members. Consequently, the M1 model was deemed to be the most suitable to represent the mixing model for spring flows, with direct precipitation, soil water, and baseflow before the monsoon season as the major end-members.Fig. 9 Projection of ionic concentrations in spring flow for selected solutes using relative contributions of endmembers determined by EMMA (the results from the M1 model as an example) and ionic concentrations in endmembers.

Fig. 9

The EMMA results were summarized in Table 4, where the simultaneous equations (1), (2), (3) were applied to all individual samples and the ‘±’ values represent one “standard deviation” of fractional contributions of all samples for respective end-members of each spring. The results indicated that spring flows in the Dungju Ri and Yude Ri catchments were dominated by direct precipitation (snowmelt and rainwater) with a mean fractional contribution of 0.59 ± 0.21 in the water year 2022. The contribution of shallow groundwater (baseflow) was second, with a mean of 0.31 ± 0.18, while subsurface water contribution was least with a mean of 0.10 ± 0.15. Its source water was dominated by precipitation-derived flow components (69 % combining both direct precipitation and soil subsurface water).Table 4 The fractional contributions of three end-members to the spring flows.

Table 4Location Code	Elevation (m)	Direct Precipitation	Soil Water	Baseflow	
YS16	3266	0.90 ± 0.01	0.07 ± 0.04	0.03 ± 0.03	
YS13	2626	0.85 ± 0.02	0.03 ± 0.05	0.12 ± 0.06	
DS8	2608	0.75 ± 0.05	0.13 ± 0.09	0.11 ± 0.07	
YS7	2361	0.84 ± 0.02	0.00 ± 0.00	0.16 ± 0.02	
YS6	2192	0.77 ± 0.01	0.00 ± 0.01	0.23 ± 0.01	
YS5	2148	0.75 ± 0.03	0.00 ± 0.00	0.25 ± 0.03	
YS21	2111	0.59 ± 0.05	0.07 ± 0.03	0.35 ± 0.06	
DS15	1886	0.57 ± 0.03	0.01 ± 0.03	0.42 ± 0.04	
DS16	1847	0.54 ± 0.04	0.01 ± 0.02	0.45 ± 0.04	
YS2	1841	0.65 ± 0.03	0.00 ± 0.00	0.35 ± 0.03	
DS18	1834	0.61 ± 0.02	0.15 ± 0.07	0.24 ± 0.07	
DS4	1793	0.34 ± 0.14	0.38 ± 0.26	0.28 ± 0.13	
DS19	1775	0.69 ± 0.04	0.06 ± 0.06	0.25 ± 0.04	
DS3	1771	0.38 ± 0.11	0.26 ± 0.18	0.36 ± 0.09	
DS2	1753	0.30 ± 0.16	0.39 ± 0.19	0.31 ± 0.19	
DS1	1719	0.47 ± 0.05	0.01 ± 0.04	0.52 ± 0.06	
DS20	1702	0.59 ± 0.11	0.04 ± 0.11	0.37 ± 0.10	
YS3	1558	0.13 ± 0.08	0.11 ± 0.05	0.76 ± 0.13	
Mean		0.59 ± 0.21	0.10 ± 0.15	0.31 ± 0.18	

The temporal variation of the end-member contributions to spring flow was noteworthy for springs DS2, DS3 and DS4 (Fig. 10 and also by their standard deviations in Table 4). In these springs, the fractional contributions of direct precipitation were higher in April and then gradually decreased through August and finally showed a gradual rise till January (Fig. 10). The opposite trend was observed for fractional contribution of shallow groundwater to these springs. Meanwhile, subsurface water was the dominant contributor to spring flows in July, August and September for DS2, DS3 and DS4, consistent with the peak monsoon season. These springs, situated within a range of 1753 m to 1793 m, are in close proximity to each other and emerge within the paddy fields downstream of the Dungju Ri catchment.Fig. 10 Spatial and temporal variation of fractional contribution of end-members to spring flows.

Fig. 10

The mean fractions from all samples collected over time were used to correlate with elevations of their sampling locations (Fig. 11). The contributions of end-members, especially direct precipitation and shallow groundwater, varied spatially within the study area. Mean fractional contribution of direct precipitation (surface runoff) to spring flow decreased significantly with a decreasing elevation from 0.90 ± 0.1 at 3266 m to 0.13 ± 0.08 at 1558 m (y = 0.88In(x) - 6.08, R2 = 0.65, p < 0.01). On contrast, mean fractional contribution of shallow groundwater to spring flow increased significantly with a decreasing elevation, with 0.03 ± 0.03 at 3266 m to 0.37 ± 0.19 at 1558 m (y = −0.71In(x) - 5.70, R2 = 0.65, p < 0.01).Fig. 11 Variation of fractional contribution of (a) direct precipitation and (b) shallow groundwater end-members to the spring flows with elevations.

Fig. 11

4 Discussion

4.1 Source waters and flow paths controlling spring flows

This study showed that spring flows were primarily generated from surface and shallow subsurface flow paths at higher elevations and by groundwater at lower elevations (Fig. 11), which fits well to the spring flow path models suggested by Frisbee et al. [28], Taucare et al. [45], and Somers and McKenzie [46]. As shown in the conceptual model (Fig. 12), at high elevation sites, the steep topography with poorly developed soil of the study area is conducive to fast runoff processes such as surface runoff and shallow subsurface flow, especially during snowmelt and intense thunderstorm events as reported in previous studies [26,28,[47], [48], [49]]. As precipitation contribution constitutes a significant component of spring flow at higher elevation (between 2148 m (YS5) to 3266 m (YS16)), geochemical transformations during the fast runoff process appear to be minimal likely due to short or local flow paths which does not allow substantial interaction with soil or underlying lithology. Frisbee et al. [29] described that the fractured bedrock outcrops with interconnected fractures and or/preferential flow paths in the soil as probable mechanism for the dominance of rainfall and fast runoff components in the spring flow. Similar findings were reported in the other parts of the monsoonal Himalayan catchments [26,50] and also in Western Andean Front of Central Chile [45]. In contrast, the importance of groundwater increases with increasing scale in drainage areas. It appears that more geochemically evolved groundwaters with relatively longer residence times, likely due to deeper soils, are discharging to the low-elevation springs as depicted in Fig. 11, Fig. 12, which is consistent with the 3-D catchment-mixing conceptual model proposed by Frisbee et al. [28] for streamflow generation at the large watershed scale. The dominance of Ca2+, Mg2+, and HCO3− ions in spring flow samples demonstrated sufficient recharge from shallow freshwaters [51,52]. Recharging water occurs as water percolates through the subsurface; it carries dissolved carbonate in the form of HCO3− and geochemically mobile Ca2+ and Mg2+ obtained from the dissolution of carbonate rocks such as limestone and dolomite [43]. Thus, groundwater recharging spring flows are shallow groundwater in nature, not deep groundwater. This result is also supported by the fact that all springs seem to share a similar origin to their hydrochemistry (Fig. 7), consistent with past studies of the Himalayan springs [[53], [54], [55]]. Furthermore, weakly evolved geochemistry of spring water samples implies that vegetation cover in the study area consumes the storage of soil water through root uptake, which in turn shortens the mean residence times of soil water and groundwater [56].Fig. 12 Conceptual diagram illustrating the source waters and flow paths controlling the spring flows.

Fig. 12

Lucianetti et al. [57] reported a spatial variability of snowmelt vs rainwater contributions to spring water in dolomitic mountain group in Italian Alps. In the Italian Alps, snowmelt contribution to the spring recharge was predominant (72 ± 29 %) in the high-altitude recharge areas above 2500 m, while springs below 2000 m were mostly recharged from rain. In the present study, though the contribution of rainwater is not ascertained separately from snowmelt, it is obvious that springs at the lower elevations were recharged by rainfall as snowfall in the study area occurs mainly above 2500 m. Jeelani et al. [54] investigated the temporal variation of contributions to spring flow form snow, glacier and rain in Western Himalaya, and found the dominance of snowmelt (55–96 %) at the start of melt season (March–June), and from glaciers and high altitude snowpacks (>3000–4000 m) from June to October (5–36 %), with much less fractional contributions from rainfall (4–34 %).

The variability of groundwater end-member contributions to spring flows in the study area clearly indicated that spring flow generation is topography-driven in which local, intermediate, and regional groundwater flow paths (Fig. 12) develop in basins as a consequence of the distribution of topographic features and energy gradients [28,58]. Meanwhile, the increases in soil water contribution to spring flow at the lower elevation catchments (DS2, DS3 and DS4 shown in Fig. 10) in the study area was apparently affected by the paddy fields, demonstrating that the subsurface was at saturation at the time due to recharge from irrigation or heavy summer monsoon rains, consistent with study of Liu et al. [32]. It is possible that fine textured clay and loam soils with well-developed organic layers and few rock fragments in the faddy fields in the study area facilitate the occurrence of infiltration- or saturation-excess over land flow due to its low hydraulic conductivity [40].

4.2 Implication in the Light of climate change

Our study highlights the significance of precipitation, including snowmelt and rainwater, as a major contributor to spring flows in the study area. Numerous studies predicted significant reductions in snow cover duration, snow cover area, and snowmelt volumes in mountainous regions like Himalayas due to climate change [16,57,[59], [60], [61], [62]]. These changes will have a profound impact on groundwater recharge and the water yield of springs [63]. For instance, under extreme RCP8.5 scenario, the climate models project reductions in annual snowfall by 30–50 % in the Indus Basin, 50–60 % in the Ganges Basin, and 50–70 % in the Brahmaputra Basin by 2071–2100 [64]. With more precipitation falling as rain instead of snow, rainwater flows downstream much faster, usually within hours or days after a storm, compared to a longer time frame for snow and even longer time frame for glaciers. Thereby, peak discharges of spring flow will shift earlier in spring, impacting the downstream water availability during the summer months when demand is high [46].

This study also demonstrates the importance of groundwater in sustaining spring flow during dry seasons with lower precipitation inputs, particularly in lower elevations. Huntington and Niswonger [65] have modelled the projected impacts of climate change on surface water-groundwater interactions in snow-dominated regions and found that groundwater discharge to streams was depleted during the summer, mainly due to earlier snowmelt and drainage of shallow aquifers, irrespective of increased precipitation or groundwater. They estimated that these processes resulted in a 30 % reduction in annual summer flow, implying that dry season water stress may become more severe even with an increase in annual precipitation and thus exerts more pressure on spring recharges. Similar model under the RCP 4.5 scenario with associated factors based on hydrological parameters, including topography, structural geology, land use practices/land cover and future precipitation patterns (till 2030), showed that annual spring discharges will be significantly reduced up to 50 % from 1975 to near future (2030) in Uttarakhand, Central Himalayan region, India, posing a threat to the drinking water supply [66]. This suggests a vulnerability of the water resources in our study area to future climate scenario and land use/land cover alteration.

The lower ionic concentrations in spring samples in our study area also indicated relatively short residence times for groundwater [29]. This indicates that if multiple low snow years were to occur simultaneously, the ability of groundwater to compensate for decreased precipitation inputs or earlier snowmelt would be attenuated and groundwater contributions to spring flows would be reduced. This also means that low flows during the pre-monsoon period will lead to heightened competition for limited resources, posing critical implications for water rights and water availability in our study area [67]. Furthermore, people may turn to groundwater as alternative sources of water diminish. While the prospects for snow and ice in the face of climate change are grim, shallow groundwater resources are also at risk of unsustainable use, creating an uncertain future for reliance on either storage system [68].

5 Conclusion

Spring flow is primarily generated from direct precipitation and event-related storm flow in the shallow soils, particularly at higher elevations, in the Bhutan Himalaya. Shallow groundwater contributions to spring flow, however, do increase with a decrease in elevations (or an increase in catchment areas), likely due to an increase in soil depth. Deep groundwater, with geochemistry evolved from fractured bedrock, does not seem to contribute to spring flow, at least not significantly. As a result, spring flow is strongly sensitive to any changes in precipitation and precipitation pattern in the future. A reduction in precipitation amounts and a shift of more rainfall to snowfall in the face of climate change will reduce the spring flow recharge and affect local communities for their drinking water and irrigation in the Bhutan Himalaya and beyond and thus their wellbeing.

Funding

This work was supported by the 10.13039/100000200 U.S. Agency for International Development (USAID) through the Partnerships for Enhanced Engagement in Research (PEER) program (Contract/Award number: AID-OAA-A-11-00012).

Data availability Statement

Data will be made available on request. The data used in figures and modeling of this study can be obtained from Tshewang Dendup through an email request (tshewang.sherubtse@rub.edu.bt or pelwanggatshel@gmail.com).

CRediT authorship contribution statement

Tshewang Dendup: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. Dendup Tshering: Writing – review & editing, Resources, Investigation, Funding acquisition, Conceptualization. Sonam Tobgay: Writing – review & editing, Visualization, Validation, Software, Resources, Funding acquisition. Fengjing Liu: Writing – review & editing, Visualization, Validation, Supervision, Software, Resources, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is the Supplementary data to this article.Multimedia component 1

Multimedia component 1

Acknowledgments

We express our gratitude to the U.S. National Academies of Sciences, Engineering, and Medicine and Bhutan Foundation for implementing the program. Furthermore, we extend our heartfelt gratitude to the Ministry of Finance, the Royal Government of Bhutan for facilitating this collaboration. We express our sincere appreciation to the Vice-Chancellor's Office and Sherubtse College, the Royal University of Bhutan, for their invaluable administrative support and provision of all necessary resources, including field staff for sample collection. Additionally, we extend our thanks to the Research Steering Committee, Department of Forest and Park Services, Bhutan, for granting approval to conduct research in the Yude Ri and Dungju Ri catchments. Authors are grateful to the Sakteng Wildlife Sanctuary for their logistical assistance during the field campaign. Many thanks to the National Center for Hydrology and Metrology of Bhutan for providing the hydrometeorological data.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.heliyon.2024.e36211.
==== Refs
References

1 Kulkarni H. Desai J. Siddique M.I. Rejuvenation of springs in the himalayan region Water, Climate Change, and Sustainability 2021 97 107
2 Tambe S. Reviving dying springs: climate change adaptation experiments from the Sikkim Himalaya Mt. Res. Dev. 32 1 2012 62 72 11
3 Kumar M. Rathod R. Mukherji A. Water security and spring conservation in the Himalaya Ojha H. Schofield N. Camkin J. Climate Risks to Water Security: Framing Effective Response in Asia and the Pacific 2023 Springer International Publishing Cham 15 36
4 Poudel D.D. Duex T.W. Vanishing springs in Nepalese mountains: assessment of water sources, farmers' perceptions, and climate change adaptation Mt. Res. Dev. 37 1 2017 35 46 12
5 Singh A.K. Pande R.K. Changes in spring activity: experiences of kumaun Himalaya, India Environmentalist 9 1 1989 25 29
6 Verma R. Jamwal P. Sustenance of Himalayan springs in an emerging water crisis Environ. Monit. Assess. 194 2 2022 87 35020087
7 Acharya G. Springs are drying up, say experts Kuensel 2021 Kuensel Corporation Ltd. Thimphu, Bhutan
8 National Environment Commission (NEC) National water resources inventory A Preliminary Assessment of Water Resources in the 20 Dzongkhags 2018 National Environment Commission, Royal Government of Bhutan Thimphu, Bhutan
9 Lee H. Synthesis Report of the IPCC Sixth Assessment Report (AR6), Longer Report 2023 IPCC.
10 Hock R. High Mountain Areas: in: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate 2019
11 Yadav M. Elevation-dependent precipitation in the Indian himalayan region Theor. Appl. Climatol. 155 2 2024 815 828
12 Yao T. The imbalance of the Asian water tower Nat. Rev. Earth Environ. 3 10 2022 618 632
13 Immerzeel W.W. Importance and vulnerability of the world's water towers Nature 577 7790 2020 364 369 31816624
14 Stigter E.E. Energy and mass balance dynamics of the seasonal snowpack at two high-altitude sites in the Himalaya Cold Reg. Sci. Technol. 183 2021 103233
15 Chandel V.S. Ghosh S. Components of himalayan river flows in a changing climate Water Resour. Res. 57 2 2021 e2020WR027589
16 Hill A.F. How important is meltwater to the chamkhar Chhu headwaters of the Brahmaputra River? Front. Earth Sci. 8 81 2020
17 Tshering D. Seasonal source water and flow path insights from a year of sampling in the Chamkhar Chhu basin of Central Bhutan Arctic Antarct. Alpine Res. 52 1 2020 146 160
18 National Statistics Bureau (NSB) Labour Force Survey Report 2022 National Statistics Bureau Royal Government of Bhutan: Thimphu, Bhutan. Available at: https://www.nsb.gov.bt/publications/labour-force-survey-report/
19 National Statistics Bureau (NSB) National Accounts Statistics 2023 National Statistics Bureau, Royal Government of Bhutan Thimphu, Bhutan Available at: https://www.nsb.gov.bt/publications/national-account-report/
20 Asian Development Bank/National Environment Commission (ADB/NEC) Water: Securing Bhutan's Future 1 ed 2016 Thimphu, Bhutan: Asian Development Bank/National Environment Commission, Royal Government of Bhutan Available at: https://www.adb.org/sites/default/files/publication/190540/water-bhutan-future.pdf
21 Sharma G. Conserving springs as climate change adaptation action: lessons from Chibo-Pashyor, Teesta River Basin, Kalimpong, West Bengal, India ICIMOD Working Paper, Kathmandu, Nepal 2019
22 Williams M.W. Using geochemical and isotopic chemistry to evaluate glacier melt contributions to the Chamkar Chhu (river), Bhutan Ann. Glaciol. 57 71 2016 339 348
23 Shrestha R.B. Protocol for Reviving Springs in the Hindu Kush Himalayas: A Practitioner's Manual 2018 International Centre for Integrated Mountain Development (ICIMOD)
24 Gupta A. Kulkarni H. Report of Working Group I Inventory and Revival of Springs in the Himalayas for Water Security 2018 NITI Aayog New Delhi
25 Sharma B. Springs, Storage Towers and Water Conservation in the Midhills of Nepal (ICIMOD Working Paper No. 2016/3) 2016 Nepal Water Conservation Foundation and International Centre for Integrated
26 Tarafdar S. Bruijnzeel L.A. Kumar B. Improved understanding of spring and stream water responses in headwaters of the Indian Lesser Himalaya using stable isotopes, conductivity and temperature as tracers Hydrol. Sci. J. 64 7 2019 757 770
27 Liu F. Williams M.W. Caine N. Source waters and flow paths in an alpine catchment, Colorado Front Range, United States Water Resour. Res. 40 9 2004
28 Frisbee M.D. Streamflow generation in a large, alpine watershed in the southern Rocky Mountains of Colorado: is streamflow generation simply the aggregation of hillslope runoff responses? Water Resour. Res. 47 6 2011
29 Frisbee M.D. Effect of source integration on the geochemical fluxes from springs Appl. Geochem. 28 2013 32 54
30 Christophersen N. Hooper R.P. Multivariate analysis of stream water chemical data: the use of principal components analysis for the end-member mixing problem Water Resour. Res. 28 1 1992 99 107
31 Hooper R.P. Diagnostic tools for mixing models of stream water chemistry Water Resour. Res. 39 3 2003
32 Liu F. Streamflow generation from snowmelt in semi‐arid, seasonally snow‐covered, forested catchments, Valles Caldera, New Mexico Water Resour. Res. 44 12 2008
33 Suecker J.K. Determination of hydrologic pathways during snowmelt for alpine/subalpine basins, Rocky Mountain National Park, Colorado Water Resour. Res. 36 1 2000 63 75
34 Flood Engineering and Management Division (FEMD) Flood hazard assessment for Trashigang dzongkhag Neopaney M. Thimphu, Bhutan 2019 Available at: https://www.moit.gov.bt/wp-content/uploads/2019/02/Trashigang_Flood-Hazard-Assessment_Report-1.pdf
35 Tenzing K. Millar J. Black R. Changes in property rights and management of high-elevation rangelands in Bhutan: implications for sustainable development of herder communities Mt. Res. Dev. 37 3 2017 353 366
36 Long S.P. A New 1:500,000-scale Geologic Map of Bhutan: a Detailed View of Eastern Himalayan Stratigraphy and Structural Geometry 2010 T43B T2176
37 Wilson A.M. Use of a hydrologic mixing model to examine the roles of meltwater, precipitation and groundwater in the Langtang River basin, Nepal Ann. Glaciol. 57 71 2016 155 168
38 Liu F. Determining hydrologic pathways of streamflow using geochemical tracers in a claypan watershed Hydrol. Process. 34 11 2020 2494 2509
39 Porter V.M. Controls on decadal, annual, and seasonal concentration-discharge relationships in the Sleepers River Research Watershed, Vermont, northeastern United States Hydrol. Process. 36 3 2022 e14559
40 Liu F. Seasonal and interannual variation of streamflow pathways and biogeochemical implications in semi-arid, forested catchments in Valles Caldera, New Mexico Ecohydrology 1 3 2008 239 252
41 Shaw G.D. Groundwater and surface water flow to the merced river, yosemite valley, California: 36Cl and Cl− evidence Water Resour. Res. 50 3 2014 1943 1959
42 Piper A.M. A graphic procedure in the geochemical interpretation of water-analyses Eos, Transactions American Geophysical Union 25 1944 914 928
43 Ravikumar P. Somashekar R.K. Principal component analysis and hydrochemical facies characterization to evaluate groundwater quality in Varahi river basin, Karnataka state, India Appl. Water Sci. 7 2 2017 745 755
44 Liu F. Hunsaker C. Bales R.C. Controls of streamflow generation in small catchments across the snow–rain transition in the Southern Sierra Nevada, California Hydrol. Process. 27 14 2013 1959 1972
45 Taucare M. Groundwater resources and recharge processes in the western andean Front of Central Chile Sci. Total Environ. 722 2020 137824
46 Somers L.D. McKenzie J.M. A review of groundwater in high mountain environments WIREs Water 7 6 2020 e1475
47 Frisbee M.D. Unraveling the mysteries of the large watershed black box: implications for the streamflow response to climate and landscape perturbations Geophys. Res. Lett. 39 1 2012
48 Uchida T. Analysis of flowpath dynamics in a steep unchannelled hollow in the Tanakami Mountains of Japan Hydrol. Process. 17 2 2003 417 430
49 Uchida T. McDonnell J.J. Asano Y. Functional intercomparison of hillslopes and small catchments by examining water source, flowpath and mean residence time J. Hydrol. 327 3 2006 627 642
50 Thakur N. Assessment of recharge source to springs in upper Beas basin of Kullu region, Himachal Pradesh, India using isotopic signatures J. Radioanal. Nucl. Chem. 323 3 2020 1217 1225
51 Kumar M. A comparative evaluation of groundwater suitability for irrigation and drinking purposes in two intensively cultivated districts of Punjab, India Environmental Geology 53 3 2007 553 574
52 Warix S.R. Groundwater geochemistry and flow in the spring mountains, NV: implications for the death valley regional flow system J. Hydrol. 580 2020 124313
53 Jeelani G. Geochemical characterization of surface water and spring water in SE Kashmir Valley, western Himalaya: implications to water–rock interaction J. Earth Syst. Sci. 120 5 2011 921 932
54 Jeelani G. Distinguishing and estimating recharge to karst springs in snow and glacier dominated mountainous basins of the western Himalaya, India J. Hydrol. 550 2017 239 252
55 Bhat N.A. Jeelani G.H. Delineation of the recharge areas and distinguishing the sources of karst springs in Bringi watershed, Kashmir Himalayas using hydrochemistry and environmental isotopes J. Earth Syst. Sci. 124 8 2015 1667 1676
56 Asano Y. Uchida T. Ohte N. Residence times and flow paths of water in steep unchannelled catchments, Tanakami, Japan J. Hydrol. 261 1 2002 173 192
57 Lucianetti G. The role of snowmelt on the spatio-temporal variability of spring recharge in a dolomitic mountain group, Italian Alps Water 12 8 2020 2256
58 Asano Y. An increase in specific discharge with catchment area implies that bedrock infiltration feeds large rather than small mountain headwater streams Water Resour. Res. 56 9 2020 e2019WR025658
59 Mishra B. Babel M.S. Tripathi N.K. Analysis of climatic variability and snow cover in the kaligandaki River basin, Himalaya, Nepal Theor. Appl. Climatol. 116 3 2014 681 694
60 Nepal S. Future snow projections in a small basin of the Western Himalaya Sci. Total Environ. 795 2021 148587
61 Jeelani G. Role of snow and glacier melt in controlling river hydrology in Liddar watershed (western Himalaya) under current and future climate Water Resour. Res. 48 12 2012
62 Chhogyel N. Kumar L. Bajgai Y. Spatio-temporal landscape changes and the impacts of climate change in mountainous Bhutan: a case of Punatsang Chhu Basin Remote Sens. Appl.: Society and Environment 18 2020 100307
63 Chiphang N. Bandyopadhyay A. Bhadra A. Assessing the effects of snowmelt dynamics on streamflow and water balance components in an eastern himalayan River basin using SWAT model Environ. Model. Assess. 25 6 2020 861 883
64 Viste E. Sorteberg A. Snowfall in the Himalayas: an uncertain future from a little-known past Cryosphere 9 3 2015 1147 1167
65 Huntington J.L. Niswonger R.G. Role of surface-water and groundwater interactions on projected summertime streamflow in snow dominated regions: an integrated modeling approach Water Resour. Res. 48 11 2012
66 Vijhani A. Assessment of diminishing discharge of springs in Central Himalayan region, India Hydrol. Process. 36 5 2022 e14582
67 Dailey K.R. Streamflow and groundwater response to precipitation variability in a snow-dominated, subalpine headwater catchment Colorado Rocky Mountains, USA 2019
68 Wilson A.M. High asia: the international dynamics of climate change and water security J. Asian Stud. 76 2 2017 457 480
