Active ice sheet conservation cannot stop the retreat of Sermeq Kujalleq glacier, Greenland

Introduction

Ice loss from Greenland ice sheet is raising global mean sea level by ~0.8 mm yr^-1 since 2002^1,2,3. The mass loss of the Greenland ice sheet has been attributed to two separate processes over the 1972–2018 period⁴ about 2/3 controlled by glacier dynamics (ice mass flux into the ocean) through its marine-terminating glaciers, and about 1/3 controlled by surface mass balance (SMB; surface accumulation minus runoff and other forms of ablation) from enhanced surface melt. Furthermore, although surface melting has increased since 2000s, ice dynamics is expected to remain the primary driver of Greenland ice sheet mass loss over 21st century ⁵. Oceanic forcing impacts ice dynamics through submarine melting and calving, while atmospheric forcing directly impacts surface melt and ice calving by filling crevasses with melt water. The two forcings are coupled, because ocean heat can warm the coastal atmosphere, on the other hand, ocean-driven melting can be heavily influenced by subglacial discharge which may be increased by surface runoff even in the absence of ocean warming. Both forcings also impact the buttressing provided by the ice mélange that fills the fjords in front of actively calving outlet glaciers^6,7.

Air temperatures respond rapidly to changes in radiative forcing, hence the acceleration of surface runoff can be slowed by following more stringent greenhouse mitigation pathways. But the rates can only be slowed relative to the present by reducing forcing below present levels, which is not expected for existing international agreements on limiting greenhouse gas emissions, before midcentury at best⁸. However, radiative forcing can be lowered by reducing incoming solar radiation (solar geoengineering), which is simulated to lower both surface runoff^9,10,11 and dynamic ice losses¹² relative to pure greenhouse gas scenarios, but ice losses and sea level contributions will still rise relative to present-day values. Since more than 90% of the excess heat accumulated since the industrial revolution is stored in the oceans, ocean heat content cannot be reduced rapidly.

Ice shelf basal melt driven by ocean thermal forcing reduces the buttressing as contact area between shelf and the topography it sits in are reduced and may trigger grounding line retreat. If a glacier retreats inland into deeper water, it will be subject to increased ocean thermal forcing as deep waters are warmer than at the surface in the polar regions due to their higher salinities, and the pressure-induced freezing point depression of ice. Initial retreat onto a reverse-sloping bed can lead to Marine Ice Sheet Instability (MISI)¹³. The retreat will continue until the terminus reaches a prograde slope unless large enough additional buttressing is provided, for instance, by lateral friction in a fjord, at the deep bed of the glacier, or sea floor high spots. High basal melt rates increase the formation of channels and crevasses, which affect the topography, roughness, and stability of ice shelves^14,15, and atmospheric warming increases surface meltwater supply that percolates through the surface crevasses. Both sets of crevasses weaken the ice and influence calving rates¹⁶.

One way to provide buttressing might be to lower the ocean thermal forcing, allowing the glacier to thicken, and potentially increase friction, such as by thickening any floating ice shelves that might reground on sea floor high spots. Wolovick and Moore¹⁷ used a suite of coupled quasi-2-D ice–ocean simulations to explore whether targeted geoengineering using either a continuous artificial sill or isolated artificial pinning points could counter a collapse of Thwaites Glacier. This was developed further in an engineering feasibility study¹⁸ and an elementary cost-benefit analysis¹⁹ for the Amundsen Sea outlets. Moreover, the local acceptability and sustainability aspects that determine the societal feasibility of the seabed anchored curtain were studied in Ilulissat Icefjord (also known as Kangia). Blocking warm water inflow to Greenland outlet glaciers has been suggested²⁰, but without any simulations of the impacts on outlet glaciers.

In this study, we explore whether blocking the warm water by a barrier at the fjord mouth could mitigate retreat and mass loss of the Sermeq Kujalleq (Jakobshavn Isbræ in Danish) glacier in western Greenland. Multiple forms of hypothetical subsea barriers have been proposed as a means to block deep warm water from reaching vulnerable glaciers, including earthen dams¹⁷, rigid steel plates²⁰, or flexible buoyant curtains¹⁸. We do not explore the details of barrier design in this study; our model simply assumes that an engineered barrier blocks exchange between the fjord and the external ocean below a specified blocking depth, and we explore the impact of this blockage on fjord properties and glacier dynamics. Sermeq Kujalleq (Fig. 1) has been the fastest and largest (by ice discharge volume) glacier in Greenland over recent decades except for a few brief intervals^4,21. In common with most outlet glaciers in Greenland, Sermeq Kujalleq terminates in a fjord (Ilulissat Isfjord). The fjord is 8 km wide and about 800 m deep⁷ but has a sill reducing the depth of the entrance to about 250 m (Fig. 1). The fjord basin water is mixed and renewed by ambient water in Disko Bay flowing over the sill, subglacial discharge plume from beneath the glacier terminus and glacier submarine melt^7,22. The circulation is also substantially modified by icebergs²³.

**Fig. 1: Surface velocity and bed elevation of our model domain in the Greenland Ice Sheet.**

The terminus position of Sermeq Kujalleq was relatively stable between 1950 and 1996 with seasonal fluctuation of about 2.5 km²⁴. Sermeq Kujalleq thickened substantially from 1991 to 1997, and rapidly thinned after 1997. The glacier flowed into a 15 km floating ice tongue before 1997, which then retreated and finally disintegrated in the early 2000s²⁵. Holland et al.⁶ suggested that warm water originating from the Irminger Sea enhanced basal melting of the ice tongue and triggered its break up in 1997. The glacier has accelerated since 2004 and reached maximum velocities of 17 km yr⁻¹ in 2012 and 2013²⁶. The acceleration may have been facilitated by anomalously high surface melting^27,28. The glacier then decelerated after 2012 but significant seasonal variability in ice velocity remained²⁹. The thinning rates slowed from 2014, and the glacier thickened and re-advanced between 2016–2019³⁰. This deceleration coincided with a period of cooler ocean waters in Disko Bay which connect to Ilulissat Icefjord, when ocean temperatures in Disko Bay’s upper 250 m cooled to those of the mid 1980s³⁰. Joughin et al.³¹ suggested that the ocean’s dominant influence on Sermeq Kujalleq is through its effect on winter mélange rigidity, rather than summer submarine melting.

