2.1.1. General setting

The Djankuat Glacier (43°12′ N, 42°46′ E) is a partly debris-covered, temperate valley glacier in the Caucasus, near the Russian Federation–Georgia border (Fig. 1). The glacier currently spans an elevation range between roughly 2750 and 3750 m a.s.l. It has been monitored annually since the late 1960s, producing one of the longest and most comprehensive datasets in the world. The dataset includes measurements related to on-glacier automatic weather station (AWS) data, glacier geometry, annual local surface mass balances (SMBs), glacier dynamics and debris cover evolution (Boyarsky, Reference Boyarsky1978; Baume and Popovnin, Reference Baume and Popovnin1994; Popovnin, Reference Popovnin1999; Popovnin and Naruse, Reference Popovnin and Naruse2005; Pastukhov, Reference Pastukhov2011; Lavrentiev and others, Reference Lavrentiev, Kutuzov, Petrakov, Popov and Popovnin2014; Rets and others, Reference Rets2019; Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020; WGMS, 2023; Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024; Popovnin and Gubanov, Reference Popovnin and Gubanov2025). Owing to this record, the Djankuat Glacier has been selected as a WGMS reference glacier for the Caucasus region (WGMS, 2023). Like other glaciers in the area, the glacier has strongly retreated over the last several decades, with a cumulative mass balance of −19.77 m w.e. since 1967/68, and its physical surface area decreasing by −0.57% yr⁻¹ over the same period (Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020; WGMS, 2023; Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024). Due to the presence of a significant debris cover on part of its ablation zone, as well as the availability of abundant observational data, the glacier presents itself as an ideal candidate for debris-covered glacier modelling studies.

Figure 1.

The Djankuat Glacier and its observational network used for surface mass-balance estimation. Shown are the glacier and debris outline on top of a Sentinel-2 L1C True Color satellite image of 5 September 2025, together with the locations of the measurement network for the 2024/25 balance year and the debris-covered ice AWS location in 2009. The inset in the upper right shows a drone picture of the (partly) debris-covered snout taken by Victor V. Popovnin in 2018.

1 Long description

Overview of the study area and field measurements on Djankuat Glacier, Caucasus. The main map shows the glacier extent in 2025 over high-resolution satellite imagery, with 100 m elevation contours. The debris-covered area is indicated by the orange hatched polygon. The locations of snow pits, debris stakes, ice stakes, an automatic weather station (AWS; 2009), and snow-depth measurements collected during the 2025 field campaign are shown. Major topographic features, including Djantugan Peak, the Djantugan firn plateau, Uya-tau Peak, the ice divide, and the neighbouring Bashkara Glacier, are labelled. The inset in the upper left locates Djankuat Glacier within the central Caucasus, while the upper-right inset presents an oblique aerial view of the debris-covered glacier tongue.

2.1.2. Supraglacial debris cover

The Djankuat Glacier receives near-continuous debris input from steep lateral moraines, exposed bedrock and headwall erosion, which is further enhanced by glacier retreat, the high relief surrounding the glacier and its moderate rock erodibility, as well as the high precipitation regime around the Main Caucasus Ridge (Popovnin and others, Reference Popovnin, Rezepkin and Tielidze2015; Popovnin and Gubanov, Reference Popovnin and Gubanov2025). Major rockfalls in the early 2000s further contributed to the debris supply (Stokes and others, Reference Stokes, Popovnin, Aleynikov, Gurney and Shahgedanova2007). Long-term measurements have revealed that the supraglacial debris cover has expanded rapidly, from ∼2% of the glacier area in 1968 (0.065 km²) to ∼20% in 2022 (0.439 km²). Detailed constructed debris thickness maps have revealed that the glacier-averaged debris thickness increased from 0.28 m in 1983 to 0.60 m in 2022 (Popovnin and others, Reference Popovnin, Rezepkin and Tielidze2015; Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024; Popovnin and Gubanov, Reference Popovnin and Gubanov2025). These debris thickness measurements were made by either manual excavation or by penetrating the debris with a steel probe until the ice was reached. The debris currently clearly affects the glacier geometry: the heavily debris-covered western part of the snout retreated significantly less than the cleaner eastern part, and thickly buried ice bodies extend well below the current terminus of the active glacier. The relatively thick debris cover (up to more than 1 m in some areas) slows the frontal retreat of the glacier, produces (near-)stagnant dead ice bodies and affects the (glacio)morphological and hydrological regime of the glacier (Popovnin and others, Reference Popovnin, Rezepkin and Tielidze2015; Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024; Verhaegen and others, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024; Popovnin and Gubanov, Reference Popovnin and Gubanov2025). Based on the extrapolation of these trends and the expected future climate warming (i.e. potentially resulting in more englacial melt out, reduced flow velocities and off-glacier debris discharge, and an increased bedrock exposure and slope instability), the debris supply onto the glacier is believed to increase even further in the future, as long as it does not retreat beyond the debris sources (Popovnin and others, Reference Popovnin, Rezepkin and Tielidze2015; Popovnin and Gubanov, Reference Popovnin and Gubanov2023; Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024). Site-specific debris-related properties of the granite/gneiss-type rocks on Djankuat were investigated by Bozhinskiy and others (Reference Bozhinskiy, Krass and Popovnin1986) and have been used, together with meteorological data from an AWS on top of debris-covered ice, by Verhaegen and others (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024) in a detailed energy-balance and mass-balance model for 1 specific balance year.

2.1.3. Meteorological data and mass-balance measurements

