Firn air content changes on Antarctic ice shelves under three future warming scenarios

The Antarctic firn layer provides pore space in which an estimated 94 % to 96 % of the surface melt refreezes or is retained as liquid water. Future depletion of firn pore space by increased surface melt, densification and formation of low-permeability ice slabs can potentially lead to meltwater ponding, hydrofracturing and ice-shelf disintegration. Here, we investigate the 21st-century evolution of total firn air content (FAC) and accessible FAC (i.e. the pore space that meltwater can reach) across Antarctic ice shelves. We use the semi-empirical IMAU Firn Densification Model (IMAU-FDM) with an updated dynamical densification expression to cope with changing climate forcing. The firn model is forced by general circulation model output of the Community Earth System Model version 2 (CESM2) for three climate emission scenarios (SSP1-2.6, SSP2-4.5 and SSP5-8.5), dynamically downscaled to a 27 km horizontal resolution by the Regional Atmospheric Climate Model version 2.3p2 (RACMO2.3p2). To estimate accessible FAC, we prescribe a relationship between ice-slab thickness and permeability. In our simulations, ice shelves on the Antarctic Peninsula and the Roi Baudouin Ice Shelf in Dronning Maud Land are particularly vulnerable to total FAC depletion (> 50 % decrease by 2100), even for low-emission (SSP1-2.6) and intermediate-emission (SSP2-4.5) scenarios. In the high-emission (SSP5-8.5) scenario in particular, the formation of ice slabs further reduces accessible FAC on ice shelves with low accumulation rates (current rates of < 500 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$), including many East Antarctic ice shelves and the Filchner–Ronne, Ross, Pine Island and Larsen C ice shelves. These results underline the potentially large vulnerability of low-accumulation ice shelves to firn air depletion through ice-slab formation.

1 Introduction

The Antarctic Ice Sheet (AIS) has been losing mass since at least 2002 (The IMBIE team, 2018; Rignot et al., 2019), contributing ∼ 10 % to global average sea level rise since 1993 (Oppenheimer et al., 2019). This mass loss is mainly driven by enhanced iceberg calving and basal melting beneath ice shelves (Smith et al., 2020). Both reduce the buttressing effect of ice shelves, allowing tributary glaciers to accelerate and thereby increasing ice discharge into the ocean. On the Antarctic Peninsula, the warmest region of Antarctica, mass loss is also driven by surface melt. Here, extensive melt has led to meltwater ponding and subsequently to ice-shelf disintegration by hydrofracturing, after which mass loss from tributary glaciers has accelerated (Rignot et al., 2004; Banwell et al., 2013).

Not all ice shelves are susceptible to meltwater-induced hydrofracturing. Firstly, surface melt typically only leads to meltwater ponding when the firn layer lacks sufficient pore space for meltwater to percolate into and refreeze. Currently, an estimated 94 % to 96 % of the surface melt on the AIS is retained within the firn (Medley et al., 2022; Van Wessem et al., 2018). Secondly, meltwater-induced hydrofracturing also requires sufficient tensile stress. Hydrofracturing, in turn, only induces mass loss if the ice shelf provides substantial buttressing. Currently, 60 % of the ice shelves (by area) both buttress upstream ice and are vulnerable to hydrofracturing if inundated by meltwater (Lai et al., 2020). Hence, to assess the future stability of ice shelves and predict mass loss from the AIS, it is important to estimate the future evolution of its firn layer.

Under future warming, we anticipate more surface melt and rainfall, faster firn densification and increased formation of low-permeability ice slabs (i.e. ice layers > 1 m thick) (Ligtenberg et al., 2014; MacFerrin et al., 2019; Kittel et al., 2021; Vignon et al., 2021). These processes decrease the firn air content (FAC) and consequently enable firn saturation and ponding by meltwater. On the other hand, snowfall is projected to increase (Kittel et al., 2021), adding pore space to the firn. Climate models have recently been used to assess the impacts of future climatic changes on Antarctica's firn. Based on the exceedance of a melt-over-accumulation ratio (MOA) of 0.7, van Wessem et al. (2023) estimate future melt ponding on ice shelves. Gilbert and Kittel (2021) use runoff from climate models as an indicator for meltwater ponding.

Offline firn models forced by the output of regional climate models are useful tools to simulate the transient evolution of the firn layer and can therefore also be used to assess meltwater ponding onset. The main advantage of using a firn model instead of a climate model is the lower computational cost, which enables a higher vertical resolution, a proper initialisation of the firn layer and more extensive sensitivity tests. The disadvantage of using an offline firn model is that interaction with the atmosphere is not possible. In contrast to diagnostic studies that use MOA thresholds (e.g. van Wessem et al., 2023), firn models simulate transient changes in the firn layer, thereby accounting for the time it takes to adjust to new climatic conditions. Firn models have been used to simulate the current (1979–present) AIS firn layer (Gardner et al., 2023; Keenan et al., 2021; Medley et al., 2022; Veldhuijsen et al., 2023a) and its evolution in response to climate change (Ligtenberg et al., 2014; Kuipers Munneke et al., 2014a).

The densification equations used in semi-empirical firn models are often based on assumptions of constant accumulation and temperature (Arthern et al., 2010), but simulating future firn evolution requires densification expressions that allow for changing climatic conditions. Moreover, using only total FAC to assess the firn's meltwater buffering capacity overlooks the impact of near-surface ice slabs formed by meltwater refreezing. Ice slabs are common in Greenland (MacFerrin et al., 2019; Culberg et al., 2021) and have been locally observed in Antarctica on the Larsen C Ice Shelf (Hubbard et al., 2016). They impede vertical meltwater percolation into deeper firn, limiting the fraction of the FAC that is accessible to meltwater. To assess the meltwater buffering capacity of firn, it is therefore important to include the impact of ice slabs (henceforth the accessible FAC).

Here, we use the IMAU Firn Densification Model (IMAU-FDM) for Antarctica, which has previously been evaluated for the contemporary climate (1979–2020) (Veldhuijsen et al., 2023a). In this study, IMAU-FDM is driven by realisations of the Community Earth System Model version 2 (CESM2) of the scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5 for the period 1950–2100, dynamically downscaled to 27 km resolution with the Regional Atmospheric Climate Model version 2.3p2 (RACMO2.3p2). To allow for changing climatic conditions, we updated and evaluated the densification equation (Sect. 2). In Sect. 3, we parameterise the accessible FAC based on ice slabs and FAC. In Sects. 4 and 5, we present and discuss the response of the AIS firn layer to future warming scenarios in IMAU-FDM and specifically focus on the evolution of total FAC, accessible FAC and runoff. Finally, we present the conclusions in Sect. 6.

2 IMAU-FDM model updates

