Multi-Objective Validation of the Runoff Formation Model in the High-Mountain River Basin of the Central Caucasus

INTRODUCTION

Assessing the spatial and temporal variability of water balance components using distributed hydrological models at a regional scale is crucial for addressing practical challenges in managing water resource, including evaluating the impact of climate variability on water resources. Such studies require accurate model calibration and validation, which are typically performed using a traditional method based on observed discharge data. Previous research has highlighted that process-based multi-objective calibration and validation effectively reduce uncertainty and improve forecasting accuracy in hydrological modeling (Efstratiadis and Koutsoyiannis 2010). Many studies have incorporated additional data sources to enhance model validation. Remotely sensed hydrological variables, such as actual evapotranspiration (Immerzeel and Droogers 2008; Zhang et al. 2009), soil moisture (Li et al. 2018; Parajka et al. 2009; Poméon et al. 2018; Wanders et al. 2014), snow cover (Duethmann et al. 2014; Han et al. 2019; Tong et al. 2021; Udnæs et al. 2007), and total water storage (Bai et al. 2018; Trautmann et al. 2018) have been frequently used to assess model performance. Stable isotope tracer data also offer valuable insights for evaluating hydrological models by providing information on water flow pathways and sources (Ala-Aho et al. 2017; Holmes et al. 2020; Holmes et al. 2023).

However, multi-objective validation is commonly used in lowland areas because of the availability of extensive data from monitoring stations, simpler terrain, and more predictable hydrological processes. Mountainous areas present challenges for multi-objective validation. The complex terrain, steep slopes, and varying land cover create heterogeneous hydrological processes that are more complicated for modeling. Additionally, data availability is often limited with fewer monitoring stations. This makes validation more challenging than in the lowland regions. Assessment of snow cover is crucial for more complex hydrological modeling and is perhaps one of the most important unresolved problems for modeling in mountainous areas (Dozier et al. 2016; Lettenmaier et al. 2015; Treichler and Kääb 2017). This problem becomes more pronounced in the case of complex distribution mechanisms for solid precipitation due to differences in moisture and wind conditions in complex high-mountainous topographies (Mikhalenko et al. 2020). Typically, spatially distributed snow information over large areas is the result of a combination of modeling or interpolation based on sparse in situ measurements, remote sensing data, or climate reanalysis (Treichler and Kääb 2017). The approach of determining snowpack thickness was presented in a recent study (N. E. Elagina et al. 2024), where the glacial runoff model A-Melt (Rets et al. 2021) in combination with geodetic estimates for Mount Elbrus (North Caucasus) was used (Kutuzov et al. 2019) (Kutuzov et al., 2019).

This study focuses on the significantly glaciated Baksan River basin in the North Caucasus. Snow and glaciers provide a stable water supply throughout the year, which supports drinking water, irrigation, and hydropower generation. Modeling of runoff processes provides detailed knowledge of the glacier and snow-fed components of runoff, which is essential for assessing water availability and planning for future climate impacts. Validating hydrological models ensures reliable forecasts of water availability, helps to manage water resources, reduces flood risks, and sustains essential sectors in this mountainous region.

One of the models that has proven to be effective in high-altitude river runoff modeling is the ECOMAG model (ECOlogical Model for Applied Geophysics) (Motovilov 1999), which has been applied to understand the behavior of hydrological systems and explain the impacts of different climate and land-use scenarios across a wide range of geographic regions (Moreido and Kalugin 2017; Kalugin and Motovilov 2018), including mountainous areas with present-day glaciers and permafrost (Gelfan et al. 2017; Kornilova et al. 2024; Kornilova et al. 2021; Lisina et al. 2023; Motovilov et al. 2017). Various datasets were used to test and validate the ECOMAG model in lowland basins to ensure its performance across multiple hydrological components. Simulation results based on data from the NOPEX experiment in Sweden (Motovilov 1999) showed that the model effectively captured spatial and temporal variations in water balance components, including soil moisture and groundwater levels. In the process of step-by-step calibration of parameters and modeling of snow cover and soil moisture fields in the Volga Basin, the actual fields of these characteristics were constructed based on the data of en-route snow measurement observations and productive moisture measurements at agrometeorological stations (Motovilov and Gelfan 2018). Using a modification of the ECOMAG-HM model for the Kama River basin, modeled copper concentrations were compared on a daily basis with episodic measurements at two hydrochemical monitoring stations (Motovilov and Fashchevskaya 2019). However, the model has been specifically adapted for the first time (Kornilova et al. 2023; Kornilova et al. 2024) to provide a detailed multi-objective validation in an upland basin.

