Evolution of the firn pack of Kaskawulsh Glacier, Yukon: meltwater effects, densification, and the development of a perennial firn aquifer

Abstract

In spring 2018, two firn cores (21 and 36 m in length) were extracted from the accumulation zone of Kaskawulsh Glacier, St. Elias Mountains, Yukon. The cores were analyzed for ice layer stratigraphy and density and compared against historical measurements made in 1964 and 2006. Deep meltwater percolation and refreezing events were evident in the cores, with a total ice content of $\mathrm{2.33}\pm \mathrm{0.26}$ m in the 36 m core and liquid water discovered below a depth of 34.5 m. Together with the observed ice content, surface energy balance and firn modelling indicate that Kaskawulsh Glacier firn retained about 86 % of its meltwater in the years 2005–2017. For an average surface ablation of 0.38 m w.e. yr${}^{-\mathrm{1}}$ over this period, an estimated 0.28 m w.e. yr${}^{-\mathrm{1}}$ refroze in the firn, 0.05 m w.e. yr${}^{-\mathrm{1}}$ was retained as liquid water, and 0.05 m w.e. yr${}^{-\mathrm{1}}$ drained or ran off. The refrozen meltwater is associated with a surface lowering of$\mathrm{0.73}\pm \mathrm{0.23}$ m between 2005 and 2017 (i.e., surface drawdown that has no associated mass loss). The firn has become denser and more ice-rich since the 1960s and contains a perennial firn aquifer (PFA), which may have developed over the past decade. This illustrates how firn may be evolving in response to climate change in the St. Elias Mountains, provides firn density information required for geodetic mass balance calculations, and is the first documented PFA in the Yukon–Alaska region.

1 Introduction

With the increasing effects of climate change and the need for understanding glacier and ice sheet melt rates, geodetic methods are useful for indirect measurements of mass balance (Cogley, 2009). Based on repeat altimetry, geodetic approaches to mass balance monitoring rely on several assumptions. Estimates must be made of the density of snow, firn, and ice at the sampling location, with the additional assumption that these densities remain unchanged between the two measurement dates. However, over multi-annual timescales in a warming climate this may not be true (Moholdt et al., 2010a). Meltwater percolation and refreezing can significantly change the firn density profile and mean density of the accumulation zone of a glacier (Gascon et al., 2013) and can introduce large uncertainties when using geodetic techniques to determine glacier mass balance if they are not properly accounted for. For example, Moholdt et al. (2010b) determined the geodetic mass balance of Svalbard glaciers to be $-\mathrm{4.3}\pm \mathrm{1.4}$ Gt yr${}^{-\mathrm{1}}$, based on ICESat laser altimetry, with the large uncertainty attributed to limited knowledge of the snow and firn density and their spatial and temporal variability. By altering the density and causing surface lowering, meltwater percolation, refreezing, and liquid water storage all complicate the interpretation of geodetic mass balance data.

Warming firn can result in increased meltwater production and altered firn densification processes. Initially, melt can round the snow grains and increase the snowpack density and then percolate into the firn and refreeze as ice layers or lenses (Sommerfeld and LaChapelle, 1970; Cuffey and Paterson, 2010). On glaciers with medium to high surface melt, and high annual snow accumulation, meltwater that percolates below the winter cold layer often will not refreeze, and may thus form a perennial firn aquifer (PFA) if this water cannot effectively drain through crevasses or moulins (Kuipers Munneke et al., 2014). These internal accumulation processes can significantly increase the firn density, and once ice layers or PFAs form they affect how meltwater percolates through the firn pack (Gascon et al., 2013). Due to the spatial heterogeneity of meltwater retention, percolation, and refreezing processes, there are still many gaps in knowledge of how to model these processes and subsequently estimate firn density in areas where these processes occur (van As et al., 2016).

Meltwater retention in firn is also important for estimating glacial run-off contributions to sea level rise. Numerous recent studies have investigated meltwater refreezing processes in northern locations such as southern Greenland (Humphrey et al., 2012; Harper et al., 2012; De La Peña et al., 2015; MacFerrin et al., 2019; Vandecrux et al., 2020), the Canadian Arctic Archipelago (Noël et al., 2018; Zdanowicz et al., 2012; Bezeau et al., 2013; Gascon et al., 2013), and Svalbard (Noël et al., 2020; Van Pelt et al., 2019; Christianson et al., 2015). In many locations with cold deep firn, short-term increases in surface melt rates may not result in proportional increases in surface run-off due to percolation and refreezing of meltwater in the firn pack (e.g., Harper et al., 2012; Koenig et al., 2014; MacFerrin et al., 2019). However, in the long term this may lead to expansion of low-permeability ice layers, causing run-off to increase and expediting the movement of water from glaciers to the ocean (MacFerrin et al., 2019; Machguth et al., 2016). Current knowledge of these processes is limited for mountain glaciers in other regions, although this information is required for improved estimates and models of glacier mass balance and associated sea level rise.

In this study two firn cores were retrieved in spring 2018 on Kaskawulsh Glacier, St. Elias Mountains, Yukon, and analyzed for density and the effects of meltwater percolation and refreezing. We also use a firn model (Samimi et al., 2020) forced by bias-adjusted ERA climate reanalyses to investigate the evolution of the firn through time at this location. Comparisons of these measurements with firn density profiles collected at a nearby site in 1964 and 2006 enable us to (i) quantify contemporary firn characteristics and densification processes, (ii) determine how the physical properties of the firn pack have changed over the past $\sim \mathrm{50}$ years, and (iii) assess the potential widespread presence of a PFA on the upper Kaskawulsh Glacier.

2 Study area

The St. Elias Mountains are located in the southwestern corner of Yukon Territory, Canada, and contain many peaks higher than 3000 m, including the highest mountain in Canada, Mount Logan, at 5959 m a.s.l. (Fig. 1). The St. Elias is home to the largest ice field outside of the polar regions, with an area of $\sim \mathrm{46}\mathrm{000}$ km${}^{\mathrm{2}}$ (Berthier et al., 2010). Measurements presented here are focused on the upper accumulation zone of Kaskawulsh Glacier (Fig. 1), which is part of an extensive ($\sim \mathrm{63}$ km${}^{\mathrm{2}}$) snowfield at an elevation of 2500–2700 m a.s.l. This plateau region has subtle topographic variations and includes the drainage divide between the Kaskawulsh and Hubbard glaciers.

Figure 1

Field locations in the St. Elias ice field, Yukon. The IRRP A site is the site of the 1964 firn core that is referenced in our study (Grew and Mellor, 1966). Base map from http://openmaptiles.org/ (last access: 17 March 2021).

[Figure omitted. See PDF]

Kaskawulsh Glacier is a large valley glacier located on the eastern side of the St. Elias Mountains within the Donjek Range and is approximately 70 km long and 3–4 km wide. Our 2018 drill site was located on the upper northern arm of the glacier in the accumulation zone (60.78${}^{\circ }$ N, 139.63${}^{\circ }$ W), at an elevation of 2640 m a.s.l. Based on satellite imagery, Foy et al. (2011) estimated an average equilibrium line altitude (ELA) for the glacier of 1958 m a.s.l. for the period 1977–2007, while Young et al. (2020) provided a mean ELA of $\mathrm{2261}\pm \mathrm{151}$ m a.s.l. for the years 2013–2019. Our core site is thus well above the ELA and has remained within the main accumulation area of the glacier. Mean annual and summer (JJA) air temperatures from 1979–2019 were $-\mathrm{10.7}$ and $-\mathrm{2.5}$ ${}^{\circ }$C, respectively, based on bias-adjusted ERA5 climate reanalyses (Hersbach et al., 2020). The main melt season occurs from June through August. Over the period 1979–2016, Williamson et al. (2020) reported that the St. Elias ice field air temperature warmed at an average rate of 0.19 ${}^{\circ }$C per decade at an elevation of 2000–2500 m a.s.l., rising to 0.28 ${}^{\circ }$C per decade at an elevation of 5500–6000 m a.s.l.

Previous studies of Kaskawulsh Glacier have included an analysis of volume change over time based on comparisons of satellite imagery and digital elevation models (Foy et al., 2011; Young et al., 2020). Several reports in the 1960s documented various glaciological characteristics and processes occurring in the St. Elias ice fields as part of the Icefield Ranges Research Project (IRRP) (Wood, 1963; Grew and Mellor, 1966; Marcus and Ragle, 1970). Firn density and temperature measurements to 15 m depth were made during this period at site IRRP A, near the Kaskawulsh–Hubbard divide and about 5 km from our core site. Additional snow accumulation data are available from the “Copland Camp” site on the upper Hubbard Glacier, located $\sim \mathrm{12}$ km southwest of our drill site and at a similar elevation (Fig. 1). A weather station located on a nunatak near Copland Camp has been in service since 2013 (60.70${}^{\circ }$ N, 139.80${}^{\circ }$ W; $\sim \mathrm{2600}$ m a.s.l.; Fig. 1). Other relevant studies in the region include ice cores collected from the Eclipse ice field, located 12 km northwest of our drill site (Yalcin et al., 2006; Zdanowicz et al., 2014) but at a higher elevation (3017 m a.s.l.).

We consider snow accumulation rates, weather conditions, and earlier firn core studies across several different locations within this broad snowfield region that constitutes the upper accumulation areas of the Kaskawulsh and Hubbard glaciers. Some caution is needed in comparing different sites, but the region is relatively flat and uniform, with the exception of some nunataks. Away from the nunataks there is negligible influence from topographic obstacles or valley walls, so we hypothesize that the upper accumulation area will be exposed to similar climate conditions and snow accumulation rates over long periods. The possibility of significant spatial variability cannot be ruled out, however, so we consider this further in the data analysis.

3 Methods

3.1 Ice core field collection

Two 8 cm diameter cores were drilled between 20 and 24 May 2018, using an ECLIPSE ice drill (Icefield Instruments, Whitehorse, Yukon). With a starting depth of 2 m below the snow surface, Core 1 was 34.6 m long and reached a depth of 36.6 m, and Core 2 was 19.6 m long and reached a depth of 21.6 m. The two cores were drilled 60 cm apart, and core stratigraphy and density were recorded in the field. At a depth of 34.5 m below the snow surface, liquid water became evident in Core 1; drilling was stopped at a depth of 36.6 m to avoid the risk of the drill freezing in the hole.