The observed historical variability of Sermeq Kujalleq correlated with ocean temperature motivates us to explore the lowering the oceanic thermal forcing by simulating a raised sill at the fjord mouth under a range of climate warming scenarios (SSP2-4.5 and SSP5-8.5) and with historical ocean and air temperatures. Here we simulate the effects of an artificial barrier on the top of the sill at the fjord mouth of various heights on the fjord water properties and glacier basal melt rate, to examine if such an intervention might mitigate retreat and mass loss of Sermeq Kujalleq.

Result

Fjord water temperature and melt rate evolution under different scenarios

We first describe the fjord water temperature evolution in the absence of any barrier. The ocean water near the fjord mouth at 150–250 m depths was relatively cool during the years 2015–2019 (Fig. 2), which led to the modeled fjord water temperature being noticeably lower at the same time (Fig. 3). There is no significant trend for fjord water temperature under SSP2-4.5 for the remainder of the century. However, under SSP5-8.5 there is clear warming trend in the fjord water, especially during 2070–2100, though with seasonal fluctuation and inter-annual variability. Surface runoff greatly increases after 2070 under SSP5-8.5, leading to more subglacial discharge at the grounding line (Fig. 4), which enhances the mixing of plume and fjord water, and results in more melt of the glacier terminus (Figs. 3, 5).

**Fig. 2: Reanalysis and projected monthly ocean temperatures, 1993–2100.**

Fig. 3: Modeled annual averaged fjord water temperature near the ice front in the experiments with barrier height of 100 m, 50 m and no barrier under SSP2-4.5 and SSP5-8.5 for the period 2009–2100.

**Fig. 4: Subglacial discharge and surface runoff under SSP2-4.5 and SSP5-8.5.**

**Fig. 5: Modelled mean melt rate at the ice front in summer and rest of the year with different barrier heights under SSP2-4.5 and SSP5-8.5, 2009–2100.**

Higher barriers yield cooler fjord water (Fig. 3) and reduced melt rates at the ice front (Fig. 5). We report results with barrier heights of 50 and 100 m. We did explore higher barriers but we do not have confidence in the fjord model simulation of the two-layer exchange flow when the fjord mouth is reduced to less half of its the present-day opening, so we only discuss the 50 and 100 m results compared with no barrier. The 0~250 m depth-averaged fjord water temperature in July 2100 under SSP2-4.5 and SSP5-8.5 drops by ~0.4 °C and ~0.6 °C respectively for every 50 m increase in barrier height. Melt rate at the ice front in summer (June-July-August, JJA) is much larger than in other seasons (Fig. 5), due to the presence of subglacial discharge and higher water temperature.

Glacier evolution under different scenarios

The glacier is subjected to both atmospheric forcing (surface runoff, SMB) and ocean forcings (basal melt rate, ice mélange buttressing effect in calving). We use the moderate SSP2-4.5 forcing (Ex1; Table 1) and the SSP5-8.5 high end forcing (Ex2) outputs from MAR¹¹ and CNRM-CM6−1-HR³². We also explore how ocean and atmospheric forcing compare by keeping ocean forcing under SSP5-8.5 but replacing elements of the atmospheric forcing (Table 1). The three extra forcings are: (i) Ex3: repeated surface runoff in 2010–2019 and SMB under SSP5-8.5, (ii) Ex4: repeated surface runoff and SMB in 2010–2019, (iii) Ex5: repeated surface runoff and SMB in 1990−1999. We also do another experiment (Ex6) using repeated ocean forcing in 1993–2002 (cool ocean) and repeated surface runoff and SMB in 1990–1999.

Table 1 Experiments under future scenarios

Full size table

Figure 6 shows that modeled ice front position in July retreats in response to different atmospheric and ocean forcings in the experiments (Table 1). More surface runoff boosting the development of surface crevasses, increases ice front calving frequency, intensifies glacier retreat, and also yields more subglacial discharge thus enhancing the basal melt rate at the ice front. The fjord water temperature increases with warmer ocean water near the fjord mouth, causing more melt at the ice front. Meanwhile, the ice mélange inside the fjord gives less buttressing to the glacier terminus as fjord water warms (Eq. (5)). In general, the modeled glacier terminus retreats more with warmer atmospheric forcing (more surface runoff; less SMB) or warmer ocean temperature forcing.

Fig. 6: Modeled cross-sectional profiles along the center flowline near the terminus in July over the period 2020–2100 (color-coded by the top scale) using α₁ = 0.6, α₂ = 0.15, β = 34.

With no barrier at the sill, the glacier retreats 19.91 km under SSP2-4.5, and 51.39 km under SSP5-8.5 from 2020 to 2100. Modeled ice front retreats by 22.11 km and 18.94 km with barrier height of 50 m and 100 m under SSP2-4.5, and 51.24 km and 50.66 km with barrier height of 50 m and 100 m from 2020 to 2100 under SSP5-8.5. Therefore, barriers with height up to 100 m make little difference to the glacier terminus retreat under these future climate scenarios.

In Ex3 (3rd column of Fig. 6) surface runoff under SSP5-8.5 is replaced with that from 2010–2019. Comparing Ex3 with Ex2, annual mean subglacial discharge during 2020–2100 has a reduction of 19%, and JJA-mean subglacial discharge during 2020–2100 has a reduction of 35%. The subglacial discharge since 2060 under SSP5-8.5 is significantly higher than that in 2010–2019 (Fig. 4). The subglacial discharge in JJA months during 2060–2100

