The impact of ice structures and ocean warming in Milne Fiord

Arctic tidewater glaciers and ice shelves are undergoing rapid attrition, with warmer ocean temperatures playing an important role. However, the relationship between ocean temperature and ice structure retreat is complex and may change as the ocean warms and as the ice structure geometry evolves. In order to explore ice–ocean interactions and the impact of retreating ice structures in a glacial fjord, we use a numerical ocean model of Milne Fiord, which features an ice shelf and a tidewater glacier with a floating glacier tongue (part of which is detached). We model past, present, and potential future ice configurations. Our results reveal that the average submarine melting is negligible across the ice shelf (<2 cm a⁻¹) but can dominate thinning rates (>20 cm a⁻¹) at specific locations where the ice is thick (>50 m) along the seaward edge. Our simulations also indicate that the temperature of water reaching the grounding line does not vary significantly when the ice shelf and glacier tongue are removed. In addition, we carry out a series of simulations with increasing ocean temperature which reveal a quasi-linear relationship between ocean temperature and submarine melting at the grounding line. Using this relationship and ocean temperature predictions for different greenhouse gas emission scenarios (2020 to 2100), we estimate that Milne Glacier will continue to retreat for at least 50 years, solely in response to ocean forcing. This study highlights the vulnerability of ice structures in the Arctic, even in a region regarded as the Last Ice Area.

1 Introduction

Ice structures, defined here as large submerged glacier ice features, are currently experiencing a global retreat in the Arctic; 85 % of marine-terminating glaciers (hereafter “tidewater glaciers”) in the Northern Hemisphere retreated between 2000 and 2020 (Kochtitzky and Copland, 2022). During this period, many glaciers also lost their ice tongues, including the Jakobshavn (Motyka et al., 2011), Zachariae Isstrøm (Mouginot et al., 2015), and Yelverton glaciers (White and Copland, 2019). A warmer ocean is often considered a major cause of this phenomenon (e.g. Holland et al., 2008; Straneo et al., 2012; Cai et al., 2017; Millan et al., 2022). In the Southern Hemisphere, the intrusion of warm modified Circumpolar Deep Water inside ice shelf cavities has also been identified as a major cause of ice shelf retreat (e.g. Azaneu et al., 2023; Si et al., 2024).

It is an axiom that increased water temperature in contact with ice leads to more melting, but the exact relationship between the temperature of water offshore of glacial fjords or ice shelf cavities and the melting of ice structures is complex. This relationship depends on multiple factors, such as the presence of a sill (Zhao et al., 2019, 2021), the presence of icebergs (Davison et al., 2020; Hager et al., 2024), ice shelf thickness and geometry (De Rydt et al., 2014; Bradley et al., 2022; Steiger et al., 2022), and the amount of subglacial discharge (Cowton et al., 2015; Slater et al., 2019). For instance, for four similar tidewater glaciers in West Greenland, Rignot et al. (2016) found that the dependency of glacier melting on water temperature can vary by more than 50 % (their Table S1). This result suggests that simple parameterizations for ice melting as a function of water temperature derived for a specific site should be used with caution for other glacial fjords or ice shelf cavities, unless validated. Moreover, in addition to the variability between sites, changes at a specific location, such as glacier retreat or a large calving event, can potentially change the derived relationship by altering the ocean circulation and hydrography. This process is important because it could create a positive feedback loop where warmer water leads to enhanced calving, which leads to more warm water intrusion. A retrograde bathymetry, such as that found in Petermann Fjord, is an example of a feature potentially leading to enhanced submarine melting following glacier retreat because of the increasing water temperature with depth (Millan et al., 2022). On the other hand, the emergence of a shallow sill during glacier retreat could obstruct the flow of warmer deep water to the grounding line, as is currently observed with the Ryder and 79N glaciers (Schaffer et al., 2020; Jakobsson et al., 2020). Two recent numerical studies on Antarctic ice shelves indicate that major calving events do not substantially change the overall average submarine melting rate (<15 % change) but can lead to significant melt increases at specific locations (>100 %) (Bradley et al., 2022; Poinelli et al., 2023).

In this study, we evaluate the impact of changing ice structures (ice shelf, glacier tongue) and increasing ocean temperature on ice–ocean interactions in Milne Fiord. Milne Fiord, located on the northern coast of Ellesmere Island in the Canadian Arctic, is an ideal site because it hosts an ice shelf at its mouth in addition to a tidewater glacier at its head (separate ice structures). Moreover, observational efforts in the last 15 years enabled the development and validation of a high-resolution ocean model, which can predict the melting/freezing of the ice structures (Bonneau et al., 2024b). The model is firstly used to evaluate the impact of removing ice structures (calving) on the hydrography and submarine melting. Then, it is employed to evaluate the impact of the forecasted ocean warming in the region on the submarine melting of the glacier.

Figure 1 Study area. (a) Earthview image of the northern coast of Ellesmere Island taken on 18 August 2024. Ice shelf extent from Mueller et al. (2017). The black box is the extent of panel (b), and the blue box is the extent of the numerical model domain. The top-right inset shows the location of panel (a). (b) PlanetScope near-infrared image of Milne Fiord from 18 August 2024 showing the extent of the Milne Ice Shelf (MIS) and the Milne Glacier tongue (MGT). (c) Along-fjord profile of ice thickness from ArcticDEM (Porter et al., 2018) and IceBridge (Paden et al., 2010) derived using hydrostatic equilibrium (for ArcticDEM data). Seabed elevation under the glacier is from IceBridge, and seabed elevation under the ocean is from manual soundings.

2 Background on the region

Milne Fiord (82.5° N, −81° W) is situated on the northernmost coast of the Canadian Arctic Archipelago, facing the Arctic Ocean (Fig. 1a). This unique site hosts an ice shelf, a detached glacier tongue, and a tidewater glacier. It is located at the centre of the coastline of Canada's newly designated Tuvaijuittuq Marine Protected Area. Tuvaijuittuq, meaning “the place where the ice never melts” in Inuktut, is projected to be one of the last areas with multiyear sea ice in the Northern Hemisphere (Jahn et al., 2024), thereby serving as a refuge for ice-dependent species (Vincent and Mueller, 2020). However, in spite of the predicted resilience of sea ice in this area, the ice structures (ice shelves, glacier tongues, glaciers) in Tuvaijuittuq are currently experiencing sustained retreat (e.g. Mueller et al., 2017; Kochtitzky and Copland, 2022; White and Copland, 2019).

2.1 Milne Ice Shelf (MIS)