Once the cores were retrieved, the presence of ice layers, ice lenses, and “melt-affected” firn was logged, and the stratigraphic character (e.g., texture, opacity), depth, and thickness were recorded. Melt-affected firn refers to any firn that displays physical characteristics indicating that there was the presence of liquid water at some point (Fig. S5). This can result in ice layers or ice lenses or can be indicated by the lack of grain boundaries, the presence of air bubbles, and opacity. When an ice horizon extended across the entire diameter of the core, it was labelled as an ice layer. If the ice horizon was of more limited lateral extent, it was labelled an ice lens. Ice lenses were occasionally wedge-shaped.

All of the density measurements for Core 1 were completed in the field. The Core 2 samples could not be measured for density in the field due to lack of time and were flown to Kluane Lake Research Station frozen, where the measurements were made within 24 h of arrival. A random assortment of 125 out of the 196 Core 2 sample bags were damaged during this transport and were not included in the measurements. This left a random reduced set of samples available to use for the density analysis, so we were able to construct a density stratigraphy (Core 2 in Fig. 2b), but uncertainties are higher for Core 2, and most of our analysis therefore focuses on Core 1.

Figure 2

Measured firn densities of (a) Core 1 and (b) Core 2 (20–24 May 2018), with uncertainties and best-fit logarithmic curves (black line). The depth scales are truncated at the location of the last summer surface at 4.2 m depth as the profile consisted of seasonal snow above this.

[Figure omitted. See PDF]

3.2 Ice core density analysis

Ice core density measurements were completed in the field. Each core was sawed into $\sim \mathrm{10}$ cm long sections and the diameter of the sections measured at each end. The sections were then double-bagged; weighed; and assessed for the quality of the core sample and its cylindrical completeness, which we denote $f$. The average diameter was used to determine the volume of the core section ($V$). Together with the mass of the core section, $m$, density was calculated following

1$$\mathit{\rho }=m/V,\text{ with }V=f\mathit{\pi }L(D/\mathrm{2}{)}^{\mathrm{2}},$$ where $\mathit{\rho }$ is the density of the firn, $D$ is the average core section diameter, $L$ is the length of the section, and $f\in$[0,1] is the subjectively assessed fraction of completeness of the core section. For example, if visual inspection indicated that about 5 % of the core was missing (e.g., due to missing ice chips caused by the core dogs of the drill head), then $f$ would be 0.95. Outliers were removed for the background firn density calculations if they were not physically possible (i.e., values $>\mathrm{917}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ or $<\mathrm{300}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ at depths below the last summer surface). Outliers from 32–36 m depth had residual liquid water in them, so these higher density values were retained.

In order to calculate the uncertainty in density, $\mathrm{d}\mathit{\rho }$, random and systematic sources of error have to be taken into account in the propagation of errors: 2$$\mathrm{d}\mathit{\rho }=\mathit{\rho }\sqrt{{\left({\displaystyle \frac{\mathrm{d}m}{m}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}V}{V}}\right)}^{\mathrm{2}}}.$$ The mass uncertainty was assumed to be 0.3 g, which is a conservative estimate given the scale's accuracy ($\pm \mathrm{0.1}$ g) but accounts for potential residual snow or water on the scale. The volume uncertainty is calculated by breaking down Eq. (1) for sample volume, $V=fAL$, where cross-sectional area $A=\mathit{\pi }(D/\mathrm{2}{)}^{\mathrm{2}}$. There is uncertainty in the measured length of the core section, $L$; the radius of the core section, $D/\mathrm{2}$; and the assessment of the completeness of the core sample, $f$. Each of these was calculated independently; and the propagation of uncertainty was calculated from 3$$\mathrm{d}V=V\sqrt{{\left({\displaystyle \frac{\mathrm{d}f}{f}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}A}{A}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}L}{L}}\right)}^{\mathrm{2}}};$$ $\mathrm{d}L$ was assumed to be 0.25 cm because the tape measure had ticks at every millimetre so it could be measured with precision, but core sections were often uneven, with crumbly edges caused by the drill cutters. The same uncertainty was assigned to the measurement of core diameter. Given two independent measurements, the uncertainty in the diameter is $\mathrm{d}D=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{(\mathrm{0.25}{)}^{\mathrm{2}}+(\mathrm{0.25}{)}^{\mathrm{2}}}=\mathrm{0.18}$ cm. For the cross-sectional area, the uncertainty $\mathrm{d}A=\mathit{\pi }D\mathrm{d}D/\mathrm{2}$.

Values of $f$ were determined by assessing the shape of the core and deciding how complete a cylinder the core section represented (e.g., accounting for missing volume due to chips from the core dogs along the edges). Three different people performed this evaluation, so there was subjectivity in each of the $f$ values, and it is best to be conservative with this estimate. We assigned this to be $\mathrm{d}f=\mathrm{0.2}$ for $f<\mathrm{0.8}$ and $\mathrm{d}f=\mathrm{0.1}$ for $f\ge \mathrm{0.8}$. The uncertainty in a higher $f$ value is lower because when a core was of good quality it was obvious. Less complete cylinders were more difficult to assess, hence the greater uncertainty when $f\le \mathrm{0.8}$. The $f$ value has the greatest effect on the overall uncertainty calculation for firn density. We did not record $f$ values for Core 2 in the field, so values are based on the measurements from Core 1. The minimum value recorded in Core 1 was $f=\mathrm{0.7}$, with a maximum of 1 and an average of 0.96. We assume a value of $f=\mathrm{0.96}\pm \mathrm{0.1}$ for all of Core 2.

The resulting uncertainty in the density was calculated from 4$$\mathrm{d}\mathit{\rho }=\mathit{\rho }\sqrt{{\left({\displaystyle \frac{\mathrm{d}m}{m}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}V}{V}}\right)}^{\mathrm{2}}}.$$ For the average densities, $\overline{\mathit{\rho }}$, the uncertainty can be calculated from the standard error of the mean, $\mathrm{d}\overline{\mathit{\rho }}=\mathrm{d}\mathit{\rho }/\sqrt{N}$, for sample size $N$. This can be estimated from the average value of $\mathrm{d}\mathit{\rho }$, but we report the more precise uncertainty calculated from the root mean square value of all point values, $\mathrm{d}{\mathit{\rho }}_k$: $\mathrm{d}\overline{\mathit{\rho }}=\frac{\mathrm{1}}{N}{\left[\sum _N\mathrm{d}{\mathit{\rho }}_k^{\mathrm{2}}\right]}^{\mathrm{1}/\mathrm{2}}$. Density can be expressed as water equivalence (w.e.) for each core section from the conversion $w=L\mathit{\rho }/{\mathit{\rho }}_{\mathrm{w}}$, where ${\mathit{\rho }}_{\mathrm{w}}$ is the density of water. For the whole core, of length $L_{\mathrm{c}}$, the water equivalence is $w_{\mathrm{c}}=L_{\mathrm{c}}\overline{\mathit{\rho }}/{\mathit{\rho }}_{\mathrm{w}}$, with units m w.e. We also include an estimate of the age of the cores, based on an estimate of the average annual net accumulation rate, $\overline{a}$, with units m w.e. yr${}^{-\mathrm{1}}$. The age of the core is then ${\mathit{\tau }}_{\mathrm{c}}=w_{\mathrm{c}}/\overline{a}$. Uncertainty is estimated by propagation of uncertainties in $w_{\mathrm{c}}$ and $\overline{a}$. We use an uncertainty of $\mathrm{d}{L}_{\mathrm{c}}=\mathrm{0.5}$ m for the total length of the core, $L_{\mathrm{c}}$, which is based on measurements during retrieval of Core 1 of 35.05 m from the drill panel, 34.59 m from the addition of core lengths, and 34.25 m from the sum of the $\sim \mathrm{10}$ cm samples. For Core 2 the length was 19.75 m from the drill panel, 19.35 m from the addition of core lengths, and 19.63 m from the sum of the $\sim \mathrm{10}$ cm samples.

Ice fraction, $F_{\mathrm{i}}\in [\mathrm{0},\mathrm{1}]$, was calculated for each 10 cm section of the firn core. Here ice was defined based on its lack of air bubbles and crystalline structure, as compared to the granular structure of firn. We refer to this as ice fraction rather than melt percent as melt percent generally assumes that the meltwater remains within the net annual accumulation layer (Koerner, 1977), which cannot be assumed here due to evidence that meltwater percolates beyond the annual accumulation layer and refreezes into previous years' accumulation. The thickness of individual ice layers was summed within each 10 cm core section. In core samples that had ice lenses, their diameter typically occupied about 50 % of the core sample; therefore their thickness was divided by 2 before being summed. For each core section, total ice content was divided by the length of the section, $L$, to give $F_{\mathrm{i}}$. These values were also summed to give the total ice core ice content.

To understand the firn densification process in the absence of refrozen meltwater, the “background” firn density is of interest. For each sample, we estimated this by subtracting the mass and volume of the ice to give the firn density in the absence of ice content. We used a 30 cm moving average of total ice content and density in order to smooth out a possible error of $\pm \mathrm{10}$ cm in assigning the location of the ice features within the stratigraphy. Each sample had a measured bulk density, ${\mathit{\rho }}_{\mathrm{b}}$, which we assume resulted from a binary mixture of ice and firn, with densities ${\mathit{\rho }}_{\mathrm{f}}$ and ${\mathit{\rho }}_{\mathrm{i}}$. Ice and firn fractions, $F_{\mathrm{i}}$ and $F_{\mathrm{f}}$, were defined with $F_{\mathrm{i}}+F_{\mathrm{f}}=\mathrm{1}$. The background firn density was then calculated following 5$${\mathit{\rho }}_{\mathrm{f}}=({\mathit{\rho }}_{\mathrm{b}}-{\mathit{\rho }}_{\mathrm{i}}{F}_{\mathrm{i}})/F_{\mathrm{f}}.$$ In cases where there was no ice fraction ($F_{\mathrm{i}}=\mathrm{0}$), ${\mathit{\rho }}_{\mathrm{f}}={\mathit{\rho }}_{\mathrm{b}}$. Ice layers and lenses were assumed to have a density of $\mathrm{874}\pm \mathrm{35}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, based on the average density of firn core sections that were 100 % ice in Greenland (873 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) and Devon ice cap (875 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) (Bezeau et al., 2013; Machguth et al., 2016). This is different from the 917 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ upper bound used in the outlier analysis because that is the theoretical limit for pure ice, whereas 874 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ is based on measured field data, which include observed ice layers and lenses which have small bubbles and imperfections in them.

