Surface and subsurface flow of a glacierised catchment in the cold-arid region of Ladakh, Trans-Himalaya

Hydrological assessments of high-altitude catchments in Trans-Himalayan Ladakh are necessary for a better understanding of water availability in the context of irrigated cultivation under conditions of insufficient quantitative information on cryospheric meltwater discharge. In this study, an integrated spatially distributed temperature index model and a coupled surface/subsurface flow model were used to simulate daily, seasonal, and annual surface and subsurface flows to assess the proportion of corresponding source contributors from the Stok catchment. Snow and glacier meltwater discharge secures irrigated agriculture of more than 300 households in this catchment. The models were forced by temperature, precipitation, ice-and snow-covered areas at daily time steps with calibration (2019; 108 days) and validation (2018; 93 days) against the observed discharge. The simulated discharge shows a good agreement with the observed discharge with R 2 and RMSE of 0.8 (p < 0.01) and 0.6 m 3 /s, respectively. The results between 2003 and 2019 show that the snowmelt contribution to the total annual discharge is largest with 65 %, followed by glacier melt and rainfall contributions of approximately 19 % and 16 %, respectively. A reduction in glacierised areas by 4.2 % was observed while snow-covered areas showed high inter-annual variation. Simulated subsurface flow makes up 62 % (mean = 37.2 × 10 6 m 3 ) of the total discharge with less inter-annual variation. The simulation suggests that while surface flow ceases during the winter period and peaks in August, the annualized mean flow amounts to ~23.7 × 10 6 m 3 . More than 50 % of the melt occurs in the summer months of June, July and August, when the intensity of snowmelt, ice melt, and rainfall reach its maximum. The findings of this study on meltwater availability and surface/subsurface flow is important for irrigated agriculture of Stok village on a local scale, and it might also help to better understand socio-hydrological dynamics and situations of water scarcity in the wider cold-arid region of Ladakh.


Introduction
The Himalayan region contains enormous freshwater reservoirs in the form of glaciers, snow, and permafrost feeding several river basins.Meltwater from these cryospheric sources provides an essential water supply to millions of people in the lowlands (Azam et al., 2021;Immerzeel et al., 2010;Nie et al., 2021;Pritchard, 2019) and in the mountain regions (Mukherji et al., 2019;Nüsser et al., 2019a).A large number of studies shows that the ongoing climate change has a massive impact on the cryosphere and regional hydrology (Azam et al., 2018;Bolch, 2019;Bolch et al., 2012;Gardelle et al., 2013;Kääb et al., 2012;Schmidt & Nüsser, 2012).In addition, population growth, socioeconomic development, and urbanisation processes over the past few decades have increased water demands in the region (Dame et al., 2019;Kaser et al., 2010;Lutz et al., 2016).One of the distinctive hydrological features in the Himalayan region is the rapid transformation of liquid precipitation into runoff.The steep slopes and rugged terrain lead to enhanced surface runoff during heavy precipitation events, often causing flash floods and landslides.However, precipitation of the Indian Summer Monsoon decreases from east to west and rarely reaches the Trans-Himalayan region of Ladakh and the Karakoram (Böhner, 2006;Immerzeel et al., 2009;Singh & Kumar, 1997) where the precipitation is dominated by the mid-latitude westerlies in winter (Banerjee & Dimri, 2019;Chevuturi et al., 2018;Mölg et al., 2014).Thus, the glaciers and snowfields at higher altitudes contribute significantly to river flow ensuring a continuous water supply in downstream regions.In addition, the higher elevations of Ladakh, being part of the upper Indus basin at the western edge of the Himalayas depend mostly on meltwater from cryospheric sources as compared to the catchments of the eastern and central Himalaya (Lutz et al., 2016).
Various hydrological models have been developed and employed to understand the complex water dynamics in the Himalayan region using data from meteorological observations, remote sensing, and groundbased measurements (Tiel et al., 2020).The choice of specific models to simulate and predict surface and subsurface hydrology depends on factors such as data availability, catchment characteristics, specific research objectives, and the level of complexity required to represent hydrological processes accurately (Sidle, 2021).Due to remoteness, extreme weather conditions, difficult terrain, and geopolitical reasons, meteorological and hydrological information from the Himalayan region is generally rare (Winiger et al., 2005) which limits quantitative studies on changes in the regional water resources (Fort, 2015;Singh et al., 2016;Tewari et al., 2017;Tse-ring et al., 2010).Generally, the interconnections between surface and subsurface hydrology in the Himalayas are complex and tightly coupled.For instance, the availability of subsurface water influences the baseflow of rivers during lean flow periods, providing crucial support to river ecosystems and downstream communities.Conversely, the presence of surface water bodies and their interactions with the subsurface influence groundwater recharge rates.Therefore, a detailed exploration of surface and subsurface flows from a small high-altitude catchment may provide insights to better understand the functioning of mountain hydrology.
In the case of the Trans-Himalayan region of Ladakh, characterized by cold-arid conditions, the socio-hydrological system largely depends on meltwater discharge from the cryosphere with regular water scarcity in spring (Nüsser et al., 2012).The ongoing climate warming together with a sharp increase in the tourism footprint associated with high water demand have pushed the region to shift towards alternative sources of water supply, including groundwater extraction (Dame et al., 2019;Müller, 2022).On the other hand, the construction of ice reservoirs in central Ladakh demonstrates the importance of adaption strategies to cope with recurrent water scarcity in spring (Nüsser, et al., 2019b).These developments call for broader socio-hydrological analyses, which integrate environmental and socio-economic dynamics (Nüsser, 2017;Sivapalan et al., 2012).Therefore, the present study aims to simulate both the surface and subsurface hydrology of a high-altitude catchment in Ladakh on daily, seasonal, and annual time scales to assess the proportion of corresponding source contributors from the Stok catchment.For this, we used an integrated spatially distributed temperature index model and a coupled surface/subsurface flow model.