IMAU-FDM is a semi-empirical 1-dimensional firn densification model that simulates the evolution of firn depth, density, temperature and liquid water content. Firn compaction is calculated based on the semi-empirical dry-snow densification equations of Arthern et al. (2010). Meltwater percolation is simulated using the bucket method, whereby each firn layer has a maximum irreducible water content that decreases with increasing density (Coléou et al., 1999). The meltwater can percolate through all layers in a single time step and (partly) refreezes when it reaches a layer with a temperature below the freezing point. Once meltwater has saturated the lowermost firn layer beyond the maximum irreducible water content, we assume that it leaves the firn column as runoff instantaneously. Standing water and lateral runoff over ice layers are currently not considered. For further details on the model setup, we refer to Veldhuijsen et al. (2023a) and Brils et al. (2022). Version 1.2A of IMAU-FDM (referred to as FDM v1.2A) has been extensively evaluated over Antarctica against in situ observations of firn density and temperature and remote sensing altimetry measurements (Veldhuijsen et al., 2023a). In this study, we update the model to version IMAU-FDM v1.2AD (referred to as FDM v1.2AD) by implementing a densification expression that is suitable for transient climate change experiments, as is described in the next section.

2.1 Densification expression

The rate of the dry-firn densification ($\mathrm{d}\mathit{\rho }/\mathrm{d}t$) in IMAU-FDM v1.2A is calculated using the semi-empirical equation of Arthern et al. (2010) in combination with a calibration factor (MO_∗):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}={\mathrm{MO}}_{\ast }D\dot{b}g({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho })e^{\left(\frac{E_{\mathrm{c}}}{RT}-\frac{E_{\mathrm{g}}}{R{T}_{\text{ave}}}\right)},\end{array}$$

in which D is a constant, $\dot{b}$ is the long-term average accumulation rate (kg m² yr⁻¹), g is the gravitational acceleration, ρ_i is the density of bubble-free ice (917 kg m⁻³), ρ is the layer density (kg m⁻³), T is the instantaneous layer temperature (K), T_ave is the long-term average surface temperature (K), R is the universal gas constant, and E_c and E_g are the activation energies for creep (60.0 kJ mol⁻¹) and grain growth (42.4 kJ mol⁻¹). The constant D has different values above (0.03) and below (0.07) the critical density level of ρ = 550 kg m⁻³ to represent two distinct densification mechanisms (Herron and Langway, 1980). The calibration factor, MO_∗, depends on annual average accumulation and is defined separately for ρ < 550 kg m⁻³ (MO₅₅₀) and for 550 < ρ < 830 kg m⁻³ (${\mathrm{MO}}_{{\mathrm{830}}^{\ast }}$). These calibration factors are based on the ratio of modelled and observed values of depths of critical density levels ρ = 550 kg m⁻³ (z₅₅₀) and ρ = 830 kg m⁻³ (z₈₃₀), where $z_{{\mathrm{830}}^{\ast }}$ = z₈₃₀ − z₅₅₀. MO₅₅₀ and ${\mathrm{MO}}_{{\mathrm{830}}^{\ast }}$ were chosen as logarithmic and power-law functions, respectively, of the long-term average accumulation rate:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{MO}}_{\mathrm{550}}=\mathit{\alpha }-\mathit{\beta }\ln \left(\dot{b}\right),\text{(3)}& {\displaystyle }& {\displaystyle }{\mathrm{MO}}_{{\mathrm{830}}^{\ast }}=\mathit{\delta }{\dot{b}}^{-\mathit{ϵ}}+\mathit{\varphi },\end{array}$$

in which α, β, δ, and ϵ and ϕ are fit coefficients. Equations (1) to (3) use the long-term annual average accumulation rate ($\dot{b}$) as a proxy for overburden pressure, and Eq. (1) uses long-term annual average temperature (T_ave) in the grain growth part of the exponential term. To include the effect of a changing climate on firn densification, previous studies with IMAU-FDM used running average accumulation and temperature over the 40 years preceding each time step (Ligtenberg et al., 2014; Kuipers Munneke et al., 2014a). Since there is a large spatial variation in firn age across the AIS (e.g. the firn age at the pore close-off depth ranges from 20 to ∼ 3200 years; Veldhuijsen et al., 2023a), the firn temperature and overburden pressure will change at different rates across the AIS in a warming future climate. In addition, the change in temperature and overburden pressure will also differ vertically within the firn column. Advection and conduction transport heat vertically; however, firn has a relatively low thermal conductivity (0.2 to 2 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$; Calonne et al., 2019). For example, Muto et al. (2011) show that firn temperature throughout a firn column differs by about 1 K due to historical temperature trends in East Antarctica of 1 to 1.5 K over 50 years. So, using a running average of the accumulation and firn temperature is a crude approximation to capture the transient response of firn.

In this work, we aim to resolve this by replacing $\dot{b}$ and T_ave in the densification equations (Eqs. 1 to 3) while staying as close as possible to FDM v1.2A. To do so, we revert to the complete transient, dynamical model of Arthern et al. (2010) for the evolution of firn density. Equation (4) is derived from this equation using several simplifying assumptions, such as steady-state accumulation. The equation is based on evolving squared grain size (r²; mm²) and overburden pressure (σ; $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and includes the MO_∗ term:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}={\mathrm{MO}}_{\ast }{k}_{\mathrm{c}}({\mathit{\rho }}_{\mathrm{i}}-\mathit{\rho })e^{\left(\frac{E_{\mathrm{c}}}{RT}\right)}{\displaystyle \frac{\mathit{\sigma }}{r^{\mathrm{2}}}},\end{array}$$

where k_c is a constant. The evolution of r² for one time step includes dry-snow metamorphism:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{r}^{\mathrm{2}}}{\mathrm{d}t}}=k_{\mathrm{g}}{e}^{\left(-\frac{E_{\mathrm{g}}}{RT}\right)},\end{array}$$

where k_g is a constant. In the case of refreezing, r² is calculated as the mass-weighted average of the solid firn grain size and refrozen liquid water grain size, under the condition that refreezing can only increase the grain size. The initial (0.03 mm) and refreezing (0.25 mm) grain sizes are taken from van Dalum et al. (2022). To determine MO_∗ (Eqs. 2 and 3), we replace $\dot{b}$ in those equations by the local (layer) long-term mean accumulation rate (${\dot{b}}_{\text{loc}}$) defined as

$$\begin{array}{}\text{(6)}& {\dot{b}}_{\text{loc}}={\displaystyle \frac{\mathit{\sigma }}{gA}},\end{array}$$

wherein A is the firn age. For updates on the layer merging-and-splitting approach, we refer to Text S1 in the Supplement.

2.2 Atmospheric forcing

