Modeling the surface mass balance of Penny Ice Cap, Baffin Island, 1959–2099

Abstract

Glaciers of Baffin Island and nearby islands of Arctic Canada have experienced rapid mass losses over recent decades. However, projections of loss rates into the 21st century have so far been limited by the availability of model calibration and validation data. In this study, we model the surface mass balance of the largest ice cap on Baffin Island, Penny Ice Cap, since 1959, using an enhanced temperature index model calibrated with in situ data from 2006–2014. Subsequently, we project changes to 2099 based on the RCP4.5 climate scenario. Since the mid-1990s, the surface mass balance over Penny Ice Cap has become increasingly negative, particularly after 2005. Using volume–area scaling to account for glacier retreat, peak net mass loss is projected to occur between ~2040 and 2080, and the ice cap is expected to lose 22% (377.4 Gt or 60 m w.e.) of its 2014 ice mass by 2099, contributing 1.0 mm to sea level rise. Our 2015–2099 projections are approximately nine times more sensitive to changes in temperature than precipitation, with an absolute cumulative difference of 566 Gt (90 m w.e.) between +2 and −2°C scenarios, and 63 Gt (10 m w.e.) between +20% and −20% precipitation scenarios.

1. Introduction

Glaciers in the Canadian Arctic Archipelago (CAA) have experienced increasing mass loss rates in recent decades, particularly since 2005 (Gardner and others, 2011, 2012; Harig and Simons, 2016; Millan and others, 2017; Serreze and others, 2017; Thomson and others, 2017; Noël and others, 2018; Ciracì and others, 2020; Hugonnet and others, 2021). The area–averaged mass loss in the southern part of the archipelago (hereafter: southern CAA; Baffin and Bylot islands) was more than double that of the more northern sectors over the period 2003–2009 (Gardner and others, 2013). Noël and others (2018) used the Regional Atmospheric Climate Model 2.3 (RACMO2.3; https://www.projects.science.uu.nl/iceclimate/models/racmo-model.php) to show that mass loss rates in the southern CAA increased from 11.8 ± 4.5 Gt a^–1 over the period 1958–1995 to 21.9 ± 4.5 Gt a^–1 over the period 1996–2015. More recently, Ciracì and others (2020) used GRACE satellite gravity measurements to determine a mean climatic mass balance of −31.8 ± 5 Gt a^–1 for this region over the period 2002–2019. While temperatures have clearly increased in the southern CAA since the 1950s, there has been no significant change in precipitation (Gardner and others, 2012; Vincent and others, 2015; Noël and others, 2018). How quickly the rapid warming will lead to peak mass loss and subsequent disappearance of ice caps across this region is therefore an important question.

The largest ice cap in the southern CAA is Penny Ice Cap (67°N, 66°W; Fig. 1). Lenaerts and others (2013) used a coupled atmosphere/snow model forced with the IPCC's moderate RCP4.5 greenhouse gas concentration scenario to model sustained and irreversible glacier mass losses across the whole of the CAA, increasing from 29 ± 6 Gt a^–1 over 2000–2011 to 62 ± 10 Gt a^–1 by the end of the 21st century. Prior to the Lenaerts and others (2013) study, Penny Ice Cap was only modeled as part of global–scale assessments of glacier losses, which can have major uncertainties due to the lack of local model calibration or validation (e.g. Radic and Hock, 2011; Marzeion and others, 2012; Slangen and others, 2012; Huss and Hock, 2015). None of these studies were supported by spatially-distributed data from Penny Ice Cap itself, and their projections for the future of this ice cap were therefore poorly-constrained.

Figure 1. Map of Penny Ice Cap showing the five surface mass balance survey lines, automatic weather stations (Summit AWS, AWS2) and NASA Airborne Topographic Mapper altimetry (ATM altimetry) lines. Background image: Landsat 5, 19 August 1985.

In this study, we model the daily surface mass balance of Penny Ice Cap from 1958 to 2099 using an enhanced temperature-index model calibrated with in situ data collected on the ice cap between 2006 and 2014. Model outputs are used to quantify contributions to sea level rise, and to determine when the ice cap will reach its peak melt output. This is locally relevant since one of the most important hydrological outlets of Penny Ice Cap is in Akshayuk Pass (Fig. 1), which is a popular hiking route in Auyuittuq National Park. This valley has experienced sudden and dramatic increases in streamflow due to increased glacier melt and rainfall, resulting in the evacuation of hikers (The Innovator, 2008). Planning for visitor safety in the valley would therefore benefit from a better understanding of glacier meltwater contributions to streamflow and how they may change in the future. Furthermore, the modeled mass balance over Penny Ice Cap will complement other regional and global studies which report increased mass loss over the last few decades (e.g. Gardner and others, 2012, 2013; Lenaerts and others, 2013; Noël and others, 2018; Ciracì and others, 2020; Hugonnet and others, 2021; Otosaka and others, 2023). This is the first study to model the surface mass balance of Penny Ice Cap using extensive spatially-distributed data obtained from the ice cap itself and provides insight into how it will respond to future climate scenarios.

2. Study site

Penny Ice Cap (67°N, 66°W) covers an area of ~6300 km², with a summit elevation of ~1930 m a.s.l. (above sea level; Fig. 1), and maximum ice thickness of ~880 m (Shi and others, 2010). The ice cap terminates in a broad, gently sloping lobe-like region to the west, while major outlet glaciers flow from the ice cap interior down deeply-incised valleys toward the north, east and south. Two of these are tidewater terminating: Coronation Glacier in the south-east sector, and an unnamed glacier in the north-central sector (Fig. 1). The majority of the ice cap moves slowly (<20 m a^–1), particularly in the interior, with faster motion limited to the upper reaches of outlet glaciers, where velocities range from ~100–250 m a^–1 (Van Wychen and others, 2015; Schaffer and others, 2017). Over the period 1985–2011 there has been a general deceleration of outlet glaciers, at a rate of 12–25% decade^–1 (Heid and Kääb, 2012; Schaffer and others, 2017).

Based on annual surface mass balance measurements using stakes between 2006 and 2014 along four survey lines (Fig. 1), the mean equilibrium line altitude was ~1646 m a.s.l., varying between ~1320 and ~1820 m in low and high mass loss years, respectively. The mean surface balance rate in the same period was −1.2 m w.e. a^–1 between 329 and 1817 m a.s.l. Thinning rates are 3–4 m a^–1 at the ice cap margin, among the highest measured in the CAA, and recent surface melt rates are comparable to those last experienced over 3000 years ago (Fisher and others, 2011; Zdanowicz and others, 2012). Infiltration of surface meltwater has resulted in increased firn density since the mid-1990s and caused 10 m firn temperatures to rise ~10°C between the mid-1990s and 2011 (Zdanowicz and others, 2012). Measurements by an automatic weather station at the ice cap summit (Summit AWS, 67.25°N, 65.85°W; Fig. 1) in 2007 and 2008 recorded mean, maximum and minimum annual air temperatures of −15.4, 2.9 and −42.9°C, respectively (Zdanowicz and others, 2012). Based on microwave satellite data acquired from 2007–2010, summer melt typically begins in late May and ends in early September (F. Dupont, personal communication, 2015).

3. Methods

In this study we follow the terminology described by Cogley and others (2011), whereby surface mass balance is defined as the sum of surface accumulation and surface ablation (i.e. refreezing in the snow/firn is not included). The balance, including the internal mass gain from refreezing, is referred to as climatic mass balance. Upper-case symbols are used to denote glacier-wide terms, while lower-case symbols refer to point values.

The glacier-wide surface mass balance (B) over Penny Ice Cap was calculated from the sum of ablation (A), accumulation (C) and frontal ablation (A _f, mass losses at the marine-terminating front, e.g., by iceberg calving and submarine melt). Mass gains are defined positively and losses negatively:

(1)$$B = A + C + A_f$$

3.1 Surface mass balance

The surface mass balance was calculated with the open-access Distributed Enhanced Temperature Index Model (DETIM; Hock, 1999; https://regine.github.io/meltmodel/). We have chosen to use a temperature index model because there is insufficient data to calibrate a full energy balance model. Temperature index models have been shown to produce realistic simulations of glacier melt and meltwater output, especially at the catchment scale and on seasonal to interannual time scales (Hock, 1999, 2003). For example, Huss and Hock (2015) found that to produce global estimates of glacier mass loss, a classical degree-day model performed better than a simplified energy balance model. Furthermore, glacier surface mass balance in the CAA is highly correlated to summer temperatures (Gardner and others, 2011; Sharp and others, 2011; Noël and others, 2018).

For each gridcell, surface ablation by melt (a) is calculated by:

(2)$$a = \left\{{\;\matrix{ {( {\;F_m + F_{r\;snow/ice}I} ) T, \;} \hfill & {T > 0} \hfill \cr {0, \;} \hfill & {T \le \;0} \hfill \cr } } \right\}$$

where F _m is a melt factor (mm d^–1°C^–1), F _{r snow/ice} is a radiation factor for snow or ice (mm m² W^–1 °C^–1 d^–1), I is the potential clear-sky direct solar radiation (W m^–2) and T is the daily mean air temperature (°C). Specific values of F _m and F _r were derived for Penny Ice Cap, as described in section 5.0.

For each day of the year, I was calculated from solar geometry and topographic shading (Hock, 1999) determined from the Canadian Digital Elevation Dataset, which has an average gridcell size of 37 × 93 m over the study area. Calculating the solar radiation is computationally intensive, so we used the same daily values for each year of the model run.

Snow accumulation at each gridcell (c) was computed from precipitation (p) and a threshold temperature (T _snow) to distinguish between solid precipitation (snow accumulation) and rainfall, with linear interpolation of the snow fraction within the range $T_{snow}$ ± 1°C:

(3)$$c = \delta \;p\;\left\{{\matrix{ {\delta = 1, \;} \hfill & {T\;\le \;T_{snow}-1} \hfill \cr {\delta = \displaystyle{{T_{snow}} \over 2}, \;} \hfill & {\;T_{snow}-1 < \;\;T < \;T_{snow} + 1} \hfill \cr {\delta = 0, \;} \hfill & {T\;\ge \;T_{snow} + 1} \hfill \cr } } \right\}$$

3.2 Frontal ablation

The A _f term accounts for all mass losses at a marine-terminating glacier front, including calving and submarine melting. A constant value of A _f of −0.02 Gt a^–1 (Van Wychen and others, 2015) was used rather than incorporating a separate model (e.g. Trussel and others, 2015). This is because frontal ablation is a minor component of the mass budget of Penny Ice Cap, with only two tidewater glaciers accounting for ~0.2% of its total net mass loss in 2011 (Van Wychen and others, 2015). It is not known how this may have varied over time.

3.3 Geometry changes

Glacier mass loss is often accompanied by a reduction in glacier extent, resulting in mass-balance feedbacks. Ice is typically lost at the lowest elevations first, so the remaining glacier mass has a higher mean elevation and smaller area, which then reduces the overall specific mass loss rate. We accounted for changes in glacier extent based on a volume–area scaling approach (Bahr and others, 2015):

(4)$$V = d( S ) ^\gamma $$

where d and γ are empirical coefficients, V is the glacier volume, and S is the glacier area. Here we use its differentiated form to relate the change in glacier volume ΔV to the change in surface area ΔS following Arendt and others (2006):

(5)$$\Delta V = d\gamma \Delta S( {S_i} ) ^{\gamma -1}$$

where S _i is the initial glacier surface area (surface area at the beginning of the DETIM model run). At the end of each mass-balance year, ΔV was derived from the modeled mass balance assuming a mean bulk density of 900 kg m³, and this was then used to calculate ΔS. For the scaling parameter γ we used a value of 1.25, as suggested for ice caps (Bahr and others, 2015). Parameter d was derived by solving Eqn (4) using the total volume obtained from Penny Ice Cap ice thicknesses provided by Huss and Farinotti (2012), and the associated ice cap area, to yield a value of 0.9608 m^3–2γ. This value is considerably lower than the value of 1.7 m^{3–2 γ} suggested by Radić and Hock (2010), but only slightly lower than values suggested by Grinsted (2013).

We compared these thicknesses to those measured with the NASA airborne Multichannel Coherent Radar Depth Sounder in 2013 during Airborne Topographic Mapper (ATM) altimetry flights. These showed a similar distribution of ice thickness with elevation, and, on average, the NASA thicknesses were 5.8 m greater than modeled. We consider this difference to be acceptable given that the ice cap is hundreds of meters thick, reaching a maximum of ~880 m (Shi and others, 2010), and that ablation rates reach up to 4 m yr^–1 at low elevations (Schaffer and others, 2020).

The glacier area is automatically adjusted in DETIM by calculating the number of grid cells represented by ΔS, then removing grid cells sequentially starting with those at the lowest elevation. If the mass balance is positive, the area is kept constant. This method of calculating B assumes that all the volume loss is compensated for by retreat (rather than elevation change) and provides a theoretical upper limit for retreat.

4. Model inputs

4.1 Climate data

DETIM was forced by daily data of near-surface air temperature and precipitation from RACMO. For the period 1959–2014 we used the output from RACMO2.3 forced by reanalysis data from ERA-interim, with a grid resolution of 11 km (Lenaerts and others, 2012; Noël and others, 2015). For the period 2015–2099, we used outputs from an earlier version of the model, RACMO2.1 (Lenaerts and others, 2013; Van Angelen and others, 2013), because the RACMO2.3-generated climate fields were not available beyond 2015. The main difference between RACMO2.3 and RACMO2.1 is that RACMO2.3 has been improved with major changes in the description of cloud microphysics, surface and boundary layer turbulence, and radiation transport (Noël and others, 2015). Precipitation outputs have also been modified to be exclusively snowfall under freezing conditions, but this update does not impact our results since we provided DETIM with the RACMO's total precipitation and with a calibrated, ice cap-specific snow/rain threshold temperature. Inputs to RACMO2.1 were from the Coupled Model Intercomparison Project Phase 5 (CMIP5) general circulation model HadGEM2-ES (Lenaerts and others, 2013), which was itself forced using the RCP4.5 scenario leading to a stabilized radiative global mean forcing of ~4.5 W m^–2 by 2100. This scenario also leads to a global mean warming of 1.8°C and a 3.6% increase in precipitation over the period 2081–2100, relative to the 1986–2005 average.

4.1.1 Validation and bias correction of RACMO2.3 outputs

The RACMO2.3 temperature outputs were highly correlated with the near-surface temperatures recorded by the Summit AWS on Penny Ice Cap over the period 2007–2014 (Pearson's correlation coefficient r = 0.98, p < 0.01; Figs 2a and b). However, the modeled temperatures were slightly cooler than observed temperatures, with an average offset of −0.28°C. An offset of +0.28°C was therefore added to each RACMO2.3 data point and the resulting mean annual air temperature at the summit over the period 2007–2014 was ~–15.3°C. This agrees well with the estimated mean annual temperature at the summit of −16 ± 1.5°C based on AWS data collected between 1992–2000 and 2007–2011 (Zdanowicz and others, 2012).

Figure 2. (a, b) Comparison between in situ daily mean air temperatures measured at the Summit AWS and the modeled 2 m air temperature from the closest RACMO2.3 gridcell. (c, d) Comparison between in situ spring snowpack measurements at two mass balance stakes closest to the Summit AWS (P000 and P101) vs RACMO2.3 cumulative winter total precipitation (snow and rain). Each year's snow pack values in (c), for both the in situ and RACMO2.3 data, refer to the end of winter (~April) and represent the accumulated snowfall since the end of the previous summer (~Sept.). For example, year 2008 refers to snow pack measurements obtained in April 2008 that represent the accumulated snowfall since Sept. 2007. The black line in (b) and (d) is the 1:1 line.

The RACMO2.3 modeled precipitation outputs for the period 2007–2014 were compared to late winter measurements of the snowpack water equivalent (SWE in mm) on Penny Ice Cap summit (stakes P000 and P101; Figs 2c and d). Here, ‘winter’ refers to the span of time from the end of the melt season (typically early September) to the time when the snowpack measurements were taken, usually mid-April. Comparisons were made by assuming that all RACMO2.3 precipitation (whether rain or snow) remained within the winter snowpack, and by summing the thickness and density of layers above the last summer surface to determine cumulative snow water equivalent. The RACMO2.3 cumulative winter precipitation captured the interannual variability in the winter SWE (r = 0.65, p = 0.021; Figs 2c and d), but underestimated winter precipitation by 21% on average. The RACMO2.3 precipitation values were therefore increased by 21% for use in DETIM.

The average RACMO2.3 precipitation between 1963 and 2011 was also compared to the average net accumulation recorded in firn cores taken near the ice cap summit over the same period. The firn cores recorded an average net accumulation of 0.40 ± 0.05 m w.e. a^–1 (Zdanowicz and others, 2012), which is equal, within error limits, to the unadjusted RACMO2.3 value of 0.45 m a^–1 w.e., but less than the adjusted precipitation of 0.54 m a^–1 w.e.

4.1.2 Validation and bias correction of RACMO2.1 outputs

The RACMO2.1 temperature outputs were well correlated with temperatures recorded by the Summit AWS between 2007 and 2014 (r = 0.79, p < 0.01; Figs S1a and b), but less correlated than RACMO2.3 and cooler, with an average offset of −2.05°C. RACMO2.1 precipitation outputs did not correlate well with measurements of snowpack water equivalence, underestimating it by 31% on average over the period 2007–2014 (Figs S1c and d). We modeled the surface mass balance with the RACMO2.1 data in two trial runs: (i) using the offsets from in situ values (2.05°C and 31% increase in precipitation), and (ii) with offsets used for RACMO2.3 (0.28°C and 21% increase in precipitation). When using the RACMO2.3 offsets, the mean surface mass balance calculated over the 2005–2013 calibration period of −4.59 Gt a^–1 (−0.73 m w.e. a^–1 on average over the entire ice cap) was much closer to the mean surface mass balance measured with the 2005–2013 ATM-altimetry (−4.56 Gt a^–1 or −0.72 m w.e. a^–1), compared to using the RACMO2.1 offsets, and nearly identical to the modeled mass balance using the RACMO2.3 data for the same time period. We therefore decided to apply the RACMO2.3 offsets (0.28°C and 21% increase in precipitation) to the RACMO2.1 dataset.

4.1.3 Extrapolation of temperature and precipitation values across Penny Ice Cap

DETIM was forced with RACMO data for the gridcell closest to the ice cap Summit AWS. From this point daily mean air temperature data were extrapolated to other grid cells using monthly lapse rates (Table 1) between the summit and ~490 m a.s.l. on Glacier 1 (Table 1), where data from another AWS were available. In July 2012, the estimated lapse rate between the two AWS was 4.69°C km^–1, which is nearly equal to the lapse rate of 4.66°C km^–1 obtained from RACMO at these grid cells. These figures are also very close to earlier estimates of lapse rates on CAA ice caps of 4.6–4.9°C km^–1 (Mair and others, 2005; Shepherd and others, 2007; Gardner and others, 2009). In the present study, daily lapse rates were held constant within each particular month (Table 1).

Table 1. Average monthly lapse rates for Penny Ice Cap derived from RACMO2.3 data using temperatures from 2007–2014 at a gridcell near the ice cap summit and at ~490 m a.s.l

The negative sign indicates a decrease in temperature with increasing elevation. These lapse rates were used to extrapolate daily mean air temperature data to each gridcell over the ice cap.

Precipitation over Penny Ice Cap was extrapolated using a gradient calculated from total annual precipitation near the ice cap summit and close to sea level. The total precipitation near the summit (~1817 m a.s.l.) was estimated from the mean net accumulation in firn cores between 1963–2011 (0.40 m a^–1 w.e.), using the assumption that at this high altitude net accumulation is approximately equal to total annual precipitation (discounting losses from sublimation and winter wind scouring). At lower elevations precipitation data from the Pangnirtung weather station (23 m a.s.l.; Fig. 1) was used. Here the mean annual precipitation was ~0.24 m a^–1 w.e. between 2008 and 2014, yielding a precipitation gradient of + 3.8% per 100 m increase in elevation over the ice cap. While somewhat crude, the estimate is necessarily constrained by the scarcity of data available for verification. Above 1900 m a.s.l., the precipitation was left constant in DETIM to account for reduced air moisture content and increased wind scouring at higher elevations.

4.2 Initial snow cover

DETIM requires an estimate of the ice-cap-wide snow extent and thickness at the start of a simulation. We used the average winter SWE measured at each stake on Penny Ice Cap from 2006–2014 to create a linear model of changes in snow cover with respect to elevation, after applying a square-root transformation to the snow cover data to obtain normally distributed residuals (coefficient of determination r ² = 0.16, F-test p ≤ 0.001, standard error of coefficient p ≤ 0.001). This model was then used to create a snow cover map over the entire ice cap using the Canadian Digital Elevation Dataset, which was in turn used for initializing the calibration process. Subsequently, the snow cover was set to zero automatically by DETIM on the first day of the winter mass balance season each year (September).

The modeled snowpack SWE has a gradient of 10.2% per 100 m, which is much larger than the total precipitation gradient of 3.8% per 100 m because the former only accounts for the fraction of the total precipitation which accumulates as snow on the ice cap between ~September and April. These snow fractions were 48% of the total precipitation at Pangnirtung and 80% at the Summit AWS, respectively, based on our snowpack model. We approximated the fraction at the Pangnirtung AWS by summing the total precipitation during winter (defined as mid-September to May) between 2008 and 2014, which yielded a fraction of 48% matching our modeled results here.

4.3 Glacier inventory and surface type data

For the historical (1959–2014) simulation, changes in glacier area over time were incorporated by updating an initial 1959 glacier outline with 1975, 2001 and 2014 versions. The initial 1959 outline was created by Evelyn Dowdeswell (University of Bristol) from aerial photographs taken in 1959 by the Royal Canadian Air Force, with minor modifications including orthorectification and the addition of peripheral ice masses. The outlines from 1975 (16 August & 2 September), 2001 (29–31 July) and 2014 (26 July) were derived from cloud-free Landsat images. Likewise, the surface type (snow, ice or firn), needed to choose the melt factor, was updated with end of summer surface type grids for 1975, 2001 and 2014, manually derived from these Landsat scenes. A Landsat image surface reflectance threshold of 180 (band 4 for 1975, band 8 for 2001 and 2014) was used to delineate snow patches, with all patches >0.5 km² mapped.

5. Model calibration

DETIM was customized for Penny Ice Cap by calibrating four parameters in the model (F _m, F _{r snow}, F _{r ice} and T _snow) using in situ point surface mass balance measurements and NASA ATM elevation change observations (Table 2). Point mass balances were measured along the four survey lines (Fig. 1) at 42 locations and over an elevation range of 71 to 1822 m a.s.l. A total of 120 measurements made between 2006 and 2014 were used to calibrate the DETIM model. The ATM altimetry measurements were carried out in spring 1995, 2000, 2005, 2013 and 2014, prior to the ablation season (Fig. 1; see Schaffer and others (2020) for further details).

Table 2. Parameter values for the optimal, maximum and minimum parameter combinations for Penny Ice Cap that met the calibration requirements

Modeled outputs were compared to in situ surface mass balance measurements to calculate the RMSE and r ². Maximum and minimum correspond to the parameter combinations with the largest and smallest ice cap-wide mass loss rate selected from a range of 13 parameter combinations (Fig. S2), respectively. The ice cap-wide mass loss rate modified for firn densification over the same time period derived from ATM altimetry data was −4.6 Gt a⁻¹ ± ~1.9 Gt a⁻¹, using a density of 900 kg m⁻³.

The model calibration was accomplished using a Monte Carlo approach and the optimal values for each parameter were found by systematically varying all parameters over a range of physically plausible values. We required that the cumulative ice-cap-wide surface mass balance simulated by DETIM over the period 2005–2013 match the ATM altimetry-derived geodetic mass balance of −4.56 Gt a^–1 along flight lines (Schaffer and others, 2020). To ensure a physically realistic model we also required that F _m be >0 mm d^–1°C^–1, that F _{r ice} be 0.2 mm m² W^–1°C^–1 d^–1 greater than F _{r snow}, and that T _snow be ≤3°C. The optimization was further constrained by the requirement that the RMSE between the modeled and in situ 2006–2014 mass balance be ≤ ± 0.25 m w.e. a^–1, which is the estimated typical error for in situ surface mass balance measurements (Cogley and others, 1996). Daily DETIM mass balance outputs were summed over each mass balance year for comparison with annual in situ measurements. The mass balance year for this purpose was defined as starting on the date that the in situ mass balance measurements were taken, which was normally in April, but varied from year to year. This calibration method forces a good fit of the model output to the decadal geodetic balance to ensure that the long-term ice cap response to climate is realistically captured.

With F _m, F _{r snow}, F _{r ice} and T _snow initial increments of 0.5 mm d^–1°C^–1, 0.2 mm m² W^–1°C^–1 d^–1, 0.2 mm m² W^–1°C^–1 d^–1 and 1°C, respectively, a total of 13 parameter combinations were identified that met the aforementioned criteria (Table 2, Fig. S2). Additional model runs were then performed to identify the combination with the lowest RMSE. This step used increments of 0.02, 0.02, 0.02 and 0.5 for F _m, F _{r snow}, F _{r ice} and T _snow, respectively, and iterative convergence was achieved when there was no further change in RMSE to three decimal places. The optimal parameter combination had an RMSE of 0.45 m w.e. and r ² of 0.77 (Fig. 3; Table 2). These 13 best-performing parameters sets were used for further analysis.

Figure 3. Observed annual point mass balances and corresponding modeled values for: (a) the optimal parameter combination; (b) the optimal, minimum and maximum parameter combinations plotted against elevation for the period 2005–2013. Of the 13 parameter combinations with an RMSE ≤0.5 m w.e. the combination which results in the most negative mass balance is referred to as ‘minimum’, the most positive ‘maximum’ and closest to in situ and NASA altimetry data as ‘optimal’. The RMSE was calculated from the difference between observed and modeled values.

6. Model validation

The simulated surface mass balance was compared against altimetry-derived geodetic balances from 1995–2000 (Abdalati and others, 2004) and 2000–2005 (Gardner and others, 2012) not previously used in the calibration process (Fig. 4). For this comparison DETIM outputs were converted from units of m w.e. to Gt. The original ATM data were converted to mass change for this plot using a density of 900 kg m⁻³ (solid lines) and 850 kg m⁻³ (dashed line; see Schaffer and others (2020) and Huss (2013) for further details). Modeled results closely replicated the ATM-derived mass changes between 2000 and 2005 but overestimated the mass loss rate between 1995 and 2000 (Fig. 4). Over the entire period 1995–2005, the range of model results lie entirely within the error limits of the ATM altimetry data.

Figure 4. Modeled and measured mass-change rates over Penny Ice Cap averaged over three multi-year periods during 1995–2013. Shaded portions represent the 95% confidence intervals for the ATM altimetry data (blue) and the range of mass loss estimates obtained from the 13 DETIM parameter combinations (red). The bold red lines show the mass change obtained with the optimal DETIM parameter combination. The dashed blue line represents the mass change inferred from the 2005–2013 ATM altimetry data adjusted to account for the change in elevation due to firn densification (Schaffer and others, 2020), which was used for the model calibration.

The mean surface mass balance measured with ATM altimetry over the 2005–2013 calibration period was −4.56 Gt a^–1 (−0.72 m w.e. a^–1; dashed blue line in Fig. 4) which is very similar to that obtained using DETIM (−4.59 Gt a^–1 or −0.73 m w.e. a^–1; solid red line in Fig. 4). The modeled surface mass balance using RACMO2.1 and RACMO2.3 input data for the same period is nearly identical, providing confidence that RACMO2.1 results for 2015–2099 will be similar to those that would be obtained with RACMO2.3 input data if the latter were available after 2015.

7. Model projections

We first ran DETIM to simulate the historical surface mass balance variations for Penny Ice Cap over the period 1959–2014, without volume–area scaling. Subsequently, the mass balance time series of the most representative outlet glacier, Glacier 1 (Fig. 1), was compared to that of the entire ice cap. The correlation between the two was strong (r = 1.00, p = <0.001), implying that the mass balance trends for Glacier 1 can be considered representative of the ice cap as a whole. We therefore restricted our future projections (2015–2099) to Glacier 1, and then upscaled these. The simulated net mass losses for Glacier 1 over the period 1959–2014 represent 10.4% of the ice-cap-wide losses, so we used this ratio to upscale the projected changes for Glacier 1 to the entire ice cap.

For the future projections, DETIM was first run with a constant glacier area from 2014, and again with an evolving glacier area over time using the volume–area scaling described in section 3.3. The 2014 surface type grid (i.e., spatial distribution of firn vs glacier ice) was used for the future projections. Glacier mass losses were converted to sea level rise equivalent by assuming an ice density of 900 kg m^–3 and a global ocean area of 3.625 × 10⁸ km² (Cogley and others, 2011).

8. Results

Simulations of past and future surface mass balance changes over Penny Ice Cap from 1959 to 2099 are shown in Figure 5. For the historical period 1959–2014, the simulations show a relatively stable, albeit negative, mass balance until the mid-1990s, followed by an increasingly negative mass balance since then. The spatial pattern of simulated mass balance changes over the entire ice cap is shown in Figure 6. As can be seen, the western sector and low-lying outlet glaciers show the largest losses, >2.5 m w.e. a^–1 by 2010–2014. There is also an obvious decrease in the size of the accumulation area through time (area outlined with black line in Fig. 6).

Figure 5. Annual modeled surface mass balance rate for Penny Ice Cap between 1959 and 2099. The solid line is the optimal model parameter combination, while the shaded portion shows the range covered by the 13 parameter combinations. A constant glacier area was assumed for the dark blue series, while volume-area scaling was applied to the green series. Calving is not shown, but would add an additional 0.02 Gt a ¹ of mass loss. The thicker lines are 10-year running means.

Figure 6. Decadal averages of the specific surface mass balance over Penny Ice Cap modeled for the period 1960 to 2014. For the last interval (2010–14), values shown are the 5-year means. All areas outlined with a black line are specific surface mass balance values above zero.

For future projections (Fig. 5), the simulation that includes volume–area scaling (green line) produces a lesser rate of mass balance decrease than without it (blue line), which results in important differences at the end of the simulation period. The mass balance of Penny Ice Cap is expected to decrease over most of the period, reaching a minimum in ~2070. The lowest 10-year mean values are −9.42 Gt a^–1 ( − 1.49 m w.e. a^–1) and −5.31 Gt a^–1 ( − 0.84 m w.e. a^–1), for the two projections, which are 2.1 and 1.2 times more negative, respectively, than the average surface mass balance of −4.56 Gt a^–1 ( − 0.72 m w.e. a^–1) over the calibration period (2005–2014).

Figure 7 shows the projected pattern of mass balance changes over the Glacier 1 catchment from 2010 to 2099. The surface mass balance minimum in ~2070 coincides with when the accumulation area on Glacier 1 essentially disappears in the 2070s, with this loss sustained moving forward in time. The most negative specific balances occur at the glacier terminus, but the greatest total mass losses occur in higher elevation bands that cover a relatively large area. To determine the variability in total mass change with altitude, the specific mass balance was multiplied by the area for each 20 m elevation band (Fig. 8a). Results show that largest mass losses for Glacier 1 are expected to occur between approximately 600–1450 m a.s.l. (0.5 Gt a^–1 or 74% of the total mass loss; Fig. 8b).

Figure 7. Modeled surface mass balance for the period 2010 to 2099 for Glacier 1 with a constant glacier area. Ten-year averages are shown. All areas outlined with a black line are specific surface mass balance values above zero.

Figure 8. (a) Hypsometry (gray bars) and average predicted surface mass balance rate of Glacier 1 over the period 2010–2099 in m w.e. a⁻¹, and (b) in Gt a⁻¹. The elevation range within which the largest mass losses are expected to occur (~600–1450 m) is highlighted with gray shading in (b).

By 2099 Penny Ice Cap is projected to lose between 22 and 35% of its 2014 volume for the volume–area scaling and constant area approaches, respectively (Table 3). These figures translate to total mass losses of 377.4 and 602.2 Gt, respectively. Adding the mass changes due to frontal ablation slightly increases the projected mass loss rate by 0.02 Gt a^–1. The projected mass losses translate to 1.0 to 1.7 mm of sea level rise. The predicted cumulative mean surface lowering over Penny Ice Cap by 2099 is 59.8 to 95.0 m. If the mass loss rates of 4.0 Gt a^–1 (volume–area scaling) or 8.8 Gt a^–1 (constant glacier area) between 2090 and 2099 are sustained in the future, the ice cap is projected to disappear entirely by the early 2200s for the constant area scenario, and by the mid–2400s for the volume–area scaling scenario.

Table 3. Projected cumulative mass loss for Penny Ice Cap between 2014 and 2099 in Gt; as a percentage of total mass lost compared to its 2014 volume; in terms of sea level rise contributions; and mean ice cap-wide surface lowering for the constant area and volume-area scaling approaches

Adding the mass changes due to frontal ablation slightly increases the volume loss by −0.02 Gt a^−t. The expected disappearance date for both approaches assumes that the average rate of mass loss between 2090–2099 continues into the future.

9. Discussion

Our simulations show that the surface mass balance of Penny Ice Cap has become increasingly negative since the mid-1990s (Figs 5 and 6), in good agreement with the modeled surface mass balance of Noël and others (2018). A more negative mass balance since the mid-1990s is also seen from in situ observations and models from other ice caps in the northern CAA (Lenaerts and others, 2013). A pronounced increase in the modeled mass loss rate occurred after 2005, which also agrees with previous historical surface mass balance modeling (Noël and others, 2018), with ATM altimetry measurements over the ice cap (Fig. 4), and with patterns of mass change over the entire CAA derived from satellite elevation changes (ICESat) and gravity measurements (GRACE) (Gardner and others, 2012; Ciracì and others, 2020). Similar patterns of mass change have also been observed for the Greenland Ice Sheet (Otosaka and others, 2023) and globally (Hugonnet and others, 2021). Ten-year running means of the modeled surface mass balance for Penny Ice Cap are negative for the entire simulation period (1959–2099), supporting projections of irreversible mass losses of CAA glaciers (Lenaerts and others, 2013). Penny Ice Cap is expected to lose 22 and 35% of its 2014 volume by 2099 for the volume–area scaling and constant area approaches, respectively. The cumulative surface mass loss between 2015 and 2099 is reduced by 224.8 Gt (Fig. 9a) when volume–area scaling is applied. In this case, mass loss results in the removal of glacier area at the lowest elevations, resulting in a reduced overall specific mass loss rate through mass-balance feedbacks (see section 3.3), which explains the striking difference between the two projections. The volume–area scaling approach is comparable in magnitude to the bulk glacier volume loss of 18% for the whole southern CAA predicted by Lenaerts and others (2013) using the same RCP4.5 scenario and RACMO2 climate dataset.

Figure 9. Cumulative modeled surface mass balance for the period 2015 to 2099 for: (a) four temperature scenarios; (b) four temperature scenarios with volume area (VA) scaling; (c) four precipitation scenarios; (d) four precipitation scenarios with VA scaling. Outputs using the unmodified RACMO2.1 data are plotted as black solid and dashed lines for the constant area and volume area scaling scenarios, respectively.

Both approaches show an eventual stabilization of mass loss rates (Fig. 5). When volume–area scaling is accounted for this stabilization begins by ~2030, and by ~2070 when it is not. We hypothesize that the stabilization is due to the aforementioned mass-balance feedback, which is amplified when volume–area scaling is included. Radic and Hock (2011) predicted that mass losses for glaciers and ice caps worldwide will peak in the 2040s, while Marzeion and others (2012) predicted a globally-averaged peak in glacier mass loss rates between ~2050 and 2060 under RCP4.5 when incorporating volume–area scaling. The 10-year running mean peak mass loss occurs between 2043 (− 5.04 Gt a^–1) and 2078 (− 5.31 Gt a^–1), which overlaps well with these globally-averaged estimates.

9.1 Mass balance sensitivity

To assess the sensitivity of the projections for Penny Ice Cap to climatic variables, simulations were performed over the period 2015 to 2099 for four different temperature and precipitation scenarios. Uniform temperature changes (+1, + 2, −1 −2°C) were added to each daily temperature in the adjusted RACMO2.1 dataset, while precipitation was unaltered. Likewise, uniform precipitation changes ( + 10, + 20, −10, −20%) were applied to the RACMO2.1 precipitation dataset without any adjustments to air temperature. Given the scenarios evaluated (Fig. 9), the mass balance projections are far more sensitive to temperature changes than to precipitation changes. The absolute difference in the cumulative surface mass balance (2015–2099) between the most extreme temperature scenarios (+2 and −2°C) is 1067 Gt (Fig. 9a), but only 104 Gt between the most extreme precipitation scenarios (+20 and −20%; Fig. 9c).

When volume–area scaling is implemented in the model, the projected cumulative surface mass loss between 2015 and 2099 is reduced by 37% under the RCP4.5 scenario (from 592.8 to 373.1 Gt), and by 45% with the additional +2°C perturbation (from 1183.0 to 644.8 Gt). The absolute difference in the cumulative surface mass balance (2015–2099) between the most extreme temperature scenarios ( + 2 and −2°C) is 566 Gt (Fig. 9b), but only 63 Gt between the most extreme precipitation scenarios ( + 20 and −20%; Fig. 9d).

9.2 Refreezing

A limitation with our DETIM-based findings is that they only provide values for surface mass balance, and not climatic mass balance, as they do not explicitly account for refreezing. Given that refreezing of meltwater has become an increasingly important mass accumulation process on Penny Ice Cap (Zdanowicz and others, 2012) and other CAA ice caps (Bezeau and others, 2013; Gascon and others, 2013) over the past decade or two, it is useful to consider the effect that it may have on the ice cap mass balance.

Borehole and shallow core measurements from the summit of Penny Ice Cap show that the firn density has been increasing since the mid-1990s due to the formation of thick infiltration ice layers. In the future it is expected that the firn zone will be replaced completely with superimposed ice. To determine when this will occur we modeled the refreezing process separately using four refreezing parameterizations and variations outlined in Reijmer and others (2012; their Eqns 5–7) and Janssens and Huybrechts (2000; their Eqn 7). All the parameterizations presume that the refrozen mass (R) is equal to the minimum of either the available water mass (W _r), assumed equal to the mean annual snowfall, or the maximum amount of water that can potentially be refrozen (P _r), which depends on the energy available to melt all the snow:

(6)$$R = min[ {P_r, \;\;W_r} ] $$

If R is limited by the energy available, only part of the snowpack is melted and refrozen. This expression includes not only the water that refreezes in snow but also the water that refreezes at depth to form superimposed ice (Reijmer and others, 2012). The parameterizations differ in their calculation of P _r, with the simplest based on the firn cold content and more complex versions additionally accounting for the filling of pore spaces within the firn.

The parameterization outputs were compared to firn cores collected near the ice cap summit (Zdanowicz and others, 2012), and the model with the best fit to this in situ data was selected. These cores provide a proxy record of annual surface melt on Penny Ice Cap over the period 1963–2010, based on the volumetric percentage of infiltration ice (melt feature ‘MF’). The period starts in 1963 as this year is identifiable by a radioactive layer in the firn (Zdanowicz and others, 2012). The best fit to the firn core proxy melt record was obtained using a simple parameterization for P _r that only accounts for the cold content given by Huybrechts and De Wolde (1999):

(7)$$P_r = \displaystyle{{C_i} \over {L_f}}\;h_a\vert {T_s} \vert $$

where C _i is the heat capacity of ice (2050 J kg^–1 K^–1), L _f is the latent heat of fusion (0.334 × 10⁶ J kg^–1), h _a is the thickness of the thermally active layer and T _s is the annual mean glacier surface temperature which was obtained from RACMO2.3. The thickness of the thermally active layer is the maximum depth of the 0°C isotherm in summer. The value of W _r was assumed to equal the mean annual snowfall, also obtained from RACMO2.3 near the ice cap summit. Values of h _a were calculated from thermistor string measurements made in boreholes near the summit of Penny Ice Cap in 1953 and 2011 (Zdanowicz and others, 2012). In 1953, h _a was ~1.4 m and the firn temperature at ~10 m depth was ~–13°C. For 2011, h _a was ~3.1 m and the 10 m temperature was −3°C.

Equation 6 was calibrated to estimate values of the summer melt percentage from firn cores by varying h _a through time within a realistic range of values informed by the thermistor string measurements. The best fit was obtained for values which varied linearly between 1.3 m in 1963 and 3.1 m in 2014. A simulation of refreezing at the summit of Penny Ice Cap for 1963–2014 using the optimal parametrization (Fig. 10a) reproduces the observed long-term trend recorded in firn cores (difference <1%, and within 95% confidence limits). While the model does not accurately capture the interannual variations, this is less important for long-term projections.

Figure 10. Refreezing parameterization outputs for Penny Ice Cap, expressed as a percentage (R/ W _r) for: (a) 1963–2014 forced with RACMO2.3 data (black line) and (b) 2005–2098 forced with RACMO2.1 data (5 year running mean in blue). Outputs in (a) are compared to the volumetric percentage of refrozen meltwater ice (MF = melt feature) in firn cores collected near the summit from 1963–2014 (gray line). The linear regression lines for the modeled (blue line) and measured firn core data (yellow dashed line) are shown.

The refreezing model was applied until 2098 using the RACMO2.1 data. The 5 year running mean (Fig. 10b) shows an obvious trend of increasing MF % over time, approaching 95% by 2098, indicating that by that time firn at the surface has nearly disappeared, to be replaced by superimposed ice.

10. Conclusions

This study is the first to model the surface mass balance of Penny Ice Cap using extensive spatially-distributed data obtained from the ice cap itself. To calibrate DETIM we used: (a) a comprehensive dataset of in situ mass balance measurements covering a representative range of elevations, slopes and aspects, and (b) ATM altimetry data. DETIM inputs and outputs were validated with several independent datasets including mass balance data, ATM altimetry, and firn cores covering the entire historical modeling period (1959–2004). The spatial resolution of the mass balance outputs (~60 m over the entire ice cap between 1959–2015 and for Glacier 1 between 2015 and 2100) is a significant improvement over the 11 km resolution provided in an earlier, regional-scale study (Lenaerts and others, 2013). Ice marginal retreat was also accounted for using a volume–area scaling approach. The model in this study performed well when compared to mass balance measurements during the calibration period (r ² = 0.77, RMSE = 0.45 m w.e.), and gave estimates that were identical, within error limits, to those inferred from the 1995–2005 ATM altimetry data.

The majority of predicted surface mass losses on Penny Ice Cap over the period 2015–2099 will occur at elevations between approximately 600 and 1450 m. The refreezing parameterization developed in this study predicts that the ice cap surface will be nearly firn-free by 2100. Considering this loss of refreezing capacity and recent trends in mass balance, Noël and others (2018) suggests the inevitable disappearance of the ice cap. If the average rate of mass loss projected for 2090–99 is assumed to be sustained into the future, we predict that Penny Ice Cap will disappear entirely sometime between the early 2200s for the constant area scaling option, and mid-2400s for the volume–area scaling option.