Modeling of Extreme Hydrological Events in the Baksan River Basin, the Central Caucasus, Russia

Abstract

High mountain areas are prone to extreme hydrological events, and their study is especially important in the context of ongoing intensive deglaciation. In this research, a model “chain” consisting of a hydrodynamic model and a runoff formation model is adopted to simulate a glacier lake outburst flood (GLOF) from Bashkara Lake (the Central Caucasus, Russia) and its effect on downstream. In addition to an actual GLOF event that occurred on 1 September 2017 and led to casualties and significant destruction in the Adylsu and Baksan Rivers valleys, possible scenarios for the re-outburst of the lake are considered. The hydrographs of the outburst and the downstream movement of the flood wave along the Adylsu River valley are estimated using STREAM_2D two-dimensional hydrodynamic model. The water discharges in the entire river network of the Baksan River are assessed using the ECOMAG (ECOlogical Model for Applied Geophysics) runoff formation model. The output flood hydrograph from the hydrodynamic model is set as additional input into the Baksan River runoff formation model in the upper reaches of the Adylsu River. As a result of the simulations, estimates for the contribution of GLOFs and precipitation to an increase in peak discharge along the Baksan River were obtained. The actual outburst flood contributed 45% and precipitation 30% to the peak flow in the Baksan River at the mouth of the Adylsu River (10 km from the outburst site). In Tyrnyauz (40 km from the outburst site), the contributions of the outburst flood and precipitation were equal and, in Zayukovo (70 km from the outburst site), the outburst flood contributed only 20% to the peak flow, whereas precipitation contributed 44%. Similar calculations were made for future potential re-outburst flood, taking into account climatic changes with an increase in air temperatures of 2 °C, an increase in precipitation of 10% in winter and a decrease of 10% in summer. The maximum discharge of the re-outburst flood in the Adylsu River mouth, according to model estimations, will be approximately three times less than the discharge of the actual outburst on 1 September 2017 and can contribute up to 18% of the peak discharge in the Baksan River at the confluence.

1. Introduction

Continuing global warming changes the frequency and magnitude for many natural hazards [1,2]. An increasing rate of glacier downwasting is reported for most mountain areas [2,3]. The same situation is observed in the Central Caucasus, Russia [4,5,6,7]. There is also high confidence that the number and area of glacier lakes will continue to increase in most regions in the coming decades, and new lakes will develop closer to steep and potentially unstable mountain walls where lake outbursts can be more easily triggered by the impact of landslides [8,9,10,11,12,13]. In contrast to the number and size of glacier lakes, trends in the number of glacier-related floods are not well known for recent decades [14,15], although a number of phases of increased and decreased flood activity have been documented for individual glaciers in North America and Greenland [16,17]. Glacier lake outburst floods (GLOFs) occur in various mountain regions in the world, including the Himalayas [18,19,20,21,22], Tien Shan [23,24,25,26], the Andes [27,28,29], the Alps [30], North American Cordillera [31] and the Caucasus [32,33,34,35,36]. An increasing GLOF activity in the catchments with large glaciers is projected for the end of the 21st century [37]. The study and prediction of flash floods is an important hydrology and socio-economic task, especially in mountainous areas in the case of a small number of gauges [38,39,40].

Important triggers for a GLOF are an air temperature increase and extreme precipitation, which requires special attention given current climatic changes [23,41,42]. One of the problems studying outburst floods triggered by extreme precipitation is that rain-induced flow could amplify the outburst peak discharge and may cause a corresponding an increase in hazard and risk in downstream areas. If a GLOF is accompanied by a heavy rainstorm, it is difficult to trace the separate impacts of the GLOF and the rainstorm. For example, as result of the Chorabari Lake (Gandhi Sarovar, India) burst, which occurred during intensive rainfalls, more than 6000 people perished and several settlements, Kedarnath, Rambara, Gaurikund and Sonpreyag in the valley of the Mandakini River below the lake, were destroyed [10]. In the Himalayas, GLOF waves attenuate rapidly and only the farthest-reaching GLOFs may surpass the estimated 100-year floods for up to 85 (+45/−65) km downstream [43]. In the Caucasus, the reaches of the main rivers, which can be affected by the GLOF waves from the tributaries, is unknown. There is still a lack of analysis as to which part of the peak discharge can be attributed exactly to the outburst, and which to the rainfall, that evidently results in a water flow increase in the whole river network, and this lack of knowledge can lead to difficulties in designing protective measures. Thus, the task to assess the outburst flood contribution to the peak water discharge increase downstream of the burst lake during an extreme rainstorms is critically important.

Hydrodynamic models or “process chain” of models, considering only hydrodynamic modeling within the valley, are most commonly used for outburst flood studies [44]. This approach allows for the reconstruction of GLOFs, considering possible scenarios of outburst floods induced by different triggers [45,46]. However, this kind of modeling requires highly detailed data, and thus hydrodynamic modeling only allows for the estimation of the main hydrodynamic characteristics of an outburst flow for particular river reaches [47,48,49]. Runoff formation models are used worldwide for various high-mountain basins for regional assessments [50,51,52] and allow, for example, the assessment of possible climate impacts on runoff change [53]. This type of model can also be used for big river valleys to estimate the GLOF influence on the river runoff [54]. The chains with different types of hydrological models (runoff formation and hydrodynamic) are widely spread in hydrological forecasting [55] as they combine the advantages of both models and allow for the consideration of runoff formation processes, and describe flow movement on the floodplains and in river channels in great detail. In this study, combined numerical simulation were used to estimate GLOF impact, using both a hydrodynamic model to assess the outburst flood wave movement and a runoff formation model to take into account meteorological factors.