IMAU-FDM is forced at the upper boundary with 3-hourly values of snowfall, sublimation, snowdrift erosion, 10 m wind speed, surface temperature, snowmelt and rainfall from RACMO2.3p2. This regional climate model is used to dynamically downscale ERA5 reanalysis data (Hersbach et al., 2020) for the contemporary climate (RACMO2.3p2–ERA5) and Community Earth System Model version 2 model output (CESM2; Danabasoglu et al., 2020) for future projections (RACMO2p3–CESM2) to a 27 km resolution. RACMO2.3p2–ERA5 aims to provide an accurate description of past AIS weather (and climate) from 1979–present and has been thoroughly evaluated (van Wessem et al., 2023).

For future projections, van Wessem et al. (2023) used RACMO2.3p2 to dynamically downscale a historical CESM2 realisation (1950–2014) and one realisation of each of the low-, middle- and high-emission future (2015–2100) projection scenarios (SSP1-2.6, SSP2-4.5 and SSP5-8.5, respectively). CESM2 simulates the coupled interactions between the atmosphere–ocean–land components of the climate system on the global scale at a 1° horizontal resolution (Danabasoglu et al., 2020) and has been thoroughly evaluated over the AIS (Gorte et al., 2020; Dunmire et al., 2022). The model has a relatively detailed representation of polar processes and is among the best CMIP6 models in representing the present Antarctic surface mass balance (SMB; Gorte et al., 2020). The projected Antarctic warming in SSP5-8.5 (+7.7 °C) in CESM2 is stronger than the mean warming in CMIP6 models (+5 °C), which enables us to assess the AIS firn layer response to strong warming.

Arthern et al. (2010)Arthern et al. (2010)

Table 1 Abbreviations and characteristics of IMAU-FDM versions used in this study.

^∗ Veldhuijsen et al. (2023a).

Figure 1 Map of Antarctica showing average firn air content (FAC) for 2014 from FDM v1.2AD-C. The circles indicate locations of in situ observations of the depth of critical density level ρ = 550 kg m⁻³ (z₅₅₀), the depth of critical density level ρ = 830 kg m⁻³ (z₈₃₀) and firn air content (FAC). The blue star indicates the location of Fig. 3b. The names indicate the ice shelves referred to in the text.

RACMO2.3p2–CESM2 time series of mean annual surface temperature, accumulation and surface melt over the AIS are shown in Fig. S1 in the Supplement. Compared to RACMO2.3p2–ERA5 for the overlapping period (1979–2014), we find a temperature bias of −1.2 °C over the AIS, an accumulation bias of −7.6 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (−4 %) and a snowmelt bias of −1.0 mm (−11 %). These biases can be explained by the cold bias in CESM2 (Dunmire et al., 2022). RACMO2.3p2–CESM2 has also been used by van Wessem et al. (2023) in their assessment of AIS meltwater ponding and compares well to meltwater lake volume observations of the Sentinel-2 satellite. Figure S2 in the Supplement shows the differences in mean annual surface temperature, accumulation and surface melt in the period 1979–2014 between RACMO2.3p2–CESM2 and RACMO2.3p2–ERA5. Accumulation and temperature in RACMO2.3p2–CESM2 are in general lower, except in Dronning Maud Land and Enderby Land. Temperatures are most notably lower on the Ross Ice Shelf and in West Antarctica. Surface melt is in general lower, except in Dronning Maud Land and on the Larsen C Ice Shelf.

2.3 Experimental setup

The four IMAU-FDM versions used in this study differ in the dry-snow densification expression (non-dynamical or dynamical; Sect. 2.1) and the indirect atmospheric forcing (ERA5 or CESM2; Sect. 2.2). Abbreviations and characteristics of the IMAU-FDM versions used are listed in Table 1. FDM v1.2AD-C, the dynamical model indirectly forced by CESM2, is used to simulate future firn evolution over the AIS. FDM v1.2A-E and FDM v1.2AD-E are used for the evaluation of the dynamical model over the current climate (Sect. 2.4), and FDM v1.2A-C is used to assess the impact of the dynamical model on future firn evolution (Sect. 4.5; tested for SSP5-8.5). FDM v1.2A-E provides firn characteristics over the AIS from 1979 to the present and has been extensively evaluated (Veldhuijsen et al., 2023a). Details about the model initialisation are given in Text S2 in the Supplement.

Figure 2 Simulated against observed (a) z₅₅₀, (b) $z_{{\mathrm{830}}^{\ast }}$ and (c) firn air content (FAC) for FDM v1.2AD-E, FDM v1.2A-E and FDM v1.2AD-C.

2.4 In situ measurements

To calibrate and evaluate the firn model, we compare the overlapping historical part (1979–2014) of the FDM v1.2A-E, FDM v1.2AD-E and FDM v1.2A-C simulations to in situ firn core density measurements (Fig. 1). We used 112 density profiles across the AIS by combining multiple published datasets (van den Broeke, 2008; Schwanck et al., 2016; Bréant et al., 2017; Fernandoy et al., 2010; Montgomery et al., 2018; Fourteau et al., 2019; Olmi et al., 2021; Winstrup et al., 2019). Detailed information about the dataset is presented in Veldhuijsen et al. (2023a).

2.5 Calibration

The densification equation of the non-dynamical model indirectly forced by ERA5 (FDM v1.2A-E) has been calibrated in Veldhuijsen et al. (2023a). To calibrate the dry-snow densification rate of the dynamical model versions, we first performed simulations for locations with firn density observations (Fig. 1) without MO corrections – i.e. in which the MO values are equal to 1. The resulting MO fits and statistics are shown in Fig. S3 in the Supplement. All model versions yield similar R² values. The MO₅₅₀ fits of the dynamical versions are steeper and higher than FDM v1.2A-E due to a relatively low ratio of overburden pressure and grain size in the upper firn (above z₅₅₀). In addition, the CESM2 forcing alters the MO₅₅₀ and ${\mathrm{MO}}_{{\mathrm{830}}^{\ast }}$ fits due to differences in modelled accumulation rates and surface temperature.

2.6 Performance of the dynamical densification model

We evaluate the performance of the dynamical model, FDM v1.2AD, by comparing it with FDM v1.2A and in situ observations. For the latter, we compare the simulated to the observed depths of the critical density levels ρ = 550 kg m⁻³ (z₅₅₀) and ρ = 830 kg m⁻³ (z₈₃₀) and FAC (Fig. 2a–c). Both versions indirectly forced by ERA5 yield a similar root mean square error (RMSE) and bias, which indicates a similar performance of the updated dynamical firn model for the same forcing. The CESM forced run with the dynamical model yields slightly higher RMSE, indicating that the CESM2 forcing, which is not constrained by observations like ERA5, results in a slightly poorer performance.