The aim of this study is to consider the possibilities and limitations of multi-objective validation of a distributed hydrological model in a high-mountain river basin. For the multi-objective validation of the ECOMAG hydrological model, we used available data on different components of the water balance for the Baksan River Basin, including MODIS snow cover data, ablation of representative glaciers, and isotope analysis data in the high-mountain Baksan River Basin of the Central Caucasus. Additionally, the results of the glacial runoff calculations using the A-Melt model were used to cross-validate both the ECOMAG and A-Melt models.

Data and methodology

Study area

The Baksan River Basin is located within the Central Caucasus region, which is part of the larger Caucasus Mountain range known for its high and inaccessible topography. This area represents one of the most important centers of glacial activity in the Caucasus. Approximately 232 km² of the region is covered by glaciers, and approximately 30% of this total area belongs to the Elbrus volcanic massif. Mount Elbrus, the highest peak in the Caucasus at 5,642 m a.s.l., has a glacier covering approximately 109 km2 (in 2017), which is approximately 10% of the total area glaciated in the whole of the Caucasus.

The overall glacier coverage within the Baksan Basin ranges from 7.4% at the Zayukovo outlet to 16.5% at the Tyrnyauz outlet (Fig. 1). The river’s hydrological regime is characterized by relatively low flows during the autumn and winter months, and significantly higher flows during the spring and summer months, complicated by sporadic peaks of rainfall floods. In this study, the gauging station “Zayukovo” was selected as the outlet point, which collected the majority of glacial runoff from the Baksan River (Fig.1). The average elevation of the river basin is 2,350 m with a total area of 2,100 km². Of this area, 64% is located in the high-altitude region of the basin above 2,000 m. According to data from hydrological gauges Tyrnyauz and Zayukovo between 1970 and 2022, mean long-term discharge values were 24.6 m³/s and 35.4 m³/s for these stations, respectively.

Runoff formation model

ECOMAG is a hydrological model designed to simulate key processes within the water cycle. It accounts for infiltration, evaporation, soil heat and water dynamics, snow accumulation and melt, freezing within the snowpack, as well as surface, subsurface, and groundwater flow - ultimately modeling river discharge. A comprehensive list of model parameters and equations is available in the study by Motovilov et al. (1999).

A simplified glacier block was used in the ECOMAG model for the study basin, based on the simple temperature index approach, assuming that ice and snow melting occurs during periods of positive temperatures (Hock 2003). The glacier block of the ECOMAG runoff formation model operates as follows: elementary watersheds containing glaciation are defined along with the percentage of glaciation within each watershed. The ice melt calculation subroutine is always activated in computational watersheds where glaciers are present. At the beginning of summer, when the watershed is covered with snow, only snow melts, using the snowmelt coefficient. As snow melts and recedes, the snow-covered area of the computational unit decreases. When the snow-covered area becomes smaller than the glaciated area, ice melting begins on the exposed surface using the ice melt coefficient. The melt coefficients may differ from basin to basin since they must implicitly take into consideration every element affecting the heat balance. Therefore, in the core version of the ECOMAG model, the possibility of calibrating the melt coefficients was retained while ensuring that the values remained within the natural range.

For this study, the basin was divided into 662 elementary watersheds, 102 of which were glaciated, 30 of which were more than 50% glaciated. For meteorological input, we used data from seven stations in the region’s meteorological network for the period 1977-2022 (Fig. 1). Information about the underlying surface, including topographic, soil, and landscape maps as well as glaciation data, was determined for every elementary watershed (Table 1).

Water discharge is a common hydrological modeling method, so it was used to calibrate and validate the model. We compared the observed and simulated values to assess the ability of the model to accurately reproduce river flow dynamics. Using this approach, we calibrated the model and evaluated its accuracy during validation to ensure realistic simulation of river behavior. We utilized criteria generally accepted in hydrology (Motovilov and Gelfan 2018), such as BIAS and Nash-Sutcliffe efficiency, to assess the quality of our modeling. The results are presented in section 3.1.

Data and methods of runoff formation model validation

We used remotely sensed snow cover data, isotope analysis data, observed ablation data from representative glaciers, and glacier melt modeling results to validate the ECOMAG model (Table 2). This section provides an overview of the data used in the multi-objective validation.