in Ex3 has a reduction of 48% compared with that under SSP5-8.5. The glacier retreats noticeably less after 2060 in Ex3 than under SSP5-8.5 in Ex2: at 2100 by 40 km, 41 km, and 30 km, with barrier heights of 0 m, 50 m, 100 m. Reduced runoff at the ice front (Fig. 4) since 2060 leads to less ice calving and also contributes to the deduction of glacier retreat in Ex3.

In Ex4 (4th column of Fig. 6) we replace both surface runoff and SMB under SSP5-8.5 for 2020–2100 with repeated data for 2010–2019 (an increase in SMB of 41%), which reduces glacier retreat to 23 km, 20.35 km and 19.77 km, with barrier heights of 0 m, 50 m and 100 m. As with the other experiments, the largest reduction in glacier retreat occurs during 2070–2100. Comparing Ex4 result with Ex3, indicates the importance of maintaining positive SMB to stabilize the glaciers.

Comparing Ex4 with Ex1 without barrier in which they have the similar SMB and runoff in a cold atmospheric scenario but different ocean forcing, we can see the glacier is more unstable and retreats at earlier years on the retrograde slope when subject to warmer ocean thermal forcing. The averaged relative temperature of the inflow ocean water in JJA with respect to freezing point at 250 m depth is 4.6 °C in Ex4 under SSP5-8.5 in the period 2070–2100, while it is 3.4 °C in Ex1 under SSP2-4.5 (Supplementary Fig. 2). This leads to more retreat in Ex4 since 2070 with ocean thermal forcing under SSP5-8.5, although the glacier sits on a prograde slope. Comparing results in Ex5 with Ex6 without barrier leads to a similar finding.

In Ex5 (5th column of Fig. 6) we replace surface runoff and SMB under SSP5-8.5 in 2020–2100 with repeated data from 1990–1999, in which the subglacial discharge is 0.4~0.88 m²·s^-1 (48.5% lower) and surface runoff at the ice front is 0.82~1.27 m w.e. month^-1 (27% lower). The glacier retreat is again further reduced to 19.49 km, 19.49 km, 18.48 km with barrier heights of 0 m, 50 m, 100 m by 2100.

In Ex6 (6th column of Fig. 6) we replace ocean temperature under SSP5-8.5 in 2020–2100 with repeated data from 1993–2002, but keep the same atmospheric forcing as in Ex5. The glacier retreat is 15.54 km, 12.47 km, 12.91 km with barrier height of 0 m, 50 m, 100 m at 2100. Cooler ocean temperature reduces retreat by a further 20.27%, 36.02%, 30.14% with barrier height of 0 m, 50 m, 100 m relative to Ex5 under the same atmospheric forcing.

The plane view of glacier terminus position evolution in different experiments is shown in Fig. 7. Under SSP2-4.5 forcing the basin retreats from the coastal margins, but under SSP5-8.5 the dramatic retreat along the sub-glacial central trough is much more apparent. This retreat is largely prevented under the conditions in Ex4, 5 and 6. The barrier makes little discernible difference to the pattern of retreat under all experiments. The modeled result using the second best set of parameters are α₁ = 0.5, α₂ = 0.11, β = 36 (Supplementary Figs. 3 and 4) is close to that using the best set of parameters with differences in retreat of typically 1–6 km by 2100.

Fig. 7: Grounding lines contoured decadally as a function of year (color bar) with barrier heights of 100 m (top row), 50 m (middle row) and no barrier (bottom row) using α₁ = 0.6, α₂ = 0.15, β = 34.

Discussion

The mean fjord water during the 2080–2100 is 0.5 °C, 1.1 °C cooler with barrier heights of 50, 100 m than that without barrier under the SSP5-8.5 scenario. We investigated the mechanism behind the sill cooling fjord water temperature. The thermal forcing in the inflowing layer over the sill during the 2080–2100 period is reduced by 23% with a barrier height of 50 m and 48% with a barrier height of 100 m relative to the period 2009–2019 with no barrier. The exchange flow changes are much smaller: increases of 2% with a barrier height of 50 m and a decrease of 2% with a 100 m barrier. Hence, the sill largely acts via cooling water temperatures behind the barrier rather than altering exchange fluxes.

Climate-ice sheet intervention aims to tackle the global-scale impact of climate warming on sea level rise. This might be done by local engineered methods at a handful of high-leverage locations, such as ice streams and outlet glaciers that are subject to MISI. Several marine-terminating Greenland outlet glaciers have their terminus in deep fjords that extend far inland of the coast. A leading cause of these outlet glacier retreats is the warm subsurface water into the fjords from the North Atlantic Ocean^6,33,34 coupled with vigorous circulation driven by subglacial discharge plumes³⁵, which results in undercutting at glacier termini, and contributed to the loss of almost all Greenland ice shelves, which then leads to widespread glacier retreat and ice flow acceleration. This motivates the exploration of an intervention we simulate here.

There are other issues that would need to be considered before any such intervention were done including social acceptability and political feasibility of the method to the Greenland authorities, and by local communities. These communities exploit the rich fishing that has increased in value as the fresh water supply from beneath the marine-terminating glacier has increased, bring more deep nutrients into the photic zone supporting more phytoplankton³⁶. The warming Atlantic waters have also introduced new commercially important species such as Greenland halibut to the Disko Bay region^36,37,38. But any such considerations depend on the question of whether the intervention may feasibly reduce global sea level rise. If the intervention is unlikely to be successful, then the other considerations are moot. Our simulations suggest that an intervention to block warm water intrusion is unlikely to prevent the further retreat of Sermeq Kujalleq. While additional uncertainty remains, and no one modeling study should ever be taken as the final word, these results do suggest that this glacier may have passed the point at which an ocean-focused intervention can save it. Future research into glacial interventions should aim to constrain when different glaciers will pass the point at which interventions will no longer help.