The average absolute difference in FAC over the AIS between FDM v1.2AD-C and FDM v1.2A-E for the period 1979–2014 (the overlapping historical period) is 5.2 % (Fig. S4b in the Supplement). This is caused by a combination of the updated dry-snow densification expression and different climatic forcing. Conversely, the absolute FAC difference over the period 1979–2020 between FDM v1.2AD-E and FDM v1.2A-E, caused solely by the updated dry-snow densification expression, averages only 1 % (Fig. S4a). Over the grounded ice, the differences are even smaller (0.6 %), while they are higher in regions with high MOA, such as on Wilkins and Larsen C ice shelves and around the grounding lines of, for example, the Roi Baudouin and Amery ice shelves. The reason is that densification is enhanced in regions with high MOA (Fig. S4c), which is consistent with theoretical considerations that a MOA > 0.7 initiates the depletion of pore space (Pfeffer et al., 1991; van Wessem et al., 2023). Overall, the results of the (non-)dynamical versions are similar over the historical period, which is as expected since there are no trends in accumulation and temperature, reducing the impact of the dynamical densification formulation. The difference in FAC evolution in response to future warming scenarios between the (non-)dynamical versions is discussed in Sect. 4.5.

3 Parameterisation of accessible firn air content

The dynamical FDM version (FDM v1.2AD-C) is used to simulate the future FAC evolution to assess the vulnerability of ice shelves to meltwater ponding. However, meltwater refreezing can form ice slabs in firn, which can impede meltwater percolation into deeper layers, thereby reducing the firn's meltwater buffering capacity (Machguth et al., 2016; MacFerrin et al., 2019; Culberg et al., 2021). That is why, in addition to FAC (henceforth referred to as total FAC), we also calculate accessible FAC. To do so, we use ice layer thickness as a measure of their permeability. It is important to note that the calculation of accessible FAC is a post-processing step and that ice layers within the firn model itself are permeable for vertical liquid water movement.

Figure 3a shows relationships between ice layer thickness and permeability from several observational studies. A small-scale field experiment conducted in Greenland shows that ice layers of 0.12 m can still be completely permeable for liquid water (Samimi et al., 2020). On the other hand, in another small-scale field experiment in Greenland, ice lenses of only 0.03–0.05 m have also been found to be partly impermeable (Clerx et al., 2022). Samimi et al. (2021) assume that ice layers thicker than 0.5 m act as impermeable barriers and prescribe a non-linear decrease in permeability between 0.1 and 0.5 m. These results also show that processes of meltwater percolation and refreezing occur on a small scale. Firn is therefore spatially heterogeneous (e.g. Samimi et al., 2020; Vandecrux et al., 2019), and on a larger horizontal scale, such as that of a FDM grid cell (27 km), ice layers can be discontinuous, allowing meltwater to still percolate through. For ice layers to be impermeable on the model-resolved spatial scale, this requires at least a larger thickness than the ones found in the small-scale field experiments. Here, we also assume that lateral connectivity of ice layers increases with ice thickness.

Radar data and firn cores show that horizontally continuous > 1 m thick ice slabs typically develop on top of refrozen ice layers after extreme melt years. A radar study on Devon Ice Cap, Canada, revealed that an initially widespread ice layer that formed during an extreme melt year thickened by between 0.5 and 4.5 m over the subsequent 5 years (Gascon et al., 2013). Similarly, over time, thickening ice layers of 1–2 m formed in Greenland following an extreme melt year (Culberg et al., 2021). While these large-scale radar observations do not give an exact relation of thickness to permeability, they do give an indication that ice layers thicker than 0.5 m are at least partly impermeable on a larger scale. In addition, firn temperature measurements in Greenland show that no percolation occurred through a 5.5 m thick ice slab even during an extreme melt year (Charalampidis et al., 2016). Based on these observations, we propose a relationship between ice layer thickness (z) and permeability factor (Pf) (Fig. 3a) using a sigmoid function:

$$\begin{array}{}\text{(7)}& \mathrm{Pf}={\displaystyle \frac{\mathrm{1}+b+\left(az\right)}{b+e^{\left(az\right)}}},\end{array}$$

in which a and b are tuning coefficients for which we propose values of a = 1.130 and b = 3.245. We use a sigmoid function because it has a characteristic S-shaped curve, representing that on a large spatial scale, thin ice lenses are permeable, while even very thick ice layers only approach complete impermeability.

Since exact observations that evaluate the permeability as a function of ice layer thickness on the regional scale are lacking, the relationship in Eq. (7) must be regarded as a rough estimate. In addition, the locations of the observations have a larger surface slope (> 0.4°) (Yi et al., 2005) than most Antarctic ice shelves (< 0.15°) (Slater et al., 2018). The small surface slopes on ice shelves result in less lateral flow on top of refreezing layers, which may impact the permeability. Considering these uncertainties, we test the sensitivity of our results to adjusting the values of a (1.119 and 1.250) and b (13.490 and 0.0594), indicated by the envelope of the black-shaded surfaces in Fig. 3a. In addition, we also include the total FAC, which represents full permeability.

To calculate the accessible FAC of a layer, its FAC is multiplied with the permeability factors of all overlying ice layers. For example, if two distinct ice layers with a Pf of 0.5 are overlying a firn layer, the FAC of that layer is multiplied twice by 0.5 to yield its accessible FAC. The sum of accessible FAC of the individual layers equals the firn accessible FAC. An ice layer can range from a single to numerous model layers. Impermeable ice layers are usually defined as having a density of > 830 kg m⁻³ (the pore close-off density). Here, we use a threshold of > 900 kg m⁻³, which corresponds to the density of refreezing ice layers in the model. This choice limits the impact on the accessible FAC of changes in high-density non-refreezing layers in the deep firn. Figure 3b shows the impact of ice layer formation on accessible FAC for an example location in a high-emission scenario (SSP5-8.5). As can be seen, ice layer formation from 2060 onwards depletes the accessible FAC compared to the total FAC. Henceforth, we refer to (a set of) ice layers that have a substantial impact on the accessible FAC as ice slabs.

Figure 3 (a) Relations between the ice layer thickness and permeability factor, including results from in situ observations (stars), a previous estimate (dotted line), observed ranges of reduced permeability from radar observations (coloured shaded surfaces) and a relationship proposed here (solid line), including tested sensitivity ranges (black-shaded surface). (b) Example plot of simulated density including total firn air content (FAC) and accessible FAC (calculated using the proposed relationship from panel a) simulated with FDM v1.2AD-C for SSP5-8.5. The location is on the Shackleton Ice Shelf and is indicated by the blue star in Fig. 1.