Study area
The Stok catchment (34 • 2' N, 77 • 26' E to 33 • 56' N, 77 • 32' E) is located on the southern bank of Indus River, opposite of Leh town, the capital of the Union Territory of Ladakh, India (Fig. 1).It covers an area of about 66 km 2 .The elevation of the SSW-NNE oriented catchment ranges from 3600 m a.s.l. at the Stok village head up to the highest peak Stok Kangri (6145 m a.s.l.).Snow cover occupies almost the entire catchment during winter and melts gradually from spring onwards ( Passang et al., 2022).Above 5200 m a.s.l.seven small glaciers (<1.5 km 2 ) cover an area of ~3.4 km 2 equivalent to ~5.3 % (in 2019) of the catchment area.All of these glaciers face the northern quadrant (NW-N-NE) like most (74 %) glaciers in the Ladakh region (Kamp et al., 2011;Schmidt & Nüsser, 2012, 2017;Soheb et al., 2022).Other cryospheric components include aufeis, which occurs during winter by successive water overflow and freezing on ice-covered surfaces along the valley floors (Brombierstäudl et al., 2021(Brombierstäudl et al., , 2023)), and permafrost.Meltwater from these cryospheric sources feeds more than 300 households of Stok village before joining the Indus River (LAHDC-Leh, 2017).
Ladakh receives most of its precipitation from July-September and January-March.The mean annual precipitation in Leh (3500 m a.s.l.) amounts to 100 mm with high inter-annual variation while the mean monthly temperature ranges from − 7.2 ℃ in January to 17.9 ℃ in July.Most of the total precipitation falls during the monsoon period with frequent torrential events (Chevuturi et al., 2018;Dimri et al., 2017;Thayyen et al., 2013).Huge tracts of Ladakh are barren with scarce dwarf scrub vegetation, and hygrophilous willows and poplars along the meltwater-fed rivers and streams.While irrigated agriculture entirely depends on regular meltwater supply, the cryosphere is essential for livelihood security and local economy (Nüsser et al., 2012;Schmidt & Nüsser, 2017;Schmidt et al., 2020;Soheb et al., 2020).

Methods
Hydrological models have advanced with greater complexity, computing power and more detailed understanding of processes (Sidle, 2021).Their applicability depends on data availability, regional topography and other complexities.While physically based models provide a comprehensive representation of underlying hydrological processes, they generally require a large amount of information, typically available in the form of reanalysis datasets with low spatial and temporal resolution.Due to their large spatial coverage, the strength of these datasets (e.g., ERA5, CRU, APHRODITE, etc.) lies on the regional scale.Bias-correction or downscaling of reanalysis datasets require extensive observed datasets spanning various parameters, especially for applications on a catchment scale (Berg et al., 2012).It is important to note that bias correction or downscaling of reanalysis temperature (T) and precipitation (P) with in-situ observations, obtained from a location near the study area (Fig. 1), does not provide additional value to the current research compared to the extrapolation method employed in this study.Given the limitations imposed by the scarcity of additional meteorological data (wind and radiation components, humidity, surface temperature, etc.) and hydrological data for the study area, any effort to improve our understanding of hydrology in these catchments requires models that operate with minimal input and necessitate fewer adjustments.Therefore, this study utilized two models, (i) a fully distributed temperature index model: to determine the melt from ice and snow surfaces; and (ii) coupled surface/subsurface flow model: to generate the discharge flow from the catchment (Fig. 2).The models were forced by daily temperature (T), precipitation (P), snow-covered area (SCA) and glacierised area (GA) at each 90 m grids of hydrologically conditioned HydroSHEDs DEM (Hydrological data and maps based on shuttle elevation derivatives at multiple scales Digital Elevation Model; Lehner et al., 2008) over the Stok catchment (Fig. 2).HydroSHEDs is specifically designed for hydrological applications, incorporating relevant corrections to enhance its quality for such purposes.Although SRTM offers higher spatial resolution, it requires additional processing, especially in narrow valleys, before becoming suitable for hydrological modelling.Consequently, HydroSHEDs was the preferred choice for this study.The models were calibrated (108 days; 16 July to 31 October 2019) and validated (93 days; 22 June to 03 October 2018) against the observed discharge at the catchment outlet.The method was further used to model the catchment hydrology from 2003 to 2019 on seasonal and annual time scales.

Discharge measurement
Discharge measurements were carried out in 2018 and 2019 at the catchment outlet upstream of Stok village at a distance of ~ 10 km from the Stok glacier terminus (Fig. 1, Fig. 3).The measurement site at 3600 m a.s.l. was chosen based on regular accessibility for continuous measurements and relatively lower turbulence in the stream water.Discharge measurements were performed using the widely applied areavelocity method where the surface velocity was measured using wooden floats and the water level was monitored on an installed graduated staff gauge (Kumar et al., 2020;Rana et al., 2020;Wulf et al., 2016) (Fig. 1, Fig. 3).The velocity and water level measurements were carried out twice a day (~07:00 and ~ 16:00 h) from 22 June to 3 October 2018 and from 16 July to 31 October 2019.Each time velocity measurement was repeated 3-5 times and the average of the measurements was taken as the final flow velocity.Cross-section of the stream bed was measured three times (early, mid, and end summer) a year using a dipstick to minimize errors due to changes in the stream bed.There is considerable uncertainty associated with the measurement of flow velocity of water at an unpaved gauging site due to technical limitations and human error (Frenierre & Mark, 2014;Mandal et al., 2020).Therefore, the combined uncertainty in discharge measurements is up to 25 % of the total discharge (Eeckman et al., 2017;Mandal et al., 2020).

Surface and sub-surface discharge
The daily total melt and runoff from the catchment were computed with the help of two models i.e., Model-1: A spatially distributed temperature index model to generate the melt rate at each 90 m grid cells, and Model-2: A coupled surface/subsurface flow model to estimate runoff from the snow cover and glaciers using the following expression: where snow and ice melt were generated by Model-1, and runoff, infiltration and exfiltration were estimated using Model-2.

