SWAT‐GL: A new glacier routine for the hydrological model SWAT

The hydrological model Soil Water Assessment Tool (SWAT) is widely used in water resources management worldwide. It is also used to simulate catchment hydrology in high‐mountainous regions where glaciers play an important role. However, SWAT considers glaciers in a simplistic way. Although some efforts were done to overcome this limitation, there is no official version available that considers glaciers adequately. This strongly impairs its applicability in glacierized catchments. In this technical note, we propose a novel version of the traditional SWAT, called SWAT‐GL, which introduces (1) a mass balance module and (2) a glacier evolution routine to represent dynamic glacier changes. Mass balance calculations are based on a conceptual degree‐day approach, similar to the snow routine implemented in SWAT. Glacier evolution is realized using the delta‐h (∆h) parameterization, which requires a minimum of data and is thus suitable in data‐scarce regions. The approach allows users to simulate spatially distributed glacier changes. Annual mass balance changes are translated to distributed ice thickness changes depending on the glacier elevation. We demonstrate how SWAT‐GL is technically integrated into SWAT and how glaciers are merged with the existing spatial units. Model code and test data are freely accessible to promote further model development efforts and a wide application. Ultimately, SWAT‐GL aims to make SWAT easily applicable in glacierized catchments without the need of additional tools.


| INTRODUC TI ON
The Soil Water Assessment Tool (SWAT) (Arnold et al., 1998) is a widely used and recognized tool in water resources management.It is used for a vast number of applications ranging from the evaluation of land management practices (as originally developed for) (Himanshu et al., 2019;Ullrich & Volk, 2009) and sediment transport (Betrie et al., 2011;de Oliveira Serrão et al., 2022) to water quality studies (Nazari-Sharabian et al., 2019;Qi et al., 2020) or the assessment of the impacts of climate change (Schaffhauser et al., 2023;Schürz et al., 2019).SWAT is a physically based and semi-distributed model mainly applied on the watershed scale.These applications include a large number of studies in mountainous catchments, which are dominated by snow and glacier processes (Andrianaki et al., 2019;Khan & Koch, 2018;Marahatta et al., 2021;Omani et al., 2017;Rahman et al., 2012;Schaffhauser et al., 2023;Shukla et al., 2021;Tuo et al., 2018;Xu et al., 2015).While snow processes are considered by a degree-day approach, glacier processes are not specifically represented in an official release yet.
However, some efforts have already been made in the past to address this topic.One of the most widely used extension of SWAT that considers glacier processes is the approach proposed by Luo et al. (2013).Here, a volume-area scaling was integrated into SWAT, and then applied several times (Gan et al., 2015;Luo et al., 2018;Ma et al., 2015;Shafeeque et al., 2019;Wang et al., 2018) and recently transferred to the successor of SWAT, namely SWAT+ (Yang et al., 2022).
Most of the existing approaches represent glacier dynamics in SWAT either by coupling it to an external glacier model, using static approaches with unrestricted glacier melt that does not allow to represent glacier retreat or by using a volume-area scaling to simulate glacier retreat.Another way is to increase the initial snow storage to compensate for the missing glaciers (Rahman et al., 2012).Besides, none of the approaches has been considered in any of the official releases, nor are these methods always easily or openly available.
To overcome these limitations, modified versions of the traditional hydrological model SWAT are required.We thus propose a modified version which we call SWAT-GL.It consists of two features to model glacier dynamics a mass balance module as well as a glacier evolution module.Like most of the hydrological models (van Tiel et al., 2020), SWAT-GL simulates mass balance variations based on a simple degree-day approach.However, in contrast to existing approaches, glacier evolution is simulated by implementing the ∆h-parameterization first proposed by Huss et al. (2008) and further developed by Huss et al. (2010).The ∆h-parameterization is an empirical approach, where glaciers are distributed in elevation zones.The elevation zones are normalized for the minimum and maximum elevation and each zone receives a specific ice thickness change.The idea is that lower elevated zones are subject to stronger ablation than higher zones.The method is mass-conserving and can be applied with a minimum of input requirements (Li et al., 2015), making it also suitable in data-scarce regions.While the method found its way already to hydrological models such as HBV (Li et al., 2015;Seibert et al., 2018), or WASA (Duethmann et al., 2015), the community did, to the best of our knowledge, not yet integrate the approach in SWAT.However, an alternative approach is, for example, shown by Wortmann et al. (2016) for SWIM, a model which is very similar to SWAT.The glaciological response units are based on elevation zones and aspect and ice flow can occur between elevation zones.
Our technical note aims to illustrate the technical implementation of the glacier routine in the SWAT code and how glaciers can generally be incorporated into hydrological models such as SWAT.Hereby, we take up the point made by Seibert et al. (2018) who encouraged the community already to implement the approach in other hydrological models.To increase the interoperability, reproducibility, and accessibility based on the FAIR principles (Wilkinson et al., 2016), our model code is well documented and openly available via a Git repository.The technical note of the glacier routine is structured in four parts.(1) We provide details on input and preprocessing requirements in combination with the spatial integration of glaciers in SWAT.In (2), the mass balance module is described, followed by (3) the implementation of the glacier evolution module.(4) Last, we provide information where users can access, download, or further develop SWAT-GL, as well as where users can obtain the data and model of our model application.