There is surface lowering associated with melting but without associated mass loss due to subsurface refreezing. This surface lowering is an “apparent ablation” in airborne or satellite altimetry signals. We calculated this for each core section using the background firn density, ${\mathit{\rho }}_{\mathrm{f}}$, and length of the section, $L$. The “thinning” or surface lowering of a given core section, $\mathrm{\Delta }L$, was estimated by reverting the ice to the density of the background firn following 6$$\mathrm{\Delta }L=L\left[\left(F_{\mathrm{f}}+{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}{F}_{\mathrm{i}}}{\mathit{\rho }}}}_{\mathrm{f}}\right)-\mathrm{1}\right].$$ Summed over the full firn column, this gives the total surface lowering associated with meltwater that percolates and refreezes, with no actual loss of mass.

3.3 Historical measurements

As part of an expedition undertaken by the IRRP, Grew and Mellor (1966) measured snow density and temperature to a depth of 15 m at the IRRP A site on 23 July 1964 (Fig. 1). The first $\sim \mathrm{4}$ m was measured in a snow pit, while the remaining $\sim \mathrm{11}$ m was based on measurements of a core drilled with a Cold Regions Research and Engineering Laboratory (CRREL) coring auger. The original data are not available, so values were reconstructed based on digitization of the density plot provided in Fig. 4 of Grew and Mellor (1966). This digitization was undertaken with WebPlot Digitizer 4.3 (Rohatgi, 2020) and has an estimated error of $\pm \mathrm{2}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for density and $\pm \mathrm{0.01}$ m for depth. Errors were calculated by clicking the same point 25 times and evaluating the variability in the points (i.e., the standard deviation).

From 14–17 July 2006, snow density and temperature measurements were recorded every 10 cm to a depth of 10.4 m at the Copland Camp as part of a University of Ottawa field class. Measurements from 0 to 5.4 m were recorded in a snow pit, while those from 5.5 to 10.4 m were based on a core recovered with a Kovacs Mark II coring system (Kovacs Enterprises, Oregon, USA). Density measurements in the snow pit were undertaken with a 250 cm${}^{\mathrm{3}}$ RIP 2 cutter (Snowmetrics, Colorado, USA) and in the ice core by measuring and weighing core sections and using Eq. (1). Errors in the density measurements were determined from Eq. (4) and verified against density values recorded in a second snow pit dug to a depth of 4.0 m, approximately 2 m away from the first. All temperature measurements were undertaken with a Thermor PS100 digital stem thermometer with an accuracy of $\pm \mathrm{0.5}$ ${}^{\circ }$C.

Annual snow accumulation at the Copland Camp was measured between 2004–2011 with a Campbell Scientific SR50 sonic ranging sensor mounted on a cross-arm on a vertical steel pole drilled into the firn. The SR50 was connected to a Campbell Scientific CR10X logger and included a correction for the change in speed of sound with air temperature. The mounting pole was raised annually to keep it above the snow surface, and densities recorded in snow pits collected during annual University of Ottawa field classes (typically in early July) were used to convert the SR50 depth measurements into w.e. values.

3.4 Energy balance and firn modelling

ERA climate reanalyses were used to examine changes in climate and annual surface melting at the study site since the 1960s, coupled with a firn model to simulate the decadal evolution of firn temperature, hydrology, ice content, and density. Daily melt rates were calculated from 1965 to 2019 using a surface energy balance model (Ebrahimi and Marshall, 2016), coupled to a subsurface model of coupled thermal and hydrological evolution in the snow and firn (Samimi et al., 2020). The model calculates the surface energy budget and snowmelt based on incoming shortwave and longwave radiation, temperature, relative humidity, wind speed, and air pressure, with internal parameterizations of surface albedo evolution and outgoing longwave radiation. Conductive heat flux to the snow surface and snow surface temperatures are simulated within the subsurface snow and firn model. Snow and firn densification are parameterized following Vionnet el al. (2012) for the firn matrix (“background firn”), with bulk density including the additional mass of any ice or water content. Details of the model are provided in the Supplement.

Meteorological inputs for the surface energy balance model were derived from the ERA5 climate reanalysis for the period 1979 to 2019 (Hersbach et al., 2020) and extended back to 1965 using the ERA 20th-century reanalysis (ERA20c; Poli et al., 2016). ERA5 outputs are at a resolution of 0.25${}^{\circ }$ latitude and longitude, and data for our analysis were averaged from ERA5 grid cells located at 60.75${}^{\circ }$ N, 139.75${}^{\circ }$ W, and 60.75${}^{\circ }$ N, 139.5${}^{\circ }$ W. ERA20c data are at 1${}^{\circ }$ latitude and longitude resolution, and we interpolated meteorological conditions to the upper Kaskawulsh Glacier from the four model grid cells at 60 to 61${}^{\circ }$ N and 139 to 140${}^{\circ }$ W. ERA20c fields were homogenized with ERA5 through bias adjustments for 2 years of overlap in the reanalyses, 1979 and 1980, with ERA5 assumed to be the more accurate reconstruction. Monthly bias adjustments based on this period of overlap were then applied to the ERA20c data from 1965 to 1978.

The reanalysis data represent the climatology over the region of the upper Kaskawulsh–Hubbard divide (i.e., a 0.25${}^{\circ }$ grid cell) and are not specific to our core site. The firn modelling is therefore taken to be generally applicable for this upper plateau region. ERA air pressure and 2 m temperature and humidity fields were bias-adjusted to the specific elevation of our core site, 2640 m (see the Supplement). ERA5 temperature fields were evaluated against Copland weather station data from 2014–2018, which indicate a small (0.6 ${}^{\circ }$C) cold bias in the ERA5 data for average summer (JJA) temperatures over this period. ERA temperatures were further bias-adjusted by this amount. Our firn core site, the Copland weather station, Copland Camp, and IRRP research sites all fall within the same ERA5 grid cell, and we make the assumption that climate conditions are similar for similar elevations and glaciological settings within this region.

Surface energy balance and melt were calculated every 30 min using mean daily meteorological forcing from ERA and a parameterization of the diurnal cycles of temperature and incoming shortwave radiation (Ebrahimi and Marshall, 2016). Snow accumulation is based on the ERA total precipitation, with a constant scaling factor of 1.6 in order to give representative annual totals. With this scaling, mean annual precipitation at the site ($\pm \mathrm{1}$ standard deviation) was $\mathrm{1.83}\pm \mathrm{0.32}$ m w.e. yr${}^{-\mathrm{1}}$. Snow is updated monthly in the numerical simulation. We neglect rainfall as we do not have a good constraint on how much summer precipitation falls as rain, and this will not be reliably predicted in the climate reanalysis. Summer temperatures are cool (mean value of $-\mathrm{2.4}$ ${}^{\circ }$C), and during our experience while working at the Copland Camp in the month of July, we have experienced numerous snow events but no rainfall. While rain must occur from time to time, we believe it to be rare at the study site.

Subsurface temperatures were modelled for a 35 m firn column, with a simple model for meltwater percolation that accounts for meltwater refreezing and the associated latent heat release where snow or firn is below 0 ${}^{\circ }$C (Samimi and Marshall, 2017; Samimi et al., 2020). For the current study, we discretize the snow and firn into 58 layers from 0.1 to 1 m in thickness, with higher resolution near the surface. The firn model is coupled with the surface energy balance model, solving for the firn thermodynamic and hydrological evolution at 30 min time steps for the period 1965 to 2019. The subsurface temperature evolution includes vertical heat conduction and latent heat release from refreezing. Heat advection associated with snow accumulation is neglected. When subsurface temperatures reach 0 ${}^{\circ }$C, liquid water is retained or percolates to depth, following a Darcian parameterization for water flux: $q_{\mathrm{w}}=-k_{\mathrm{h}}\mathrm{\nabla }\mathit{\phi }$, for hydraulic conductivity $k_{\mathrm{h}}$ and hydraulic potential $\mathit{\phi }$ (Samimi and Marshall, 2017). For the numerical experiments in this study we set $k_{\mathrm{h}}={\mathrm{10}}^{-\mathrm{5}}$ m s${}^{-\mathrm{1}}$ in snow and 10${}^{-\mathrm{6}}$ m s${}^{-\mathrm{1}}$ for snow and firn, respectively. Capillary water retention is calculated following Coléou and Lesaffre (1998). The default model parameters are based on calibration at DYE-2, Greenland, in the percolation zone of the southern Greenland Ice Sheet (Samimi et al., 2020). A broader range of model parameters are explored in sensitivity analyses presented in the Supplement.

The model is “spun up” through a 30-year simulation with perpetual 1965 climate forcing (i.e., running through 30 annual cycles with 1965 climate conditions). This provides the initial temperature, density, and ice layer structure within the firn column. Ideally, a spin-up simulation forced by the historical meteorological conditions (e.g., the years 1935–1965) would be preferable to assuming a perpetual climatology from a single year. Mean annual and mean summer temperatures were close to standard climatological values in 1965, so the model results are not strongly sensitive to this assumption (discussed in the results), but we explore a range of numerical experiments to examine the model sensitivity to these initial conditions and the spin-up assumptions.

Table 1

Total ice content, ice fraction ($F_{\mathrm{i}}$), bulk density (${\mathit{\rho }}_{\mathrm{b}}$), background density (${\mathit{\rho }}_{\mathrm{f}}$), and total water equivalent ($w$) for the portion of each core below the last summer surface. Depths are reported from the May 2018 snow surface and the firn portion of the core started at the 2017 summer surface, at 4.2 m depth.

	Depth below	Total Ice	$F_{\mathrm{i}}$	$F_{\mathrm{i}}$	${\mathit{\rho }}_{\mathrm{b}}$	${\mathit{\rho }}_{\mathrm{f}}$	$w$
	surface (m)	content (m)	(% vol)	(% mass)	($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$)	($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$)	w.e. (m)