Snow cover

We used data on snow cover derived from Moderate Resolution Imaging Spectroradiometer (MODIS) images, which are instruments used by the Earth Observing System (EOS) satellites Terra and Aqua (Justice et al. 2002), to help check our model’s accuracy. The spatial resolution for interpreting snow cover information was 500 m, with a revisit time of one day per area. Snow cover masks derived from MODIS satellite data with daily temporal resolution were obtained using the MODISNOW tool developed by Gafurov et al. (2016). The method was developed to estimate the actual pixel coverage for cloud-covered areas by using spatial and temporal information about cloud-covered pixels and processing the portion of input data by removing a fraction of cloud-covered pixels. The ECOMAG model consists of a module that generates the outputs of individual characteristics that determine the runoff formation for each small catchment. This information was used to create maps showing the spatial distribution of snow cover on any selected day within the study area. In this study, we explored the potential of using additional data from MODIS Aqua and Terra images to validate the temporal and spatial distributions of snow cover. Daily data from 2000 to 2019 were used to validate the simulation of basin snow cover (%) using the ECOMAG model. Additionally, the spatial distribution of snow cover was validated using the example of the outstanding water year 2017.

Isotopic hydrograph separation

The results of isotopic analysis of natural waters can serve as an additional source for validating the runoff formation model. The distinct differences in the isotopic composition of glacial meltwater, snow water, and liquid precipitation allow for the estimation of the proportion of these components in total river runoff.

Dissecting a hydrograph into its genetic components using natural tracers is a widely used method for glacier-fed basins (Hoeg et al. 2000; Vasil’chuk et al. 2016). Tracers can be either total mineralization or specific ions, atoms, or molecules.

The results of the isotopic hydrograph separation performed for the Baksan River by Rets et al. (2024) were used in this study. Event-based isotope sampling was conducted along the Baksan River during the 2020-2021 period, approximately 5-6 times per season (Table 3, Fig. 2). A heavy oxygen isotope with an atomic mass of 18 (δ18O) was used as the tracer. After the sample analysis, the hydrograph of the Baksan River was divided into two components: meltwater runoff (seasonal snow and glacier meltwater) and rainfall runoff (Rets et al. 2024). The study does not outline subsurface/groundwater flow as a runoff source; instead, it considers it as a mixture of liquid precipitation, seasonal snow, and glacier ice melt navigating through the subsurface portion of the watershed. Since water stable isotopes are conservative tracers, their concentration remains unaffected by the flow pathway unless it is associated with substantial kinetic fractionation, such as evaporation from water reservoirs or watershed surfaces, or snow sublimation (Oshun et al. 2016). Genetic partitioning of a hydrograph using natural tracers is based on the water balance equation (Rets et al. 2024).

The ECOMAG model of runoff formation, based on the available set of model outputs, allows for the estimation of runoff components only for a month. To estimate the genetic components, the monthly average outputs of the water yield from snow and ice, precipitation on soil, and evaporation from soil for each elementary catchment were used. These data were compared with the meltwater and rainfall runoff components (m3/s, %) based on the results of the isotope analysis. The contribution of the glacial component was considered by increasing runoff when the model’s glacier block was activated. From all sampling dates, the mid-month dates were chosen for comparison with the model results.

Mass balance glacier observations

The mass balance time series from the World Glacier Monitoring Service (WGMS) dataset can serve as a valuable reference for validating ice-melt simulations. By comparing modeled glacier ablation with observed data, researchers can assess the accuracy of simulated ice-melt processes (Eis et al. 2021; Schaefli et al. 2005). The outcomes of the mass balance studies of the reference glaciers of the North Caucasus were employed to validate the runoff formation model. The continuous series of observations on the small valley glacier Djankuat (2750-3670 m a.s.l, 2.5 km²) started in 1967 and continues to date (Popovnin et al. 2024). The Garabashi glacier is a part of the Mt. Elbrus glacier volcanic massif and is located on the southern slope of the glacier system (3400-4900 m a.s.l., 3.9 km²). Monitoring began in 1982 and has continued to date (Rototaeva et al. 2019). Both glaciers are included in the observation system of the World Glacier Monitoring Service (WGMS) and are now considered representative of the Caucasus region. The mass balance observations included annual, winter, and summer measurements of mass balance according to ablation stakes, radar and manual snow measurements, and snow pits.