| S WAT-G L
SWAT-GL is a revised version of SWAT considering glacier mass balance as well as glacier evolution.The core of the model is the discretization of glaciers in the aforementioned increments to allow for spatially distributed glacier changes according to the ∆h-parameterization.From now on, we use the term elevation section (ES) for these increments in which glaciers are divided.ES differs from the traditional elevation bands

Research Impact Statement
We introduce a novel version of the hydrological SWAT, called SWAT-GL, which considers spatially distributed glacier changes.
(EB) of SWAT in which snow processes are calculated and climatic inputs can be adjusted.ESs are integrated in the existing approach and only used for glacier-specific computations.Climate information here are also received from the individual EBs.
Figure 1 gives an overview on the general workflow of SWAT-GL and how it fits in the actual SWAT code.

| Preprocessing and spatial integration in SWAT
To apply the glacier model, two additional input datasets (compared to the standard version of SWAT) are required: initial glacier outlines and initial ice thickness values (from which glacier mass can be inferred).Several openly available datasets can be acquired for that task, for example, the ice thickness estimates from Farinotti et al. (2019) or Millan et al. (2022) and the glacier outlines from the RGI (RGI Consortium, 2017) might serve as appropriate for the model setup.A further approach to determine initial glacier thicknesses is to use geodetic approaches such F I G U R E 1 Flowchart of the new glacier routine implemented in SWAT-G.Hydrological response units (HRUs) are differentiated in glacierized and non-glacierized, its initialization is based on the intersection of the land use map and the glacier outlines.Glacier-specific calculations, except accumulation, only take place when the corresponding HRU is not fully snow-covered (indicated by snow cover < glacier cover).Accumulation takes place when snow is present in an HRU.HRUs are distinctly assigned to a specific elevation section (ES).Daily glacier calculations on the HRU scale are aggregated to annual mass balance estimations at the end of the ablation period upon which the ∆h-parameterization takes place.After the annual mass balance calculations are distributed over the ES, the HRU values are updated accordingly.SWAT, Soil Water Assessment Tool.
as GlabTop (Linsbauer et al., 2012).The approach makes use of an empirical relationship between the elevation range of a glacier and average basal shear stress.However, this assumes an adequate DEM (Digital Elevation Model) and glacier outlines to be available for the purpose and in best case represents the model starting year, which is probably seldom the case.
As a first preprocessing step before setting up the SWAT model, users have to define ES spacing on the basin scale.A narrow spacing allows for a more detailed representation of glacier evolution.Then, users have to merge their original land use map with the glacier outlines and prepare a land use class for each ES rather than one land use class for glaciers.In the current version of SWAT-GL, the land use class is static, which means glacierized hydrological response units (HRUs) keep its land use class regardless whether the glacier receded or not.The following descriptive example emphasizes the preprocessing step: If we assume an example in which an ES increment of 50 m would result in 10 ES, the user would need to define 10 separate land use classes.This would refer to an adaption of the land use map, its lookup table, and the SWAT crop database, in which all 10 classes must be defined accordingly.Typical values for the increment can be 50-200 m (Duethmann et al., 2015;Seibert et al., 2018).Section 3 illustrates the step for a real-world example.This step assures that glaciers and especially the corresponding glacierized HRUs get a distinct location on the glacier scale.An HRU with a land use class that corresponds to the lowest ES thus receives an indirect altitude information (as an ES covers an elevation range).In other words, there is a distinct assignment of HRUs to ESs which ensures the applicability of the ∆h-parameterization.
Area and ice thickness information are provided in a separate initialization file called swat_gles_full.txt.
The file contains the initial values for each ES and subbasin.This is important as glaciers are defined on the subbasin scale in SWAT-G.In turn, this means that users define the initial information for all glacierized subbasin in the input file.ES is defined glacier-wide, however, if we, for example, assume three glacierized subbasins and 10 ES, one needs to determine ice thickness and glacier area for each ES separately per subbasin.Detailed information on how to format the new files so that they are correctly read by SWAT-GL, as well as an example, can be found in the manual of our GitLab repository (see Section 4).Besides, a second input file with the initial parameterization for glacierized HRUs has to be provided.We hereby follow the approach to put the information of all HRUs in one file rather than having one file per HRU (as standard in SWAT).Parameters refer to glacier melt, accumulation, and sublimation, which are introduced in the next chapter.When the files are prepared and the model is set up as for the standard SWAT, it can be run.An overview of all new input and output files as well as the new parameters introduced can be found in Table 1.