Some uncertainty of the model result comes from the choice of climate forcing data. We only used modeled ocean temperature from CNRM-CM6-1-HR due to its high resolution, after bias-corrected using reanalysis data. Some uncertainty also comes from the use of calving parameterizations, but the comparison between the best and second best set of parameters (Figs. 6 and 7 with Supplementary Figs. 3 and 4) shows qualitatively similar results, but with some difference in the duration of stability on the small prograde bedrock slopes during retreat. This suggests that the results are not particularly sensitive to the selection of the model parameters used. For the treatment of calving in BISICLES, we use a calving parameterization depending on the crevasse penetration depth, which relies on surface runoff. After tuning the parameters, the model did not do a particularly good job of fitting the observed calving front position and motion. But calving is a hard problem for all ice sheet models. Many models struggle to reproduce steady-state calving front positions, let alone the dynamic history of ice front position. Several calving parameterizations have been proposed, including a calving law depending on crevasse depth³⁹, along- and across-flow strain rates (eigencalving)⁴⁰, and the von Mises tensile stress⁴¹. Choi et al.⁴² evaluated four calving laws over nine tidewater glaciers of Greenland, and found von Mises calving law generally produces terminus change closest to observed, and eigencalving the furthest. More simulations with other calving laws, and other mechanisms for coupling calving to ocean conditions, are needed to improve predictions of glacier dynamics and to improve our understanding of the responsiveness (or lack thereof) of glaciers to interventions.

Some uncertainty also comes from the submarine melting parameterization which we calculate at the vertical face of ice cliff in the fjord model. The plume model result depends on entrainment coefficients, drag coefficient, and turbulent transfer coefficients for heat and salt, etc. Almost all tidewater glacier studies with plume model use the standard set of values (e.g^43,45,46,46.). But these coefficients are empirically derived, and have not been validated. For instance, the frictional drag at the ice-ocean interface is parameterized as a quadratic function of the plume speed with a drag coefficient. And the standard drag coefficient value comes from MacAyeal (1985)⁴⁷ for simulations in which he ignored the Coriolis force. To match observations, Jackson et al.⁴⁸ suggests that if the basic physics of plume-melt theory are correct, then critical parameters in the plume parameterization must be adjusted. Davison et al.⁴⁹ varied the turbulence parameters for iceberg melt and showed that larger values of turbulent transfer enhance the impacts of icebergs. In our plume model, the subglacial discharge is assumed to be uniform across the width of water sheet released at the grounding line, hence it does not include entrainment in the across-fjord direction. But for a discrete localized source of subglacial discharge input at the grounding line, the half-conical plume model would be more appropriate based on localized runoff input at the grounding line⁴⁴ using a different entrainment coefficient.

We find that the maximal monthly melt rate modeled by the fjord model in the summer months (JJA) can be approximated reasonably well by a function that has a nonlinear dependence on inflow water temperature with a power exponent of 0.46–0.55, and a nonlinear dependence on subglacial discharge with a power exponent of 0.23–0.24. The power exponent for dependence on inflow water temperature is about half that is found by Jenkins (2011)⁴³, but he used ambient temperature which will be different from the fjord inflow water temperature that we use. The power exponent of dependence on subglacial discharge is slightly smaller than the 1/3 power found by Jenkins (2011)⁴³. The approximation function we find for melt rate is

$$m={e}^{6.9379}cdot {{Q}_{sg}}^{0.2362}cdot {({T}_{a}-{T}_{af})}^{0.4635}$$

(1)

under SSP2-4.5 (R² = 0.87) and

$$m={e}^{6.8622}cdot {{Q}_{sg}}^{0.2279}cdot {({T}_{a}-{T}_{af})}^{0.5538}$$

(2)

under SSP5-8.5 (R² = 0.96), where e is the Euler number, m is the maximal basal melt rate at the vertical face of the ice cliff, Q_sg represents subglacial discharge, T_a is the inflow temperature water at 250 m depth at the fjord mouth, and T_af is the freezing point of water at depth of 250 m, depending on the salinity of the inflow water. Using the above approximation, a reduction of 48% in JJA-mean subglacial discharge during 2060–2100 in Ex3 compared with Ex2 leads to a reduction of 14% in melt rate. A large fraction of glacial melt rate variability can be thought of as a function of subglacial discharge (e.g^35,50.), which comes from surface meltwater. Hence, atmospheric forcing would still be a dominant factor for glacier retreat in the future. However, recent work⁵¹ also suggests that seawater intrusions underneath the floating portion Sermeq Kujalleq might contribute considerably to basal melt rates, so perhaps the use of a plume parameterization itself represents a large uncertainty. In addition, we do not use the melt rate at the vertical face of ice cliff in the ice sheet model. Instead, we use equivalent basal melt rate beneath the very short ice shelf. More simulations that implement melt rate at the vertical face of ice front are required to test if this makes significant differences to the modeled glacier dynamics.

The bedrock slope plays an important role in glacier evolution. The modeled glacier retreats in all the experiments we design under either cooler or warmer climates than at present. The fast retreat of the grounding line in the first few decades is partially because the glacier sits on a retrograde slope bedrock, which is consistent with previous studies of MISI (e.g^13,52.). The glacier is also reacting to the loss of its ice shelf in 2003, and the warmer atmospheric conditions that appear to remove the winter short section of ice shelf that forms seasonally. The retreat under MISI is characterized by insensitivity to atmospheric forcing, and can only be stabilized by addition of buttressing, or when the glacier rests on a prograde slope at the terminus. This effect can be observed in the simulations where grounding line retreat is slowed on the prograde slope at about 40 km in Fig. 6. Under relatively cool atmospheric forcing (Fig. 6), the glacier does not retreat past this bump in the bed.

