Revealing the influence of topography and vegetation on hydrological processes using a stepwise modelling approach in cold alpine basins of the Mongolian Plateau

Topography and vegetation are critical factors influencing catchment hydrology; however, their individual contributions are often underestimated in hydrological models. This limitation is particularly evident in cold, mountainous regions such as the Mongolian Plateau, where observational data are sparse. To address this, we employed a stepwise, top-down modelling strategy based on a flexible modelling framework to systematically assess the influence of topography and vegetation on hydrological processes in the Bogd Uliastai and Zavkhan Guulin river basins. Beginning with a lumped model (FLEX^L), we successively integrated snow processes (FLEX^L-S), topographic distribution (FLEX^D), and finally, a landscape-based parameterization accounting for vegetation heterogeneity (FLEX^T). Both FLEX^D and FLEX^T outperformed the lumped models in simulating runoff and snow water equivalent (SWE). Interestingly, FLEX^T showed similar performance to FLEX^D – likely due to limited vegetation heterogeneity – it offers more physically realistic parameterization by explicitly representing landscape units, suggesting its potential in more complex basins. The ratio of snowmelt runoff to streamflow was quantified as 23.6 % ± 0.7 % and 15.9 % ± 1.3 % in the Bogd Uliastai and Zavkhan Guulin river basins, respectively, with peaks in spring and a clear increase with elevation. At high elevations, runoff is primarily snowmelt-driven, resulting in delayed and gradual runoff, whereas lower elevations dominated by rainfall generate rapid runoff. Controlled by distinct dominant hydrological mechanisms, different landscape units contribute unequally to streamflow. This study underscores the pivotal roles of topography and vegetation in runoff generation and demonstrates the effectiveness of a stepwise modelling framework for improving hydrological understanding in cryospheric and data-scarce regions.

1 Introduction

Understanding and accurately simulating hydrological processes are fundamental for elucidating basin hydrological patterns and supporting water resource management and ecological protection under global environmental change (Oki and Kanae, 2006; Gomes et al., 2023). Topography and vegetation play essential roles as drivers of hydrological processes, influencing precipitation, interception (Dwarakish and Ganasri, 2015), snowmelt (Hammond et al., 2019), evaporation (Jiao et al., 2017), and runoff generation (Qin et al., 2025). Their combined influence underpins landscape organization, forms the basis for Hydrological Response Units (HRUs), and shapes spatial heterogeneity and dominant hydrological mechanisms (Savenije, 2010; Sivapalan, 2009). Limited observational data in cold-arid regions often results in oversimplified hydrological models, emphasizing the need to represent landscape controls more explicitly (Ragettli et al., 2014; Tarasova et al., 2016).

Topography governs water redistribution across landscapes by shaping soil moisture patterns, modulating precipitation and evaporation, and controlling runoff pathways (Wicki et al., 2023). In mountainous basins, strong relief introduces substantial challenges for hydrological modelling and increase predictive uncertainty (Seibert and Mcdonnell, 2002). Steeper slopes typically lead to more rapid runoff, while gentler slopes promote greater infiltration and moisture retention, thereby affecting the spatiotemporal distribution of water resources (Ye et al., 2023). In addition, topographic features such as slope, aspect, and elevation control snow accumulation and melt, resulting in highly heterogeneous hydrological responses across basins (Broxton et al., 2020).

Vegetation also regulates hydrological processes by intercepting snow and rainfall, facilitating water infiltration and storage in the root zone, and absorbing and transpiring water through its roots. Canopies reduce effective precipitation reaching the ground and can mitigate surface runoff (Cheng et al., 2020), while root systems shape soil moisture dynamics and water redistribution (Volpe et al., 2013; Luo et al., 2022). Consequently, different vegetation types (e.g., forests, grassland, and abandoned farmland) often exhibit distinct hydrological behaviors (Chen et al., 2023). In cold mountainous regions, vegetation further alters snow processes through shading and moderating wind redistribution, thereby influencing the meltwater timing and magnitude (Sun et al., 2022). Over long timescales, snow-vegetation interactions can modulate runoff sensitivity to precipitation, evapotranspiration, and vegetation dynamics, especially in regions such as Central Asia (Feng et al., 2025) and the Tibetan Plateau (Ni et al., 2025).

Despite the recognized importance of topography and vegetation, their influences remain insufficiently represented in many hydrological modelling frameworks (Stephens et al., 2021). This gap is especially critical in cryospheric regions, where snow accumulation and melt dominate runoff generation and are strongly shaped by terrain and vegetation controls (Zhong et al., 2021; Dharmadasa et al., 2023; Immerzeel et al., 2010). Quantifying and integrating these controls into hydrological models is therefore essential for advancing cold-region hydrology and improving model realism.

Existing hydrological models often struggle to adequately capture the complexities introduced by topography and vegetation heterogeneity. Early lumped models, relying on basin-averaged inputs, oversimplify spatial heterogeneity within catchments (Beven, 2012). Distributed hydrological models allow spatially explicit simulations but are highly dependent on the availability and quality of input data (Fenicia et al., 2016), which remains a significant limitation in cold, high-mountain regions. Remote sensing has become an invaluable tool for providing spatially continuous data on topography (e.g., Digital Elevation Models, DEMs), vegetation (e.g., Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI)) (Xiong et al., 2023), and snow dynamics (e.g., SWE) (Duethmann et al., 2014), enabling a more realistic representation of landscape-dependent controls in hydrological models (Gao et al., 2014).

In the absence of direct measurements of individual hydrological processes, the top-down modelling approach offers a powerful and practical means of exploring the internal dynamics of basin behavior (Fenicia et al., 2008b). Originally proposed by Klemeš (Klemeš, 1983) and later reformulated by Sivapalan et al. (Sivapalan et al., 2003), the top-down approach is rooted in a deductive philosophy that infers the underlying “causes” from the overall observed “effect” of a system. In hydrological modelling, this method begins with a simple structure that is progressively refined to address limitations in reproducing observed catchment behavior (Fenicia et al., 2008a).

