Coupling the Glacio-hydrological Degree-day Model (GDM) with PCRaster for spatial dynamic modeling of Himalayan river basins

We have seen a surge in glacio-hydrological modeling efforts in the past few decades. This form of modeling is also being carried out in the Himalayan river basins, but a comprehensive high-resolution simulation software that can be effective with a limited number of hydrometeorological data is recommended. In this regard, an open-source, scalable, flexible, and distributed modeling system called PyGDM has been developed by fully coupling the Glacio-hydrological Degree-day Model with PCRaster Python software. To evaluate the potential of using this model in the Himalayan river basins, we calibrated and then validated the model for Trishuli River basin using geographical data and the existing hydrometeorological data. The tests showed a promising result with respect to the effective application of the model in the entire Himalayan region. The PyGDM source code was optimized and adapted to the process models of glacier melting and hydrological processes in Himalayan basins. It increased the speed of the simulation, made the model highly scalable to accommodate new submodels, and enhanced the flexibility of the model to ingest various types of input data and parameters. Hence, the PyGDM model strives to simulate the glacio-hydrological processes of the entire Himalayan region.

1 Introduction

Modeling glacio-hydrological processes started with different modeling approaches (e.g., Buch et al., 1993; Kayastha et al., 1999) and has gained momentum, with hundreds of studies in Asia and other parts of the world focusing on glacial melting process and the related calibration and validation processes of modeling. (see the review in van Tiel et al., 2020; Y. Chen et al., 2017).

There are mainly stochastic, conceptual, and physics-based models used in semi-distributed and distributed modes (Hock and Jansson, 2006). We see some novel use of artificial neural network (ANN; Buch et al., 1993) machine learning techniques, such as long short-term memory (LSTM) and extreme gradient boosting (XGBoost) (Ji et al., 2021) in Central Asia and gated recurrent units (GRUs) integrated with the revised version of the Soil and Water Assessment Tool (SWAT+) that considers glacial hydrological processes (GSWAT+; Yang et al., 2023). Nevertheless, the ease of model development and the need for only a few types of data have given us more conceptual and semi-distributed models.

To give just a few examples of the variety of glacio-hydrological models that have been developed for various river basins around the world, the variable infiltration capacity glacier (VIC-glacier) model was used in the Eastern Pamirs (Ren et al., 2018) and China (Zhao et al., 2013); an energy-balance model (Zhang et al., 2016), the Coupled Routing and Excess Storage (CREST-glacier) model (X. Chen et al., 2017), a tracer-aided hydrological model (Nan et al., 2022), and the SWAT+ degree-day model were employed on the Tibetan Plateau (Adnan et al., 2019); the SWAT+ degree-day model (Rahman et al., 2013) and the HBV-light model (Finger et al., 2015) were used in Switzerland; the PRMSglacier model was employed in Alaska (Van Beusekom and Viger, 2016, 2018); glacier melt models were used in the Andes (Kinouchi et al., 2013), Northwestern China (Ji and Luo, 2013), and the Upper Indus Basin (Baig et al., 2021); and the WRF-Hydro/Glacier model was employed in southern central Norway (Eidhammer et al., 2021).

The Himalayan region has also required its own set of modeling efforts (Azam et al., 2021). This has involved a mix of different modeling approaches, such as ANN models (Buch et al., 1993), the wasm overflow detector (WASMOD) model with precipitation from the Weather Research and Forecasting (WRF) model in the Beas River basin (Li et al., 2019), the VIC-glacier model in five river basins (Chandel and Ghosh, 2021), the coupled VIC–glacier melt (GM2) model (VIC-GM2) in the Tarim River (Wang et al., 2021), an energy-balance-based distributed glacier model in the Hunza River basin (Shrestha et al., 2015), a snow and glacier melt model in the Upper Indus Basin (Gupta et al., 2021), the SWAT model in the Western Himalayas (Shukla et al., 2021) and Hunza (Garee et al., 2017), a distributed hydrological model in the Western Himalayas (Wulf et al., 2016), and the snowmelt model (SNOWMOD) for the Gangotri Glacier (Arora et al., 2016).

In Nepalese Himalayan river basins, an integrated hydro-glaciological modeling system is proposed so that a common platform is available. This system has the advantage of a distributed model but only requires a few input data and parameters. This paper describes the development and testing of this model that integrates the Glacio-hydrological Degree-day Model (GDM; Kayastha et al., 2020) with PCRaster.

2 Model development

There is the need for a distributed modeling approach that is simple yet powerful and flexible to handle the vast number of spatial and temporal data generated. We chose PCRaster Python as the data processing and simulation software, as it has been successfully used in other modeling frameworks such as the Spatial Processes in Hydrology (SPHY) model (Terink et al., 2015), the PCRaster Global Water Balance 2 (PCR-GLOBWB 2) model (Sutanudjaja et al., 2017; Wiersma et al., 2022), and other similar hydrological models (Karssenberg et al., 2015; López et al., 2019).