Temperature index model (Model-1)
To determine the amount of seasonal snow and glacial melt, a combination of ground-based and satellite data was used to force the distributed temperature index model on a daily time step (Hock, 2003(Hock, , 2005;;Kumar et al., 2016;R. Singh et al., 2018;Zhang et al., 2018).The model estimates melt rates from glacierised and non-glacierised areas.Snowmelt, glacier melt, and rainfall contribute from the glacierised area whereas only snowmelt and liquid precipitation contribute from the non-glacierised area.
The temperature index model is based on the strong correlation between air temperature and melting rates.The simplest and most common expression relating daily melt to the temperature index is, where M is the melt (mm) for a given period, DDF is the degree-day factor (mm ⁰C -1 day − 1 ), and T and T m are the mean air temperature and threshold temperature for melt (i.e., 0 ⁰C).The air temperature threshold for melting of snow and ice depends on various factors including the presence of impurities (e.g., debris, sediment, dust) and the impact of local conditions (albedo, wind, radiation, cloud cover, topography).Therefore, it is difficult to determine different air temperature thresholds for each case.Thus, the commonly used threshold of 0 • C was adopted according to the literature (Hock, 2003, Lutz et al., 2014, Azam et al., 2019).3.2.1.1. Temperature andprecipitation. Temperature data (2003-2019) were obtained from the Indian Meteorological Department (IMD) Station at Leh, located approximately 8 km north of Stok village.Daily temperature data from Leh station were extrapolated and distributed to HydroSHEDs DEM 90 m resolution grids using monthly lapse rates (Supplement Table 1) obtained by Thayyen & Dimri (2014).They used daily temperatures from three stations located in Leh valley at different elevations (Leh station at 3500 m a.s.l., South Pullu at 4700 m a.s.l., Phutse glacier at 5600 m a.s.l.).Temperature extrapolation performed on the Leh station was compared to the daily in-situ temperature from Lato station, located ~30 South-east of Stok valley at 5050 m a.s.l.(Fig. 1).The comparison between the extrapolated temperatures with in-situ temperatures at Lato station showed good agreement with R 2 = 0.9 at p < 0.01 (Soheb et al., 2020) for the period from 1 July 2018 to 27 September 2019.
The daily precipitation data obtained from the Leh station was observed manually by IMD.In regions like Ladakh where both solid and liquid precipitation occur, IMD measures precipitation with snow gauge, rain gauge and snow poles (IMD, 2023).In high altitude regions, precipitation undercatch is common and significant (Immerzeel et al., 2009).To minimize undercatch in solid precipitation due to wind, the snow gauge employed at Leh station is supplemented with proper wind shield and undergoes regular maintenance.IMD converts solid precipitation to liquid by adding a known amount of hot water, and this conversion is reversed to determine the water equivalent.During heavy snowfall events, snow depth is measured at the snow pole, using a 10:1 ratio for snow to water equivalent.These datasets undergo scrutiny and quality checks at the National Data Centre (NDC) Pune, Govt. of India, following WMO guidelines (Jaswal et al., 2014(Jaswal et al., ,2015)).

Precipitation gradient.
Precipitation distribution over a mountainous region is more complex than temperature because precipitation amounts are not spatially uniform and have strong vertical dependence (Immerzeel et al., 2014;Barry, 2013;Winiger et al., 2005).Here, a gradient of 0.1 m km − 1 (Immerzeel et al., 2014;Soheb et al., 2020) was used to estimate the amount of precipitation.It is crucial to separate solid and liquid forms of precipitation to estimate the contribution to the runoff from these sources.Previous studies (Kayastha et al., 2005;Hagg et al., 2004) have shown that using a temperature threshold to separate the precipitation type is acceptable for modelling hydrological processes in glacierised areas.Therefore, a threshold temperature (T p ) of 1 ℃ (Lejeune et al., 2007) was used.

Degree day factor.
In this study, a degree day factor of 3.1 mm ⁰C -1 day − 1 and 5.9 mm ⁰C -1 day − 1 was used for snow and ice, respectively (Soheb et al., 2020).Degree-day factors for snow and ice were calculated with the help of six ablation stakes, installed on each elevation band of Stok glacier, during the summer (August-September) of 2015, 2016 and 2019 as no significant fresh snowfall has been observed during this period.For each stake location, corresponding cumulated positive degree-days (PDD) were estimated by extrapolating temperature from the Leh station using monthly lapse rates (Thayyen & Dimri, 2014).

Glacierized and snow-covered area.
To incorporate the dynamics of glacier and snow cover changes into the model, manual mapping of glaciers was conducted, and daily snow cover data were obtained from Muhammad and Thapa (2020).Glacier outlines were mapped for three distinct years (2003,2009,2016) as annual changes in glacierized areas were relatively small.The mapping process involved cloud-free ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) Level 1 Precision Terrain Corrected V031, Landsat orthorectified Level-1 T Thematic Mapper (TM), and Planet-Scope imageries, acquired at the end of the ablation period.To minimize errors at high elevations due to perennial snow, the extent of glaciers in the accumulation zone were kept unchanged (Bhambri et al., 2013;Bolch et al., 2010;Garg et al., 2019).The model extrapolates the glacierized areas from the three mapped years for the intervening years.Additionally, the vertical component of glacier change was not considered as the model specifically focuses on generating melt solely from the glacier surface area.Furthermore, the improved daily MODIS (Moderate Resolution Imaging Spectroradiometer) Terra/Aqua snow cover product (Muhammad & Thapa, 2021) from 2003 to 2019 was used to obtain the daily snow cover extent.The Temperature Index Model (Model-1) provides data on snow and ice melt as input for the flow model (Model-2).Details of the data used in this study are presented in Table 1.
The melt contribution from the glacierised area constitutes snow and ice melt, melt induced by heat transfer from rain to snow and ice surfaces, (Kumar et al., 2016;Singh et al., 2018).Firn, which properties lie between snow and glacial ice, was included in the snow-covered area.This inclusion was based on its smaller surface area compared to the seasonal snow cover (roughly around 6-8 % of the study area) predominantly located at higher altitudes (>5500 m a.s.l.).Therefore, its contribution to the overall melt is not anticipated to be substantial to be analysed in this study as a separate component.
Whereas meltwater from the non-glacierised part of the catchment comes only from three sources, i.e., melt contribution from the snow cover, melt induced by the liquid precipitation on the SCA, and direct rainfall.The total melt contribution from the glacierised and nonglacierised areas can be computed using the following equations: where, Q ice and Q snow are the melt contributed by ice and snow (m 3 /day), Q rain on ice and Q rain on snow are the melt induced by rain on ice and snow surface (m 3 /s), Q rain is the direct rainfall contribution (m 3 /s).Moreover, the rainfall-induced melt on snow and ice surfaces can be computed using the following equations: where a ice/snow are the snow and ice cover areas (m 2 ), Q p is the energy supplied to the surface by rain (kJ m − 2 day − 1 ), ρ is the density of water (1000 kg m − 3 ), h f is the latent heat of fusion of water (335 kJ kg − 1 ), B is the thermal quality of surface [B = 0.95-0.97(snow), B = 1 (ice)], C p is the specific heat of water (4.2 kJ kg − 1 ⁰C -1 ), P r is the total rain (mm), and T r and T s are the rain temperature (here mean air temperature was taken as rain temperature) and surface temperature (0 ⁰C for ice and snow surface).

