Firn changes at Colle Gnifetti revealed with a high-resolution process-based physical model approach

. Our changing climate is expected to affect ice core records as cold ﬁrn progressively transitions to a temperate state. Thus there is a need to improve understanding and further develop quantitative process modeling, to better predict cold ﬁrn evolution under a range of climate scenarios. Here we present the application of a distributed, fully coupled energy balance model, to simulate high-alpine cold ﬁrn at Colle Gnifetti over the period 2003–2018. For the ﬁrst time, we force such a model with high-resolution, long-term and extensively quality-checked meteorological data measured in closest vicinity of the ﬁrn 5 site, at the Capanna Margherita (4560 m a.s.l.). The model incorporates the spatial variability of snow accumulation rates, and is calibrated using several, partly unpublished high-altitude measurements from the Monte Rosa area. The simulation reveals a very good overall agreement in the comparison with a large archive of ﬁrn temperature proﬁles. The rate of ﬁrn warming at 20 m depth is estimated at 0.44 ◦ C per decade. Our results show that surface melt over the glaciated saddle is increasing by 3–4 mm w.e. yr -2 depending on the location (29–36 % in 16 years), although with large inter-annual variability. Analysis of modeled 10 melt indicates a marked tendency towards small melt events (< 4 mm w.e.), which collectively represent a signiﬁcant fraction of the melt totals. Atmospheric humidity is found to be a prominent control over melt occurrence, with considerable amounts of sublimation happening in dry conditions. For future developments, we recommend implementing a physical simulation of meltwater inﬁltration, to be calibrated with more measurements of percolation

for more information). In both panels X and Y are metric CH1903/LV03 coordinates. Orthophoto and topographic map source: Federal Office of Topography swisstopo. The consistent availability (> 95 %) of long-term data (since mid-2002) makes the series a valuable data-set for highalpine research. Still, an extensive evaluation and processing of the data is essential due to the extreme measuring conditions (Martorina et al., 2003). In particular, freezing of the anemometer -as well as snow and ice accumulation interfering with the pyranometer -are relatively common occurrences. Therefore, the hourly data from eight other high-altitude AWS located in the region (Fig. 1b, id 2-9 in Table 2) were used to perform quality checks, to fill gaps in the CM time series, and to provide 95 data for parameters not measured at CM. The selection of stations to include was determined by the availability of different parameters measured at each site, by the effort to provide an unbiased geographic coverage in all directions from the CM AWS, and by the station tendency towards concurrent failures in challenging conditions, such as winter storms.
All AWS series were pre-processed to remove clock errors, detected with cross-correlation analysis. Then each series was entirely reconstructed from the data of the best correlated others using quantile mapping (Cannon et al., 2015;Feigenwinter 100 et al., 2018), to provide a reference for robust outlier detection. High-resolution reanalysis series of the parameters measured at CM were also collected as an additional basis for comparison (COSMO-REA2 and COSMO-REA6: Wahl et al., 2017;Bollmeyer et al., 2015;Frank et al., 2018). Due to the large volume of data involved, an automated pre-filtering routine was implemented, based on objective criteria (absolute values, rates of change, comparison with reconstructed series and reanalysis) to mark single values as potential outliers, which were then manually checked. The CM AWS was always processed last in 105 order to be compared to the highest-quality data. After quality check, all gaps in the hourly CM series were filled with the corresponding reconstructed values. The close match of the reconstructed series (Table 3) despite the large elevation differences involved (Table 2) confirms the benefit of the extensive processing described.
The CM AWS does not record series of relative humidity (RH) or cloudiness, which are required within the model input.
Thus RH was entirely derived from the other quality-checked time series; cloudiness was reconstructed from the incoming 110 SW radiation at CM when available, and otherwise computed from the incoming LW radiation and observed cloud cover, respectively at Stockhorn and Plateau Rosa. Lapse rates of temperature and pressure -used to extend the AWS record to the gridded model domain -were computed from the difference between CM and the other stations. Further details on the series  processing can be found in Mattea (2020). Figure 2 shows the final time series of all meteorological parameters assigned to CM and used as model input.

115
To determine some of the model calibration parameters (Table 4), we also re-analysed data from two additional weather and energy-balance stations (10 and 11 in Table 2), temporarily installed above 4000 m a.s.l. within the ALPCLIM project (Auer et al., 2001). The Seserjoch station operated from September 1998 to October 2000 for the detailed SEB investigations of Suter et al. (2004); the Colle del Lys station (Rossi et al., 2000a,b) was active between 1996 and 2000.