Figure 4 Relative change in total firn air content (FAC) (a, b, c) and accessible FAC (d, e, f) by 2100 compared to 2014 for (a, d) SSP1-2.6, (b, e) SSP2-4.5 and (c, f) SSP5-8.5. Please note the different scales for decreasing (red) and increasing (blue) FAC. (g, h, i) Difference between the total FAC change by 2100 (a, b, c) and accessible FAC change by 2100 (d, e, f).

4 Results

4.1 Total and accessible firn air content in 2100 CE

Using the updated firn model FDM v1.2AD, we project the firn evolution over the AIS for the period 1950–2100 driven by RACMO2.3p2–CESM2 for climate scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5. Figure 4a–c show the relative change in total FAC by 2100 compared to 2014. On average, FAC over the grounded ice shifts by −0.7 %, −1.1 % and −2.3 % for scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5, respectively. This implies that the effect of enhanced precipitation is weaker than the effect of firn warming. For ice shelves, FAC shifts by −16 %, −20 % and −42 % for scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5, respectively. For SSP1-2.6 and SSP2-4.5, we find a substantial FAC decrease on Larsen C (−70 % and −75 %), Wilkins (−59 % and −58 %), Roi Baudouin (−67 % and −81 %), George VI (−47 % and −56 %) and Bach (−59 % and −66 %) ice shelves. For SSP5-8.5, we additionally find a substantial FAC decrease for ice shelves in Dronning Maud Land, such as the Fimbul Ice Shelf (−91 %), and elsewhere for the Abbot (−76 %), Pine Island (−74 %), West (−83 %) and Shackleton (−76 %) ice shelves. FAC decrease on the Ross and Filchner–Ronne ice shelves is more modest (−21 % and −19 %).

Figure 5 Relative change in total firn air content (FAC) by 2100 for SSP5-8.5 compared to 2014 (a) as a function of temperature change and accumulation change by 2100 compared to 2014 for locations that do not experience melt by the end of the century in SSP5-8.5 and (b) as a function of annual average accumulation and temperature (1950–2014) for the entire AIS. (c) Average difference between accessible FAC and total FAC (accessible FAC minus total FAC) in 2060–2100 for SSP5-8.5 as a function of annual average accumulation and temperature (1950–2014) for the entire AIS. Contour lines in panel (b) indicate the number of pixels per bin. Please note the different scales for decreasing (red) and increasing (blue) FAC in panel (b).

Figure 4d–f show the relative change in accessible FAC by 2100 compared to 2014, and Fig. 4g–i show the difference between the total FAC change and accessible FAC change by 2100. The figures show that enhanced depletion of accessible FAC by ice layer formation mainly occurs over the ice shelves and is limited over the grounded ice sheet, except for the region east of the Ross Ice Shelf under SSP5-8.5. For SSP1-2.6, the decrease in accessible FAC is accelerated compared to the total FAC on Larsen D ice shelf (from −46 % to −65 %), but, on average, the difference over all ice shelves is limited (−16 % vs. −18 %), indicating limited ice layer formation. For SSP2-4.5, the depletion is accelerated on e.g. the Amery (from 30 % to −49 %), Larsen D (from −49 % to −68 %) and Fimbul (from −43 % to −55 %) ice shelves. The average difference over all ice shelves under this scenario is limited as well (−19 % vs. −23 %). In contrast, for SSP5-8.5, the accelerated accessible FAC depletion is substantial (on average, from −42 % to −53 %) and especially pronounced on e.g. Amery (from −54 % to −92 %), Shackleton (from −76 % to −94 %), Brunt (from −73 % to −93 %) and Filchner–Ronne (from −18 % to −33 %).

4.2 Climatic drivers of changes in firn air content

In this section, we link the changes in FAC under SSP5-8.5 to various climate variables. Due to low temperatures in the interior ice sheet, 76 % of the AIS does not experience melt by the end of the century even in this strong warming scenario. Here, changes in FAC are solely driven by increasing firn temperature and changing accumulation rates. In 42 % of melt-free locations, FAC decreases by 2100, while in the remaining 58 %, FAC increases. In Fig. 5a, we show the relative change in total FAC by 2100 compared to 2014 for non-melt locations as a function of temperature and accumulation change. For 4–6 °C warming, an increase of at least 30 % in accumulation is needed to compensate for the increased densification, whereas for > 8 °C warming, an increase of at least 70 % in accumulation is required. Because of these compensating mechanisms, for most of the AIS, the FAC remains relatively stable.

When we compare the change in FAC of the entire AIS to current climate conditions (1950–2014; Fig. 5b), we observe the largest relative decrease in FAC (−55 %) in currently warm regions (> −22 °C) receiving less than 1000 mm accumulation annually. For high-accumulation regions (> 1000 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$), such as parts of the Getz Ice Shelf and the northwestern part of the Antarctic Peninsula, the decrease in FAC is notably smaller (−23 %). FAC also slightly decreases (−14 %) in colder regions (−34 and −24 °C) with low accumulation (< 100 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$). Large differences over the period 2060–2100 between total FAC and accessible FAC mainly occur in current intermediate warm (−28 to −16 °C) and relatively dry (< 500 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) locations (−0.85 m) (Fig. 5c) and most prominently in regions with temperatures between −22 and −18 °C and accumulation rates between 100 and 400 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (−1.6 m). We select this period as this is when total FAC and accessible FAC start to diverge. For colder initial conditions, the projected melt is too weak to allow for ice lens formation, while for warmer locations, the firn layer is projected to completely disappear in 2090–2100. In the wettest locations (> 600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$), the average difference is smaller (−0.59 m). The patterns in Fig. 5 are, in general, similar for SSP1-2.6 and SSP2-4.5, albeit with smaller magnitudes (Fig. S5 in the Supplement).

Figure 6 Time series of (a, b) total firn air content (FAC), (c, d) accessible FAC and (e, f) runoff extent of 12 ice shelves for (a, c, e) SSP2-4.5 and (b, d, f) SSP5-8.5 for the period 2015–2100. The shaded areas indicate the sensitivity to the relationship between ice layer thickness and permeability factor shown in Fig. 3a.

4.3 Temporal evolution of total and accessible firn air content and runoff over ice shelves

Time series of total FAC, accessible FAC and runoff extent for 12 major ice shelves under SSP2-4.5 and SSP5-8.5 are shown in Fig. 6 (see Fig. S6 in the Supplement for SSP1-2.6). From 2020 onwards there is a gradual decrease in FAC on all ice shelves for all scenarios. An accelerated FAC depletion occurs around 2030 on the Antarctic Peninsula (Larsen C, Wilkins and George VI) ice shelves and Roi Baudouin ice shelf, irrespective of the scenario. Further ahead, for SSP5-8.5, enhanced FAC depletion on the Fimbul, Abbot, Pine Island and Shackleton ice shelves is projected around 2050–2060. Under SSP1-2.6 and SSP2-4.5 the Fimbul Ice Shelf also experiences enhanced FAC depletion around 2050–2060. By 2100, the Antarctic Peninsula ice shelves and the Roi Baudouin Ice Shelf are projected to retain less than 6 m FAC under SSP2-4.5. Under SSP5-8.5, FAC levels are anticipated to decrease by over 4 m across all ice shelves by 2100, resulting in less than 3 m of FAC remaining on the Roi Baudouin, Larsen C, Wilkins, Fimbul and Shackleton ice shelves.