| Glacier routine
When SWAT-GL is executed, first the mass balance module is called to estimate melt, accumulation, and sublimation on the daily scale.The calculations are performed on the HRU scale.Due to the unique allocation of HRUs to ESs within a subbasin, the initial values of an ES (specified in the swat_gles_full.txt)are uniformly assigned to the HRUs that belong to this ES.For this purpose, a new subroutine readglsubs_t.fwas created.Daily mass balance estimates are aggregated for each glaciological year (01.10. to 30.09.) before they are used in the ∆h subroutine, which refers to the second step of the routine.

TA B L E 1
Overview of new input files required by SWAT-GL, as well as new or modified output files generated when applying SWAT-G.The degree-day factor of ice is calculated analogously to the one of snow, which is already included in SWAT and is time varying.The Accumulation is formulated as with SWE t being the snow water equivalent of the current day [mm H 2 O] and f acc being the accumulation rate coefficient [-], which varies seasonally according to Luo et al. (2013).The described accumulation process is a very simple representation of the transition from snow to firn and ice (Luo et al., 2013;Seibert et al., 2018).Moreover, we want to mention that we introduced a simple snow redistribution scheme, which is standard in many other hydrological models such as WASA (Duethmann et al., 2015).If the SWE of an EB reaches a critical threshold, the difference between the actual SWE and the threshold SWE is shifted to the next lower EB.
Lastly, sublimation is described as function of potential evapotranspiration where ETP t is the potential evapotranspiration rate [mm H 2 O day −1 ] and represents a sublimation coefficient [-].The sublimation coefficient varies seasonally.