Core 1	4.2–14.2	0.67 $\pm$ 0.07	6.7 $\pm$ 0.7	13.0 $\pm$ 1.3	588 $\pm$ 8	565 $\pm$ 9	5.88 $\pm$ 0.08
	4.2–21.6	1.51 $\pm$ 0.15	8.7 $\pm$ 0.9	15.6 $\pm$ 1.6	640 $\pm$ 6	613 $\pm$ 7	11.08 $\pm$ 0.11
	4.2–36.6	2.33 $\pm$ 0.26	7.2 $\pm$ 0.7	11.9 $\pm$ 1.2	698 $\pm$ 5	676 $\pm$ 6	22.49 $\pm$ 0.15
Core 2	4.2–14.2	0.42 $\pm$ 0.04	4.2 $\pm$ 0.4	8.4 $\pm$ 0.8	572 $\pm$ 7	556 $\pm$ 7	5.72 $\pm$ 0.07
	4.2–21.6	0.81 $\pm$ 0.08	4.7 $\pm$ 0.5	8.5 $\pm$ 0.9	624 $\pm$ 5	609 $\pm$ 6	10.85 $\pm$ 0.09
Average	4.2–14.2	1.18	4.0	564	580 $\pm$ 5	560 $\pm$ 5	5.80 $\pm$ 0.05
	4.2–21.6				632 $\pm$ 4	611 $\pm$ 4	10.97 $\pm$ 0.07

4 Results

4.1 Ice core density

The density data are plotted in Fig. 2, fitted with a logarithmic curve to quantitatively compare our two cores. The first 4.2 m of both 2018 cores was dry and had an average density of $\mathrm{450}\pm \mathrm{21}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, with no ice content. At 4.2 m there was a significant ice crust, with large crystal size, rounded grains, and high impurity content, which was assumed to represent the last summer surface (LSS) from 2017. The snow above this LSS layer was therefore classified as seasonal snow. In this section we focus on the firn characteristics below the LSS, so our discussion is centred on the core recovered between 4.2 and 36.6 m below the surface for Core 1 (i.e., total firn length of 32.4 m) and between 4.2 and 21.6 m below the surface for Core 2 (i.e., total firn length of 17.4 m). For consistency, we reference all depths to the seasonal snow surface throughout this paper.

In the upper 10 m of firn (4.2 to 14.2 m below the surface; Table 1), Cores 1 and 2 had average densities of $\mathrm{588}\pm \mathrm{8}$ and $\mathrm{572}\pm \mathrm{7}$ kg m${}^{-\mathrm{3}}$, respectively, giving an overall average density of $\mathrm{580}\pm \mathrm{5}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Over the upper 17.4 m of firn in each core (4.2 to 21.6 m below the surface, the depth to the bottom of Core 2), Kaskawulsh firn had an average density of $\mathrm{632}\pm \mathrm{4}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. The full 32.4 m of firn at Core 1 (4.2 to 36.6 m below the surface) had an average density of $\mathrm{698}\pm \mathrm{5}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Ice content generally increased with depth in the upper $\sim \mathrm{25}$ m of the core, but deeper sections were less icy (Table 1). The bottom 5 m of firn in Core 1 had an average density of $\mathrm{826}\pm \mathrm{13}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ but with no identified ice layers. Based on the high density and texture of this deep firn, along with the presence of liquid water in the deepest sections of the core, we believe that we drilled to near the base of the firn at the core site but cannot confirm this as we halted drilling before reaching glacier ice.

Total ice content in the 32.2 m firn portion of Core 1 (4.2 to 36.6 m below the surface) was $\mathrm{2.33}\pm \mathrm{0.26}$ m of ice, or $\mathrm{2.67}\pm \mathrm{0.24}$ m w.e. This is equivalent to 7.2 % by volume and 11.9 % by mass (Table 1). Using Eq. (5) and the values for ice content in Core 1, we estimate a background firn density of $\mathrm{676}\pm \mathrm{6}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for the full column of firn, 3.2 % less than the bulk density of the firn (Table 1). The two cores had very similar bulk and background densities over the upper 10 m of firn (4.2 to 14.2 m below the surface) and 17.4 m (4.2 to 21.6 m below the surface), where a direct comparison was possible. The total water equivalent of firn in Core 1 was calculated to be $w_{\mathrm{c}}=\mathrm{22.5}\pm \mathrm{0.2}$ m w.e.

Figure 3

Stratigraphy of the cores collected in May 2018. LSS is the last summer surface at 4.2 m, the boundary between seasonal snow above and firn below. Ice layer thicknesses were classified in the legend by thickness distribution. Note that the ice layers in the first several metres of the core are interpreted as wind crusts.

[Figure omitted. See PDF]

4.2 Ice core stratigraphy

The stratigraphy of the 2018 cores indicates numerous ice layers as well as melt-affected firn, distinguished by a lack of grain boundaries or opaque, bubbly firn. The first 4.2 m comprised the seasonal snowpack, with firn below. Within the first 6 m below the surface, there were several small ice layers ($<\mathrm{2.5}$ cm thick), interpreted as wind crusts (Fig. 3). Several thick ($>\mathrm{10}$ cm) ice layers were found between 6 and 26 m depth (1.8 to 21.8 m in the firn). The largest ice layer in Core 1 was 22 cm thick, found at 14.1 m (9.9 m in the firn). At 26.4 m (22.2 m in the firn) the ice layers and lenses disappeared. Below this the firn was almost entirely meltwater-affected, based on its appearance and texture, but without the quantity of ice lenses or ice layers that were present in the first 25 m. We interpret this section of the core as infiltration ice, consisting of water-saturated firn that has experienced refreezing. At 30 m depth (25.8 m in the firn), the meltwater effects were absent, and there were two small ice layers and an ice lens. At 30.6 m depth the firn was melt-affected again. From 34.5 to 36.6 m (30.3 to 32.4 m in firn) the core sections expelled liquid water as they were extracted from the core barrel. The expulsion of liquid water can be viewed in Video Supplement 1 and 2 (10.5446/50918 and 10.5446/50919, Ochwat, 2021a, b).

In Core 2 there were numerous ice layers starting at a depth of 3.8 m, and below 4.4 m (0.2 m in the firn) the core was meltwater-affected. There was a thick ice layer at 6.6 m (2.4 m in the firn) that was 30 cm lower than a similar ice layer in Core 1 at 6.3 m. There were numerous melt-affected layers between ice lenses much closer to the surface in Core 2 than Core 1. In Core 1 there were several ice layers at $\sim \mathrm{10}$ m depth (5.8 m in the firn), but these layers were not present in Core 2. At 14.4 m (10.2 m in the firn) another section of the firn had numerous ice layers ($\sim \mathrm{20}$–30 cm deeper than recorded in Core 1), and at 14.6 m the thickest ice layer was encountered (12 cm), corresponding well with the thickest layer in Core 1. Between 16 and 21.5 m (11.8 to 17.3 m in the firn) the core was melt-affected. We attribute differences between Core 1 and Core 2 stratigraphy to uncertainty in the depth of features (as discussed in Sect. 3.2) and horizontal variability in meltwater infiltration, which is known to occur at length scales less than 1 m (Parry et al., 2007; Harper et al., 2011). Spatial heterogeneity in firn is common in areas with high surface melt due to differential melting and percolation that is complex due to the presence of sastrugi and wind crusts, different permeability of the snow and firn, and vertical piping mechanisms (Marchenko et al., 2017; Parry et al., 2007).

Figure 4

(a) Comparison of densities averaged over 1 m segments at IRRP A on 23 July 1964 (Grew and Mellor, 1966; blue), at Copland Camp on 14–17 July 2006 (red), and at Core 1 on 20–24 May 2018 (green). Depth of LSS (i.e., boundary between seasonal snow above and firn below) was 3.28 m in 1964, 3.50 m in 2006, and 4.22 m in 2018; the density data for 2018 begin at the LSS due to the difference in time of year of the measurements compared to the others. (b) Comparison between cumulative w.e. content in the 1964, 2006, and 2018 profiles, starting at the LSS.

[Figure omitted. See PDF]

4.3 Changes in firn characteristics over time

The firn in the accumulation area of Kaskawulsh Glacier has become denser since 1964 (Fig. 4a). The mean density of the upper 7 m of firn was 516 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in 1964 (3.3 to 10.3 m below the surface), 590 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in 2006 (3.5 to 10.5 m below the surface), and 549 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in 2018 (4.4 to 11.4 m below the surface). The difference between the average densities from the upper 7 m of the 1964 and 2018 core is 33 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, which is an increase of $\sim \mathrm{7}$ %. It is difficult to assess whether firn temperatures have changed over this time as limited data are available from below the depth of the annual temperature wave ($\sim \mathrm{10}$ m for heat diffusion and deeper than this with the effects of subsurface meltwater infiltration and latent heat release). Borehole temperature records from Grew and Mellor (1966) indicate temperate (0 ${}^{\circ }$C) conditions at 15 m depth in the summer of 1964, which suggests that deep temperate firn may have existed at this site in the 1960s. This supports the assumption that Kaskawulsh Glacier is temperate (Foy et al., 2011), despite mean annual air temperatures of about $-\mathrm{11}$ ${}^{\circ }$C on the upper glacier.

Accumulation data from the IRRP A site, Copland Camp, and our 2018 measurements do not show any evidence for a significant change over time, although there can be high interannual variability. At IRRP A, Wagner (1969) reported values between 1.3 and 1.9 m w.e. yr${}^{-\mathrm{1}}$ for 1963. Marcus and Ragle (1970) measured a winter snow accumulation of 1.6 m w.e. from 1964–1965. Holdsworth (1965) reported an estimated mean annual accumulation rate of 1.8 m w.e. yr${}^{-\mathrm{1}}$ in the early 1960s (year not specified) (Holdsworth, 1965). Yearly snow accumulation data from 2004–2011 collected with the SR50 at Copland Camp indicate a mean annual accumulation rate of 1.77 m w.e. yr${}^{-\mathrm{1}}$, with variations between 1.3 and 2.4 m w.e. yr${}^{-\mathrm{1}}$. The seasonal snowpack at our drill site was 4.2 m in May 2018, with an average snow density of 440 kg m${}^{-\mathrm{1}}$, giving a total accumulation of 1.85 m w.e. for 2017–2018.

