3.1. Field observations

Air temperature and SMB data for this study were collected through field observations. The Geological Survey of Canada (GSC) maintains an annual mass balance monitoring program on the northwest basins of the Devon Ice Cap (Koerner, Reference Koerner2005; Burgess, Reference Burgess2017). The northwest measurement transect (NW transect) extends from the ice cap summit to near sea level.

We utilized 2 m air temperature data obtained from three AWS located along the NW transect (Fig. 1b): AWS D1-1H (75.3691°N, 82.6702°W, and 1829 m a.s.l.), AWS D2-4F (75.3863°N, 82.76450°W, and 1737 m a.s.l.), and AWS D3-7D (75.4219°N, 82.9473°W, and 1638 m a.s.l.). The air temperature data were collected using Campbell Scientific 44212 temperature probes with an accuracy of ±0.1°C, housed inside a ventilated Gill radiation shield. Air temperature data from all three AWSs were available from 1999 to 2022 averaged an hourly time scale (with some gaps). This data was used to adjust biases in the reanalysis data applied in this study (see Section 3.2 for details).

Field observations along the NW Transect (Fig. 1b) were obtained from the monitoring program of the GSC (Koerner, Reference Koerner2005). Accurately accounting for refreezing processes within field-based SMB monitoring remains a recognized challenge. In the GSC program, SMB is determined using the stratigraphic approach described by Cogley and others (Reference Cogley2011). This method involves measuring the annual vertical displacement between successive end-of-summer surfaces relative to a fixed survey pole. The measured height change is then converted to net annual SMB in water equivalent by applying the average density of the material accumulated or lost at the measurement site. Although internal refreezing is not directly incorporated in the calculations, it is assumed that a portion of the apparent summer mass loss above the firn line is retained within the ice cap, proportional to the firn density in the vicinity of the pole. Based on this methodology, we estimate a conservative uncertainty of approximately 14% in the GSC SMB pole measurements. We utilized net annual SMB data from 14 locations along the NW transect of the study area (Fig. 1b) for each mass balance year from 2004 to 2021. The annual SMB data is utilized only for validating the model results. In this study, the mass balance year is defined as starting on September 1st of a given calendar year and ending on August 31st of the following year (e.g. Tedesco and others, Reference Tedesco2013; Strand and others, Reference Strand, Christiansen and Gilbert2022). Henceforth, all annual quantities will adhere to this definition of the mass balance year.

We used field data on vertical firn density profiles derived from shallow firn cores that were collected at three sites: Colgan E (75.4494° N, 82.5305° W, and 1525 m a.s.l.), Colgan C (75.3399° N, 82.6763° W, and 1825 m a.s.l.), and Colgan B (75.3427° N, 83.3525° W, and 1630 m a.s.l.), across the Devon Ice Cap (Colgan and Sharp, Reference Colgan and Sharp2008; Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013). These data were collected in 2004/05 and 2012, and we used them to evaluate the corresponding model results. Additionally, we used the mean bulk firn density of the top 2.5 m of the firn layer from 10 locations along the CryoSat line (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a) (Fig. 1b) to validate our model results. This data was available for the years 2004, 2006, 2008, 2011, and 2012 (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013). We also used in-situ mean bulk firn densities from six shallow firn cores collected in 2012 (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013) and 2022 (Hallé and others, Reference Hallé2025) in southwest Devon Ice Cap and compared them with the corresponding model simulations.

3.2. Reanalysis data and bias corrections

We used 2 m air temperature, incoming shortwave radiation, cloud cover and 2 m dew point temperature from the ERA5-Land dataset (Muñoz-Sabater and others, Reference Muñoz-Sabater2021) from 1999 to 2022, with hourly time resolution. The relative humidity (RH) is not directly available in the ERA5-Land dataset. We compute the hourly RH by determining the ratio of actual vapour pressure to saturation vapour pressure, using 2 m air temperature (T_air) and dew point temperature (T_d) from the ERA5-Land dataset, according to the following relation (Te Chow and others, Reference Te Chow, Maidment and Mays1988):

(1)\begin{equation}RH = 100 \times \frac{{0.6108 \times {\text{exp}}\left( {\frac{{17.625\,{T_d}}}{{243.04 + {T_d}}}} \right)}}{{0.6108{ } \times {\text{exp}}\left( {\frac{{17.625\,{T_{air}}}}{{243.04 + {T_{air}}}}} \right)}}.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\end{equation}

As cloud cover data is not provided in the ERA5-Land, we used the total hourly cloud cover data from the ERA5 dataset (Hersbach and others, Reference Hersbach2020). Additionally, the precipitation data from the ERA5-Land dataset shows substantial biases, with differences in the order of magnitude when compared to the typical observed precipitation over the Devon Ice Cap (Koerner, Reference Koerner1966; Shepherd and others, Reference Shepherd, Du, Benham, Dowdeswell and Morris2007; Noël and others, Reference Noël, van de Berg, Lhermitte, Wouters, Schaffer and van den Broeke2018). Therefore, we utilized precipitation data from the North American Regional Reanalysis (NARR) (Mesinger and others, Reference Mesinger2006), following the previous literature (e.g. Gardner and Sharp, Reference Gardner and Sharp2009; Gascon and others, Reference Gascon2014). However, we chose to retain ERA5-Land for 2 m air temperature, humidity, and other meteorological variables due to its higher spatial resolution (9 km vs. NARR’s ∼32 km) and hourly temporal resolution, which better match the requirements of our modelling framework (Section 3.3).

We compared the hourly air temperature data from ERA5-Land for three AWS locations (Fig. 1b) to estimate the biases over the period from 1999 to 2022. To perform the temperature comparison, we corrected for elevation between the corresponding ERA5-Land grid and AWS locations (Section 3.1) by applying a constant lapse rate of −5.7°C km⁻¹ (Ashmore and others, Reference Ashmore, Mair and Burgess2020). At each AWS location, a time series of annual temperature biases for ERA5-Land was estimated from 1999 to 2022, with some gaps due to the availability of observed data (Fig. S1). The mean annual temperature biases for the three AWS locations were calculated and assumed to remain constant across the entire study area during the model simulation period.