Sub-surface data
For model validation (Sect. 4.1) we used archived firn temperature profiles (Fig. 1a), measured between 2003 and 2018. Beside the CG saddle, we included two profiles from the Grenzgletscher slopes and four from Seserjoch. In total, we considered 25 temperature profiles from 18 boreholes, some locations having been measured more than once. Detailed description of the profiles is reported in Hoelzle et al. (2011) andGLAMOS (2017). Earlier information on the deep temperatures (Haeberli and Funk, 1991) was also used as reference during model set-up (Sect. 3.4).

125
Moreover, we compiled point measurements of annual accumulation (Fig 3a), derived from layer thicknesses in GPR profiles (Konrad et al., 2013) and from the stake network of Suter and Hoelzle (2002), as well as from archived ice core measurements (densities and annual dating) acquired at CG between 1982 and 2019. In total, 14 core profiles were available; only one (core Zumsteinkern, from 1991) is located on the high-accumulation, south-facing Zumsteinspitze slope (Fig. 3a). Description of the individual cores can be found in Licciulli (2018) and Lier (2018).

130
Finally, we hand-drilled a 5.5 m core (unifr-2019, Fig. 1a) near the CG saddle point on June 25, 2019, analyzing density and stratigraphy in the field. Core description is presented in the appendix.

Topography
The and long-wave (LW ) radiation, sensible (Q SH ) and latent (Q LH ) turbulent fluxes, heat advection from rainfall (Q rain ) and heat conduction into the snow or ice (Q g ). Then the SEB (Eq. 1) is solved for surface temperature and melt amounts: these, together with the lower boundary condition of geothermal heat flux, drive the sub-surface evolution.
Simulation of surface processes is developed along the lines of Klok and Oerlemans (2002), while the multi-layer sub-surface snow model is based on the SOMARS approach (Simulation Of glacier surface Mass balance And Related sub-surface processes, Greuell and Konzelmann, 1994). In this work, the model version described by van Pelt et al. (2019) was used, with a parametrized water percolation routine simulating preferential flow (Marchenko et al., 2017), and an updated scheme for albedo decay based on Bougamont et al. (2005). This model participated in the firn meltwater Retention Model Intercomparison Project (RetMIP) under the designation "UppsalaUniDeepPerc" (Vandecrux et al., 2020).
In the following we highlight the main EBFM routines and their respective adaptations to the CG setting. The model was originally developed over large, polythermal Arctic glaciers: set-up for a high-alpine cold firn saddle was possible thanks to a rich archive of energy-balance measurements from the Monte Rosa area.

Radiative fluxes 155
The model computes incoming SW radiation as where T OA shaded (x, y) is the unattenuated top-of-atmosphere radiation corrected for topographic shading, t rg , t w , t a and t cl gaseous, water vapor, aerosol and cloud transmissivities, n fractional cloud cover, and a and b calibration parameters. Values for 160 a and b (Table 4) were derived from Greuell et al. (1997) as calibrated in high-alpine terrain, compared to the EBFM defaults which were tuned from measurements in the Arctic. The choice was supported by a distribution analysis of the incoming SW flux, measured at the CM AWS under overcast sky conditions.
Reflected SW radiation is controlled by a broadband, isotropic surface albedo (Eq. 4). Albedo evolution is modeled after Oerlemans and Knap (1998) as an exponentially decaying function of time since last significant snowfall (Eq. 5), bounded 165 by constant values for fresh snow and firn. Decay time-scale is a function of snow surface temperature (Eq. 6) to account for slower metamorphism in cold conditions (Bougamont et al., 2005;van Pelt et al., 2019):

170
where t * is the computed albedo decay time-scale, t * wet (t * dry ) time-scale for a melting (dry) surface at 0 • C, K a calibration parameter, T snow surface temperature and T max,t * a temperature cut-off value for decay slow-down. Values for these parameters (Table 4) were estimated from a new analysis of albedo evolution, measured at the Seserjoch and Colle del Lys stations respectively in 1999 and between 1996 and 2000. Incoming LW radiation is computed with the Stefan-Boltzmann law for grey-body radiation (Eq. 7); sky emissivity is mod-175 eled after Konzelmann et al. (1994) as a function of cloud cover, air temperature and humidity (Eqs. 8 and 9): e = e cs (1 − n 2 ) + e cl n 2 (8) where e cs is clear-sky emissivity, b a calibration parameter, V P vapor pressure, T air temperature, e sky emissivity, n fractional 180 cloud cover and e cl cloud emissivity. Parameters b and e cl were estimated for CG from radiation measurements at Seserjoch.