Based on the above review, we adopt an estimate of $\overline{a}=\mathrm{1.8}\pm \mathrm{0.2}$ m w.e. yr${}^{-\mathrm{1}}$ for the net accumulation from 2005 to 2018. Using this value, the firn layer of Core 1 represents $\mathrm{12.5}\pm \mathrm{1.4}$ years of net accumulation (i.e., 2005–2017) or $\mathrm{13.2}\pm \mathrm{1.4}$ years if the seasonal snowpack on top is counted. Over 12.5 years, the total measured ice content of 2.67 m w.e. in the firn equates to an average meltwater refreezing rate of 0.22 m w.e. yr${}^{-\mathrm{1}}$.

Table 2

Climate, surface energy balance, and firn conditions, 1965 to 2019, based on the ERA meteorological forcing at the core site. Decadal trends are reported from linear fits to the data. The period 1965–1975 represents the historical baseline period, when much of the work of the IRRP was completed. The years 2005–2017 represent the period of record of Core 1, and 2013 was an exceptional year, which potentially marked the initial development of the firn aquifer at this site. Melt and refreeze refer to the total annual melting and refreezing in the 35 m snow and firn column, “drainage” is the total annual melt minus refreezing, and “net melt” is the surface mass loss or drawdown associated with summer melting. Freeze–thaw cycles in the surface layer of the snow mean that some fraction of the net energy that is available for melt is directed to refrozen (i.e., recycled) meltwater. As a result, the net melt – meltwater that is available to percolate into the deeper snowpack and firn pack – is less than the total summer melt; note that rainfall is neglected in this study and is assumed to be negligible.

	1965–2019	Trend	1965–1975	2005–2017	2013
	Mean ($\pm \mathrm{1}\mathit{\sigma }$)	(decade${}^{-\mathrm{1}})$	Mean ($\pm \mathrm{1}\mathit{\sigma }$)	Mean ($\pm \mathrm{1}\mathit{\sigma }$)
Meteorological conditions
$T_{\mathrm{ann}}$ (${}^{\circ }$C)	$-\mathrm{10.7}\pm \mathrm{0.9}$	$+\mathrm{0.16}$	$-\mathrm{11.2}\pm \mathrm{0.7}$	$-\mathrm{10.4}\pm \mathrm{1.0}$	$-\mathrm{9.6}$
$T_{\mathrm{JJA}}$ (${}^{\circ }$C)	$-\mathrm{2.4}\pm \mathrm{0.8}$	$+\mathrm{0.07}$	$-\mathrm{2.2}\pm \mathrm{0.9}$	$-\mathrm{2.1}\pm \mathrm{0.7}$	$-\mathrm{0.7}$
$T_{\mathrm{SJJA}}$ (${}^{\circ }$C)	$-\mathrm{2.3}\pm \mathrm{0.8}$	$+\mathrm{0.29}$	$-\mathrm{2.9}\pm \mathrm{0.9}$	$-\mathrm{1.8}\pm \mathrm{0.6}$	$-\mathrm{0.8}$
Positive degree days (${}^{\circ }$C d)	54 $\pm$ 23	$+\mathrm{3.60}$	49 $\pm$ 16	69 $\pm$ 31	123
$q_{\mathrm{v}}$ (g kg${}^{-\mathrm{1}}$)	3.7 $\pm$ 0.2	$+\mathrm{0.10}$	3.5 $\pm$ 0.2	3.9 $\pm$ 0.2	4.2
Surface energy balance (JJA values)
$Q*$ (W m${}^{-\mathrm{2}}$)	18 $\pm$ 11	$+\mathrm{3.7}$	8 $\pm$ 3	26 $\pm$ 13	45
$Q_N$ (W m${}^{-\mathrm{2}}$)	10 $\pm$ 9	$+\mathrm{2.6}$	4 $\pm$ 3	16 $\pm$ 13	37
Net melt (mm w.e. yr${}^{-\mathrm{1}}$)	230 $\pm$ 210	$+\mathrm{62}$	100 $\pm$ 80	380 $\pm$ 310	895
Melt (mm w.e. yr${}^{-\mathrm{1}}$)	520 $\pm$ 270	$+\mathrm{81}$	360 $\pm$ 130	720 $\pm$ 375	1360
Refreeze (mm w.e. yr${}^{-\mathrm{1}}$)	500 $\pm$ 195	$+\mathrm{48}$	360 $\pm$ 130	615 $\pm$ 205	1100
Drainage (mm w.e. yr${}^{-\mathrm{1}}$)	20 $\pm$ 120	$+\mathrm{32}$	0 $\pm$ 0	105 $\pm$ 215	260
Firn conditions
$T_{\mathrm{1}}$ (${}^{\circ }$C)	$-\mathrm{12.8}\pm \mathrm{0.9}$	$+\mathrm{0.2}$	$-\mathrm{13.3}\pm \mathrm{0.8}$	$-\mathrm{12.4}\pm \mathrm{0.9}$	$-\mathrm{11.5}$
$T_{\mathrm{10}}$ (${}^{\circ }$C)	$-\mathrm{7.3}\pm \mathrm{3.4}$	$+\mathrm{1.8}$	$-\mathrm{11.3}\pm \mathrm{0.8}$	$-\mathrm{2.9}\pm \mathrm{2.4}$	$-\mathrm{3.0}$
$T_{\mathrm{20}}$ (${}^{\circ }$C)	$-\mathrm{7.2}\pm \mathrm{3.6}$	$+\mathrm{2.1}$	$-\mathrm{12.2}\pm \mathrm{0.5}$	$-\mathrm{3.7}\pm \mathrm{2.5}$	$-\mathrm{4.5}$
$T_{\mathrm{35}}$ (${}^{\circ }$C)	$-\mathrm{8.0}\pm \mathrm{3.5}$	$+\mathrm{2.1}$	$-\mathrm{12.7}\pm \mathrm{0.4}$	$-\mathrm{4.8}\pm \mathrm{1.9}$	$-\mathrm{5.2}$
$z_{\mathrm{thaw}}$ (m)	6.8 $\pm$ 9.4	$+\mathrm{3.6}$	1.2 $\pm$ 1.0	13.1 $\pm$ 12.5	18
$E_{\mathrm{lat}}$ (MJ m${}^{-\mathrm{2}}$)	126 $\pm$ 41	$+\mathrm{9.3}$	$\mathrm{98}\pm \mathrm{30}$	147 $\pm$ 43	258
${\mathit{\rho }}_{\mathrm{b}}$ ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$)	655 $\pm$ 10	$+\mathrm{4.6}$	645 $\pm$ 3	663 $\pm$ 12	671
Ice content (m)	2.0 $\pm$ 0.6	$+\mathrm{0.2}$	1.1 $\pm$ 0.3	2.3 $\pm$ 0.4	2.6

Figure 5

Modelled meteorological, surface mass balance, and firn conditions from 1965 to 2019: (a) summer (JJA) air and snow surface temperatures (${}^{\circ }$C); (b) annual melting and “drainage” (melting minus refreezing; mm w.e. yr${}^{-\mathrm{1}}$); (c) annual mean snow and firn temperature at the surface (0.1 m) and at depths of 10, 20, and 35 m (${}^{\circ }$C); (d) modelled maximum depths of the summer wetting and thawing fronts (m); (e) average firn density for the full firn column and in the upper 20 m ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$); (f) total firn ice content (metres).

[Figure omitted. See PDF]

4.4 Surface energy balance and firn modelling

Reconstructed air temperature, melt, and firn trends from 1965–2019 are shown in Fig. 5. Summer air temperature from the reanalysis (Fig. 5a) shows a modest but statistically significant increase over the study period, with a trend of $+\mathrm{0.07}$ ${}^{\circ }$C per decade. Table 2 reports changes in meteorological, energy balance, and modelled firn conditions over this time. Specific humidity and incoming longwave radiation increase markedly over the 55 years, with trends of $+\mathrm{0.1}$ g kg${}^{-\mathrm{1}}$ per decade and $+\mathrm{3.5}$ W m${}^{-\mathrm{2}}$ per decade, respectively. This echoes the findings of Williamson et al. (2020), who report decadal-scale, high-elevation warming in the St. Elias Mountains in association with increases in atmospheric water vapour and longwave radiation. These trends augment the net energy available for melt through increases in both the net radiation and latent heat flux. Modelled annual melt averaged $\mathrm{230}\pm \mathrm{210}$ mm w.e. yr${}^{-\mathrm{1}}$ from 1965 to 2019 and $\mathrm{380}\pm \mathrm{310}$ mm w.e. yr${}^{-\mathrm{1}}$ from 2005 to 2017, 70 % higher than the long-term average. The latter interval represents the approximate period of record of the firn core. The trend in surface melting is $+\mathrm{62}$ mm w.e. yr${}^{-\mathrm{1}}$ per decade from 1965 to 2019 (Fig. 5b). The summer of 2013 was exceptional; it had the warmest summer temperatures on record, $T_{\mathrm{JJA}}=-\mathrm{0.7}$ ${}^{\circ }$C, with 895 mm w.e. of meltwater (Table 2).