Using explorer reports from expeditions in 1875–1876 and 1906, Vincent et al. (2001) estimated that the north coast of Ellesmere Island was fringed by an ∼8400 km² ice shelf at the beginning of the 20th century. This “Ellesmere Ice Shelf” (unofficial name) is believed to have formed around 5500 years ago (England et al., 2017). Over the past century, the Ellesmere Ice Shelf has almost completely disintegrated; it had an extent of 535 km² in 2015 (Mueller et al., 2017). The last relatively stable fragment of the Ellesmere Ice Shelf was in Milne Fiord, until 2020, when the Milne Ice Shelf (MIS) calved and lost 43 % (80 km²) of its area, including its thickest ice (∼90 m thick). Before the calving, the mean thickness of the MIS was ∼44 m (Fig. 1c; Hamilton, 2016). Ice cores and remote sensing have revealed that the MIS is a composite ice shelf, formed by the floating extension of glaciers and the accumulation of sea ice, nourished by basal accretion and snow accumulation (Jeffries, 1992; Mortimer et al., 2012; Richer-McCallum, 2015; White, 2019). Since 2011, the only glacier feeding the MIS is Glacier 2 (Fig. 1b; Mortimer et al., 2012)). Until 2020, the MIS acted as a dam for surface runoff from the watershed, forming a freshwater layer floating above the seawater known as an epishelf lake (Hamilton et al., 2017; Bonneau et al., 2021; Veillette et al., 2011). The Milne Fiord epishelf lake drained through a basal channel under the MIS, where the ice draft was approximately 8 m thick (Bonneau et al., 2021). Although the epishelf lake has dramatically thinned since the 2020 calving event, it is unclear at the moment if it remains perennial or has transitioned to a seasonal feature.

2.2 Milne Glacier

Around 1959, the floating extension of Milne Glacier disconnected from the MIS and retreated up-fjord, creating a gap between the ice shelf and the Milne Glacier tongue (MGT) (Jeffries, 1986). This gap, filled by the epishelf lake and covered by perennial lake ice, has been expanding ever since (Mortimer et al., 2012). Due to the presence of the ice shelf and perennial ice, icebergs do not exit the fjord, resulting in a fractured glacier tongue, similar to the now-disintegrated tongues of the Tracy, Steensby, and C.H. Ostenfeld glaciers (Hill et al., 2017) and to the current Thwaites Western Ice Tongue (e.g. Benn et al., 2022). Between 2009 and 2011, a large rift formed <1 km from the grounding line (Mueller et al., 2017; Antropova et al., 2024), separating the MGT from the glacier. In 2018, the MGT covered ∼66 km². A recent study combining remote sensing and radar observations from 2011 to 2023 estimated that the grounding line of Milne Glacier retreated between 30 and 124 m a⁻¹, depending on the location across the glacier (Fig. 1b; Antropova et al., 2024). The Milne Glacier grounding line has a maximum depth around 190 m and was located on a retrograde slope in 2011 (Fig. 1c). The estimated ice flux from the glacier varied from 0.03 to 0.14 Gt a⁻¹ between 2011 and 2019. The associated glacier velocities ranged from 20 to 160 m a⁻¹ between 2011 and 2020, with the acceleration of up to ∼160 m yr⁻¹ between 2016 and 2020 (Van Wychen et al., 2016; Millan et al., 2017; Wychen et al., 2020; Antropova et al., 2024).

2.3 Milne Fiord oceanography

The water column in Milne Fiord is composed of different water masses. At the surface is the ∼10 m thick epishelf lake, which is warm (>0 °C) and fresh (<1 g kg⁻¹) (Hamilton et al., 2017; Bonneau et al., 2021). It is ice-covered year round and therefore isolates the fjord's seawater from atmospheric forcing. Below the epishelf lake is an ∼100 m thick layer of polar water near the freezing point. Below the polar water is the warmer (>0 °C) and saltier Atlantic Water. Because of the presence of the MIS, freshwater from surface runoff, subglacial discharge, and submarine melting accumulates in the upper 50 m of the water column (Hamilton et al., 2021). An ∼230 m sill restricts the renewal of deep water, but it has limited impact on water reaching the glacier because it is deeper than the grounding line (Hamilton et al., 2021). The results from a multiyear numerical simulation of the circulation in Milne Fiord, validated with observations, have revealed that the circulation is weak but highly unsteady: three circulation modes with different patterns alternate in response to density variations on the coastal shelf (Bonneau et al., 2024b). The water on the coastal shelf, which flows westward, originates mainly from the Canada Basin (Hamilton et al., 2021; Timmermans and Marshall, 2020). Except for intermittent leads along the coast in July, August, and September (Fig. 1a), the offshore region is covered by multiyear land-fast sea ice.

3 Methods

3.1 Numerical model

This study relies on a numerical ocean circulation model (MITgcm; Marshall et al., 1997) of Milne Fiord developed and evaluated in Bonneau et al. (2024b, a). The model covers the extent of the fjord using a 150 m by 150 m horizontal (Fig. 1a). The numerical domain extends 25 km east, west, and north offshore from the fjord's mouth, with a telescoping grid. In the vertical, the cells are 2 m for the first 50 m and then linearly extend to 30 m at 600 m depth. Vertical eddy viscosity and diffusivity are set to 10⁻⁵ and 10 m² s⁻¹, respectively. Horizontal eddy viscosity is determined using the Smagorinsky method (Smagorinsky, 1963) (coefficient = 2.0), and the horizontal diffusivities are set to 3.0 m² s⁻¹ (Table 1). To model the heat, salt, and momentum flux at the boundaries with the ice structures, we take advantage of the separation of the glacier tongue from the glacier in 2009–2011. The ice shelf and the glacier tongue are modelled using the ShelfIce package (Losch, 2008) with a drag coefficient (C_d) automatically determined so that u and v are zero at the interface (corresponding to a no-slip boundary condition). At the glacier face, the fluxes are determined with the IcePlume package (Cowton et al., 2015) (C_d=0.01, heat transfer coefficient Γ_Θ=0.044, salt transfer coefficient Γ_S=0.00124, minimum background velocity v_hor=0.035 m s⁻¹), modified for a glacier with a constant overhanging slope of 4:1 (Bonneau et al., 2024b). The employed values for C_d, Γ, and v_hor are in line with Jackson et al. (2020) and yield glacier retreat estimates consistent with observations (Sect. 4.2). The amount of subglacial discharge is determined by integrating the negative surface mass balance from RACMO2.3 (Noël et al., 2018). Two subglacial discharge outlets are used (in the IcePlume package), both coincident with depressions in the bathymetry (Sect. 4.2). While the exact number of subglacial discharge outlets is unknown, observations showed one outlet in a bathymetry depression on the west side of the fjord (Hamilton et al., 2021). Another one was therefore added in a depression on the other side of the fjord. Subglacial discharge is released equally in the two outlets. Offshore, the model is forced with temperature, salinity, and velocities from ORAS5 reanalysis (Zuo et al., 2019). To reproduce the epishelf lake, which recharges every summer because of surface runoff, the upper three grid cells were relaxed towards mooring observations (temperature and salinity).

Table 1 Parameters used in the MITgcm ocean model.

Figure 2 Model configuration. (a) Ice draft and ice extent for the different ice structure scenarios. (b) Mean offshore ocean temperature for the different ocean warming scenarios. (c) Time series of amount of subglacial discharge (same for all scenarios).

3.2 Numerical experiments