We selected for our study the most glaciated and populated part of the Central Caucasus, the Mt. Elbrus region (Russia), where high mountain resorts and corresponding infrastructure are located. One of the biggest recent GLOF events in the region was an outburst of Bashkara Lake and Lapa Lake on 1 September 2017 [56]. The lake outbursts with a total water volume of 0.8 mln. m³ involved in the flood formation were triggered by an extreme rainstorm accompanied by preceding rainy days, causing an overwetting of the dam and increasing of the lakes’ level. As a result, significant damage was observed both in the Adylsu River valley, in the right tributary of the Baksan River with the total length of 11.5 km, and in the Baksan River valley at a distance of more than 50 km downstream from the mouth of the Adylsu River. In particular, the gas pipeline and the federal highway leading to the main resorts in the Elbrus region were destroyed, the transport connection was interrupted for a two-week period, and three people perished. The total damage is estimated at about 810 mln. rubles (about 11.5 mln. EUR) [56]. The water volume of Bashkara Lake after the outburst decreased by more than ¾, from more than 1 million m³ to approximately 290 thousand m³ [56]. The intensive downwasting of the Bashkara Glacier took place during the last two years: two small lakes near Bashkara Lake appeared and increased considerably (Figure 1). Thus, the lake system is unstable and still potentially hazardous.

The lake system has been the object of annual monitoring by Lomonosov Moscow State University since 1999 [33], and extensive field data were collected and analyzed in 2017 during and immediately after the 1st September GLOF event [56]. Previous outburst flood modeling for the Adylsu River valley, performed before and after the GLOF on 1 September 2017 [56], concerned mostly the GLOF moving downstream from the Adylsu River; the impact of the GLOF on the Baksan River runoff has not been studied before.

The aim of this study is to set up a model system to simulate the hydrological routing under the scenarios of a different future climates and Bashkara Lake dam failure and to quantify the impact on the flood characteristics along the Adylsu and Baksan Rivers.

The following steps were performed in the research: (1) Modeling the actual outburst flood on 1 September 2017; (2) potential outburst flood modeling in case of rockfalls into the lake or in case of lake moraine dam destruction; (3) runoff modeling in the Baksan River basin for the historical period 2000–2017 and taking into account different future weather scenarios (according to IPCC (The Intergovernmental Panel on Climate Change) scenarios, an increase in background air temperatures of 2 °C in combination with an increase in precipitation of 10% in winter and a decrease of 10% in summer); and (4) assessment of the outburst flood impact on water discharges downstream from the Baksan River under various scenarios.

2. Background

2.1. The Baksan River Basin

The Baksan River basin is located within the Russian Central Caucasus, the highest part of the Greater Caucasus. This region is one of the major centers of glaciation in the Caucasus. Glaciers cover 232 km² [57], and approximately 60% are located at the southeast slopes of Mt. Elbrus.

The Baksan River belongs to the Terek River basin and is its second-order tributary. The Baksan River source is located in the southern slopes of Elbrus; the river starts from the Bolshoy and Maly Azau Glaciers at an altitude of 2680 m, then it crosses several high mountain ridges and flows into the Malka River near the town of Prokhladnyj at the foothills. The Baksan River basin is significantly glaciated; 9.3% of the modeled catchment area is covered by glaciers. The rest of the basin is covered by alpine and subalpine meadows and mixed pine-birch forests.

In the upper reaches of the Baksan River basin there are numerous mountain resorts with a total capacity of ca. 5000 people. The population alongside the Baksan River, including the towns of Tyrnyauz and Baksan, is more than 100,000 people.

The flow regime of the Baksan River is characterized by an October–April low-water period and a May–September high flow period, with sharply changing flow rates and a general maximum in July–August. The main hydrologic gauging station on the Baksan River is the Zayukovo station, with a catchment area of 2100 km² and a mean catchment elevation of 2350 m. In addition, the Tyrnyauz gauging station with a basin area of 838 km² is located on the Baksan River (Figure 2). The distance of the Tyrnyauz and Zayukovo stations from the Baksan River mouth is 125 and 82 km respectively. According to data analysed for the period 1977–2016 for the Tyrnyauz and Zayukovo gauge stations, the mean annual water discharge s are correspondingly 23.9 and 34.9 m³/s, with maximum daily water discharges of 1% reoccurrence ‒218 and 306 m³/s) and 10% reoccurrence ‒157 and 172 m³/s. A negative trend is characteristic for the annual peak flow of the Baksan River, which is quite common for the Central North Caucasus. The mean value of annual peak discharge is 6% lower at Zayukovo gauging station in 1981–2016 compared to 1931–1980 [58]. A slight increase of less than 10% in 1978–2010 compared to 1933–1977 is observed in the mean annual runoff in the Baksan River basin [59].

The Baksan River basin is a debris flow hazard area, especially taking into account intensive glacier melting and an increase in extreme precipitation [60,61]. Over the 68-year observation period (1951–2018) in the basin, 45 years were noted with debris flows, during which 440 debris flow event cases took place [62].

At the end of the 21st century, rising temperatures in the North Caucasus of 2–3 °C are predicted under the RCP4.5 scenario, and of more than 3 °C for the scenarios RCP6.0 and RCP 8.5. Precipitation changes are seasonal in all scenarios, for example, precipitation will increase by 10% in winter and decrease by 10% in summer according to scenario RCP4.5 [63]. Such an increase in summer air temperature is not compensated by an increase in winter precipitation and can lead to continuing degradation of glaciers, formation, expansion, and bursts of glacier lakes and change of GLOF hazards.

2.2. Lakes Near the Bashkara Glacier

Bashkara Lake is located in the upper reaches of the Baksan River, in the headwaters of the Adylsu River between the right lateral moraine and the valley Bashkara Glacier. The lake has appeared in the early 1930s; for a long time it was surrounded by ice shores on three sides. In August 1958 and 1959, as well as in October 1960, there were three lake outbursts through the grottoes in the ice wall, causing high-magnitude debris floods downstream from the Adylsu valley [64].