Differences between the accessible and total FAC mainly occur under SSP5-8.5 and are most pronounced on Shackleton, Fimbul, Pine Island, Roi Baudouin, Amery and Filchner–Ronne ice shelves, where depletion is accelerated by up to 20 years and up to 5 m. What these ice shelves have in common is that they are currently relatively dry (< 500 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$). On most ice shelves, the difference increases gradually. However, at the Fimbul and Shackleton ice shelves, ice layer formation around 2060 results in a rapid depletion of more than 5 m (> 50 %) of accessible FAC within 5 years. This coincides with an episode of high melt rates (25 % higher on Fimbul and 38 % higher on Shackleton compared to the previous and following 5 years). On the other hand, on the currently relatively warm (> −19 °C) and wet (> 600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) Wilkins, Getz, Abbot and George VI ice shelves, there remains little difference between accessible and total FAC, which corresponds with Fig. 5c. For SSP2-4.5, there is generally little difference between total FAC and accessible FAC, except for Roi Baudouin, Fimbul and Amery ice shelves (>1 m reduction in accessible FAC compared to total FAC after 2070).

Figure 7 (a, b) The maximum (a) absolute and (b) corresponding relative difference between total firn air content (FAC) and accessible FAC (total FAC minus accessible FAC) from 1950–2100 for SSP5-8.5, shown for locations with at least 25 % accessible FAC depletion in 2100 compared to 2014. (c–h) Time series of total FAC (solid red line), accessible FAC (dashed red line), annual surface melt, annual accumulation and 5-year running average melt-over-accumulation (MOA) ratio for individual grid points. The locations of the grid points are indicated in panel (a). The grey shaded areas indicate extreme melt seasons that cause persistent reduction in accessible FAC.

The runoff extent is the areal fraction of the ice shelf where runoff is generated in a given year. Since ice layers are fully permeable for meltwater percolation in IMAU-FDM, the runoff time series are closely related to the total FAC time series, with some differences (Fig. 6e and f). For the Wilkins, Abbot and George VI ice shelves, the runoff extent is large compared to the FAC depletion (>50 % runoff extent with 5 m FAC) in contrast to the Shackleton and Fimbul ice shelves, where the runoff extents are about 45 %, with 2.5 m of FAC. The reason for the high runoff extent on Wilkins, Abbot and George VI ice shelves is the combination of high accumulation and high melt rates. The high melt allows for saturation of a thick firn layer, which is maintained by high accumulation. In addition, runoff here also occurs year-round from firn aquifers, which are perennial subsurface bodies of liquid water that become more ubiquitous in a warmer Antarctica (Bell et al., 2018).

For the Wilkins Ice Shelf, we see a quick increase from 0 % to > 90 % runoff extent for both scenarios, which indicates limited spatial variation in the firn state across the ice shelf. On the other hand, Larsen C, a larger ice shelf with large latitudinal extent and climate gradients, shows a gradual southward migrating increase in runoff extent. On average, 6 % and 25 % of the entire Antarctic Ice Shelf area experience runoff in 2100 under SSP2-4.5 and SSP5-8.5, respectively, indicated by the black lines in Fig. 6e and f.

Figure 8 (a) Average total firn air content (FAC) difference between FDM v1.2AD-C and FDM v1.2A-C (FDM v1.2AD-C subtracted from FDM v1.2A-C) against average FAC for FDM v1.2A-C for the period 2080–2100. The colour indicates the maximum firn age (taken as the age at the critical density level ρ = 830 kg m⁻³) from FDM v1.2A-C (averaged over 1950–2014). (b) Average FAC difference between FDM v1.2AD-C and FDM v1.2A-C for the period 2080–2100.

4.4 Ice layer formation and its climatic drivers

In the previous sections, we showed that ice layer formation on some ice shelves causes enhanced depletion of accessible FAC. Figure 7a and b, respectively, show the maximum absolute difference and the corresponding relative difference between total FAC and accessible FAC over the period 1950–2100 for SSP5-8.5. The differences are the highest on ice shelves in Dronning Maud Land and on the Amery, West, Shackleton, Ross and Filchner–Ronne ice shelves. In contrast, for the Ross Ice Shelf, the maximum difference is found near the grounding line and, for the Filchner–Ronne Ice Shelf, near the seaward edge. The differences are lowest in the Bellingshausen Sea region, on ice shelves such as Wilkins, George VI, Stange and Abbot and on the Getz and Crosson ice shelves (< −3.9 m and < −47 %), which are among the warmest and wettest ice shelves of the AIS (> −19 °C and > 600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$). The absence of ice slabs under these conditions is also depicted in Fig. 5c.

In Fig. 7c to h, we zoom in on FAC and accessible FAC time series for selected locations plotted together with accumulation, melt and associated MOA. At the Amery, Shackleton and Filchner–Ronne ice shelves, we see that extreme melt seasons can cause a persistent reduction in accessible FAC (indicated by the grey-shaded area). We define a melt season as extreme when the melt exceeds the 93rd quantile of the detrended time series. The 5-year running average MOA values for which this occurs on Amery, Shackleton and Filchner–Ronne are between 0.59 and 0.71. On Larsen C, where melt is more constant over the years, this occurs for a MOA of about 1.04. On the wetter Wilkins location, ice layer formation is limited, even though the firn layer becomes completely depleted and a MOA of 1 is exceeded. In addition, on the wetter Getz Ice Shelf location, ice layer formation is also limited as a MOA of 0.6 has not been reached, even though > 50 % of the FAC has depleted. The general pattern of these figures is that meltwater-blocking ice lenses are primarily formed in drier locations with considerable interannual variability in melt.

4.5 Impact of the dynamical densification model on future firn evolution

The change to a dynamical densification expression changes the temporal evolution of the firn. In this section and in Fig. 8, we quantify this effect by comparing the FAC by 2080–2100 under SSP5-8.5 simulated by FDM v1.2A and by the dynamical model (FDM v1.2AD). We select the period 2080–2100 as the firn layer disappears in 2090 on some ice shelves (Fig. 6b). As shown in Sect. 2.1, both models include the effect of warmer snow densifying faster. However, in FDM v1.2A, the compensating effects of enhanced grain growth, which makes snow stiffer, and enhanced accumulation, which increases the overburden pressure, are parameterised using 40-year running averages of temperature and accumulation.