Valley cross-section and fill.
To model the flow of water on the surface and through the subsurface, the sediment or valley fill thickness of the catchment is required.However, direct measurements of sediment thickness are expensive, time-consuming, and spatially not uniform thus such information is always local and requires significant interpolation to get the complete coverage of a region.For the study area of Stok catchment, no direct observations of sediment thickness are available.
Studies have largely used area-volume scaling (Blöthe & Korup, 2013;Straumann & Korup, 2009) or morphometric methods combined with geophysical measurements (Hinderer, 2001;Otto et al., 2009;Schrott et al., 2003) to determine the volume of unconsolidated material in a valley.These methods either do not yield the spatial distribution as in the case of area-volume scaling or are limited to the areas where direct measurements are readily available.Therefore, a relationship between valley width and unconsolidated material thickness was developed following the method used by Mey et al., 2015 (Fig. 4).This approach was based on the geometrical argument that the width of the valley fill correlates positively to the valley fill thickness.It follows a simple geometric projection of hillslopes, along the flow accumulation as the centre line, into the subsurface up to the point of intersection along the valley profile (Fig. 4a).The fill depth at the point of intersection was considered as the maximum depth of the valley fill where the depth decreases with the distance from the centre line (Fig. 4a).The method was applied to 18 different cross-sections at different elevations (Fig. 4b) and different fill width settings.The valley fill width was manually mapped using high-resolution (3 m) PlanetScope images from 2018.The obtained information was further used to develop a regression (Fig. 4d) between valley fill widths and fill maximum depths.A good correlation (R 2 = 0.9 at p < 0.01) was found between fill width and fill depth.This regression was then used to estimate the valley fill depth in the entire catchment at each 90 m grid cells (Fig. 4e).This method determines the valley fill depth based on the slope of each grid cell and the obtained relationship between valley fill width and depth (see section 3.2.2.3.1).A slope threshold of 14 • was set to determine the thickness of a sediment surface, above which there is no fill sediment was considered.