In addition to Bashkara Lake, there is Lapa Lake at the glacier terminus, which appeared in the 1990s. During the last three decades there was noticeable surface lowering at the Bashkara Glacier snout adjacent to the lakes. The area and volume of Lapa Lake significantly increased in 1999–2017, and Bashkara Lake’s size was quasi-stable [33,56]. In 2008, overtopping of the Bashkara Lake dam was registered for the first time during the period of observations (from 1990th) and an outburst hazard was declared. During the period 2009–2014, the water level in Baskara Lake was lower, and overflow over the moraine dam was not observed. However, in 2015–2017, the water level in the lake rose sharply; overflow was registered annually and its duration increased from one week in 2015 to two months in 2017. During field and aerovisual inspections, seepage through the moraine dam was also observed [56].

During the night from 31 August to 1 September 2017, there was an overflow and the consequent outburst of Bashkara Lake occured. The immediate trigger of the lake outburst was an extremely heavy rainstorm accompanied by strong overwetting of the moraine dam. A four-day period of rainy weather before the 31st of August ended with an extreme rain storm. The total precipitation during the night from 31 August to 1 September close to the lake Djankuat station was approximately 100 mm/8 h (Figure 3).

The pre-outburst level of Bashkara Lake was about 2672 m. According to data from the water level logger, the level in Lake Bashkara rose 55 cm in 5 h during this rainstorm. As a result, the crest of the lake dam was eroded during outflow. According to the data from the field survey in the days after the event, the water level in Bashkara Lake lowered by 17 m relative to the pre-outburst one. The GLOF wave moved first from the upper Bashkara Lake across the surface of the Bashkara Glacier snout, then partially drained Lapa Lake, and afterwards began to move downstream from the Adylsu River valley. Before the outburst, according to the volume curves, the volume of Lake Bashkara was 1,024,000 m³ and the volume of Lake Lapa was 274,000 m³. After the outburst, the water volume in the lakes had decreased to 287,000 m³ in Bashkara and 222,000 m³ in Lapa Lake (Figure 4).

Bashkara Lake dam is located below the local permafrost level, which is above 2700 m in this part of the Central Caucasus. However, buried ice is still preserved in Little Ice Age lateral moraines. Small lenses of ice were found in the dam during field inspections before the event. Clear GPR (ground penetration radar) reflections at the dam section were not received during surveys in 2009 [65], hence it can be concluded that there was no significant ice content at the dam. There were no observation of dam crest lowering in 2010–2017, therefore we conclude that buried ice melting did not play a significant role in GLOF initiation.

3. Data and Methods

3.1. Models Chain Composition

We used a composition of two models (”model chain”) to simulate the Bashkara Lake outburst flood and its impact on the Baksan River hydrograph. Two-dimensional hydrodynamic model STREAM_2D was applied to simulate the outburst process and the flood wave movement downstream from the Adylsu River (river length 11.5 km, basin area 101 km²).

To take into account the runoff process formation in all the high mountain parts of the Baksan River basin, from the upper reaches to the Zaykovo gauge (Baksan River length—87 km, basin area—2100 km²), we used the ECOMAG runoff formation model. The outburst flood discharges from Bashkara Lake, simulated on the basis of the hydrodynamic model, were set as additional contributions into the runoff formation model in the upper reaches of the Adylsu River downstream from Lapa Lake. Separate datasets were prepared for each type of model. Data and models are described in the next sections, following the workflow chart (Figure 5).

3.2. Data for Hydrodynamic Modelling

Digital terrain models (DTM) and the results of bathymetric surveys for the Bashkara and Lapa lakes were used as initial data input to run the STREAM_2D model. The bathymetry of the lakes before the outburst was obtained during surveys in 2008–2014 using a boat-mounted Lowrance 525 CF echo sounder; and corresponding area-volume curves were constructed [56]. A bathymetric survey was performed in September 2020 to estimate the current volume of the lakes.

We used the DTM generated from the SPOT 7 stereo pair (resolution—3.2 m, survey date 01.08.2017) to simulate the actual outburst flood. The morphometry of the dam break section was set according to the field investigation as follows: the width of the breach 15–20 m and the depth 30–40 m. The DTM generated from the Pleiades stereo pair (resolution—1 m, survey date 03.09.2017) were used for scenario modeling of the recurrent flood. Water channel relief before and after the outburst was reconstructed using a morphometry of cross-sections, obtained during the field surveys, and were incorporated into the corresponding DTMs.

Roughness coefficients in the Adylsu valley were chosen according to different types of surfaces; these types were defined using automated classification of satellite image on the base of ERDAS Imagine software by the method of k-means [66]. Values of the roughness coefficients were initially set using published data [67] as 0.035 for river channels and 0.05 for floodplains, and then adjusted during model calibration. To set ice channel roughness, we used our data from a special field survey at the Dzhikiyugankez Glacier (Mt. Elbrus), which included measuring of flow velocities, water surfaces slopes, and channels morphometry. Roughness coefficients of 0.024 for supraglacial channels at debris covered glaciers were calculated on the base of these characteristics using the Chezy formula [67].

For the analysis of water level fluctuations in Bashkara Lake, our data from waterlogger Keller for 2014–2018 years were used. The initial lake water level for actual outburst modeling was set corresponding to the pre-outburst level. For the simulation period, rain intensity of 20 mm/h, which was set up homogeneously during the period of hydrodynamic simulation, was taken into account.

3.3. Hydrodynamic Modeling of the Outburst Floods from the Lakes Near the Bashkara Glacier

We applied the STREAM_2D software [68] based on the numerical solution of the non-stationary two-dimensional Saint-Venant equations in the shallow water approximation for hydrodynamic modeling [69]. The shallow water equations consider the main forces in a stream with a free surface (gravity, friction, pressure, and inertia; Coriolis’s force and wind influence can be considered in addition), and the three-dimensional orography of the channels and floodplain. The main equations and the descriptions of model computational schemes can be found in [69]. STREAM_2D is based on hybrid triangular-quadrangular meshes with a variable grid size and has a sediment transport and ice block [70,71]. STREAM_2D software is widely used in Russia for river valleys of different scale [70,71]. Previous studies of potential GLOF from Bashkara Lake were performed using the same model [33,72], and we chose this software for better comparison of results. Input data for hydrodynamic modeling are a digital terrain models and initial water levels; output data includes water levels, flow depths, velocities for each model cell, inundation zone, and water discharges in chosen cross sections.

