Hybrid hydrological modeling for large alpine basins: a semi-distributed approach

Alpine basins are important water sources for human life, and reliable hydrological modeling can enhance the water resource management in alpine basins. Recently, hybrid hydrological models, coupling process-based models and deep learning (DL), have exhibited considerable promise in hydrological simulations. However, a notable limitation of existing hybrid models lies in their failure to incorporate spatial information within the basin and describe alpine hydrological processes, which restricts their applicability in hydrological modeling in large alpine basins. To address this issue, we develop a set of hybrid semi-distributed hydrological models by employing a process-based model as the backbone and utilizing embedded neural networks (ENNs) to parameterize and replace different internal modules. The proposed models are tested on three large alpine basins on the Tibetan Plateau. A climate perturbation method is further used to test the applicability of the hybrid models to analyze the hydrological sensitivities to climate change in large alpine basins. Results indicate that proposed hybrid hydrological models can perform well in predicting runoff processes and simulating runoff component contributions in large alpine basins. The optimal hybrid model with Nash–Sutcliffe efficiencies (NSEs) higher than 0.87 shows comparable performance to state-of-the-art DL models. The hybrid model also exhibits remarkable capability in simulating hydrological processes at ungauged sites within the basin, markedly surpassing traditional distributed models. In addition, the results also show reasonable patterns in the analysis of the hydrological sensitivities to climate change. Overall, this study provides a high-performance tool enriched with explicit hydrological knowledge for hydrological prediction and improves our understanding about the hydrological sensitivities to climate change in large alpine basins.

1 Introduction

Alpine basins are important water sources, playing a crucial role in various aspects of human life and the environment, such as domestic water supply, irrigation, hydropower generation and climate regulation (Cui et al., 2023; Huss et al., 2017; Viviroli et al., 2011). Developing reliable hydrological models is crucial for managing floods and improving water use efficiency under climate change (Blöschl et al., 2019).

Process-based hydrological models, such as EXP-Hydro (Patil and Stieglitz, 2014), CRHM (DeBeer and Pomeroy, 2017) and THREW (Nan et al., 2021), are widely used approaches for hydrological simulation in large alpine basins. These models depend on physical laws and empirical knowledge to describe physical processes and are grounded in well-defined physical mechanisms. They can be used to advance scientific understanding about the hydrological systems and provide the insight into the response of hydrological processes to climate changes (Cui et al., 2023; Li et al., 2021). However, the performance of these models is constrained by several factors, including an incomplete understanding of alpine hydrological processes, errors in the model structure and uncertainties in parameterization (Kuppel et al., 2018; Beven, 2006). These deficiencies also give rise to equifinality, making it challenging to accurately represent hydrological processes. This diminishes the credibility of process-based models in the context of climate change assessment.

Deep learning (DL) hydrological models are distinguished by their remarkable data mining capabilities, operating independently of hydrological knowledge. They have showcased exceptional model performance across diverse hydrological domains, including streamflow/discharge forecasting (Kratzert et al., 2018; Lees et al., 2021; Liu et al., 2021), snow water equivalent modeling (Duan and Ullrich, 2021) and groundwater level mapping (Solgi et al., 2021; Nourani et al., 2022). Most of these studies disregard the effect of spatial information from meteorological data on hydrological modeling. Li et al. (2023a) introduced an innovative spatiotemporal DL hydrological model, demonstrating that integrating spatial information can significantly improve the performance of DL models in hydrological modeling. Nonetheless, despite their remarkable capabilities, DL hydrological models still face scrutiny within the hydrological modeling community, primarily due to their “black-box” nature. Furthermore, DL models rely on the assumption that the dataset's distribution during the prediction period remains consistent with that of the training period. This assumption cannot be met when using DL models to assess the effects of climate change on hydrological modeling (Nearing et al., 2021; Zhong et al., 2023).

Hybrid hydrological models that combine process-based and DL approaches are anticipated to harness their respective strengths to achieve both impressive performance and a well-defined understanding of hydrological processes (Tsai et al., 2021; Shen et al., 2023). Previous studies have introduced various hybrid model configurations and demonstrated satisfactory outcomes (Feigl et al., 2022; Frame et al., 2021; Kashinath et al., 2021; Quilty et al., 2022; Bhasme et al., 2022; Kumanlioglu and Fistikoglu, 2019; Xie et al., 2021; Lu et al., 2021), while the underlying concept in many of these hybrid models remains centered on either pure DL models or process-based models. For instance, Frame et al. (2021) utilized long short-term memory (LSTM) models as post-processors for the United States National Water Model, highlighting that integrating DL models can improve performance by rectifying errors in the outcomes of process-based models. Xie et al. (2021) introduced a physically guided LSTM model by incorporating synthetic samples during model training to capture underlying physical mechanisms. Recently, some studies attempted to implement differentiable models to facilitate a bidirectional integration between process-based models and DL models (Shen et al., 2023; Baydin et al., 2018; Höge et al., 2022). Feng et al. (2022) introduced hybrid hydrological models that integrated a lumped hydrological model HBV as the foundation and incorporated embedded neural networks (ENNs) to parameterize, enhance or replace internal components without prior training. The proposed models demonstrated comparable performance to DL models and can output untrained physical variables. Our earlier work further developed hybrid models by employing ENNs to replace the internal modules of the lumped model EXP-Hydro and systematically test the impact of replacing different internal modules with ENNs (Li et al., 2023b). The findings suggest that substituting any internal component with ENNs can enhance model performance, but increasing the number of internal component replacements does not guarantee improved outcomes. Achieving optimal performance requires a delicate equilibrium between the quantity of ENNs and the process constraints inherent in the process-based model. However, Feng et al. (2022) and Li et al. (2023b) have predominantly employed lumped hydrological models as the foundational framework in hybrid models. They have not adequately accounted for the spatial information of meteorological inputs and underlying surfaces within the basin, which limits their applicability in large basins. Additionally, the effectiveness of hybrid models in the Tibetan Plateau's large alpine basins, particularly in assessing hydrological sensitivities to climate change, is yet to be clearly established. Therefore, there is a need to evolve hybrid models from lumped to distributed to adequately capture the spatial information within the basin. Moreover, it is also essential to incorporate alpine hydrological processes in hybrid models for adapting them to alpine basins and evaluate the adaptability of these hybrid models in analyzing the hydrological sensitivities to climate change in large alpine basins.