The local annual SMBs on the Djankuat Glacier are determined from direct stake measurements, repeated snow pit observations and winter snow depth surveys all over the accessible/safe glacier area. The observational network typically includes ca. 30 stakes, with their quantity varying from year to year (Fig. 1, Table 1). In addition, 2–4 reference snow pits are excavated annually and monitored until the seasonal snow cover has completely melted (if so). Usually, only 0–1 of these pits persists through the summer, in which case they can also provide a direct local estimate of the positive annual mass balance by using the measured snow density values (Popovnin and others, Reference Popovnin, Gubanov, Lisak and Toropov2024). Snow depth surveys, based on direct manual measurements of the seasonal snow cover by metallic rods, consist of 250–450 measurement points, varying from year to year. Such abundant measurements make the Djankuat Glacier one of the few Caucasian glaciers with a multi-decadal, high-quality mass-balance record. In this study, we use the direct stake readings from both the clean ice and debris-covered ice areas during the 2009–25 period for tuning of the energy/mass-balance models. Over the course of those balance years, this resulted in more than 500 stake readings, of which ∼12.5% were situated in the (partly) debris-covered area of the glacier (Table 1).

Table 1.

Stake network and observed glacier-wide mass balance of the Djankuat Glacier.

Meteorological data (air temperature and precipitation) are available from the Terskol weather station (2141 m a.s.l., at 16 km NW of the glacier) at 3 hourly temporal resolution (Fig. 2). Data from this station have been proven to show a high correlation with the meteorological regime on the glacier itself (Rybak and Rybak, Reference Rybak and Rybak2018; Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020) and hence serve as an ideal forcing dataset for the mass-balance model. We transfer the data from the Terskol weather station to any given surface elevation on the Djankuat Glacier by using elevation-dependent temperature and precipitation gradients (Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020). Insolation data at the top of the atmosphere for the latitude of the Djankuat Glacier can be further calculated from basic astronomical formulas.

Figure 2.

The climatic data used to force the energy- and mass-balance model in this study during the 2009–25 period. Subplots represent (a) daily minimum, maximum and average air temperature ( ${T_{\text{a}}}$), (b) daily precipitation sums ( $P$) and (c) theoretical average and maximum daily incoming solar insolation ( ${S_ \downarrow }$) at the top of the atmosphere (TOA) calculated from basic astronomical formulas. Temperature and precipitation data are representative for the Terskol weather station location (2141 m a.s.l.) but can be extrapolated to any surface elevation on the glacier using elevation-dependent temperature and precipitation gradients (Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020). The solar radiation data are further influenced by atmospheric transmissivity, surface albedo, topographic shading and surface geometry. Note that in the actual model, 3 hourly data are used instead of daily values.

Figure 2 Long description

ime series of the meteorological forcing used to drive the glacier model from the Terskol weather station for the period 2009–2025. (a) Daily mean (black), minimum (blue), and maximum (red) air temperatures, showing a pronounced seasonal cycle with winter minima below −20 °C and summer maxima exceeding 20 °C. (b) Daily precipitation totals, illustrating the highly variable occurrence and intensity of precipitation events throughout the record, with occasional daily totals approaching 100 mm. (c) Daily top-of-atmosphere (TOA) insolation, showing the seasonal variation in daily mean (orange) and maximum (red) incoming solar radiation controlled by the annual cycle in solar geometry.

Observed surface temperature data over debris are available from an AWS that was present on the debris-covered area of the glacier during the 2009 summer season (Rets and others, Reference Rets2019). The AWS (43.20° N, 42.76° E) was placed at an elevation of 3050 m a.s.l. on top of 65 cm of supraglacial debris and started fully operating on 1 July 2009. It recorded hourly values of air relative humidity, air temperature, shortwave and longwave radiation, wind speed and direction, and atmospheric pressure at 2 m above the surface until the end of September. The measured outgoing longwave radiation was used to calculate debris surface temperatures using the Stefan–Boltzmann law and a typical value of the debris surface emissivity of 0.95 (Reid and Brock, Reference Reid and Brock2014; Fontrodona-Bach and others, Reference Fontrodona-Bach2025). Since these data are only available for one summer season, we nevertheless select the data from the Terskol weather station rather than the on-glacier site as forcing data for our model.

2.2.1. The surface energy balance

The model used in this study consists of surface energy balance (SEB), vertical heat conduction (HC) and surface mass balance (SMB) components. The full energy-balance equation at a given surface can be written as:

(1)\begin{align}{\psi _{\text{s}}} & = {S_ \downarrow }\left( {1 - \alpha } \right)\tau + {\varepsilon _{\text{a}}}\sigma T_{\text{a}}^4 - {\varepsilon _{\text{s}}}\sigma T_{\text{s}}^4 + {c_{\text{a}}}{\rho _{\text{a}}}{C_{\text{E}}}u\left( {{T_{\text{a}}} - {T_{\text{s}}}} \right) \nonumber \\ &\quad + {L_{\text{v}}}{\rho _{\text{a}}}{C_{\text{E}}}u\left( {{q_{\text{a}}} - {q_{\text{s}}}} \right) + {\rho _{\text{w}}}{c_{\text{w}}}P\left( {{T_{\text{a}}} - {T_{\text{s}}}} \right) + k\frac{{\partial T}}{{\partial z}} = 0.\end{align}