Within the model, 91 % of the surface meltwater refreezes in the firn over the period 1965–2019, with 100 % of it refreezing in cool summers, when meltwater generation is limited. Meltwater that does not refreeze percolates to depth in the firn model. Figure 5b plots the annual melting minus refreezing, with positive values indicating deep percolation. If the firn is temperate (0 ${}^{\circ }$C), meltwater can percolate through the entire depth of the firn column (35 m), where it is permitted to “drain” through the lowest layer; this water leaves the system and is considered to be run-off. Porewater in the firn can also refreeze in the subsequent winter, to the depth of the winter cold wave ($\sim \mathrm{8}$ m), accounting for the negative values in Fig. 5b. This is percolated meltwater that refreezes within the firn column in the following calendar year. Complete meltwater retention is typical for most of the period from 1965 to the early 2010s, but there is a marked increase in modelled run-off over the last decade (Fig. 5b), indicating drainage through the full 35 m firn column. Only 72 % of the surface melt refroze during the period 2005–2017, with a fivefold increase in meltwater drainage, from an average of $-\mathrm{20}\pm \mathrm{120}$ mm w.e. yr${}^{-\mathrm{1}}$ from 1965–2019 to $-\mathrm{105}\pm \mathrm{220}$ mm w.e. yr${}^{-\mathrm{1}}$ from 2005–2017. Meltwater that drains to the deep firn is portioned between porewater storage (meltwater retention) and deep drainage (mass loss). This partitioning was almost equal in the firn model over the period 2005–2017, with an average of 52 mm w.e. yr${}^{-\mathrm{1}}$ stored as liquid water in the deep firn and 54 mm w.e. yr${}^{-\mathrm{1}}$ of run-off: water that drains out through the bottom layer of the firn, leaving the system. This equates to a total meltwater retention of 86 % as either liquid water or refrozen ice, with a mass loss representing 14 % of the summer melt.

Figure 6

Modelled subsurface temperature evolution for the reference model climatology and parameter settings: (a) 1965–2019, full 35 m firn column; (b) 2000–2019, upper 20 m. Deep temperate conditions conducive to a firn aquifer developed from 2013 to 2017, in response to several subsequent summers of high melting and deep meltwater infiltration.

[Figure omitted. See PDF]

Summers with high amounts of surface melt produce greater refreezing and warming of the snow and firn, eventually overwhelming the content and enabling deep percolation and drainage. Figure 5c and d plot the modelled evolution of the firn temperatures and the wetting and melting fronts, which closely coincide. Snow and firn temperatures in Fig. 5c are mean annual values at the snow surface (the upper 0.1 m) and at 10, 20, and 35 m depth. For a purely conductive environment, $\sim \mathrm{10}$ m represents the depth of the annual temperature wave (Cuffey and Paterson, 2010), but latent heat release from meltwater refreezing warms the subsurface and causes a deeper influence of surface conditions, such that 10 m temperatures are highly variable (Table 2). The modelled wetting and melting fronts in Fig. 5d suggest dramatic recent developments in firn thermal and hydrological structure at the Kaskawulsh drill site, with a regime shift in the firn structure over the period 2013–2017. This is consistent with the birth of a deep PFA at this time. Figure 6 plots the full subsurface temperature evolution over the period 1965–2019, showing the typical seasonal evolution of firn temperatures and the unusual nature of the hydrological breakthrough event that began in 2013 and persists through 2019. Figure 5e and f plot the modelled increases in average firn density and total firn ice content from 1965–2019. The average firn density in the model is 682 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in 2018, compared with $\mathrm{698}\pm \mathrm{5}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ measured in Core 1. The modelled firn densification since the 1960s roughly matches the observed density trend.

The model results in Figs. 5 and 6 are for the “reference” 1965–2019 ERA climatological forcing and firn model parameters. These are the direct ERA climate fields, bias-adjusted to represent the elevation of the core site and to give consistency with the local Copland weather station data (2014–2018). The weather station has a similar elevation and topoclimatic environment and is about 11 km from the core site, falling within the same ERA5 grid cell. Firn model settings are based on calibrations against field data at DYE-2, Greenland, within the percolation zone of the southern Greenland Ice Sheet (Samimi et al., 2020), but we have no local field calibration of these model parameters. There are therefore uncertainties within both the climate forcing and the model parameters and assumptions. The Supplement examines the sensitivity of model results to several important meteorological inputs and model parameters as well as the strategy adopted for the model spin-up.

Figure 7

Sensitivity of the model simulations to (a, c) meteorological forcing and firn model parameters and (b, d) initial conditions, through different model spin-up settings. (a) Mean annual 20 m temperatures and (c) seasonal thaw depths from 1965–2019 for the reference model and for sensitivity experiments with $\pm \mathrm{1}$ ${}^{\circ }$C and for irreducible water contents of 0.02 ($\mathit{\theta }$w2) and 0.04 ($\mathit{\theta }$w4). The line colours in (a) also apply to (c). An extended set of sensitivity tests is presented in the Supplement. (b) Temperatures at 20 m and (d) thaw depths from 1965–2019 after a 30-year spin-up with perpetual 1965 climatology (the reference model) and imposed temperature anomalies of 1, 1.5, 2, and 2.5 ${}^{\circ }$C for the spin-up. The colour legend for (c) and (d) is indicated in (b).

[Figure omitted. See PDF]

Selected results of the sensitivity tests are plotted in Fig. 7, indicating the wide range of model behaviour that is possible with perturbations to the model inputs, parameter settings, and spin-up assumptions. An air temperature anomaly of $\pm \mathrm{1}$ ${}^{\circ }$C applied to the reference ERA climatology gives very different firn evolutions from 1965–2019, with warmer temperatures driving a shift to temperate firn conditions in the late 1980s (Fig. 7a and c). Warming of 2 ${}^{\circ }$C gives temperate firn for the entire period. In the other direction, a temperature anomaly of $-\mathrm{1}$ ${}^{\circ }$C is sufficient to maintain polythermal conditions at the site, precluding the development of deep temperate firn or a PFA. Similar results are obtained with perturbations of $\pm \mathrm{10}$ W m${}^{-\mathrm{2}}$ to the incoming longwave radiation (Supplement). Increases in meltwater infiltration that are stimulated by lower values of the irreducible water content (${\mathit{\theta }}_{\mathrm{wi}}<\mathrm{0.025}$) have a similar effect to warming, promoting meltwater infiltration, firn warming, and the earlier development of temperate firn.

The simulations are also sensitive to the initial conditions (Fig. 7b and d). Given evidence from Grew and Mellor (1966) that firn at 15 m depth was temperate in the mid-1960s near our core site, we introduce temperature anomalies from $+\mathrm{0.5}$ to $+\mathrm{2}$ ${}^{\circ }$C to the spin-up climatology. A perturbation of $+\mathrm{1.5}$ ${}^{\circ }$C creates temperate conditions to 12 m depth, and $+\mathrm{2}$ ${}^{\circ }$C is sufficient to create deep temperate firn that persists for several years (Fig. 7d). Firn refreezes in the 1970s within the model and eventually follows a similar path to the reference simulation but with a memory of warmer initial firn temperatures. This permits a more rapid transition (or return) to deep temperate conditions spurred by the heavy melt season in 2013. Overall, the model sensitivities in Fig. 7 indicate that a wide range of model solutions are possible at this site, indicating that Kaskawulsh Glacier firn is very close to the threshold for either temperate or polythermal conditions. We discuss this further below.

The temperature forcing required to produce temperate firn in the 1960s is relatively strong. The year 1965 that is used for the model initialization is representative of the long-term mean climatology of the site, with mean summer and annual air temperatures of $-\mathrm{2.5}$ and $-\mathrm{10.8}$ ${}^{\circ }$C. This compares with averages of $-\mathrm{2.4}\pm \mathrm{0.8}$ ${}^{\circ }$C and $-\mathrm{10.7}\pm \mathrm{0.9}$ ${}^{\circ }$C for the period 1965–2019 (Table 2). Incoming solar and longwave radiation in summer 1965 averaged 304 and 240 W m${}^{-\mathrm{2}}$, compared with long-term averages of 298 and 255 W m${}^{-\mathrm{2}}$. Net energy and melt were slightly lower than the long-term average due to low incoming longwave radiation, but overall, 1965 was a typical year. A warm anomaly of $+\mathrm{2}$ ${}^{\circ }$C represents 2.5 standard deviations above normal, giving a mean summer temperature of $-\mathrm{0.5}$ ${}^{\circ }$C; this would represent the warmest summer on record.

The initial firn density and ice content are relatively high when we force the model to produce temperate firn conditions in the mid-1960s through the $+\mathrm{2}$ ${}^{\circ }$C air temperature anomaly in the model spin-up. Values in 1965 are 679 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and 2.8 m, compared with reference model values of 641 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and 0.7 m. Figure 8 plots the subsequent firn temperature and density evolution if the $+\mathrm{2}$ ${}^{\circ }$C temperature anomaly is maintained from 1965 to 2019 and in the case where the model forcing is restored to the reference ERA climatology from 1965 to 2019. Subsurface temperature and density evolutions in the latter case parallel those of the reference model after a transient adjustment period of about a decade, while the perpetual $+\mathrm{2}$ ${}^{\circ }$C anomaly maintains dense and temperate firn. The decadal adjustment of firn density (Fig. 8b) is the “overturning” time of the firn core for downward advection of new snow and firn to 35 m depth. The temperature adjustment (Fig. 8a, c) does not follow this as it is governed by thermal diffusion timescales in the deep firn, giving a longer memory of the initial conditions.

Figure 8

Modelled (a) 10 m firn temperature and (b) average firn density for the reference model, for a 2 ${}^{\circ }$C temperature anomaly for the spin-up, and for a sustained temperature anomaly of $+\mathrm{2}$ ${}^{\circ }$C. (c, d) Firn temperature evolution for (c) the warm spin-up, followed by the reference climatology, and (d) sustained 2 ${}^{\circ }$C temperature anomalies.

[Figure omitted. See PDF]

5 Discussion

5.1 Firn characteristics and changes over time

The accumulation area of Kaskawulsh Glacier currently has indications of widespread meltwater percolation and refreezing. Meltwater is stored within the firn as ice, as indicated by the presence of ice layers and infiltration ice, and there is liquid water at a depth of $\sim \mathrm{35}$ m below the surface. The density of the firn has increased by about 15 % since 1964 in the first 7 m of firn due to the increased presence of ice layers. However, the firn in 1964 was not without meltwater percolation and refreezing; Grew and Mellor (1966) note the presence of refrozen ice lenses and glands and report evidence for meltwater infiltration and refreezing at depths of $\sim \mathrm{5}$ m. Nevertheless, the quantity and thickness of ice layers and lenses have increased towards present day, as reflected in the changes in the stratigraphy and the density (Fig. 4). The firn modelling also indicates decadal-scale increases in firn ice content and density (Table 2, Fig. 5e). For the reference model parameter settings and ERA climate forcing, the model predicts a significant increase in melting (Fig. 5b), driving increases in the depth of the melting and wetting fronts, meltwater percolation and run-off, and latent heat release associated with refreezing since the 1960s. This fundamentally changes the way the firn contributes to the mass balance of the glacier and englacial hydrological dynamics, as discussed further in Sect. 5.3. There are significant decadal firn warming trends in the model (Figs. 5 and 6), driven by the increases in melting and meltwater percolation. The modelling is not observationally constrained, however (Fig. 7 and Supplement), so the simulated firn warming is uncertain.