3.3. The model

We employed a 1-D physically-based mass-energy balance model (Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014) over the accumulation area of the Devon Ice Cap (elevation ≥ 1300 m a.s.l., Fig. 1b) to simulate the evolution of subsurface snow and firn. The simulation used a horizontal spatial resolution of 5 km × 5 km and a vertical resolution of 1 cm, with a temporal resolution of 15 minutes over the period from 1999 to 2022. The high resolution of the vertical grid allows us to continuously model the growth and decay of ice layers ranging from 1 cm to several metres thick, a capability that is not possible with the coarser vertical resolutions typical of many existing models in the literature (Vandecrux and others, Reference Vandecrux2020). Each spatial grid represents a 1-D vertical multilayer structure of snow and firn, characterized by its corresponding temperature and density profiles. The modelling scheme consists of two primary components: a surface energy balance module and a vertical percolation–refreezing module. The model uses air temperature, precipitation, incoming shortwave radiation, relative humidity (RH), and cloud cover data as forcing variables to parameterize (Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014) the components of the surface energy balance (Hock, Reference Hock2005; Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014). The energy balance equation is expressed as the sum of its constituent components (Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014):

(2)\begin{equation}{Q_t} = G\left( {1 - \alpha } \right) + L + {Q_h} + \,{Q_l}.\,\end{equation}

here G, L, and α denote the incoming shortwave radiation (W m⁻²), the net longwave radiation balance (W m⁻²), and the surface albedo, respectively. Q_h and Q_l represent the sensible and latent heat fluxes (W m⁻²), respectively. Q_t corresponds to the total surface heat flux (W m⁻²), which contributes to surface temperature changes, melting, and the conductive heat flux into the subsurface. Following Greuell and Konzelmann (Reference Greuell and Konzelmann1994) and Morris and others (Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014), the longwave radiation balance and the sensible and latent heat fluxes were calculated using analogous parameterizations. The surface albedo (α) was parameterized as a linear function of the density of the uppermost 10 cm of the snowpack, following Morris and others (Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014):

(3)\begin{equation}\alpha = {\alpha _i} + \left( {{\rho _{10}} - {\rho _i}} \right)\left( {\frac{{{\alpha _s} - {\alpha _i}}}{{{\rho _s} - {\rho _i}}}} \right)\end{equation}