High accumulation rates result in firn mostly younger than 40 to 100 years, and in such regions, FDM v1.2A provides lower 2080–2100 FAC estimates than the dynamical model. This indicates that the effect of grain size growth, which slows down densification, is underestimated in FDM v1.2A in a transient climate for these locations. Figure 8b shows that these young firn locations are mostly found along the coasts of the Antarctic Peninsula and West Antarctica, such as the Wilkins, Getz and George VI ice shelves. Average end-of-century (2080–2100) FAC estimates of the dynamical model are 3.8 and 1.1 m higher in regions with a maximum firn age lower than 40 and 100 years, respectively.

For locations with older firn (> 100 years), the differences are smaller. This does not prove that the parameterisations of transient behaviour are correct, but at least their errors largely balance out. Clear spatial patterns are visible in Fig. 8b. At the currently relatively dry Ross and Filchner–Ronne ice shelves (105 and 170 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ annual average accumulation, respectively), the dynamical model results in a quicker FAC depletion (−0.45 and −0.41 m, respectively) compared to FDM v1.2A due to a dominating effect of the increase in temperature (+8.1 and +8.4 °C, respectively) compared to accumulation (+10 % and +61 %, respectively) by 2080–2100 compared to 1950–2014.

5 Discussion

5.1 Firn air depletion

We use a firn model with an updated dynamical densification expression to simulate FAC over the AIS under future warming scenarios. We also estimate the impact of ice slabs on the accessible FAC. For the contemporary climate, the performance of the updated model is comparable to the non-dynamical model (Fig. 2). In a warming transient climate (here tested for SSP5-8.5), the Ross and Ronne-Filcher ice shelves in the updated model are more vulnerable to firn air depletion. On the other hand, high-accumulation ice shelves, such as Getz, Wilkins and George VI, are less vulnerable (Fig. 8).

The high vulnerability of the Antarctic Peninsula ice shelves and the Roi Baudouin Ice Shelf to total FAC depletion (> 50 % decrease in FAC by 2100 for all scenarios) aligns with previous studies (Gilbert and Kittel, 2021; Kuipers Munneke et al., 2014a; van Wessem et al., 2023). According to our results, enhanced FAC depletion on these ice shelves is expected to start around 2030 for all scenarios (Figs. 6 and S6). For SSP5-8.5, we also find a substantial FAC depletion (> 74 % decrease) for ice shelves in Dronning Maud Land, such as Fimbul, and elsewhere for Abbot, Pine Island, West and Shackleton ice shelves. Enhanced FAC depletion on these ice shelves is expected to occur around 2050–2060 for SSP5-8.5.

By including the effect of reduced permeability of ice slabs, our results reveal enhanced accessible FAC depletion under SSP5-8.5 on ice shelves with current accumulation rates of < 500 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$, mean annual temperatures of < −16 °C and interannual variability in melt. Such ice shelves are mainly situated in East Antarctica but also include the Ross, Filchner–Ronne, Larsen C and Pine Island ice shelves. Enhanced depletion starts around 2060–2070 for Shackleton, Fimbul and Pine Island shelves (∼ 4 °C AIS-wide average warming) and around 2080 for the colder Ross and Filchner–Ronne ice shelves (∼ 6 °C AIS-wide average warming). For SSP1-2.6 and SSP2-4.5 the accelerated accessible FAC depletion is limited to Larsen D, Fimbul, Amery and Roi Baudouin ice shelves. The formation of ice slabs thus increases the divergence between firn air depletion for the high-emission scenario on the one hand and for low- and intermediate-emission scenarios on the other hand.

Our results suggest that extreme melt seasons can initiate ice-slab formation (Fig. 7), as observed in Greenland (Culberg et al., 2021). When ice layers that formed during an extreme melt season remain in contact with the surface hydrology, they concentrate new refreezing above their horizon, amplifying the contribution of even average subsequent melt seasons to ice-slab formation (Jullien et al., 2023). However, standing water on top of ice layers is not considered by IMAU-FDM, but an alternative positive feedback mechanism is captured: ice has a higher thermal conductivity, which contributes to more efficient conductive cooling and thereby promotes the growth of ice layers and slabs (Vandecrux et al., 2020). The future accessible FAC on these ice shelves will depend on the frequency of future extreme melt events, and the timing of depletion is therefore less certain compared to non-ice-slab regions, where the depletion is mainly forced by the mean climate state and therefore more gradual.

Ice shelves where ice-slab formation under SSP5-8.5 is prevented have current accumulation rates > 600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ and mean annual temperatures > −19 °C and are mainly situated on the Antarctic Peninsula and West Antarctica, such as Wilkins, Stange, Abbot, Crosson and Getz ice shelves. High accumulation allows new pore space to regenerate FAC above the most recent refreezing layers. In addition, mild temperatures and more low-density fresh snow, both reducing the refreezing capacity per m³ of firn, prevent firn from reaching the ice density. The projected increasing snowfall in Antarctica is thus important not only to counteract ice sheet mass loss, but also to prevent ice-slab formation. This is in line with MacFerrin et al. (2019), who found that ice slabs appear to be absent in regions of high accumulation (> 572 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) on the Greenland Ice Sheet. In contrast, high accumulation (> 1000 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) and mild (> −19 °C) conditions may lead to the formation of aquifers (Kuipers Munneke et al., 2014b; van Wessem et al., 2021), which aligns with the absence of extensive refreezing. Although firn aquifers do not have a depleted (accessible) FAC, the runoff from aquifers has the potential to cause hydrofracturing. Therefore, we will explore the future expansion of AIS aquifers in a forthcoming study.

The AIS ice-shelf-wide runoff extent in Fig. 6, which is 25 % by 2100 for SSP5-8.5, is substantially lower than reported by Gilbert and Kittel (2021), who estimate a runoff extent of 98 % by 2100 using the regional climate model MAR forced by the same CESM2 realisation. In MAR, liquid water reaching a firn layer with a density of > 830 kg m⁻³ is converted to runoff, whereas in IMAU-FDM, ice layers are completely permeable for meltwater percolation. Both models thus represent extremes in permeability assumptions, which explains the large difference. This also underlines the potentially large impact of ice slabs when estimating the onset of meltwater ponding. Besides this, MAR also yields higher melt rates compared to RACMO2 for the same forcing (Carter et al., 2022), which is potentially caused by a more active snowmelt–albedo feedback. Differences in the initialisation of the firn layer could also contribute to the discrepancy.