Surface/subsurface runoff modelling.
A coupled surface/subsurface flow model is used to convert total melt, as discussed in subsection 3.2.2.1, into stream discharge at the head of Stok village.The model is based on a shallow flow assumption for the three-dimensional velocity fields in the surface and subsurface compartments.The shallow flow assumption states that i) the vertical velocity component is zero, ii) the horizontal velocity components do not deviate much from their vertical average and iii) there exists a unique vertical position of the water surface at every horizontal position.For the surface, these assumptions allow one to rigorously derive the shallow water equations (also dynamic wave equation) for the vertically averaged horizontal velocity components and water depth from the three-dimensional Navier-Stokes equations with a free surface (Weiyan, 1992).Assuming further that gravitational and bottom friction terms are dominant, the shallow water equations can be simplified to the diffusive wave approximation equation (Ponce and Simons, 1977;Bolster & Saiers, 2002) for the position of the surface water table w s (x, t):  (Yu and Lane, 2006;Bauer et al., 2006).water column at (x, t) and a flat water surface ∇w s = 0 results in a welldefined no flow situation since η > 0 (lake at rest).A comparison of dynamic, diffusive and kinematic wave approximation for surface flow is given by (Ponce and Simons, 1977;Cea et al., 2010).
For groundwater, the starting point is the groundwater flow equation consisting of mass conservation and Darcy's law.Vertical averaging based on the shallow flow assumption replaces Darcy's law with the Dupuit-Forchheimer equation (Bear, 1979): with φ(x) the porosity, w a (x, t) the vertical position of the groundwater table, q a (x, t) the volumetric horizontal water flux in the subsurface, r(x, w s , w a ) is the exchange with the surface water, K( x where r i is infiltration into groundwater from surface water and r e is exfiltration from groundwater to surface water.Exfiltration happens when w a > w s and is simply modelled by a first-order exchange term with a constant L e .For infiltration two conditions are necessary: w s > w a and w s > b s (i.e.water is present at the surface).The fraction (with C = 1) is zero when w s ≤ w a and approaches 1 for w s -w a becoming large.The infiltration flux is driven by w s -b s > 0 and L i is one of the fitting parameters in the calibration.The differential equation model needs to be complemented by initial and boundary conditions.Boundary conditions are taken of Dirichlet type for water height: x ∈ ∂Ω Since the computational domain is chosen much larger than the Stok catchment the boundary conditions do not affect the computed discharge.Initial conditions are explained in the section on calibration.3.2.2.2.1.Numerical simulation.The coupled system of surface and subsurface flow is discretized with a cell-centred finite volume method (FVM) with two-point flux approximation (TPFA) in space (Droniou, 2014) and the (fully) implicit Euler method in time.The FVM computes the water height w T (k) at time t k in cell T of size 90 m.The TPFA of q s from a cell T with centre x T to a neighbouring cell T′ centred at x T' is then computed (for the surface flow) as where the first max ensures non-negativity of the water height, the second max is upwinding and the third max ensures continuity of the flux independent of the bathymetry.The numerical flux for the subsurface is defined accordingly.
In each time step the arising large nonlinear algebraic system is solved with Newton's method, using numerically computed Jacobian systems and is solved in parallel with the BiCGStab method preconditioned with an overlapping Schwarz method (Cai et al., 1998).The implementation is using the Distributed and Unified Numerics Environment (Bastian et al., 2008(Bastian et al., , 2020)).The complete simulation (2003-2019, 17 years) required 45,403 adaptively chosen time steps for 6205 days, i.e., an average time step size of 3.3 h.

Model calibration and uncertainty quantification
3.2.2.3.1.Calibration.The numerical model is calibrated to reproduce the daily discharge measured in the year 2019.Let y j be the measured discharge (in m 3 /s) for day j ∈ I 2019 , with the set I 2019 holds the numbers of the M = 108 days in 2019 from July 16th to October 30th when measurements were taken.Then let m j (k s , b 0 , L i ) be the computed discharge (model output) on day j ∈ I 2019 depending on the three parameters surface conductivity k s , aquifer thickness constant b 0 and infiltration coefficient L i .b 0 is used in the computation of the aquifer thickness b s (x T ) − b a (x T ) = b 0 + 12.3 ⋅ distance T where x T is the position of an aquifer cell and distance T is the distance of the cell from the valley boundary in number of pixels.With the objective function defined as the RMS error in daily discharges the parameter estimation problem is the nonlinear least squares optimization problem The problem has been solved approximately by evaluating the objective function for 480 different combinations of the three parameters in total, with all remaining parameters of the model fixed as given in Table 2.Each calibration run started with zero water height in the surface and subsurface, simulated one year of forcing with the data from 2018 to obtain the initial value for the subsurface, and then simulated the years 2018 and 2019 again with this initial value.The best value was obtained for the following parameters:

Table 2
Parameter values used in the simulation.
Parameter Symbol Value Unit The objective function is, however, quite flat in the vicinity of the minimum and the aquifer thickness b 0 = 67.5 appears to be quite high for the narrowest sections of the valley.Therefore, we have chosen the following parameter set.
Generally, hydraulic conductivity is highly heterogeneous in space and depth.Moreover, the valley fill is composed of varying layers of rocks, gravel, sand, and silt.Thus, it is difficult to determine exact values for the catchment.Therefore, an average value of 0.001 m s − 1 was chosen according to Freeze & Cherry (1979), which is valid for coarse sand and fine gravel.Like hydraulic conductivity, porosity is spatially heterogeneous, but to a much lesser extent.A value of 0.2, characteristic for the lower end for gravel (Geotechdata, 2013), was chosen.Porosity, together with the aquifer thickness, determines the underground storage capacity for water.Since the aquifer thickness parameter b_0 was one of the fitting parameters, one can interpret this as the water capacity being fitted to the measurement data.

Uncertainty estimation.
In this section, we provide a linearized uncertainty estimation of the used parameters.For convenience we organize the M= 108 measurements y = (y 1 ,…, y M ) T and model outputs m (p) = (m 1 (p), …, M (p))T as well as the N = 3 parameters p = (p 1 , p 2 , p 3 ) T = (K s , b 0 , L i ) T in vectors.Assuming the measurements y i are drawn from M statistically independent and normally distributed random variables Y i ~(y i , σi 2 ) with given mean yi and variance σi 2 , the objective function J(p) has the statistical interpretation as the conditional probability or likelihood Where D is a diagonal matrix with the standard deviations D ii = σ i .
Solving the parameter estimation problem can then be interpreted as maximizing the likelihood or, equivalently, as minimizing the negative logarithm of the likelihood yielding: Linearizing the model output with respect to the parameters at a parameter p 0 we obtain m(p) ≈ m(p 0 ) + ∇m(p 0 )(p) with the M × N sensitivity matrix ∇m(p 0 ).For this linearized model, we can solve the least squares minimization problem and obtain In this linearized model the maximum likelihood parameter p* is obtained from the measurements y by an affine linear map.Since the affine linear map of a normally distributed random vector is again normally distributed, the maximum likelihood parameter p* is also normally distributed as with covariance matrix The (in our case) 3 × 3 covariance matrix C describes how the uncertainty in the data given by D, translates into uncertainty of the parameters given by the covariance matrix C. The standard deviations of the measurements, i.e., D, were determined from six independent measurements per day and the resulting covariance matrix is then determined as The diagonal elements represent the variances of the parameters and the corresponding standard deviations are σ ks = 0.481, σ b0 = 6.90, σ Li = 2.18⋅10 − 6 (25) Considering the parameter values k* s = 4, b* 0 = 37.5, L* i = 2.4·10 -5 we conclude that the parameters are reasonably well defined.The offdiagonal entries of the covariance matrix give rise to the correlation √ with the numerical values: ρ ksb0 = − 0.6862, ρ ks Li = 0.8130, ρ boLi = − 0.9534

Air temperature and precipitation
The mean annual air temperature and total precipitation of the Stok catchment showed an inter-annual variation with temperature ranging between -1.4 • C (2019) and 2.6 • C (2016) and precipitation between 375 mm (2007) and 431 mm ( 2006) over the study period of 17 years (2003-2019) (Fig. 5).The coldest and warmest month were found to be January and August with a mean temperature of -12.3 • C and 14.1 • C, respectively.Whereas July received maximum precipitation of 41 mm and November received the least (27 mm).Though the region receives relatively less total annual precipitation while solid precipitation mostly falls during winter.

Glacier change and snow cover variability
The glacierised area of Stok valley decreased from 3.42 km 2 in 2003 to 3.28 km 2 in 2016, revealing an overall change of -0.14 km 2 (-4.22 %) in 13 years with a mean change of -0.1 km 2 /yr (-0.3 % yr − 1 ).The change in total glacierised area found in the present study differs slightly from the results by Schmidt & Nüsser (2017), where reported changes of -0.9 % yr − 1 were reported between 2000 and 2016.These differences are probably due to the use of images from different years and associated snow cover settings and image interpretations.However, the change was found to be less than 1 % yr − 1 in both cases.Individually, the glaciers of the catchment witnessed a reduction in the area over the 13 years but with varying magnitudes ranging between 1 % (0.06 % yr − 1 ) and 14 % (1.1 % yr − 1 ).Though the glaciers (n = 7) in the catchment are all small (<1 km 2 ), a reduction of > 5 % was associated with very small glaciers (n = 4; area < 0.35 km 2 ), while the remaining glaciers (area > 0.6 km 2 ) showed a reduction of < 3.5 % in surface area over the observation period (Table 3).Compared to the larger glaciers of Zanskar in southern Ladakh (Kamp et al., 2011), the glaciers of Stok catchment showed only minor frontal retreat.
The mean annual SCA of the catchment over the 17-year period amounts to 52 % ranging from 35 % in 2016 to 60 % in 2009 and 2019.Seasonal snow coverage varied from year to year with the highest coverage observed in the winter months followed by a gradual decline as the year progresses (Table 4).The mean seasonal snow in the Stok catchment showed a maximum coverage in February (~93 %) while the minimum coverage was found in August (~15 %).The results are in agreement with the regional snow cover analysis of Ladakh where the maximum (~71 %) and minimum (~13 %) snow cover was also found in February and August, respectively (Passang et al., 2022).

Observed and modelled discharge from the catchment
Calibration and validation of the modelled discharge from the Stok catchment were carried out against the observed daily discharge (n = 201 days).For the entire validation period (n = 93 days, 22 June to 03 October 2018), the daily modelled and observed hydrographs are in good agreement (Fig. 6, c) with R 2 and RMSE of 0.75 (p < 0.01) and 0.61 m 3 /s, respectively.However, during the validation period in 2018, the modelled hydrograph fails to match the observed discharge during certain days.For example, the model overestimated the discharge from 21 to 24 July 2018 and underestimated the discharge in the period 28-31 July 2018.As the meteorological data used in the model are based on Leh station.Reduction in the local temperature due to cloud cover above the catchment may not be reflected in the measured temperature, thus providing higher melt rates in the modelled results than in the observed discharge rates.Additionally, the different precipitation patterns due to local wind conditions in both locations (Leh station and Stok catchment) also cause such differences.Slope wind circulation may have resulted in significant underestimation of precipitation measured (Schmidt & Nüsser, 2017;Winiger et al., 2005).The precipitation appears to have a direct representation in the modelled discharge of 2018 where the hydrograph peaks are overestimated (Fig. 6, a), while peaks are not found to be so steep in the observed hydrograph.These limitations cannot be avoided in regions where in-situ measurements are rare (Ali et al., 2022).Furthermore, the distance between the source of melt water and the catchment outlet regulates the flows inducing a bias between the observed and modelled discharge.

Annual and seasonal discharge
The daily simulated surface and subsurface discharge (Supplement video 1; Supplement Fig. 1) in the Stok catchment was further analyzed for annual and seasonal variations to get an in-depth understanding of the discharge patterns.

Annual discharge
The mean simulated total surface and subsurface discharge over the 17 years was 23.66 × 10 6 m 3 and 37.16 × 10 6 m 3 , respectively.The corresponding mean PDD, precipitation and SCA over the same period were 1580 • C, 402 mm, and 52.17 %, respectively (Fig. 7).The catchment-wide mean annual surface discharge shows a large interannual variability between 2003 and 2019, with discharge ranging from 9.64 × 10 6 m 3 in 2018 to 32.76 × 10 6 m 3 in 2006.While relatively less inter-annual variation was observed in simulated subsurface discharge where the annual values ranged between 25.91 × 10 6 m 3 in 2018 to 41.74 × 10 6 m 3 in 2013.The low simulated subsurface discharge in 2017, 2018 and 2019 was the result of less surface runoff due to low SCA and PDD.The surface runoff is one of the main contributors to the subsurface flow other than the assumed permafrost contribution.Due to the existence of significant permafrost areas at higher elevations in Ladakh (Khan et al., 2021;Wani et al., 2020Wani et al., , 2023)), permafrost occurrence can be assumed in Stok catchment at elevations above 5000 m a.s.l., witnessed in some areas by surface deformation resulting from thawing processes (Supplement Fig. 2).Subsurface discharge also depends on the duration of melt, as higher melt rates over a shorter period result in a relatively higher discharge velocity and shorter residence time of surface water.This results in less subsurface flow, whereas gradual and continuous melting of snow and ice feeds the subsurface reservoir.While in the case of low melt rates, the contribution to subsurface reservoirs increases due to the relatively longer residence time of surface water.Such cases are often seen during winter when surface flow ceases to exist, but the simulated subsurface flow remains significant (Figs. 9 & 10).Even in a high-altitude catchment like Stok, melt may occur during daytime, when sub-daily temperatures rise above 0 • C while daily mean temperatures are below the freezing point (Soheb et al., 2018).Subsurface water at lower elevations (unaffected by permafrost) is protected from freezing by thick sediment layers and seasonal snow cover, hence subsurface flow remains throughout the year feeding the stream as baseflow and through springs (Supplement Fig. 3).
The results show that snowmelt is the largest contributor to total discharge with values ranging from 49 to 71 % (mean = 65 %), followed by ice melt and rainfall ranging from 13 to 25 % (mean = 19 %) and 13-29 % (mean = 16 %), respectively (Fig. 8).Between 2016 and 2019, it is evident that the annual contribution of snowmelt was lower than in other observation years, which reduced surface/subsurface discharge, while contributions from ice melt and rainfall remained similar (Figs. 7,  8).In the years 2010, 2017, 2018, and 2019, the contribution from rainfall to discharge was higher than from ice melt (Fig. 8), probably due to frequent cloud bursts and intense rainfall events observed in Ladakh (Banerjee & Dimri, 2019;Thayyen et al., 2013).The decline in discharge after 2015 is attributed to a combination of factors, with snow-covered areas (SCA) and positive degree days (PDD) playing pivotal roles.Although the PDD was notably high in 2016, the SCA was smaller than