Turbulent heat fluxes
In the EBFM, turbulent heat exchange is modeled with the glacier katabatic wind parametrization of Oerlemans and Grisogono (2002). This was developed with a focus on large valley glaciers; notably, it computes heat fluxes which are independent on the ambient wind field, since wind speeds are estimated from the katabatic flow model (Oerlemans and Grisogono, 2002). This 185 was deemed inadequate for a high-alpine, wind-exposed saddle: thus the EBFM computation of turbulent heat exchange was re-implemented, following the bulk aerodynamic equations of Essery and Etchevers (2004) computed as where ρ a is air density, c p specific heat of dry air, V wind speed, T air temperature, L s,v latent heat of sublimation or vaporization (chosen depending on the modeled surface temperature), q specific humidity, and the s and 1 subscripts refer respectively to the snow surface and the measurement level (2 m). Exchange coefficient C h is defined as

195
where C hn is the value under neutral conditions and f h a correction for atmospheric stability, expressed in terms of the bulk Richardson number Ri B : Here, k is the von Kármán constant, z 1 the measurement level, z 0 the surface roughness length, g the gravity acceleration, and the ratio of molecular weights for water and dry air.

Precipitation model
The model precipitation routine was adapted to reproduce the extreme spatial gradient of snow accumulation distinctive of 205 the site (Sect. 1). Since the EBFM does not include a blowing snow routine, the simple model of linear precipitation rates with altitude (van Pelt et al., 2019) was replaced by a gridded precipitation time series, already corrected for snow lost to wind scouring. This was computed with a three-phase anomaly method inspired from New et al. (2000), by combining a fixed climatological grid, an annual anomaly series, and a temporal down-scaling coefficient. Specifically, for grid cell (x, y), at simulation time-step t, in year i, precipitation was expressed as where C is the long-term annual accumulation climatology, A the domain-wide annual anomaly, and D the down-scaling coefficient. The main assumption of the method is that spatial patterns of relative accumulation do not change over time. The climatological grid C was assembled by interpolating point values of long-term net accumulation, estimated from the snow mass of dated firn cores and from the mean layer thickness in GPR profiles (Fig. 3a). To extend data coverage and reduce the 215 occurrence of extrapolation, stake measurements from Suter and Hoelzle (2002) were also used in the western and southern domain regions. Accumulation measured at each stake over a single year was re-scaled to a mean annual estimate by using the overlaps with firn cores and GPR points.
For annual anomalies A, the multiplicative snow mass anomaly of the KCC deep core (Bohleber et al., 2018) was found to be moderately anti-correlated with wind speed measured at the CM AWS (Fig. 3b). Linear fit over 9 annual data points yielded 220 the formulation where A is core snow mass anomaly and V the median wind speed observed at CM over the corresponding year. Anomaly from the KCC core was then applied to the whole domain.
Finally, the down-scaling coefficients D were computed by normalizing to a unit sum each year of hourly precipitation, 225 averaged over the three closest rain gauges (Passo Monte Moro, Bocchetta delle Pisse and Rifugio Zamboni: Fig. 1b, Table 2).
Equation 17 produces an hourly series of gridded precipitation (already corrected for wind erosion) which was used to force the EBFM.

Water infiltration
The EBFM features a parametrized water infiltration routine, which routes any liquid water from the snow surface towards 230 depth, instantly and according to a constant vertical distribution (Marchenko et al., 2017), as long as the modeled near-surface layers are not impermeable. Water at depth can then refreeze, be held by capillary forces as irreducible water, or form slush layers which drain gradually (van Pelt et al., 2012). In this study, the water vertical distribution was chosen to follow the normal law; the maximum depth z lim was tuned to a value of 4 m, to match the firn temperature of the CG saddle point at a depth of 20 m (Haeberli and Funk, 1991).

Results
The EBFM computes and logs a wide variety of surface and sub-surface variables. In the following, we focus on firn tempera-245 tures and melt amounts as relevant descriptors of current cold firn evolution. For these variables, a conspicuous archive of field measurements and model estimates exists at CG, allowing verification and comparison of our results.