| Glacier evolution model
The ∆h relationship between normalized ice thickness change and normalized glacier evolution was first determined for 48 Swiss glaciers from which three general relationships have been inferred based on different glacier sizes (Huss et al., 2010).The three general relationships are illustrated in Figure 2 and serve as the core of the glacier evolution module.The parameterization is based on the assumption that ice thickness changes are stronger in lower elevated zones (ablation zone) and small in the accumulation zone.The normalization of the ES is based on the minimum and maximum glacier elevation in the basin: where E norm,i is the normalized elevation of ES i [-]; E max and E min refer to the maximum and minimum glacier elevation [m]; and E i is the actual The ∆h-parameterization uses annual mass balance changes on the subbasin scale (EW values from Equation 1 are multiplied with subbasin area to obtain volume changes).Technically, to be consistent with the original SWAT code, it is implemented in the virtual.fsubroutine, where other subbasin-specific aggregations also take place.Annual hereby refers to the glaciological year.At the end of each glaciological year, the mass balance change serves as input to the ∆h module.Therefore, a corresponding deltah_t.fsubroutine is called from the virtual.ffile, on the 1st of October of each year.
The mass balance change is distributed over the ESs of a subbasin to simulate either ablation or accumulation (of existing ESs).Glacier advance is not yet considered, similar to the original description by Huss et al. (2010).Accumulation is, therefore, restricted to glacierized areas of the corresponding time step.However, upcoming releases of SWAT-GL will consider advance, for example, following the approach of Seibert et al. (2018).
The ice thickness change follows one of the three equations of Figure 2, which in its general form is written as with a, b, c, y being coefficients differing for the glacier size classes and Δh i being the normalized ice thickness change for ES i. Users should note that specific coefficients, different from the ones proposed by Huss et al. (2010), can theoretically be determined to describe local conditions more adequately.However, the required information to derive the coefficients must be available for a specific glacier, usually two DEMs from different points in time.In the next step, a scaling factor f s [m] which converts the dimensionless ice thickness change ∆h i into an actual change is determined.
The scaling factor, in combination with the sum of the individual ES areas multiplied with the respective ∆h values, must equal the calculated annual mass balance change of a subbasin.To obtain the annual varying scaling factor, the following equation is used.
where V a refers to the annual glacier volume change expressed as water equivalent [m 3 ], determined by multiplication of the annual aggregated EW values (Equation 1) with the respective subbasin area.A i refers to the area of ES i and n is the total number of ES.
The scaling factor is then used to update the actual ice thickness of each ES h i,1 hereby refers to the updated ice thickness [m water equivalent] after a glaciological year in ES i , while h i,0 represents the ice thickness [m water equivalent] in the same ES before the ∆h-parameterization was applied.If h i,1 becomes zero, the ES is assumed to be ice-free and the glacier extent is updated accordingly.
After the ice thickness and the glacier area are updated, the HRU ice water equivalents are updated accordingly.
Empirical relationship between ice thickness change and elevation based on Huss et al. (2010).The mathematical description of the relationship is based on the glacier size and three cases can be distinguished: large glaciers-red; medium size glaciers-blue; small glaciers-purple.