Glacier melting modeling

The additional method used for validation was glacial runoff calculated using the A-melt energy-balance model. One of the most important hydrological variables is glacial runoff, which is a component of the water–ice balance. More complex methods of calculating the glacier melting amount based on energy-balance methods are commonly developed for glaciological purposes (Gabbi et al. 2014; Reveillet et al. 2017; Sakai and Fujita 2017; Shannon et al. 2019). These models require detailed meteorological data as inputs, such as long and shortwave radiation balance, wind speed, humidity, and pressure. At the same time, it is complicated to transfer energy-balance calculations to a regional scale because of the large amount of input data required for such models. Simpler temperature index approaches are commonly used for regional runoff modeling in large basins. In this study, we used the results of runoff modeling from the Elbrus glaciers in the Baksan River basin using an energy-balance A-Melt model of snow and ice melt in the high mountain zone and the ECOMAG runoff formation model results for the cross-validation of both models.

The A-Melt model is based on the heat balance equation and calculates the mass balance of the glacier and its components based on meteorological data. The model simulates the vertical distribution of density, temperature, and water content in the snowpack, considering water seepage, refreezing in the snow thickness, and the snow–ice and firn thickness boundary and runoff. The A-Melt model has been successfully adopted for glaciers of the North Caucasus (Rets et al. 2011), Central Tien-Shan (Rets et al. 2021), and Svalbard Archipelago (Elagina et al. 2021). A more precise description of the A-Melt model has been provided in previous studies (Elagina et al. 2021; Rets et al. 2011; Rets et al. 2021). The A-Melt model experiment procedure and the results used in the current study for ECOMAG model validation are fully described step-by-step in (Elagina et al. 2024) and will be briefly described further.

Glacial runoff from the southern slope glaciers of Mt. Elbrus contributes to the formation of river runoff in the Baksan River Basin. After the model was applied for the entire Mt. Elbrus glacier system, the model output was extracted for the glaciers Bolshoy Azau, Maliy Azau, Garabashi, Terskol, Basin 25, Irik, and Irikchat, which supply the Baksan River basin with meltwater.

The model computes the heat balance equation for each point of the 250-m spatial grid. As meteorological input, we used data from the ERA5 reanalysis (Hersbach et al. 2020), which allows the reliable reconstruction of the meteorological conditions on Elbrus in accordance with the experiments given in the study (Mikhalenko et al. 2020). For the incoming radiation data, we used reanalysis data for a single level. Reanalysis data for air temperature, relative humidity, and wind speed were obtained for various pressure levels (500, 550, 600, 650, and 700 hPa), corresponding to an altitude range of 3000–5500 m above sea level. The values of the aforementioned meteorological parameters that fell between the input pressure levels were interpolated. The input data related to topographical features such as glacier area and outlines and surface elevation were obtained by analyzing available space data. The data sources include Spot 7, Landsat imagery, 2017 Pléiades DEM by the Pléiades Glacier Observatory initiative of the French Space Agency (CNES) and the 1997 0.5m-resolution DEM by the Faculty of Geography, Moscow State University (Zolotarev and Kharkovets 2000). To identify the surface type (ice/firn/snow), we used the Google Earth Engine automated algorithm proposed by Li et al. (2022), which allows classification of bare ice and snow cover on glaciers using Landsat images. The model does not simulate ice flow, so the topography features, surface type, and glacier boundaries were updated in the input data in 1984, 1996, 2006, 2011, 2016.

Because the model does not include the orographic component of precipitation distribution, a spatial coefficient distributed over the 250-m spatial grid was applied to distribute the solid precipitation amount over the entire area of Mt. Elbrus. The procedure is as follows:

1. We calculated the average summer mass balance and obtained a map of its average distribution for Mt. Elbrus from 1997 to 2017
2. The average winter mass balance map for 1997-2017 was determined by subtracting the geodetic mass balance map of the Elbrus glacier system for 1997-2017, obtained from Kutuzov et al. (2019) and the modeled average summer balance map for the same period. This subtraction provides an understanding of solid precipitation patterns across the glacier system.
3. The average winter mass balance at every grid point was then normalized by the amount of precipitation according to the ERA5 reanalysis data for the future calculation of winter mass balance for each balance year.