Building upon our earlier work about hybrid lumped models (Li et al., 2023b), this study aims to propose hybrid semi-distributed models that employ a hydrological model as the backbone and employ ENNs to parameterize and replace different internal modules within the sub-basin scale. The proposed models are then comprehensively assessed across three large mountainous basins on the Tibetan Plateau. A climate perturbation method is further used to analyze the hydrological sensitivities to climate change in large alpine basins. The remainder of this paper is organized as follows: Sect. 2 outlines the proposed hybrid models, study area and data; Sect. 3 shows the evaluation results of the proposed models; Sect. 4 provides details about the hydrological sensitivities to climate change; and we conclude in Sect. 5.

2 Methods and materials

2.1 Model development

This study develops hybrid semi-distributed hydrological models by integrating the process-based model and embedded neural networks (ENNs; Fig. 1). Specifically, the proposed models use a semi-distributed EXP-Hydro model as the backbone, with ENNs parameterizing and replacing different internal modules. The differential programming framework is utilized to achieve a bidirectional integration between the process-based model and ENNs, enabling simultaneous parameter training of both entities.

Figure 1 The schematic diagram of hybrid semi-distributed models. (a) The basin is first divided into many sub-basins and (b) all meteorological and hydrological processes included in the EXP-Hydro model are calculated in each sub-basin. The precipitation partition, snowmelt and runoff modules can be optionally replaced by embedded neural networks. For detailed formulations of these processes, refer to the main text.

2.1.1 The semi-distributed EXP-Hydro model

In this study, the hybrid semi-distributed models are built upon the foundation of the semi-distributed EXP-Hydro model (Patil et al., 2014). The originally lumped EXP-Hydro model, proposed by Patil and Stieglitz (2014), treats each basin as a singular areal unit, disregarding the spatial information within the basin. The EXP-Hydro model encompasses a snow accumulation bucket and a basin bucket represented by snow storage (S₀) and basin water storage (S₁), respectively. Within the model, four processes are represented: precipitation partition (rainfall P_r or snowfall P_s), evapotranspiration (ET), snowmelt (M) and runoff (Q). For detailed equations, refer to Appendix A1 and Patil and Stieglitz (2014). The semi-distributed EXP-Hydro model was subsequently extended to incorporate the spatial heterogeneity within the basin (Patil et al., 2014). Initially, the study basin is divided into multiple sub-basins using a digital elevation model (DEM). The EXP-Hydro model is run independently within each sub-basin, and the overall basin runoff is derived by summing the runoff outputs from all sub-basins (Eq. A12). Patil and Stieglitz (2015) and Patil et al. (2014) showcased the efficacy of the semi-distributed EXP-Hydro model in hydrological modeling across 295 basins spanning the continental United States. Their studies indicated that this model outperforms the original EXP-Hydro model.

2.1.2 The hybrid semi-distributed models

Using the semi-distributed EXP-Hydro model as the backbone, the hybrid models integrate ENNs to parameterize and replace various internal modules within the differential programming framework (Baydin et al., 2018). This configuration enables the model to comply with basic physical principles while enhancing its representational capability of the corresponding meteorological and hydrological modules, thus increasing the accuracy of hydrological simulations. ENNs utilize both static attributes (Table A1) and dynamic meteorological time series from each sub-basin as inputs. These inputs are employed to characterize the disparities in physical mechanisms among sub-basins and to drive the precipitation–runoff processes. The hybrid models are realized via four steps:

1. Data pre-processing. A DEM is employed to partition the study basin into multiple sub-basins, guided by a drainage area threshold (Grieve et al., 2016; Noël et al., 2014). The static attributes (Table A1) and daily meteorological time series for each sub-basin are derived by calculating the areal averages from the original dataset.

2. Semi-distributed model development within the differential programming framework. All equations within the hydrological model are formulated to be differentiable to ensure operation within the differential programming framework (Shen et al., 2023; Li et al., 2023b; Levine et al., 2016). This framework facilitates the computation of derivatives from model outputs to inputs and intermediate variables, thus enabling an “end-to-end” training approach. The hybrid model achieves simultaneous training of both the semi-distributed hydrological models and ENNs. Only runoff data are employed as the training target, eliminating the need for observed data for ENN outputs. Furthermore, a physical recurrent neural network (P-RNN) is established to simulate hydrological dynamic processes and retain the memory of past basin storage sequences (Li et al., 2023b; Jiang et al., 2020).

3. ENN parameterization and replacement. The calibration parameters of all sub-basins within the basin are assumed to be the same in the semi-distributed EXP-Hydro model, while many of them related to sub-basin attributes should be different (Feng et al., 2022). To capture the spatial diversity of these calibration parameters at the sub-basin scale, we build an ENN to derive calibration parameters only using static attributes as inputs. Additionally, ENNs are employed to potentially substitute distinct internal modules of the EXP-Hydro model, utilizing static attributes and corresponding dynamic time series as inputs. Specifically, three ENNs are designed for simulating runoff, precipitation partition and snowmelt processes in this study.