This study is divided into two main parts. The first part explores the impact of the ice structures in Milne Fiord on water properties and melting/freezing of the remaining ice structures. For this, the numerical model is run with four different ice structure configurations (Fig. 2a). The first ice configuration is the same as in Bonneau et al. (2024b, a) and represents Milne Fiord before the July 2020 MIS calving event; it is denoted “pre2020”. The second part represents the fjord as of 2024; it is denoted “now2024”. The third ice configuration represents the fjord without an ice shelf; it is denoted “nois”. The fourth configuration represents the fjord without an ice shelf and without a glacier tongue; it is denoted “nogt”. nois and nogt are configurations expected to arise in the coming decades, as has happened for other fjords along the northern coast of Ellesmere Island over the last 2 decades (Copland et al., 2007; Mueller et al., 2003; White and Copland, 2019; Mueller et al., 2017). For the nois and nogt runs, there is no temperature and salinity relaxation in the three uppermost grid cells to reproduce the epishelf lake. For each configuration, the model is spun up for 0.5 years starting in May 2011, and the last 3 years (November 2011 to November 2014) are used for analysis. A 3-year simulation is sufficient to obtain accurate (<20 % error) average along-fjord heat and salt fluxes and accurate mean circulation patterns (Bonneau et al., 2024b).

The second part of this study explores the impact of warming ocean temperatures on the melting of the glacier face. For this purpose, the model is run with the nogt ice configuration and with four different temperature increases (added to the baseline nogt ORAS5 offshore forcing). The temperature increases were determined using a multimodel mean from the fifth and sixth Coupled Model Intercomparison Project (CMIP; Taylor et al., 2012; Eyring et al., 2016). The four ocean warming simulations (T03, T09, T16, and T30) are intermediate steps in different greenhouse gas emission scenarios (RCP2.6/SSP126, RCP4.5/SSP245, and RCP8.5/SSP585; S1). The CMIP multimodel mean ocean temperature increase corresponding to a warming of 0.3, 0.9, 1.6, and 3.0 °C at 220 m depth was added to the ORAS5 baseline offshore forcing (nogt). 220 m corresponds to the grid cell just above the sill and just below the grounding line of Milne Glacier. Above 220 m, the temperature increase decreases with height for all scenarios, as CMIP models do not show a temperature increase at the surface, except for extreme scenarios (Fig. 2b). The amount of imposed subglacial discharge at the grounding line of Milne Glacier is the same for all cases (Fig. 2c). For each of the ocean warming simulations, the model is run 2.5 years and the last 2 years are used for analysis. A time of 2 years is sufficient because the melting of the glacier face does not show high interannual variability (Sect. 4.2); therefore this duration is sufficient to obtain an accurate average.

3.3 Water property metrics

To characterize the water properties, the water column inside the fjord is divided into four distinct layers and the water property anomalies, and the differences between the different ice structure configurations in the simulation are analysed. Temperature anomalies are defined as the temperature inside the fjord minus the temperature offshore (Θ−Θ_off), and salinity anomalies are defined likewise (S−S_off). We do not consider the depth of the sill here, as the grounding line is shallower than the sill's depth. Unless stated or plotted otherwise, the water properties at the centre (18 km down-fjord from the grounding line) are used for $\mathrm{\Theta }/S$ when calculating these anomalies. We also calculate the volume exchange between the coastal shelf and the fjord (Q_ex), the pycnocline upwelling inside the fjord (Δz_ρ), the fraction of submarine meltwater (f_SM), and the total amount of freshwater (f_SM+SD). Q_ex is calculated using along-fjord velocities (U):

$$\begin{array}{}\text{(1)}& Q_{\mathrm{ex}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\int \left|U\right|\mathrm{d}A,\end{array}$$

where A is the cross-section at the mouth of the fjord (X=33 km). The upwelling, Δz_ρ, is calculated as the average depth difference between isopycnals offshore and at the centre of the fjord (X=20 km) for a specific layer. Below the settling depth of the subglacial discharge plume (∼40 m), the fraction of submarine meltwater (f_SM) is estimated using the temperature difference between the fjord (Θ) and offshore along isopycnals, such that ice melting produces a cold anomaly represented by the mixing of ambient offshore water (Θ_off,ρ) with water at −92.5 °C (Θ_SM). The temperature of −92.5 °C accounts for the latent heat of fusion (Hamilton et al., 2021):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{f}_{\mathrm{off},\mathit{\rho }}+f_{\mathrm{SM}}& =\mathrm{1}\cup f_{\mathrm{off},\mathit{\rho }}{\mathrm{\Theta }}_{\mathrm{off},\mathit{\rho }}+f_{\mathrm{SM}}{\mathrm{\Theta }}_{\mathrm{SM}}=\mathrm{\Theta }\\ & \to f_{\mathrm{SM}}={\displaystyle \frac{\mathrm{\Theta }-{\mathrm{\Theta }}_{\mathrm{off},\mathit{\rho }}}{{\mathrm{\Theta }}_{\mathrm{SM}}-{\mathrm{\Theta }}_{\mathrm{off},\mathit{\rho }}}},\end{array}\end{array}$$

where f_off,ρ is the fraction of offshore water. Θ_off,ρ is taken at the same density level as Θ but offshore. Above the minimal settling depth of the subglacial discharge plume, the combined fraction of submarine meltwater and subglacial discharge water (f_SM+SD) is estimated using the fact that both subglacial discharge and submarine melting waters have a negligible salinity ($S_{\mathrm{SM}+\mathrm{SD}}\approx \mathrm{0}$ g kg⁻¹). Therefore, f_SM+SD can be determined using the salinity inside the fjord (S_obs):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{f}_{\mathrm{off},z}+f_{\mathrm{SM}+\mathrm{SD}}& =\mathrm{1}\cup f_{\mathrm{off},z}{S}_{\mathrm{off},z}+f_{\mathrm{SM}+\mathrm{SD}}{S}_{\mathrm{SM}+\mathrm{SD}}=S\\ & \to f_{\mathrm{SM}+\mathrm{SD}}=\mathrm{1}-{\displaystyle \frac{S}{S_{\mathrm{off},z}}}.\end{array}\end{array}$$

The offshore salinity (S_off,z) is taken at the same depth as S. Since subglacial discharge entrains saltier deep waters and transport them upwards, Eq. (3) underestimates f_SD+SM. Similarly, Eq. (2) also underestimates f_SM below 40 m since the neglected subglacial discharge water and associated upwelled deep water result in warm anomalies, offsetting part of the cold anomaly used to determine f_SM. Nevertheless, these estimates satisfy the present objective of comparing water properties for the different ice structure configurations. The amount of f_SM and f_SM+SD is reported in water height (∫fdz). For example, an average fraction of 0.01 over a 10 m layer is 10 cm.

3.4 Glacier face melting parameterization

To compare the results from our model to other modelled glacial fjords, the submarine melting parameterization of Rignot et al. (2016) is employed:

$$\begin{array}{}\text{(4)}& M_{\mathrm{GL}}=(Ah{q}_{\mathrm{sd}}^{\mathit{\alpha }}+B){\text{TF}}^{\mathit{\beta }},\end{array}$$

where M_GL (in units of m d⁻¹) is the horizontal submarine melting rate at the grounding line (undercutting), h is the thickness of the glacier at the grounding line (200 m), and q_sd is the rate of subglacial discharge volume (in units of m³ d⁻¹) divided by the submerged area of the glacier face. TF is the temperature above the local freezing point (Θ_f), commonly referred to as the thermal forcing (Θ−Θ_f). In this study, the TF value used in Eq. (4) is the average TF value from the sea surface to the grounding line depth, following Morlighem et al. (2019). Equation (4) is not used in our model but on the output. Using remote sensing observations to calibrate Eq. (4) for five tidewater glaciers in West Greenland, Rignot et al. (2016) arrived at mean values for the coefficients A, α, B, and β of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$, 0.39, 0.15, and 1.18, respectively. These coefficients have not been evaluated in other glacial fjords, and they differed by up to 50 % between fjords in the Rignot et al. (2016) study. The effects of the fjord and the glacier geometries and of the ocean forcing mechanisms (e.g. tides, coastal trapped waves, buoyancy-driven circulation) are encapsulated in these coefficients. The more vigorous the fjord's circulation and the exchange with the coastal shelf, the higher these coefficients should be. In this study, melting rates presented in metres and cubic metres (m and m³) per unit of time are based on an ice density of 920 kg m⁻³.

3.5 Glacier retreat

To evaluate our glacier melt estimates and to predict the future retreat of Milne Glacier due to submarine melting, we focus on the grounding line position at the centre of the fjord. A thinning of the glacier at the grounding line (δh, negative for thinning, in m a⁻¹) will result in a retreat (δx, negative for retreat, in m a⁻¹) (Thomas and Bentley, 1978; Rignot, 1998; Wood et al., 2021):

$$\begin{array}{}\text{(5)}& \mathit{\delta }h=\mathit{\delta }x\left[\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\mathit{\alpha }}_{\mathrm{b}}-{\mathit{\alpha }}_{\mathrm{s}}\right],\end{array}$$