Increased firn meltwater and ice content as well as potential firn warming in recent decades will affect firn densification processes. Melting rounds snow grains and increases the rate of the first stage of densification. With enough melt to drive meltwater percolation through the snow and firn layer, meltwater can fill in air pockets and refreeze, further accelerating the transition from snow to ice (Cuffey and Paterson, 2010). The overall pattern of density measurements from 2018 resembles a logarithmic densification curve (Fig. 2), as is typical for Sorge's law of densification in dry snow (Sorge, 1935; Bader, 1954). However, with increasing meltwater percolation and refreezing effects, higher densities are common in the upper portions of the firn, as observed in our cores. Bezeau et al. (2013) report similar findings from Devon ice cap, where they found a depth–density reversal and suggest that Sorge's law no longer holds in areas of significant warming. To account for this, firn densification models are being revised to address the effects of ice layers and warming temperatures on the rate of densification (Reeh, 2008; Ligtenberg et al., 2011), and other studies are revising mass balance estimates based on dynamic densification rates (e.g., Schaffer et al., 2020).

5.2 Perennial firn aquifer

We found unequivocal evidence for a deep perennial firn aquifer on the upper Kaskawulsh Glacier, with excess water in the firn pore space below about 32 m depth. Some of this water drained during firn core acquisition (Video Supplement 1 and 2: 10.5446/50918 and 10.5446/50919, respectively). We cannot tell whether this PFA is a new feature at this site. Borehole temperature measurements from 1964 at a site close to our cores indicate temperate conditions at 15 m depth at this time (Grew and Mellor, 1966), and it is possible that firn has been temperate since that time, conducive to a PFA below the depth of the annual winter cold wave. There are no historical temperature measurements from greater firn depths at the site, and earlier coring efforts and radar surveys from the upper Kaskawulsh, Divide, or Eclipse sites make no comment or inference about the presence of liquid water, so we cannot attest to the age or origins of the PFA. It may well be a new feature.

The modelling results suggest that there are significant decadal increases in melting and refreezing since the 1960s at this site, driving firn warming, increased ice content, and densification (Table 2). The firn model predicts the development of wet, temperate conditions in the deep firn following the 2013 melt season, although it takes 4 years to fully develop (Fig. 6). This was triggered by meltwater penetration to 11 m depth in 2013, which is below the depth of penetration of the winter cold wave. Temperate conditions propagated downwards in the following years and persisted to 2019, supported by several more years with above-average melting. Deep meltwater percolation during these years would support the development and recharge of a PFA or perched water table at the glacier ice–firn interface. This agrees with the stratigraphy found in the field. The presence of firn that has not been visibly affected by meltwater overlying the PFA implies that deep meltwater infiltration through vertical piping may be an important process here and may allow the PFA to be recharged in a heavy melt season. In the model, deep recharge does not occur every summer after the establishment of a temperate firn column; the summer melt still needs to break through the winter cold layer, which typically extends to 6–7 m depth (Fig. 6). Also of interest in Fig. 6 is a large melt event in 2007, which led to meltwater infiltration and warming to about 9 m depth. This was similar to the 2013 melt event, but the summers of 2008 to 2010 were relatively cool (average JJA temperatures of $-\mathrm{2.8}$ ${}^{\circ }$C and melting of 111 mm w.e.), leading to refreezing in the upper 9 m of firn. Thawing of the full 35 m firn column to shift it from polythermal to temperate conditions requires several years of sustained melt forcing in the model.

There are significant uncertainties in the modelling associated with the climatological forcing, surface energy balance and firn model parameterizations, and initial conditions. The Supplement explores these in detail, while Fig. 7 provides an illustration of the range of simulated behaviour for different model settings. The “reference model” results presented in Figs. 5 and 6 should be seen as just one scenario, corresponding to our best estimate of the parameter settings. We lack local calibration and validation studies, so we cannot preclude different firn temperature and melt evolutions at this site, particularly given the inference of Grew and Mellor (1966) that firn at 15 m depth was temperate in the mid-1960s. The default model parameters and spin-up settings do not produce this; augmented warming or incoming radiation fluxes need to be introduced to the ERA climatology to produce temperate firn at this time. It is possible that strong melt seasons in the early 1960s created temporary temperate conditions in the upper firn column. Alternatively, the surface energy balance and firn hydrological models may underestimate the amount of melting and meltwater infiltration. The one firm conclusion is that the climatological and glaciological conditions on the upper Kaskawulsh Glacier are very close to the tipping point between polythermal and temperate conditions. A slight nudge to either side of the reference model settings can give either persistently sub-zero or persistently thawed conditions in the deep firn at this site (Figs. 7 and S1).

The presence of the deep PFA in 2018 indicates that it is currently temperate, despite mean annual air temperatures of about $-\mathrm{11}$ ${}^{\circ }$C. Meltwater refreezing releases enough latent heat to bring the firn to 0 ${}^{\circ }$C. All model simulations concur on this, although the long-term evolution is uncertain. We do not know the fate of the water that drains through the firn, but the reference model predicts a total drainage of 1.13 m w.e. over the 55-year simulation, most of this over the last decade. Some of this is retained within the PFA, but some can be expected to run off. The water in the PFA on Kaskawulsh Glacier is likely to be flowing, redistributing mass. The drill site was located high in the glacier's accumulation zone, with a gently sloping surface ($<\mathrm{0.6}$${}^{\circ }$) resulting in a subtle hydraulic gradient. We likely drilled into the top of the water table of the PFA, and with densities near the pore close-off density, it is likely the PFA does not extend much deeper. There may be downslope flow along the firn–ice interface as well as possible Darcian flow within the PFA itself (e.g., Christianson et al., 2015).

The liquid-phase meltwater retention on Kaskawulsh is similar to the PFAs found in the high-accumulation areas of southern Greenland and Svalbard (e.g., Miège et al., 2016; Christianson et al., 2015) and different than the water-saturated layers commonly found on temperate glaciers. PFAs that have been studied on temperate mountain glaciers typically have a saturated layer close to the surface (for example, 5 m below the surface at Storglaciären); have active discharge and recharge processes (Fountain and Walder, 1998; Schneider, 1999; Glazyrin et al., 1977); and appear to experience seasonal drainage over the winter months (Fountain, 1989, 1996; Jansson et al., 2003), likely due to high hydraulic gradients. Active water flow in the firn has been observed in 19 and 25 m pits at Abramov Glacier (Glazyrin et al., 1977) as well as Austfonna ice cap in 1985 at 7 m depth, where they also found sub-horizontal melt channels at 7, 15, and 30 m (Zagorodnov et al., 2006). In 2012, “water-saturated” firn was found at 40 m depth in an ice core from Mt. Waddington, British Columbia (Neff et al., 2012). However, they reported no significant alteration of chemistry from the melt above this layer, and no additional analysis of this layer was discussed (Neff et al., 2012). In 2015, a PFA was found on Holtedahlfonna ice field in northwestern Svalbard (Christianson et al., 2015), and in 2019 a PFA was investigated at Lomonosovfonna ice cap approximately 100 km to the southeast (Hawrylak and Nilsson, 2019).

According to Kuipers Munneke et al. (2014), PFA formation in Greenland is contingent upon a high annual snow accumulation, which helps to insulate the underlying firn from the winter cold wave. Mean annual temperatures in Greenland are well below 0 ${}^{\circ }$C, and PFAs require latent heat release from meltwater refreezing to warm the snow and firn to 0 ${}^{\circ }$C, along with meltwater penetration to depths of at least 10 m to evade the winter cold wave (Kuipers Munneke et al., 2014). The firn modelling suggests that meltwater penetration to depths of 10 m is rare at Kaskawulsh Glacier but can occur in strong melt seasons. Based on our measurements and earlier reports from the IRRP (Wood, 1963; Grew and Mellor, 1966), the estimated accumulation rate at our study site is 1.8 m w.e. yr${}^{-\mathrm{1}}$. This is similar to reported accumulation rates where PFAs have been identified in southeastern Greenland (e.g., Miège et al., 2016). Melt rates at southeastern Greenland PFA locations are also comparable to those on the upper Kaskawulsh. Miège et al. (2016) report 0.73 m w.e. yr${}^{-\mathrm{1}}$ over the time period 1979–2014, while Miller et al. (2020) estimated annual melt rates from 0.24–0.50 m w.e. yr${}^{-\mathrm{1}}$ in a PFA field study at 1700 m elevation in the Helheim Glacier catchment. Modelled melt rates on the upper Kaskawulsh Glacier are estimated at $\mathrm{0.52}\pm \mathrm{0.27}$ m w.e. yr${}^{-\mathrm{1}}$ from 1965–2019 (Table 2). Recent (2005–2017) Kaskawulsh melt rates increased to $\mathrm{0.72}\pm \mathrm{0.38}$ m w.e. yr${}^{-\mathrm{1}}$, similar to the long-term estimate of Miège et al. (2016) in southeastern Greenland and perhaps close to the threshold for PFA formation and recharge.

5.3 Implications for geodetic mass balance

Liquid water is commonly found in the temperate firn of low- and mid-latitude mountain glaciers and plays an important role in meltwater storage and glacier hydrology and mass balance (Fountain and Walder, 1998; Schneider, 1999). For example, storage of meltwater in a PFA accounts for as much as 64 % of internal accumulation for glaciers in Alaska and Sweden (Trabant and Mayo, 1985; Schneider and Jansson, 2004). In general, melt can result in net surface lowering in four main ways but with differing impacts on mass balance: (i) melt which runs off results in direct mass loss; (ii) melt which percolates and refreezes internally can result in surface lowering, with little to no mass loss; (iii) melt which makes it into a PFA likely contributes to mass loss, but the storage of liquid water at the firn–ice interface delays run-off from hours to weeks or longer (Jansson et al., 2003); and (iv) accelerated compaction of warming firn can result in accelerated surface lowering, without any mass loss. These components can be interrelated, and their relative importance depends on many factors including spatial and temporal variations in melt, PFA thickness, the presence of ice lenses, and firn temperature, so their effects are difficult to disentangle. We provide further information about these components below, including an assessment of their relative importance for Kaskawulsh Glacier.