In case there is a sufficient amount of initial information, including the bathymetry of the lake and morphometry of the dam break, the STREAM_2D software package allows us to calculate the hydrograph of the outflow directly from the hydrodynamic model. This approach was used to assess the outburst hydrograph from Tangjiashan Lake [73] and showed good concurrence between the modeled and observed discharges and water levels. Due to the availability of the necessary information, a similar approach was applied in our study.

For a detailed study of the GLOF, the hydrodynamic modeling area was divided in two parts, and two computational grids were built: separate grids for the lake system and the Bashkara Glacier area (upper reach model—hydrodynamic model No. 1) and for the Adylsu River valley (lower reach model—hydrodynamic model No. 2) (Figure 6). This division was done to simplify the simulations and calibration process. Using hydrodynamic model No. 1, we set a detailed grid onto the dam break section and simulated numerous outburst scenarios, model No. 2 was used to only estimate floods for the chosen scenarios with maximum outburst discharges. Taking into account the resolution of input relief data for the river channels and the dam section, a quadrangular grid with a step of approximately 2 * 5 m was set, and an irregular triangular grid with variable size 5–10 m was created for the floodplain.

On the basis of hydrodynamic model No. 1, the outflow hydrographs were simulated at the section of the dam breach (Figure 6, Section 1), as well as in the cross-section of the Bashkara River below Lapa Lake (Section 2). The hydrograph obtained at the outlet of Lapa Lake was used as an input hydrograph for the lower model for the Adylsu River valley. In addition to the input hydrograph, a water discharge of 15 m³/s was also set for the tributaries of the Adylsu River. This discharge was estimated on the basis of the observed water discharge in the Djankuat River during the night of 31 August to 1 September 2017, one of the tributaries where daily hydrological observations are carried out within the framework of the Djankuat glaciological station program [74]. The input discharge of the Baksan River at the confluence with the Adylsu River was set based on the results of the ECOMAG runoff formation model and assessed as 40 m³/s. On the basis of hydrodynamic model No. 2, flow velocities, flooding depths, and inundation zones in the Adylsu valley, as well as water discharges at the key sections at the Djantugan mountain hotel (Figure 6, Section 3), the Shkhelda mountain hotel (Section 4), and the Adylsu River mouth (Section 5) were estimated. Information about the timing of flood movement was used to calibrate the model No. 2.

The roughness coefficients were adjusted assuming that the modeled time for GLOF wave travel would fit the values obtained from interviews with local residents. According to witnesses, the head wave of the flood reached the Shkhelda mountain hotel (6.4 km below the source of Bashkara Lake) at approximately 1:20 a.m. September 1st 2017, the Elbrus settlement (10 km below the source)—no later than 1:30 a.m. Several options for specifying the roughness coefficients showed that when the channel roughness coefficient is changed from 0.02 to 0.05, the flood peak travel time changes by 5 min. For further simulations, a set of parameters No. 3 was used, where the travel time matched well with the witnessed information (

Table 1. Hydrodynamic model STREAM_2D calibration results based on local’s poll about flood timing.

		Local Time, 1/09/2017		Roughness Coefficient
		Shkhelda Camp	Elbrus Settl.	Channels	Floodplain	Forested Floodplain
Observed (locals poll)		1:20	1:30
Model parameters set	1	1:16	1:25	0.020	0.030	0.035
	2	1:21	1:32	0.050	0.060	0.070
	3	1:20	1:30	0.045	0.050	0.055

).

For verification of the upper and the lower hydrodynamic model, the data on flooding zones within the snout of the Bashkara Glacier and the Adylsu River valley, obtained from satellite images (Pleiades, resolution—1 m, survey date 03.09.2017) and from the results of modeling, were used. A comparison of the simulated and actual flood zones generally showed satisfactory visual agreement (Figure 7); the difference in the inundated area between the satellite image and modeling results was from 5% up to 20% depending on the section of the valley. However, the difference in the boundaries of the flooded zones along the banks can reach 30 m.

3.4. Data for Runoff Formation Modeling

The input data for the runoff formation model for the upper part of the Baksan River basin to the Zayukovo gauge included digitized landscape and soil maps with a scale of 1: 750,000, namely the void-filled shuttle radar topography mission (SRTM) DEM from 2000 with a spatial resolution 30 m., which is available at the website of the Consultative Group for International Agriculture Research Consortium for Spatial Information (CGIAR-CSI) via http://srtm.csi.cgiar.org/. We used the glacier outlines from the Randolph Glacier Inventory (RGI-6) to account for the parts of the basin covered with glaciers. RGI-6 is a data set of glacier outlines that combines previously existing outlines from different sources (including the Global Land Ice Measurements from Space (GLIMS) initiative data) with new data from various contributors [57].

As meteorological forcing for the runoff formation model, we used the information from seven weather stations: data from three stations (Terskol, Cheget, Nalchik) were available for the period 2000–2019 and earlier, and data from the full set of seven meteorological stations were available from 2008 to 2019 years. Daily runoff data from the Tyrnyauz and Zayukovo hydrological gauge stations for the period 2000–2017 years were used for calibration and verification of the runoff formation model.

3.5. Runoff Formation Model for the Baksan River Basin

In order to trace how the Bashkara GLOF affects the runoff of the Baksan River, ECOMAG (ECOlogical Model for Applied Geophysics) software, developed under the leadership of Yu.G. Motovilov [75,76,77,78], was used. The ECOMAG model has shown its high efficiency during studies in different regions worldwide, from large river basins (Lena, Makenzie, Yenisei, Northern Dvina, etc.) to small basins with different nature conditions, including present-day glaciers and permafrost [79,80].

ECOMAG is a distributed model where the basin surface is divided into elementary watersheds, taking into account topography, soil types, and landscapes. The model describes the main processes of the hydrological cycle: infiltration, evaporation, thermal and water regime of soils, the formation of snow cover and snow melting, the formation of surface, subsurface, and ground and river runoff. One dimensional kinematic wave approximation is used in the ECOMAG model for the simulation of surface and river runoff movement. The input meteorological data for the model are air temperature, precipitation, and air humidity deficit with daily time steps. The output data of the model are the water discharges in the river network and all components of water balance of each elementary watershed.