here ρ ₁₀ denotes the mean density of the uppermost 10 cm (kg m⁻³). The constants in the Eqn (3) take the following values: bare ice albedo (α_i) = 0.50, albedo of freshly fallen snow (α_s) = 0.81, ice density (ρ_i) = 910 kg m⁻³, and density of freshly fallen snow (ρ_s) = 410 kg m⁻³. We also compared the simplified albedo parameterization (Eqn (3)) with the high-resolution Copernicus Arctic Regional Reanalysis (CARRA) dataset (accessible at https://climate.copernicus.eu/copernicus-arctic-regional-reanalysis-service). The modelled interannual variability in albedo shows a good match with that obtained from the CARRA data over the study region (Fig. S2).

Snowfall or rainfall was determined using a 1°C threshold on 2 m air temperature. Snow was added to the top, increasing its thickness and adjusting its bulk density with fresh snow density ρ_s; rainfall percolated immediately. Meltwater generation was computed from the integral of Q_t (Eqn (2)) over time, with Q_t first warming the surface to 0°C before inducing melt. Note that the model excludes snow/firn compaction and meltwater flow to simplify complexity and enhance computational efficiency. We implemented a commonly used bucket scheme to simulate vertical meltwater percolation and refreezing (e.g. Simonsen and others, Reference Simonsen, Stenseng, Ađalgeirsdóttir, Fausto, Hvidberg and Lucas-Picher2013; Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014; Verjans and others, Reference Verjans, Leeson, MacFerrin, Noël and van den Broeke2019). Here, meltwater and rainfall are instantaneously infiltrated from the surface, percolating through successive layers and refreezing at depth according to the available cold content. The latent heat released during refreezing can sometimes increase the layer temperature to the melting point, initiating deeper percolation and prompting the model to update the layer temperature accordingly. Percolation continued until all liquid refroze or encountered an impermeable ice layer, at which point the excess became runoff (e.g. Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014; Ashmore and others, Reference Ashmore, Mair and Burgess2020). The novel development to the model used in this study is our application and evaluation of three different parameterizations to represent the permeability of subsurface ice layers (Section 3.5). In each model grid, the net annual SMB (unit m w.e. a ⁻¹) was estimated using the following relation,

(4)\begin{equation}SMB=Precipitation-Surface\;melt+Total\,refreezing.\,\,\,\,\end{equation}

3.4. Model implementation

We note that the horizontal resolution of the model differs from that of the forcing reanalysis data. Thus, all forcing variables for each model grid were extracted from the corresponding reanalysis data (Section 3.2) using a nearest-neighbour approach. The air temperature shows a strong dependence on elevation compared to the other forcing variables, and the mean elevation of a model grid differs from that of the nearest-neighbour reanalysis data grid. Therefore, the air temperature data were corrected for elevation differences between the model grid and the ERA5-Land grid using a constant lapse rate of −5.7°C km⁻¹ throughout the simulation period (Ashmore and others, Reference Ashmore, Mair and Burgess2020). The mean elevation of each model grid was derived from the ArcticDEM (version 4.1, Mosaic 2 m) (Porter and others, Reference Porter2023). Additionally, all hourly forcing data were linearly interpolated to the model’s time step of 15 minutes. In this study, we defined the vertical model domain to span from +15 m to −15 m relative to the initial surface position. The model updates the surface position after each mass balance year, based on the balance between annual accumulation and melting. Here, we ignored the surface elevation changes arising from ice dynamics, since we have calculated that this contribution is negligible (see the supplementary section S1 for details). We initialize the vertical density and temperature profiles at each model grid by following the procedure described by Ashmore and others (Reference Ashmore, Mair and Burgess2020). The model’s initial firn density profile is derived from measurements collected in 2000/01 (Mair and others, Reference Mair, Burgess and Sharp2005) by fitting an exponential function to the observed density–elevation relationship. The initial temperature profile is similarly parameterized by fitting an exponential curve between the snow surface and 10 m depth. Here, the 10 m temperature is assumed to represent the mean air temperature over the preceding 12 months, while the surface temperature is set to the mean May air temperature. Sensitivity tests (Ashmore and others, Reference Ashmore, Mair and Burgess2020) indicate that these assumptions have only a minor influence on the resulting density stratigraphy. Then, we used five years of forcing data from 1994 to 1998 for the model spin-up, followed by running the model from 1999 to 2022. To verify the robustness of the results, we perturbed the vertical density and temperature profiles by 25% during the spin-up runs and confirmed that the final model outcomes from 1999 to 2022 were not influenced by the initial conditions.

3.5. Methods for defining ice layer permeability

Although field observations and laboratory experiments indicate that ice layers can be permeable under certain conditions (e.g. Harper and others, Reference Harper, Humphrey, Pfeffer, Brown and Fettweis2012; Sommers and others, Reference Sommers, Rajaram, Weber, MacFerrin, Colgan and Stevens2017; Fowler and Iverson, Reference Fowler and Iverson2023), modelling studies often assume them to be impermeable (e.g. Sørensen and others, Reference Sørensen2011; Kuipers Munneke and others, Reference Kuipers Munneke, Ligtenberg, van den Broeke, van Angelen and Forster2014; Vandecrux and others, Reference Vandecrux2020). These modelling approaches commonly employ uneven grid spacing with depth and various meltwater percolation–refreezing schemes to represent ice layers within firn and snow. Ashmore and others (Reference Ashmore, Mair and Burgess2020) were the first to introduce variable thicknesses of impermeable ice layers (0.01–5 m) to evaluate model performance against field observations. A greater thicknesses of impermeable ice layers (>1 m) allows more meltwater to percolate and refreeze within the firn before runoff occurs, resulting in higher retention and reduced runoff. This model criterion was qualitatively able to reproduce the observed near-surface thin ice layers in the accumulation zone of the Devon Ice Cap (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a). Although varying impermeable ice layer thickness values between 0.25 and 1 m had little impact on the SMB comparison between observed and modelled values.

Our new prioritisation of criteria that determine ice layer permeability follows from a suite of cold laboratory experiments that have undertaken (Shiggins and others, Reference Shiggins Connor, Mair Douglas, Ashmore David, Nias Isabel, Lea James and Brown2025) to observed and monitored the effective permeability of predefined ice layers within snowpacks injected with dyed meltwater to simulate surface melt within two defined thermal regimes: a) cold and b) temperate (Shiggins and others, Reference Shiggins Connor, Mair Douglas, Ashmore David, Nias Isabel, Lea James and Brown2025). When the thermal regime is cold (i.e. snow/firn temperature immediately below the ice layer is <−1°C) all ice layers were impermeable and forced meltwater to pond, runoff and refreeze. Even the thinnest ice layer (5 mm), exposed to extreme meltwater fluxes, provided an impermeable barrier to meltwater and was resilient to degradation and fracturing. When the snowpack was temperate (i.e. snow/firn temperature immediately below the ice layer is >−0.15°C), ice layers of all thicknesses tested (10–60 mm) were permeable and allowed meltwater to percolate to deeper depths. For the thinnest ice layer (10 mm), permeability was due to small fractures forming in the ice layer beneath the ponded meltwater, while all the thicker layers experienced meltwater seeping into the ice and being transported through the layer via water veins (Fowler and Iverson, Reference Fowler and Iverson2023). Temperate ice layers became permeable in as little as ∼6–8 hours (e.g. through layers 40–50 mm thickness), significant within a diurnal cycle (Shiggins and others, Reference Shiggins Connor, Mair Douglas, Ashmore David, Nias Isabel, Lea James and Brown2025). Consequently, Shiggins and others (Reference Shiggins Connor, Mair Douglas, Ashmore David, Nias Isabel, Lea James and Brown2025) suggest that temperature, not thickness, of an ice layer is a first-order control on its permeability.

In this study, we tested three different methods for ice layer permeability, which regulate meltwater percolation and runoff availability. Method 1 stipulates that if the temperature immediately below the ice layer (T_th) is warmer than −0.15°C, the layer is considered permeable up to an ice layer thickness of 1 m (H_imp). Here, the subscripts ‘th’ and ‘imp’ denote the threshold temperature value and the impermeable ice layer, respectively. Method 2 specifies that if the temperature below the ice layer (T_th) is warmer than −0.15°C, the layer is permeable irrespective of the ice layer thickness. For both Methods 1 and 2, if T_th is cooler than −0.15°C the ice layer is impermeable. Method 3 follows the same approach as described by (Ashmore and others, Reference Ashmore, Mair and Burgess2020) with H_imp = 1 m and no thermal criteria for permeability. The rationale for adopting a standard H_imp value of 1 m was informed by the standard definition of an ‘ice slab’ in previously observed GPR data (e.g. Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023), and following previous comparisons of modelling outcomes to observations over the Devon Ice Cap (Ashmore and others, Reference Ashmore, Mair and Burgess2020) and across Greenland (Brils and others, Reference Brils2024). The above framework enabled us to identify which permeability criteria compared most favourably with measurements of near-surface density and SMB, and therefore which most effectively simulate percolation, refreezing and runoff, yielding more realistic estimates.