Table 3
Change in individual glaciers (Fig. 1) and total glacierised area of Stok catchment.the 17-year average, leading to a decline in discharge.However, relatively large contributions from ice melt and rainfall were also observed.
In 2018, a small SCA and low PDD had a cumulative effect on the resulting low discharge.In contrast, despite the above average SCA, a decline in discharge was primarily due to a low PDD in 2017 and 2019, emphasizing the intricate relationship between snow cover, temperature, and subsequent hydrological impacts.
The simulation results indicate that 62 % of the mean annual melt discharged through the subsurface over the entire observation period of 17 years.However, these estimates vary between different years, where the subsurface discharge ranged between 73 % (2018) and 54 % (2006) of the total discharge.Higher subsurface flows in Stok catchment help to rejuvenate groundwater bodies, while exploitation of these resources is currently increasing as around 90 % of households in Stok village own private borewells (interview with village head in 2020).In the adjoining areas especially in the urban agglomeration of Leh, groundwater resources are currently under extreme stress due to high and increasing water demand (Dolma et al., 2020;Müller et al., 2020).

Seasonal discharge
Seasonal surface discharge in the Stok catchment shows a similar variation in almost all years between 2003 and 2019, where it normally starts in April, peaks in August, and ends in October (Figs. 9 & 10).Such variation in seasonal discharge is quite common in the western Himalayan region (Engelhardt et al., 2017;Kumar et al., 2016;Lutz et al., 2014).For example, in 2010, surface runoff appeared 45 days earlier than in other years because of the enhanced snowmelt contribution in downstream areas of the catchment (Figs. 9 & 10).Whereas, a delay in the onset of surface discharge was observed in May 2016 and 2017 as well as in June 2018 and 2019 when the snow coverage and PDD were below the average.The peak surface discharge in August ranges between 2 m 3 /s (2018) and 6 m 3 /s (2006), due to the additional input from glaciers above 5300 m a.s.l.The maximum seasonal contribution of snowmelt and rainfall varies inter-annually but generally peaks between June and August (Fig. 10).
Simulated subsurface discharge shows a similar trend each year where it peaks a month later than the surface discharge (i.e., September of each year) due to the slow infiltration process of meltwater to subsurface layers.Unlike surface discharge, which ceases during winter (November-March), subsurface discharge occurs throughout the year at varying rates with mean simulated values ranging between 0.45 m 3 /s in March and 2.05 m 3 /s in September.The subsurface flow is well insulated from freezing and receives regular input from melt that occurs during daytime when temperatures are above the freezing point.However, the estimation of sub-daily melt input to the subsurface is beyond the scope of this study as the present model operates on a daily time step.