4. Model training. Through the aforementioned steps, all parameters of the hybrid models, encompassing the EXP-Hydro model and ENNs, can be jointly trained using observed runoff data as the training target. The Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) is utilized as the loss function during training.

Our previous study has shown that the utilization of ENNs to substitute internal components of lumped hydrological models can elevate model performance in hydrological modeling (Li et al., 2023b). ENNs possess the flexibility to optionally replace any single or multiple internal modules of the hydrological model. Similar to Li et al. (2023b), the ENN dedicated to precipitation partition employs precipitation and air temperature as inputs to compute the snowfall ratio. Rainfall is then determined by subtracting snowfall from the precipitation. The snowmelt ratio is determined through an ENN that takes air temperature as the input. The ENN related to the runoff process is developed using basin water storage, the combined value of rainfall and snowmelt, and air temperature as inputs. The inclusion of air temperature serves to depict the influence of soil freeze–thaw dynamics on the runoff process in alpine basins within the Tibetan Plateau (Zhong et al., 2023). Apart from the dynamic driving time series, all ENNs utilized for replacing internal components also incorporate static attributes as inputs, aiming to differentiate disparities among various sub-basins. The detailed ENNs inputs refer to Appendix B. ENN_θ is used to represent the ENNs that parameterize the process-based model. ENN_Q, ENN_S and ENN_M are utilized hereinafter to denote the ENN that replace runoff, precipitation partition and snowmelt processes, respectively.

In this study, we develop and evaluate five hybrid models denoted as DM_θ, DM_θ−Q, DM${}_{\mathit{\theta }-Q-T}$, DM_θ−QSM and DM${}_{\mathit{\theta }-QSM-T}$ (Table 1). The DM_θ model solely employs the ENN_θ for parameterizing calibration parameters across sub-basins. The DM_θ−Q and DM${}_{\mathit{\theta }-Q-T}$ models go a step further by incorporating ENN_Q to replace the runoff process. Expanding upon this, the precipitation partition and snowmelt processes are substituted by corresponding ENNs in DM_θ−QSM and DM${}_{\mathit{\theta }-QSM-T}$ models. Notably, the inputs for the ENN_Q include air temperature in DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models, while DM_θ−Q and DM_θ−QSM models do not consider it.

Table 1 Design details of different hybrid models. “✓” represents that the model employs the corresponding ENNs, while “×” means that it does not.

2.1.3 Comparison models

We also compare our proposed models with the state-of-the-art distributed hydrological model THREW (Tsinghua Representative Elementary Watershed) and deep learning models LSTM and CNN-LSTM. The THREW model, originally proposed by Tian et al. (2006), operates by delineating the basin into representative elementary watersheds (REWs) through DEM calculation. Furthermore, each REW is subdivided into sub-zones, which serve as the fundamental units for hydrological modeling. The THREW model has demonstrated successful applications across diverse basins, including representative ones within the Tibetan Plateau, Alps and Tianshan (Cui et al., 2023; He et al., 2014). To establish a fair comparison of model performance between the THREW model and the proposed hybrid models, the THREW model in this study is subjected to the same spatial discretization utilized by the hybrid models. LSTM models (Hochreiter and Schmidhuber, 1997) have recently shown excellent capabilities in hydrological simulation all over the world (Kratzert et al., 2019; Lees et al., 2021; Li et al., 2023a). To benchmark against our proposed hybrid models, we have sourced the LSTM and CNN-LSTM model results from Li et al. (2023a). These models are renowned for their superior accuracy in existing deep learning research studies within the study basins. Furthermore, we also include the hybrid lumped hydrological models EXP_Q and EXP_QSM, proposed by Li et al. (2023b), for comparative evaluation. Their backbone model is the lumped hydrological model EXP-Hydro. This allows us to assess the effect of spatial information on hydrological modeling within hybrid frameworks. Notably, the EXP_Q and EXP_QSM employ the same dynamic time series inputs of ENNs for module replacement as the DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models, respectively. In addition, DM and EXP are utilized hereinafter to denote semi-distributed and lumped EXP-Hydro models if not specified otherwise.

2.2 Study area and data

2.2.1 Study area

The Tibetan Plateau (TP; Fig. 2), acclaimed as the “Third Pole” and the “water tower of Asia”, stands as the world's highest plateau. The TP provides a significant source of abundant water resources, crucial for the sustenance of downstream communities. To evaluate the performance of proposed hybrid models in large alpine basins, this study focuses on the source regions of three major river basins: the Yellow River, the Yangtze River and the Lancang River. These basins are recognized as extensive mountainous regions within the TP (Fig. 2). Each of these study basins spans an area exceeding 90 000 km², characterized by diverse topography, with elevation fluctuations exceeding 3000 m. The significant topographic variations within the basin lead to notable spatial heterogeneity in meteorological elements such as precipitation and temperature (Fig. A1). To accurately capture this heterogeneity in hydrological modeling, it is necessary to divide the basin into different computational units. In addition, previous studies have shown that the glacier process has a minimal impact on runoff modeling in the three study basins, and it is neglected in this study (Cui et al., 2023). Hereinafter, Yellow, Yangtze and Lancang are used to denote the corresponding source regions in this study.

Figure 2 The terrain of the Tibetan Plateau and the location of the four study basins.

2.2.2 Data used

This study utilized the reanalysis and remote sensing datasets for input variables of hybrid models and the THREW model as follows:

1. for precipitation, the China Meteorological Forcing Dataset (CMFD) with 0.1° spatial and 3 h temporal resolution (Yang et al., 2010);

2. the air temperature at 2 m a.g.l. (T2) from the fifth generation of ECMWF atmospheric reanalysis of the global climate (ERA5) reanalysis dataset with 0.1° spatial and 1 h temporal resolution (Hersbach et al., 2020);

3. the potential evaporation from the ERA5 reanalysis dataset with 0.1° spatial and 1 h temporal resolution (Hersbach et al., 2020);

4. for the DEM, the Shuttle Radar Topography Mission (SRTM) with 90 m spatial resolution, provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences (http://www.gscloud.cn, last access: 12 May 2022);

5. for the leaf area index (LAI), the MOD15A2H dataset from MODIS product with 500 m spatial and 8 d temporal resolution (Myneni et al., 2015);

6. for the normalized difference vegetation index (NDVI), the MOD13A3 dataset from MODIS product with 1 km spatial and 1-month temporal resolution (Didan, 2015).

The daily observed runoff data at hydrological stations (Fig. 2) are used for the model calibration/training and evaluation. The dataset is provided by local water agencies.

2.3 Experimental design

2.3.1 Model evaluation schemes

We conduct two sets of experiments to comprehensively evaluate the performance of proposed hybrid semi-distributed hydrological models in this study.

1. Model performance in trained sites. All proposed hybrid semi-distributed models are developed, trained and evaluated in three study basins. The comparison models are then utilized for a range of purposes: comparing the performance of the proposed models against state-of-the-art DL and distributed hydrological models, examining the effects of ENNs parameterization and replacement on hydrological modeling, and appraising the impact of spatial information on model performance. Due to the limitation of the observed runoff data, TNH in Yellow, ZMD in Yangtze and JZ in Lancang are utilized as the evaluation stations in this experiment. For Yellow and Yangtze, the training and evaluation periods are, respectively, designated as 1982–2004 and 2007–2014. In the case of the Lancang, these periods span 1988–2003 and 2005–2010.

2. Model performance in untrained sites within the basin. By capturing the spatial heterogeneity within the basin, hybrid semi-distributed models provide the opportunity to predict hydrological processes at any untrained sites within the basin. To assess the proficiency of hybrid semi-distributed models in ungauged sites within the basin, the MT, MQ and JG stations, situated upstream of the TNH station in the Yellow (Fig. 2), are simulated using Yellow (TNH) hydrological models in this section. The evaluation phase encompasses the years 2009 to 2014 for all hydrological stations.

2.3.2 The climate perturbation method

This study uses the climate perturbation method to test the applicability of the hybrid models to analyze the hydrological sensitivities to climate change in three large alpine basins. Using precipitation and temperature data from the reanalysis dataset (Sect. 2.2.2) as the reference, the additional perturbation sequences are added to represent the potential climate changes. Perturbed precipitation sequences are extracted by multiplying the reference precipitation data from 80 % to 120 % with an increment of 10 % (Su et al., 2023). Perturbed temperature sequences are generated by adding from 0.5 to 2 °C with an increment of 0.5 °C to the reference temperature input (Cui et al., 2023). The impact of increased temperature on the potential evapotranspiration is calculated by the regression between observed temperature and potential evapotranspiration in each sub-basin (Cui et al., 2023; van Pelt et al., 2009; Xu et al., 2019). In total, one reference, four perturbed temperature and four perturbed precipitation sequences are conducted to assess the influence of precipitation and temperature change on hydrological processes. The changes of other underlying surfaces are not considered in this study.

2.3.3 Evaluation metrics

Three common hydrological metrics – including NSE, modified NSE (mNSE; Legates and McCabe , 1999) and the absolute value of peak flow bias (PFAB; Yilmaz et al., 2008) – are employed to evaluate the model performance. They can be defined as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^T{\left(Q_{\mathrm{obs},i}-Q_{\mathrm{sim},i}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^T{\left(Q_{\mathrm{obs},i}-\stackrel{\mathrm{‾}}{Q_{\mathrm{obs}}}\right)}^{\mathrm{2}}}}\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{mNSE}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^T\left|Q_{\mathrm{obs},i}-Q_{\mathrm{sim},i}\right|}{\sum _{i=\mathrm{1}}^T\left|Q_{\mathrm{obs},i}-\stackrel{\mathrm{‾}}{Q_{\mathrm{obs}}}\right|}}\end{array}$$

$$\begin{array}{}\text{(3)}& \mathrm{PFAB}=\mathrm{100}\times \left|{\displaystyle \frac{\sum _{l=\mathrm{1}}^L\left(Q_{\mathrm{sim}:l}-Q_{\mathrm{obs}:l}\right)}{\sum _{l=\mathrm{1}}^L{Q}_{\mathrm{obs}:l}}}\right|,\end{array}$$

where Q_obs,i and Q_sim,i are the observed and simulated values, T is the length of the evaluation period, and $\stackrel{\mathrm{‾}}{Q_{\mathrm{obs}}}$ is the averaged observed values. Q_sim:l and Q_obs:l are the observed and simulated runoff sorted in descending order, respectively. L is the number of flow values that are in the top 2 % of all flows. Both NSE and mNSE measure the overall goodness of fit of simulated and observed data, while mNSE gives less weight to high values than NSE and thus focuses on the baseflow. A NSE and a mNSE of 1 indicate a perfect fit, and a NSE of 0.55 is the threshold for good performance (Newman et al., 2015; Knoben et al., 2019). PFAB emphasizes the performance for peak values, and a value closer to zero indicates a smaller peak bias.