3.6. Model validation and parameter sensitivity

We validated the modelled annual SMB, firn density profiles, and firn bulk density (upper 2.5 m) from the three different methods of defining ice layer permeability (Section 3.5) against corresponding field observations (Section 3.1). The performance of each method was evaluated using root mean squared error (RMSE), R-squared (R²), and Pearson correlation (r) values. This validation procedure was used to identify the method with the overall most effective performance.

In this study, we did not calibrate the energy-mass balance model parameters. Instead, we use the previously field-based calibrated model parameters from the CryoSat line (Fig. 1b) on the Devon Ice Cap (Morris and others, Reference Morris, Mair, Nienow, Bell, Burgess and Wright2014; Ashmore and others, Reference Ashmore, Mair and Burgess2020). We conducted a model sensitivity test on the parameters (T_th and H_imp) that control ice layer permeability for the selected the most effective method. We performed the sensitivity test only at the locations for which where we had field data for validation (Section 3.1).

3.7. Climate response of mass balance and summer runoff

We utilized the interannual variability of modelled specific summer runoff (m w.e. a⁻¹), annual mass balance (m w.e. a⁻¹), annual precipitation (m w.e. a⁻¹), and summer air temperature (°C) to assess the corresponding climate response at each model grid. Since runoff is mostly available only during the summer (>90%), we considered only summer runoff for the analysis. In each model grid, the climate response of annual SMB and summer runoff (Q) is related to annual precipitation (P) and summer temperature (T) by (e.g. Zheng and others, Reference Zheng, Zhang, Zhu, Liu, Sato and Fukushima2009; Banerjee, Reference Banerjee2022; Laha and others, Reference Laha, Banerjee, Singh, Sharma and Thamban2023):

(5)\begin{equation}\delta SMB = s_{SMB}^P\delta P + s_{SMB}^T\delta T,\end{equation}

(6)\begin{equation}\delta Q = \,s_Q^P\delta P\, + \,s_Q^T\,\delta T.\,\,\,\,\,\end{equation}

here the δ notation represents the anomaly of the corresponding variables. The variables s^P and s^T represent the precipitation sensitivity (m a⁻¹ °C⁻¹) and temperature sensitivity (m a⁻¹ °C⁻¹) of summer runoff (Q) and annual SMB. We define precipitation sensitivity as the change in SMB or runoff resulting from a 10% change in precipitation, consistent with conventions used in the literature (e.g. Oerlemans and Reichert, Reference Oerlemans and Reichert2000; Woul and Hock, Reference Hock2005; McGrath and others, Reference McGrath, Sass, O’Neel, Arendt and Kienholz2017). The potential bilinear interaction term δP × δT (Lang, Reference Lang, Kraijenhoff and Moll1986) in the linear climate response model (Eqs. (5) and (6)) was not included because the correlation between δP and δT was not statistically significant (r = 0.07, p > 0.01) for all the model grids over the Devon Ice Cap.

4.1. Parameter sensitivity

We adjusted the H_imp value to 0.25 m, 1.5 m, and 2 m in Method 1, with T_th set to −0.15°C, and ran simulations across the study area to assess the model’s sensitivity to H_imp. These adjustments increased the RMSE between the observed and modelled SMB data by 8%, 11%, and 18%, respectively, compared to the results with H_imp set at 1 m. Furthermore, the RMSE for the firn density profiles increased by more than 40% across all three sites, while the RMSE for the top 2.5 m mean bulk firn densities increased by over 60%. These results suggest a high sensitivity of the model to H_imp. Consequently, we opted to keep the parameter’s optimal value at 1 m. We further confirmed that the sensitivity of H_imp in Methods 2 and 3 fell within a similar range as reported above.

We further tested T_th values of −0.5°C, −1°C, and −2°C, while maintaining H_imp at 1 m in Method 1. In these cases, the maximum increase in RMSE for SMB, firn density profiles, and top 2.5 m mean bulk firn density was less than 7% when compared to the results at T_th = −0.15°C. This suggests the model has a relatively low sensitivity to variations in T_th values. Thus, we retained the optimal T_th value of −0.15°C, as supported by the laboratory experiments (Shiggins and others, Reference Shiggins Connor, Mair Douglas, Ashmore David, Nias Isabel, Lea James and Brown2025).

4.2. Model validation

For model validation (as with Model implementation in Section 3.4), we utilized a nearest-neighbour model grid for each point-scale near-surface firn bulk density (top 2.5 m), firn density profile, and SMB measurement across the study area for three different methods. Here, Method 1 assumes the ice layer is permeable if the temperature below it (T_th) warmer than −0.15 °C and permeable up to an ice layer thickness (H_imp) 1 m. Method 2 applies the same temperature threshold but considers the ice permeable regardless of thickness. Method 3 adopts the approach of Ashmore and others (Reference Ashmore, Mair and Burgess2020) with no thermal criteria and setting H_imp at a value of 1 m. We then identified the most effective method based on validation performance.

4.2.1. Near-surface mean bulk firn density

The observed average bulk firn density measurements for the top 2.5 meters were collected from 10 different locations along the CryoSat line, covering an elevation range of 1400 to 1685 m a.s.l. (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013). These point-scale measurements were available for the years 2004, 2006, 2008, 2011, and 2012, with some gaps in the data. The average bulk firn densities for the top 2.5 meters exhibit an overall decreasing trend with elevation, as observed in both the models and the observations (Fig. 2). This trend highlights the densification of the near-surface firn layer at lower elevations. We found that the model density values using Method 3 were overestimated compared to the other two methods (Fig. 2 and S3). However, we did not observe any significant (at p < 0.01 level) trend in the modelled density biases relative to the corresponding elevations (Fig. S3). For all data points, the RMSE values between the modelled densities and the corresponding observations were 47.6, 47.8, and 60.8 kg m⁻³ for Method 1, Method 2, and Method 3, respectively. The corresponding R ² values were 0.8, 0.8, and 0.6, while the r values were 0.9, 0.9, and 0.8, all statistically significant at the p < 0.01 level. These results indicate that the performance of Methods 1 and 2 was similar and that both relatively outperformed Method 3.