Extensive meltwater ponding is expected when accessible FAC is strongly depleted. Figure 6 shows that runoff is enhanced when total FAC (averaged over the ice shelf) roughly approaches 2.5 m for low-accumulation ice shelves, such as the Roi Baudouin Ice Shelf, and 5 m for high-accumulation ice shelves, such as the Wilkins Ice Shelf. We expect that the accessible FAC thresholds for runoff generation are approximately similar. Based on these assumptions, Fig. 6 implies an extensive melt ponding onset around 2060 on the Roi Baudouin Ice Shelf for all scenarios. On Larsen C and Wilkins ice shelves, meltwater ponding is expected around 2060 for SSP5-8.5 and around 2080 for SSP1-2.6 and SSP2-4.5. For these three ice shelves, meltwater ponding is thus expected to occur in the 21st century irrespective of the emission scenario. For SSP5-8.5, extensive meltwater ponding is expected to occur on George VI around 2070; on Fimbul around 2075; and on the Shackleton, Amery and Pine Island ice shelves around 2090. However, extensive meltwater ponding has already been observed on the northern George VI Ice Shelf (van Wessem et al., 2023), which suggests that our initial FAC estimates (Fig. 1) may be overestimated. On the other hand, locally observed meltwater ponding along the grounding lines of Amery and Roi Baudouin ice shelves is in line with our low initial FAC estimates there.

5.2 Limitations

FDM v1.2AD uses the bucket scheme to simulate vertical meltwater movement, which does not allow for standing water over impermeable ice slabs, lateral meltwater movement or preferential flow. Neglecting these processes has impacts on densification and ice layer formation. Verjans et al. (2019) found that a firn model using the bucket scheme produces similar density results as the more physics-based Richards equation in a single domain for four locations on the Greenland ice sheet. Using the Richards equation in a dual domain, which accounts for partitioning between matrix and fast preferential flow, simulated more ice layers and underestimated the FAC to a greater extent but is better at reproducing density variability with depth. Using two firn models that use the bucket scheme, Thompson-Munson et al. (2023) simulated ice slabs in Greenland and found that most ice slabs detected by IceBridge accumulation radar data in 2014 (MacFerrin et al., 2019) overlap with simulated ice slabs and ablation zones. On average, the observed ice slabs are located at slightly higher altitudes than the modelled ones. However, the results are also strongly influenced by the densification equation and climatic forcing used (Verjans et al., 2019).

The densification scheme used in FDM v1.2AD (Eq. 4) is developed for dry-snow densification. The presence of liquid water may also impact the densification rate of firn; but, due to a lack of physical understanding and available measurements, this effect has not been included. In line with this, the in situ observations used for the calibration and evaluation of the model are mainly (90 %) situated in non-melt regions. Equation (4) depends on initial grain size (taken here as 0.03 mm), which is a poorly constrained parameter over the range of climates on ice sheets. Fortunately, the model grain size is not very sensitive to the initial value, since the growth rate is independent of the initial grain size, and the grain size quickly becomes magnitudes larger (e.g. a growth rate of 0.0056 mm d⁻¹ for a firn layer of 250 K). Observed and modelled albedo evolution through melt–refreezing cycles shows that refreezing strongly increases the grain size (e.g. van Dalum et al., 2022). However, tests with a different refreezing grain size (e.g. 0.4 mm instead of 0.25 mm) did not lead to very different results (Fig. S7 in the Supplement). Again, this process is poorly constrained by direct observations, and hence the impact of refreezing on grain size and subsequently on firn compaction remains uncertain. Firn compaction by horizontal divergence and strain softening is not included in the model; therefore, densification in high-strain and horizontal stretching areas, such as Pine Island and Thwaites ice shelves, is likely underestimated (Horlings et al., 2021; Oraschewski and Grinsted, 2022).

Our outcomes of accessible FAC depend on ice-slab formation, which is a complex process to model. The current approach (as a diagnostic post-processing step rather than including it in the model) should be regarded as an exploratory study. We present these results together with the total FAC and test the sensitivity to several ice layer thickness–permeability relations. We found that the timing and magnitude of accessible FAC depletion averaged over ice shelves are not very sensitive (< 10 years and < 1 m) to the thickness–permeability relationships used (Figs. 3 and 6). Including the impermeability of ice layers interactively within IMAU-FDM will be tested in future work.

The current view is that meltwater ponding weakens ice shelves, potentially leading to ice-shelf disintegration by hydrofracturing (Banwell et al., 2013). However, Bell et al. (2017) found that when the ice surface slope is sufficiently steep, rivers can form that transport liquid water laterally to the ocean, thereby preventing its destructive effects on ice shelves. FDM does not include a liquid water routing scheme to transport water laterally by streams and rivers.

Our results are furthermore limited by only a single general circulation model (CESM2) and one model for dynamical downscaling (RACMO2). Future work should address this by an ensemble of Earth system models and other regional climate models such as MAR (Kittel et al., 2021) and MetUM (Orr et al., 2021), which have different outcomes for melt (Carter et al., 2022). The exact timing and amount of (accessible) FAC depletion will vary when using a different forcing dataset; however, the broad spatial variability among climatic regions and the dependency of ice layer formation on accumulation should not be fundamentally different.

6 Conclusions

In this study, we explored future (accessible) FAC evolution over the AIS under various emission scenarios. Our main tool is a firn model with an updated densification expression that allows for changing climatic conditions. Compared to the previous model version, the Ross and Filcher–Ronne ice shelves emerge as being more vulnerable to FAC depletion, while high-accumulation ice shelves, such as Getz, Wilkins and George VI, are less vulnerable. In our simulations, ice shelves on the Antarctic Peninsula and the Roi Baudouin Ice Shelf in Dronning Maud Land are particularly vulnerable to FAC depletion, also under low-emission (SSP1-2.6) and intermediate-emission (SSP2-4.5) scenarios (> 50 % decrease in FAC by 2100). In the high-emission (SSP5-8.5) scenario, all ice shelves experience considerable FAC depletion by 2100 (> 19 % decrease). The formation of ice slabs further reduces accessible FAC under SSP5-8.5, especially on low-accumulation ice shelves (current rates < 500 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) with a considerable interannual variability in surface melt. This includes many East Antarctic ice shelves and the Filchner–Ronne, Ross, Pine Island and Larsen C ice shelves. Under emission scenarios SSP1-2.6 and SSP2-4.5, ice-slab formation is limited to the Larsen D, Roi Baudouin, Fimbul and Amery ice shelves. High accumulation rates (currently > 600 $\mathrm{mm}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$) on the Antarctic Peninsula and West Antarctic ice shelves, such as the Wilkins, George VI, Getz, Abbot, Stange and Crosson ice shelves, prevent the formation of ice slabs under all scenarios. These results highlight the potentially large vulnerability of low-accumulation ice shelves to firn air depletion through ice-slab formation.