Modelling the evolution of Djankuat Glacier, North Caucasus, from 1752 until 2100 AD

. We use a numerical flow line model to simulate the behaviour of the Djankuat Glacier, a WGMS reference glacier situated in the North Caucasus (Republic of Kabardino-Balkaria, Russian Federation), in response to past, present and future 10 climate conditions (1752−2100 AD). The model consists of a coupled ice flow−mass balance model that also takes into account the evolution of a supraglacial debris cover. After simulation of the past retreat by applying a dynamic calibration procedure, the model was forced with data for the future period under different scenarios regarding temperature, precipitation and debris input. The main results show that the glacier length and surface area have decreased by ca. 1.4 km (ca. -29.5 %) and ca. 1.6 km² (-35.2 %) respectively between the initial state in 1752 AD and present-day conditions. Some 15 minor stabilization and/or readvancements of the glacier have occurred, but the general trend shows an almost continuous retreat since the 1850s. Future projections using CMIP5 temperature and precipitation data exhibit a further decline of the glacier. Under constant present-day climate conditions, its length and surface area will further shrink by ca. 30 % by 2100 AD. However, even under the most extreme RCP 8.5 scenario, the glacier will not have disappeared completely by the end of the modelling period. The presence of an increasingly widespread supraglacial debris cover is shown to significantly delay 20 glacier retreat, depending on the interaction between the prevailing climatic conditions, the debris input location, the debris mass flux magnitude and the time of release of debris sources from the surrounding topography.


Introduction
Recently, a lot of attention has been given to modelling mountain glaciers, in particular due to their worldwide observed 25 shrinkage and important role within the current changing climate (e.g. Shannon et al., 2019;Zekollari et al., 2019;Hock et al., 2019). The observed warming trend is a significant matter of concern to scientists and all other people (in)directly involved in the behaviour of these glacial systems, as projected scenarios point towards an even further increase of the global mean temperature in the future, especially if no efficient mitigation strategies are implemented (Vaughan et al., 2013;Rasul and Molden, 2019;Hock et al., 2019). Being consistent with this global trend, the accelerated retreat of Caucasian glaciers 30 during the last several decades has been clearly noticed (e.g. Shahgedanova et al., 2014;Zemp et al., 2015;Tielidze, 2016). Accordingly, total glaciated area has decreased from 691.5 ± 29.0 km² to 590.0 ± 25.8 km² (-0.52 % yr⁻¹) in the period between 1986 and 2014 (Tielidze et al., 2020). Further degradation of Caucasian glaciers may affect the supply of water used for drinking, irrigation and hydroelectric energy generation, whereas it may also pose a threat for downstream communities via flooding, glacier collapses, avalanches, debris flows and glacial lake outbursts (e.g. Volodicheva, 2002;Ahouissoussi et 35 al., 2014;Taillant, 2015;Chernomorets et al., 2018). Furthermore, the presence of glaciers in the Caucasus can be considered important for paleoclimatic research, tourism, cultural heritage and biodiversity (e.g. Popovnin, 1999;Shahgedanova et al., 2005;Hagg et al., 2010;Makowska et al., 2016;Tielidze and Wheate, 2018;Rets et al., 2019). Despite these rising concerns, however, modelling of Caucasian glaciers is scarce and has only been attempted in a few studies (e.g. Rezepkin and Popovnin, 2018;Belozerov et al., 2020). 40 In a warming climate, debris coverage onto the glacier's surface is believed to increase drastically due to the build-up of more englacial melt-out material, lower flow velocities and increased slope instability, which favour the occurrence of rock slides and mass movements from the surrounding topography (Østrem, 1959;Kirkbride, 2000;Stokes et al., 2007;Jouvet et al., 2011;Carenzo et al., 2016). During the last decades, a sharp increase of debris-covered glacier surfaces has been observed over the Caucasus region, owing to the combined effects of steep terrain, a wet climate, small average glacier size, 45 large lateral moraines and the presence of local easily erodible sedimentary rock outcrops. Accordingly, debris coverage has expanded at a rate of ca. +0.22 % yr⁻¹ between 1986 and 2014 when the entire Caucasus region is considered (Tielidze et al., 2020). Scherler et al. (2018) estimate the supraglacial debris cover on Caucasian glaciers to be 26.2 % (ca. 155 ± 6.7 km²) at present-day, hence enabling the area to hold the world's most abundant share of debris-covered glacier surfaces in relative terms. Evidently, the presence of such supraglacial debris can influence the evolution of mountain glaciers in a variety of 50 ways, depending on its thickness, properties and spatial/temporal distribution (Nicholson and Benn, 2006;Anderson and Anderson, 2016). Apart from a slight melt enhancement for a very thin debris layer, thick debris has been shown to reduce runoff volumes and reverse mass balance gradients due to its melt-reducing effect (e.g. Østrem, 1959;Bozhinskiy et al., 1986;Anderson and Anderson, 2016). If a thick supraglacial debris cover is present over a large portion of a glacier's ablation zone, surface melting and terminal retreat can be drastically suppressed, even under a warming climate (e.g. 55 Scherler et al., 2011;Benn et al., 2012). In such cases, debris-covered glaciers are shown to lose mass by lowering the surface in their ablation zone (downwasting), rather than by terminus retreat (e.g. Hambrey et al., 2008;Rowan et al., 2015).
The pronounced effect of debris should therefore not be ignored in numerical experiments to determine the future evolution of mountain glaciers, yet only few studies have included this complex process in time-dependent models (e.g. Jouvet et al., 2011;Rowan et al., 2015;Huss and Fischer, 2016;Kienholz et al., 2017;Rezepkin and Popovnin, 2018;Wirbel et al., 2018). 60 In this paper, we focus on modelling the Djankuat Glacier (North Caucasus, Russian Federation), a WGMS (World Glacier Monitoring Service) reference glacier which has a broad observational network in both space and time. However, despite abundant field data availability and increasing interest concerning its future behaviour, the Djankuat Glacier has not yet been modelled extensively. Here, we present a numerical flow line model to simulate its response to past, present and future of stakes) and indirect (stereophotogrammetrical) measurements, of which the resulting maps are reported in Aleynikov et al. (1999) and Pastukhov (2011). Glacier-wide ice thickness maps have also been constructed by Lavrentiev et al. (2014), using ground-based radio-echo measurements. However, direct and reliable observations lack at the higher elevations (> 3600 m) and the Djantugan Plateau due to difficult accessibility. In these areas, ice thickness values have been derived indirectly using surface velocity and slope (Aleynikov et al., 2002b;Pastukhov, 2011). The current ice thickness has been found to go 100 up to ca. 100 m in the central part of the main glacier body, and to more than 200 m at the Djantugan Plateau. Furthermore, the glacier's cumulative surface mass balance during the 1967/68−2016/17 period exhibited a strongly negative value of -14.33 m w.e., with a mean equilibrium line altitude (ELA) of 3213 m (WGMS, 2018). Moreover, the mass balance profile in the upper areas is significantly modified (at 3600 m by ca. -76 % of the value that the specific mass balance would have if it were extrapolated according to the mass balance gradient found below) by snow redistribution processes (Pastukhov, 2011). 105 Both glacier-averaged debris thickness (from 0.28 m in 1983 to 0.54 m in 2010) and total debris-covered area (from ca. 0.10 km² or 3.5 % in 1968 to ca. 0.34 km² or 12.7 % of the glacier in 2010 AD) have increased largely. However, at the debriscovered left side of the snout (when seen in the downstream direction), debris thickness increased exponentially over the years, resulting in mean values of 1 m at the glacier front in 2010, compared to 0.29 m in 1983 and 0.45 m in 1994 AD . Recent observations have shown the importance of the debris cover on the Djankuat Glacier, as the 110 debris-covered left side of the front clearly retreated slower than the less affected right part. As of 2010 AD, the length difference between both sides was ca. 180 m ( Fig. 1) but this has increased to ca. 250 m by 2017 AD (Rets et al., 2019).
The climate around the glacier can be inferred from nearby weather stations, such as Terskol (elevation 2141 m, approx. 20 km northwest of the glacier) and Mestia (approx. 16 km southwest from Djankuat Glacier, in Georgia, at 1441 m elevation), see Fig. 1 and Table 2. The average mean annual temperatures in Terskol and Mestia are 2.8 °C and 6.0 °C respectively for 115 the 1981−2010 reference climate. For the summer half-year from April to September (AMJJAS), the corresponding mean temperatures are 8.7 °C and 12.0 °C. Precipitation, on the other hand, is rather complex in the region due to variations of atmospheric circulation patterns, orographic uplift and convective precipitation in the summer season (Boyarsky, 1978;Shahgedanova et al., 2007;Hagg et al., 2010;Popovnin and Pylayeva, 2015). At Terskol and Mestia, total annual precipitation amounts equal 1001.1 and 1035.1 mm yr⁻¹ respectively for the 1981−2010 climate. During the accumulation 120 season (October to March, ONDJFM), the corresponding precipitation values are 418.4 and 490.0 mm yr⁻¹ w.e. respectively.
In 2007, two automatic weather stations (AWS) were additionally installed, one in the Adylsu Valley at ca. 2640 m elevation (AWS 1 in Fig. 1) and one in the ablation zone of the glacier at ca. 2960 m on a sparsely debris-covered ice surface (AWS 2 in Fig. 1). During the summer seasons (June to September, JJAS) of 2007−2017, a wide range of additional meteorological variables have therefore been acquired by both AWSs (air temperature, dew point temperature, incoming and outgoing 125 shortwave/longwave radiation, relative humidity, wind speed and direction, air pressure and for AWS 1 also precipitation amounts). The AWSs did, however, not operate outside the JJAS period (Rets et al., 2019).