Figure 2. A comparison of the modelled mean bulk firn density of the top 2.5 m layer for three different methods against observations along the CryoSat line (Figure 1) at 10 different locations (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013). The symbols open circles, crosses, and upper triangles denote Method 1, Method 2, and Method 3, respectively. Each point is coloured according to the corresponding elevation. The solid magenta, violet, and red lines represent the best-fit lines for Methods 1, 2, and 3, respectively. The corresponding RMSE values are 47.6, 47.8, and 60.8 kg m⁻³, with R ² values of 0.8, 0.8, and 0.6, respectively. The dotted grey line represents the 1:1 reference line.

4.2.2. Vertical firn density profiles

We compared the modelled firn density profiles up to a 10-meter depth at three Colgan sites: C (1825 m a.s.l.), B (1630 m a.s.l.), and E (1525 m a.s.l.) on the Devon Ice Cap (Colgan and Sharp, Reference Colgan and Sharp2008; Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013) for the years 2012 (Fig. 3) and 2004/05 (Fig. S4). The vertical density profiles produced by Methods 1 and 2 were identical at all three Colgan sites (Figs 3 and S4). This similarity arises because no ice slabs exceeding 1 m in thickness developed at these locations, resulting in comparable performance of Methods 1 and 2 (Section 3.5). Note that 2012 was a warm year and 2004/05 was a relatively cool year. For brevity, herein we primarily report comparisons between modelled and observed measurements for the warm year 2012, while additional figure for 2004/05 is provided in the supplementary materials (Fig. S4). In 2012, Method 3 produced an ice layer more than 2 m thick at approximately 2 m depth, particularly at the Colgan B and E sites (Fig. 3b and 3c). Such a thick ice layer was absent in both the observations and in the modelled results from Methods 1 and 2. At the higher-elevation Colgan C site, Method 3 produced an ice layer nearly 2 meters thick at around 3 meters depth, whereas the observed data showed density variations ranging from 450 to 500 kg m⁻³ (Fig. 3a). In both 2012 and 2004/05, the modelled densities do not exactly match the observed values (Figs 3 and S4), which is to be expected given the comparison between point-scale core measurements and spatially averaged modelled values over areas of approximately 25 km². However, Methods 1 and 2 perform better than Method 3, particularly at the Colgan B and C sites in 2012. We also found that the modelled density profiles from Methods 1 and 2 matched better in the warm year 2012 than in the cool year 2004/05. In some cases, the model exhibited a slight positional lag in capturing the density peaks, which may be due to the spatial resolution of the model versus the resolution at which density was measured along the firn cores.

Figure 3. Comparison between the observed firn density profiles at three sites—Colgan C (a), Colgan B (b), and Colgan E (c) (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013), and the corresponding modelled profiles for three different methods during 2012. Note that Methods 1 and 2 performed identically across all three Colgan sites (see Section 4.2.2).

Across all three sites, Methods 1 and 2 yielded better RMSE values compared to the Method 3 (Table 1), with the best match and lowest RMSE observed at the Colgan B site for the year 2012. At Colgan C, B, and E, the RMSE values for Methods 1 and 2 improved by more than 33%, 38%, and 21%, respectively, compared to Method 3 (Table 1). We also verified our model performance by applying a 20 cm moving average to the observed and modelled firn density profiles (Fig. S5) to reduce the impact of existing noise in the data (Fig. 3) on the RMSE values. We found that Methods 1 and 2 outperformed the other method, achieving RMSE improvements of 38%, 41%, and 32% for the Colgan C, B, and E sites, respectively. In 2004/05, at the three Colgan sites, Methods 1 and 2 yielded the same RMSE, which was lower than that of Method 3 (Table S1). This highlights that, independent of the noise in the data, the modelled density profiles generated by Methods 1 and 2 provide better results at all three Colgan firn density measurement sites. Modelling firn density profiles sometimes exhibits substantial positive biases near the surface and negative biases in deeper sections when compared to observations (e.g. Gascon and others, Reference Gascon2014). In contrast, our current model mitigates these biases through higher vertical resolution (1 cm) and a temperature-based ice layer permeability criterion. These advancements facilitate the simulation of a more realistic evolution of ice layers and their subsequent impacts on runoff.

Table 1. RMSE values comparing the modelled and observed firn density profiles at three Colgan sites during 2012 (Figure 3) using three different methods.

4.2.3. Annual SMB data

We compared the observed annual SMB data collected from 14 locations along the NW transect between 2004 and 2022 (Fig. 1b) with the corresponding model results from the three different methods (Fig. 4). We found that Method 1 outperformed the other two methods, with R² values of 0.48, 0.21, and 0.24 for Method 1, Method 2, and Method 3, respectively (Fig. 4). The corresponding RMSE values were 0.09, 0.19, and 0.18 m w.e. a⁻¹. Furthermore, the r values were 0.69, 0.42, and 0.46, statistically significant at p < 0.01 level. Therefore, our newly introduced Method 1 substantially improves the R ², RMSE, and r values compared to the previously employed Method 3 (Ashmore and others, Reference Ashmore, Mair and Burgess2020). We observed a positive bias in the modelled SMB estimates for both Method 2 and Method 3, particularly for the observed negative SMB values (Fig. 4b and 4c) A similar positive bias in the modelled SMB for Method 3 was also observed by Ashmore and others (Reference Ashmore, Mair and Burgess2020). This bias was primarily attributed to substantial refreezing of percolating meltwater within the snow/firn layer in Methods 2 and 3, leading to the formation of a relatively thicker ice layers near the surface (Fig. S6). Method 1 more accurately captures the spatio-temporal variability of SMB across the NW transect (Fig. 1b).