The Bogd Uliastai and Zavkhan Guulin river basins, in the headwaters of the Zavkhan River on the western Mongolian Plateau, exhibit a cold, arid continental climate with sparse vegetation (Baasanmunkh et al., 2019). Strong cryospheric influences control hydrological processes, making these basins key areas for regional ecological security and downstream water-resource resilience. However, monitoring networks are sparse and long-term observational data are limited, leaving runoff generation mechanisms and their interactions with topography, vegetation and the cryosphere poorly understood (Dorjsuren et al., 2024).

To address these challenges, we employ a top-down modelling framework that begins with a lumped model to assess runoff dynamics and progressively incorporates distributed representations of snowmelt, topography, and vegetation. This framework allows assessment of hydrological responses across different landscape units and addresses three research questions:

1. How can runoff be effectively simulated in data-scarce, cold mountainous regions using a top-down modelling approach?

2. How can the contribution of snowmelt runoff to streamflow be quantified using a landscape-based hydrological model?

3. How do topography and vegetation influence runoff generation processes?

The study aligns with the IAHS HELPING Decade (Hydrology Engaging Local People IN one Global world, 2023–2032), which calls for interdisciplinary research in data-scarce and ecologically vulnerable regions. By developing a modelling framework for cold basins with limited observational data and improving understanding of hydrological processes on the Mongolian Plateau, this work provides an important case study contributing to the goals of the HELPING initiative.

2 Study site

2.1 Bogd Uliastai river basin

The Bogd Uliastai river basin (47°30^′–48°10^′ N, 96°45^′–97°45^′ E) is located in the northern part of the Zavkhan river headwaters, along the southern foothills of the central Khangaï Mountains in Mongolia (Fig. 1). The basin spans an area of 1610 km² and is predominantly mountainous, with elevations ranging from 1753 to 3972 m a.s.l. The region receives an average annual precipitation of approximately 200 mm, with more than 80 % of rainfall occurring between June and September. The average annual temperature is −1 °C, while winter temperatures frequently fall below −30 °C, reflecting a typical alpine climate. Runoff displays strong seasonal variability, with distinct peaks during the spring and summer and almost no flow in winter, resulting in extreme hydrological conditions (Dorjsuren et al., 2024). The vegetation exhibits clear altitudinal zonation: alpine meadows and tundra, dominated by mosses and lichens, prevail at higher elevations, whereas needlegrass steppe and low shrublands are common in mid and low elevation areas (Baasanmunkh et al., 2019).

Figure 1 Location, landscape (a) and topography (b) of the Bogd Uliastai and Zavkhan Guulin river basins on the Mongolian Plateau.

2.2 Zavkhan Guulin river basin

The Zavkhan Guulin river basin (46°30^′–47°50^′ N, 96°45^′–97°00^′ E), located in the central and southern parts of Zavkhan Province, lies within the transitional zone of the southern Khangaï Mountains (Fig. 1). The basin covers an area of approximately 12258 km² and is predominantly composed of low mountains and hills, with elevations ranging from 1785 to 3980 m a.s.l. The basin's annual average precipitation is about 160 mm, with most precipitation concentrated in summer. The annual average temperature is approximately −3 °C, with summer temperatures exceeding 20 °C and winter temperatures dropping as low as −50 °C, characteristic of a temperate continental climate (Dorjsuren et al., 2023). Vegetation in the region is sparse, primarily dominated by drought-tolerant Artemisia species, with scattered distributions of grass and shrubs. At higher elevations, the landscape is characterized by alpine meadows, exposed rock surfaces, and cold desert environments. Soils are nutrient-poor, and the ecological environment is fragile, facing severe challenges such as soil erosion (Baasanmunkh et al., 2019).

3 Data

3.1 Data set