Ice dynamic model
The ice dynamic model is implemented as a numerical flow line model, in which the prognostic continuity equation for ice thickness change is solved. We choose to only model ice flow along a central axis in the x-direction and not upgrade the 130 model to 3D due to the abundant amount of experiments that were conducted. However, the y-dimension is implicitly taken into account due to inclusion of glacier width along this central axis. One central flow line is considered with a total length of 5 km, stretching from the glacier top near the Djantugan peak down to the current snout and further into the Adylsu Valley (Fig. 1). The flowline is constructed perpendicular to the surface elevation isolines, generally close to the location where the cross-sectional ice thickness and ice velocity are maximal as determined from ice thickness and surface velocity 135 maps. The model treats ice flow as a non-linear diffusion problem in a vertically integrated approach (e.g. Oerlemans, 2001) where is the ice thickness, the time, µ the slope of the lateral valley walls, ? the glacier bed width, ()* the glacier surface width, -*. the ice volume flux, the horizontal distance, 3 the local annual surface mass balance, -the ice 140 density, the gravitational acceleration, G the flow parameter related to internal deformation, ( the flow parameter related to basal sliding and ℎ the surface elevation. The vertically integrated velocity is calculated by assuming that the 1D Shallow Ice Approximation is applicable to derive driving stresses on a xz plane and that ice is treated as a homogenous, incompressible and isothermal non-Newtonian fluid in Glen's flow law. For basal sliding, a simplified Weertman-type flow law is used where the basal water pressure is proportional to the ice thickness and the basal shear stress equals the driving 145 stress (e.g. Oerlemans, 1992;Oerlemans, 2001;Leclercq et al. 2012): Here, is the vertically averaged horizontal velocity, and G and ( are the velocity components related to internal deformation and basal sliding respectively. Equation (1) is then solved on a staggered grid with a spatial resolution ∆ of 10 m. Integration over time is achieved with a forward in time, centered in space (FTCS) numerical scheme using a time step ∆ 150 of 0.0005 years, as determined by the CFL-condition for diffusion problems.

Mass balance model
The mass balance model is based upon the difference between accumulation and runoff over the balance year (1 October -30 September), so that the local surface mass balance 3 is defined as: (3) 155 each point along the flow line is only dependent on the part of the total precipitation that is solid ( (\]-G ), which only takes place if precipitation occurs below a certain threshold temperature _Y.(I : Air temperatures 3-Y from Terskol weather station were interpolated to any height on the Djankuat Glacier by applying 160 vertical temperature lapse rates γ b (Table 1). A direct comparison of measured air temperatures between AWS 2 on Djankuat and the Terskol weather station was found to exhibit a strong correlation (R² = 0.81), generating a summer season lapse rate of -0.0067 °C m⁻¹ between 2007−2017 AD. Due to lack of AWS data outside of the JJAS period, a temperature lapse rate of -0.0049°C m⁻¹ was used for the winter half-year (ONDJFM), in accordance with a mean annual ELA temperature of -3.75 °C for Djankuat Glacier (WGMS, 2018). The term b.Y(c\] * . represents the precipitation in the 165 Adylsu Valley, calculated by multiplying the precipitation in Terskol with a horizontal precipitation enhancement factor . to account for horizontal precipitation variations. In this study, a value for . of 1.5 between Terskol and the Adylsu Valley was found after a comparison of precipitation amounts from AWS 1 in the glacier valley. The factor (*3]. is used to scale the obtained precipitation amounts to the entire glacier from the Adylsu Valley to any surface elevation ℎ, by making use of a vertical precipitation gradient γ j , where the latter is used as a tuning parameter due to a lack of data (see Sect. 3.1): 170 At last, the factor Y.G represents a snow redistribution factor which corrects the solid precipitation for redistribution by wind and/or avalanches. Here, a topographic characteristic is used to parameterize snow addition or removal from the glacier surface. It was quantified by dividing the linear accumulation profile (without the redistribution factor) with the observed profile and correlating these anomalies to the laterally averaged surface slope along the flow line (e.g. Huss et al., 2009).
As such, a polynomial fit was found. For slopes steeper than the threshold, removal of snow can occur, and is assumed to be 175 influenced by the surface slope itself: The critical slope *Y-_ distinguishes between slopes that either favour snow addition or snow removal (Table 1). We do acknowledge that the Y.G parameterization is solely used for curve fitting of the accumulation profile.
Melt production , on the other hand, only takes place when the net energy flux per unit area at the surface ? is positive 180 (e.g. Oerlemans, 2001;Nemec et al., 2009): where } is the water density and • the latent heat of fusion. As discussed further in section 2.5, the melt term is further modified by the debris cover and meltwater retention in the snowpack to obtain the total runoff . The net energy flux is parameterized as (Oerlemans, 2001;Giesen and Oerlemans, 2010;Leclercq et al., 2012): 185 Here, is the atmospheric transmissivity, is the surface albedo, while ? and 5 are constants to describe the air temperature-dependent fluxes (i.e. the net longwave, latent heat and sensible heat fluxes). Hence, for air temperatures below the threshold …Y.3c , ? has a constant value. For higher temperatures, however, ? increases linearly with 3-Y , where the rate of increase is determined by 5 . The downward incoming solar radiation at the surface ↓ 190 incident on an inclined surface with a certain surface slope and aspect, is calculated as (e.g. Oerlemans, 2001): where ↓(b ‡ˆ) is the incoming instantaneous extraterrestrial shortwave radiation on a horizontal plane at the top of the atmosphere, . and • are the solar elevation and zenith angle respectively calculated using basic astronomical formulas (e.g. Iqbal, 1983;Allen et al., 2006;Duffie and Beckman, 2006) and is the angle of incidence, taking into account the surface 195 geometry. Furthermore, G-Y and G-) are the fraction of direct and diffuse solar radiation, which are derived from parameterizations used by Oerlemans (1992Oerlemans ( , 2001Oerlemans ( , 2010 and Voloshina (2002) that use the fractional cloud cover *] : At last, surface albedo is parameterized as (e.g. Oerlemans and Knap, 1998;Nemec et al., 2009): where ("\} is the snow albedo, -*. the ice albedo and ("\} * a characteristic snow depth. Measurements of the incoming solar radiation from the AWS 2 were used to derive atmospheric transmissivity. These data were therefore compared to the theoretical maximum incoming solar radiation at the top of the atmosphere, calculated with standard astronomical formulas (e.g. Iqbal, 1983;Allen et al., 2006;Duffie and Beckman, 2006). Consequently, the overall atmospheric transmissivity in the summer season over the Djankuat Glacier could be deduced as an average of 0.53 (Table  205 1). The ice albedo -*. can, according to raw data from the AWS 2, vary between 0.15 and 0.40 depending on the presence of water, moraine cover and other impurities and has an average value of 0.22, corresponding to moderately debris-loaded ice. Sparse snow-covered conditions during the ablation season causes ("\} to increase to the 0.40−0.90 range (mean 0.79).
Next, values for G-Y and G-) are derived from the parameterization of the fractional cloud cover *] over the Djankuat (Voloshina, 2002), of which the latter was derived from measurements by AWS 2 on the glacier surface. The analysis points out that direct and diffuse solar radiation are approximately equally important for the glacier (Table 1). The constants ? , 5 and …Y.3c , describing the air-temperature dependent fluxes and their relationship with the air temperature 3-Y itself, are at last derived from measurements of AWS 2 of the net longwave radiation, as well as from a parameterization of the sensible and latent heat fluxes via Kuzmin's method (Kuzmin, 1961;Toropov et al., 2017). Here, these fluxes are added up and 215 analyzed against air temperature following the method of Giesen and Oerlemans (2010) and Leclercq et al. (2012), as can be seen from Eq. (8).

Debris cover model
The supraglacial debris cover on the Djankuat Glacier was parameterized in order to account for the effects of melt reduction under debris-covered ice. The debris thickness was approached with a steady deposit model adopted from Anderson and 220 Anderson (2016), where debris input onto the glacier is generated from a fixed point on the flow line. In the model, debris thickness then changes according to either melt-out from debris-loaded ice (first term), the downstream advection of supraglacial debris (second term) and the input or removal of supraglacial debris on the glacier surface (third term): Here, G.…Y-( is the debris thickness, the time, G.…Y-( the englacial debris concentration, G.…Y-( the debris cover porosity, 225 G.…Y-( the debris rock density, 3 the specific surface mass balance, ()* the glacier surface velocity and G.…Y-( the input or removal of debris from the glacier surface. The advection equation (Eq. 12) is solved using a first order upwind scheme with ∆ = 0.01 years, in accordance with the CFL-condition for advection problems. In the model, the factors G.…Y-( and G.…Y-( are constants in space and time and taken at 0.43 and 2600 kg m −3 respectively (Bozhinskiy et al., 1986). For G.…Y-( , we use a value of 1.05 kg m −3 , referring to a bulk debris concentration inside the ice of 0.12 % as found by the same authors for the 230 Djankuat Glacier in the 1980s (Table 1). Also here, a constant value in space and time is assumed. Incorporating englacial debris pathways or the spatial distribution of englacial debris concentration would add more detail than warranted by the lack of reliable data regarding this value.
Next, at the debris input location G.…Y-( , a steady debris flux per unit area G.…Y-( -"ª¥_ transmits material from the surrounding topography to the glacier by means of a debris deposition rate (m yr⁻¹), starting from G.…Y-( onwards. Here, G.…Y-( is defined 235 as the time at which the topographic debris source firstly starts to release its mass flux towards the glacier surface. We set the debris input location G.…Y-( at 1680 m from the highest point (just below the ELA, at 88% of the distance between the terminus { and the ELA «{ˆ) , since it is the furthest point up-glacier for which observed debris thickness values are reported in . It was chosen to keep the debris input location at a fixed position due to the general absence of direct observations regarding past (static or moving) topographic debris sources. However, a comparison of 240 present-day satellite imagery with those from the 1970s (Pasthukov, 2011) points out that the debris patches exhibited only minor up-glacier migration on the main glacier tributary and the debris-covered part of the snout, lending some support to this assumption.
To avoid the buildup of unrealistically high debris thickness in low flow velocity zones in the future, we furthermore choose to let the debris mass flux stop when the surface width at point G.…Y-( has reached a value lower than 90 % ( 6()*-5? % ) of 245 its original value at time G.…Y-( . This is considered a reasonable value, as the current observed debris-covered area is ca. 10 % at this specific point (Fig. 3). Connectivity issues between the topographic source and the main glacier are forwarded as the main reason to justify this modification of the Anderson and Anderson (2016) model. Consequently, by then the glacier has laterally shrunk too much to ensure that debris fluxes could still reach its surface. At the terminus (the last non-zero ice thickness grid point), debris is removed into the foreland by a debris flux per unit area G.…Y-( =-{ (Anderson and Anderson, 250 where { is the terminus position and { is a constant describing the strength of debris removal from the terminus into the foreland, for which we used the same value as suggested in Anderson and Anderson (2016), i.e. { = 1 ( Table 1) Here, G.…Y-( * is a characteristic debris thickness (i.e. the debris thickness at which the melt rate is e⁻¹ or ~ 37% of the clean ice melt rate). It must be noted that the melt enhancement that may occur for a very thin debris cover was not implemented in 260 the debris model. However, values in literature of the debris thickness for which a maximum amount of melt enhancement occurs on the Djankuat Glacier vary from 0.02 m to 0.07 m (Bozhinskiy et al., 1986;Popovnin and Rozova, 2002;Lambrecht et al., 2011), and the areal fraction of Djankuat Glacier that holds these thin thickness values is very small (Popovnin and Rozova, 2002;. It is therefore not believed to have a significant influence on the ablation of Djankuat Glacier. 265 Next, the fractional debris covered area along the flow line is parameterized based upon the distance from the terminus b , for which an exponential relationship was found from observations that can, of course, not exceed 1: is furthermore worth noting that the debris model also neglects other processes that may potentially play a role in the spatial 270 and temporal distribution of debris, such as the formation and thickening of medial moraines, ice cliffs and surface ponds (Anderson and Anderson, 2016).
As such, in the case that snow is present at the glacier surface, runoff is calculated as the meltwater outflow from a saturated snowpack ("\} , following the principles applied in Schaefli and Huss (2011). On the other hand, in case of snow-free conditions, runoff is affected by the presence of a debris cover on the glacier ice (e.g. Lambrecht et al., 2011): 275 where is the melt production (see Sect. 2.4), ("\} is the water outflow from the saturated snowpack, ("\} the liquid snow store, η ( the water holding capacity of the snowpack, G.…Y-( the melt-reduction factor from debris, G.…Y-( the debris covered area and ("\} the snow depth.

Mass balance model
We used the 1967/68−2006/07 period to calibrate the mass balance model, as this time frame holds both specific (elevationdependent) and mean specific (glacier-wide) surface mass balance measurements (WGMS, 2018). Accordingly, 3-hourly temperature and precipitation data of the corresponding period were used from the Terskol weather station. For geometric data that serve as input for solar geometry calculations, we use laterally averaged values for slope and aspect, calculated by 285 averaging all intra-glacier values along a line perpendicular to the flow line. Surface elevations were directly extracted from a DEM for 2009/10 AD conditions. We hereby take into account the same spatial spacing of 10 m that is used in the flow model. Afterwards, geometric input data were smoothed using a window size of ± 100 m around every grid point.
Calibration of the mass balance model further assumes the geometry (slope, aspect, glacier length and surface area) to be fixed over the 1967/68−2006/07 period, whereas in fact length and surface area decreased by 113 m and 0.346 km² 290 respectively.
For the accumulation part, the vertical precipitation gradient γ ¾ was used as a tuning parameter by fitting the accumulation profile of the glacier. In literature, several values for this parameter have been proposed, varying between 0.0005 and 0.0046 m yr⁻¹ w.e. m⁻¹ (e.g. Boyarsky, 1978;Hagg et al., 2010;Giesen and Oerlemans, 2012;WGMS, 2018). To ensure successful calibration, a precipitation gradient of 0.0023 m yr⁻¹ w.e. m⁻¹ was derived to extract these data over the entire glacier surface. 295 At last, the snow redistribution factor Y.G was used for curve fitting of the accumulation profile, as discussed before.
holding capacity of snow η ( , it was used to calibrate the ablation in the accumulation area. Additionally, the intercept of the air temperature-dependent fluxes c ? was chosen as a second tuning parameter, again due to the lack of reliable and/or sufficient data during the observational period (Table 1). For the factor G.…Y-( * , which controls the strength of the melt-300 reducing effect of debris, several values have already been proposed for the Djankuat Glacier (e.g. Bozhinskiy et al., 1986;Popovnin and Rozova, 2002;Lambrecht et al., 2011). Due to the large uncertainty, it was used as a third tuning parameter, this time for the lower elevation areas. Here, a value of 1.15 m was found to exhibit the best fit with the observations. We acknowledge that this value implies that the gradient of the exponential decay in Eq. 14 is somewhat out of range with respect to earlier studies for other glaciers (e.g. Anderson and Anderson, 2016). This rather atypical value can however be 305 linked to the relatively high thermal conductivity of the granite-type debris cover on the glacier (2.8 W m −1 °C −1 ) and the high debris cover porosity (0.43 for Djankuat Glacier, Bozhinskiy et al., 1986). Also, the relatively low water saturation and large particle size, as suggested by Lambrecht et al. (2011), may imply that heat conduction towards the debris-ice interface seems to occur efficiently on the Djankuat Glacier.
With the calibrated surface energy balance model, the multiyear mean mass balance profile of the Djankuat Glacier during 310 the 1967/68-2006/07 period is successfully reproduced, as the calculated mass balance vs. elevation profile matches nicely with the observations (Fig. 4). This profile reflects the determining processes affecting the Djankuat Glacier's mass balance: As a remark, it must be noted that the calibration dataset for the mass balance model is quite long (39 years from 1967/68 to 2006/07 AD), making it credible to assume that the parameters calibrated to this period have some validity for past and 320 future conditions as well. Apart from the high-elevation areas, where data availability is limited and snow redistribution processes create complex conditions (> 3600 m, of which the areal fraction is only ca. 3 % of the glacier area in 2010 AD), it can be expected that the environmental setting within the calibration window also holds for periods prior to and after the observational period. It must furthermore be noted that there are only few independent data to validate our model results with a sufficient degree of certainty. 325

Debris cover model
For the debris model calibration, we matched the temporal evolution of the average debris thickness at the front (i.e. the first 30 grid points) as well as the debris covered area, using G.…Y-( , G.…Y-( -"ª¥_ and ˆ as tuning parameters. Values for the observed debris cover at different elevation bands from the survey year 1968 AD (only for debris area) as well as for 1983, 1994 and 2010 AD (for both debris area and thickness) are available from . Moreover, to obtain more detailed 330 information concerning the current debris covered area on a spatial scale, the debris cover extent was manually digitized based on satellite imagery of the year 2010 (see Fig. 1).
Accordingly, the observed debris thickness evolution was found to be best reproduced by setting G.…Y-( to 1958 and G.…Y-( -"ª¥_ to 1.60 m yr⁻¹ (Table 1). At last, a power relation (R² = 0.85) was found between the growth factor ˆ and the modelled mean debris thickness at the glacier front as obtained in the previous step: 335 Where G.…Y-( )Y\"_ is the modelled debris thickness at the front (i.e. the first 30 grid points) as obtained before. As such, the RMSE between modelled and observed values between 1967/68 and 2009/10 AD was reduced to 0.07 m (R² = 0.83) for debris thickness at the front and 0.9 % (R² = 0.95) for the fractional debris-covered area respectively (Fig. 5).

Ice dynamics model 340
To calibrate the flow model, it was initially run from zero ice thickness until a steady state situation was reached, which is achieved when the glacier has less than 0.002 % change in its total volume per year. The steady state situation of the ice flow model was then tested by comparing the ice flux with the integrated upstream mass balance, by ensuring that the integrated surface mass balance over the entire glacier approaches 0 to within an acceptable accuracy (0.006 m yr⁻¹ w.e.), and by calculating the volume change with time. As expected for the model setup, all results exhibited an appropriate steady state 345 situation for the glacier. The parameters G and ( were adopted to minimize the RMSE between observed and modelled ice thickness for 2010 AD conditions, assuming a steady state. Geometric input data for the flow model were extracted from a DEM for 2010 AD conditions. Hence, bedrock elevation was derived in combination with ice thickness maps from Pastukhov (2011) andLavrentiev et al. (2014). Surface width was extracted by measuring the intra-glacier distance of 10-m spaced lines perpendicular to the orientation of the flow line. After extracting the lateral valley slopes, the width at the bed 350 was calculated assuming a trapezoidal valley shape (e.g. Oerlemans, 1992;Gantayat et al., 2017). All data were finally joined to the closest point on the flow line for every 10 m and smoothed with a window of ± 100 m around every grid point.
Additionally, the bed width for the assumed trapezoidal-shaped cross section was slightly adjusted to ensure that the parameterization fits the observed area-elevation distribution for a total surface area of 2.688 km². The full set of parameter 355 values used in the model is given in Table 1.
The flow model for the Djankuat Glacier was able to produce a steady state glacier profile with a length of 3.26 km after 200 years (Fig. 6). The model approaches the observed ice thickness as it minimizes the RMSE to 14.27 m (R² = 0.90). Despite minimized RMSE, the mismatch near the snout and steep slopes of the Djantugan peak increase the error of the model. presence of a supraglacial debris cover at the front, or the lack of reliable and direct ice thickness observations at the highest elevations of the glacier. As with the mass balance and debris cover model, there are no, or only few, independent data to validate our model results with a sufficient degree of certainty.
Modelled current surface velocity for the Djankuat Glacier goes up to 79.7 m yr⁻¹ near the ice falls of the Djantugan Plateau and also peak in the middle section of the glacier, which fits well with observations of maximum velocities in the 60−80 m 365 yr⁻¹ range (Aleynikov et al., 1999;Pastukhov, 2011). Moreover, the modelled deformational and basal sliding components comprise respectively 45 % and 55 % of the vertically averaged ice flow velocity along the flow line.

Basic sensitivity experiments
With the calibrated submodels, some basic sensitivity tests were conducted which all initially started from a steady state glacier resembling the present-day geometry. Perturbed mass balance profiles (in steps of 0.25 m yr −1 w.e.) were 370 subsequently used as forcing into the flow model, until a new steady state was reached. As such, a relationship with a slight deviation from linear was found between the steady state length and the mass balance perturbations ∆ 3 , exhibiting a value of ca. 1100 and 1355 m (m yr −1 w.e.) −1 for negative and positive perturbations respectively (Fig. 7a). On the other hand, the e-folding length response time (i.e. the time needed to achieve 1 -e -1 or ~63% of the total length change) of Djankuat is in the order of 31 ± 3 years. Additional sensitivity experiments with the mass balance model show that the Djankuat Glacier, 375 when its 2010 AD geometry and other parameters are considered fixed, is quite sensitive to both temperature (-0.70 m yr⁻¹ w.e. °C −1 ) and precipitation changes (0.20 m yr⁻¹ w.e. 10 % −1 ). As such, a 1 °C annual temperature change for the Djankuat Glacier is only compensated when the precipitation change is in the order of ca. 35 %. Mass balance sensitivity to temperature changes shows a non-linear behavior, whereas the relationship is linear for precipitation changes (Fig. 7b).
To assess the climate and glacier sensitivity for equilibrium conditions, mass balance profiles were furthermore altered by 380 temperature and precipitation perturbations within the -3 to +3 °C and -25 % to +25 % range respectively (as compared to the 1967/68−2006/07 AD reference values). Sensitivity of steady state length to temperature changes was found to exhibit a linear behaviour (815 m °C −1 ) for perturbations between -1.4 and +0.7 °C, but is modelled to vary between 400 and 1400 m °C −1 when assessed over the entire range (Fig. 7c). The glacier sensitivity depends largely upon geometry and increases (decreases) for more negative (positive) mass balance perturbations, predominantly due to the flatter (steeper) terrain. The 385 sensitivity also peaks around a temperature perturbation of +1 °C, i.e. when the glacier front is positioned at the transition between the broad accumulation area and the narrower snout (ca. x = 2300 m on the flow line). Also, the non-linear nature of the temperature-mass balance relationship (Fig. 7b) triggers a deviation from linear behaviour. Consequently, the change in forcing needed for a retreat from 2 to 1 km is nearly twice as large as for a retreat from 4 to 3 km. For precipitation the sensitivity is more or less constant for a value of 250 m 10 % −1 (Fig. 7d). A temperature increase of +3.4 °C compared to the 390 Djantugan Plateau melts away 470 years after the induced perturbation.

LIA extent of the glacier
All three submodels (ice flow, mass balance and debris cover) are finally coupled to determine the past and future evolution 395 of the Djankuat Glacier. Here, the mass balance model and debris cover model calculate annual surface mass balance profiles, which are then used as input for the continuity equation in the ice flow model after conversion to ice equivalents.
Glacier length is calculated by multiplying the number of non-zero ice thickness grid points by ∆ . is thus not necessarily equal to the glacier terminus position, as the glacier may disintegrate in several sections during retreat. As a first step, the model is initialized with a spin up run in which a steady state glacier, as well as a steady state debris cover, are 400 produced for the balance year 1752/53 AD. Although we have no clear indication to suspect steady state behaviour at this time due to lack of reliable data on debris cover, mass balance and length change, it was imposed to start the simulations without unwanted behavior at the initial stage.
We choose to let the glacier grow until the length indicated by the end moraine of the 19ᵗʰ century (4.62 km), as determined by lichenometric dating in the paleovalley (Boyarsky, 1978;Zolotarev, 1998;Petrakov et al., 2012), cf.  (Fig. 6d). 410 We furthermore chose to include a supraglacial debris cover in the initialization procedure. However, as can be deduced from the large lateral moraines in the Adylsu Valley (Fig. 1), the Djankuat Glacier used to export most of its debris to the margins in the historic period, rather than developing a supraglacial debris cover. Furthermore, debris sources from surrounding topography and melt-out processes were likely less widespread in the historic period because of the colder climate (i.e. the current exposed slopes were covered by the glacier itself and were more stable). Also, the fast-flowing 415 nature of the paleo-glacier tongue in the valley (up to 100 m yr -1 around 1752 AD, Fig. 6d) disfavours the accumulation of thick debris on the glacier surface. For this reason, supraglacial debris is believed to have been much less widespread prior to the observational period of 1967/68 AD, implying that the glacier was not very much influenced by debris cover in the historic period. Nevertheless, there is also indirect evidence for at least some supraglacial debris in the historic period from the presence of moraines in the valley (Fig. 1) and a photograph taken around 1930 showing some debris patches on the 420 snout (Aleynikov et al., 2002b). It would be furthermore unrealistic to only introduce a debris cover in the model once the model approaches the start of the observations. This would contradict the presence of moraines and the observation that there already was an expanding debris cover during the first data collection in 1967/68 AD . Because there is no direct evidence for the origin of the debris cover, it was chosen to include only melt-out processes in the model initialization, which implies that debris mass fluxes from surrounding topography are not incorporated in the initialization 425 procedure (i.e. G.…Y-( -"ª¥_ = 0 m yr⁻¹ and G.…Y-( = 1.05 kg m −3 ). With these values, the LIA steady state debris cover had a thickness of 0.64 m at the front and occupied a fractional area of ca. 8 % (ca. 0.331 km 2 of the 1752 AD glacier).

Evolution of the glacier from 1752 AD to present
To force the model in the historic period, climatic data at 3-hourly intervals are needed. Historic climatic datasets for Terskol weather station were therefore constructed using a multiproxy approach, including information from various weather stations 430 in the area, such as Mestia, Pyatigorsk (approximately 100 km northeast from the glacier at 512 m elevation) and Mineralnye Vody (approximately 115 km northeast from the glacier at 321 m elevation). Additionally, historic data from the CRUTEM4 and CRU TS datasets, as well as from tree ring reconstructions for the broader Caucasus area, were used for the remaining uncovered data gaps since 1752 AD (D'Arrigo et al., 2001;Toucham et al., 2003;Griggs et al., 2007;Köse et al., 2011;Jones et al., 2012;Harris et al., 2014;Holobâcă et al., 2015;Martin-Benito et al., 435 2016;Dolgova, 2016). Data from the pre-observational period outside the Terskol time series were therefore averaged over all the available datasets to create a multiproxy mean time series, for which mean monthly temperatures and total precipitation amounts were derived by matching the mean (corrected additively for temperature and multiplicatively for precipitation) and standard deviation of the overlapping part in the observed Terskol dataset (Table 2, e.g. Huss and Hock, 2015;Zekollari et al., 2019). To obtain a record with a 3-hourly temporal resolution, the data sequence for Terskol over 440 which measurements with a 3-hourly interval are available is repeated into the past and future in order to maintain intra-daily and intra-annual variability in the data. These data were afterwards corrected for the monthly mean temperature and precipitation amounts obtained in the previous step ( Table 2). The reconstruction of temperature and precipitation clearly indicates a shift in the climatic conditions after 1752 AD. Especially during the last few decades, an accelerated warming trend has occurred, as the latest 10-year climatic interval exhibits a mean annual temperature anomaly of +0.5 °C compared 445 to the 1981−2010 mean. This makes it the warmest period in the reconstructed time series. For temperature, a clear sequence of colder and warmer intervals can be seen. Changes in precipitation show a sequence of drier and wetter periods (Fig. 8).
After using the steady state glacier of 1752AD as an initial input feature for the time-dependent model, dynamic calibration is applied by incorporating artificial mass balance perturbations ( (_) ) into the model. This factor was not explicitly calculated but was instead derived and adjusted iteratively by a trial and error procedure. The obtained perturbations were 450 then superimposed on the mass balance profile that was simulated with the climatic input, until the reconstructed glacier length sufficiently matched with the observed values (e.g. Oerlemans, 1997;Zekollari et al., 2014): perturbation that was applied in the dynamic calibration procedure. Such a procedure is needed to counteract imperfections 455 in the flow model, mass balance model and the climate forcing. The added value of this procedure is to ensure a current glacier state that matches the observed one, as the glacier is still responding to changes in past climate, geometry and dynamics. The procedure required a maximum additional mass balance perturbation of +0.5 m yr⁻¹ w.e. but which varies over time (Fig. 9b). Nevertheless, since the balance year 1967/68 AD, i.e. the year from which the mass balance model was calibrated, no additional perturbations were needed. It can thus be stated that the model performs well when forced with the 460 observed Terskol climatic data, and that credibility can be assigned to the dynamic calibration procedure. It furthermore implies that future projections are no longer influenced by the corresponding artificial mass balance corrections, keeping in mind an e-folding length response time of ca. 31 years for the Djankuat Glacier (see Sect. 4).
The resulting mass balance series shows clear peaks around the 1870−80s, early 1900s, late 1910s, 1940s, 1970s and early 2000s AD, hereby coinciding with slightly colder and/or wetter periods in the climatic datasets (Fig. 9c). Clear minima in the 465 mass balance series can be noted in the 1860s, 1890s, early 1910s, 1920s, late 1940s and in the 21ᵗʰ century, which agrees fairly well with earlier mass balance reconstructions of Djankuat (Dyurgerov and Popovnin, 1988;Fyodorov and Zalikhanov, 2018) and Garabashi glaciers on the Elbrus massif (Rototaeva et al., 2003;Dolgova et al., 2013). As the Djankuat Glacier reacted to these climatic perturbations, an almost continuous retreat since the 1850s AD had occurred, exhibiting some minor readvances or steady states as well. As was already discussed earlier, the past behaviour of the 470 Djankuat Glacier is in line with the general observed trend for other Caucasian glaciers (Fig. 2). During the last several decades, however, the addition of a thickening layer of supraglacial debris on the snout aided to temporarily postpone rapid retreat and more or less maintain steady state conditions. Still, the glacier has lost a total length of 1.39 km at present-day compared to the start of the reconstruction in 1752 AD (-29.4 %). The reconstruction also shows that the total glacier area around 1752 AD decreased by 35.2 % when compared to the 2010 AD situation (an area of 4.147 against 2.688 km 2 , see 475 Figs. 6 and 9a). Moreover, evolution of glacier surface area matches nicely with observed values except for the outlier around 1983, which has to do with a migrating ice divide on the Djantugan Plateau (Fig. 9a).
A historic model run conducted with a 100 % clean-ice glacier, shown as an inset in Fig. 9a, revealed that debris played only a minor role prior to ca. 1980 AD, with length differences of only 20 to 40 m. By 2010 AD, however, the modelled length difference between a debris-free and debris-covered glacier already increased to 160 m (Fig. 9a). 480 6 Future glacier evolution to 2100 AD

Response to future climate forcing
Future projections of temperature and precipitation were obtained by a multi-model approach, using output from the Coupled Model Intercomparison Project Phase 5 (CMIP5) simulations (Taylor et al., 2012) for the grid cell closest to the Djankuat 485 Glacier. Mean temperature and total precipitation amount at monthly resolution from 21 Global Circulation Models (GCMs) for the RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5 scenarios were used, based upon their availability (Table 2 and 3). The data were downloaded for both historical runs (from 1981 AD) and for projections (until 2100 AD). Although the choice of ensemble member can largely influence the eventual results (e.g. Huss and Hock, 2015), we solely focus on the first realization, i.e. ensemble member r1i1p1. As with the historic climate datasets, climate data were scaled to match the mean 490 and standard deviation of the Terskol meteorological station. Absolute GCM data were therefore at first scaled to anomalies with respect to the 1981−2010 reference values for each respective model, so that additive (temperature) and multiplicative (precipitation) biases could be removed when matching to the past forcing. For each RCP, the monthly temperature and precipitation data were then averaged over all models, resulting in a multi-model mean time series. To account for year-toyear variability, the CMIP5 data were rescaled with respect to the standard deviation of the overlapping period for the 495 observed Terskol data (e.g. Huss and Hock, 2015;Zekollari et al., 2019). As with the past, the observed 3-hourly Terskol data sequence was finally used to downscale the monthly data to the temporal resolution that suits the mass balance model.
Concerning debris cover evolution, the debris input location G.…Y-( and flux magnitude G.…Y-( -"ª¥_ were left unchanged. Consequently, once the contribution from G.…Y-( stops, either due to shrinkage of the surface width or rapid retreat beyond the input location, no additional debris source is released. Hence, only melt out from debris-loaded ice and supraglacial 500 debris advection contribute to the evolution of the supraglacial debris cover. Later on, we will, however, conduct several experiments to determine the impact of potential additional debris sources from the surrounding topography on the future glacier evolution (Sect. 6.2).
All scenarios exhibit a further increase of the temperature, which is most pronounced in the summer season. Projected precipitation, on the other hand, shows slightly decreasing values at annual resolution, but shows a tendency for a drier 505 summer half year (April to September, AMJJAS) and a wetter winter half year (October to March, ONDJFM). By 2071−2010 AD, the mean AMJJAS temperature (total ONDJFM precipitation) anomalies with respect to the 1981−2010 period are +1.4°C (+0.1 %), +2.3°C (+3.7 %), +2.7°C, (+11.2 %) and +4.5°C (+11.7 %) for the RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5 scenarios respectively (Figs. 10a and b). Additionally, also a future projection is made under a no change scenario, in which the last observed 10-year climatic interval (2009−2018 AD) is repeated with respect to its mean 510 (corresponding to a AMJJAS mean temperature and a total ONDJFM precipitation amount anomaly of +0.5 °C and -11.0 % mm yr −1 w.e. respectively).
All future scenarios agree to a rapid decline of the glacier length and surface area in the following decade, as a response to the significant warming since the late 1990s AD. By 2100 AD in the no change scenario, the total length and surface area of the glacier are projected to be 2370 m (-29.3 %) and 2.01 km² (-27.3 %), whereas the glacier front will be positioned at an 515 elevation of 2844 m. It is thus clear that, at present day, the Djankuat Glacier is not in equilibrium with the current climatic conditions and hence will strive towards a new steady state with a much smaller surface area in the future. For the RCP 2.6, 4.5, 6.0 and 8.5 scenarios, the total glacier length further decreases to 1560 m (-52.2 %), 1250 m (-61.7 %), 1070 m (-67.2 %) and 510 m (-84.4 %) by 2100 AD respectively. Meanwhile, total glacier surface area decreases to 1.17 km² (-56.5 %), 0.71 km² (-73.6 %), 0.49 km² (-81.8 %) and 0.20 km² (-92.6 %) by 2100 AD respectively ( Fig. 11a and b). As such, for the 520 RCP 6.0 and 8.5 scenarios, the glacier retreats back as far as into the bedrock depression of the Djantugan Plateau.
With respect to total runoff volume changes and water resources management, the Djankuat Glacier is close to surpassing its peak water discharge point, as the modelled annual glacier runoff reaches its maximum around 2020 AD (Fig. 11c). Hence, all RCP scenarios exhibit a further decline of the produced runoff volume into the future, which is in accordance with earlier work for this area (Huss and Hock, 2018;Hock et al., 2019). The actual course of runoff changes, however, is dependent 525 upon the trade-off between remaining glacier surface area and magnitude of melt. As such, the RCP 8.5 scenario initially produces the highest melt and corresponding runoff volume. Later on, however, the 'no change' scenario yields the highest runoff volumes due to the larger remaining glaciated area. It must also be noted that near the end of the modelling period, runoff volume temporarily stabilizes for the RCP 6.0 and RCP 8.5 scenarios. This process is related to the melting of the ice on the Djantugan Plateau, which then reinforces itself due to the mass balance-elevation feedback. 530 However, even under the most extreme RCP 8.5 scenario, the glacier would not completely disappear by the end of the modelling period. Despite accelerated melting of the high-elevation plateau because of the mass balance-elevation feedback, a decreased climate sensitivity due to the steeper laterally averaged slopes in the upper glacier part, as well as the large ice thickness on the Djantugan Plateau (up to 200 m at present-day), prevent a complete disappearance by the end of the modelling period. It must furthermore be noted that the averaging of the future climatic data implies a reduction of spread. 535 When, for example, the model was forced with the highest warming scenario of all CMIP5 models (i.e. the RCP 8.5 scenario of the GFDL-CM3 model, with mean AMJJAS temperature increase of +7.9 °C by 2071−2100 AD), the glacier will cease to exist by 2086 AD.

Impact of supraglacial debris cover on glacier evolution
Despite present-day areas of visible clean ice on the tongue, a relatively steep slope below the ELA, relatively high ice 540 velocities and a short response time, observations also show that the supraglacial debris cover on the Djankuat Glacier has significantly affected glacier geometry during the last several decades, as evident from the differential retreat of the snout ( Figs. 1 and 9a). Its importance for this specific glacier has also been demonstrated by e.g. Rezepkin and Popovnin (2018), who showed that the debris cover is believed to drastically affect the Djankuat Glacier in terms of its geometry and melting patterns. Debris input onto the Djankuat Glacier's surface due to mass fluxes from surrounding topography are furthermore 545 expected to increase even further in the future Rezepkin and Popovnin, 2018). To determine the potential effect of these additional debris sources onto the glacier surface, we performed additional experiments with varying debris input location, debris input magnitude and time of the release of the debris source from the surrounding topography.
We repeated the procedure used in Sect. 2.5, but indicate a 'debris reference scenario', in which a second debris mass flux is initiated from G.…Y-( = «{ˆ at G.…Y-( = 2035 with a magnitude of G.…Y-( -"ª¥_ = 1.5 m yr⁻¹. For «{ˆ, the average position of the 550 ELA was calculated during a window of ±15 years surrounding G.…Y-( in the 'no additional debris scenario' (Sect. 6.1), which hence varies for each climatic scenario. We therefore choose to not initiate debris fluxes from positions above the ELA, due to the neglect of englacial pathways in our debris model (see Sect. 2.5). We then let one of these three variables change, while keeping the other two at their original value of the 'reference situation'. As such, the debris input location G.…Y-( was changed to 80 %, 60 % and 40 % of the distance between «{ˆ and { (further downstream), the time of release 555 G.…Y-( to 2045, 2055 and 2065, and at last the magnitude of the debris flux G.…Y-( -"ª¥_ to 0.75, 2.25 and 3.0 m yr⁻¹. It must be noted that the values of these parameters are arbitrary, as the exact location, time and magnitude of future debris sources cannot be predicted. By assessing a range of possible values for each of these parameters, we encompass various potential future scenarios in order to account for the high uncertainty regarding these parameters. Figure 12 shows the impact of these variables (rows) on the future length of the Djankuat Glacier under different climatic scenarios (columns). The black lines 560 indicate the scenario where no additional debris source is released in the future. The other lines are for experiments that include an additional future debris source from the surrounding topography for varying values of the earlier mentioned debris-related parameters. It is clear that the addition of an increasingly widespread debris cover dampens glacier retreat. It should however be noted that the effects on glacier length are not immediate, as it takes some time for the debris to be advected to the terminus after its initiation at time G.…Y-( . 565 The effect of the timing of the source release is straightforward: the earlier the debris mass flux is released, the larger the extension of the glacier by the year 2100 AD, as the melt-reducing effect starts earlier in time. The main decisive factor here is the efficient debris advection towards the terminus, because flow velocities are larger in 2035 AD compared to 2050, 2065 and 2080 AD (Fig. 12). The magnitude of the debris input flux G.…Y-( -"ª¥_ is another crucial parameter determining the length extension of the Djankuat Glacier in the future period. It is, hence, obvious that a higher flux magnitude will contribute more 570 efficiently to a higher debris growth rate. This enhanced effect is a direct consequence of the implementation of Eq. 14, where the debris-related melt reduction depends on the debris thickness (Fig. 12). Concerning the debris input location, results suggest that the closer the input source is located to the terminus, the longer the extension of the glacier will be compared to the situation without an additional debris source. This makes sense, as the time that it takes for the supraglacial debris to be advected to the front is shorter for down-glacier input locations. Hence, the debris cover will be able to apply its 575 melt-reducing effect much earlier in time, as well as much further down-glacier in space on a still relatively long glacier.
Again, the effects on glacier length are not immediate, as it takes some time for the debris to be advected to the terminus.
The effect of climatic conditions on debris-related melt reduction and its impact on glacier geometry is twofold. Initially, the melt-reducing effect increases with higher temperature, as can be seen in the case of the no change, RCP 2.6 and RCP 4.5 scenarios. This can be related to the fact that a higher temperature will increase the melt-out of material from debris-loaded 580 ice, whereas also decreased flow velocities prevent sufficient discharge and allow the debris to thicken quickly up-glacier ( Fig. 12). Moreover, the distances between the input point and the glacier front at the time of source release decrease with increasing temperature, whereas also retreat rates are relatively larger for higher temperatures. This allows the relatively thick debris to encounter the glacier front much earlier in time. At last, it is important to note that for the same melt reduction However, for the RCP 6.0 and RCP 8.5 scenarios, the impact of the supraglacial debris cover on the glacier decreases again.
Here, a counteracting effect occurs as temperatures rise even further, because the risk of rapid loss of debris-covered area increases. This can be related to either the breaking of the glacier into several fragments where areas of 'dead ice' prevent proper connectivity between the main glacier body and the glacier front, or because the front is too close to (or has already passed) the debris source by the time it is released. Finally, the accelerated shrinkage also favors foreland deposition instead 590 of debris accumulation due to frontal retreat, as well as the loss of proper connectivity between the debris source and the main glacier body at the debris input location (Eq. 13, Fig. 12).

Conclusion
In this study, a coupled ice flow−mass balance−supraglacial debris cover model was used to simulate the response of the Djankuat Glacier to past, present and future climatic changes between 1752 and 2100 AD. We conducted, for the first time, 595 explicit time-dependent modelling of a Caucasian glacier, including an extended and physically based subroutine related to supraglacial debris cover evolution that was not yet integrated in time-dependent numerical flow line models. As it turns out, the Djankuat Glacier has been retreating almost continuously since the 1850s AD, with some minor steady states or readvances during periods with clusters of colder and/or wetter conditions. The model reconstructed the observed retreat fairly well but required additional mass balance perturbations up to a maximum of +0.5 m yr⁻¹ w.e., which were applied 600 iteratively via dynamic calibration. However, since the start of the calibration period in the balance year 1967/68 AD, no artificial mass balance perturbations were needed, ensuring proper model calibration and credibility.
The future behaviour of the glacier is determined by corresponding changes in air temperature, precipitation and supraglacial debris cover. A temperature increase of 1 °C can only be compensated by a precipitation increase of ca. 35 %, which is not indicated by future climatic projections in the study area. Hence, all scenarios agree to a rapid decline during the following 605 decade, as a response to the accelerating warming since the 1990s AD. Even after considering constant present-day climatic conditions, the glacier will shrink drastically to ca. 30 % of its current length and surface area by 2100 AD, indicating the imbalance between the current glacier geometry and the present climate. However, none of the future scenarios cause a total disappearance by the end of the modelling period. Nevertheless, the glacier will retreat most drastically (ca. -93 % of its current surface area) under the RCP 8.5 scenario, as even the thick ice on the high elevations of the Djantugan Plateau will 610 be affected by significant melting. Although the glacier is close to surpassing its peak water discharge point, the modelled temporal evolution of total runoff volumes indicates that, in particular the melting of ice on these higher parts of the glacier in higher-temperature scenarios, temporarily stabilizes runoff near the end of the modelling period due to the mass balanceelevation feedback.
The presence of a supraglacial debris cover is shown to significantly affect glacier geometry during the modelling period. 615 Hence, the effect of debris-related melt reduction on the eventual glacier length by 2100 AD is dependent upon the trade-off between the growth rate of the total supraglacial debris mass, the efficiency of down-glacier advection of supraglacial debris, the glacier retreat rate, the connectivity between the debris source and the main glacier, and finally the distance between the front and the input location at the time of source release. It turns out that debris-related effects are highest when either debris thickness and area are large, or when melt-reducing effects start earlier in time and/or more down-glacier in space in a 620 relatively warm climate. However, it must be noted that for some of the conducted experiments, the addition of an extra debris source did not (significantly) influence the glacier's geometry. As such, when temperatures increase even further, potential inhibiting effects of too rapid shrinkage are to be considered. Hence, accelerated frontal retreat, disrupted debris discharge and/or connectivity issues at the debris input location may prevent the establishment of a proper melt-reducing effect. 625 Code availability. The model code was written in MATLAB_R2019a. A coupled ice flow-supraglacial debris cover model for the Djankuat Glacier, that was used as the basis for this research, can be found and downloaded from:    (Solomina et al., 2016;WGMS, 2018). Approximate distances and direction to the Djankuat Glacier are indicated. 880 1850 1875 1900 1925 1950 1975 2000 Year (   . We refer to the text and Table 2 for more details. total annual mass balance of the Djankuat Glacier with changing geometry. Observed length variations are derived from lichenometric dating of moraines in the paleovalley, historic documents, field measurements and/or recent satellite imagery (Boyarsky, 1978;Zolotarev, 1998;Petrakov et al., 2012;WGMS, 2018). An additional model run for a 100% clean ice glacier was conducted, which is shown in the box in (a).      (D'Arrigo et al., 2001;Toucham et al., 2003;Griggs et al., 2007;Köse et al., 2011;Martin-Benito et al., 2016) 1752−present (a) Use multi-proxy / multi-model mean approach.
(b) Convert to 3-hourly values by using the observed Terskol data sequence as base but corrected for monthly amounts derived before.
(b) Convert to 3-hourly values by using the observed Terskol data sequence as base but corrected for monthly amounts derived before.