| Model and test setup
SWAT-GL was tested and evaluated in the small and glacierized Martelltal catchment in South Tyrol in Italy.Our test catchment has an area of approx.77 km 2 with a glaciated area of 11.7 km 2 (approx.15%) based on the Randolph Glacier Inventory V6 (RGI Consortium, 2017).
The elevation ranges from 1800 m to over 3700 m.Ten ESs were defined based on a spacing of 100 m with the first ES starting from roughly 2666 m.Ten of 32 delineated subbasins were (partly) glaciated.We used discharge data of two gauges as well as glacier mass balance and area data to calibrate and validate SWAT-G.While daily discharge data of the Zufall Hut were mostly available during spring and summer, the discharge data at the basin outlet (Zuftritt Reservoir, see Figure 3) were artificially generated from the elaboration of the hydropower station data.This discharge time series was already used in a previous study (Puspitarini et al., 2020)  Glacier mass balance measurements were available for the Langenferner Glacier from 2004 to 2021 which were here used to exemplify the capabilities of the glacier routine (Galos et al., 2017).Mass balance observations from further glaciers were not available.However, the model setup included all glacierized parts of the basin and was not restricted to the Langenferner.See Figure 3 for the location of the Langenferner and the other glacierized parts in the basin.Meteorological data were acquired from the Autonomous Province of Bozen/ Bolzano-South Tyrol.
Calibration in the example application was performed in a multi-objective way for glacier mass balance of the Langenferner and discharge at the two gauges.An automatic scheme using multiple iterations of Latin-Hypercube samples was used to identify reasonable parameter ranges.The automatic procedure was partly extended by manual calibration runs based on expert knowledge.Used parameters varied between the calibrated subbasins, but included all snow and (introduced) glacier-related parameters.Performance metrics provided in the following paragraphs refer to one of the best solutions of the last iteration and serve only as an example.
Glaciers were initialized using the ice thickness estimates from Farinotti et al. (2019) and glacier outlines from RGI (RGI Consortium, 2017).
The assumption is that this combination might be representative in our application that spans a relatively short period in the 2000s.However, for simulations starting decades ago, this assumption is likely to not represent former conditions adequately and approaches such as from Linsbauer et al. (2012) can offer valuable solutions.

| Test results
Figure 4 summarizes the results of an example simulation in which the model was calibrated and validated for discharge and glacier mass balances.In our example, the model was evaluated by means of the NSE (Nash-Sutcliffe Efficiency), KGE (Kling-Gupta efficiency), and PBIAS (Percentage Bias).Due to differences in data availability, the two gauges were evaluated for different periods.SWAT-GL calculates mass balances on the subbasin scale, which means that actual observed mass balance values of a glacier can only be used directly if a subbasin almost perfectly matches the outline of that glacier, which is difficult in practice.In our watershed delineation, we ended up with a subbasin which is larger than the outline of the Langenferner and thus also included adjacent glaciers.We thus assessed glacier mass balance changes using annual mass balance anomalies rather than absolute mass balance values.As our subbasin includes adjacent glaciers, the mass balance observations that only cover the Langenferner would mismatch.The assumption is therefore that the relative changes of the Langenferner observations and the corresponding subbasin simulations should be in the same order of magnitude.Anomalies were calculated with respect to the period of 2008-2012.
Discharge calibration and validation results (Figure 4a-d) for KGE are 0.90 (Zuftritt) and 0.79 (Zufall Hut), 0.89 (Zuftritt), and 0.84 (Zufall Hut) for NSE and the PBIAS ranged from −11.6% to 5.5% at both gauges (calibration and validation), respectively.In general, the results were slightly better for the Zufall Hut compared to the Zuftritt Reservoir.The degraded results at the reservoir outlet stem mainly from an overestimation in simulated discharge at the beginning of the high flow period in May.However, as the observed time series was artificially reconstructed, mismatches are inherent.We further can see that the glacier routine is able to reproduce the general catchment behavior, with Spring flow being dominated by snowmelt while glacier melt is dominant in late Summer (Figure 4f) at the Zufall Hut.At the Zufall Hut, the mean annual contribution of the cumulative snow and glacier melt to the total runoff is around 70%.
With respect to the glacier mass balance anomalies of the Langenferner (Figure 4g), the general pattern was relatively well represented in our routine.Besides, the model was able to slow down the present recession trend from 2012 to 2014 (Figure 4g,h).However, the glacier routine nearly produced an equal cumulative mass balance loss over the simulation period as shown in the observations (Figure 4h).The correlation (Pearson) of the anomalies is around 0.99.The core of the ∆h-parameterization in SWAT-GL is visible in Figure 4e.It is shown that the lower elevated glaciated parts of the basin are subject to a stronger glacier retreat while upper parts of the basin remain stable.It can be seen that lower ES faced a decrease in ice volume of more than 80% in this example.

