1. Introduction
Svalbard is one of the fastest warming regions on Earth (Nordli and others, Reference Nordli, Przybylak, Ogilvie and Isaksen2014; Isaksen and others, Reference Isaksen2016). Glaciers, which cover 57% of the archipelago (Nuth and others, Reference Nuth2013), are significantly impacted by increasing surface melt due to rising temperatures. During 2000–19, Svalbard glaciers lost mass at a rate of 
$-8 \pm 6$ Gt a−1 (Schuler and others, Reference Schuler2020), and mass loss is projected to increase in the coming decades (Hock and others, Reference Hock2019; van Pelt and others, Reference van Pelt WJJ, Schuler, Pohjola and Pettersson2021; Geyman and others, Reference Geyman, van Pelt, Maloof, Aas and Kohler2022).
The impact of increasing melt rates on total mass loss and associated sea level rise depends on the amount of water that runs off from the glacial system. Here, firn plays an important role by retaining meltwater and reducing runoff in the accumulation zones (Pfeffer and others, Reference Pfeffer, Meier and Illangasekare1991; Harper and others, Reference Harper, Humphrey, Pfeffer, Brown and Fettweis2012). In Svalbard glaciers, refreezing is a significant component of the climatic mass balance – the sum of the surface and internal mass balance (Cogley and others, Reference Cogley2011) – accounting for around 25% of the annual accumulation (Aas and others, Reference Aas2016; Østby and others, Reference Østby, Schuler, Hagen, Hock, Kohler and Reijmer2017; van Pelt and others, Reference van Pelt2019; Schuler and others, Reference Schuler2020). However, the reduction in firn cover extent and concomitant firn densification have led to a decrease in firn pore space in Svalbard glaciers, thereby diminishing their firn buffering capacity (Østby and others, Reference Østby, Schuler, Hagen, Hock, Kohler and Reijmer2017; van Pelt and others, Reference van Pelt2019; Noël and others, Reference Noël2020). This negative trend is projected to continue during the 21st century (van Pelt and others, Reference van Pelt WJJ, Schuler, Pohjola and Pettersson2021).
Numerous studies have reported a drastic reduction in firn retention capacity in some polar regions due to the formation of ice slabs (> 1 m thick ice layers) near the surface, which prevent meltwater flow to greater depths (Machguth and others, Reference Machguth2016; MacFerrin and others, Reference MacFerrin2019). Ice slabs have been discovered in the western and northern Greenland (MacFerrin and others, Reference MacFerrin2019), and in the Canadian Arctic (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013). Furthermore, firn aquifers, liquid water reservoirs in firn pore space that do not refreeze over winter (Forster and others, Reference Forster2014), have been detected along the southeastern margin and northwest of Greenland (Forster and others, Reference Forster2014; Miège and others, Reference Miège2016), in central (van den Akker and others, Reference van den Akker, van Pelt, Petterson and Pohjola2025) and northwestern Svalbard (Christianson and others, Reference Christianson, Kohler, Alley, Nuth and van Pelt2015), and in Antarctica (Montgomery and others, Reference Montgomery2020). In the past two decades, an expansion of both ice slabs and firn aquifers has been observed on the Greenland Ice Sheet (Horlings and others, Reference Horlings, Christianson and Miège2022; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). The presence of ice slabs or firn aquifers predominantly depends on the rates of surface water input (melt and rain) and snow accumulation (Brils and others, Reference Brils2024).
How the stratigraphic and thermal properties of firn will change in the warming climate is a question of high importance due to the impact on meltwater retention capacity, glacier hydrology, and dynamics. By sealing off pore space, ice slabs can increase runoff and accelerate mass loss (Machguth and others, Reference Machguth2016; MacFerrin and others, Reference MacFerrin2019), while perennial firn aquifers can delay runoff (Forster and others, Reference Forster2014). Moreover, ice velocity may be largely affected by water potentially draining from aquifers to the glacier bed through crevasses (Poinar and others, Reference Poinar2017; Cicero and others, Reference Cicero2023), and by liquid water storage in firn altering the thermal regime of the ice (Irvine-Fynn and others, Reference Irvine-Fynn, Hodson, Moorman, Vatne and Hubbard2011; Wilson and Flowers, Reference Wilson and Flowers2013).
Previous observational studies on long-term changes in firn properties have largely focused on Greenland and the Canadian Arctic. This work extends those studies to Svalbard by focusing on Austfonna – an Arctic ice cap where changes in firn due to recent atmospheric warming have not yet been documented. We combine field observations and modeling to gain a better understanding of the stratigraphy and thermal state of the firn, as well as the meltwater retention processes on the ice cap. We present data from 16 firn cores drilled on Austfonna between 1958 and 2022, exploring spatial and temporal changes in density and ice content. Additionally, we provide a firn temperature record spanning from 2009 to 2022. Furthermore, we apply the CryoGrid community energy balance model to simulate changes in firn temperature, water content, and density. In situ surface mass balance, temperature, and density data are used to validate the model. The results are relevant to other glaciers and ice caps in the Arctic, as the rapid warming in Svalbard may signal broader changes that could affect other regions in a similar way.
2. Study area
Austfonna is located on Nordaustlandet in northeastern Svalbard, centered at 79.7 ∘N and 24 ∘E (Fig. 1). With an area of 
$ \gt 7800\,\mathrm{km}^2$ (Moholdt and Kääb, Reference Moholdt and Kääb2012), Austfonna is the largest ice cap in the Svalbard archipelago, accounting for a third of its ice volume (Fürst and others, Reference Fürst2018). The ice cap ranges in elevation from sea level up to around 800 m above sea level (a.s.l.), where the thickest ice of about 500 m is found (Dowdeswell and others, Reference Dowdeswell1986). Three of the marine-terminating basins have been observed to surge in the past (Lefauconnier and Hagen, Reference Lefauconnier and Hagen1991), and two ongoing surges have been occurring in Basin-3 and Etonbreen since 2012 and 2023, respectively (Dunse and others, Reference Dunse2015; Mannerfelt and others, Reference Mannerfelt, Schellenberger and Kääb2025).
(a) Study area Austfonna ice cap, (b) sites of newly and previously drilled firn cores presented in this study, and (c) location on the Svalbard archipelago. 50 m contour lines and glacier basins are derived from a digital elevation model (DEM) based on aerial photos from 2011 (Norwegian Polar Institute, 2014b). Elevation contours are in m a.s.l.. The basemap is derived from an orthophoto from 2012 (Norwegian Polar Institute, 2014a).
The mean annual air temperature measured at a weather station on Austfonna at 370 m a.s.l. over the period 2004–12 was 
$-8.3^\circ\mathrm{C}$, with hourly values ranging from 
$-43.2^\circ\mathrm{C}$ to 
$9.1^\circ\mathrm{C}$ (Schuler and others, Reference Schuler, Dunse, Østby and Hagen2014). Temporal variability is high; even in winter, the temperature can occasionally rise to above 0
$^\circ\mathrm{C}$, while summer snowfalls may also occur. The melt season typically lasts from mid-June to the end of August or beginning of September (Schuler and others, Reference Schuler2007, Reference Schuler, Dunse, Østby and Hagen2014). The primary moisture source for Austfonna’s precipitation is the Barents Sea east-southeast of Svalbard (Pinglot and others, Reference Pinglot, Hagen, Melvold, Eiken and Vincent2001; Taurisano and others, Reference Taurisano2007).
Surface mass balance measurements using the glaciological method are available for Etonbreen, a 640 km2 glacier draining the ice cap to the west. From the mass balance years 2003/04 to 2021/22, the winter balance was 0.40 m w.e. and the summer balance was −0.54 m w.e. resulting in an annual balance of −0.14 m w.e. (Fig. 2). While the winter balance is relatively stable, the summer balance, and consequently also the annual balance, vary considerably from year to year. Particularly high melt occurred in the summers of 2004, 2013, 2020 and 2022 with annual mass balances of −0.56 m w.e., −0.59 m w.e., −0.76 m w.e. and −0.82 m w.e., respectively. The equilibrium line altitude (ELA) of Etonbreen was in average 556 m a.s.l. in 2004–22 (Fig. 2). Accumulation rates are generally higher in the southeast of the ice cap compared to the northwest, as shown by observations from shallow ice cores (Pinglot and others, Reference Pinglot, Hagen, Melvold, Eiken and Vincent2001) and ground-penetrating radar (GPR) (Taurisano and others, Reference Taurisano2007; Dunse and others, Reference Dunse2009). This is consistent with the primary moisture source originating from the Barents Sea. Mapped glacier facies based on GPR measurements indicate that even the summit area has been subject to melting during the summer, with no dry snow zone present (Dunse and others, Reference Dunse2009). Superimposed ice is typically found in the lower accumulation area.
(a) Equilibrium-line altitude (ELA) and (b) annual, winter, and summer glacier-wide surface mass balance of Etonbreen derived from stake and snow pit measurements for the mass-balance years 2003/04 to 2021/22. Note that the ELA is derived from the mass-balance gradient and therefore exceeds the maximum elevation of the ice cap in 2019/20.
3. Data and methods
3.1. Firn cores
We use data from five new firn cores retrieved in 2022 and 11 earlier cores from various sources (Table 1) and locations (Fig. 1). To identify where and when each core was drilled, the following nomenclature is used: site_year, where “site” refers to the location and “year” to the core retrieval year.
Firn cores used in this study. The last number in the core name refers to the core retrieval year. The elevations are taken from the original source, or, if not available, derived from the DEM by Norwegian Polar Institute (2014b). The drilling depth refers to the depth below the snow surface. The segment length is the length over which density was measured. The sources refer to (1) this study, (2) Langhammer Reference Langhammer(2014), (3) Watanabe and others Reference Watanabe(2001), (4) Schytt Reference Schytt(1964), and (5) Pinglot and others Reference Pinglot, Hagen, Melvold, Eiken and Vincent(2001).
3.1.1. Firn cores drilled in 2022
Between 11 and 16 May 2022, a total of five shallow firn cores were collected at Etonbreen. Four cores (Eton-4_22, Eton-5_22, Eton-6_22 and Camp_22) were retrieved along the mass balance stake transect, with one core (RR-3_22) being retrieved slightly outside of it. The cores were extracted using a KOVACS Mark III ice coring system, which consists of a 1 m core barrel with a diameter of 7.25 cm, extension rods, and an adapter for an electric drill machine. Cores RR-3_22 and Eton-5_22 were collected from the snow surface, while the three other cores were collected from the snow-firn interface of a snowpit.
The drilling procedure and measurement techniques closely followed the methods of Rennermalm and others Reference Rennermalm(2022). After each core section was extracted, its length was measured, and it was photographed in a half-pipe against an aligned ruler. Core section lengths varied substantially (
$0.32 \pm 0.22$ m, mean ± standard deviation). Sections often broke naturally into smaller pieces or were cut manually so that their lengths varied roughly between 
$0.1-0.2$ m. These shorter pieces (hereinafter called segments following Rennermalm and others Reference Rennermalm(2022)) were weighed using an OHAUS CS Series scale (capacity 2000 g, accuracy 1 g), and their length and diameter were measured with a stiff measuring tape and a caliper, respectively. Additionally, the segment intactness (in percent) was visually estimated, as segments were often not complete (e.g., some material was lost at the edges). The thickness and depth of ice layers were logged with a resolution of 0.005 m and thinner ice layers were recorded with a thickness of 0.0025 m. Coring became more challenging when approaching the glacier ice, as the barrel tended to get stuck and thus, the total lengths of the cores vary considerably (
$2.8-9.7$ m). We refer to glacier ice as ice that we assume is not underlain by firn, which distinguishes it from an ice layer.
We calculate the density of each core segment from the measurements of diameter, length, weight, and intactness. Density for 15% of core segments exceeded the density of pure ice (917 kg m−3) and was corrected as described in the Supplementary Text. The mean density (± standard deviation) of all pure ice segments 
$853 \pm 52\,\mathrm{kg\,m}^{-3}$ is consistent with previous studies that found ice densities ranging from 
$843\pm36\,$ kg m−3 to 873 ± 25 kg m−3 (Harper and others, Reference Harper, Humphrey, Pfeffer, Brown and Fettweis2012; Machguth and others, Reference Machguth2016; Rennermalm and others, Reference Rennermalm2022).
3.1.2. Previously drilled cores
We used data from 11 cores collected near the summit area between 1958 and 2015 (Table 1). Four cores were drilled between 2004 and 2015 by the University of Oslo and the Norwegian Polar Institute using similar methods as described above. Data from seven cores have previously been published (Table 1). One of them was drilled at Camp XV in 1958 (Schytt, Reference Schytt1964), two at site J in 1998 and 1999 (Watanabe and others, Reference Watanabe2001), two at site X in 1999 (Pinglot and others, Reference Pinglot, Hagen, Melvold, Eiken and Vincent2001) and 2010 (Langhammer, Reference Langhammer2014), and two at Summit in 2012 and 2013 (Langhammer, Reference Langhammer2014). Among the various cores presented by Pinglot and others Reference Pinglot, Hagen, Melvold, Eiken and Vincent(2001), we use only core X_99 due to its proximity to the summit area and its comparable elevation to the other drill sites.
The density data of all cores were corrected where necessary (Supplementary Text). Some earlier cores were drilled from the snow surface, and the depth of the snow-firn interface had to be estimated. This was done based on snow pit measurements from the same year or by scrutinizing density profiles (typically a clear spike in density was found below the snow-firn interface).
3.2. Firn temperature
Firn temperatures have been monitored near the summit at site Camp (Fig. 1) by the University of Oslo and the Norwegian Polar Institute since 2009. A thermistor string with 10 temperature sensors was inserted into a borehole on 6 May 2009, which was then backfilled with snow. The sensors were spaced at intervals of 1 to 3 m, with the lowest sensor located at 13 m below the snow surface. Data were recorded until 29 April 2014. A second thermistor string was installed, and logging started on 3 May 2014, with the lowest sensor positioned at a depth of 12 m. An 8-month gap from August 2018 to April 2019 occurred due to full memory of the storage unit. A third thermistor string was installed on 3 May 2015 logging data until 2022. The lowest sensor was located at a depth of 11.2 m at installation.
Since continuous surface height measurements at site Camp are not available, we used surface heights simulated by the CryoGrid community model (Schmidt and others, Reference Schmidt, Schuler, Thomas and Westermann2023; Westermann and others, Reference Westermann2023) (described in the next section) to track the temporal evolution of sensor depths. For all three thermistor strings, the sensor depths were initialized from the depths recorded at installation. We did not account for movement of the sensors due to firn compaction when assessing sensor depths over time.
We corrected the sensor readings for systematic biases by shifting each thermistor string’s time series by the observed bias at 
$0^\circ\mathrm{C}$, by using temperature measurements during melting conditions. The average bias for the three thermistor strings was 
$0.25^\circ\mathrm{C}$, 
$-0.10^\circ\mathrm{C}$, and 
$-0.16^\circ\mathrm{C}$, respectively. Data from the first two months after the installation of each thermistor were discarded to ensure records were not affected by disturbances from the drilling procedure.
3.3. Modeling
3.3.1. Model description
To simulate subsurface firn temperature, density, and water content at site Camp, we used the multi-layer, energy balance model CryoGrid (Westermann and others, Reference Westermann2023). While originally developed for permafrost, it has recently been adjusted for glacier applications (Schmidt and others, Reference Schmidt, Schuler, Thomas and Westermann2023). Our model set-up closely follows Schmidt and others Reference Schmidt, Schuler, Thomas and Westermann(2023). Below, we present a brief description of the model while further details and model parameter values used in this study are given in the Supplementary Text and Supplementary Tables S1-S3.
The model calculates the energy flux toward the surface from the surface energy balance. When the surface temperature is below the melting point, the excess energy is used to remove the cold content. Once the surface becomes temperate, excess energy is used for melting. The subsurface scheme accounts for snow and firn compaction, heat transport, refreezing, water percolation, and runoff. The parameterizations of subsurface processes are based on the snowpack model Crocus (Vionnet and others, Reference Vionnet2012), with a few modifications.
CryoGrid uses the so-called tipping bucket approach (Westermann and others, Reference Westermann2023) to simulate the water balance; water infiltrates by gravity to the next, underlying layer, when the field capacity (0.05 m3 m−3) of the current layer is exceeded.
Runoff may occur in layers above an impermeable layer (defined by density > 830 kg m−3) if the amount of water surpasses the field capacity (Westermann and others, Reference Westermann2023). Runoff is delayed by a runoff time scale, allowing water to be stored in a layer until it refreezes or runs off. In contrast to previous applications, here some water is allowed to flow through individual “impermeable” layers if they overlay permeable firn. The amount of water depends on the hydraulic conductivity of the layer above it, which is influenced by the snow grain diameter, snow density, and saturation (see Schmidt and others Reference Schmidt, Schuler, Thomas and Westermann(2023)). This modification was necessary to prevent pooling of water near the surface, as this was not supported by observations. Water input to the surface layer originates from snow/firn melt or rainfall. Lateral water movement is not considered.
Calculations are performed on a vertical grid with user-prescribed initial vertical spacing. Here, the top layer was assigned a thickness of 0.01 m w.e., and layer thickness doubles every 20 layers. Snow and firn layers are merged or split using a layer thickness criterion (Zweigel and others, Reference Zweigel2021; Schmidt and others, Reference Schmidt, Schuler, Thomas and Westermann2023): when the layer exceeds the user-defined snow water equivalent threshold by more than 50%, it is divided in half. A layer is merged with the underlying one, when the current layer is 50% below the threshold. Hence, the surface snow or firn layer thickness ranges from 0.005 m w.e. to 0.015 m w.e.. Below the firn-ice interface, we used 38 ice layers with a total thickness of 30 m, which remain constant throughout the simulation period.
The model was forced by 2 m air temperature, wind speed, specific humidity, incoming shortwave and longwave radiation, precipitation, and air pressure from the Copernicus Arctic Regional ReAnalysis (CARRA) (Schyberg and others, Reference Schyberg2020). The temporal resolution is 3 hours, and the spatial resolution is 2.5 km. We use the grid point closest to site Camp’s elevation.
3.3.2. Model initialization and application
The model was run over the period from 1 January 2009 to 31 December 2022 to encompass the entire period of firn temperature observations. Since thermistor measurements were not available before May 2009, we used the data from 1 January 2010 to represent the thermal conditions at the start of the simulation (Supplementary Table S3). Firn density was initialized with firn core observations from Summit_2012, as this is the firn core closest in time to our simulation start date and has a long enough density record. Both temperature and density observations were linearly interpolated to the model layers. In the absence of direct observations of the depth of the firn-ice interface at site Camp in 2009, we assumed a depth of 10 m, based on available observations nearby (Core X_10, GPR measurements from Dunse and others Reference Dunse2013, and average seasonal snow depth (Dunse and others, Reference Dunse2009)).
4. Results
4.1. Firn density and stratigraphy
Figure 3 illustrates the stratigraphy and density profiles of the five new cores drilled in May 2022. Core Eton-4_22 (506 m a.s.l.) and core Eton-5_22 (637 m a.s.l.) are located slightly below and above the average ELA for the period 2004–22 (556 m a.s.l.), respectively. Therefore, we assume that glacier ice is found at depths of −1.3 m and −3.1 m in Eton-4_22 and Eton-5_22, respectively. However, for cores Eton-6_22, Camp_22, and RR-3_22, the depth of the firn-ice interface is less clear; we assume the bottom layers to be ice layers rather than part of the continuous ice body underlying the firn. All cores display numerous ice layers throughout their entire length. The thickest ice layer was found at Eton-6_22 and was at least 0.79 m thick. Firn thickness varies substantially with elevation, increasing from 0.12 m (Eton-4_22) to > 10 m at the highest core (RR-3_22) 18 km away. Density and ice content increase with depth (Fig. S1). While density and ice content are higher in the first few meters at Eton-5_22, there is no discernible decreasing trend with increasing elevation among the three highest cores (
$736-784$ m a.s.l.) (Fig. S1).
Figure 4 shows the firn stratigraphy and density of the 13 cores drilled in close proximity to the summit between 1958 and 2022 (Fig. 1). The drilling sites are up to 7 km apart from each other, but the elevation differences between them are < 41 m and thus, spatial variations in climatic conditions can be assumed small. The thickest ice layer among these cores is found in core J_99 (0.81 m located below −6 m depth). The depth of the firn-glacier ice transition varies among the cores and over time. At site J, this depth was −8.3 m in 1998, −8.7 m in 1999, and −6.8 m in 2004. Core J in 2004 displays the shallowest firn column among all cores, while the deepest firn column is found at Summit in 2013 (12.3 m). Despite spanning more than 60 years, the cores drilled in the summit area show no apparent trend in the abundance or distribution of ice layers.
Firn stratigraphy and density for cores drilled in the summit area (sites in Figure 1) between 1958 and 2022. Elevations are in m a.s.l from the DEM by Norwegian Polar Institute (2014b) or from the original sources (Table 1) ranging from 750 m a.s.l. to 791 m a.s.l.. Black lines show the mean density for each core segment. Only the upper 16 m is shown for cores that are longer. Dashed blue line indicates the start of glacier ice for cores with no stratigraphy data based on densities and/or other field observations.
4.2. Firn temperature
Measured subsurface temperatures at site Camp from May 2009 to December 2022 show seasonal variations between winter cooling and summer melt (Fig. 5a). However, what stands out is a clear transition in the firn thermal state. Between 2009 and 2013, firn is seasonally temperate during summer and fall, warmed by meltwater percolation and refreezing, while temperatures otherwise remain below 
$0^\circ\mathrm{C}$. No temperate firn persists throughout the entire year, suggesting that all meltwater refreezes in the subsequent winters. However, after 2015, firn becomes temperate year-round at depth, with sub-freezing temperatures found only near the surface during winter. This thermal regime shift occurred between 2013 and 2015, although precise timing is unclear due to data gaps and the temperature sensors not reaching sufficient depth. Field observations in May 2014 at sites Camp and RR-3 support the presence of liquid water (and temperate firn): liquid water was encountered in the lower part of the borehole at Camp during the installation of a new thermistor string, and at −10.2 m depth when drilling a firn core at RR-3. Based on the temperature record and field observations from the pre-melt season in 2014, we infer that the shift to year-round temperate firn at depth likely occurred in 2013.
(a) Observed and (b) modeled firn temperatures at site Camp, and (c) their differences, (d) modeled liquid water content (LWC) and (e) density from 2009 to 2022. The black solid line indicates the surface as derived from the simulation, and the dashed line shows the modeled firn-glacier ice interface (derived from density ≥ 830 kg m−3 with no permeable layers below). The depth is relative to the surface on 6 May 2009. The black dots in (a) show the temperature sensor locations at installation. The white line outlines the presence of liquid water which is inferred from temperatures 
$\ge -0.1^\circ\mathrm{C}$ (panel a) and liquid water content ≥ 0.001 kg kg−1 (panel b). The pink line in (d) outlines fully saturated firn (saturation ≥ 99%).
Meltwater percolation depth, inferred from temperatures 
$\geq -0.1^\circ\mathrm{C}$, increased between 2009 and 2013 (Fig. 5a). In summer 2009, meltwater percolated down to approximately 7 m below the surface, while by summer 2013, it reached depths exceeding 10 m below the surface. After the shift to temperate conditions at depth in 2013, meltwater percolated considerably deeper, reaching at least 13 m below the surface by the end of the record. Interestingly, between 2009 and 2013, meltwater did not percolate below a depth of approximately 8 m below the May 2009 surface. This indicates either an ice layer or glacier ice blocking the percolation.
4.3. Model validation
The model was validated by comparing the simulation results with observations of surface mass balance, firn temperature, and firn density. The observed and modeled cumulative surface mass balance over the mass-balance years 2009/10 to 2018/19 align well (Fig. 6).
Cumulative surface mass balance derived from stake and snow pit measurements at site Summit and simulated for site Camp for the balance years 2009/10 to 2018/19. The model assumes that the start of the mass-balance year is 1 September.
The simulated subsurface temperature closely follows the measurements. The correlation coefficient between daily modeled and observed temperatures is 0.70, and the root mean square error is 
$2.6^\circ\mathrm{C}$ (Fig. 5c). Except when the firn is temperate, the simulated temperature profile is colder with a mean difference between the simulated and observed temperature of 
$-1.2^\circ\mathrm{C}$. Especially from summer 2010 to spring 2013, the simulated firn is notably colder and meltwater percolation is shallower than observed (Fig. 5a,b). For example, in summer 2010, meltwater percolates to 7.5 m below the surface, whereas the simulated percolation reaches 2 m below the surface. After summer 2013, both the simulation and observations generally indicate temperate conditions at depth, with temperature differences concentrated near the surface. In some years, such as winters 2010 and 2016, the simulated firn is too warm, with differences of up to 
$2^\circ\mathrm{C}$ (Fig. 5c).
Simulated and observed densities in six firn cores drilled in the summit area of Austfonna generally match well, despite small-scale discrepancies due to stratigraphic features such as ice layers (Fig. 7). For three cores, we find negative model biases of −2 kg m−3 (RR3_15), −22 kg m−3 (Summit_12) and −44 kg m−3 (RR3_14). For three other cores, we find positive model biases of 4 kg m−3 (Summit_13), 29 kg m−3 (RR3_22), and 74 kg m−3 (Camp_22). In core Camp_22, where the model overestimates density the most, the largest difference occurs below −4 m depth. Overall, the model validation provides confidence in the model’s performance, despite some shortcomings in accurately predicting firn temperature and density.
Observed and simulated density in the upper 10 m of six cores drilled in the summit area between 2012 and 2022. Observed and simulated densities refer to averages over each measured segment and each simulated layer thickness, respectively.
4.4. Simulated firn conditions
Modeled firn temperature and liquid water content indicate that between 2009 and 2013, the firn is seasonally temperate, with summer melt percolating to depths of 
$1.7 - 5.0$ m below the surface (Fig. 5b,d). However, in 2013, meltwater is simulated to percolate down to 13.5 m below surface, resulting in some firn to remain temperate throughout the year as winter cold temperatures no longer penetrate the entire firn column. In the summers of 2014 and 2017, simulated meltwater does not reach deeper than approximately 
$3\!-\!5$ m below the surface, but temperate firn is preserved below 
$\sim10$ m from the surface.
According to the simulation results, seasonal aquifers form in 2013, 2016-17 and 2020-22, when water pools above the glacier ice, and the firn is saturated. In the summers of 2020 and 2022, the top of the saturated firn is shallowest at 8 m and 5 m below the surface, respectively. The depth of the water table is highest in summer and decreases toward winter.
Densities range from 298 kg m−3 to 917 kg m−3, increasing with depth (Fig. 5e). The annual mean density averaged over the top 10 m increases significantly between 2009 and 2022 (p-value < 0.01), rising from 524 kg m−3 to 711 kg m−3.
Cold content peaks in 2011, followed by a subsequent decline; however, no statistically significant trend is observed (Fig. 8a). The year 2011 is the coldest in the simulation, with a mean annual temperature of 
$-5.7^\circ\mathrm{C}$ in the top 10 m, and the firn column is cold throughout below the summer melt. Cold content reaches a minimum in 2016, which is the warmest year in the simulation, with a mean annual temperature of 
$-1.8^\circ\mathrm{C}$ in the top 10 m, associated with temperate conditions below the penetration depth of seasonal cold.
The highest annual amounts of refreezing in the top 10 m occur in 2013 and 2017, with values of 0.48 m w.e., and 0.42 m w.e., respectively (Fig. 8b). However, no statistically significant trend is detected. Sharp increases in refreezing are related to more available meltwater. In contrast, simulated pore space decreases significantly (p-value < 0.01) during this period (Fig. 8c). The decrease in pore space is connected to increases in density within the same period (Fig. 5e).
Modeled annual (a) cold content, (b) refreezing, and (c) firn pore space in the top 10 m from 2009 to 2022. Annual values are derived from averaging (a, c) and summing (b) daily model output. The dashed line is the linear fit (only shown when significant; p-value < 0.05). The p-value is the probability value at a 95% confidence level, and r 2 is the coefficient of determination.
5. Discussion
5.1. Temporal changes in firn density and stratigraphy
To investigate temporal changes in firn stratigraphy, we examine variations in average density and ice content across three sections (0–3 m, 3–6 m and 6–9 m from the surface) among all available firn cores near the summit (Fig. 9). No clear temporal trends can be detected. However, in all sections, density appears lower between 2008 and 2015 compared to before and after, with ice content in the middle sections following this pattern. Core J_99 exhibits exceptionally high ice content in the bottom section, reaching 62%. Comparing cores drilled over the years by various field teams is challenging due to the lack of international standardized measurement protocols for retrieving and measuring firn core properties. Additionally, ice content and density can vary considerably on small scales (within meters) (Marchenko and others, Reference Marchenko2017a; Xiao and others, Reference Xiao2022), which may impact the representativeness of a single core from a given year. Nevertheless, we find that our observed variations in firn stratigraphy are consistent with variations in melt.
Mean density (a-c) and ice content (d-f) from 1998 to 2022 for the top (0-3 m, 11 and 7 cores for density and ice content, respectively), middle (3-6 m, 11 and 7 cores) and bottom sections (3-6 m, 8 and 6 cores) of the firn cores near the summit area (Figure 1). Core CampXV_58 is not included in the figure: the mean density for the top section is 448 kg m−3, and the ice content for the top and middle sections are 3.1% and 15.8%, respectively. The sections are depths below the snow surface, i.e. snow is included in the top sections. Glacier ice is not included in any of the cores; the bottom sections of cores J_99 and X_99 only include depths until −8.7 and −8.8 m, respectively.
To explore mechanisms behind the observed and simulated variations in firn properties, we calculate annual positive degree day (PDD) sums, defined as the sum of mean daily air temperatures above the melting point (Cogley and others, Reference Cogley2011), as a proxy for melt (Fig. 10). We use the CARRA dataset for site Camp from 1993 to 2022, which is the same dataset used to force the model simulations. The mean PDD sum over this 30-year period is 37 K d. Between 2005 and 2012, the PDD sum is approximately at or below average, ranging from 7 to 39 K d, indicating an 8-year period with little melt. From 2013 onward, several years – such as 2013, 2015, 2020 and 2022 – exhibit PDD sums substantially higher than the 30-year average.
Positive degree day sum from 1993 to 2022 based on mean daily air temperature from the CARRA forcing at site Camp. The dashed line indicates the mean of the 30-year period.
We find significant correlations between firn core density in the middle (3-6 m) and bottom (6-9 m) sections and the PDD sum of the year preceding core retrieval (p-values < 0.01, Fig. 11). This indicates deep meltwater percolation and refreezing, which is further supported by the temperature data. Intensive melt years lead to decreased pore space; conversely, the period of approximately average or below average melt from 2005 to 2012 may have caused the decrease in firn density and ice content in cores collected from 2008 to 2013 compared to those from 1998 to 2004. This indicates that firn pore space can replenish when less meltwater infiltrates and refreezes in the firn. The increased density in 2022 may be related to the high melt year of 2020. Our findings that firn density is affected by changes in annual melt are consistent with firn core studies from other regions (e.g. Braithwaite and others, Reference Braithwaite, Laternser and Pfeffer1994; Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013; Machguth and others, Reference Machguth2016; Rennermalm and others, Reference Rennermalm2022). To explore the longer-term impact of melt on firn density, we also correlate density for each section with mean PDD sums from two to five years prior to core retrieval (Fig. S2). Results suggest that for the middle section, melt conditions up to three years prior to core retrieval also affect the evolution of the stratigraphy. However, further data are needed to confirm these findings.
Positive degree day sum and mean density of the (a) top (0-3 m, 11 cores), (b) middle (3-6 m, 11 cores), and (c) bottom (6-9 m, 8 cores) sections. The PDD sum is from the year prior to core retrieval. The dashed line is the linear fit (only shown when significant; p-value < 0.05). The p-value is the probability value, and r 2 is the coefficient of determination. The significance is evaluated at a 95% confidence level.
Our examination of the firn stratigraphy reveals an absence of ice slabs (> 1 m ice layers) in the summit area of Austfonna, in contrast to recent observations from Greenland (e.g. Machguth and others, Reference Machguth2016; MacFerrin and others, Reference MacFerrin2019; Rennermalm and others, Reference Rennermalm2022) and the Canadian Arctic (e.g. Zdanowicz and others, Reference Zdanowicz2012; Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013; Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013). In Greenland, ice slabs have expanded and thickened since 2000 (MacFerrin and others, Reference MacFerrin2019; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023), and they have been shown to form rapidly after an extreme melt season (Culberg and others, Reference Culberg, Schroeder and Chu2021). Here, even extreme melt years like 2020 (Figs. 2, 10) do not considerably alter the firn stratigraphy so that thick impermeable layers form. Hence, the effect of ice slabs in accelerating runoff and mass loss by limiting percolation has likely not been relevant at our sites.
5.2. Firn thermal regime
Our firn temperature observations and modeling results indicate that between 2009 and 2022, the thermal regime of the firn at depth shifted from cold to temperate conditions throughout the year. Both observations and model simulations indicate that this thermal transition most likely occurred in 2013, suggesting that a single extreme melt year can induce a lasting change in the firn thermal regime. We note that in some years, simulated meltwater percolation is considerably shallower than the observed depths, and subsurface temperatures generally have a cold bias (Fig. 5c). The absence of a parameterization for deep percolation due to preferential flow paths in firn (Marchenko and others, Reference Marchenko2017b) could at least partially account for the observed cold bias below the wetting front. Nevertheless, the model captures the observed regime shift, indicating that this shift is a robust feature.
Similar shifts in the thermal state of deeper firn, transitioning from cold to temperate year-round, have been inferred based on temperature measurements in central Svalbard comparing data from 1965 and 1976 to data from 2012 and 2013 (Marchenko and others, Reference Marchenko2017a), as well as through simulations covering period 1957–2014 (Østby and others, Reference Østby, Schuler, Hagen, Hock, Kohler and Reijmer2017). However, our study is the first in Svalbard to document a thermal shift based on continuous multi-year observations further supported by firn modeling.
Firn temperature observations in the summit area of Austfonna are scarce before 2009. Schytt Reference Schytt(1964) reported a temperature of 
$-4.2^\circ\mathrm{C}$ at a depth of 7 m below the surface at Camp XV in 1958. Watanabe and others Reference Watanabe(2001) found firn temperatures of 
$-1^\circ\mathrm{C}$ and 
$-2.8^\circ\mathrm{C}$ at −10 m depth at site J in 1998 and 1999, respectively. Zagorodnov and others Reference Zagorodnov, Sinkevich and Arkhipov(1989) observed predominantly negative temperature profiles in 1987; however, in two boreholes, temperatures increased with depth, approaching the melting point at approximately 
$6\!-\!7$ m below the surface.
5.3. Firn aquifer
The model simulations indicate that, following the shift in thermal regime, firn aquifers formed above the simulated firn-ice interface in 2013, 2016-17 and 2020-22. We attribute the aquifer formation to year-round temperate firn and deep meltwater percolation, the latter facilitated by the absence of ice slabs. The aquifer appears to be increasing in duration and thickness over time. The intermittent depletion and thinning of the aquifers may result from refreezing and/or drainage. Refreezing appears to have contributed to the depletion of the thin aquifer in 2013, with refreezing peaking notably in 2013 (Fig. 8). However, when the firn above the ice-firn interface remains temperate from 2016 onward, thinning and depletion are mainly attributed to lateral drainage, which is governed by the runoff time scale (see Supplementary equation S10). We note that the model does not account for englacial water transport, such as firn aquifer water draining through cracks or crevasses (Poinar and others, Reference Poinar2017; Cicero and others, Reference Cicero2023), which would contribute to decreased aquifer thickness.
Arkhipov and others Reference Arkhipov(1987) and Zagorodnov and others Reference Zagorodnov, Sinkevich and Arkhipov(1989) reported evidence of year-around water in firn at the summit of Austfonna in the 1980s, inferred from borehole drillings in the summers of 1984, 1985 and 1987. However, high spatial variability was observed, with cold firn detected near boreholes with water incursions, indicating water pooling in local depressions. It remains unclear whether these local firn aquifers persisted over time or were episodic, short-lived phenomena facilitated by the abnormally warm summers in 1985 and 1986 (Zagorodnov and others, Reference Zagorodnov, Sinkevich and Arkhipov1989). However, at least at site Camp, our observations and simulations from 2009 to 2022 indicate a general transition from cold to temperate conditions, along with the emergence of firn aquifers.
The formation of perennial firn aquifers is highly dependent on climatic drivers, particularly accumulation rate and liquid water input. Previous studies simulating conditions conducive to the formation of perennial firn aquifers have primarily focused on Greenland (Forster and others, Reference Forster2014; Munneke and others, Reference Munneke, SRM, van den Broeke, van Angelen and Forster2014; Miège and others, Reference Miège2016; Brils and others, Reference Brils2024). They conclude that accumulation must be high enough to provide the pore space needed for water storage while preventing refreezing, and that melt and rain must be sufficient to enable aquifer formation. It has been suggested that an aquifer region could turn into an ice slab region when liquid water input increases (Brils and others, Reference Brils2024). Further research is needed to better understand the conditions controlling the formation and extent of firn aquifers in Svalbard glaciers, where high-elevation accumulation zones typically have temperate deep firn (van Pelt and others, Reference van Pelt2019). A strong easterly-westerly climatic gradient in Svalbard (Schuler and Østby, Reference Schuler and Østby2020) and, for example, the largely asymmetric accumulation pattern on Austfonna (e.g. Dunse and others, Reference Dunse2009) would be interesting study targets for examining the impact of climatic conditions on the extent of aquifers.
5.4. Refreezing
Previous Svalbard-wide modeling studies report decreasing trends in refreezing, mainly attributed to a reduction in firn area (Østby and others, Reference Østby, Schuler, Hagen, Hock, Kohler and Reijmer2017; van Pelt and others, Reference van Pelt2019; Noël and others, Reference Noël2020). In contrast, our study simulates refreezing at a single location, where changes reflect available pore space, cold content, and surface water input, independent of glacier-wide firn area reductions. We do not find a statistically significant trend in refreezing in the top 10 m firn layer between 2009 and 2022, despite a decreasing trend in pore space (Fig. 8). This suggests an increase in meltwater and rain that can refreeze due to a sufficient build-up of cold content during winter. However, refreezing is expected to decrease in the future at our simulation site as firn thickness, cold content, and available pore space decline due to increased melt and negative annual mass balances (van Pelt and others, Reference van Pelt WJJ, Schuler, Pohjola and Pettersson2021). Our results indicate that the turning point at which refreezing decreases has not yet been reached in the summit area of Austfonna. However, interannual variability in refreezing is large (Noël and others, Reference Noël2020), and Schuler and others Reference Schuler(2020) suggest that the region-wide turning point might not be as clearly discernible in Svalbard glaciers as in the peripheral glaciers of Greenland (Noël and others, Reference Noël2017).
6. Conclusions
Sixteen firn cores drilled on the Austfonna ice cap between 1958 and 2022 are characterized by frequent ice layers, even in the summit area, but ice slabs (> 1 m ice layers) that are found in other regions, such as western and northern Greenland (MacFerrin and others, Reference MacFerrin2019) and the Canadian Arctic (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013) are absent. Ice content and density appear generally lower in core sections 3-6 m below the surface between 2010 and 2015 compared to earlier and later, consistent with lower melt amounts inferred from positive degree-day calculations. However, no clear trends in time or elevation are observed.
Firn temperature observations near the summit between 2009 and 2022 uncover a shift in the thermal regime of the firn. In 2013, the firn above the firn-ice interface transitioned from cold to temperate year-round. This suggests enhanced warming of the deeper firn layers due to increased latent heat release from refreezing, with insufficient winter cooling to offset it. We used the CryoGrid community model, validated against in situ surface mass balance, temperature, and density data, to simulate firn temperature, water content, and density. The model simulations support the observed thermal shift in 2013 and also the repeated formation of firn aquifers since then.
Direct observations such as firn drilling and ground-penetrating radar surveys would be beneficial for confirming the inferred presence of a firn aquifer at site Camp, as well as determining its thickness and extent. Additionally, direct observations of firn thickness are necessary to improve model simulations. Further investigations are also needed to assess the impact of these firn aquifers on the mass balance of the ice cap.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10072.
Author contributions
SI gathered, compiled, quality-controlled and curated the datasets, performed the data analyses and simulations, and generated the figures. SI, RH and TVS initiated and designed the study. TVS, and, in particular RH, provided overall guidance throughout this work. SI interpreted the results with significant input from RH, FC and TVS, and wrote the manuscript with contributions from all co-authors, in particular RH. LSS, FC and TVS assisted with figure visualization. LSS adapted the CryoGrid model code in direct collaboration with SI and assisted considerably with the simulations. SI led the firn core fieldwork in 2022 with support from GM. TVS and GM collected and provided other field data. All authors commented on and contributed to the clarity of the manuscript.
Acknowledgements
Fieldwork on Austfonna has been supported by ESA project CryoVEx (PEA no 4000110725), the Research Council of Norway projects AltiMelt (grant 251355/F50) and SIOS-InfraNor (grant 269927). Further work has been supported by the Research Council of Norway through the Nansen Legacy project (grant 276730).
SI was supported by Maa- ja Vesitekniikan Tuki Ry (scholarship 4438). RH was supported by the Research Council of Norway (project 324131) and ERC-2022-ADG grant 101096057 (GLACMASS).
We thank the Norwegian Polar Institute for the logistic support and Eirik Hellerud for assistance in the field in 2022. We are grateful to Takao Kameda and Hideaki Motoyama for providing previously published firn core data in digital form. Jon Ove Hagen, Elisabeth Isaksson, and Lisbeth Langhammer provided valuable information in locating previously unpublished firn data. Veijo Pohjola gave insightful comments on this work. We thank two anonymous reviewers and the editor Dr. William Colgan for their valuable contributions to the manuscript.
Data and software availability
The version of the CryoGrid code used for this study is available at https://doi.org/10.5281/zenodo.15838031. The firn core data are available at https://doi.org/10.21334/NPOLAR.2025.9A2A014B. The firn temperature data are available at https://doi.org/10.21334/NPOLAR.2025.10DB6623. The surface mass balance data are available at https://mosj.no/en/indikator/climate/land/mass-balance-for-glaciers-in-svalbard/.
Open access