The climate reanalysis suggests that the effects of meltwater storage through refreezing or liquid retention are increasing as the climate is warming. Geodetic mass balance measurements are compromised by climate-change-induced densification changes that are not accounted for when interpreting surface lowering of the accumulation zone (Reeh, 2008; Huss, 2013). Mass balance studies in Greenland indicate that changing melt regimes, meltwater refreezing, and the unknown density and pore capacity of snow and firn pose significant uncertainties when modelling the surface mass balance of ice sheets (Lenaerts et al., 2019). Meltwater retention as porewater or refrozen ice will delay surface run-off, dependent on the water storage characteristics of firn (e.g., pore space availability, water at interstitial grain boundaries) (Fountain and Walder, 1998; Schneider, 1999). If ice layers become too extensive or thick, they can form an “ice slab”, a thick impermeable barrier that leads to enhanced surface run-off (MacFerrin et al., 2019). The thickness of ice layers that prevents percolation is not well understood. For example, in Greenland 12 cm thick ice layers were still permeable (Samimi et al., 2020), whereas Bell et al. (2008) reported that a 1–2 cm ice layer prevented percolation at Devon ice cap, Canada. These phenomena and effects are not limited to Greenland and the high Arctic. This study demonstrates that Kaskawulsh Glacier also experiences meltwater storage in the form of ice layers and liquid water retention (as a PFA), with potentially significant recent changes in firn structure and meltwater retention capacity. The increases in firn density and ice content found on Kaskawulsh Glacier appear to be similar to other high-accumulation Arctic regions (Pohjola et al., 2002; De La Peña et al., 2015; Bezeau et al., 2013).

The surface energy balance model is not observationally constrained at this site, so we do not have quantitative confidence in the modelled mass balance and melt rates, but the reconstructed trends indicate a $\sim \mathrm{70}$ % increase in summer meltwater production at this site since the 1960s, leading to increased rates of refreezing and also the onset of meltwater run-off in recent years. We neglected the potential influence of summer rainfall in this study as we believe that rain events remain rare at this site. However, they likely happen from time to time and could become more prevalent in a warming climate. Rain events would add sensible heat and liquid water to the snow and firn, further increasing snow and firn temperatures; water infiltration into the firn; and meltwater run-off, where there is inadequate cold content or pore space to retain this water (as in recent years at the site). In situ field studies are needed to confirm and constrain the meteorological and energy balance conditions on the upper Kaskawulsh Glacier to inform the mass balance processes for both the glacier and the PFA.

Melt totals are much less than the annual accumulation ($\sim \mathrm{1.8}$ m w.e.), so the site remains within the accumulation area of the glacier, with most of the meltwater refreezing. Increases in meltwater refreezing have driven decadal-scale firn warming along with increases in ice content and firn density at this site. The modelling suggests a $\sim \mathrm{5}$ % increase in firn density since the 1960s and a doubling of the ice content, from 1.1 to 2.3 m over the full 35 m snow and firn column. Increases in summer melting over the last decade are associated with meltwater infiltration and firn warming in the deep firn, with the likelihood of meltwater run-off from the site in recent years. Within the model, 91 % of surface melt refreezes over the 55-year simulation, but this declines to 73 % for the period 2005–2017. The remaining 27 % drains to the deep firn through this period, where it is either retained within the PFA, or it may drain from the system. In the firn model, a total of $\sim \mathrm{0.7}$ m of meltwater is stored as liquid water in the deep firn, and $\sim \mathrm{0.7}$ m w.e. “runs off” through the period 2005–2017, draining through the bottom layer and leaving the system. In reality, this meltwater may drain through lateral transport in the PFA or at the ice–firn interface.

The modelled 2018 firn core has an ice content of 2.6 m compared to a total ice content of $\mathrm{2.33}\pm \mathrm{0.26}$ m measured in Core 1. The modelled ice content is a completely independent estimate but is in reasonable agreement with the firn core, giving some confidence in the modelled refreezing. The model may slightly overestimate the melt or the meltwater infiltration, given that the modelled ice content is about 15 % too high. However, that inference is not consistent with the apparent cold bias in the model spin-up. Alternatively, firn in the model may be too cold through much of the simulation, causing an overestimate of the modelled meltwater refreezing and retention capacity. If this is the case, run-off (summer mass losses) from the site will be higher than our estimates, with negative implications for Kaskawulsh Glacier mass balance.

The accumulation zone of Kaskawulsh Glacier is estimated to have experienced a minimum of $\mathrm{0.73}\pm \mathrm{0.23}$ m of surface lowering due to internal refreezing over the period represented by Core 1, which we estimate to be 12.5 years, or $\mathrm{0.06}\pm \mathrm{0.02}$ m yr${}^{-\mathrm{1}}$ from 2005–2017. This estimate of thinning is likely low because neither the meltwater retention due to the infiltration ice nor the presence of the PFA is included in this estimate. In previous measurements of surface elevation changes on Kaskawulsh Glacier, Foy et al. (2011) found that the accumulation zone thinned by an average of 0.04–0.11 m yr${}^{-\mathrm{1}}$ from 1977–2007, with a total thinning of 1–3 m over this period. Larsen et al. (2015) reported mean elevation losses of 0–1 m yr${}^{-\mathrm{1}}$ towards the head of the glacier from 1995–2000. The thinning signal due to meltwater percolation and refreezing is within the estimates of Foy et al. (2011) and Larsen et al. (2015), suggesting that some or all of the reported lowering could be due to mass redistribution and not mass loss. The density of the firn has increased from 1964–2018 due to meltwater percolation and refreezing. It is therefore likely that the surface has lowered since 1964 because of this increased densification.

6 Conclusion

The upper accumulation zone of Kaskawulsh Glacier firn has undergone significant changes since 1964, becoming denser and more ice-rich. The mean density of the first 32 m of firn (4.2 to 36.2 m below the surface) was $\mathrm{698}\pm \mathrm{5}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. Analysis of historical density data indicates that the firn of Kaskawulsh Glacier has become up to 15 % denser since the early 1960s due to the increased ice content and melt-affected firn. Increases in firn density due to meltwater refreezing over the 13-year period represented by the firn core (2005–2018) are equivalent to a surface lowering of $\mathrm{0.73}\pm \mathrm{0.23}$ m ($\mathrm{0.06}\pm \mathrm{0.02}$ m yr${}^{-\mathrm{1}}$), and this rate of surface lowering is likely increasing in association with the overall densification, connected to increases in firn temperature as well as ice content. Though not observationally constrained, and therefore uncertain, the modelling results suggest the likelihood of significant increases in melting and refreezing since the 1960s at this site, driving decadal-scale increases in firn temperature, ice content, and density. The estimates of firn density and the evidence for densification can help to inform geodetic mass balance measurements from this region.

Our study also illustrates a high-elevation accumulation area in the St. Elias Mountains that is undergoing a transformation in response to climate change. The firn on upper Kaskawulsh Glacier now contains a PFA, which has likely developed over the past decade. Firn modelling suggests that the PFA has developed in response to increased summer melting, meltwater infiltration, and firn warming from the associated latent heat release. The Kaskawulsh Glacier PFA needs to be more widely studied as the spatial extent and depth of the aquifer are not yet known. Ground-penetrating radar measurements may provide a method to investigate the spatial extent of the feature. Use of an electrothermal drill that can drill through water-saturated firn may allow estimations of the depth of the firn aquifer as well as subsequent studies on the potential flow of the water within the aquifer. The firn modelling suggests that this site is experiencing meltwater run-off over the past decade in relation to the development of the PFA. A better understanding of this feature is needed to quantify the extent of meltwater retention and the mass balance of Kaskawulsh Glacier. This region will likely continue to experience increasing amounts of surface melt and refreezing within the snowpack and firn, extending to higher elevations, so there is some urgency to obtain climate records from this region.

Code availability

The MATLAB code used for the firn modelling is publicly available at 10.5683/SP2/WRWJAZ (Marshall, 2021). This is research code under active development; the code is not supported and is not thoroughly documented. Users are free to adapt and apply it with appropriate acknowledgement of the source, and questions about the code can be addressed to Shawn J. Marshall.

Data availability

Raw density data are available by contacting the corresponding author.

Video supplement

The video supplements related to this article are available at 10.5446/50918 (Ochwat, 2021a) and 10.5446/50919 (Ochwat, 2021b). They are videos filmed during the fieldwork expedition and capture the moment when the team drilled into the perennial firn aquifer.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-15-2021-2021-supplement.

Author contributions

NEO and ASC collected field data. ASC ran ion analyses, supervised the field campaign, and helped with figures. SJM contributed to the design and funding of the study and was responsible for the firn modelling. BJM and SJM provided supervision during the project. LC provided weather station data and contributed to the collection and interpretation of data. NEO analysed the data and wrote the manuscript, to which all co-authors contributed.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Parks Canada for permission to conduct this research in Kluane National Park, under research and collection permit KLU-2018-28117. We are grateful for the field crew members Étienne Gros and Peter Moraal, Icefield Instruments Inc., for assistance with coring and members of the University of Ottawa Glaciology and Northern Field Research classes, particularly Jean Bjornson, for assistance with the snow pit measurements. The Arctic Institute of North America, Kluane Lake Research Station, and Icefield Discovery supported fieldwork logistics. We thank Kristina Miller of the University of Calgary for field support and countless glaciological discussions and Eduard Khachatrian for translating the Glazyrin et al. (1977) paper.

Financial support

This research has been supported by Polar Knowledge Canada with the grant titled “Cryosphere-Climate Monitoring at Kluane Lake Research Station, Yukon Territory, Canada” and the Natural Sciences and Engineering Research Council (NSERC) of Canada with the Discovery Grant titled “Modelling glacier and ice sheet response to climate change”.

Review statement

This paper was edited by Etienne Berthier and reviewed by three anonymous referees.