where ρ_w and ρ_i are the density of seawater (1030 g kg⁻¹) and glacier ice (920 g kg⁻¹). α_b and α_s are the slopes at the surface (fixed at −0.02) and at the base (bed) of the glacier at the grounding line (0.02 in 2020, varies along glacier). Both α_b and α_s are negative when sloping upward up-glacier (Fig. 1c). The thinning at the grounding line (δh, negative for thinning, in m a⁻¹) is a function of the submarine melting at depth (undercutting, M_GL, positive), the velocity of the glacier (u_G), the slope of the glacier face (4:1 in this case), and the surface mass balance (δh_s, negative for ablation, in m a⁻¹):

$$\begin{array}{}\text{(6)}& \mathit{\delta }h=(u_{\mathrm{G}}-M_{\mathrm{GL}})/\mathrm{4}+\mathit{\delta }{h}_{\mathrm{s}}.\end{array}$$

This leads to an expression for the retreat of the glacier (δx) as a function of submarine melting at the grounding line (M_GL):

$$\begin{array}{}\text{(7)}& \mathit{\delta }x=\left[{\displaystyle \frac{(u_{\mathrm{G}}-M_{\mathrm{GL}})}{\mathrm{4}}}+\mathit{\delta }{h}_{\mathrm{s}}\right]/\left[\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}\right){\mathit{\alpha }}_{\mathrm{b}}-{\mathit{\alpha }}_{\mathrm{s}}\right].\end{array}$$

In this study, M_GL is the average melting at 187 m depth (closest grid cell above the grounding line). For a zero surface mass balance, if the undercutting is higher than the glacier velocity, the glacier will thin and the grounding line will retreat, and vice versa.

4 Results: the impact of ice structures

4.1 Water properties

The impact of removing ice structures on the hydrography varies greatly with depth (Fig. 3). To characterize these differences, the water column is divided into four layers: (1) the epishelf lake (0–10 m), (2) the subglacial discharge accumulation layer (10–40 m), (3) an intermediate layer affected by the ice shelf and glacier tongue (40–90 m), and (4) a deep water layer between the intermediate layer and the sill (90–220 m).

4.1.1 0–10 m: epishelf lake

The ∼10 m deep epishelf lake at the top of the water column is an obvious layer with fresh (<10 g kg⁻¹) and warm ($>-\mathrm{1}$ °C) anomalies (Θ−Θ_off, S−S_off; Figs. 3 and 4e). The 2020 MIS calving event results in a shorter basal channel, allowing a doubling of the outflow from the epishelf lake (Table 2). This results in a 0.12 m thinning of the epishelf lake accompanied by a 0.5 g kg⁻¹ salinity increase. This small difference is because the increased exchange only occurs below S_A=5 g kg⁻¹. When the ice shelf is completely removed, the epishelf lake disappears: salinity in this layer increases abruptly to seawater values (>29 g kg⁻¹), and temperature decreases below −1.3 °C.

4.1.2 10–40 m: subglacial discharge accumulation layer

For the two runs with an ice shelf at the mouth of the fjord (pre2020 and now2024), the subglacial discharge plume in Milne Fiord generally settles between 10 and 40 m. This results in warm and fresh anomalies within this layer (Figs. 3 and 4a, i). The influence of subglacial discharge (0 °C, 0 g kg⁻¹) is noticeable on Θ-S diagrams (Fig. 3e, f) by water properties parallel to the freezing line below 31 g kg⁻¹. Within this layer, the integrated amount of subglacial discharge and submarine melting is 2.9 m for the pre2020 run and 2.4 m for the now2024 run. This decrease in f_SD+SM is attributed to enhanced volume exchange of this layer through the basal channel following the 2020 calving event (120 m³ s⁻¹ after vs. 90 m³ s⁻¹ before). The decreased f_SD+SM leads to saltier (0.5 g kg⁻¹) water in this layer. The complete removal of the ice shelf (nois, nogt) leads to unrestricted exchange with the coastal shelf resulting in water properties and stratification near offshore values (Fig. 3). The unrestricted exchange leads to minimal amounts of subglacial discharge and submarine melting waters remaining inside the fjord (<20 cm) and an associated salinity increase (>1 g kg⁻¹). The temperature also decreases by more than 0.5 °C following the calving event, indicating that subglacial discharge and deep water entrained in the plume (both warm anomalies), not submarine melting (negative temperature anomaly), dominate the temperature anomaly signal for the pre2020 and now2024 cases.

4.1.3 40–90 m: intermediate layer