The water that formed during the ablation process, excluding the water that refreezes in the snow and firn layer and at the boundary with ice, that is, gravitationally free water ready to flow off the glacier, is considered glacial runoff in these calculations. The contribution of seasonal snow covering the surrounding slopes was excluded from the calculation of glacial runoff by Elbrus. Meltwater infiltration into the soil was not considered in the calculations. The values of glacial runoff from Elbrus were converted to units of volume per unit time. These data were then compared to the glacier runoff calculations obtained using the ECOMAG model. The glacier runoff in the ECOMAG model was calculated as glacier and snow ablation, excluding water refreezing in the snowpack.

Results

River discharge

The study (Kornilova et al. 2023) provides a step-by-step description of the model calibration and validation process. Briefly, data from the Tyrnyauz and Zayukovo gauging stations on the Baksan River were used, with a calibration period of 2000–2008 and a validation period of 2009–2017. Key runoff-influencing parameters were adjusted, including evaporation coefficients, active layer thickness, melting coefficient, snow water-holding capacity, and temperature gradients. The model demonstrated high accuracy (NSE > 0.8, annual runoff deviation < 6%), as illustrated in Fig. 3 and Table 4.

Snow cover

MODIS daily snow cover data from 2000 to 2019 were used to validate the ECOMAG model. The model effectively captured snow cover dynamics (Fig. 4a). The coefficient of determination R² for the average monthly snow cover values was 0.85, with a relative error of 20% (Fig. 4b). However, the model tended to overestimate snow cover values, particularly during the spring and winter (Fig. 4c).

To validate the runoff formation model for the spatial distribution of snow cover, the calculation results for the first day of each month in 2017 were used and compared with the snow cover data from the MODIS Terra and Aqua satellites (Fig. 4d). Overall, the model accurately simulated the annual course of snow cover changes. The peak was observed in February, followed by a gradual decrease until a minimum was reached in August and then an increase until December. However, during the spring months, the modeling results showed higher snow coverage than that observed in the MODIS data. Nevertheless, the modeling results for snow cover at the end of August agreed well with the MODIS data and generally corresponded to the mean annual position of the climatic snowline. It should be noted that the model data were presented for specific catchments with an average area of 9.9 km². Furthermore, within a specific catchment area, all model parameters were averaged. The pixel area of the satellite image was 0.25 km², which could be one reason for the ECOMAG snow cover model’s overestimation of snow cover model.

Isotopic hydrograph separation

From all sampling dates, mid-month dates were chosen for comparison of the isotopic analysis and model results (Fig. 5a-d). Based on the comparison, it can be seen that the model follows the intra-annual variability in the runoff sources, specifically a decrease in meltwater runoff and an increase in rain runoff during the summer months.

The data suggest that the proportion of snowmelt runoff decreased during the summer months, whereas the proportion of rain runoff increased. According to the model data, snowmelt accounts for approximately 70% of the total runoff in June but only approximately 35% in August. Additionally, the amount of snowmelt naturally decreases from the source to the downstream areas, while rain runoff increases. At the Elbrus station, the proportion of snowmelt ranged from 40% in August to 90% in June. It varied between 20% and 45% at Islamey station.

The average relative error between the model and the actual values for all measurements (four sampling points) was 19 %. The coefficient of determination for the meltwater component is 0.81 and for the rainfall component it is 0.86. The model tends to overestimate snowmelt values in June and underestimate them in July and August, on average.

Mass balance glacier observations

The actual values of ablation for reference glaciers for the period 2000-2009, obtained from the data of the World Glacier Monitoring Service (Zemp et al. 2015) and (Popovnin et al. 2024), were compared with the model values of snow and ice melting for the period of field observations (May-September) in the elementary catchment where the glacier is located (Fig. 6).

For the Garabashi Glacier, there was a good agreement between the modeled and actual values, with a discrepancy of -12%. However, for the Djankuat Glacier, the relative error was much higher, averaging -55% over the study period. Coefficient of determination R² is 0.20 for the Garabashi Glacier and 0.18 for the Djankuat Glacier, respectively.

Glacier melting modeling