3 Model evaluation

To adopt the hybrid semi-distributed models and THREW models in three basins, the Yellow, Yangtze and Lancang basins are delineated into 83, 99 and 63 sub-basins (Fig. 2) based on the actual river network and the divided sub-basin numbers in other relevant studies (Cui et al., 2023). The performance of all proposed models in gauged and ungauged sites is evaluated as follow.

3.1 Hybrid model evaluation in trained sites

3.1.1 The effect of ENNs on runoff modeling

In general, all hybrid semi-distributed models exhibit notable performance, adeptly capturing the runoff peaks with appropriate magnitudes and timings across three study basins (Fig. 3 and Table 2). Specifically, the comparison results show that the DM_θ model exhibits a closed but slightly better performance than the DM model in overall runoff modeling, with a slight increase in NSE and mNSE of 0.01–0.03 in all three basins. Additionally, lower PFAB results imply that the DM_θ model contributes to an improved performance in peak runoff modeling. The incorporation of ENNs to represent the spatial heterogeneity of calibration parameters can reduce the peak simulation biases and slightly improve the overall performance.

Figure 3 The comparison of simulated (DM, DM_θ, DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models) and observed runoff processes in the evaluation period at the trained TNH, ZMD and JZ station in Yellow, Yangtze and Lancang, respectively.

Table 2 The results of three hydrological metrics for different hybrid semi-distributed models in the three study basins.

The notably enhanced performance in DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models indicates that the inclusion of ENNs for replacing internal modules yields further improvements in model performance (Fig. 3 and Table 2). First, the results between DM${}_{\mathit{\theta }-Q-T}$ and DM_θ models show the significant improvement in runoff modeling brought by the incorporation of ENN_Q. This enhancement is illustrated by an increase in NSE and mNSE values, ranging from 0.06 to 0.09 in Yellow and Yangtze. Since the DM_θ model already exhibits commendable performance in Lancang, the advancements achieved by the DM${}_{\mathit{\theta }-Q-T}$ model are relatively marginal in comparison. PFAB results suggest that the ENN_Q does not lead to substantial improvements in peak flow performance. In addition, evaluation findings for the DM${}_{\mathit{\theta }-QSM-T}$ model show that replacing precipitation partition and snowmelt modules by ENNs can further improve the model performance with an increase NSE of 0.01–0.05. It also does not translate into better peak runoff modeling, as evidenced by comparable PFAB scores across all three basins. ENNs employed for replacement in hybrid hydrological models have proven to be effective in enhancing the model performance in runoff modeling. Among them, the ENN_Q leads to the most substantial improvements in runoff prediction performance. The replacement of ENNs for snow-related processes (ENN_S and ENN_M) results in comparatively minor enhancements. These findings align with our hydrological understanding as the runoff module directly generates runoff and thus plays a central role in runoff modeling. It thus contributes the most to the overall performance of runoff prediction. Conversely, the influence of snow-related processes on runoff modeling performance improvements is indirect and thus relatively modest (Li et al., 2023b).

The air temperature is employed as the additional input of the ENN_Q to implicitly represent the soil freeze–thaw process in this study (Zhong et al., 2023; Gao et al., 2021). Results indicate that DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models exhibit improved performance in peaking runoff modeling compared to the DM_θ−Q and DM_θ−QSM models, respectively. This enhancement in peaking runoff modeling is evident through closed NSE and mNSE and lower PFAB values in all three basins. Moreover, the enhancement observed due to the inclusion of air temperature is notably more pronounced in Yellow and Yangtze compared to Lancang. This pattern aligns with expectations because Lancang features a smaller extent of permafrost regions, resulting in a smaller influence of the soil freeze–thaw process on runoff modeling in this region.

3.1.2 The impact of spatial information on runoff modeling

Hybrid lumped models proposed by Li et al. (2023b) are similar to our proposed hybrid semi-distributed models but did not consider the spatial heterogeneity. Hybrid lumped and semi-distributed models are used to test the effect of spatial information on hydrological modeling. It is important to highlight that while the ENNs of the hybrid lumped models utilize the same dynamic time series inputs as those of the distributed models, they do not include the static attributes of the basin. Results show that both hybrid lumped models, EXP_Q and EXP_QSM, exhibit strong performance in runoff modeling with NSE more than 0.74 in all three basins (Fig. 4 and Table 3). It demonstrated the suitability of hybrid lumped models for hydrological modeling on the TP. In comparison to EXP_Q and EXP_QSM models, the DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models show more impressive performance in runoff modeling with the increase NSE and mNSE of 0.01–0.14 in three basins. PFAB results affirm that DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models excel in simulating peak flow processes, achieving PFAB values of less than 10 % across all three basins. Consequently, the incorporation of spatial heterogeneity within the basin into hybrid models leads to improved performance in both overall and peak runoff modeling. This finding is seamlessly consistent with our hydrological understanding and is also corroborated by related studies in the case of distributed process-based hydrological models and DL hydrological models (Li et al., 2023a; Patil et al., 2014). In practice, we recommend the utilization of hybrid semi-distributed models for hydrological modeling, particularly in the context of large basins, to attain enhanced performance outcomes.

Figure 4 (a–c) The comparison of simulated and observed runoff processes in the evaluation period in Yellow, Yangtze and Lancang, respectively. DM_T and EXP are denoted to hybrid semi-distributed and lumped models, while DM represents the hybrid semi-distributed models without inclusion of air temperature in ENN_Q. Circles, squares and triangles refer to NSE, mNSE and PFAB. (d) The model comparison with state-of-the-art models.