Parameters of the model are physically meaningful and can be related to the measurable characteristics of river basins. Most parts of the soil and landscape parameters are initially included in the databases of the ECOMAG model. During calibration, multipliers of each parameter, uniform for all types of soils and landscapes, can be changed. The list of parameters, as well as the model equations, can be found in [75,76,77,78].

The simplified glacier block in the ECOMAG model, which is based on the assumption that thawing from glaciated parts of the basin is unlimited during periods of positive temperatures, was applied for the study basin. From 662 elementary watersheds included in the runoff formation model schematization, 102 watersheds were glaciated, and 30 of these by more than 50%.

For calibration and verification, data from two hydrological gauging stations located on the Baksan River were used: the Tyrnyauz city and Zayukovo settlement. The calibration period was nine years (from 2000 to 2008), and the verification period was nine years from 2009 to 2017 (

Table 2. ECOMAG model calibration and validation periods performance criteria.

Period, Years	NSE, Daily		PBIAS, %
	Tyrnyauz	Zayukovo	Tyrnyauz	Zayukovo
2000–2008 (calibration)	0.82	0.85	−1.71	+3.24
2009–2017 (validation)	0.90	0.90	+5.98	+3.58

).

Parameters that have a significant impact on the runoff in the high-mountain Baksan River basin were identified and adjusted in the process of calibration. Among the more sensitive parameters for the study basin were multipliers for the evaporation coefficient, thickness of the active layer, melting coefficient, water-holding capacity of snow, critical temperature for the phase of atmospheric precipitation, critical temperature of snow melting, correction coefficient for density of freshly fallen snow, temperature, and precipitation gradients (Appendix A).

Good correspondence of simulated and observed daily discharges were obtained for both the calibration and verification periods for the Tyrnyauz and Zayukovo gauges. Percent bias of the annual runoff volumes did not exceed 6.5%; Nash-Sutcliffe criteria (NSE) [81] for daily discharges was higher than 0.80 for both periods and stations.

To take into account the GLOF flood wave, an additional input section located below Lapa Lake was introduced into the ECOMAG runoff formation model, in which the average daily water discharge for the period of the outburst flood obtained from the data of the hydrodynamic model was set. As a result, the runoff formation model reproduced peak flow well. For the outburst period in 2017 the percent bias of the peak water discharge estimation was 1% according to the data of the Tyrnyauz hydrological gauge (Figure 8).

3.6. Simulation Scenarios

3.6.1. Scenarios of the Bahkara Lake Outburst

The initial scenario for all types of modeling was the GLOF that had already occurred in 2017 (scenario hd0). Further, possible future scenarios of the lake re-outburst were analyzed. Two different options were taken into account as possible mechanisms for a next outburst.

The first one consists of a large rock fall into the lakes, forming displacement waves (scenario hd1). Outburst floods caused by similar processes are observed in various regions of the world [18,82,83]. This scenario was considered due to observed large rockfalls from the Bashkara Peak towards the Baskara Glacier that had taken place in the summer 2018 and in spring 2019, when rock avalanches had almost reached Bashkara Lake. In future, the possibility of repeated collapses is not excluded, considering the instability of the rocks on the mountain slopes around the lakes [84].

The second scenario assumes a re-outburst caused by water level increase (due to heavy precipitation or intense melting), accompanied by the expansion of the existing dam break. There was an increase in the lake level by 1.5 m in 2020 compared to 2017 (up to the mark 2673.5 m). By our estimations, in the case of a lake re-outburst, the dam break can deepen up to 8 m as a result of erosion. These parameters were assumed for the modeling of the second scenario (scenario hd2).

3.6.2. Scenarios of Different Meteorological Conditions of Runoff Formation in the Baksan River Basin during the Outburst Flood

To assess the impact of weather conditions on the Baksan River runoff, several scenarios were compared. The first one was a scenario with the actual meteorological conditions of 2016–2017, for which in autumn, on 1 September 2017, an outburst flood was superimposed (scenario rf1a+hd0). Another scenario that coupled GLOF with the maximum dischargers for during summer 2017 was also considered (scenario rf1s+hd0). Similar scenarios have been done for a potential future re-outburst (rf1a+hd2, rf1s+hd2) (

Table 3. Scenario descriptions.

№	Scenario	Description
hd0	Occurred outburst, hydrodynamic modelling only	Modeling the actual outburst flood on 1 September 2017
hd1	Possible next outburst in case of rock avalanche, hydrodynamic modelling only	Modeling the potential re-outburst in case of rock avalanche
hd2	Possible next outburst, hydrodynamic modelling only	Modeling the potential re-outburst in case of dam failure
rf1a+hd0	Occurred outburst, runoff formation modeling based on observed meteorology, autumn	Modeling the actual outburst flood on 1 September 2017 taking into account observed autumn meteorological conditions
rf1s+hd0	Occurred outburst, runoff formation modeling based on observed meteorology, summer	Modeling the actual outburst flood taking into account observed summer meteorological conditions
rf1a+hd2	Possible next outburst, runoff formation modeling based on observed meteorology, autumn	Modeling the potential re-outburst taking into account observed autumn meteorological conditions
rf1s+hd2	Possible next outburst, runoff formation modeling based on observed meteorology, summer	Modeling the potential re-outburst taking into account observed summer meteorological conditions
rf2a+hd0	Occurred outburst, runoff formation modeling based onprojected meteorology, autumn	Modeling the actual outburst flood on 1 September 2017 taking into account projected autumn meteorological conditions
rf2s+hd0	Occurred outburst, runoff formation modeling based onprojected meteorology, summer	Modeling the actual outburst flood on 1 September 2017 taking into account projected summer meteorological conditions
rf2a+hd2	Possible next outburst, runoff formation modeling based on projected meteorology, autumn	Modeling the potential re-outburst taking into account projected autumn meteorological conditions
rf2s+hd2	Possible next outburst, runoff formation modeling based on projected meteorology, summer	Modeling the potential re-outburst taking into account projected summer meteorological conditions

).