Firn temperatures
Comparison of modeled firn temperatures to measured borehole profiles (Table 5, Fig. 4 and 5) shows that most model deviations are similar in magnitude to the spatial variability of firn temperatures. Indeed, profiles CG08-1/08 and CG08-2/08 north-facing slopes, and too warm in the flat or south-facing regions of CG and Seserjoch. Moreover, model biases appear to be maintained over time for boreholes with repeated measurements (e.g., CG05-1 and CG13-1). The time series of modeled firn temperatures (Fig. 6) shows large differences over rather small distances at the CG saddle,

Melt amounts and dynamics
Modeled mean annual melt amounts have an extreme spatial variability (Fig. 7), broadly reflecting surface elevation, slope and aspect. Values increase from less than 1 cm w.e. yr -1 on the steepest slopes of the Signalkuppe, to 17 cm w.e. yr -1 at the saddle point, and about 23 cm w.e. yr -1 on the Zumsteinspitze slope. Even higher melt amounts, exceeding 30 cm w.e. yr -1 , are simulated for the lower elevation Grenzgletscher slopes (towards the western border of the domain) and for Seserjoch. Grid 270 average is 21 cm w.e. yr -1 . Within the overall energy balance, melt represents a relatively minor component (Fig. 8): the largest mean monthly contribution, in August, is well below 10 %, and in every month sublimation is a more effective energy sink than melt. Still, in the NE domain region -where wind scouring is strongest -annual melt amounts correspond to a significant fraction of net accumulation (Fig. 3a). With regard to temporal patterns, the entire surface was found to always refreeze at night over the modeled period. Moreover, no melt is simulated between November and March, with only minor amounts in April 275 and October (Fig. 8). Despite the large spatial heterogeneity, and an inter-annual variability exceeding 50 %, a common trend of melt increase could be detected in the annual time series (Fig. 9): the fitted slope ranges between (3 ± 2) and (4 ± 3) mm w.e. yr -2 across the saddle. While the trend is somewhat masked by the 2003 extreme melt year at the very beginning, it becomes statistically significant over the rest of the period.