We developed the PCRaster Python Glacio-hydrological Degree-day Model (PyGDM) by coupling the PCRaster Python framework for spatial dynamic modeling (Karssenberg et al., 2010) with the Glacio-hydrological Degree-day Model (originally developed in Fortran and then in C#) (Kayastha et al., 2005, 2020).

Our PyGDM model is a spatially distributed and parsimonious glacio-hydrological model especially designed for application in Himalayan river basins. The simplicity of the GDM is coupled with the power and ease of the dynamic modeling of PCRaster by embedding the GDM code into the initial and dynamic sections of the PCRaster script.

2.1 Modeling framework

The modeling method is shown in Fig. 1. Kayastha et al. (2003) proposed positive degree-day factors for the Himalayan river basins, and this method has become a simple and convenient approach to simulate glacierized basin discharge in the Himalayan basins (e.g., Pradhananga et al., 2014; Silwal et al., 2016; Kayastha et al., 2020).

Figure 1 PyGDM modeling framework. The abbreviations used in the figure are as follows: GCM – global circulation model; RCM – regional climate model; and SSP – Shared Socioeconomic Pathway.

Precipitation is calculated as rain or snow according the threshold temperature (Eq. 1):

$$\begin{array}{}\text{(1)}& \mathrm{Precipitation}=\left\{\begin{array}{ll}\mathrm{rain}& T\ge T_{\mathrm{T}}\\ \mathrm{snow}& T<T_{\mathrm{T}},\end{array}\right.\end{array}$$

where T is temperature and T_T is the critical/threshold temperature.

The snowmelt and ice melt value, M (mm d⁻¹), is calculated using degree-day factors for snow, bare ice, and debris-covered ice:

$$\begin{array}{}\text{(2)}& M=\left\{\begin{array}{ll}{K}_{\mathrm{d}}\text{ or }{K}_{\mathrm{s}}\text{ or }{K}_{\mathrm{b}}\times T& T>\mathrm{0}\\ \mathrm{0}& T\le \mathrm{0},\end{array}\right.\end{array}$$

where T is air temperature (°C), K_s is the degree-day factor for snow (mm °C d⁻¹), K_b is the degree-day factor for bare ice (mm °C d⁻¹), and K_d is the degree-day factor for debris-covered ice (mm °C d⁻¹).

The Hargreaves method is a very simple empirical method for estimating potential evapotranspiration (PET) based on temperature data. The advantage of using this method is that we need only temperature data as input; this enables us to generalize the method for the entire Nepalese domain. The Hargreaves method estimates PET based on Eq. (3):

$$\begin{array}{}\text{(3)}& \mathrm{PET}=\mathrm{0.0023}(T_{\mathrm{mean}}+\mathrm{17.8})R_a\mathrm{0.0408}\sqrt{T_{\max }-T_{\min }},\end{array}$$

where PET is potential evapotranspiration (mm d⁻¹); T_mean is average daily temperature (°C); T_max is maximum daily temperature (°C); T_min is minimum daily temperature (°C); and R_a is extraterrestrial radiation (MJ m² d⁻¹), which can be calculated based on location, date, and time using established formulas.

The Hargreaves method is also easy to use due to its dependence on only temperature and radiation, although we may consider other variables such as humidity and wind in future versions of the model. A simple interception model was used with monthly threshold values provided in the model.

The two-reservoir-based modeling approach of SWAT (Luo et al., 2012) is adopted to simulate the hydrological response of the baseflow and rainfall runoff contribution to river discharge. Then, the discharge due to rainfall and snowfall are calculated with the use of runoff coefficients, and the daily accumulated amount of water flowing through the river basin is simulated using the partial contribution from the past-day discharge (Eqs. 4, 5). The details of the discharge simulation are explained in more precisely in Kayastha et al. (2020).

We decided to utilize the static and dynamic modeling capabilities of PCRaster software to simulate the hydrological response of the Himalayan river basins. PCRaster is a highly capable environmental modeling software with a Python version that is open-source and free. The deterministic modeling required for the distributed hydrological modeling starts with the static modeling framework, as summarized by Eq. (6). In Eq. (6), Z is the model state variable, I represents the inputs, P represents the parameters, and f defines the model structure. In the dynamic modeling section of the PCRaster script, a time component is added in the static model. In the dynamic modeling framework, there are several sections such as the “initial” section, the “dynamic” section, the dynamic database, the timer, and so on. The initial section is used to initialize static values and parameters used in the hydrological model, to store lookup tables of various parameters, and to define the major geographical databases, such as digital elevation model, land use and land cover, and basin spatial maps. The dynamic section is used to apply the temporal behavior of the hydrological model by running the operations on maps for every time step in sequence. PCRaster also has a multicore module to facilitate the parallel running mode. We used the multicore capability in the workstation to decrease the simulation time.

$$\begin{array}{}\text{(4)}& {\displaystyle }{Q}_G=Q_r\cdot C_r+Q_s\cdot C_s\text{(5)}& {\displaystyle }{Q}_d=Q_R\left(\mathrm{1}-k\right)+Q_{R(d-\mathrm{1})}k+Q_B\text{(6)}& {\displaystyle }Z=f(Z,I,P)\end{array}$$

The only spatially and temporally varying input data required for our model are the precipitation and temperature. Our model is able to interpolate the station data as well as use the gridded NetCDF data files that may be obtained from global and regional models and reanalysis data.

2.2 Input data

In the present version of our model, land use and land cover are classified into eight categories: forest, grassland, agriculture, barren land, waterbodies, built-up area, debris-covered glacier, and clean/bare glacier.

To calibrate and validate PyGDM, we applied the model to the Trishuli River basin, which covers part of China and Nepal. We used the gridded land use and land cover data (Fig. 2), which were obtained by modifying and reclassifying the Esri LULC (land use and land cover) derived from Sentinel-2 imagery at a 10 m resolution (Karra et al., 2021) (https://livingatlas.arcgis.com/landcoverexplorer, last access: 10 June 2024). To reclassify land classes for glacier, the Randolph Glacier Inventory Version 6.0 (RGI v6.0) glacier outlines were used. A 30 m resolution digital elevation model (DEM) was obtained from the Shuttle Radar Topography Mission (SRTM) 1 arcsec Global 30 m DEM from the United States Geological Survey (https://earthexplorer.usgs.gov, last access: 10 June 2024).

Figure 2 Land use and land cover (LULC) of Trishuli River basin. We used the LULC from Esri and modified it as per the requirements of our model. The Esri LULC was produced by Impact Observatory, Microsoft, and Esri (Karra et al., 2021). Publisher's remark: please note that the above figure contains disputed territories.

2.3 Parameters

Table 1 shows the main parameters used in the discharge simulation, baseflow calculation of shallow and deep aquifers, and critical temperature. All of the parameters have been identified from our past research (Kayastha et al., 2020).

Table 1 Parameters and coefficients for modeling.

Table 2 Degree-day factors and runoff coefficients.

Table 2 shows the degree-day factors, runoff coefficients, and interception thresholds for each of the land use categories used in the model. These factors and coefficients are taken from our research as well as past publications (Kayastha et al., 2005, 2020).

The interception is dependent on the land use and is modeled as a threshold value. The interception thresholds vary per land use type and have been adopted from the guidelines provided by Gerrits (2010).

3 Results and discussion

The PyGDM model ran successfully in multicore mode with highly acceptable results. The simulation time was below 10 min for a 1-year simulation at a 30 m grid resolution. The simulation time period may be altered by the actual computational power, RAM, and the CPU core characteristics.

A total of 4 years of simulation was used for the calibration and validation of the model, respectively. The calibration and validation results of the model are shown in Fig. 3. The Nash–Sutcliffe efficiency and the volume difference were the statistical metrics used in the evaluation. Both of the metrics are within the acceptable limits of the performance of the models. PyGDM is also capable of evaluating the component contribution of hydrological processes, such as baseflow, rain, snowmelt, and ice melt, on the total discharge of the basin. The contribution of process components to the discharge of Trishuli River basin is shown in Fig. 4. Although the spatial variation is not presented in this paper, the PyGDM modeling system is able to plot the spatially distributed maps of the output variables as well as animations depicting the temporal development of the hydrological response of the basin. The calibration and validation results are within the acceptable regime. The speed of simulation in multicore mode also promises very effective application of this model in multiyear and decadal simulations of the hydrological response and climate change scenarios in the Himalayan river basins. PyGDM is also capable of ingesting a diverse set of input data, and the use of the Python language gives an easy platform to upscale the model with respect to the use of new subprocesses.

Figure 3 Calibration and validation of discharge for the Trishuli River basin.

Figure 4 Components of discharge for the Trishuli River basin.

4 Conclusions

A glacio-hydrological model built on positive degree-day factors for glacier melt modeling was successfully integrated with PCRaster Python software to enable a distributed approach and enhance the scalability of the conceptual model for application in the entire Himalayan region. The dynamic modeling capability of PCRaster has been coupled with the simple yet highly effective glacio-hydrological model so that the spatiotemporal resolution of the simulation could be enhanced. To test this newly developed PyGDM model, calibration and validation phases of model development were successfully completed for the Trishuli River basin. PyGDM is ready for the comprehensive modeling and simulation of glacio-hydrological processes in the Himalayan basins. The flexibility, scalability, and the simplicity of using available environmental data in the Himalayan region are the main advantages of PyGDM.