Comparison between highest and lowest discharge years
To highlight the interannual variability over the observed period from 2003 to 2019, the two most extreme years were compared based on the daily discharge simulations.While 2006 was the year with the highest surface discharge (46 %; equivalent to 32.76 × 10 6 m 3 ), 2018 showed the lowest surface discharge (27 %; equivalent to 9.64 × 10 6 m 3 ).In 2006, the winter snow coverage and possibly the snow thickness were high resulting in long lasting snow coverage above elevations of 5000 m a.s.l., which even persisted until the next winter (Fig. 11 a, b).In contrast, the winter of 2018 was characterized by less snowfall and quickly shrinking snow coverage, which completely melted away by the end of July 2018.In both years the onset of temperature rise in spring (April-May) was around the same time while the SCA in 2018 was already reduced to around 50 % of the catchment area by the end of April resulting in a striking delay of surface discharge.At that time, melt water (if any) was mainly feeding the subsurface due to the high residence time.The contribution of snowmelt amounted to 49 % (18 × 10 6 m 3 ) only, while that of ice-melt and rainfall increased relatively to 22 % (8 × 10 6 m 3 ) and 29 % (10 × 10 6 m 3 ), respectively (Fig. 11c).In contrast, in 2006 the snow cover gradually melted from April until October, continuously contributing to the surface and subsurface flow of Stok catchment.Thus, the contribution of snowmelt to total discharge was 69 % (50 × 10 6 m 3 ) while that of ice-melt and rainfall amounted to 17 % (13 × 10 6 m 3 ) and 14 % (10 × 10 6 m 3 ), respectively.Furthermore, the rainfall contribution was similar in both years, while that of ice was slightly below the 17-year average in 2018.The inter-annual snow cover variability is important for the onset of surface discharge and melt water security in irrigated agriculture.
4.6.Limitation and applicability of this approach on a regional scale Hydrological models are valuable tools to understand processes and behaviour of the hydrological cycle through mathematical representations.Such simplifications cannot capture the complexity of real-world systems entirely leading to potential inaccuracies and uncertainties (Sidle, 2021).Therefore, model calibration and validation of the results are essential components for improving the accuracy and reliability of hydrological models.Our study involves for the first-time simulated surface and subsurface flow components for the Trans-Himalayan region of Ladakh.In general, observed data on surface/subsurface flows are scarce for the Himalayan region (Singh et al., 2016).In the present study, surface flow was calibrated with in-situ measurements at the catchment outlet, however, direct subsurface flow measurements were not available, which may result in higher uncertainty and hamper the model performance.Thus, subsurface flow measurements are crucial and warrant attention in future studies.Subsurface flow information can be obtained by employing direct and indirect approaches.Direct access to the subsurface can be obtained using boreholes and wells for water level recordings.While tracer-based studies (Cascarano et al., 2021;Anderson et al., 2015;Baraer et al., 2015), employing substances like dyes, and geophysical methods (Zuo et al., 2023;Christensen et al., 2020;Cassidy et al., 2014;McClymont et al., 2011) such as electrical resistivity and ground-penetrating radar, contribute collectively to a comprehensive understanding of subsurface hydrology.
Hydrological models operate at specific spatial scales and temporal resolutions, based on data availability and computing ability.Although the present study utilises the best available data, the model lacks the ability to explain local heterogeneities and processes on a finer spatial and temporal scale.One example is the case of aufeis, which is regularly formed along the valley floor of the Stok catchment during winter.Although it covers an area of about 0.9 km 2 between 3650 m and 4700 m a.s.l.(Brombierstäudl et al., 2021), the detection of narrow aufeis fields (average width: 88 m) is not possible using MODIS snow product due to the low spatial resolution.For future hydrological studies it is recommended to include aufeis as a separate cryospheric component.In addition, accurate estimations of parameters representing the characteristics of a hydrological system, are also challenging in Ladakh, due to the lack of data, topographical heterogeneity, and temporal variability.Furthermore, limited process representation also hampers the model performance at different scales.The present method focuses on major hydrological components while other processes, such as evaporation, transpiration and sublimation have been neglected due to unavailability of data and modelling requirements.Since the study area of Stok catchment is a north-facing valley in the cold-arid region of Ladakh, evaporation and transpiration rates are expected to be very low and nonsignificant.While sublimation is expected to be high in the catchment due to the arid conditions, it was not estimated as the present study focused on simulated surface and subsurface flows rather than on the hydrological budget or hydrological balance.
Despite these limitations, the study not only enhances our understanding of the surface and subsurface flow in Stok catchment, but it also holds potential for creating extreme event scenarios.While our study does not specifically address the interactions between extreme climatic events and hydrology, the presented findings and approach could serve as valuable tools to simulate and understand the impact of extreme events on meltwater dynamics and surface flow.This information is crucial for identifying vulnerable areas and potential threats in the future, especially in regions where detailed meteorological and hydrological data are lacking.Such comprehensive analyses are essential for socio-hydrological studies on local and regional scales, integrating hydrological and environmental dynamics with socio-economic development paths to improve sustainable water management, as highlighted by Nüsser (2017).In the case of Stok catchment, examples include assessing water availability for irrigated crop cultivation and addressing ongoing groundwater depletion linked to the observed increase in bore wells in Stok village.
In future hydrological investigations of Ladakh, where meltwater significantly influences streamflow, an emphasis on physically based models and reanalysis products is essential.Further, the complexity of predicting hydrology in these areas arises from factors such as snow redistribution, ablation, infiltration into frozen soils, bidirectional phase changes, and episodic flow processes.Models like CRHM (Cold Regions Hydrological Modelling; Pomeroy et al., 2022), VIC (Variable Infiltration Capacity; Liang et al., 1994), SWAT (Soil and Water Assessment Tool) or like need to be explored.While CRHM, because of its comprehensive capabilities in simulating various hydrological processes of a cold region, may possibly be well-suited for the particular environmental conditions of Ladakh.It considers essential factors like radiation,    sublimation, wind, energy for melt, hill shades, frozen soil, and permafrost which is obtained through reanalysis products (e.g.ERA5, ERA-L, HAR) and remote sensing (Aubry-Wake and Pomeroy, 2023;Pradhananga and Pomeroy, 2022).The reanalysis products as evidenced by their successful application in the Himalayas (e.g., Bhattacharya et al., 2020;Bannister et al., 2019;Li et al., 2018), when bias-corrected and downscaled, offer promising solutions for data scarcity in different regions of High Mountain Asia, including Ladakh.Despite their global availability and openness, uncertainties in correction processes and limitations in capturing local variations should be considered.Such an exploration holds significant promise for expanding our understanding of the total hydrology and climate of Ladakh including diverse ecological settings.
Moreover, observations from Chandra-Bhaga and Chhota-Shigri catchments of the western Himalayan region (Azam et al., 2019;Srivastava and Azam Mohd, 2022;Supplement Fig. 4) correspond with the present study regarding the dominance of snowmelt and timing of peak discharge between June and August.The annual rainfall contribution differs between 16 % in Stok, 10 % in Chhota-Shigri and 22 % in Chandra-Bhaga catchments, due to the higher impact of the monsoon in the Chandra-Bhaga catchment.