| CON CLUS I ON S AND OUTLOOK
We present a new and openly available version of the hydrological model SWAT, which is not only capable to simulate glacier mass balances on the catchment scale, but considers glacier evolution based on the ∆h-parameterization.
The mass-conserving approach allows users to represent spatially distributed glacier changes.It overcomes limitations of hydrological models where glaciers are either not integrated at all, or areal glacier changes are not taken into account (e.g., glacier evolution).The latter case can lead to unrealistic runoff contributions, especially when substantial glacier retreat phases are simulated.In contrast to volume-area scaling, the ∆h-parameterization explicitly considers basin-wide ice thickness distributions and enables elevation-dependent ice thickness changes which are translated into distributed area changes.Volume-area scaling only provides the average ice thickness for each glacier and cannot account for a number of parameters that have a large influence on glacier volume, such as surface slope or climatic conditions that affect mass balance (Helfricht et al., 2019).Furthermore, ∆h-parameterization can be applied to individual glaciers, which is not recommended for volume-area scaling (Bahr et al., 2015).
Moreover, our application demonstrates the capabilities of SWAT-GL, although the calibration scheme does not explicitly consider snow processes, for example, in the form of snow cover maps.Even though our small test showed generally good agreements in the mass balance patterns as well as the cumulative changes of the Langenferner, weaknesses were discovered in the accumulation phase of 2012-2014 where the model slowed down the recession but without reversing it.The general evaluation and benchmarking of SWAT-GL need clearly more and comprehensive investigations which is an ongoing task, while we here wanted to focus on the technical details.We highly recommend the application of any of the existing solutions to represent glacier processes in glaciated catchments, when SWAT is the favorable model to be used.In recognition of the barriers to use SWAT in these catchments, such as the inaccessibility of existing solutions or the requirement to use and manually couple additional (glacier-specific) software, our work can contribute to circumvent these obstacles and expand the applicability of SWAT.We also want to encourage the community to make model code available to promote model development, so feedback, reuse, or further development of our own code is highly appreciated.Generally, improved access to existing and new software would facilitate reproducibility and promote benchmarking in hydrological studies.
However, compared to Luo et al. (2013), glacier advance is not yet included but will be included in one of the next releases.A promising approach is to follow the suggestion of Seibert et al. (2018), which allows to include glacier advance within the initial glacier outline.Care should thus be taken if SWAT-GL should be used for long-term simulations in the past where glacier advance might be significant.Besides, the current study focuses only on one example in Italy.A comprehensive benchmarking and verification of SWAT-GL is thus missing so far.
As the focus in this technical note is on the presentation and implementation of the new module, future studies will focus explicitly on the question of SWAT-GL's applicability in and transferability to other glaciated catchments across the world.A current limitation of the approach is that, particularly for narrow ES spacing, the land use thresholds in the HRU definition are recommended to be close to 0 to allow a best possible representation of all ESs and to avoid an accidental removal of ES.However, this in combination with the general approach to define one land use class per ES could result in relatively high numbers of HRUs, especially in large catchments.Besides, land use classes are static in our current version, whereas changing land use classes in periods of full glacier retreat would be more realistic.An approach could be further developed to define two land use classes for glacier areas, one for periods of glaciation and one for ice-free times in future studies.However, parameter stationarity assumption is a general issue that should be given more consideration in hydrology, but which is beyond the scope of this work.Further difficulties refer to the data availability for the initialization as that information is rather scarce and usually not available at multiple points in time for a specific region.However, global available datasets such as the ice thickness dataset from Farinotti et al. (2019) or Millan et al. (2022) provide reasonable solutions.Besides, fragmented glaciers within a subbasin are treated as one construct, a limiting but necessary simplification in semi-distributed hydrological models.One solution for the consideration of debris-covered glaciers would be the adjustment of the melt factor of ice.Different studies (e.g., Muhammad et al., 2020) concluded that the degree-day-factor for ice melt is reduced through debris cover.However, a sound sensitivity study for glaciers with different degrees of debris cover is necessary for defining an appropriate parameter range.This was (also) not within the scope of this study, but we encourage the scientific community to test the new glacier routine under different conditions.
We also recognize that future model development efforts should be put into the newer successor model SWAT+ (Bieger et al., 2017).
Future versions of SWAT-GL that address current shortcomings will thus increasingly focus on SWAT+.However, the spatial integration of glaciers in SWAT+ will likely differ from our current implementation, due to the revised spatial structure of SWAT+.
corresponding equation is where b gmlt,mx and b gmlt,mn correspond to the maximum and minimum degree-day factors for June 21 and December 21 [mm H 2 O °C −1 day −1 ].In theory, degree-day factors of ice are usually higher than those of snow.An overview of the magnitude of degree-day factors can be found inSingh et al. (2000).If degree-day factors of snow exceed those of ice, SWAT-GL automatically corrects the latter and uses the factors of snow for ice as well in order to preserve the physical relationship.
to calibrate and validate a hydrological model in Martelltal.Details on the discharge data of the Zuftritt Reservoir can be found in Puspitarini et al. (2020).F I G U R E 3 Overview of the study area, (a) elevation distribution, (b) land use map showing an example where 10 ES (indicated as 341-350 in the legend) are separately considered as distinct land use classes.Originally, there was only one global land use class for glaciers and no separate land use class for a specific elevation range was considered.