Below the subglacial plume settling depth, from 40 to 90 m (maximum MIS thickness before the MIS calving event), the pre2020 run shows a temperature anomaly of −0.08 °C. The temperature in this layer increases as ice structures are removed (Fig. 4a–d), indicating that trapped submarine meltwater is responsible for this cold anomaly. The amount of submarine meltwater at the centre of the fjord is estimated at 4 cm for the pre2020 run. This decreases to 2.5 cm after the 2020 calving event and further decreases to ∼1 cm for the nois and nogt runs, further confirming trapped submarine meltwater is responsible for this anomaly. The removal of the ice shelf and the glacier tongue has the double effect of decreasing the amount of submarine meltwater (less ice available to undergo melting) and increasing the exchange with the coastal shelf, with both acting to decrease the fraction of submarine meltwater, thereby generating the temperature difference between the runs. The positive ∼ 0.2 g kg⁻¹ salinity anomaly (Fig. 4i), on the other hand, cannot be caused by trapped submarine meltwater, which would decrease salinity. Examination of along-fjord isopycnals reveals that the positive salinity anomaly is caused by upwelling. For the pre2020 run, the mean upwelling at the centre of the fjord is 5.9 m. Together with the average vertical salinity gradient of 0.038 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$, this upwelling yields a salinity anomaly of 0.22 g kg⁻¹, which is consistent with the salinity anomaly observed. The same analysis reveals a decrease in the amount of upwelling with 5.3, 4.3, and 4.6 m for the now2024, nois, and nogt runs, respectively, explaining the freshening observed when ice structures are removed. Upwelling should also result in a temperature increase, but this signal is masked by the accumulation of submarine meltwater (cold anomaly). Consistent with the change in temperature, the reduced upwelling is the consequence of enhanced exchange with the coastal shelf when the ice shelf is removed (Table 2).

4.1.4 90–220 m: deep layer

Between 90 m and the sill, cold temperature anomalies from submarine meltwater are present near the Milne Glacier grounding line for all four configurations (Fig. 4a–d). The presence of submarine meltwater is also evidenced by a temperature depression, which moves the water closer toward the freezing temperature, moving parallel to the melting line on T–S diagrams (Fig. 3e, f). Moving down-fjord, submarine meltwater is diluted with offshore water, increasing the temperature of this layer and allowing the upwelling signal (positive temperature anomaly) to prevail. For example, the submarine meltwater content in this layer decreases from 37 cm 2 km away from the grounding line to 11 cm 18 km away from the grounding line. As the ice structures are removed, the amount of upwelling decreases from 3.5 m (pre2020) to 3.4 m (now2024) to 2.1 m (nois) to 1.7 m (nogt). This decrease results in a cooling (Fig. 4f–h) and freshening (Fig. 4n–p) of the 90–220 m layer.

Table 2 Depth of the 5 g kg⁻¹ halocline, Volume exchange (Q_ex), isopycnal upwelling (Δz_ρ), subglacial discharge fraction (f_SD), and submarine meltwater fraction (f_SM) for the four layers and four ice structure configurations. The number following ± indicates 1 standard deviation from the 3-year average.

4.1.5 Comparison to observations and variability in the water properties

Observations from field campaigns carried out in July 2012 to 2019, 2022, and 2023 generally agree well with model values of f_SM, f_SD+SM, and Δz_ρ (90–220 m layer) (Fig. S2). On the other hand, observed values for Δz_ρ (40–90 m layer) are generally lower than modelled values. This difference is in part because the observed subglacial discharge settling depth is ∼15 m lower than simulated. The most significant cause of this discrepancy, however, is likely the high temporal and spatial variability in the observations.

Overall, our modelling reveals that hydrographic changes in Milne Fiord following the 2020 calving event are small (Fig. 3) and well within the range of observed and modelled variability (Figs. S2 and S3). Without the exact same boundary conditions, these differences would be impossible to detect. Accordingly, observations after the MIS calving event do not reveal significant hydrographic changes compared to before (Fig. S3), except at the very surface (the thinning of epishelf lake).

Figure 3 Water properties from observations (left column) and model simulations (right column). (a) July 2013 temperature profiles from offshore (dark-blue line), at the centre of Milne Fiord (cyan line), and near the grounding line (dotted cyan line). (b) Modelled 3-year average temperature profiles at the centre of Milne Fiord (solid lines) and near the grounding line (dotted lines) and the average offshore temperature profile (black line). (c) Same as panel (a) but for salinity. (d) Same as panel (b) but for salinity. (e) Temperature–salinity plot from the three profiles in panels (a) and (c). The dot-dashed line is the melting line; water involved in melting ice cools down parallel to this line. The negative slope is the freezing point at the surface. The horizontal grey lines denote the four layers identified in Sect. 4.1. (f) Same as panel (e) but for profiles in panels (b) and (d). (g) Density stratification for the three profiles in panels (a) and (c). (h) Same as panel (g) but for profiles in panels (b) and (d).

Figure 4 Along-fjord temperature and salinity anomalies (relative to offshore) and differences (relative to the pre2020 run). (a–d) Time-averaged temperature anomalies (Θ-Θ_off) for the four ice configurations. (e) Time-averaged temperature for the pre2020 run. (f–h) Time-averaged temperature difference between the now2024, nois, and nogt runs and the pre2020 run. (i–l) Time-averaged salinity anomalies (S-S_off) for the four ice configurations. (m) Time-averaged salinity for the pre2020 run. (n–p) Time-averaged salinity difference between the now2024, nois, and nogt runs and the pre2020 run.

4.2 Melting of the ice structures

4.2.1 Milne Ice Shelf (MIS)

The model shows that the MIS experiences both basal melting and freezing before and after the 2020 calving event (Fig. 5a, b). Melting/freezing rates are highly variable in both space (Fig. 5a, b) and time (Fig. 6a); they vary from 15×10³ m³ d⁻¹ of freezing (equivalent to an area-average thickening of 0.06 m a⁻¹) to 75×10³ m³ d⁻¹ of melting (equivalent to 0.46 m a⁻¹ for the pre2020 run and to 0.80 m a⁻¹ for the now2024 run). Generally, melting occurs over areas where the ice draft is greater than ∼40 m, while freezing occurs in shallower areas. This pattern is consistent with water temperature which tends to be near the freezing point (within 0.1 °C) from 10 to 40 m (Fig. 3).

While depth predominantly controls whether melting or freezing occurs because of the temperature profile, the magnitude of the currents under the MIS (|U|) is responsible for periods of enhanced melt (M) (Fig. 6b; R² of 0.89 between spatially averaged M and |U|). Consequently, the highest time-averaged melting rates (0.4 m a⁻¹) are found under the thicker seaward areas of the ice shelf where currents are stronger. Conversely, the maximum freezing rates are located on the west side of the ice shelf in thinner areas where water is trapped and is therefore prone to enhanced vertical heat loss through the ice shelf as lateral exchange is not possible. The high temporal variability in melt rates under the MIS driven by currents is consistent with the unsteady circulation described in Bonneau et al. (2024b).

The 2020 calving event results in a marked increase in melting along the new ice shelf front and around the thicker ice on the west side at X=27 km where Glacier 2 feeds the MIS. However, it does not impact areas of freezing (ice draft <40 m; Fig. 5d). This increase in melting is due to stronger currents now reaching the remaining ice shelf (Fig. 6b) and the warming of the 40–90 m layer inside the fjord (Figs. 4f, 6c), suggesting a weak positive feedback where calving leads to more basal melting, thereby further weakening the ice shelf.