Table 3 The results of three hydrological metrics for different hybrid semi-distributed and lumped models in the three study basins.

3.1.3 The comparison to the state-of-the-art models

We further use the optimal hybrid semi-distributed model DM${}_{\mathit{\theta }-QSM-T}$ to compare with state-of-the-art models: the distributed hydrological model THREW and the DL models LSTM and CNN-LSTM (Li et al., 2023a). Results show that the DM${}_{\mathit{\theta }-QSM-T}$ model outperforms the THREW model by a substantial margin and holds comparable performance to the LSTM and CNN-LSTM models (Fig. 4). This reveals that our hybrid semi-distributed model can effectively harness the advantages of both process-based models and DL models. Specifically, it attains the high performance characteristic of DL models while adhering to the physical mechanism constraints inherent in process-based models, creating a synergy not entirely realized in other models.

3.2 Hybrid semi-distributed model evaluation in untrained sites within the basin

As proposed hybrid models operate in a semi-distributed manner, it is imperative to further investigate whether models trained using the basin outlet point can effectively simulate hydrological processes in any untrained sites within the same basin. In this study, runoff processes at three hydrological stations (JG, MQ and MT), situated upstream of TNH in Yellow, are simulated using our proposed hybrid models trained by TNH data (Figs. 2 and 5).

Figure 5 The NSE results between simulated (different hydrological models in TNH) and observed runoff processes at JG, MQ and MT. e represents the static attributes of sub-basins normalized by the maximum–minimum in TNH.

Results reveal that all models trained on TNH data exhibit impressive performance in simulating runoff processes at JG and MQ stations, with NSE values exceeding 0.71. The DM${}_{\mathit{\theta }-QSM-T}$ model achieves an especially high NSE of 0.84. However, the models demonstrates lower accuracy at the station that is most upstream, MT (Fig. 2). This is because the alpine hydrological processes in the basin above the MT station, such as soil freeze–thaw and snow and ice processes, play a more significant role in runoff processes (Fig. 5e). This increases the difficulty of hydrological simulation, leading to reduced model accuracy. Among them, DM and DM_θ models show the most significant reduction in accuracy due to their insufficient representation of alpine hydrological processes. On the other hand, hybrid semi-distributed models with ENN replacement, including DM${}_{\mathit{\theta }-Q-T}$ and DM${}_{\mathit{\theta }-QSM-T}$ models, exhibit notably enhanced abilities in runoff modeling compared to DM and DM_θ models, resulting in NSE improvements ranging from 0.09 to 0.58. The DM${}_{\mathit{\theta }-QSM-T}$ model demonstrates the strongest performance in runoff modeling across all three stations, particularly in MT where its NSE reaches 0.54, whereas the other three models yield NSE values lower than 0.22 (Fig. 5). The findings show that the proposed hybrid semi-distributed models exhibit strong performance in hydrological modeling for untrained sites within the basin. It is also demonstrated that the hydrological relationships established by ENNs are credible and robust.

4 The applicability of hybrid models for hydrological sensitivities to climate change

Perturbed precipitation and air temperature dataset are input to trained DM${}_{\mathit{\theta }-QSM-T}$ models to test the applicability of the hybrid models to analyze the hydrological sensitivities to climate change in three large alpine basins.

Figure 6 Runoff responses to altered precipitation in the Yellow, Yangtze and Lancang basins (a–c for annual; d–f for monthly). The error bars in panels (a)–(c) and the shaded areas in panels (d)–(f) denote the range of simulated runoff.

Figure 7 Relative change of annual (grey background) and monthly(yellow background) runoff response to the perturbed precipitation (a–c) and air temperature (d–f) in Yellow, Yangtze and Lancang, respectively.

4.1 Sensitivities of runoff to perturbed precipitation

Figures 6 and 7a–c depict the runoff sensitivities to various altered precipitation scenarios within three study basins. The findings suggest a consistent trend in the relationship between runoff and precipitation: runoff rises (decreases) as precipitation increases (decreases). Specifically, the annual runoff increases at rates of approximately 33.8, 18.1 and 44.9 mm per 10 % with the increase of precipitation within Yellow, Yangtze and Lancang, respectively. The relative change in runoff surpasses that of precipitation in all three study basins: a 10 % increase in precipitation leads to a 15 % to 20 % increase in runoff in all three study basins. In addition, annual runoff exhibits greater sensitivity to increases in precipitation compared to decreases (Fig. 7a–c). As an illustration, an increase of 20 % in precipitation results in a substantial 40 % increase in annual runoff, whereas a 20 % decrease in precipitation leads to a notable 30 % reduction in annual runoff in Yellow. It is indicated that runoff exhibits an amplification effect in response to precipitation changes due to the increase in the runoff coefficient with rising precipitation. Figure 6a–c also illustrate that the inter-annual variation in runoff follows a pattern consistent with the annual runoff: there is a greater (lesser) variation in inter-annual runoff when there is an increase (decrease) in precipitation.

Figure 8 Runoff responses to altered temperature in the Yellow, Yangtze and Lancang basins (a–c for annual; d–f for monthly). The error bars in panels (a)–(c) and the shaded areas in panels (d)–(f) denote the range of simulated runoff.

Moreover, the monthly runoff across all months shows a consistent response to perturbed precipitation, yet the extent of change varied among different months (Figs. 6d–f and 7a–c). Notably, the alterations during the wet seasons (June to October) are more pronounced compared to those in the dry seasons. This indicates that increased precipitation contributes to a more concentrated distribution of runoff. Figure 6d–f also demonstrate that intra-annual runoff variation becomes more pronounced with higher levels of precipitation. These findings can be attributed to the fact that the augmented precipitation primarily occurs during the wet seasons, and the primary runoff components during these periods consist of direct rainfall runoff.