F I G U R E 4
Example application of SWAT-GL in the Martelltal: (a) Daily discharge simulations at the Zufall Hut for both the calibration and validation period.Gaps in the time series indicate missing observations; (b) mean seasonal discharge simulations based on the combined calibration and validation period; (c) same as (a) but for the basin outlet; (d) same as (b) but for the basin outlet; (e) elevation-dependent glacier volume change for the whole simulation period.Changes are shown as the volume at the beginning (gray) and end (cyan) of the simulation for the different ESs.The initial volume was derived from the ice thickness estimates of Farinotti et al. (2019); (f) mean monthly simulated runoff contributions of snow and glacier melt on the basin scale (no observations); (g) model performance of the glacier routine evaluated as glacier volume changes for the Langenferner glacier.Results are shown as annual volume anomalies relative to the period 2008-2012, as the subbasin covering the Langenferner exhibits a larger glacier extent as the Langenferner glacier alone.(h) The cumulative mass balance change expressed as surface elevation change for the Langenferner and the covering subbasin.
Jiang et al. (2010)lcuBonekamp et al., 2019;h the existing snow-based computations (snom.f),withanexceptionforsublimationwhich is estimated in the original etact.froutine,wherealsosnowsublimation is considered.The mass balance can be formulated as where EW t is the water equivalent of ice [mm H 2 O] at day t, M t is the melt rate [mm H 2 O day −1 ], S t is the sublimation rate [mm H 2 O day −1 ], C t is the accumulation rate [mm H 2 O day −1 ], and f refers to a refreezing factor of melt, for example also introduced inLuo et al. (2013)where a value of 0.2 fromJiang et al. (2010)was stated.In our case, f is a variable calibration parameter that can be used to reduce high melt rates.It has been reported that refreezing of melt water in glaciers or ice sheets can have a considerable contribution(Abrahim et al., 2023;Bonekamp et al., 2019; GLMFMXInput Melt factor for ice on June 21 [mm H 2 O/(°C•day)] GLMFMN Input Melt factor for ice on December 21 [mm H 2 O/(°C•day)] f Input Refreezing factor of glacier melt [-] f acc Input Conversion factor of snow to firn and ice [-] 2.2.1 | Mass balance model with T mx,t being the daily maximum temperature [°C], T gmlt being the threshold temperature [°C] glacier melt to occur, b gmlt refers to the degree-day factor of ice [mm H 2 O °C −1 day −1 ], and A sc and A gc refer to the snow and glacier covered fractions, respectively.Glacier melt is only initiated when the daily maximum temperature reaches the critical threshold temperature and a glacier HRU is not fully snow covered.