The temporal and spatial average melting rate over the ice shelf is 0.01 m a⁻¹ for the pre2020 run and 0.005 m a⁻¹ for the now2024 run. These values are negligible compared to the temporal and spatial average surface mass balance of 0.30 m a⁻¹ from 2008 to 2017 (White, 2019). However, similarly to the basal mass balance, the surface mass balance of the MIS is spatially variable and has a strong along-fjord gradient, with average ablations of 0.16 m a⁻¹ over the down-fjord half of the ice shelf (X>32 km) and 0.55 m a⁻¹ over the up-fjord section (White, 2019). Therefore, basal melting can be an important process where the ice shelf is thick (>50 m). On the down-fjord section of the ice shelf, it can be the dominating thinning mechanism.

The submarine melting/freezing rate estimates from the pre2020 run have a similar spatial pattern to estimates using a simplified two-equation parameterization and hydrographic observations (Hamilton, 2016). However, melt rates in this realistic numerical simulation are, on average, 10 times less than those estimated by Hamilton (2016) and include negative values (not permitted by the Hamilton (2016) model). By combining the surface mass balance gradient from White (2019) with the basal melt rate estimates from the pre2020 run, we obtain an average ice shelf thinning rate of 31 cm a⁻¹, which closely matches the 29 cm a⁻¹ obtained by Mortimer et al. (2012) for the 1981–2009 period (Fig. 5e). This correspondence supports the low rates of submarine melting obtained in this study and suggests that Hamilton (2016) employed drag and transfer coefficients values that were likely too high and overestimated |U|. The currents they employed were 2–8 times higher than those modelled here. Furthermore, the spatial distribution of thinning over the MIS obtained by Mortimer et al. (2012) aligns well with our results (Fig. 5e), providing additional confidence in the modelled basal melt rates.

Figure 5 Time-averaged basal melting and freezing of the Milne Ice Shelf (MIS) and the Milne Glacier tongue (MGT) from November 2011 to November 2014. (a) Average melting (positive) and freezing (negative) rates for the pre2020 run. Colour scale in panel (c). (b) Same as panel (a) but for the now2024 run. (c) Same as panel (a) but for the nois run. (d) Melt rate difference between the now2024 and pre2020 runs for the MIS and between the nois and pre2020 runs for the MGT. The colour scale is shown in panel (c). (e) Thinning rate for the MIS and MGT using the surface mass balance from White (2019) and submarine melting from the pre2020 run in panel (a). Note the different colour scale for panel (e) (up to 0.6 m a⁻¹).

Figure 6 Time series of total melt rates (M), spatially averaged current speed (|U|), and spatially averaged thermal forcing ($\text{TF}=\mathrm{\Theta }-{\mathrm{\Theta }}_{\mathrm{f}}$). For the MIS and MGT, |U| and TF are from the cells directly under the ice. For the glacier face, |U| and TF are from an across-fjord cross-section 1 km down-fjord from the glacier. Panels (a)–(c) are for the Milne Ice Shelf (MIS), panels (d)–(f) are for the Milne Glacier tongue (MGT), and panels (g)–(i) are for the Milne Glacier face (MG). The grey shading denotes the periods with subglacial discharge.

4.2.2 Milne Glacier tongue (MGT)

The basal melting of the MGT shows a similar spatial pattern to the MIS, with freezing above an ice draft of 40 m and melting below (Fig. 5a). On the other hand, contrary to the MIS, the melting of the MGT has a clear seasonal pattern, driven by velocity and temperature increase during and following periods of subglacial discharge (Fig. 6d–f). Nonetheless, the melting/freezing rates are very small (−0.1 m a⁻¹ $<M<\mathrm{0.1}$ m a⁻¹). The spatially averaged rate for the entire simulation is equivalent to 0.02 m a⁻¹ of freezing for the pre2020 run. This value is negligible compared to the surface mass balance (melting rate of ∼0.55 m a⁻¹; Fig. 5a, e; White, 2019). The negligible basal melting rates obtained over most of the MGT are in agreement with remote sensing estimates from 2011 to 2015 (Hamilton, 2016). However, these melting rates at the bottom of the keel are perhaps too low, likely because of the very low currents there (Fig. 6e) and because ShelfIce does not take into account the ice slopes, which are known to enhance submarine melting (Rosevear et al., 2025).

The partial removal of the ice shelf (now2024) has no significant effect on the basal melting rates of the MGT, current magnitude, and thermal forcing (Fig. 6d–f). The complete removal of the ice shelf (nois) results in an increase in current magnitude by 20 % and a decrease in thermal forcing by 23 % because of stronger exchange with the coastal shelf (Sect. 4.1). The net effect is a 70 % increase of the freezing rates, which overall are still negligible (3 cm a⁻¹).

4.2.3 Glacier face

The melting of the glacier face in the pre2020 run was already discussed in Bonneau et al. (2024a). The main takeaways are that, although melting increases by ∼15 % in summer due to subglacial discharge, most of the melting (∼ 85 %) occurs outside of the subglacial discharge plumes and does not exhibit a strong seasonality. The average melting rate of the glacier face is 33 m a⁻¹, with rates of ∼100 m a⁻¹ at the depth of the grounding line (Fig. 7a).

With a surface slope of −0.02, a basal slope of 0.02 (Fig. 1c), and a grounding line retreat of 53 m a⁻¹ (centreline; Antropova et al., 2024), Eq. (5) yields a thinning rate δ_h of 0.93 m a⁻¹. Combined with ablation stake measurements closest to the grounding line showing a negative surface mass balance (δh_s) of $\sim -\mathrm{60}$ cm a⁻¹(White, 2019) and an undercutting rate (M_GL) of 87.2 m a⁻¹, Eq. (6) yields an average glacier velocity of 85.9 m a⁻¹, which compares well with the observed surface velocities (20 to 160 m a⁻¹; Van Wychen et al., 2016; Millan et al., 2017; Wychen et al., 2020; Antropova et al., 2024), supporting the submarine melt rates obtained from this numerical model.

The results from the four model runs with different ice structure configurations show that the presence of the MIS and MGT has a very limited impact on the melting of the Milne Glacier face (Fig. 7). Melting rates from the now2024 run are the same as the melting rates from the pre2020 run (within 0.5 %). The complete removal of the ice shelf (nois run) leads to a 9 % decrease in melting. This decrease is due to the cooling and freshening of the 90–220 m layer (Figs. 3 and 4), which decreases the thermal forcing by ∼0.05 °C. The further loss of the glacier tongue partly offsets this cooling, resulting in an increase in melting rate for the nogt run compared to the nois run, although, in both cases, the melting rate is lower than before the calving event (pre2020 run). A noteworthy change arising from the removal of the glacier tongue is the surfacing of the two subglacial discharge plumes (Fig. 7d). This is due to the lower density stratification in front of the glacier (Fig. 3h) because of the increased exchange with the rest of the fjord (fresher water not trapped behind the glacier tongue; Fig. 3d).

Figure 7 Time-averaged melting rates over the glacier face for the pre2020 (a), now2024 (b), nois (c), and nogt (d) simulations.