In addition, a scenarios with an air temperature increase of 2 °C in combination with an increase in precipitation of 10% in winter and a decrease of 10% in summer were modelled (scenarios rf2+hd2). These scenarios correspond to a possible change in meteorological conditions according to moderate IPCC scenario RCP4.5. To simulate the change in meteorological forcing, the so-called “delta-change method” [79] was applied. Observed daily air temperatures were increased by 2 °C for this scenario. For precipitation, separate change coefficients were applied for cold (October–March) and warm (April–September) seasons. As a result, a transformed series of daily values of precipitation and air temperature was obtained for each of the meteorological stations. To calculate the transformed values of the air humidity deficit, it was assumed that the dew point temperature changes by the same amount as the air temperature (i.e., by 2 degrees).

4. Results

4.1. Hydrodynamic Modeling of Glacial Lake Outburst Flood (GLOF) on 1 September 2017

According to the simulations for the scenario hd0, the outburst flood wave was flattened while moving downstream of the Adylsu River valley; the value of maximum water discharge decreased, the triangular shape of the hydrograph flattens, the time between the beginning of the increase in flow rates and the peak of the hydrograph in the section increased, and also the duration of the flood also increases. The modelled maximum flow discharge at the outflow of the Bashkara Lake was 710 m³/s and at the outlet of Lake Lapa—440 m³/s three minutes after the initial outburst. At the mouth of the Adylsu River, the peak of the flood was observed 40 min after the outburst, and the maximum water discharge was 320 m³/s.

The maximum flow velocities were observed in the dam breach section (16–18 m/s), at the inflow to Lake Lapa (12 m/s), at the Shkhelda mountain hotel (11–15 m/s), in the upper reaches of the Adylsu River (10–13 m/s) and below the confluence with the Djankuat River (12 m/s). It should be noted that the distribution of the maximum velocities coincides with the zones of intense vertical deformations, which were mapped during the field survey.

4.2. Scenario Modeling of Potential GLOF

The scenario of a rock avalanche impact on the lake (scenario hd1) did not lead to a significant response of the model for different variants of falling into the lake and volumes of an avalanche. This means, no displacement waves was modelled and the modeled water levels at the different control points remained close to initial values.

Based on the simulation for the second scenario (scenario hd2), which implies an increase in the existing dam break, hydrographs of flood wave movement in various sections of the Adylsu valley were obtained (Figure 9b). According to the modelling results for this scenario, the flood overcomes the system of lakes after six minutes. A subsequent flood wave flattening can be observed, and the flood wave movement is slower than the one of the 1 September 2017 event. The recession lasts for about 1 h. The maximum water discharge from Bashkara Lake for scenario hd2 is 298 m³/s, at the outlet from Lake Lapa 138 m³/s and at the mouth of the Adylsu River due to wave spreading, the maximum discharge is 101 m³/s (which is approximately 3 times less than for the actual outburst on 1 September 2017).

The maximum flow velocities are observed at the snout of the Bashkara Glacier (>15 m/s), in the dam break (10 m/s), and near the Shkhelda mountain hotel (7–8 m/s). Everywhere, in the upper reaches of the Adylsu River, the flow velocity is about 6 m/s, except for the wider part of the valley near the Djantugan mountain hotel. The maximum flow depths observed near the Shkhelda and the Djantugan mountain hotels are approximately 4–5 m.

4.3. Projected GLOF Impact on the Baksan River Peak Flow According to Scenario Modeling

Several scenarios were considered while modeling the runoff formation in the Baksan River basin: the superposition of the actual outburst flood on the autumn and summer runoff (scenarios rf1a+hd0, rf1s+hd0), the overlap of a possible repeated outburst flood on the contemporary autumn and summer runoff (scenarios rf1a+hd2, rf1s+hd2), as well as modified runoff regime according to climatic scenarios (scenarios rf2a+hd2, rf2s+hd2).

Based on scenario simulations, the impact of the outburst flood on the Baksan River discharge significantly decreases downstream. The modeling results for the 1 September 2017 outburst (scenario rf1a+hd0), in the Elbrus settlement on the Baksan River, at the mouth of the Adylsu River show that the outburst flood contribution into the total Baksan River water discharge accounts for 45% of the peak discharge, in the Tyrnyauz city (40 km downstream) it accounted for 30% and in the Zayukovo village (70 km downstream) it accounted for 20% (Figure 10a).

If the outburst of the lakes in 2017 overlap with summer annual peak flow (scenario rf1s+hd0), then both the base and total water discharges would be higher, while the GLOF contribution to the total water discharge at the indicated sections would decrease by 38%, 26%, and 17%, respectively (Figure 10b and Figure 11).

An assessment of the contribution for a potential GLOF to the discharge in the Baksan River showed that its influence seems not to be so significant as for the occurred GLOF. The contribution of the outburst flood varies from 18% near the Elbrus settlement to 5% near Zayukovo under the “autumn” scenario rf1a+hd2 (Figure 10c) and from 13% to 5% under the “summer” scenario rf1s+hd2, depending on the distance from the outburst site (Figure 10d).

Contributions of the GLOF, base runoff, and precipitation to peak discharge in the Adylsu and Baksan Rivers is shown in Figure 10. The GLOF accounts for almost all peak discharge downstream from the lakes, drops rapidly after the confluence with Baksan and becomes equal to the contribution of the precipitation input at a distance of 40 km from the lakes (Figure 11).