Our model appears more sensitive to atmospheric forcing than it is to ocean forcing. But keep in mind that ocean and atmospheric forcing are coupled in our model. Submarine melting of the glacier is influenced not only by ocean temperature coming through the fjord mouth, but also by the subglacial discharge which comes from surface runoff. Ice calving is also impacted by surface runoff in our crevasse depth dependent calving criterion. We may exercise some skepticism about over-interpreting results from this model as the simulations of the historical terminus position was not particularly good. The real glacier is sensitive to ocean forcing as well as atmospheric forcing; but it is not clear whether this sensitivity will continue in the future, even though we know that the central trough is deep and much of the ice discharge in the basin will remain in direct contact with the ocean even though the margins of the basin will retreat away from the ocean inland. So, the terminus will continue to be exposed to deep warm waters even after a large retreat.

Our model is an offline, asynchronous coupled ice sheet-plume model. To better understand the ice-ocean interaction at the glacier terminus, asynchronous or synchronous coupled ice sheet-plume models (e.g^46,53,54.) and ice sheet-ocean models (e.g^55,56,57.) have been developed over the last decade. Synchronously coupling the ice sheet model with a 3D ocean model is more computationally expensive and more challenging, since ocean models have difficulties accounting for continuously changing ice-margin geometry. But they can simulate both small- and large-scale spatial and temporal evolution of the basal melting beneath ice shelves. However, they require more oceanographic observations which are scarce, and more fjord bathymetry data which are mostly missing or inaccurate than the plume model, which questions the benefit of using computationally expensive 3D models for future predictions.

Conclusion

The barrier scheme we evaluated cannot stop Sermeq Kujalleq retreat in the next few decades under either cooler or warmer climates than present, because the glacier now sits on a submarine retrograde slope bedrock. Our simulations cannot stabilize the glacier even by using the repeated historical runoff and SMB during 2010–2019 or 1990–1999 – which is not very different from the 1960s forcings. However, use of historical atmospheric forcing can largely reduce the retreat especially during the 2070–2100 period, when there is a large increase of surface runoff under SSP5-8.5. This finding aligns with the assessment of the indigenous knowledge holders in Greenland: already before modeling was made, they expressed doubts on the effectivity of the curtain in Ilulissat Icefjord, because of surface melt⁵⁸.

This result suggests that some glaciers undergoing rapid retreat at present cannot be stabilized under geoengineering scenarios, either by cooling the atmosphere such as with stratospheric aerosols geoengineering¹², or by the targeted interventions explored by the barrier simulations here. This conclusion is commensurate with a study of local acceptability of an intervention in Ilulissat⁵⁸. The locals fear that erecting a barrier would harm their traditional indigenous way of life and economy where fishing and hunting play important role. Other reasons include the distrust in the efficiency and technical feasibility of the curtain, the special value of the Icefjord locally and globally, and the fear of colonialism if a potentially harmful curtain was placed in Ilulissat for preventing global sea level rise. Hence, based on workshops, interviews and informal discussions with over 50 participants in Ilulissat and Nuuk between 2021–2023, there appears to be neither a local or global benefit from making an intervention in Ilulissat. However, there are several other fast flowing and notable outlets for global sea level in Greenland⁵⁹. While human settlements in Greenland are typically located close to productive fjords³⁸, it is plausible that there are other fjords that would be suitable for field tests, with the consent of Greenland authorities, where there are no local fishing communities relying on the productivity of fjords resulting from the enhanced nutrients supplied by accelerating melt and iceberg discharge. Several of these are in north and northeast Greenland where both oceans and atmospheric thermal forcing is less than for Sermeq Kujalleq. It may be possible to prevent retreat in a MISI regime, such as forecast for Zacharia glacier⁶⁰. Recent observations of the 79 N glacier⁶¹, located next to Zacharia glacier suggest that its fringing ice shelf is being rapidly thinned by warming Atlantic waters. Thus, a barrier may have better potential in stabilizing that glacier than others in Greenland. In any case it appears that Sermeq Kujalleq already passed the tipping point of MISI, and the time for an intervention to make a difference to Sermeq Kujalleq is gone.

Method

Glacier data and Climate Forcing data

The bedrock elevation comes from IceBridge BedMachine Greenland, version 3 at a resolution of 500 m (Fig. 1; Table 2). We take the year 2009 as our initial conditions for simulations. The ice thickness comes from IceBridge MCoRDS L3⁶² at a resolution of 500 m for the year 2009. The observed surface velocity in the year 2009 is from MEaSUREs⁶³ at a resolution of 200 m. The observed terminus positions for the period 2009–2021 are from TermPicks⁶⁴ and Loebel et al.⁶⁵ The observed surface velocity for the period 2009–2023 is from MEaSUREs Greenland Ice Velocity, version 4⁶⁶.

Table 2 Data used in this study

Full size table

The modeled evolution of Sermeq Kujalleq over the 21st century is driven by surface runoff, SMB and ocean thermal forcing. Surface runoff and SMB for the period 1990–2100 are from monthly outputs of regional surface mass and energy balance model MARv3.11.3 with 20 km resolution¹¹. MARv3.11.3 is driven by historical forcing data for the period 2009–2014, and CNRM-ESM2-1 for the period 2015–2100⁶⁷ under SSP2-4.5 and SSP5-8.5 scenarios. Figure 4 shows that almost all surface runoff and meltwater discharge at the grounding line occurs in summer (JJA).