Figure 4. Comparison of observed annual SMB data from 2004 to 2021 along the NW transect with the corresponding modelled values for Method 1 (a), Method 2 (b), and Method 3 (c). Here, each point is coloured according to the corresponding mass balance year, and the dotted grey line represents the 1:1 reference line. The solid red line in each plot denotes the best-fit line.

Our analysis of model validation using near-surface bulk firn densities and vertical firn density profiles indicates that Methods 1 and 2 performed better than the previously used Method 3 (Ashmore and others, Reference Ashmore, Mair and Burgess2020). Point-scale SMB validation further demonstrates that Method 1 outperformed the other two methods, yielding more consistent modelled SMB values relative to the corresponding field observations across the study area. This suggest that the temperature-based criteria of ice layer permeability improves the existing method in the literature (Ashmore and others, Reference Ashmore, Mair and Burgess2020). Consequently, Method 1 was selected for further analysis in this study.

4.3. Development of near-surface ice layers

In the model results, surface meltwater percolation and refreezing are represented by the formation of near-surface ice layers (within the top 10 m; Fig. 5), predominantly in the lower accumulation zone where warmer summer temperatures create wetter snow conditions. These processes are inferred from changes in firn density (Fig. 5), as the current model configuration does not explicitly simulate liquid water content profiles. Vertical near-surface density stratigraphy was observed using GPR data along the CryoSat line on the Devon Ice Cap during spring 2007, 2009, 2010, 2011, and 2012, covering elevations from 1440 to 1690 m a.s.l. (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a). Accordingly, we discuss the modelled near-surface density stratigraphy along this CryoSat line (Fig. 1b) to further evaluate the model results. The observed density stratigraphy indicated an increased thickness of near-surface ice layers in the lower accumulation zone (elevation range 1400–1500 m a.s.l.) compared to the upper accumulation zone (elevation > 1600 m a.s.l.) (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a). A similar pattern was also observed in our model simulations (Fig. 5).

Figure 5. Modelled depth-density profiles for the model grids along the CryoSat line (Figure 1b) during the spring seasons of 2007 (a), 2012 (b), 2017 (c), and 2022 (d), illustrating changes over time. The colour scale indicates the corresponding density (kg m⁻³) values.

Compared to the observed density stratigraphy, the model simulations allow for a more detailed discussion of the evolution of ice layers within the firn layer (Fig. 5). Our model results suggest that the observed thin near-surface ice layer at an elevation of approximately 1400 m a.s.l. during 2005 (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a) formed during the relatively warm summer of 2001 (Fig. S7). This ice layer was present throughout the CryoSat line up to 1800 m a.s.l. During spring 2007, the modelled depth of the above ice layer was between 6–7 m from the surface at elevations of around 1400–1600 m a.s.l., respectively (Fig. 5a), which was similar to the field GPR data (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a). By the year 2012, the near-surface ice layer had become thicker compared to year 2007, with GPR-based observed thickness measurements of 5.7 m, 3.8 m, and 0.7 m at elevations of 1400, 1490, and 1610 m a.s.l., respectively (Gascon and others, Reference Gascon, Sharp, Burgess, Bezeau and Bush2013a). The corresponding modelled thicknesses at elevations of approximately 1400–1500 m and 1600 m a.s.l. were about 3 m and 0.8 m, respectively (Fig. 5b). This indicated that our model estimates of ice layer thickness aligned well with observations at relatively higher elevations. However, the model tended to underestimate ice layer thickness at lower elevations. As the near-surface ice layer developed from 2005 onward, it prevented meltwater from percolating deeper, resulting in more meltwater being available for runoff. This was evident from the time series of the mean annual runoff over the study area, where we observed a slight increasing trend in runoff between the years 2005 and 2012 (Fig. S8).

In 2012, the number of thin ice layers increased at higher elevations (>1600 m a.s.l.) compared to 2007 (Fig. 5a and 5b), due to a relatively warm summer. However, the initial ice layer formed during the summer of 2001 was no longer present near the surface by 2012 (Fig. 5b), as it moved deeper (below 9 m) into the subsurface due to continuous accumulation. At elevations > 1600 m a.s.l., we observed patches of multiple thin ice layers ranging from 2 to 6 meters in depth in 2012 (Fig. 5b), which advanced to greater depths by 2017 (Fig. 5c). Additionally, new thin near-surface ice layers formed between 2018 and 2021, which were advected to approximately 4 meters depth by 2022 (Fig. 5d). However, the number of near-surface thin ice layers in 2022 was reduced compared to 2012.

Our model results revealed that the near-surface ice layer at lower elevations (∼1400 m a.s.l.) gradually thickened from 2007 to 2022 (Fig. 5). This trend was also reflected in the time series of mean annual runoff over the study region, where the runoff was notably lower (0.07 m a⁻¹) during 2000–2006 (Fig. S8), a period when the near-surface thick ice layers had not yet formed. The presence of thicker near-surface ice layers subsequently contributed to an increase in mean annual runoff (0.14 m a⁻¹) and the corresponding interannual variability during 2007–2022 (Fig. S8).