Under scenarios with changed meteorology, a rise in runoff in winter and a drop in summer is projected. The contribution of an outburst flood to the total peak discharge decreases (Figure 12). For example, for September 1st, 2017 GLOF, under changed meteorological conditions, the contribution of the outburst near the Elbrus settlement would be 42%, and would decrease to 30% near Tyrnauz in autumn (scenario (scenario rf2a+hd0), and correspondingly from 41% to 30% in summer (scenario rf2s+hd0). The contribution of the repeated outburst flood from the total discharge decreases in the village of Elbrus to 16%, in Tyrnyauz at 40 km to 10%, and in Zayukovo at 70 km to 6% in the “autumn” scenario rf2a+hd2 (Figure 12a) and “summer” scenario rf2s+hd2 respectively (Figure 12b).

In all cases, the contribution of the GLOF to the peak discharges of the Baksan River is less in summer compared to autumn conditions. In the case of projected changes in meteorological conditions, the role of the base runoff in total peak flow of the Baksan River increases under the “autumn” scenario, and decreases under the “summer” scenario.

5. Discussion

In the presented study, an actual GLOF occurred in the Baksan River basin on 1 September 2017, possible scenarios of a potential re-outburst in the case of rock avalanche impact and expansion of the existing slit, as well as calculations of the outburst flood propagation downstream from the Adylsu valley were implemented using only the STREAM_2D software package. A number of hydrodynamic models, or “process chains” are used frequently for modeling outburst floods, as they allow for a complete description of the process of the outburst, from the beginning to the end, including assessing the potential hazard [85,86]. The results of each process are used to simulate subsequent processes in the chain. Various software and models are often used at each stage. For example, in [87], potential avalanches were modeled using RAMMS [88], displacement wave dynamics using FLOW3D [89], moraine erosion using BASEMENT [90], and flood wave propagation downstream and inundation using FLO-2D [91]. Differences between two-dimensional hydrodynamic models, such as Delft 3D [92], HEC-RAS [93,94], TELEMAC [95], and BASEMENT [96], mainly consist in the methods of schematization of the computational domain and the composition of additional blocks of models. The STREAM_2D software, based on hybrid triangular-quadrangular grids with a variable grid size, gives the opportunity to describe flow movement in great detail, but at present it does not have a debris flow module developed, as, for example, software packages FLO-2D [97].

Previous scenario modeling of the outburst flood in the Adylsu River valley was performed before the GLOF on 1 September 2017 [33]. The modeling was done on the basis of the “River” model, the predecessor of the modern version of STREAM_2D [98] and FLO-2D [99] software. After the GLOF, the comparative modeling was also carried out with simulations based on RAMMS software [72]. A comparison of the different models with and without consideration of debris transport has shown close values of the travel time, flow velocities, depths, and the zone of maximum inundation. Therefore, in our work, only simulations of the water flow movement were carried out on the basis of STREAM_2D without considering debris entrainment and transport.

Earlier studies carried out before the actual Bashkara Lake outburst were based on the assumption that the lake will break through an englacial tunnel in the Bashkara Glacier; such an outburst mechanism implies a lower maximum water discharge, a longer outflow duration, and an asymmetry of outburst flood hydrograph [24,33,100]. To calculate a possible hydrograph of the lake outburst, the model developed by Yu.B. Vinogradov [101] was used. The maximum outburst discharge was estimated as 123.5 m³/s and the total volume of the outburst flood as 764 * 10³ m³ [100]. In reality, from leaning against the Bashkara Glacier in the early 21st century, a moraine dam became less stable due to the significant lowering of the snout surface in 2008–2020. Causing the entrance to the englacial channel to be moved approximately 30 m from the lake bank. The moraine dam has been significantly destroyed during the GLOF and the narrow entrance into the englacial channel could not bypass the flood wave. The discharge, much higher than expected, was a result of overflow through a moraine dam instead of a burst through a tunnel. A comparison of our modeling results with estimates, which were obtained before the 2017 outburst, using the same modeling software, has shown that the flood reached the Adylsu River mouth less than in 1 h in both cases. However, flood depths and velocities for the September, 1, 2017 GLOF were 1.2–1.5 times higher than simulated earlier. It might be explained by another form of outburst hydrograph.

The characteristics of the GLOF modeled in our work, on the basis of STREAM_2D, are in general consistent with both the estimates obtained during the field surveys after the 2017 GLOF, and the results based on the RAMMS simulation performed after the outburst [72]. The asymmetric shape of the outflow hydrograph from the lake obtained as a result of hydrodynamic modeling corresponds to the shape of the outburst floods hydrograph, with an abrupt and rapid destruction of the moraine dam, and is used in other works [44,102].

The modeling was carried out without considering lateral erosion, thus the difference between the boundaries of the flooded zones along the banks reached 30 m.

According to the assessments from the field surveys in early September 2017, the maximum discharge of the frontal flood wave can be estimated to be in the range of 500–700 m³/s [56,72]. The simulation results based on the STREAM_2D modeling (710 m³/s) are at the upper boundary of the field assessment. Some estimates of the outburst flood discharges in the Adylsu valley [56], for example, a section of the Shkelda bridge, were made taking into account debris transport (800 m³/s) and exceeded the modeled values, since the modeling did not take into account the debris component of the discharge. Considering the tasks mentioned in the Introduction, we provide only estimates for the water component. In the future, it will possible to expand the work, taking also into account debris entrainment, transport, and accumulation.

Two possible scenarios of the lake re-outburst were considered, and according to the first one, connected with the impact of a rock avalanche onto the lake (scenario hd1), we did not receive a significant model response. Such results can be explained by the rather small size of the lake. For example, when we applied the same methodology to the Tangjiashan Lake, where the potential rock fall impacted a deep and narrow part of the lake, a peak outflow discharge 5000 m³/s was received [73].

According to the second scenario of the lake re-outburst (scenario hd2), the maximum discharge is 58% lower than the discharges of the 2017 GLOF. However, a recurring outburst still poses a threat to infrastructure in the Adylsu and Baksan River valleys.

This work has demonstrated the possibility of using a chain of two inherently different models (hydrodynamic and runoff formation model) to assess the impact of outburst floods on the Baksan River flow. As far as we know, in earlier studies of outburst floods, either use the “process chain”, considering only hydrodynamic modeling within the valley [44], or model runoff formation, which was used to assess the effect of the GLOF on river runoff based on a comparison of simulated and actual values of water discharges [54]. In this work [43], the authors refined a semi-analytical flood-wave propagation model that analytically approximates the kinematic wave equation for simulating downstream wave attenuation, and estimates local peak discharge mainly as a function of downstream distance, channel gradient, hydrograph volume, and flood-wave length. The approach proposed in this study provides an opportunity to assess the contribution of GLOFs to peak discharges in major rivers on a regional scale.

Our estimates indicate that the impact of a large GLOF (about 1 million m³) from the tributary to the discharge of the receiving river downstream is significant over a distance of 40–50 km, and can be detected at a distance of up to 70 km from the mouth of the tributary. A repeated outburst flood of a smaller volume, according to our estimates, will affect only the water discharge at the place where the Adylsu River flows into the Baksan, while downstream of the Baksan River its contribution is insignificant.

The involvement of a runoff formation model also allows the assessment of the joint impact of outburst floods and possible climatic changes on river runoff and peak discharge. Here, we assess how delta-changes in meteorological conditions based on data from 2016–2017 impact the runoff of the Baksan River and the resulting peak discharge of a joint rainfall–runoff and GLOF event. Our calculations showed that the contribution of changing meteorological conditions to the Baksan River runoff is seasonal. In winter, the runoff increases from 10% at Elbrus to 7% at Zayukovo. In summer, mean runoff, on the contrary, decreases by 10% at Elbrus and by 20% at Zayukovo. The maximum water discharges can decrease by up to 20% compared to current ones with an increase in temperature of 2 °C and a decrease in precipitation of 10% in summer. Early studies of runoff changes in the Baksan River basin for long-term period [51] were based on assumption that air temperature and precipitation will increase in a future, when modern climate model estimations [1,63] have indicated a possible decrease in summer precipitation for the North Caucasus in future. At the same time, in the study [51], different scenarios of deglaciation were considered, and it was shown that further warming leads to a stronger runoff increase for the scenarios with modern glacier cover. Without glaciers, summer runoff drops below the baseline hydrographs; the water shortages in July, August, and September equal −10%, −20%, and −18%, respectively [51]. Therefore, it is very important for the next stage of our research to take into account long-term changes to glaciers in the basin, because even current trends already indicate a 6% decrease in the maximum annual flow of the Baksan River since 1981 [58].

In the scenario, with a repeated outburst flood from Bashkara Lake overlapping with heavy rainfall in the Baksan River basin, the rainfall runoff will play the main role in the lower reaches of the Baksan River valley. However, taking into account other outburst-hazardous lakes in the basin a similar research approach and a more detailed analysis of changes in the runoff in the Baksan River basin are required, taking into account climatic changes, reported [4] and projected glacier downwasting, consequent lakes formation, and expansion [103].

6. Conclusions

In the present study, a new version of the two-dimensional hydrodynamic model for the Adylsu River valley was prepared based on STREAM_2D software, which makes it possible to simulate the Bashkara Lake GLOF on 1 September 2017. According to the modeling, the maximum GLOF discharge was 710 m³/s. At the mouth of the Adylsu River, the peak of the flood was observed 40 min after the outburst and the maximum water discharge at the mouth was 320 m³/s. The GLOF travel time to different reaches in the Adylsu valley obtained on the basis of modeling coincides with the observed one.

Two main scenarios were considered as possible mechanisms for a re-outburst at Bashkara Lake: a rock avalanche impact, forming displacement waves, and a lake re-outburst as a result of an increase in the water level accompanied by erosion and expansion of the existing dam break. Under the rock avalanche scenario, there was no significant model response. Based on the results of the modeling of the second scenario, hydrographs of the outburst flood in various sections were obtained: the maximum discharge of the outflow from the lake was 298 m³/s, and the flood peak in the mouth section of the Adylsu River was observed 55 min after the start of the outburst with a maximum water discharge of 101 m³/s.

The ECOMAG runoff formation model was adapted for the Baksan River basin for the first time. ECOMAG calibration and verification using the data from the Tyrnyauz and Zayukovo hydrological stations showed good results, the Nash-Sutcliffe efficiency criteria (NSE) for the verification period exceeded 0.8 for both hydrological stations, and the percent bias of annual runoff volumes was less than 6.5%. The synthesis of the STREAM_2D hydrodynamic model and the ECOMAG runoff formation model in the model chain showed the possibility of combining them to assess the scale of extreme hydrological events in mountainous areas.

The range of possible impacts of the outburst flood on the peak water discharges in the receiving Baksan River was assessed. The greatest influence of the outburst flood is manifested at a distance of up to 50–60 km from the GLOF origin site. The outburst flood contribution to the increase in the Baksan River water discharge on 1 September 2017, according to the modeling, was 45% at the confluence of the Adylsu River and the Baksan River, and 30% near Tyrnyauz (40 km below the confluence). Further, up to 70 km, the contribution of the outburst flood to the peak discharge reduced to less than 20%, and the influence of the rainfall runoff event prevailed.

According to the modeling results, the re-outburst flood from the lakes near the Bashkara Glacier mainly affects the Adylsu River valley and the area of its confluence with the Baksan River; the contribution of the re-outburst flood to the peak discharge at the confluence may be up to 18% and downstream, the influence of the re-outburst flood will be insignificant.

An assessment of the sensitivity of water discharge to changes in meteorological forcing showed that an increase in air temperatures of 2 °C and decrease in summer precipitation of 10% lead to a decrease in maximum water discharge of 20%. Hence, further research to assess the long-term climate change impacts on river runoff in the high-mountainous part of the Baksan River, taking into account changes in glaciation in the basin, is a necessary further step.

The results obtained in this work indicate the need for a comprehensive assessment of both the flood wave movement and the accompanying meteorological situation in the mountain river basins, which should be taken into account both in the analysis of events that have already occurred and in the development of forecast and warning systems.

The obtained results and the adapted model complexes can be used to assess hazard and exposure of extreme hydrological events such as rain-induced and outburst floods in high-mountain glaciated catchments.