The water in Disko Bay can flow over the 250 m deep sill at the mouth of Ilulissat Fjord. The ocean temperature and salinity above 250 m depth and close to the fjord mouth comes from the monthly reanalysis data CMEMS (https://marine.copernicus.eu/) at spatial resolution of 0.083° × 0.083° for the period 2009–2019, and bias corrected modeled monthly outputs of CNRM-CM6-1-HR for the period 2020–2100. We select CNRM-CM6-1-HR because this has very high resolution for the ocean (25 km) compared with other climate models around Greenland. We do bias correction to the modeled monthly ocean temperature and salinity using the equidistant cumulative distribution functions (CDF) matching (EDCDFm) method. This follows standard procedures in bias correction which tries to account for probability distribution changes, as well changes to mean, between the projection and baseline periods for the modeled variable⁶⁸. There is an apparent 3-year ocean water cooling period 2016–18 which can be most easily seen at depths greater than 100 m shown by the reanalysis data CMEMS. And there is abrupt ocean water warming, most obvious in January and April, starting in 2019 (Fig. 2), due to switching the ocean forcing data.

Ice sheet model

We model Sermeq Kujalleq using the ice sheet model BISICLES^69,70. BISICLES employs a vertically integrated ice flow model based on Schoof and Hindmarsh⁷¹ which includes longitudinal and lateral stresses and a simplified treatment of vertical shear stress which is best suited to ice shelves and fast-flowing ice streams.

The ice thickness h and horizontal velocity u satisfy a two dimensional mass transport equation

$$frac{partial h}{partial t}+nabla cdot ({{bf{u}}}h)=a-M,$$

(3)

and two dimensional stress-balance equation

$$nabla cdot left[varphi hbar{mu }left(2dot{{{boldsymbol{epsilon }}}}+2{tr}left(dot{{{boldsymbol{epsilon }}}}right){{bf{I}}}right)right]-{{{boldsymbol{tau }}}}^{{{boldsymbol{b}}}}={rho }_{i}{gh}nabla s$$

(4)

where a is surface mass balance (SMB), M is the oceanic met rate only applied to the ice shelf bottom, φ is the stiffening factor, (hbar{mu }) is the vertically integrated effective viscosity, (dot{{{boldsymbol{epsilon }}}}) is the horizontal strain-rate tensor, I is the identity tensor, ({rho }_{i}) is the ice density, g is gravity acceleration, s is the surface elevation⁶⁹. The submarine melt rate in Eq. (3) beneath the ice shelf of Sermeq Kujalleq is not estimated by a simple parameterization, but using the output of a fjord model described later.

BISICLES employs adaptive mesh refinement to obtain fine spatial resolution in fast flow region and coarse resolution in low flow region. We use 500 m in the region of grounding line, and 1 km upstream and 2 km far inland (Supplementary Fig. 5).

We use a crevasse-based calving criterion that calves ice where the crevasse penetration depth, d_s, is greater than upper surface elevation, which is based on the idea that the downglacier velocity gradient and the ice elevation above the waterline are the primary controls on glacier terminus position⁷². We also include the effects of surface runoff in the opening of crevasses and the rigid ice mélange in suppressing calving on Jakobshavn Isbræ^73,74. Then d_s is estimated as

$${d}_{s}=frac{alpha cdot S}{gcdot {rho }_{i}}+frac{{rho }_{w}}{{rho }_{i}}cdot Rcdot beta$$

(5)

where S is the magnitude of tensile stress and R is surface runoff; α is a tuning coefficient, which is linearly dependent on ocean temperature to represent the buttressing of ice mélange in the fjord:

$$alpha ={alpha }_{1}+{alpha }_{2}cdot frac{({T}_{o}-mu )}{{sigma }^{2}}$$

(6)

where T_o is the average water temperature above 250-m-depth in the Ilulissat Fjord, μ and ({sigma }^{2}) are mean and variance values of T_o during the simulation period respectively; β is a tuning parameter of runoff. The parameters α and β are tuned to best reproduce the observed calving front position for the period 2009–2019.

A linear Weertman sliding relation is used for grounded ice to define the basal friction:

$${{{boldsymbol{tau }}}}_{b}=-C{{{bf{u}}}}_{b}$$

(7)

where ({{{boldsymbol{tau }}}}_{b}) is the basal shear stress, ({{{bf{u}}}}_{b}) is the basal velocity, and C is the basal friction coefficient.

The basal friction coefficient C and the stiffening factor φ (Eq. 4) are fixed throughout the simulations and obtained simultaneously by solving the inverse problem⁷⁰ with bedrock topography, ice thickness and ice surface velocity in the year 2009 (Table 2).

Fjord model

We consider an elongated fjord roughly representative of the geometry at Illulissat Icefjord in front of Sermeq Kujalleq (Fig. 1). This picture (Fig. 8) is inspired by the detailed observational and modeling description of Ilullissat circulation in Gladish et al.⁷ That paper described a fjord whose properties and internal overturning circulation are essentially determined by two competing plumes: at the seaward end, dense inflowing water cascades over a shallow sill to fill the deep basin of the fjord, while at the glacier end, buoyant meltwater creates a buoyant upwelling plume. Our model is customized to capture these specific dynamics.

In our model, the seaward end of the fjord consists of a shallow sill, while the landward end is the ice sheet terminus. At the seaward end, a two-layer exchange flow with the surrounding ocean results in dense water entering the fjord at relatively shallow depths. This inflow is denser than the ambient fjord properties and thus it cascades down the spillway slope at the back of the sill as a negatively buoyant turbulent plume. When this plume reaches its level of neutral buoyancy, its outflow spreads horizontally into the ambient fjord waters. At the glacier end of the fjord, fresh subglacial discharge triggers another plume, this one positively buoyant. The glacier plume entrains ambient deep water from the body of the fjord; the heat from this deep water triggers melting at the ice face, and the mixture of entrained water and ice melt is upwelled and then discharged into the ambient shallow layers of the fjord. As the two outflows spread out horizontally at their levels of neutral buoyancy within the fjord, they displace ambient water masses up and down to maintain a stable stratification. The vertical stratification of the fjord is thus determined by the balance between entrainment and outflow of these two plumes, and our model therefore consists of three coupled 1D models: two plume models (a positively buoyant one at the glacier face and a negatively buoyant one on the sill spillway), that communicate with each other via a vertical advection-diffusion equation for keeping track of the ambient fjord properties in the middle.

In our simplified model, the ambient vertical profiles of water properties within the fjord are determined by the solution to the 1D vertical advection/diffusion equation, with source and sink terms corresponding to the entrainment and outflow of the two plumes, and an additional term to account for parameterized iceberg melt. The evolution equations for conservative temperature ({{Theta }}) and absolute salinity S in depth z and time t are thus:

$$frac{{{partial }}{{Theta }}}{{{partial }}t}=-{w}_{f}frac{{{partial }}{{Theta }}}{{{partial }}z}+frac{{u}_{s}}{l}({{{Theta }}}_{s}-{{Theta }})+frac{{u}_{g}}{l}({{{Theta }}}_{g}-{{Theta }})+frac{{m}_{i}}{l}({{{Theta }}}_{i}-{{Theta }})+kappa frac{{{partial }}{{{Theta }}}^{2}}{{{{partial }}}^{2}z}$$

(8)

$$frac{{{partial }}S}{{{partial }}t}=-{w}_{f}frac{{{partial }}S}{{{partial }}z}+frac{{u}_{s}}{l}({S}_{s}-S)+frac{{u}_{g}}{l}({S}_{g}-S)-frac{{m}_{i}}{l}S+kappa frac{{{partial }}{S}^{2}}{{{{partial }}}^{2}z}$$

(9)

where w_f is the fjord upwelling velocity; u_s and u_g are the outflow velocities from the sill and glacier plumes, respectively, with associated outflow temperature and salinity given by Θ_s, Θ_g, S_s, and S_g; m_i is the iceberg melt rate, described below; Θ_i=−L/c_p is the equivalent water temperature of frozen ice, where L is the latent heat of melting and c_p is the specific heat of seawater (the salinity of frozen glacier ice is assumed to be 0); l is the fjord length from glacier to sill; and κ = 10^-3 m²s^-1 is a vertical diffusivity added to stabilize the problem.

The fjord upwelling velocity is computed by integrating the net entrainment and outflow of the two plumes:

$${w}_{f}(z)=frac{1}{l}{int }_{{z}_{b}}^{z}({u}_{s}+{u}_{g}-{e}_{s}-{e}_{g}){dz}^{{prime}}$$

(10)

where e_s and e_g are the entrainment rates of the sill and glacier plumes, respectively, and z_b is the depth of the fjord bottom. Note that, while we do not explicitly impose the boundary condition that w=0 at the sea surface, we implicitly enforce this condition by the fact that both plumes conserve mass, such that the integrated entrainment and outflow terms balance.

The two plumes follow the turbulent buoyant plume model of refs. ^43,75. They are governed by four equations describing the conservation of mass, momentum, heat, and salt. At the ice face the plume model uses the simplified two-equation melt model from Jenkins⁴³. We use the parameter values from Table 1 of Jenkins⁴³ for both plume models. For the sill plume we disable melting and integrate the plume model downwards along the spillway, with an assumed slope of 3.4°. The initial conditions for the sill plume are derived from the sill exchange model, described below. For the glacier plume, we integrate the plume equations upwards along the vertical ice face starting from an initial subglacial discharge. The plume model continues integrating beyond the neutral buoyancy depth so as to capture the magnitude of dynamic overshoot; this overshoot is used to compute a vertical spreading function which determines the depth distribution of the plume outflow into the rest of the fjord (see Supplementary Information).

The final component needed by our model is an expression for the exchange flow over the sill that connects the fjord with the surrounding ocean. This two-layer critical flow provides the starting boundary conditions for the dense sill plume. For this purpose, we use a theoretical treatment for critical exchange flows inspired by Farmer and Armi⁷⁶. The geometry of the two-layer exchange flow is shown in Fig. 4b. For any given value of the interface height, we can define layer thicknesses h₁ and h₂, along with average densities ρ₁ and ρ₂ given by the properties of either the fjord profile (for the outflow layer 1) or the external ocean (for the inflow layer 2), averaged over the appropriate depth range. The densities can be used to compute a reduced gravity, g’=Δρ/ρ₀. Our unknowns are u₁ and u₂, the exchange velocities. To conserve the overall mass of water in the fjord, the exchange velocities are related by,

$${u}_{1}{h}_{1}-{u}_{2}{h}_{2}={q}_{0}$$

(11)

where q₀ is the net barotropic flux out of the fjord (this term is necessary to balance the inflowing ice, but in practice it is quite small in comparison to the exchange fluxes). The condition that the flow over the sill must be a hydraulically controlled critical flow, as observed by Gladish et al.⁷ can be expressed as:

$${G}^{2}={F}_{1}^{2}+{F}_{2}^{2}={left(frac{{u}_{1}}{sqrt{{g}^{prime} {h}_{1}}}right)}^{2}+{left(frac{{u}_{2}}{sqrt{{g}^{prime}{h}_{2}}}right)}^{2}=1$$

(12)

where G is the composite Froude number of the exchange flow, composed of the Froude numbers for the individual layers F₁ and F₂. Using mass conservation to express u₁ in terms of u₂, we can then rearrange the above equation to form a quadratic equation for u₂:

$$({h}_{1}^{3}+{h}_{2}^{3}){u}_{2}^{2}+2{h}_{2}^{2}{q}_{0}{u}_{2}+{q}_{0}^{2}{h}_{2}-{g}^{prime} {h}_{1}^{3}{h}_{2}=0$$

(13)

We test all possible interface heights between the top of the sill or barrier and the sea surface, and compute u₂ (and thus u₁) for each of them. Finally, by making the assumption that the system will maximize entropy production when the magnitude of the exchange flow is maximized, we can select the interface depth that separates the inflow and outflow layers, thus completing the problem. If a barrier is placed at the sill the calculation of exchange flow is the same, except confined to a shallower region above the barrier.

We assume the runoff at the glacier surface is all discharged at the grounding line, which drives the plume near the ice front. We assume the subglacial discharge comes out from the grounding line through a 5 km wide water sheet. We use the subglacial discharge per kilometer (Fig. 4). The modeled melt rate at the ice front given by the fjord model is at the vertical cliff face of the ice terminus. It primarily depends on fjord water temperature and subglacial discharge which drives circulation. Melt rate can be increased by an increase in either variable, although it increases linearly with fjord water temperature and only sub-linearly with subglacial discharge (e.g³⁵.). We convert the maximum melt rate at the vertical ice front to an equivalent melt rate beneath the seasonal ice shelf using the grid cell size and ice thickness ratio.

The fjord model we use is adapted specifically for the conditions in Ilulissat Icefjord: the glacier end of the fjord is a nearly vertical ice cliff, justifying the use of a vertical plume model; the fjord itself is narrow and choked with icebergs, inhibiting the formation of eddies or large-scale cross-fjord structure; and the entrance has a relatively shallow sill, with two-layer hydraulically controlled exchange flow occurring in the space above the sill⁷. This model would be inappropriate to use for fjords where these conditions are not met, for example if the sill is deep enough to permit three-layer or more complex exchange flows (e.g⁷⁷.). We acknowledge that this fjord model is simple and neglects important fjord parameters including the effects of mixing at sills the earth’s rotation, and fjord-wide circulation. Slater et al.⁴⁵ suggests plumes can drive an energetic fjord-wide circulation which effectively doubles the glacier-wide melt rate compared to estimates of melting within plumes alone, and melting driven by fjord-scale circulation should be considered. However, the primary goal of the fjord model is to modify an ocean temperature boundary to produce a thermal forcing on the glacier and not to recreate fjord circulation perfectly. This is like other glacial fjord box models and parameterizations (e.g^78,79,80.). Nevertheless, we would welcome 2D and 3D ocean general circulation models, off- or on-line coupled to an ice sheet model to repeat the work and compare with our result.

Initialization, calibration and coupled model scheme

We assume that the bedrock topography beneath the ice sheet is unchanged in decades. We used bedrock topography in the year 2017, ice thickness and velocity in 2009 (Table 2) to provide a state in 2009. The initialization and calibration was performed as follows.

1.

We solved the inverse problem for basal conditions (Eq. 7) and stiffening factor following Cornford et al.⁷⁰ Supplementary Fig. 6 shows the discrepancy between observed velocity field⁶³ and the velocity derived from the inversion.
2.

Using φ and C from the inversion of Step 1, we performed a surface relaxation for 0.1 years to reduce errors caused by inconsistencies in data coverage and observation time.
3.

We ran the fjord model over the 2009–2019 period to produce water temperatures inside Ilulissat Fjord needed for Eq. (5), and the monthly basal melt rate at the ice front, forced by the water properties and velocity at the fjord mouth, and subglacial discharge through the grounding line.
4.

We used many individual decade simulations to determine the parameters α₁, α₂ and β, which represent mélange buttressing and crevasse depth sensitivity to surface runoff. We minimized that mismatch between observed and modeled monthly calving front position during the period 2009–2019.

The best set of parameters are α₁ = 0.6, α₂ = 0.15, β = 34 (Fig. 9). The second-best set of parameters are α₁ = 0.5, α₂ = 0.11, β = 36 and these we use in a sensitivity analysis of glacier response to parameter selection. The modeled timings of maximum extent and minimum extent each year are in satisfactory agreement with observations, and the error distribution shows a single minimum (Supplementary Fig. 7). The glacier terminus reaches maximum extent with slowest ice velocity in December and January, and minimum extent with highest velocity between June and August. The modeled calving front shows net retreat with seasonal fluctuations, and retreats 3.5 km from summer of 2009 to the summer of 2015, which is consistent with observations. Ocean temperature cooling from 2015 to 2017 (Fig. 2) triggered substantial advance of calving front (Fig. 9), demonstrating the sensitivity of our ice mélange buttress parameterization to ocean forcing. Data from after 2019 were not included in the parameter fitting, but the model captures the calving front position and velocity reasonably well (Fig. 9), and better than the second-best parameter set (Supplementary Fig. 8).

Fig. 9: Comparison of modeled and observed terminus position from July 2009 to December 2021 and ice flow velocity near the ice front from July 2009 to August 2023 along the central ice flowline (Fig. 1c) using α₁ = 0.6, α₂ = 0.15, β = 34.

We use asynchronous coupling between the ice sheet model and the fjord model for the projections from 2020 to 2100. The water temperature inside Ilulissat Fjord and basal melt rate under the ice shelf from 2009 to 2100 come from the fjord model. The time step of the fjord model is 1/4 day, while that of BISICLES is one month. The monthly basal melt rate is used for forcing the ice sheet model. The ice sheet model yields the glacier terminus position change, which corresponds to the fjord length change. Fjord length essentially controls how quickly fjord properties can change in response to changes in the incoming fluxes. Small changes of fjord length have little influence on the modeled fjord water properties. We chose to update the fjord length every 10 years to balance computation cost and accuracy of the evolution of glacier terminus, but we use the modeled monthly basal melt rate from the fjord model to drive the ice sheet model.

We designed the following scheme for coupling the ice sheet model BISICLES and the fjord model for the period 2020–2100:

1.

We repeated the year 2020 drivers for 3 years to spin-up the fjord model. We then run the fjord model for the first decade: 2020–2029.
2.

We applied the monthly average of modeled water temperature and basal melt rate from the fjord model to the ice sheet model BISICLES, and run BISICLES for 10 years. This produced the modeled glacier velocity and the evolution of glacier geometry including glacier terminus position.
3.

We updated the length of fjord with the last glacier terminus position from step 2, then repeated steps 1 and 2 for the next decade. Each new decadal simulation is preceded by a 3-year spin up of fjord model with repeated forcing from the last year of the last decade.