Here, ${S_ \downarrow }$ is the downward solar radiation, $\alpha $ is the surface albedo, $\tau $ is the sky transmissivity, $\varepsilon $ is the emissivity, $\sigma $ is the Stefan–Boltzmann constant, $c$ is the specific heat capacity, $\rho $ is the density, $T$ is the temperature, $u$ is the wind speed, ${C_{\text{E}}}$ is the turbulent exchange coefficient, $P$ is the precipitation rate, ${L_{\text{v}}}$ is the latent heat of evaporation, $q$ is the specific humidity (which can be derived from the relative humidity RH_a, temperature T_a and air pressure p) and $k$ is the thermal conductivity. Subscripts ‘a’, ‘s’ and ‘w’ stand for air, surface and water, respectively. The first term on the right-hand side denotes the net shortwave radiation ( ${Q_{\text{S}}}$), the second and third terms, the net longwave radiation ( ${Q_{\text{L}}}$), the fourth and fifth terms, the sensible ( ${Q_{{\text{SH}}}}$) and latent heat fluxes ( ${Q_{{\text{LH}}}}$), respectively. The sixth term represents the heat added due to rainfall ( ${Q_{\text{R}}}$), and the last term is the surface conductive heat flux ( ${Q_{\text{C}}}$). Fluxes are considered positive when directed towards the surface. Below, we denote the temperature-dependent fluxes (the sum of the net longwave, rain, sensible and latent heat fluxes) with symbol ${Q_{\Delta T}}$. The full SEB equation as depicted in Eqn (1) was used in the sophisticated mass-balance model of Verhaegen et al. (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024) for the Djankuat Glacier during 1 balance year. In practice, however, several terms in Eqn (1) are difficult to constrain accurately. This motivates the development of simplified parameterizations that reduce the number of required inputs while retaining the dominant physical controls on sub-debris melt.

2.2.2 A relationship for the ‘temperature-dependent’ fluxes

Due to the presence of hard-to-constrain parameters in the full energy-balance equation, a relationship between the air temperature ${T_{\text{a}}}$ and the temperature-dependent fluxes ${Q_{\Delta T}}$ for snow and clean ice was already established by Oerlemans (Reference Oerlemans2001) and applied to a large set of glaciers worldwide by Giesen and Oerlemans (Reference Giesen and Oerlemans2012), as mentioned earlier. They found an approximately constant value of the temperature-dependent fluxes for air temperatures below a certain threshold, the latter often taken to be the melting point. This constant value was denoted ${c_0}$. For temperatures higher than the threshold, a linear increase of the temperature-dependent flux with air temperature was noted. In that case, ${T_a}$ is scaled with a factor ${c_1}$. The net energy flux at a snow/ice surface ${\psi _{\text{s}}}$ (and hence the energy available for melting of ice, below denoted with subscript ‘i’, and snow, with subscript ‘s’) is thus parameterized as (Oerlemans, Reference Oerlemans2001; Giesen and Oerlemans, Reference Giesen and Oerlemans2012):