Ice layer fraction (IF), defined as the percentage of ice layer content over the length of a firn core, was reported from six shallow firn cores (up to 3 m depth) in the study area by Bezeau and others (Reference Bezeau, Sharp, Burgess and Gascon2013) based on 2012 measurements. These sites were revisited by Hallé and others (Reference Hallé2025) in 2022 (Table 2), providing an opportunity to assess observed decadal changes in IF. We compared these observed changes with simulated decadal variations in IF from the corresponding nearest-neighbour model grids (Fig. 6, Table 2). Although the modelled IF does not exactly match the observed IF in either 2012 or 2022, the downward advection of ice layers between 2012 and 2022 is evident in both the observations and our model simulations (Fig. 6). However, the spatial pattern of modelled IF during 2012 (Fig. 6) shows a similar decreasing pattern with elevation compared to the averaged IF valued during 1989–2003 (Fig. 8 in Colgan and Sharp, Reference Colgan and Sharp2008). The observed data (Hallé and others, Reference Hallé2025) indicate that the largest changes in IF occurred below 1575 m a.s.l., a pattern that is also reproduced in our model results (Fig. 6). This suggests that the model qualitatively captures the observed decadal changes in IF. The lower IF observed in both the observations and model simulations for 2022 was attributed to cooler summer conditions over the Devon Ice Cap in that year compared to 2012 (Hallé and others, Reference Hallé2025). This suggests that a reduction in IF within the firn layer can occur during cooler summers, highlighting the sensitivity of IF to interannual variability in surface air temperature. Importantly, changes in IF on multi-year timescales are not necessarily irreversible. A decrease in IF associated with a cooler summer may subsequently be offset by the formation of new ice layers during the next warmer summer (Thompson-Munson and others, Reference Thompson-Munson, Kay and Markle2024). This phenomenon is also evident in the simulated density stratigraphy shown in Fig. 5.

Figure 6. Comparison of observed decadal changes in IF at six shallow firn core sites with corresponding model simulations.

Table 2. Details of the six shallow firn core locations used to estimate decadal observed (Bezeau and others, Reference Bezeau, Sharp, Burgess and Gascon2013; Hallé and others, Reference Hallé2025) and modelled changes in IF.

4.4. Spatio-temporal variability of near-surface mean bulk firn density

We calculated the mean bulk firn density of the top 3 m across all model grids in the study region to examine its effect on runoff generation from each corresponding grid (Fig. 7). In 2007, we observed minimal spatial variability in bulk firn density across the study area (Fig. 7a). Our results indicated that bulk firn density was substantially high in 2012, particularly in the western part of the study area (Fig. 7b), due to development of ice layers within the top 3 m. In this region, the development of near-surface ice layers extended into the interior. Additionally, high bulk firn density was evident in the lower accumulation zone of the southeastern part of the study area in 2012. However, by 2017, these ice layers had advected to depths below 3 m (Fig. 7c), leading to a reduction in the spatial variability of bulk firn densities compared to 2012, particularly in the interior regions of the study area. In 2022, new ice layers formed within the top 3 meters across most parts of the study area, resulting in relatively high bulk firn density values (Fig. 7d).

Figure 7. Modelled top 3 m mean bulk firn density across the study area during the spring seasons of 2007 (a), 2012 (b), 2017 (c), and 2022 (d), illustrating changes over time. The colour scale indicates the corresponding density (kg m⁻³) values. The black dotted line represents the position of the CryoSat line (Figure 1b), which is referenced in the discussion of the vertical density profiles (Figure 5) in Section 4.3. The grey dotted lines represent the elevation contours.

Throughout the simulation period, we observed that mean bulk firn density values were consistently higher in the western part of the study area, where annual runoff was also somewhat greater (Fig. 10b), in comparison to the eastern part. Additionally, we found a positive correlation (significant at the p < 0.05 level) between annual runoff and the corresponding mean near-surface density in the lower accumulation zone of the western part of the study area during the simulation period (Fig. S9). This suggests that annual runoff was partly intensified by the presence of near-surface ice layers, which impacted the meltwater retention process (Section 4.3).

4.5. Meltwater retention

Meltwater retention, defined as the fraction of meltwater contributing to internal accumulation within the annual snowpack (Ashmore and others, Reference Ashmore, Mair and Burgess2020), exhibited mean values ranging from 0.34 to 0.81 across the study area between 1999 and 2022 (Fig. 8a). Retention was primarily controlled by mean summer temperatures (r = 0.85, p < 0.01), with cooler summers promoting higher retention and warmer summers enhancing surface melt and near-surface ice layer formation, which reduced the snowpack’s retention capacity (Fig. 8a).

Figure 8. The spatial variability of the mean fraction of annual meltwater retention (a) and the corresponding interannual variability (b) across the study area from 1999 to 2022. The grey dotted lines represent the elevation contours.

We used the coefficient of variation (C_v), defined as the ratio of the standard deviation to the mean, to assess interannual variability in meltwater retention across the study region (Fig. 8b). From 1999 to 2022, variability was low in the upper accumulation areas (>1600 m a.s.l.), where mean retention was higher (Fig. 8a), but relatively high in the upper sections of marine-terminating outlet glaciers over the Devon Ice Cap, driven mainly by annual melt fluctuations. Temperature played a key role in regulating runoff versus retention in these zones (Fig. 10f). While total melt and refreezing showed a slight increasing trend in the upper accumulation areas (Fig. S10), no significant trend (at the p < 0.05 level) in meltwater retention was detected across the study area.

4.6. Simulated SMB and runoff