Figure 8 Melting rate of the glacier face for the five ocean warming scenarios (a–e; see subpanel title). The spatial and temporal average from November 2011 to November 2013 is in the top-left corner of each panel. (f) Time series of spatially integrated melting over the Milne Glacier face for the different temperature scenarios (see legend).

5 Results: consequences of ocean warming

According to our simulation with different ice structure configurations, the presence of the MIS and MGT has a negligible impact on the melting rate of the glacier face (≤9 %; Fig. 7). Therefore, the results of this study of the impact of ocean warming on the melt rate of Milne Glacier are applicable regardless of the state of the ice shelf and glacier tongue.

The time-averaged TF (Θ−Θ_f) is 0.76, 0.84, 1.2, 1.5, and 3.1 °C for the nogt T03, T09, T16, and T30 runs, respectively. The results show that a TF increase of 0.08 °C (T03) has limited impact (total melt over the glacier face increases by 8 %) but that further warming has a more substantial impact: 38 % increase for T09, 141 % increase for T16, and 420 % increase for T30 (Fig. 8). For the five simulations, melt rates increase by a similar proportion during summer: up to ∼25 % for summer 2012 and up to 10 % for summer 2013, i.e. linear response to subglacial discharge. We note a possible residual effect of subglacial discharge on melting as elevated values of submarine melting starting at the onset of subglacial discharge persist ∼3 months after subglacial discharge stops.

Calculating the submarine melting rate at the grounding line (M_GL) as the average melting at 187 m and comparing it to the parameterization of Rignot et al. (2016) Eq. (4) reveals substantial discrepancies between this parameterization and our results (> factor 2; Fig. S4). Therefore, the coefficients of Eq. (4) were determined using the five ocean warming simulations. The expression obtained for M_GL for Milne Fiord is

$$\begin{array}{}\text{(8)}& M_{\mathrm{GL}}=\mathrm{0.31}{\text{TF}}^{\mathrm{0.90}}.\end{array}$$

The main difference between our parameterization (Eq. 8) and the one from Rignot et al. (2016) (Eq. 4) is the exclusion of the subglacial discharge term (q_sd). This exclusion is justified by our optimization (least-squared method), which yields a better correlation between modelled and parameterized M_GL when A is set to 0. This is explained by the negligible increase in undercutting (M_GL) during summer (Fig. S4); unlike the total amount of melting over the whole glacier face (e.g. Figs. 6g, 8f), M_GL does not exhibit a seasonal cycle. We attribute this major difference in the parameterization proposed by Rignot et al. (2016) to the different numerical setup: Rignot et al. (2016) (using the model from Xu et al., 2013) modelled a 150 m wide glacier with two subglacial discharge outflows, while we modelled a 4000 m wide glacier with two outflows. Therefore, the proportion of the glacier face affected by subglacial discharge plumes at depth is ∼27 times lower in our case, resulting in a negligible dependence of M_GL on subglacial discharge when averaged over the width of the glacier.

Figure 9 Undercutting from submarine melting (M_GL) as a function of thermal forcing (TF).

6 Discussion

6.1 Ice cover

This is the first study that we are aware of that examines the impact of ice shelf or glacier tongue calving on water properties and submarine melting rates in a glacial fjord. Our results show that the removal of ice structures in Milne Fiord has a limited impact on water properties in the fjord. The exception is the surface water (above 15 m depth), which becomes considerably saltier following the complete removal of the ice shelf and drainage of the epishelf lake. This drainage will lead to sea ice formation inside the fjord (instead of the formation of lake ice) resulting in occasional ice-free periods in summer, as is seen in neighbouring fjords (e.g. Veillette et al., 2008). For example, since the calving of Ward Hunt Ice Shelf (80 km east; Fig. 1a) in 2001 and the subsequent drainage of its epishelf lake, open water was visible in Disraeli Fiord in 17 out of the 20 following summers. In stark contrast, open water in Milne Fiord has occurred only once during the same period.

6.2 Ice shelf and glacier tongue melting

The low simulated submarine melting rates under the MIS and MGT (space and time average <0.02 m a⁻¹) are consistent with studies which identify atmospheric warming as the main cause for the loss of more than 75 % of the ice shelf extent along Ellesmere Island over the last 60 years (White and Copland, 2019; Copland et al., 2017, 2007; Mueller et al., 2003). Warmer air temperatures are also thought to have played a role in the recent retreat of Hunt Fjord Ice Shelf in northern Greenland (Ochwat et al., 2023). However, the comparison of the spatially varying basal melt rates obtained in this study with surface mass balance estimates reveals that submarine melting can be the dominating thinning mechanism where ice is thick (>50 m) near the mouth of the fjord. At this location, the surface mass balance is also less negative because of lower summer temperature, higher precipitation, and summer fog (White, 2019). We note observations of surface mass balance even closer to equilibrium near the seaward edge of Ward Hunt Ice Shelf (0.07 m a⁻¹ of melting; Braun et al., 2004). According to ORAS5 reanalysis (Copernicus Climate Change Services, 2021), the thermal forcing of the ocean above 100 m in front of Milne Fiord did not change substantially between 1958 and 2019 (Fig. S5). Considering the seaward portion of Ayles Ice Shelf (15 km to the east of Milne Fiord, calved in 2005) was ∼44 m thick (Copland et al., 2007) and that the seaward portion of Ward Hunt is ∼50 m thick (Braun et al., 2004), it appears that the Ellesmere Island Ice Shelf was sufficiently thick for submarine melting to be an important component of the mass balance. Thus, we hypothesize that submarine melting has substantially contributed to the disintegration of the Ellesmere Island Ice Shelf by weakening its seaward edge. Since the ocean along the north coast of Greenland has a similar water column structure than along Ellesmere Island (Johnson et al., 2011; Schaffer et al., 2017), it is not surprising that ocean-induced melting is thought to have played a role in the more recent retreat of ice shelves in northern Greenland (Hill et al., 2017; Millan et al., 2023).

The limited sensitivity of submarine melting rates to past (now2024) and future (nois) calving events in Milne Fiord is similar to the results obtained in a recent numerical study of Pine Island Ice Shelf (Bradley et al., 2022), where it was estimated that the retreat of the ice front by up to 50 km increased the spatial and temporal average ice shelf melt rate by less than 10 %. Similar results were also found for Larsen C Ice Shelf for which numerical simulations indicate that the calving of iceberg A-68 (5800 km²) did not substantially impact the spatial and temporal average ice shelf melt rate (Poinelli et al., 2023). However, similarly to our results from the now2024 run, Poinelli et al. (2023) showed a local increase in melting along the new ice shelf front (Fig. 4b, c). Moreover, even though the spatial and temporal average melt rates did not vary significantly following the calving of A-68, Poinelli et al. (2023) calculated a doubling of melt rates around an important ice shelf pinning point and noted that localized change in submarine melt rates can have an overall destabilizing impact. As the MIS is only grounded on its sides, the pinning points are lateral and the most important one is likely where Glacier 2 (west side, X=27 km; Fig. 1b) feeds the ice shelf. Our results show this is the location where submarine melt rates increase the second most following the 2020 calving event (∼75 % increase), suggesting that enhanced local submarine melting may lead to further structural weakening.