(2)\begin{equation}{Q_{{M_i}}} = {Q_{{M_s}}} = {{{\psi }}_s} = \left\{ \begin{array}{*{20}{c}} {{S_ \downarrow }\left( {1 - \alpha } \right)\tau + {c_0} if {T_a} \lt 0^\circ C} \\ {{S_ \downarrow }\left( {1 - \alpha } \right)\tau + {c_0} + {c_1}{T_a} if {T_a} \geqslant 0^\circ C} \end{array}\right.\end{equation}

For snow and clean ice, the surface conductive heat flux is assumed to be 0. The parameters ${c_0}$ and ${c_1}$ represent the scaling factors for the ‘temperature-dependent fluxes’, which are linearized around the melting point. In the model, these parameters are used for tuning the clean ice mass-balance model (see Section 2.2.5). We also check the applicability of this relationship for the Djankuat Glacier specifically, using the detailed 2D energy and mass-balance model of Verhaegen and others (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024), see Section 3.1.

Apart from a snow and ice surface, a relationship for the ‘temperature-dependent’ fluxes ${Q_{\Delta T}}$ is also apparent for a debris-covered ice surface. In that case, a scaling with ${T_{\text{a}}}$ alone, however, does not seem to work well. Instead, a relationship between ${T_a}$ −  ${T_{\text{s}}}$ and the temperature-dependent fluxes is to be expected here, as the temperature-dependent fluxes all respond to changes in the air-surface temperature difference (McCarthy, Reference McCarthy2025; Pellicciotti and others, Reference Pellicciotti2026). Interestingly, such a relationship would create an attractive option for a rather simple but physically motivated parameterization of energy- (and hence mass-) balance models for debris-covered glaciers, where it can help to counteract the large amount of (poorly constrained) parameters in the energy-balance equation. Accordingly, we adopt an Oerlemans-type formulation of the SEB over debris, with an intercept ( ${d_0}$) and a slope ( ${d_1}$):

(3)\begin{equation}{\psi _{\text{s}}} = {S_ \downarrow }\left( {1 - \alpha } \right)\tau + {d_0} + {d_1}({T_{\text{a}}} - {T_{\text{s}}}) + {k_{\text{d}}}\frac{{\partial {T_{\text{d}}}}}{{\partial z}} = 0.\end{equation}

In this case, the subscript ‘d’ refers to the debris pack. The constants ${d_0}$ and ${d_1}$ quantify the temperature-dependent fluxes ${Q_{\Delta T}}$ (the net longwave flux, the sensible and latent heat fluxes, and the heat flux added by rain), similarly to the parameterization for snow and clean ice (Fig. 3). This approach therefore parameterizes the explicit representation of sky/surface emissivity, wind speed, humidity gradients and heat added by rain. In other words, we have replaced these hard-to-constrain inputs with simplified and tunable parameters, rather than having removed those parameters. Again, we check whether this relationship applies for the Djankuat Glacier specifically, using the detailed 2D energy- and mass-balance model of Verhaegen and others (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024) (Section 3.1).

Figure 3.

Simplified representation of the heat conduction and surface energy- and mass-balance model used in this study.

Figure 3 Long description

illustration of the conceptual framework of the debris-covered glacier model. A debris layer overlies glacier ice and may be covered by a transient snow layer. The debris surface exchanges energy with the atmosphere through shortwave and longwave radiation, sensible and latent heat fluxes, and rainfall heat input. The resulting surface energy flux drives heat conduction through the debris, whose thermal response depends on its thermal conductivity, density, and heat capacity. The transient temperature profile within the debris determines the conductive heat flux reaching the debris–ice interface, where it is used to melt the underlying ice.

2.2.3. SMB model

The SMB model, to which the above-mentioned parameterizations will be applied, is largely based upon the models presented in Verhaegen and others (Reference Verhaegen, Huybrechts, Rybak and Popovnin2020, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024). It has a spatial resolution of 25 by 25 m and is based upon the difference between accumulation (ACC) and runoff (RO) over the balance year, so that the local SMB is defined as:

(4)\begin{equation}{b_{\text{s}}} = \smallint\limits_{{\text{yr}}} \left( {{\text{ACC}} - {\text{RO}}} \right){\text{d}}t.\end{equation}

The balance year is here considered between the period of 1 October and 30 September. The model is forced by 3 hourly insolation and air temperature/precipitation data from the Terskol weather station.

2.2.3.2. RO from snow and clean ice

For snow and clean ice, the RO is calculated as:

(5)\begin{equation}RO = \left\{ \begin{array}{*{20}{c}} {{W_s} if {h_s} \gt 0} \\ {{M_i} if {h_s} = 0 \& {h_d} = 0} \end{array}\right. \end{equation}

In the case of snow on the surface ( ${h_{\text{s}}} \gt 0$), RO is calculated as the meltwater outflow from a saturated snowpack ${W_{\text{s}}}$, according to Schaefli and Huss (Reference Schaefli and Huss2011). For snow-free conditions ( ${h_{\text{s}}} = 0$) over clean ice ( ${h_{\text{d}}} = 0$), RO is equal to the ice melt ${M_{\text{i}}}$, which is calculated as:

(6)\begin{equation}{M_{\text{i}}} = {\text{max}}\left( {0,\frac{{{Q_{{{\text{M}}_{\text{i}}}}}{{\Delta }}t}}{{{\rho _{\text{w}}}{L_{\text{m}}}}}} \right),\end{equation}

where ${\rho _{\text{w}}}$ is the density of water, ${L_{\text{m}}}$ the latent heat of fusion, ${Q_{{{\text{M}}_{\text{i}}}}}$ the energy available for ice melt and ${{\Delta }}t$ the model time step. The parameters ${c_0}$ and ${c_1}$ in the energy-balance equation (Eqn (2)) are used for tuning the clean ice mass-balance model (see Section 2.2.5), while all other parameters are taken from the models of Verhaegen and others (Reference Verhaegen, Huybrechts, Rybak and Popovnin2020, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024).

2.2.3.3. RO from debris-covered ice

For debris-covered ice, we use the simplified formulation of the SEB discussed in Section 2.2.2 (Eqn (3)), where the temperature-dependent fluxes ${Q_{\Delta T}}$ have been parameterized. The iterative numerical Newton–Raphson method is then used to calculate surface temperatures from the SEB equation, which are used as an upper boundary condition for the vertical heat-transfer model. The air temperature is used as a first guess for the initial time step of the balance year and for the first time step after complete snow cover disappearance. In those cases, a simple linear vertical temperature gradient within the debris is assumed. For all other time steps, the initial guess for the surface temperature is defined as the surface temperature from the previous time step, while the change in the internal energy within the debris pack is then defined by the thermodynamic heat equation. For simplicity, our model will only incorporate the heat-conduction term and thus discards any other nonconductive processes that may contribute to heat transfer from the surface towards the ice (Fig. 3):

(7)\begin{equation}{\rho _{\text{d}}}{c_{\text{d}}}\frac{{\partial {T_{\text{d}}}}}{{\partial t}} = \frac{\partial }{{\partial z}}\left( {{k_{\text{d}}}\frac{{\partial {T_{\text{d}}}}}{{\partial z}}} \right)\end{equation}

with ρ is the density (kg m⁻³) and c the heat capacity (J kg⁻¹ K⁻¹) and where $\Delta z$ is represented by a fixed number of ${n_{\text{d}}}$ = 10 layers (i.e. $\Delta z$ =  ${h_{\text{d}}}/{n_{\text{d}}}$). The subscript ‘d’ refers to values that are applicable to the debris pack and can be thought of as a combination of the thermal properties of the rock (subscript ‘r’) and the contents of the pore spaces (subscript ‘a’), i.e. ${k_d} \approx {k_r}\left( {1 - {\phi _d}} \right) + {k_a}{\phi _d}$ and ${\rho _d}{c_d} \approx {\rho _r}{c_r}\left( {1 - {\phi _d}} \right) + {\rho _a}{c_a}{\phi _d}$, with ${\phi _{\text{d}}}$ the debris cover porosity. The pore spaces can be filled with air, water or ice, and the value of the debris cover porosity can vary with depth. In our model, however, we parameterize all these processes in a single depth-independent value for ${\rho _{\text{d}}}{c_{\text{d}}}$ and ${k_{\text{d}}}$, but we perform some model sensitivity experiments to evaluate how including a changing value for ${\phi _d}$ affects the model outcome.

To solve for the internal debris temperatures ${T_{\text{d}}}$, we use a numerical Crank–Nicholson scheme (Reid and Brock, Reference Reid and Brock2010; Rounce and others, Reference Rounce, King, McCarthy, Shean and Salerno2018). In essence, this equation propagates the surface perturbation from the SEB equation downward as a damped and delayed thermal wave, resulting in nonlinear vertical temperature profiles. We furthermore assume that the ice is at its melting point at the ice–debris interface, maintaining a constant temperature of 0°C. The sub-debris ice melt ${M_{\text{d}}}$ is calculated using the downward conductive flux at the vertical ice–debris interface $Q_{\text{C}}^ \downarrow $, which itself depends on the temperature gradient ( $\partial {T_{\text{d}}}$/ $\partial z$) and the effective thermal conductivity ( ${k_{\text{d}}}$) at the base of the debris layer (Reid and Brock, Reference Reid and Brock2010; Verhaegen and others, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024):

(8)\begin{equation}{Q_{{{\text{M}}_{\text{d}}}}} = Q_{\text{C}}^ \downarrow = {\left( {{k_{\text{d}}}\frac{{\partial {T_{\text{d}}}}}{{\partial z}}} \right)_{\text{b}}},\end{equation}

where subscript ‘b’ implies the base of the debris layer. The sub-debris melt rate is then eventually calculated as (with ${Q_{{{\text{M}}_{\text{d}}}}}$ the energy available for sub-debris melt):

(9)\begin{equation}{M_{\text{d}}} = {\text{max}}\left( {0,\frac{{{Q_{{{\text{M}}_{\text{d}}}}}{{\Delta }}t}}{{{\rho _{\text{w}}}{L_{\text{m}}}}}} \right)\end{equation}

and the eventual RO is similar to that of the clean ice mass-balance model:

(10)\begin{equation}RO = \left\{ \begin{array}{*{20}{c}} {{W_s} if {h_s} \gt 0 } \\ {{M_d} if {h_s} = 0 \& {h_d} \gt 0 } \end{array}\right. \end{equation}

The parameters ${d_0}$, ${d_1}$ and ${k_{\text{d}}}$ in the energy-balance equation (Eqn (3)) were used for tuning the debris-covered ice mass-balance model (see Section 2.2.5). Note that we assume fully debris-covered pixels (i.e. ${A_{\text{d}}}$/ $A$ = 1) for all debris thicknesses in the model unless indicated otherwise (see, e.g., Appendix A).

With ${M_i}$ and ${M_d}$ known, the Østrem curve can be constructed by calculating relative melt modification factors ${f_{\text{d}}}$ =  ${M_d}$/ ${M_i}$ and plot them with respect to the corresponding debris thickness ${h_{\text{d}}}$. Within the Østrem curve, we then also define some notable points to quantify its shape (Popovnin and Rozova, Reference Popovnin and Rozova2002; Kirkbride and Dugmore, Reference Kirkbride and Dugmore2003; Verhaegen and others, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024). These are the effective debris thicknesses (i.e. the debris thickness at which relative melt is maximally enhanced), the critical debris thicknesses (i.e. the debris thickness at which relative melt is equal to that of clean ice) and the characteristic debris thicknesses (i.e. the debris thickness at which relative melt is a fraction ${e^{ - 1}}$ or ∼37% of that of clean ice).

2.2.4. Model configuration and experiment design

The primary model calibration uses the Oerlemans-type surface energy balance (SEB) approach, coupled to a one-dimensional heat-conduction (HC) model that includes the temporal internal heat storage (IHS) term within the debris pack, as discussed above. Afterwards, we conduct some sensitivity experiments by varying key parameters in the model. However, in addition to this, we also perform a series of experiments with different levels of model complexity to evaluate the robustness of our sub-debris melt model approach (Table 2). Hence, we verify the results against two alternative model configurations. First, we compare the simplified SEB parameterization with a sophisticated full energy-balance model, similar to the complexity of the Verhaegen et al. (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024) model. For this experiment, summer season incoming shortwave radiation ( ${S_ \downarrow }\tau $), near-surface wind speed ( $u$), surface pressure ( $p$), sky emissivity ( ${\varepsilon _{\text{a}}}$) and relative air humidity ( ${\text{R}}{{\text{H}}_{\text{a}}}$) for the 2009–25 period were obtained from ERA5-Land (Muñoz-Sabater and others, Reference Muñoz-Sabater2021) for the grid cell closest to Djankuat Glacier. These were subsequently bias-corrected using the 2009 AWS observations, following the methodology of Verhaegen and others (Reference Verhaegen, Rybak, Popovnin and Huybrechts2024). Second, we perform an experiment that excludes the temporal IHS term within the debris, thereby assuming only linear vertical temperature gradients within the debris layer. These additional experiments should provide insight into the sensitivity of the modelled sub-debris melt to the degree of physical complexity included in the debris energy-balance and heat-conduction formulation and help assess the robustness of the simplified SEB approach.

Table 2.

Main model configurations with varying degrees of complexity in this study. The abbreviations ‘SEB’ and ‘IHS’ stand for ‘surface energy balance’ and ‘internal heat storage’ respectively. See main text for the meaning of the symbols.

Table 2 Long description

Three model configurations for debris-covered glacier melt are compared by their surface energy balance approach, internal heat storage treatment, required atmospheric forcing, and debris input variables. The simplest option uses a parameterized surface energy balance and a simplified heat storage assumption based on linear vertical temperature gradients, forced only by air temperature, precipitation, and incoming shortwave radiation, and it needs a small set of debris properties including albedo and thermal conductivity. The main study configuration keeps the parameterized surface energy balance but adds full internal heat storage within the debris, using the same limited atmospheric forcing while requiring additional debris properties such as density and heat capacity. The most complex option computes the full debris surface energy balance and full internal heat storage, optionally allowing a vertical porosity gradient, and it requires the broadest atmospheric forcing set, adding wind, air emissivity, humidity, and air pressure to the basic inputs. It also expands debris inputs to include surface emissivity, an exchange coefficient, a surface moisture or humidity-related term, and the debris cover porosity. Overall, complexity increases from by adding more physical processes and therefore more forcing data and debris parameters. Variable names are listed as categories of inputs rather than defining measurement units or data sources, so comparability depends on consistent parameter definitions across implementations.

›In Appendix A, we further explore a potential and optional model extension in which the relative melt enhancement for thin/patchy debris is included in the model by parameterizing debris cover patchiness.

2.2.5. Model calibration

(a) Parameter values and input data

The Digital Elevation Model (DEM) from Morozova and Rybak (Reference Morozova and Rybak2017) is used to calculate the insolation based on the local slope and surface aspect using basic astronomical formulas, as well as to incorporate the effects of topographic shadow and self-shading (Verhaegen and others, Reference Verhaegen, Rybak, Popovnin and Huybrechts2024). We also use the DEM to spatialize the air temperature and precipitation fields, using elevation-dependent gradients (Verhaegen and others, Reference Verhaegen, Huybrechts, Rybak and Popovnin2020). For the atmospheric transmissivity and debris surface albedo, an AWS-derived value of 0.58 and 0.10 is assigned, respectively.

For the volumetric heat capacity ${\rho _{\text{d}}}{c_{\text{d}}}$ in Eqn (7), values are taken over from Bozhinskiy and others (Reference Bozhinskiy, Krass and Popovnin1986) for the Djankuat Glacier. However, in other instances, standard values can be used depending on the rock lithology. For the gneiss/granite-type material that is found on the Djankuat Glacier, for example, typical values for the rock density and specific heat capacity are 2650 kg m⁻³ and 800 J kg⁻¹ °C⁻¹, respectively (Robertson, Reference Robertson1988). A value of 0.3 is also often assumed to be a good representative for the overall porosity of a debris layer (Nicholson and Benn, Reference Nicholson and Benn2013; Fontrodona-Bach and others, Reference Fontrodona-Bach2025). For our model, this implies an effective volumetric heat capacity of $1.5 \times {10^6}$ J m⁻³ °C⁻¹.

(b) Snow and clean ice model

Model calibration of the multiyear energy/mass-balance model over the 2009–25 period starts with the tuning of the snow and clean ice mass-balance model. For that, we calibrate the parameters ${c_0}$ and ${c_1}$ in Eqn (2) until the root-mean-square error (RMSE) between the observed and modelled local annual SMBs from the stake measurements in the clean ice area is minimized (Fig. 4):

(11)\begin{equation}{\text{RMS}}{{\text{E}}_{\text{M}}} = \sqrt {\frac{1}{{{n_{\text{M}}}}}\sum\limits_{i = 1}^{{n_{\text{M}}}} {{\left( {{b_{{\text{s}},{\text{obs}}}} - {b_{{\text{s}},{\text{mod}}}}} \right)}^2}} ,\end{equation}

where ${n_{\text{M}}}$ is the number of local mass-balance observations in the clean ice area. We vary ${c_0}$ in steps of 5 W m⁻² between −110 and 0 W m⁻² and ${c_0}$ in steps of 1 W m⁻² °C⁻¹ between 2 and 25 W m⁻² C⁻¹.

Figure 4.

Schematic visualization of the calibration procedure for the simplified surface energy-balance approach (‘simple SEB + full IHS’) in this study.

Figure 4 Long description

The flowchart illustrates the calibration process for snow/clean ice and debris-covered ice surface mass balance (SMB) models. It begins with varying parameters c_0 and c_1. The next step is to find the minimized root mean square error (RMSE) between modeled and observed b_s for snow/clean ice. Then, c_0 is used as an initial guess for d_0. The process continues by varying d_1 and k_. The next step is to find the minimized composite RMSE_CN between modeled and observed b_s and T_s for debris-covered ice. If the minimum RMSE _CN for d_1 and k_d has not yet reached its lowest value, the process loops back to adjust d_0. Finally, the process concludes with determining the final values for d_0, d_1 and k_d.

(c) Debris-covered ice model

To tune the debris-covered ice model, we follow a slightly different and more sophisticated approach. Here, we use three parameters: ${d_0}$, ${d_1}$ and ${k_{\text{d}}}$ (Eqns (3) and (7)). However, to reduce the risk of parameter equifinality, we can restrict the parameter range of ${d_0}$ based on the best-fit value of ${c_0}$. Indeed, the value of ${c_0}$ (the intercept for snow/clean ice) applies for a value of ${T_a}$ = 0°C and ${d_0}$ (the intercept for debris-covered ice) for ${T_a}$ −  ${T_{\text{s}}}$ = 0°C. If we assume a melting snow/ice surface (i.e. a surface temperature fixed at 0°C), the sensible heat flux is absent in both cases, and any difference in the net longwave radiation is only caused by a different surface emissivity of a snow/ice and debris surface (which, in fact, are relatively similar). For mid-latitude glaciers, the latent heat flux is furthermore of minor importance when compared to other energy fluxes most of the time. This should imply that ${c_0} \approx {d_0}$ and that the best-fit value for ${c_0}$ can thus be used as a first guess for ${d_0}$ in the tuning procedure for the debris-covered ice mass-balance model. This leaves two remaining parameters to be varied for a single value of ${d_0}$: ${d_1}$ and ${k_d}$. In a subsequent step, ${d_0}$ can be slightly varied to check whether this further improves the agreement between observed and modelled sub-debris melt rates and debris surface temperatures (Fig. 4).

We thus start the tuning by setting ${d_0}$ to the calibrated value of ${c_0}$ and vary ${d_1}$ in steps of 2 W m⁻² between 4 and 60 W m⁻²  $^\circ $C⁻¹ and ${k_d}$ in steps of 0.2 W m⁻¹  $^\circ $C⁻¹ between 0.4 and 3.4 W m⁻¹  $^\circ $C⁻¹. The latter is a sufficiently large interval to capture commonly reported values in the literature for ${k_d}$ (Rowan and others, Reference Rowan2021; Fontrodona-Bach and others, Reference Fontrodona-Bach2025). To calibrate the sub-debris melt model, we minimize a multi-objective ‘combined normalized RMSE’ (hereafter denoted $RMS{E_{CN}}$) using two observational datasets: (i) sub-debris melt measurements from stakes in the debris-covered area (2009–25), and (ii) debris surface temperatures measured by the AWS during the 2009 summer season:

(12)\begin{equation}{\text{RMS}}{{\text{E}}_{{\text{CN}}}} = \sqrt {\begin{array}{l}{w_{\text{M}}}\left[ {\frac{1}{{{n_{\text{M}}}}} \sum \limits_{i = 1}^{{n_{\text{M}}}} {{\left( {\frac{{{b_{{\text{s}},{\text{obs}}}} - {b_{{\text{s}},{\text{mod}}}}}}{{{\sigma _{{b_{{\text{s}},{\text{obs}}}}}}}}} \right)}^2}} \right]\\ + {w_T}\left[ {\frac{1}{{{n_T}}}\sum \limits_{i = 1}^{{n_T}} {{\left( {\frac{{{T_{{\text{s}},{\text{obs}}}} - {T_{{\text{s}},{\text{mod}}}}}}{{{\sigma _{{T_{{\text{s}},{\text{obs}}}}}}}}} \right)}^2}} \right]\end{array}} ,\end{equation}

where the individual terms below the square root represent ${\text{RMS}}{{\text{E}}_{\text{M}}}$ and ${\text{RMS}}{{\text{E}}_{\text{T}}}$, which are the RMSE values of local mass balance and surface temperature observations, respectively, over debris-covered terrain. To make the total error dimensionless, the data are, however, normalized by the standard deviation $\sigma $ of the corresponding observations. Moreover, because the spread in surface temperature errors differs from that of sub-debris melt, we use weights such that both objectives affect the ${\text{RMS}}{{\text{E}}_{{\text{CN}}}}$ similarly:

(13)\begin{align}{w_M} & = \frac{{{{\overline {RMSE} }_{T\left\{ {{d_1},{k_d}} \right\}}}}}{{{{\overline {RMSE} }_{T\left\{ {{d_1},{k_d}} \right\}}} + {{\overline {RMSE} }_{M\left\{ {{d_1},{k_d}} \right\}}}}}, \nonumber\\ {w_T} & = \frac{{{{\overline {RMSE} }_{M\left\{ {{d_1},{k_d}} \right\}}}}}{{{{\overline {RMSE} }_{T\left\{ {{d_1},{k_d}} \right\}}} + {{\overline {RMSE} }_{M\left\{ {{d_1},{k_d}} \right\}}}}}\end{align}

The weights are defined by the average RMSE for the surface temperature and sub-debris melt rates over the entire parameter space of ${d_1}$ −  ${k_{\text{d}}}$ combinations. Although alternatives to our approach exist (e.g. tuning for the surface temperature and sub-debris melt rates separately) that may produce similar best-fit values (McCarthy, Reference McCarthy2025; Pellicciotti and others, Reference Pellicciotti2026), these steps ensure that both variables contribute equally to the composite objective, independent of differences in units or absolute magnitude of the errors. We note that $RMS{E_{CN}}$ represents a composite objective function rather than a classical RMSE, as it combines two datasets with different units and error structures (Deb, Reference Deb2001). In this sense, the formulation should be interpreted within a multi-objective optimization framework, where the weights determine the relative importance of sub-debris melt and surface temperatures. After each iteration of varying ${d_1}$ and ${k_d}$, we changed the offset ${d_0}$ with 5 W m⁻² to check whether the minimum $RMS{E_{CN}}$ across the parameter space further decreases. The calibration procedure stops whenever a minimum for min ( ${\text{RMS}}{{\text{E}}_{{\text{CN}}}}\left\{ {{d_1},{k_{\text{d}}}} \right\}$) as a function of ${d_0}$ is found (Figs 4 and 7d).

3.3.1. Sub-debris melt sensitivity to parameter variations

To evaluate how uncertainties in debris characteristics influence modelled sub-debris melt, we varied six key debris parameters (the effective thermal conductivity ${k_d}$, effective volumetric heat capacity ${\rho _{\text{d}}}{c_{\text{d}}}$, debris albedo ${\alpha _{\text{d}}}$, debris cover porosity ${\phi _{\text{d}}}$ and the SEB tuning parameters ${d_0}$ and ${d_1}$) and examined the resulting changes in the Østrem curves (melt-modification factor versus debris thickness). To have a uniform evaluation, each parameter was varied within an interval of −50 and +50% of its original value (Fig. 9). This also ensured consistency with most reported parameter ranges within the DebDaB database (Fontrodona-Bach and others, Reference Fontrodona-Bach2025).

Figure 9.

Sensitivity of the fitted Østrem curve to variations of debris-related key variables, including the effective thermal conductivity, effective volumetric heat capacity, debris albedo, debris cover porosity and the tuning parameters of the SEB equation ${d_0}$ and ${{\text{d}}_1}$. The inset shows the corresponding effect on the characteristic debris thickness (i.e. the thickness for which sub-debris melt is a fraction ${e^{ - 1}}$ or ∼37% of that of clean ice). The colours represent changes from −50% (red) to +50% (blue) compared to the value used in the best-fit model.

Figure 9 Long description

Sensitivity of the melt modification factor to ±50% perturbations in the principal debris and model parameters. Each panel shows the melt modification factor as a function of debris thickness (0–2 m), with red and blue curves representing negative and positive parameter perturbations, respectively. The analysed parameters are: (a) effective thermal conductivity, k_d (W m^-1 °C^-1), (b) effective volumetric heat capacity, ρ_d*c_d (J m^-3 °C^-1), (c) debris albedo, α_d (-), (d) debris porosity, ϕ_d (-), (e) heat exchange intercept parameter, d_0 (W m^-2), and (f) heat exchange slope parameter, d_1 (W m^-2 °C^-1). Insets show the corresponding characteristic debris thickness, hc (cm), at which debris melt equals clean-ice melt, as a function of the relative parameter perturbation (Δ parameter, %). Thermal conductivity produces the largest changes in both the melt modification factor and the characteristic debris thickness, whereas the remaining parameters primarily affect the thin-debris regime, with progressively smaller sensitivities for volumetric heat capacity, albedo, porosity, d_0, and d_1.

Overall, the sensitivity experiments reveal that most parameters primarily affect the Østrem curve in the thin to moderately thick debris regime. Hence, for debris thicknesses larger than ∼1.5 m, parameter sensitivity becomes increasingly negligible across all experiments and the Østrem curves converge towards a common insulating asymptote. Among all parameters, thermal conductivity ${k_d}$ has the strongest influence. Higher ${k_d}$ values are modelled to support higher melt rates under the debris by enhancing conductive heat transfer towards the ice, while lower ${k_d}$ values suppress melt more rapidly, producing earlier transitions to insulating behaviour. Differences gradually diminish for debris thicker than ∼1.5 m, where all values converge towards strong melt reduction. Also the debris cover porosity ${\phi _{\text{d}}}$ exerts a more noticeable influence by altering the effective thermal properties. Higher porosity values increase the insulation effect and reduce melt under the debris, but again the effect gradually fades for an increasingly thicker cover. The tuning parameters ${d_0}$ and ${d_1}$ affect the Østrem curve to a moderate degree and have less impact than the thermal conductivity. This is to be expected since the surface boundary condition reaches the sub-debris ice only as a dampened signal, especially for thicker debris. Changes in volumetric heat capacity ${\rho _{\text{d}}}{c_{\text{d}}}$ and debris albedo ${\alpha _{\text{d}}}$ affect the curve to a lesser degree and appear second-order or even negligible controls. Taken together, the sensitivity analysis shows that model performance is most responsive to variations in the effective thermal conductivity ${k_d}$ (including the effects of the porosity), moderately responsive to the temperature-dependent tuning parameters ${d_0}$ and ${d_1}$ and less responsive to surface albedo and debris pack volumetric heat capacity variations.

3.3.2. Sub-debris melt sensitivity to model complexity

As stated before, we also compare our model approach to two alternative configurations (Table 2): (i) one including the full SEB formulation (‘full SEB + full IHS’) and (ii) the other being without IHS for the heat-conduction submodule (‘simple SEB + simple IHS’). Figure 10a shows that all model configurations reproduce a similar Østrem curve, despite substantial differences in model complexity. These relatively small differences indicate that a simplified representation of the SEB (and heat-conduction process) can capture the overall magnitude and shape of the melt-debris thickness relationship with good accuracy. The error boxplots, displaying the corresponding errors (observed vs modelled values) with respect to both sub-debris melt and debris surface temperature observations for the best fit, further show that differences in sub-debris melt errors are generally small (Fig. 10b). Nevertheless, with respect to the sub-debris mass balances, we observe that the ‘full SEB + full IHS’ configuration yields the best overall results (median error: 0.102 m yr⁻¹ w.e.) and the ‘simple SEB + simple IHS’ configuration exhibits the least accuracy (median error: 0.183 m yr⁻¹ w.e.). The ‘simple SEB + full IHS’ approach, which is the one used in this study, produces intermediate results (median error: 0.158 m yr⁻¹ w.e.). Differences between the degree of model complexity are more pronounced for debris surface temperatures (Fig. 10c), where the simplified configuration produces the largest errors. Again, the ‘full SEB + full IHS’ configuration yields the best results (median error: −0.389 $^\circ $C) and the ‘simple SEB + simple IHS’ configuration exhibits the largest overall errors (median error: −2.776 $^\circ $C). The ‘simple SEB + full IHS’ approach used in this study sits in between those (median error: −1.610 $^\circ $C).

Figure 10.

Comparison (a) of best-fit Østrem curve reproduced by models with the full SEB (surface energy balance) and the simplified IHS (internal heat storage) for the heat-conduction model. The ‘simple SEB + full IHS’ approach corresponds to the model configuration used in this study. Subplots (b) and (c) show the boxplots of the corresponding errors (modelled values − observations) for the surface mass balances (b_s) and debris surface temperatures (T_s), where the colours correspond to model configurations mentioned in the legend of subplot (a).

Figure 10 Long description

Comparison of model performance for three debris melt model configurations with increasing levels of physical complexity: simple surface energy balance (SEB) with full internal heat storage (IHS), simple SEB with simplified IHS, and full SEB with full IHS. (a) Simulated Østrem curves showing the melt modification factor as a function of debris thickness (m). All three model configurations reproduce the characteristic decrease in melt with increasing debris thickness and yield nearly identical melt modification factors over most of the debris thickness range. (b) Boxplots of the sub-debris mass balance errors, bs (m yr^-1 w.e.), for each model configuration. The full SEB with full IHS produces the smallest error spread, while the simplified configurations exhibit only slightly larger errors. (c) Boxplots of debris surface temperature errors, Ts (°C), showing comparable performance among the three model configurations, with the full SEB and full IHS providing the best overall agreement with observations. Overall, the results indicate that incorporating internal heat storage has a greater influence on model performance than employing a fully resolved surface energy balance.