The results of the runoff simulation for the Mt. Elbrus glaciers using the ECOMAG and A-Melt models are shown in Fig. 7. Fig. 7a presents the annual runoff values from the ECOMAG model extracted for the Elbrus Basin and the annual glacial runoff series according to the A-Melt model, which includes ice and snow meltwater. The series obtained using the A-Melt energy-balance model reflects the values of the average annual flow, considering meltwater, freezing in the snow and firn layers, and at the border of snow and ice. The model results were compared with the model output extracted for the glaciers Bolshoy Azau, Maliy Azau, Garabashi, Terskol, Basin 25, Irik, and Irikchat, which supply the Baksan River Basin with meltwater. Fig. 8b shows a comparison of the annual runoff values. The results of simulations indicate that the average runoff values are in good agreement (2.0 m³/sec by ECOMAG and 2.4 m³/sec by A-Melt, pBias = 18%). The average snow and ice glacial runoff components assessed by A-Melt model are 1.0 and 1.5 m³/sec, and 1.4 and 0.6 m³/sec assessed by ECOMAG model correspondingly. The higher values of the snow component of runoff and the lower values of the ice component of runoff according to the ECOMAG model compared to the results of modeling according to the A-Melt model are due to the difference in approaches to determine the type of surface for which the calculation is made (see Discussion).

Discussions

In mountainous areas, access to observational hydrometeorological data can be limited because of the difficult terrain, harsh weather conditions, and high costs of monitoring stations. Considering these challenges, it seems unrealistic to expect an increase in observational network density in the study regions in the near future. Therefore, relying solely on traditional methods for validating discharge may not provide a complete assessment of hydrological model performance.

Multi-objective validation offers a way to evaluate the model’s performance in multiple aspects of the water balance, enhancing the reliability of the assessment by incorporating various hydrological components and providing a more comprehensive understanding of the accuracy and ability to simulate complex mountain hydrology (Table 5).

In general, the ECOMAG snow cover model showed good agreement with MODIS data in terms of inter- and intra-annual variability. However, the runoff formation model indicated that snow cover during the spring months was higher than that observed from the remote sensing data. This finding suggests that meteorological data from stations located in the Baksan River Valley may not be sufficient to accurately represent the snow cover in adjacent side valleys with complex topographic conditions.

To address this issue, the use of high-resolution meteorological data, reanalysis, mesoscale modeling, and orographic precipitation modeling (Toropov et al. 2023) could improve the model’s ability to simulate snowmelt dynamics and overall hydrological processes in the Baksan River Basin.

MODIS data provide information only on the extent of snow cover but not on snow depth. In high-mountain regions, the accurate monitoring of snow depth is challenging because of the extremely limited number of snow measurement stations. Therefore, when validating the model based on snow cover, we primarily assessed the spatiotemporal dynamics of its distribution, whereas the contribution of snowmelt to total runoff can only be inferred indirectly through isotope analysis.

The similarity between the daily values of snowmelt and rainfall runoff, estimated through isotope partitioning, and the corresponding monthly average values derived from the model increases confidence in the accurate representation of runoff formation processes in the Baksan River Basin. Such evidence is essential for the physical justification of the long-term runoff change forecasts based on this model.

To accurately compare the absolute values of runoff components from the ECOMAG model simulation with isotope analysis results, continuous sampling at selected river cross-sections is necessary.

The comparison between simulated glacial runoff and ablation produced less than optimal results. The outcome suggests potential inaccuracies in the parameterization of the model or limitations in its ability to represent the complex glacial topography. The underestimation of the model’s ablation can be explained by the fact that the catchment where the glacier is located has a dissected relief with steep slopes, resulting in uneven distribution and melting of the snow cover. In the model, the entire catchment is given a single elevation corresponding to its geometric center, which does not accurately reflect the complex processes of snow and ice melting in such a terrain.

In reality, these processes are much more complex, leading to a systematic error in the model results, which are underestimated compared with the actual conditions. The large relative error in the ablation model for the Djankuat glacier may be due to the significant proportion of avalanche-fed ice, for which the runoff formation model did not account. In contrast, avalanche activity is absent on the Garabashi glacier, making it a more suitable location for testing these models. Detailed models of snow and ice melting in the nival-glacial zone for individual glaciers, such as those proposed by other authors (Fyffe et al. 2021; Kinnard et al. 2021; Rets et al. 2021; Van Pelt et al. 2019), provide the opportunity to take into account the topographic features and describe the melting process in greater detail.