6.3 Glacier face melting

Using two-dimensional numerical simulations of Petermann Glacier, Cai et al. (2017) argued that the melt rates at the grounding line should increase following the removal of the ice shelf, as this would lead to a steeper under-ice slope. We cannot confirm that this would be the case for Milne Glacier, as the slope near the grounding line was kept constant in our numerical experiments. However, as our results show, water reaching the grounding line has similar temperature (within 0.15 °C) and salinity (within 0.1 g kg⁻¹) regardless of the ice structure configurations. Therefore, it is likely that local melt rates will increase if the slope becomes steeper.

Comparing the two northern Greenland fjords where the Petermann and Ryder glaciers terminate, Jakobsson et al. (2020) found warmer water in Petermann fjord because of a deeper sill. This difference is perhaps why the grounding line of Ryder is stable (Holmes et al., 2021), while the grounding line of Petermann is retreating (Millan et al., 2022), highlighting the role of submarine melting at depth. Milne Glacier, which is currently retreating and grounded above the sill, likely shares a similar sensitivity to ocean conditions to Petermann Glacier.

Consistent with our results, numerical simulations of Pine Island Ice Shelf (Bradley et al., 2022) and Larsen C Ice Shelf (Poinelli et al., 2023) with constant ice morphology also suggest that major calving events have limited impact on submarine melting along the grounding line, with more pronounced impacts near the seaward edge or along the main ocean intrusion pathways.

6.4 Future submarine-melting-induced retreat of Milne Glacier

Using the M_GL parameterization (Eq. 8), a constant surface mass balance ($\mathit{\delta }{h}_{\mathrm{s}}=-\mathrm{0.60}$), and a constant glacier velocity (u_G=85.9 m a⁻¹) calibrated with the observed grounding line retreat (53 m a⁻¹), we calculate grounding line retreat for three different greenhouse gas emission scenarios (SSP126, SSP245, and SSP585; S1). For each of these scenarios, the ocean temperature (therefore TF) is determined by using the (time- and depth-varying) multimodel mean temperature increase. For each scenario, M_GL is calculated using the time- and depth-varying TF (Eq. 8), then the glacier retreat is calculated for 1 year (Eq. 7). The grounding line position and depth are then updated, and the retreat for the next year is calculated.

Figure 10 Submarine-melting-induced retreat of Milne Glacier according to different carbon emission scenarios (see legend). Time series of (a) thermal forcing, (b) depth of the grounding line, (c) undercutting rate, (d) submarine-melting-induced retreat rate of the grounding line, (e) cumulative retreat of Milne Glacier grounding line.

These estimates suggest that the grounding line of Milne Glacier will experience an accelerated retreat until 2030 and a continual near-constant rate of retreat of ∼200 m a⁻¹ during the following 30 years, regardless of the greenhouse gas emission scenario (Fig. 10). After 2065, the low carbon emission scenario (RCP2.6/SSP126) leads to a stabilization of the glacier, while the moderate and high greenhouse gas emission scenarios result in an uninterrupted retreat. The rate of retreat of 200 m a⁻¹ obtained from 2030 to 2035 is similar to what was observed for Umiamako Isbrae in West Greenland (180 m a⁻¹ from 1989 to 2015; Rignot et al., 2016) and Petermann Glacier in northern Greenland (130–230 m a⁻¹ from 1992–2021; Millan et al., 2022).

Varying the basal slope (α_b) and the surface slope (α_s) by ±50 % and increasing the surface mass balance (δh_s) by a factor of 4 does not significantly alter these results (Fig. S6). However, these estimates do not include glacier dynamics. For instance, the MIS, MGT, and perennial ice likely provided some buttressing to the glacier presently and in the past (Scambos et al., 2004; Mouginot et al., 2015, 2019). The loss of these ice structures could therefore lead to a speed-up of the glacier, resulting in thinning and further grounding line retreat (Joughin et al., 2021; Millan et al., 2022). The MIS, MGT, and perennial ice also protect the glacier from ocean swell and currents (Glasser and Scambos, 2008; Massom et al., 2018). These processes could also lead to faster rates of retreat than what is estimated here. Finally, not accounted for in this study, air temperature will also increase (e.g. Cai et al., 2021), resulting in a more negative surface mass balance and additional subglacial discharge. These two mechanisms will likely lead to additional glacier retreat. Numerical modelling of Milne Glacier would allow the investigation of dynamical ice feedback and enable the validation of the grounding line retreat parameterization.

The calculated and predicted retreat of Milne Glacier is in line with the ongoing retreat of >80 % of marine-terminating glaciers in the Arctic (Kochtitzky and Copland, 2022) and the continual retreat of ice shelves (Hill et al., 2017; Carr et al., 2017; White and Copland, 2019). The retreat of marine-terminating glaciers directly leads to sea level rise (Jakob and Gourmelen, 2023; Greene et al., 2024) but can also transform ecosystems that depend on glacier-induced upwelling for nutrient delivery (Bhatia et al., 2021; White et al., 2025).

7 Conclusion

In this study, we employed a numerical model evaluated against observations to explore ice–ocean interactions in Milne Fiord. To investigate the rapid changes currently happening in this system, we examined the impact of the loss of ice structures and ocean warming on the hydrography and submarine melting rates. The main results are as follows:

1. Although the spatial and temporal average submarine melt rate under Milne Ice Shelf (<0.02 m a⁻¹) is negligible compared to the surface mass balance (0.3 m a⁻¹), submarine melting is the dominant thinning mechanism locally where the ice shelf is >50 m thick near the ocean.

2. The loss of ice structures impacts the hydrography of Milne Fiord differently with depth. Above 15 m, the removal of the ice shelf results in a strong salinity increase (>15 g kg⁻¹). Below 15 m, the impact of the removal of ice structures is limited to a temperature change of 0.15 °C and salinity change of 0.2 g kg⁻¹, resulting in small changes in spatial and temporal average submarine melting/freezing rates. However, two specific locations now experience enhanced melting following the 2024 MIS calving event: the new ice shelf front and the basal surface of Glacier 2. These localized regions, where melt rates show enhanced sensitivity to change in ice structure configuration, suggest a weak positive feedback loop in which ice shelf calving leads to more melting which could promote more ice shelf calving.

3. The melting of the Milne Glacier face responds quasi-linearly to ocean warming. According to the analysed CMIP5 and CMIP6 predictions, water offshore of Milne Fiord will warm by at least ∼0.2 °C by 2040 regardless of the carbon emission scenario. Predictions based solely on changing submarine melting indicate a retreat of the Milne Glacier grounding line by 8 km by 2065. These predictions also suggest a subsequent stabilization for low greenhouse gas emission scenarios and further grounding line retreating for higher emission scenarios.

Situated at the centre of Tuvaijuittuq's coastline, Milne Fiord is rapidly changing. This study highlights the vulnerability of this newly designated Marine Protected Area: the ice shelf and glacier tongue will disappear, non-perennial sea ice will replace the permanent lake ice, and glacier retreat will likely accelerate in the next decade.