4.2 Sensitivities of runoff to perturbed temperature

The sensitivities of runoff to changing temperature follow a more intricate pattern: runoff tends to decrease as temperatures rise. This decrease is particularly pronounced during the flood season, while in the dry season, there is a slight increase in runoff (Figs. 7d–f, 8 and 9). This shift also leads to a reduction in the intra-annual variability of runoff. Taking the temperature increase of 2 °C as an example, the annual runoff in the three study basins decreases by less than 15 %. When examining monthly runoff, the most significant increase occurs in April, while the most notable decrease is observed in June. These phenomena can be explained by the fact that changes in temperature affect the evaporation capacity, the redistribution of rainfall and snowfall, and the timing of snowmelt. Higher temperature leads to increased evaporation capability, which results in more actual evaporation and less total runoff when precipitation remains constant. During the winter and spring, the increased rainfall and earlier snowmelt, along with higher actual evaporation, tend to balance each other, resulting in a minor increase or decrease in runoff. However, in the summer, reduced snowmelt and higher evaporation significantly reduce runoff.

Figure 9 The runoff components with different perturbed temperature scenarios in Yellow, Yangtze and Lancang, respectively.

To enhance the reliability of our model and validate our findings of hydrological sensitivities to climate change, we conducted an analysis of runoff component contributions in all three study basins across scenarios with varying temperature perturbations. It is essential to highlight that the glacier module has been excluded from this model due to structural limitations. Previous studies in the study basins have demonstrated that glaciers have a negligible impact on runoff (Cui et al., 2023; Su et al., 2023). As a result, this limitation does not significantly affect the accuracy of the simulation results. In the reference scenario, rainfall runoff emerged as the primary component, contributing approximately 81.5 %, 73.1 % and 84.0 % to the total runoff in Yellow, Yangtze and Lancang, respectively. Notably, these results align with findings from other studies (Cui et al., 2023; Su et al., 2023), underscoring that our hybrid model not only excels in simulating the runoff process but also accurately represents untrained hydrological processes. Furthermore, the contribution of snowfall runoff diminishes as the perturbed temperature increases. With a 2 °C temperature rise, the contribution of snowfall runoff decreases by 5.8 %, 8.9 % and 5.0 % in the Yellow, Yangtze and Lancang basins, respectively. These results strongly support the credibility of our analysis.

5 Conclusions and limitations

In this study, we propose hybrid semi-distributed hydrological models that synergize the semi-distributed process-based model with embedded neural networks (ENNs). The hybrid models use the semi-distributed process-based model as the backbone, with ENNs parameterizing and replacing internal modules. Taking three large alpine basins on the Tibetan Plateau as the study basins, the proposed models are tested and compared with state-of-the-art models. The climate perturbation method is further carried out to test the applicability of the hybrid models to analyze the hydrological sensitivities to climate change in large alpine basins. Our main findings are as follows:

1. The optimal hybrid semi-distributed model achieves superior performance in runoff modeling, with an NSE higher than 0.87, approaching the state-of-the-art DL models and outperforming traditional process-based models. The optimal hybrid semi-distributed model also demonstrates remarkable prowess in hydrological modeling at ungauged sites within the basin.

2. Further experiments reveal that the inclusions of ENNs for parameterizing and replacing modules can lead to higher model accuracy. Considering spatial information within the basin and introducing temperature in ENN_Q to represent the soil freeze–thaw process also show enhanced predictive capabilities in hybrid models.

3. The results about hydrological sensitivities to climate change show reasonable patterns: runoff exhibits an amplification effect in response to precipitation changes, with a 10 % precipitation change resulting in a 15 %–20 % runoff change in large alpine basins. Annual runoff exhibits greater sensitivity to increases in precipitation compared to decreases. The increase in temperature enhances evaporation capacity and reduces the contributions of snowfall runoff, leading to a decrease in the total runoff and a reduction in the intra-annual variability of runoff. With a 2 °C temperature rise, the contribution of snowfall runoff decreases by 5.8 %, 8.9 % and 5.0 % in the Yellow, Yangtze and Lancang basins, respectively.

In summary, we provide an effective and easily interpretable hybrid semi-distributed hydrological model and enhance our understanding about hydrological sensitivities to climate change in large alpine basins. However, being promising in modeling hydrological processes, this study also has several limitations. First, the routing method is important for hydrological modeling, especially in large basins. The technical requirements of the differential programming framework limit the consideration of routing methods in our hybrid hydrological models. We calculate the river length from each sub-basin to the basin outlet and employ this static attribute as the inputs of ENNs to implicitly characterize the routing process within the basin. In addition, this study is limited to only using three large alpine basins on the Tibetan Plateau to evaluate proposed hybrid models due to the limitation of computational resources. Third, although numerous studies have used a climate perturbation method to calculate the response of hydrological processes to climate change, this approach has difficulty capturing the true characteristics of meteorological and hydrological changes, making it hard to validate the reasonableness of the results. In this study, we compared our findings with those of related research to demonstrate the validity of our results, thereby proving the effectiveness of our proposed coupled model in analyzing the response of hydrological processes to climate change. Future research will focus on developing hybrid distributed hydrological models, including routing processes and extending the evaluation of the hybrid model to encompass a broader range of basins.

Appendix A

A1 Distributed EXP-Hydro model equations