Hydrometeorological data: Daily precipitation, runoff, and temperature data for the Bogd Uliastai river basin (2007–2015) and the Zavkhan Guulin river basin (2000–2020) were obtained from the Information and Research Institute of Meteorology, Hydrology, and Environment (IRIMHE) via its official website (http://irimhe.namem.gov.mn, last access: 26 March 2026). For each basin, one meteorological station and one hydrological station served as the primary sources of observational data. The Arctic Snow Water Equivalent (SWE) Grid Dataset (2003–2016) was obtained from National Tibetan Plateau/Third Pole Environment Data Center (https://cstr.cn/18406.11.Snow.tpdc.271556, last access: 26 March 2026). The SWE product has a daily temporal resolution and a spatial resolution of 10 km, covering latitudes from 45 to 90° N and longitudes from 180° W to 180° E.

Topographic data: The Shuttle Radar Topography Mission Digital Elevation Model (SRTM-DEM), with a spatial resolution of 90 m, was acquired from the official website of the CGIAR Consortium for Spatial Information (CGIAR-CSI) (http://srtm.csi.cgiar.org, last access: 26 March 2026).

Land cover data: The Sentinel-2 10 m Land Use/Land Cover was accessed via ESRI's official platform (https://livingatlas.arcgis.com/landcover/, last access: 26 March 2026).

NDVI data: The normalized difference vegetation index (NDVI) data (2013–2020) were derived from the Landsat 8 Operational Land Imager (OLI) Level-2 surface reflectance products. NDVI was calculated as (NIR − Red) $/$ (NIR + Red) using bands 5 (NIR) and 4 (Red). The dataset has a spatial resolution of 30 m and a temporal resolution of 16 days. Landsat data were obtained from the United States Geological Survey (USGS) EarthExplorer platform (https://earthexplorer.usgs.gov/, last access: 26 March 2026).

3.2 Distribution of forcing data

Mountainous terrain is complex, and meteorological stations are typically located at lower elevations. Directly using point-based measurement in basin-scale simulations without accounting for elevation effects can introduce biases (Klemeš, 1989). In cold mountainous regions, higher elevations typically experience lower temperatures and greater precipitation, often in the form of snow (Lundquist et al., 2010; Stahl et al., 2006). In the study, the FLEX^D and FLEX^T models divide catchment into elevation bands and adjusts temperature and precipitation for each band using a precipitation increase rate and temperature lapse rate (Fig. 2). This distributed input approach effectively mitigates simulation bias by better capturing altitudinal variability in meteorological conditions. In this study, due to the remoteness of the region and the sparse distribution of meteorological stations, available ground observations were limited. Satellite and reanalysis products (e.g., ERA5-Land) exhibit notable biases over complex terrain and fail to capture local climatic variability. We therefore used the nearby Tsagan-Turutuin-gol catchment as a reference, which also originates in the Khangaï Mountains and shares similar topographic and climatic characteristics with the study basins (Klimek and Starkel, 1980). The model employed a precipitation increase rate of 4.2 % per 100 m and a temperature lapse rate of 0.6 °C per 100 m .

Figure 2 Area of different elevation and landscape in Bogd Uliastai and Zavkhan Guulin river basins.

4 Modelling approach

4.1 Model description

To assess the impact of topography and vegetation on hydrological processes, this study designed and tested four conceptual models with increasing complexity: FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T. The model structure and variables are shown in Fig. 3 and Table 1, while the water balance and constitutive equations are summarized in Table 2.

Figure 3 Stepwise modelling and the model structure of four models. (a) FLEX^L is a lumped model without snow module; (b) FLEX^L-S is a lumped model with snow module, and FLEX^D is a semi-distributed model with the same structure as FLEX^L-S. (c) FLEX^T is a landscape-driven semi-distributed model.

Table 1 The variables of four models. In FLEX^T model, variables associated with various landscape categories are differentiated using specific suffixes, e.g., E_i,F, represent the interception evaporation from forest.

Table 2 The water balance and constitutive equations used in four models. Note: FLEX^L model lacks the snow module, resulting in different water balance and structural equations compared to other models. For FLEX^T model, in Eqs. (7) and (9), the S_i and S_i,max represent the interception reservoir and its interception capacity, respectively, for different landscape units, including forest (S_i,maxF), grassland (S_i,maxG) and riparian area (S_i,maxR). No interception reservoir is defined for bare soil/rock areas. Similarly, the S_u and S_u,max denote the root zone reservoirs and their corresponding storage capacities for different landscapes, while K_f represents the recession coefficient of the fast response reservoir for each landscape type.

4.1.1 FLEX^L

FLEX^L is a lumped conceptual hydrological model composed of four reservoirs (Fig. 3a): an interception reservoir (S_i), a root zone reservoir (S_u), a fast response reservoir (S_f), and a slow response reservoir (S_s). A lag function is used to represent the lag time from storm to peak flow (T_lag). FLEX^L includes a total of 8 free calibration parameters (Table 3).

The interception reservoir was designed to simulate the process of precipitation interception by vegetation canopies or the ground surface (Eq. 1). Interception evaporation (E_i) was calculated by potential evaporation (E_p) and S_i, considering the interception storage capacity (S_i,max) (Eq. 2). When precipitation (P) exceeds S_i,max, the excess precipitation is routed as effective precipitation (P_e) into the root zone reservoir (Eq. 3).

In the root zone reservoir, actual evaporation (E_a) was estimated based on E_p−E_i and root zone soil moisture ($S_{\mathrm{u}}/S_{\mathrm{u},\max }$) (Eq. 13). The parameter C_e represents the threshold value controlling evaporation from the root zone soil moisture, and S_u,max is root zone storage capacity. The water retention curve from the Xin'anjiang model was used to partition P_e into stored water in S_u and runoff generated from the root zone (R_u) (Eqs. 14 and 15) (Zhao, 1992).

In the response reservoir, a splitter D was applied to divide the R_u into two fluxes (R_f andR_s) (Eqs. 16 and 17), and Eqs. (18) and (19) were used to describe the lag time between storm and peak flow. $R_{\mathrm{f}}(t-i+$1) represents the fast runoff generated in the root zone at time $t-i+\mathrm{1}$, T_lag represents the time lag between the storm and fast runoff generation. c(i) is the weight of the flow in i−1 d before and R_fl(t) is the discharge into the fast response reservoir after convolution. We used two linear reservoirs to represent the response process of subsurface storm flow (Q_f) and groundwater flow (Q_s) (Eqs. 21 and 23).

4.1.2 FLEX^L-S

FLEX^L-S builds upon the FLEX^L model by incorporating a snow reservoir (S_w) to simulate the snow accumulation and melt processes (Fig. 3b). When the daily air temperature exceeds the threshold temperature (T_t) and there is no snowpack (typically in summer), the interception process governs the initial partitioning of precipitation (Eq. 7). In contrast, when the daily mean temperature is below T_t (normally occurs in winter), precipitation is stored as snow (Eq. 5). When there is snowpack and the daily air temperature is above T_t (normally prevailing in early spring and early autumn), effective precipitation (P_e) is equal to the sum of effective rainfall after interception (P_tf) and snowmelt (M) (Eq. 11). M was calculated by the snow degree day factor (F_dd) and the threshold temperature for melting (T_tm) (Eq. 6) (Gao et al., 2017). In this study, T_tm was set to the same value as T_t. It is important to note that meltwater is conceptualized as directly infiltrating into the soil, thereby bypassing the interception reservoir.

4.1.3 FLEX^D

FLEX^D is a semi-distributed model with the same structure and parameters as FLEX^L-S (Fig. 3b). Using DEM data, the Bogd Uliastai river basin was divided into 45 elevation bands with 50 m interval, while the Zavkhan Guulin river basin was divided into 44 elevation bands as shown in Fig. 2. The FLEX^D model was operated with semi-distributed input data (see Sect. 3.2), ensuring the integration of spatial variability into the model's processes.

4.1.4 FLEX^T

The FLEX^T model classified the Bogd Uliastai river basin into four landscape elements – bare soil/rock, forest, grassland, and riparian area – based on vegetation characteristics. In contrast, the Zavkhan Guulin river basin, which has no forest, was categorized into three landscape elements. By integrating both landscape types and elevation bands, the Bogd Uliastai river basin was further subdivided into 132 HRUs, while the Zavkhan Guulin river basin consisted of 117 HRUs (Fig. 2).

The FLEX^T model's structure comprised four parallel components, representing the distinct hydrological functions associated with landscape elements (Gao et al., 2014). To capture the diverse rainfall-runoff processes in different landscape types, landscape elements were assigned model structures (Fig. 3c) and parameters (Table 3) tailored to their dominant hydrological processes.

For bare soil and rock areas, vegetation is absent; therefore, no interception occurs. Due to the low infiltration capacity of these surfaces, Hortonian overland flow (R_hof,B) dominates runoff generation when daily effective precipitation (P_e,B) exceeds the threshold parameter (P_t) (Eq. 24):

$$\begin{array}{}\text{(1)}& R_{\mathrm{hof},\mathrm{B}}=max\left(P_{\mathrm{e},\mathrm{B}}-P_{\mathrm{t}},\mathrm{0}\right)\end{array}$$

Saturation overland flow (R_sof,B), caused by the limited storage capacity of shallow soils, occurs when the unsaturated zone soil moisture reaches its maximum storage capacity (S_u,maxB). Deep percolation into groundwater (R_p,B) is governed by relative soil moisture ($S_{\mathrm{u},\mathrm{B}}/S_{\mathrm{u},\mathrm{maxB}}$) and the maximum percolation (P_max; Eq. 25). Surface runoff is further partitioned into re-infiltrating water (R_r,B) over high-permeability debris slopes and water routed directly to the channel (R_f,B) through a separator (D_B).

$$\begin{array}{}\text{(2)}& R_{\mathrm{p},\mathrm{B}}=P_{\max }{\displaystyle \frac{S_{\mathrm{u},\mathrm{B}}}{S_{\mathrm{u},\mathrm{maxB}}}}\end{array}$$

Forest and grassland hillslopes follow the same FLEX^D model structure, with differences represented through parameters governing interception and root zone processes.

For the riparian area, which is prone to saturation due to its location, we also constrained a shallower S_u,maxR, with the effect of capillary rise (C_r) taken into account. C_r is represented by a parameter (C_r,max) indicating a constant amount of capillary rise. The lag time from storm to peak flow was not considered in riparian area.

Interception thresholds were parameterized to be highest for forest (S_i,maxF), followed by grassland (S_i,maxG) and riparian area (S_i,maxR), reflecting differences in canopy interception capacity.

The root zone storage capacity was assigned deepest in forests (S_u,maxF), shallower in grasslands (S_u,maxG), and smallest in riparian areas (S_u,maxR). Bare soil and rock (S_u,maxB) were represented as belonging to the unsaturated zone rather than the root zone, with storage capacity similar to that of riparian areas. Response timescales were parameterized to be longest for groundwater (K_s), shortest for saturation-excess runoff in riparian area (K_fR) and surface runoff in bare soil/rock (K_fB), with subsurface flow from forest and grassland (K_fF and K_fG) parameterized to operate at intermediate timescales.

Table 3 Uniform prior parameter distributions of four models.

4.2 The ratio of snowmelt runoff to streamflow

This study quantifies the ratio of snowmelt runoff to streamflow using the FLEX^T model. The model assumes that snowmelt and rainfall mix rapidly and completely upon entering the model's conceptual reservoirs, thereby altering its internal composition ratios. The composition ratio of the water exiting the reservoir is identical to that within the reservoir (Liu et al., 2023). The method enables the tracking of ratio of snowmelt runoff to streamflow (f_Q,snow) at each time step by the following equation:

$$\begin{array}{}\text{(3)}& {\displaystyle }{f}_{Q,\mathrm{snow}}={\displaystyle \frac{Q_{\mathrm{M}}}{Q}}={\displaystyle \frac{Q_{\mathrm{f},\mathrm{M}}+Q_{\mathrm{s},\mathrm{M}}}{Q}}\text{(4)}& {\displaystyle }{Q}_{\mathrm{f},\mathrm{M}}={\displaystyle \frac{\left(\frac{M}{P_{\mathrm{e}}}\right)S_{\mathrm{f}}}{K_{\mathrm{f}}}}\text{(5)}& {\displaystyle }{Q}_{\mathrm{s},\mathrm{M}}={\displaystyle \frac{\left(\frac{M}{P_{\mathrm{e}}}\right)S_{\mathrm{s}}}{K_{\mathrm{s}}}}\end{array}$$

where Q is total runoff in the river channel; Q_M is snowmelt runoff in the river channel; Q_f,M is surface flow or subsurface storm flow generated by snowmelt; Q_s,M is groundwater flow generated by snowmelt.

4.3 Model calibration and uncertainty estimation

For the Bogd Uliastai river basin, the model was warmed up using data from 2007, calibrated during 2008–2011, and validated during 2012–2015. For the Zavkhan Guulin river basin, 2000 was used as the warm-up year, with 2001–2010 selected for calibration and 2011–2020 for validation.

The MOSCEM-UA (Multi-objective Shuffled Complex Evolution Metropolis Algorithm) integrates multi-objective optimization and Bayesian uncertainty analysis, featuring global search capabilities that facilitate the generation of multiple Pareto optimal solutions and provide an assessment of uncertainty (Vrugt et al., 2003). The MOSCEM_UA was run for the optimization of parameters, with 40 000 iterations for each of the four model structures. The model parameters and their prior ranges for calibration are listed in Table 3.

The Kling-Gupta Efficiency (KGE) and its logarithmic form (KGL) were used as objective functions to evaluate the simulation of daily discharge (Gupta et al., 2009). These two metrics were chosen because each emphasizes a different portion of the hydrograph: KGE is more sensitive to high-flow dynamics, while KGL better captures low-flow conditions. In this study, to accommodate minimization-based optimization algorithms, the runoff objective functions L₁ (Eq. 29) and L₂ (Eq. 30) were formulated as one minus their respective efficiency metrics. The two objective functions were assigned equal weights during model calibration to ensure a balanced representation of both high- and low-flow regimes.

$$\begin{array}{}\text{(6)}& {\displaystyle }{L}_{\mathrm{1}}=\mathrm{1}-\mathrm{KGE}=\sqrt{{\left(\mathrm{1}-\mathit{\gamma }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-\mathit{\alpha }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-\mathit{\beta }\right)}^{\mathrm{2}}}\text{(7)}& {\displaystyle }{L}_{\mathrm{2}}=\mathrm{1}-\mathrm{KGL}=\sqrt{{\left(\mathrm{1}-{\mathit{\gamma }}_{\log }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\mathit{\alpha }}_{\log }\right)}^{\mathrm{2}}+{\left(\mathrm{1}-{\mathit{\beta }}_{\log }\right)}^{\mathrm{2}}}\end{array}$$

where, γ is the correlation coefficient between simulated and observed flows, and γ_log is the correlation coefficient between their logarithmic values; α is the ratio of the standard deviations of the simulated and observed flows, and α_log is the ratio of the standard deviations of their logarithmic transformations; β is the ratio of the mean values of the simulated and observed flows, and β_log is the ratio of the mean values of their logarithmic transformations.

5 Results

5.1 Model calibration and validation

Fig.4 shows the performance of the four models during the calibration period. The Pareto-optimal front shifts progressively toward the origin, indicating that model structural modifications enhance the model's ability to capture basin runoff dynamics. The FLEX^L-S model (KGE: 0.65 and 0.65, with the former representing the Bogd Uliastai river basin and the latter representing the Zavkhan Guulin river basin, hereinafter referred to as the same; KGL: 0.68 and 0.66) outperforms the baseline FLEX^L model (KGE: 0.53 and 0.52; KGL: 0.62 and 0.48) (Table 4). This improvement highlights the importance of explicitly representing snow processes in cryospheric regions. Without accounting for snow accumulation and ablation, the model tends to overestimate small peak flow events in winter, as shown in Fig.5, underscoring the critical role of snow dynamics in shaping hydrological responses.

The FLEX^D model (KGE: 0.77 and 0.68; KGL: 0.74 and 0.74) outperforms the FLEX^L-S, with the distributed precipitation and temperature inputs significantly improving the simulation of peak flow. FLEX^D does not require a more complex model structure or additional parameters compared to FLEX^L-S. For hydrograph simulation, FLEX^T (KGE: 0.78 and 0.68; KGL: 0.75 and 0.75) performs comparably to FLEX^D.

Figure 4 Performance of the FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T models during the calibration period.

Table 4 Comparison of simulation performance among different hydrological models in the two study catchments.

Some interesting rain/snowmelt-runoff events also suggest that semi-distributed models (FLEX^D and FLEX^T) better capture basin hydrological processes. Two such events in the Bogd Uliastai river basin in April 2009 and September 2011 provide compelling evidence (Fig. 5). On 14 April 2009, despite minimal precipitation, temperatures exceeded the melting threshold, producing only a relatively insignificant peak flow. In 12 September 2011, despite a higher daily precipitation of 12.7 mm, no runoff peak was observed within the basin. Lumped models failed to reproduce these dynamics accurately, instead simulating much larger peak flows.

Figure 5 The daily observed and simulated hydrographs of the FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T models in the calibration period. The dashed boxes represent the rainfall/snowfall-runoff events on 14 April 2009 and 12 September 2011 in the Bogd Uliastai river basin.

The performance and results of the four models during the validation period are shown in Figs. 6 and 7. The results confirm the stepwise improvement in model performance, as evidenced by the points corresponding to different model structures progressively shifting toward the origin. With the gradual optimization of model structure, the model's fitness has significantly improved. Unlike during calibration, the points in the validation period do not maintain the arc shape (Fig. 4).

Figure 6 Performance of the FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T models in validation mode.

Figure 7 The daily observed and simulated hydrographs of the FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T models in the validation period.

5.2 Model test by snowpack dynamics

Snow water equivalent is a crucial indicator of snowmelt dynamics and plays an essential role in hydrological modelling, serving as an additional metric for evaluating model performance and realism (Fig. 8). In the Bogd Uliastai river basin, the FLEX^D and FLEX^T models achieved KGE values of 0.61 and 0.63, respectively, for SWE simulation, indicating their ability to capture seasonal patterns and interannual variability, particularly peak values during winter and spring. The FLEX^T model, which incorporates vegetation effects, further improved SWE simulation accuracy and enhanced responsiveness to hydrological processes. In contrast, the FLEX^L-S model yielded a KGE of only 0.37, reflecting its limitations in capturing snowpack dynamics within the basin.

Figure 8 The observed and simulated daily snow water equivalent of the FLEX^L-S, FLEX^D, and FLEX^T models.

In the Zavkhan Guulin river basin, the FLEX^L-S model demonstrated relatively stable performance, achieving a KGE of 0.50. Although lumped models struggle to capture spatial heterogeneity, they effectively reflect seasonal precipitation and snowmelt trends. FLEX^D and FLEX^T achieved KGE values of 0.55 and 0.56, respectively, showed slight improvements. Overall, the semi-distributed model exhibited consistently higher SWE simulation skill across both basins.

5.3 Model parameters composition

Figure 9 shows that model parameters exhibit distinct sensitivity across different structural configurations within the stepwise modelling framework. In models that do not account for vegetation effects, single parameter values are used to approximate basin-average hydrological behavior, whereas FLEX^T produces spatially differentiated parameter values across landscape units.

Interception capacity (S_i,max) and root zone storage capacity (S_u,max) vary across landscapes. In the Bogd Uliastai river basin, forest (present only in this basin) shows the highest S_i,max (1.19 mm), followed by riparian area (0.55 mm) and grassland (0.05 mm), while bare soil/rock has no interception storage. S_u,max decreases from forest (145 mm) to grassland (54 mm), riparian area (29 mm), and bare soil/rock (14 mm). By comparison, in the Zavkhan Guulin river basin, S_i,max is lower (grassland: 0.02 mm; riparian area: 0.41 mm), whereas S_u,max is higher (grassland: 283 mm; riparian area: 105 mm; bare soil/rock: 33 mm), reflecting the more arid conditions and associated hydrological characteristics.

The recession coefficient of the fast response reservoir (K_f) reflects runoff dynamics and storage across landscapes. In the Bogd Uliastai river basin, forest has the largest K_f (9.4 d), followed by grassland (7.7 d), while riparian areas (3.0 d) and bare soil/rock (2.6 d) have the smallest K_f, indicating faster runoff. In the Zavkhan Guulin river basin, K_f shows a similar pattern, with grassland (9.7 d), riparian areas (2.9 d), and bare soil/rock (2.5 d). The low K_f in bare soil/rock reflects Hortonian overland flow, whereas runoff in riparian areas is mainly governed by saturation overland flow.

Figure 9 The changes of averaged behavioral parameters of the FLEX^L, FLEX^L-S, FLEX^D, and FLEX^T models.

Other parameters are also refined alongside improvements in the model structure. The parameter D, which partitions generated runoff between fast and slow response reservoirs, tends to be close to 1 – indicating that most runoff is routed to the fast reservoir. Parameters related to energy processes, such as the T_tm and F_dd, exhibit a clear compensatory relationship: a higher T_tm is typically associated with a lower F_dd, and vice versa.

5.4 The ratio of snowmelt runoff to streamflow

Figure10 shows the annual and seasonal ratios of snowmelt runoff to streamflow in the Bogd Uliastai and Zavkhan Guulin river basins, as simulated by the FLEX^T model. On an annual scale, snowmelt runoff accounts for 23.6 % ± 0.7 % and 15.9 % ± 1.3 % of streamflow in the Bogd Uliastai and Zavkhan Guulin river basins, respectively. Seasonally, snowmelt plays a dominant role in sustaining spring flows, with snowmelt runoff accounting for 70.1 % ± 1.7 % and 76.9 % ±  3.2 % of streamflow in the two basins, while its contribution is considerably lower in other seasons.

Figure 10 The ratio of snowmelt runoff to streamflow based on the FLEX^T model.

Figure 11 shows the snowfall-to-precipitation ratio ($P_{\mathrm{s}}/P)$ and the ratio of snowmelt runoff to streamflow ($Q_M/Q$) across different elevations. The $P_{\mathrm{s}}/P$ increases with elevation ranging from 14.7 % to 70.3 % in the Bogd Uliastai river basin and from 12.7 % to 63.2 % in the Zavkhan Guulin river basin. Consistent with this pattern, the Q_M/Q also increases with elevation, from 11.8 % to 57.8 % and from 10.3 % to 58.4 % in the two basins, respectively. Overall, higher $P_{\mathrm{s}}/P$ are associated with higher $Q_{\mathrm{M}}/Q$ across elevations.

Figure 11 The ratio of snowmelt runoff to streamflow and the snowfall-to-precipitation ratio at different elevations based on FLEX^T model.

5.5 Runoff generation at different elevations and across landscapes

Figure12 shows significant differences in runoff contributions across 5 equal area elevation bands. In the Bogd Uliastai and Zavkhan Guulin river basins, high elevation areas (above 2900 m in the Bogd Uliastai basin and 2825 m in the Zavkhan Guulin basin) play a dominant role in runoff generation, contributing approximately 30 % of total basin runoff. This dominance is primarily attributed to orographic effect, which increase precipitation and enhance the proportion of snowmelt-derived runoff (Fig. 11). In contrast, low elevation areas rely primarily on rainfall-induced runoff. Due to limited precipitation and higher evaporative losses, their contributions to total runoff are comparatively smaller, with the lowest elevation band accounting for only about 15 % of basin runoff.

Figure 12 Runoff contribution from 5 equal area elevation bands (each representing 20 % of the total catchment area) based on FLEX^T model.

Different landscape units also exhibit differences in their contributions to runoff (Fig. 13). Grasslands occupy the largest portion of both basins (82.8 % and 85.4 % in the Bogd Uliastai and Zavkhan Guulin river basins, respectively) and contribute correspondingly high shares of total runoff (80.0 % and 82.0 %). Forested areas, found only in the Bogd Uliastai basin, account for 4.7 % of the basin area but contribute only 2.9 % of the total runoff. Riparian areas occupy a relatively small fraction of the two basins (10.3 % and 13.0 %) but contribute a comparable or higher proportion of runoff (12.8 % and 14.7 %). Bare soil/rock areas generate a relatively high runoff contribution per unit area; however, their limited spatial extent (2.2 % and 1.6 %) constrains their total contribution to streamflow (4.3 % and 3.3 %).

Figure 13 Runoff contribution from different landscapes based on FLEX^T model.

6 Discussion

6.1 The stepwise modelling framework improves the simulation of hydrological processes in cold regions

Accurately representing dominant hydrological processes and spatial heterogeneity is essential for simulating runoff dynamics in cold regions. The enhanced performance of FLEX^L-S relative to the FLEX^L model underscores the necessity of explicitly representing snow accumulation and melt, which are integral to runoff generation in cryospheric environments. Without such representation, winter runoff responses tend to be systematically mischaracterized.

The better performance of the distributed FLEX^D model highlights the importance of elevation-dependent variability in temperature and precipitation (Fenicia et al., 2008b). By maintaining distinct storage states across hydrological response units, FLEX^D captures spatial contrasts in precipitation phase and melt timing that lumped models fail to resolve, particularly during rain–snow transition events (Bormann et al., 2009).

FLEX^T systematically integrates topography, vegetation, and hydrological processes through the explicit representation of hydrological response units, thereby enhancing the physical consistency of the model structure and the interpretability of its parameters (Savenije, 2010). Under the basin conditions considered in this study, FLEX^T exhibits a high degree of consistency with FLEX^D in simulating both runoff and snowpack dynamics. Although process-based discretization approaches have been shown to improve basin-scale hydrological simulations in heterogeneous cold region basins – such as the Heihe (Gao et al., 2014) and Yigong Zangbu river basins (Gao et al., 2020) – when topographic and land-cover variability are explicitly represented, the magnitude of these benefits appears to depend strongly on the degree of landscape heterogeneity.

In the studied basins, vegetation heterogeneity is limited, with grasslands covering more than 80 % of the area (Fig. 13), which substantially reduces contrasts in vegetation properties across the landscape. Although the lumped model applies a spatially uniform NDVI and therefore cannot explicitly resolve intra-basin vegetation variability, NDVI exhibits a strong correspondence with elevation, particularly in grassland-dominated regions. NDVI values across elevation bands closely resemble those of the corresponding grassland zones, indicating that vegetation distribution is largely structured by topographic gradients. While forested and bare land areas display distinct NDVI signatures, they occupy only a minor fraction of the basin and exert a limited influence on overall runoff generation (Fig. 14).

Figure 14 Multi-year average NDVI variation across landscapes and its relationship with elevation in two study basins.

Under these conditions, key vegetation-related hydrological attributes – such as root zone storage and interception capacity – can be effectively represented through hydroclimatic conditions and topographic gradients (Roebroek et al., 2020; Antonelli et al., 2018), allowing elevation to function, to some extent, as a proxy for vegetation structure. Consequently, the comparable performance of FLEX^T and FLEX^D in simulating runoff and snowpack processes primarily reflects the constraining influence of landscape structure in the studied basins, rather than differences in model design. Within this stepwise modelling framework, FLEX^T employs a more physically grounded parameterization scheme and therefore offers a flexible and extensible platform for diagnosing hydrological processes in cold mountainous regions. This suggests a particular advantage for exploring hydrological responses to vegetation and climate change in basins with higher ecological or topographic complexity.

6.2 Process-based tracking of the ratio of snowmelt runoff to streamflow

Although direct observational data (e.g., stable water isotopes) for quantifying the ratio of snowmelt runoff to streamflow are unavailable in this study, previous research provides indirect support. Wu et al. (Wu et al., 2021) applied a similar snowmelt tracking method in the Altai Mountains and reported that snowmelt accounted for 29.3 % of annual streamflow. This result exceeds those observed in our study, largely due to regional differences in snowfall. In the Kayiertesi river basin of the Altai Mountains, annual average precipitation for one hydrological year (September to August) was 409.8 mm from 2011 to 2015 (observed at the Kuwei snow station), with snowfall from November to March comprising about 31 % of that annual precipitation (Zhang et al., 2017). In contrast, annual precipitation in the Bogd Uliastai and Zavkhan Guulin basins does not exceed 200 mm, and snowfall represents less than 15 % of the total observed precipitation.

This study also compared model-based snowmelt tracking with traditional indirect methods, which estimate the ratio of snowmelt runoff to streamflow by calculating the ratio of snowfall or snowmelt to runoff over a given period (Barnett et al., 2005; Immerzeel et al., 2010). While these indirect approaches are computationally efficient and require limited data, they implicitly assume that all meltwater is instantaneously converted to runoff, neglecting interactions with rainfall and losses due to infiltration, evaporation, and subsurface storage (Li et al., 2017).

Using the traditional indirect approach, we calculated snowmelt-to-runoff ratios of 41.3 % ± 1.6 % and 153 % ± 15.5 % in the two basins, respectively. These estimates are substantially higher than those obtained from model-based snowmelt tracking, with some values exceeding 100 %, indicating physical implausibility. This discrepancy highlights the methodological limitations of traditional approaches.

The overestimation inherent in traditional methods arises from their failure to represent snowmelt that infiltrates into the soil profile, percolates to deeper layers, or is lost through evaporation, thereby reducing the fraction of meltwater reaching the channel (Liu et al., 2023), particularly in the arid Zavkhan Guulin river basin. These findings underscore the importance of using physically based models to trace water source pathways, particularly in data-scarce and hydrologically complex regions.

6.3 Runoff generation mechanisms at different elevations and across landscapes

Elevation is a fundamental control on runoff generation and seasonal hydrological variability in mountainous basins, as it governs precipitation phase, snow accumulation, and melt timing (Jenicek and Ledvinka, 2020). In the study basins, runoff contributions increase progressively with elevation. This pattern contrasts with glacierized catchments, where the largest runoff contributions do not necessarily originate from the highest elevation bands but are often dominated by glacier melt processes (Gao et al., 2020). These elevation-dependent contrasts highlight the importance of explicitly accounting for vertical heterogeneity and basin-specific cryospheric conditions in cold-region hydrological modelling.

The lag effect in snowmelt runoff across different elevations is a defining feature of hydrological processes in mountainous basins, arising from the temporal variability of snowmelt among elevation areas. In the Bogd Uliastai and Zavkhan Guulin river basins, low elevation areas contribute rapidly to runoff during the initial stages, while high elevation areas gradually augment runoff as the season progresses, reflecting the delayed response of snowmelt-dominated areas (Hajika et al., 2024). This lag is primarily governed by elevation-dependent snow accumulation and melt dynamics: under cold conditions at higher elevations, precipitation is stored as snow, and meltwater release typically occurs several weeks to a few months later than at lower elevations, with the lag being most pronounced at the onset of the melt season (Gillan et al., 2010) (Fig. 12).

The lag in snowmelt runoff across different elevations has profound implications for water resources management. Snowmelt from high elevation areas prolongs the recession of basin runoff, maintaining downstream water supply and providing buffering capacity during droughts or low-flow periods. In contrast, low elevation areas respond rapidly to extreme precipitation, generating short-term flood peaks, while delayed high elevation snowmelt can extend high flow conditions, complicating flood management (Gu et al., 2023). Effective management in mountainous basins therefore requires integrating both immediate and delayed runoff contributions across elevations to optimize water allocation, improve flood early-warning accuracy, and safeguard downstream ecological and socio-economic water demands (Li et al., 2019).

Beyond elevation controls, vegetation influences runoff generation through its direct physiological properties (e.g., interception and root zone storage). Grasslands are the dominant land-cover type in both basins and are widely distributed across different elevation zones. Their runoff-generation characteristics are relatively spatially uniform in a basin-averaged sense and thus exert a primary control on the overall runoff response. Under the semi-arid and arid climatic setting, soils in grassland areas are generally persistently dry, with a large initial soil-moisture deficit. This condition thereby enhances the infiltration demand for rainfall and snowmelt, delaying runoff generation (Xue et al., 2025). Forested areas exhibit strong regulatory functions – intercepting precipitation, enhancing infiltration, and reducing quick flow generation (Stocker et al., 2023; Liu et al., 2018). Riparian areas function as hydrologically connected source areas, rapidly transmitting rainfall and snowmelt to the channel network (Leibowitz et al., 2023). The runoff generation capacity of bare soil/rock is high due to the lack of vegetation and low soil permeability. After rainfall or snowmelt, water infiltrates poorly and rapidly forms overland flow (Zeng et al., 2024). Although their basin-scale contribution is constrained by areal extent, these surfaces can locally intensify runoff responses and sediment transport.

These findings support the view that vegetation functions as a spatially variable regulator of runoff generation. This regulatory effect is particularly sensitive to future changes in vegetation cover and distribution under climate and land use change scenarios. For instance, overgrazing may reduce aboveground biomass and root zone storage capacity, thereby increasing runoff and erosion risks (Donovan and Monaghan, 2021), while shifts in vegetation type (e.g., shrub encroachment or forest decline) could alter hydrological partitioning along elevation gradients (Zhou et al., 2023; Hsu et al., 2025).

6.4 Applicability and limitations of the modelling framework

Previous applications of the FLEX framework in alpine basins, such as the Heihe river basin, provide an important reference for interpreting the results of this study. In the Heihe river basin, FLEX^T was configured using terrain-driven functional units, whereby the catchment was discretized into four topographic sub-units representing distinct dominant hydrological processes. This modelling strategy primarily focused on assessing model realism and transferability across sub-basins by demonstrating the ability of FLEX^T to reproduce contrasting functional hydrological behaviors (Gao et al., 2014). Building on these earlier studies, the present work further demonstrates that the FLEX framework can also be adapted to colder and more semi-arid conditions using a stepwise strategy. Importantly, our emphasis extends beyond discharge performance to diagnosing runoff-generation mechanisms—specifically the snowmelt contribution to streamflow and the relative roles of elevation bands and vegetation units—thereby clarifying how topography and land cover jointly shape hydrological responses in cold alpine basins.

While the proposed modelling framework captures key hydrological processes, several limitations should be acknowledged. First, frozen soil processes were not explicitly represented in this study due to data constraints. Frozen soil is widespread across the region and can exert strong controls on hydrological processes by restricting moisture exchange between deeper soil layers and the atmosphere, reducing soil permeability, and modifying near-surface flow pathways (Li et al., 2025). These effects may influence the partitioning between surface runoff and subsurface flow, particularly in high-elevation permafrost regions with sparse vegetation (Hu et al., 2025). In this study, the elevation-dependent runoff generation patterns identified may partly reflect the integrated influence of frozen soil conditions in these regions. Future studies that explicitly incorporate frozen soil dynamics would help to further disentangle the interacting roles of topography, vegetation, and frozen ground in shaping hydrological responses in alpine basins.

Second, regarding the assumption of complete mixing in snowmelt runoff tracking, the applied conceptual model represents storages as lumped and homogeneous reservoirs, within which incoming rainfall and snowmelt are assumed to mix instantaneously with pre-event water. Despite its limited ability to explicitly represent rapid snowmelt runoff responses and flow-path heterogeneity at the event scale (Weiler et al., 2018), this assumption is still widely applied in catchment-scale estimates of snowmelt runoff contributions and is regarded as robust at long time scales (Liu et al., 2023; Wu et al., 2021; Li et al., 2017). Importantly, the snowmelt tracker is designed to estimate the long-term average of the fraction of total streamflow that originates from snowmelt, rather than event-scale responses.

Finally, the use of two study basins entails both advantages and limitations. Focusing on a limited number of basins allows for detailed process interpretation and reduces confounding effects associated with excessive climatic or physiographic diversity, but it also constrains the generalization of the findings (Hrachowitz et al., 2013). The identified runoff generation mechanisms and model performance characteristics should therefore be interpreted as representative of cold, semi-arid alpine basins with similar hydroclimatic and geomorphic settings, rather than as universally applicable. Expanding the analysis to a larger set of basins would be necessary to assess model robustness and transferability more rigorously.

7 Conclusions

Hydrological modelling in cold and high-latitude regions remains challenges due to complex cryospheric processes and limited observational data. To address these issues, this study developed and tested a stepwise modelling framework that progressively incorporates key hydrological processes and landscape characteristics, with the aim of improving both predictive performance and process interpretability.

The results demonstrate that model structure matters. The lumped models (FLEX^L and FLEX^L-S) are insufficient for capturing runoff dynamics in basins with complex topography and heterogeneous vegetation. Although FLEX^D improves simulation by using spatially distributed inputs, it still lacks full physical interpretability. In contrast, the landscape-based FLEX^T model integrates snowpack, topography, and vegetation characteristics, providing more physical meaningful parameterization and a more mechanistic representation of hydrological processes. Validation using SWE further confirms the FLEX^T's ability to capture seasonal dynamics, interannual variability, and key hydrological mechanisms in cryospheric environments. These findings underscore the potential advantages of FLEX^T, particularly in basins with greater ecological or topographic complexity.

Results from the FLEX^T model indicate that snowmelt runoff accounts for 23.6 % ± 0.7 % and 15.9 % ± 1.3 % of streamflow in the Bogd Uliastai and Zavkhan Guulin river basins, respectively. Snowmelt runoff dominates streamflow during spring and increases with elevation, underscoring the critical role of topography in shaping runoff generation. In high elevation areas, the lagged snowmelt response leads to a sustained and gradual release of runoff, whereas low elevation areas dominated by rainfall generate rapid runoff. Controlled by distinct dominant hydrological mechanisms, different landscape units contribute unequally to streamflow.

These findings offer valuable insights into hydrological response mechanisms in cold alpine basins on the Mongolian Plateau. The stepwise modelling framework enhances model realism in cryospheric environments and provides practical support for water resource management, ecological conservation, and climate adaptation in high-latitude and data-scarce regions.