According to the results of this study, the simulated glacial runoff of the Elbrus glaciers using the A-Melt model is in good agreement with the values obtained by the ECOMAG model. At the same time, the ECOMAG model showed higher values of the snow component of runoff and lower values of the ice component of runoff compared to the results of modeling according to the A-Melt model. This variation can be explained by the differences in approaches to glacial runoff calculations. The glacier block of the ECOMAG runoff formation model operates as follows: elementary catchments with glaciation are delineated, as well as the fractional extent of glaciation within each elementary catchment. At the onset of summer, when the catchment is snow-covered, only the snow melts. As snowmelt proceeded, the fractional snow cover of the computational element decreased. When the fractional snow cover became less than the fractional extent of glaciation, ice began to melt in the area exposed to snow. The ice melt calculation depended on the ice melt coefficient set in the model. The values for the ECOMAG model were obtained at the central point of the catchment. The snow and ice delineation in the A-Melt model is based on process modeling, and a surface type is determined for each point of the spatially distributed grid. The differences in the results are primarily due to the greater detail of the description of the processes in the A-melt model. In the ECOMAG model, the quantitative assessment of ice and snow melting and all related processes is described mainly by the snow and ice melting coefficients. Therefore, the resulting differences in the estimates of the ice and snow components of runoff are likely due to an underestimation of the ice melt coefficient in the ECOMAG model or an underestimation of the percentage of the elementary catchment area occupied by glaciers. This conclusion is supported by the fact that the ECOMAG model underestimates the melt rates of the Garabashi and Djankuat glaciers.

Refining glacier-related inputs, such as surface type, melt coefficients, and meteorological data, can enhance the performance of the model in simulating glacial runoff dynamics. Research has long focused on improving the glacial component of hydrological models (Addor et al. 2014; Etter et al. 2017; Finger et al. 2015; Horton et al. 2006). However, glaciers are not well represented in large-scale models. Most models use the temperature index method, which applies a melting coefficient to simplify complex glacial processes. To improve accuracy, it is essential to couple glacial and hydrological models that can help predict future discharge changes more reliably, especially in large glacier-fed rivers, where meltwater plays a crucial role. A better representation of glaciers is particularly important for assessing the impact of climate change, especially during summer and extreme years.

Conclusions

For the first time, the ECOMAG model was tested and validated using a range of resources, including MODIS snow cover data, isotope hydrograph separation, glacier mass balance observations, and glacial runoff simulations from the A-Melt model. This study is the first to compare the results from two independent models (ECOMAG and A-Melt) to assess glacial runoff and refine the distribution of snow and ice melt contributions. Additionally, a quantitative comparison of modeled runoff sources with isotope analysis data was performed for the first time to enhance the accuracy of seasonal runoff simulations in glacier-snow-fed basins.

The findings indicate that the ECOMAG simulations performed well in the sparsely gauged Baksan River basin. A good correlation was observed between the simulated and observed daily discharge values for both the calibration and validation periods at Tyrnyauz and Zayukovo stations. The percentage bias in annual runoff volumes was less than 6.5%, and the Nash-Sutcliffe efficiency (NSE) criterion for daily discharge exceeded 0.8 for both periods and locations.

To validate the model, data form the MODIS satellite spectroradiometer were used to measure snow cover in the catchment, confirming that the model accurately simulated seasonal variations in snow cover. The coefficient of determination for the monthly mean snow cover values was 0.85, with a relative error of approximately 20%.

The validation of the snowmelt runoff model based on isotope analysis showed that it accurately reflected the seasonal variation of the input sources. Furthermore, the proportion of meltwater runoff decreased naturally when moving downstream, whereas the proportion of rainfall runoff increased. The similarity between the daily meltwater and liquid precipitation runoff values obtained from isotope hydrograph separation and the monthly averages of the model adds confidence to the accuracy of the representation of snowmelt processes in the Baksan River Basin.

The validation of the model against mass balance measurements of reference glaciers in the central Caucasus (Djankuat and Garabashi) revealed that the processes of snow and ice melting and redistribution in such dissected terrain are more complex in reality, leading to systematic errors. The best agreement was achieved for the Garabashi glacier, with an average relative error of 12%.

Based on a comparison of the model outputs from ECOMAG and A-Melt, the simulated glacial runoff values were in good agreement. However, the results of the modeling indicated differences in the proportion of snow and ice components in the runoff, which can be attributed to differences in the snow and ice delineation methods and methods used to calculate the melting layer of ice and snow within the computational blocks of the models.

The successful calibration and validation processes further demonstrated the flexibility of the ECOMAG model in conducting comprehensive water resource assessments in complex mountainous regions.