The semi-distributed EXP-Hydro model firstly delineates the basin into many sub-basins. In each sub-basin, the lumped EXP-Hydro is run independently (Eqs. A1–A12) to obtain the respective runoff. The runoff from all sub-basins is then aggregated to calculate the basin runoff (Eq. A12). The detailed equations are as follows (Patil et al., 2014).

1. Water balance is

   $$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{0}}}{\mathrm{d}t}}=P_{\mathrm{s}}-M\text{(A2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{1}}}{\mathrm{d}t}}=P_{\mathrm{r}}+M-\mathrm{ET}-Q,\end{array}$$

   where S₀, S₁, P_s, P_r, M, ET and Q are snow storage, basin water storage, snowfall, rainfall, snowmelt, evaporation and runoff, respectively.

2. Precipitation partition is

   $$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{s}}=\left\{\begin{array}{cl}\mathrm{0}& T>T_{\min }\\ P& T\le T_{\min }\end{array}\right.\text{(A4)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{r}}=\left\{\begin{array}{lc}P& T>T_{\min }\\ \mathrm{0}& T\le T_{\min }\end{array}\right.,\end{array}$$

   where P and T are precipitation and air temperature.

3. Snowmelt is

   $$\begin{array}{}\text{(A5)}& M=\left\{\begin{array}{ll}\min \left\{S_{\mathrm{0}},D_{\mathrm{f}}\cdot \left(T-T_{\max }\right)\right\}& T>T_{\max }\\ \mathrm{0}& T\le T_{\max }\end{array}\right..\end{array}$$

4. Evapotranspiration is

   $$\begin{array}{}\text{(A6)}& {\displaystyle }& {\displaystyle }\mathrm{ET}=\left\{\begin{array}{ll}\mathrm{0}& S_{\mathrm{1}}<\mathrm{0}\\ \mathrm{PET}\cdot \left(\frac{S_{\mathrm{1}}}{S_{\max }}\right)& \mathrm{0}\le S_{\mathrm{1}}\le S_{\max }\\ \mathrm{PET}& S_{\mathrm{1}}>S_{\max }\end{array}\right.\text{(A7)}& {\displaystyle }& {\displaystyle }\mathrm{PET}=\mathrm{29.8}{L}_{\mathrm{day}}{\displaystyle \frac{e_{\mathrm{sat}}\left(T\right)}{T+\mathrm{237.3}}}\text{(A8)}& {\displaystyle }& {\displaystyle }{e}_{\mathrm{sat}}\left(T\right)=\mathrm{0.611}\times \exp \left({\displaystyle \frac{\mathrm{17.3}T}{T+\mathrm{237.3}}}\right),\end{array}$$

   where PET, L_day and e_sat(T) represent the potential evaporation, the day length and the saturation vapor pressure.

5. Runoff and baseflow are

   $$\begin{array}{}\text{(A9)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{b}}=\left\{\begin{array}{ll}\mathrm{0}& S_{\mathrm{1}}<\mathrm{0}\\ Q_{\max }\cdot e^{-f\cdot \left(S_{\max }-S_{\mathrm{1}}\right)}& \mathrm{0}\le S_{\mathrm{1}}\le S_{\max }\\ Q_{\max }& S_{\mathrm{1}}>S_{\max }\end{array}\right.\text{(A10)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{s}}=\left\{\begin{array}{ll}\mathrm{0}& S_{\mathrm{1}}\le S_{\max }\\ S_{\mathrm{1}}-S_{\max }& S_{\mathrm{1}}>S_{\max }\end{array}\right.\text{(A11)}& {\displaystyle }& {\displaystyle }Q=Q_{\mathrm{b}}+Q_{\mathrm{s}},\end{array}$$

   where Q_b and Q_s are the baseflow generated depending on the available storage in the basin bucket and the capacity-excess runoff generated when the basin bucket is saturated. All the above undefined variables are calibration parameters. For details, please refer to Patil and Stieglitz (2014).

6. Basin runoff is

   $$\begin{array}{}\text{(A12)}& Q_{\mathrm{basin}}={\displaystyle \frac{\sum _{i=\mathrm{1}}^N{Q}_i\cdot A_i}{\sum _{i=\mathrm{1}}^N{A}_i}},\end{array}$$

   where Q_basin is the runoff at basin outlet. Q_i and A_i are the runoff and area of sub-basin i. N is the total number of sub-basins within the basin.

A2 Hybrid semi-distributed model equations

In all hybrid semi-distributed models, four ENNs are constructed to parameterize (NN_θ) and replace runoff (NN_Q), precipitation partition (NN_S) and snowmelt processes (NN_M). The detailed equations are as follows.

$$\begin{array}{}\text{(A13)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{d}}={\mathrm{NN}}_{\mathit{\theta }}\left(A_{\mathrm{s}}\right)\text{(A14)}& {\displaystyle }& {\displaystyle }Q={\mathrm{NN}}_Q\left(M+P_{\mathrm{r}},S_{\mathrm{1}},T,A_{\mathrm{s}}\right)\text{(A15)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{s}}=P\times {\mathrm{NN}}_S\left(P,T,A_{\mathrm{s}}\right)\text{(A16)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{r}}=P-P_{\mathrm{s}}\text{(A17)}& {\displaystyle }& {\displaystyle }M=S_{\mathrm{0}}\times {\mathrm{NN}}_M\left(T,A_{\mathrm{s}}\right),\end{array}$$

where θ_d and A_s represent calibration parameters and static basin attributes, respectively. For detailed static basin attributes, refer to Table A1.

Figure A1 The spatial heterogeneity of precipitation and air temperature in Yellow, Yangtze and Lancang.

Table A1 The summary of static basin attributes for the inputs of ENNs.