The simulated spatial mean annual SMB for the period 1999–2022 ranged from −0.09 to 0.26 m w.e. a⁻¹ across the study area, with a mean value of 0.13 m w.e. a⁻¹ (Fig. 9a). These simulated values are in reasonable agreement with the observed spatial mean SMB along the NW transect (Fig. 1b), which varied between −0.02 and 0.27 m w.e. a⁻¹ during 2004–2022. It is noteworthy that all the runoff was generated during the summer months, implying that the magnitude of annual runoff exactly matches that of the summer runoff in this study. The simulated spatial mean annual runoff for 1999–2022 ranged from 0.03 to 0.31 m w.e. a⁻¹, with an average of 0.11 m w.e. a⁻¹ (Fig. 9b). In comparison, Ashmore and others (Reference Ashmore, Mair and Burgess2020) reported spatially averaged annual SMB and runoff of 0.21 m w.e. a⁻¹ and 0.05 m w.e. a⁻¹, respectively, over the same study area for the period 2001–2016 (Fig. 9). This indicates that, when averaged over the entire study area, our model produced a 61% lower positive SMB and a 54% higher runoff compared to Ashmore and others (Reference Ashmore, Mair and Burgess2020). We also extracted SMB and runoff data from the regional climate model MARv3.11.5 (Maure and others, Reference Maure, Kittel, Lambin, Delhasse and Fettweis2023) for our study area (Fig. 9). The spatial mean annual SMB ranged from −0.01 to 0.32 m w.e. a⁻¹, and runoff ranged from 0.01 to 0.29 m w.e. a⁻¹, with mean values of 0.15 m w.e. a⁻¹ and 0.14 m w.e. a⁻¹, respectively, over the period 2000–2020 (Fig. 9). Consequently, our modelled SMB and runoff are in close agreement with the MARv3.11.5 outputs (Maure and others, Reference Maure, Kittel, Lambin, Delhasse and Fettweis2023), while also more effectively simulating the depth-density evolution (Section 4.3) compared to MAR (Machguth and others, Reference Machguth2024). In addition, observed mass balance data from point scale shallow ice core measurements at eight locations across the Devon Ice Cap, reported by Mair and others (Reference Mair, Burgess and Sharp2005), ranged from 0.13 to 0.27 m w.e. a⁻¹ during 1963–2000. We were unaware of any field measurements of point-scale runoff in the selected study region to compare with the model estimates.

Figure 9. Modelled spatial mean annual time series of (a) SMB and (b) runoff across the study area, shown in cyan colour. The modelled SMB and runoff are compared with MARv3.11.5 data (Maure and others, Reference Maure, Kittel, Lambin, Delhasse and Fettweis2023; Ashmore and others, Reference Ashmore, Mair and Burgess2020), shown in golden yellow and salmon colour, respectively.

We observed pronounced spatial variability in mean annual SMB between the eastern and western regions of the study area for the period 2000–2022 (Fig. 10a). Overall, our model results show a weak positive correlation between SMB and elevation (r = 0.31, p < 0.01), though this relationship was significant (at p < 0.01 level) only in the northwest (r = 0.62) and southwest (r = 0.54) sectors (Fig. S11). No significant linear relationship was found in the northeast and southeast sectors, likely due to the influence of open ocean and enhanced cyclonic activity over Baffin Bay, which affects SMB in the eastern region (Koerner, Reference Koerner1966; Mair and others, Reference Mair, Burgess and Sharp2005). The non-linearity in the southeast is also supported by point-scale ice core mass balance measurements (Mair and others, Reference Mair, Burgess and Sharp2005).

Figure 10. Mean annual SMB (a) and summer runoff (b) across the study area for the period 1999–2022. Distribution of climate sensitivities of SMB to annual precipitation (c) and summer temperature (e), as well as the climate sensitivities of summer runoff to annual precipitation (d) and summer temperature (f). The unit of precipitation sensitivity was m a⁻¹, representing changes in SMB (c) or summer runoff (d) due to a 10% change in precipitation. The grey dotted lines represent the elevation contours.

4.7. Climate response of annual SMB and summer runoff

The model estimated sensitivities of annual SMB to precipitation (s^P_SMB) and temperature (s^T_SMB) (Eq. (5)) ranged from 0.010 to 0.032 m a⁻¹ for a 10% precipitation change and −0.011 to −0.101 m a⁻¹ °C⁻¹, respectively, and were significant at the p < 0.05 level (Fig. 10c and 10e). The corresponding mean values of s^P_SMB and s^T_SMB were 0.017 m a⁻¹ for a 10% precipitation change and −0.042 m a⁻¹ °C⁻¹, respectively, across the study region. Our results suggest that the annual SMB is more sensitive to changes in summer temperature than to changes in annual precipitation. In literature, precipitation sensitivity of SMB values ranged from 0.03 to 0.36 m a⁻¹ for a 10% change in precipitation, based on data from 42 different locations across the Arctic region using a temperature-index method (Woul and Hock, Reference Hock2005). Also, the corresponding temperature sensitivity values ranging from −0.2 to −2.01 m a⁻¹ °C⁻¹. When compared to other Arctic regions, the Devon Ice Cap’s precipitation and temperature sensitivities were relatively small. This is due to the Devon Ice Cap’s low annual precipitation and limited interannual variability, as well as the study area being primarily focused on the accumulation zones, which may contribute to lower temperature sensitivity.

In the interior of the Devon Ice Cap (elevations > 1600 m a.s.l.), SMB showed low sensitivity to temperature and precipitation (Fig. 10c and 10e). Sensitivities were substantially higher near the equilibrium line altitude (ELA), where SMB approaches zero (Fig. 10a), and mean annual runoff was elevated due to the absence of firn in the superimposed ice zone, causing most meltwater to become runoff (Fig. S12). This highlights the vulnerability of lower accumulation areas to warming. In contrast, sensitivities decline with increasing mean SMB due to firn retention of surface melt (Fig. S12). These trends align with previous studies (Braithwaite and Zhang, Reference Braithwaite and Zhang1999; Woul and Hock, Reference Hock2005). At SMB values > 0.2 m w.e. a⁻¹, precipitation sensitivities plateau and temperature sensitivities become negligible, reflecting the role of firn and refreezing at higher elevations.

The estimated sensitivities of summer runoff to precipitation (s^P_Q) and temperature (s^T_Q) (Eq. (6)) ranged from −0.022 to −0.000 m a⁻¹ for 10% precipitation change and 0.011 to 0.101 m a⁻¹ °C⁻¹, respectively, and were significant at the p < 0.05 level (Fig. 10d and 10f). The corresponding mean values of s^P_Q and s^T_Q were −0.008 m a⁻¹ and 0.042 m a⁻¹ °C⁻¹, respectively, across the study region. We are not aware of any previous estimates of climate sensitivity for runoff in the literature to compare with our findings in the study region.