280
The EBFM shows a marked tendency towards small melt amounts: frequency of modeled melt events decays exponentially with their magnitude (Fig. 10a), and a significant fraction of total melt amounts is contributed by micro-melt events under 4 mm w.e. (Fig. 10b).  Investigation of the weather conditions leading to melt occurrence (Fig. 11) reveals the relationship between weather variables and surface melt. Air temperature provides a critical control over melt rates and total amounts, unlike cloud cover which 285 appears to have almost no effect ( Fig. 11a/      The spatial distribution of model residuals (Fig. 5) could be affected by the lack in the modeled SEB of radiation reflected 300 from the surrounding terrain. This process is expected to induce a net energy transfer from the more sun-exposed cells towards the more shaded ones, which could partially re-balance the aspect dependence of model residuals.
In addition, boreholes with positive model bias are often simulated with too strong near-surface temperature gradients, leading to sharp positive deviations in the profile near the surface (e.g., Fig. 4f/k/p/q). Such a behavior suggests too deep refreezing is occurring in the simulation, which could be linked to the parametrized preferential infiltration routine (Sect. 3.4).

305
Indeed, below a depth of 0-4 m (where refreezing is occurring) the simulation appears to be unbiased at most locations.
Model under-estimation of firn temperatures appears to correlate with the lowest mean annual accumulation values ( Fig. 3a and 5). At these locations wind scouring is strongest, and melt amounts represent a significant fraction of net annual accumulation. Simple sensitivity tests showed that the simulation of profiles at borehole location CG05-1 improves substantially by forcing increased precipitation amounts. This suggests that the accumulation model -which computes precipitation amounts 310 already corrected for losses from wind scouring -may be problematic whenever melt totals approach the value of net accumulation. Indeed, accumulation at CG results from summer precipitation events (Sect. 1), but modeled precipitation is distributed more evenly throughout the year, as it is based on weather station measurements (Sect. 3.3). Thus we expect modeled precipitation to be under-estimated in summer compared to the actual accumulation. As such, in summer modeled melt can temporarily approach (or even locally exceed) the low accumulated snow amounts. Then the parametrized percolation would generate 315 dense, thick firn layers: their limited refreezing capacity could reduce the heat release at depth during subsequent melt events, inducing under-estimation of the englacial temperatures.
Among past modeling efforts at CG, both Lüthi and Funk (2001) and Suter (2002) formulated independent predictions for firn temperatures evolution by 2020. As the target time frame for verification is reached, a major limitation of their modeled scenarios is found in the expected magnitude of atmospheric warming: while they assumed linear air temperature increases  firn warming is consistent with the results of the two studies, but lower when compared to the assumed air temperature increase.

325
A contributing factor could be the EBFM spin-up: temperatures were initialized with repeated model runs over 2004-2011, thus the initial grid at all depths is in equilibrium with the mean forcing over that period. By contrast, an adjustment time of 2-4 years is to be expected at the considered depth (Hoelzle et al., 2011). Therefore at the beginning of the simulation the EBFM may be slightly over-estimating deep firn temperatures ( Fig. 4a/b), resulting in a lower trend for 2003-2018.

Meltwater infiltration 330
A key parameter used for model set-up is the percolation depth z lim , which defines the maximum depth reached by infiltrated meltwater through preferential flow. In cold firn such a parameter is crucial: simulated firn temperatures (Fig. 6) indicate that all meltwater can be expected to refreeze not far from the initial location (controlled by the parametrized vertical distribution: Sect. 3.4). Thus z lim effectively determines not only the initial meltwater distribution, but also the depth of refreezing heat release.
The calibrated value of 4 m, together with the Gaussian vertical distribution, correspond to a mean preferential percolation 335 depth of 1.06 m. Unfortunately, in situ quantitative measurements of infiltration depths are scarce and mostly indirect. Evidence from winter snow packs and glacier accumulation areas shows that percolation and refreezing are strongly dependent on the meltwater supply amounts and rates (Marchenko et al., 2017), the temperature of the firn matrix (Koerner, 1970) and its stratigraphy (Illangasekare et al., 1990), notably affected by the previous history of infiltration and refreezing. An obvious consequence is the increase of percolation depths over the melting season (Marchenko et al., 2017). At CG, Alean et al. (1983) 340 observed slight melting but no meltwater percolation during a warm spell in summer 1981. Suter (2002) attempted to directly track meltwater refreezing by continuously logging a temperature profile at Seserjoch in 1999, but the setup failed before the onset of summer melt. Hoelzle et al. (2011) reported evidence for increasing infiltration depths, indicating (at Seserjoch and on the Grenzgletscher slopes) a transition to a percolation regime spanning several annual firn layers. By contrast, the stratigraphy of core unifr-2019 (Fig. A1) revealed dry firn layers interspersed with thin ice crusts, suggesting small amounts of infiltration 345 and refreezing, except for a thick, ice-rich layer at 4.5 m depth. Therefore in our simulation we consider z lim more as a tuning parameter than a realistic percolation depth.
Stratigraphy of the firn core points to a shortcoming of the subsurface model: meltwater distributed in cold firn along the parametrized vertical profile usually refreezes in place, producing a diluted density increase in the simulation. Instead, distinct ice layers are known to form at depth after preferential percolation, affecting the mechanical, hydrological and thermodynam-350 ical properties of the snow pack (Quéno et al., 2020). In the present EBFM formulation, such ice layers would not appear even with a very fine model grid. Refreezing after parametrized infiltration also distributes heat instantly over a fixed, large vertical extent, resulting in unrealistic, frequent warm pulses at depth after each melt event -no matter how small (inset in Fig. 6).
These observations indicate the need for a physically based percolation routine in the EBFM, accounting for the time evolution of z lim and its dependence on snow density and stratigraphy.

Melt amounts
Modeled melt amounts can be compared to the amount of refrozen ice observed in core unifr-2019 (Fig. A1). According to the computed climatology (Fig. 3a), the core location has a mean long-term accumulation of about 50 cm w.e. yr -1 : thus the 5.5 m core should span an estimated period of about 5-6 years (with a fairly large uncertainty due to the inter-annual accumulation variability). The model predicts about 85-100 cm w.e. of melt over such a time span (Fig. 7), while the core was found to 360 contain only 31 cm of ice layers. Some observations can be brought forward to put the apparent discrepancy into context.
First, refrozen ice amounts recorded in the core are affected by repeated cycles of melt-refreeze. Indeed, the very small amounts of meltwater produced during less intense (but rather frequent: Fig. 10) melt events can be expected to refreeze in the very first snow centimeters. Then any subsequent melt occurring before the next snowfall would affect the same ice surface, contributing to total melt amounts but without significant increases in ice layer thickness. Such surface crusts of relatively 365 impermeable ice have already been observed at CG (e.g., Lier, 2018). Since melt mostly happens in clustered patterns (almost only in summer and within a specific set of weather conditions: Fig. 8 and 11), contribution of repeated melt-refreeze cycles could potentially be very large.
With daily melt amounts often close to the size of single crystals (Fig. 10a), it is suggested that detection of some refrozen forms would require a resolution not achieved during our field analysis. Investigations of such minimal melt processes are very overlooked vertical ice gland embedded in our core also cannot be ruled out (measured core sections were not broken up after 375 analysis). In fact, the presence of undetected ice would be consistent with the high density variability encountered in the core profile (Fig. A1). Still, this observation could also be linked to wind compaction, and its interplay with wind erosion exposing older, denser snow.
An estimation of refreezing amounts from profile density anomalies was proposed by Lier (2018), using several new and archived cores in the CG area. The method involved sizable uncertainties due to a poorly constrained surface density value.

380
At the Sattelkern and Zumsteinkern cores (respectively at the saddle point and on the south-facing slope), estimated refreezing rates are respectively in the 1-13 and 3-33 cm ice yr -1 ranges (Lier, 2018). The EBFM predicts approximately 19 and 25 cm ice yr -1 of melt. On the shaded Signalkuppe flank, values span the range 0-15 cm ice yr -1 from multiple cores; the EBFM result over the same region is between 5 and 12 cm ice yr -1 of melt. Considering the very strong inter-annual variability of melt amounts (Fig. 9), the results are largely compatible.

385
The significant melt amounts modeled at sub-freezing temperatures suggest reconsideration of degree-day models for simulating melt at high-alpine (and possibly high-latitude) locations. This is consistent with the findings of Suter et al. (2004) at Seserjoch, who observed several surface melt events but no days with positive mean temperatures over the whole 1999 summer. Indeed, laboratory experiments by Beck et al. (1988)  incoming SW radiation. At the CM AWS values in excess of 1000 W m -2 are a common summer occurrence, supporting the plausibility of melt under even colder conditions (Fig. 11). From a theoretical perspective, Kuhn (1987) analytically explored a standard SEB equation (neglecting sub-surface heat conduction), in relation to common weather situations on an alpine glacier.
The conclusion was that melt onset likely happens at air temperatures between -10 and +10 • C, very much comparable to our results.

395
This work marked a first effort to apply a coupled, high-resolution EBFM to alpine cold firn, within a multi-year simulation forced with extensively processed meteorological data, acquired at high altitude and in closest vicinity of the study site. In both cold and near-temperate conditions, the model can achieve promising results for firn temperatures, comparable in accuracy to the spatial variability of the measurements. Therefore the EBFM can be deemed suitable for further investigations of future cold firn evolution, based on localized climate scenarios.

400
At CG, our results corroborate earlier observations on the spatial patterns of surface melt and the trends of deep temperatures.
A previously unreported feature in melt dynamics is the occurrence of micro-melt events, with daily amounts below 4 mm w.e.: these remain difficult to detect, but our analysis hints at a possibly significant impact on melt totals and firn temperatures, hence on calibration and ground-truthing of model results. More field observations are needed to verify the occurrence and improve the understanding of such events, also assessing their potential effect on subsequent water infiltration.

405
Our model results point to significant melt happening at negative temperatures, confirming earlier field observations: this would re-affirm the importance of using a full energy-balance model over a parametrized melt approach in cold conditions.
Additional field investigations of the meteorological conditions at the onset of melt would be a valuable development.
Moreover, we found a novel trend of increasing surface melt amounts, currently on the edge of statistical significance.
Further evolution of this tendency should be closely monitored, since in the near future it could affect the site suitability for 410 retrieving climate records from ice cores. Estimations of local melt amounts may also contribute to the localization of future core drilling efforts: this provides motivation to attempt model deployment at other cold firn/ice sites and potentially on a larger scale.
For such an application, it would be beneficial to better constrain patterns and depths of meltwater percolation, by implementing a physical infiltration routine and acquiring vital calibration data. In addition to thermal tracking, recent developments