Conclusion
Simulated surface and subsurface flow components were estimated for the Trans-Himalayan region of Ladakh for the first time.Such information is of utmost importance for a sustainable water management and for an improved understanding of socio-hydrological settings on the local scale.This study simulates annual and seasonal surface and subsurface flows in the Stok catchment using an integrated fully distributed temperature index model and coupled surface/subsurface flow model.The models were forced by temperature, liquid precipitation, SCA and glacierised area on a 90 m hydrologically conditioned HydroSHEDs DEM on a daily basis.They were calibrated (108 days) and validated (93 days) using discharge measurements from the catchment outlet.To estimate the infiltration rates and subsurface flow, the sediment thickness of the valley fill was estimated using a linear projection of the valley slope angles.
The study shows that snow melt is the largest contributor (65 %) to the total annual discharge followed by ice melt (19 %) and rainfall (16 %).Simulated subsurface flow indicates that most (62 %) of the available meltwater infiltrates to the ground.More than 50 % of the melt occurs during summer (June-August).Surface flow usually starts in April, peaks in August, and ends in October, while the subsurface flow remains active throughout the year.The present approach can be used in other catchments of the Trans-Himalayan region of Ladakh and in similar cold-arid environments with limited data on water sources and unknown relative contribution of different cryospheric components to

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.

Fig. 1 .
Fig. 1.Study area: Location of Stok catchment and the meteorological stations (Leh, South Pullu, Phutse and Lato) (left); and detailed view of Stok catchment depicting glaciers, villages, and the location of the discharge measurement site (catchment outlet) (right).

Fig. 2 .
Fig. 2. (a) A conceptual diagram of hydrological processes in a mountain catchment of the Trans-Himalayan region of Ladakh and (b) A flow diagram detailing the methodology used in this study.

Fig. 3 .
Fig. 3. Field photographs containing (a) hydrometric station showing the staff gauges and flow direction, (b, c) autumn and winter photographs overlooking upstream, and (d, e) autumn and winter photographs overlooking downstream from the discharge measurement site (a).
Fig. 4. Schematic illustration of the estimation of valley fill thickness; (a) The idealized valley cross section using the geometry of the valley; (b) The 18 valley profiles used in the method on 90 m HydroSHEDs DEM; (c) An example of a profile on high-resolution PlanetScope imagery of 2018; (d) regression between the valley fill width and the sediment thickness; (e) the estimated sediment thickness after applying the method to the entire catchment.
) is the hydraulic conductivity and b a (x) is the position of the impermeable bedrock (groundwater bathymetry).Note the similarity of surface and subsurface flux when setting α = 1 and η = 1.The thickness of the aquifer b a (x) − b s (x) is determined as discussed in subsection 3.2.3.The min function in the storage term ensures that storage and height of the water column are decoupled when w a (x, t) > b s (x), i.e. the aquifer is confined, which can happen below a lake or stream.The exchange between surface and groundwater is modelled byr(x, w s , w a ) = r i (x, w s , w a ) − r e (x, w s , w a )(11)r e (x, w s , w a ) = L e max(w a − w s , 0) (12)

Fig. 5 .
Fig. 5. Mean (a) annual and (b) monthly temperature and total precipitation of the Stok catchment.

Fig. 6 .
Fig. 6.Observed and modelled surface discharge of 2018 (a) and 2019 (b), and associated simulated subsurface discharge, PDD, precipitation, and SCA.Observed vs modelled discharge for 93 days in 2018 (c) and 108 days in 2019 (d) of validation and calibration period, respectively.Snow and rain bars represent the amount of precipitation at the centre of Stok catchment (~4800 m a.s.l.).

Fig. 7 .
Fig. 7. Annual simulated total surface and subsurface discharge from the Stok catchment, and the PDD, precipitation (snow and rain) and SCA between 2003 and 2019.

Fig. 8 .
Fig. 8. Annualised contribution from snowmelt, ice melt and rainfall to total discharge in the Stok catchment between 2003 and 2019.M.Soheb et al.

Fig. 10 .
Fig. 10.Total seasonal contribution from snowmelt, ice melt and rainfall to simulated surface and subsurface discharge of the Stok catchment from 2003 to 2019.

Fig. 11 .
Fig. 11.Comparison between the years with highest (2006) and lowest (2018) annual discharge.a) daily simulated surface and subsurface discharge with associated snow-covered area (SCA), precipitation (P) and degree days; b) spatiotemporal changes in SCA; c) contribution from snowmelt, ice melt and rainfall to simulated surface and subsurface discharge.

Table 1
Details of the data used in the study.
Own, in-situ discharge measurements M.Soheb et al.

Table 4
Monthly mean snow coverage (SCA in %) of the Stok catchment derived from the improved MODIS snow cover dataset of Muhammad &