Quantifying the buttressing contribution of landfast sea ice and melange to Crane Glacier, Antarctic Peninsula

The January 2022 disintegration of multi-year landfast sea ice in the Larsen B embayment, Antarctic Peninsula, was closely followed by a significant acceleration of ice flow and ice-front retreat of numerous outlet glaciers. Crane Glacier was a notable example of this, with 6 km of its floating ice shelf lost to calving in the first month following the disintegration and a 3.4 % increase in terminus flow speeds over the same time period. In this study we quantify for the first time the buttressing stresses that were transmitted to Crane by the ice melange at the glacier outlet using the ice-flow model Úa. We constrained our model with high-resolution surface elevation profiles of the glacier and ambient melange and reconstructed the observed flow velocities by optimising the rheology rate factor throughout our model domain. This allowed us to quantify the stress regime across both the glacier and ice melange. Results showed that resistive backstresses were imparted to Crane by the ice melange with a mean buttressing ratio of Θ_N=0.68 calculated at the glacier terminus (Θ_N=1 implies no buttressing). In addition, diagnostic modelling showed an expected 19.2 kPa mean increase in extensional stress at the ice front following the disintegration of the ice melange. This perturbation in stress likely triggered the observed rapid calving over the near-terminus region, leading to the periodic loss of sections of Crane's buttressing ice shelf and thus further acceleration of ice flow in the subsequent months.

1 Introduction

Landfast sea ice (fast ice) is a type of stationary sea ice which is fastened to coastlines, ice shelves and outlet glaciers. Despite its potential to regulate the dynamic behaviour of adjoining glaciers (Massom et al., 2010, 2015; Arthur et al., 2021; Christie et al., 2022), gaps in knowledge relating to the distribution and extent of fast ice have led to it being overlooked in many studies (Fraser et al., 2023). Recent observations have significantly improved our understanding of trends in fast-ice extent throughout Antarctica (Fraser et al., 2021), changes to which could impact upon the stability of adjoining ice shelves and outlet glaciers, therefore potentially impacting future rates of sea level rise (Massom et al., 2018; Miles et al., 2017).

Whilst the buttressing capacity of ice shelves has been assessed quantitatively in numerous studies (Fürst et al., 2016; Reese et al., 2018; Gudmundsson et al., 2023), buttressing provided by fast ice has primarily been discussed in terms of the correlation between perturbations in fast-ice coverage and observed changes in ice-flow velocity or migration of ice tongue calving fronts. The seasonal acceleration of the Totton Ice Shelf (Greene et al., 2018) and the Parker Ice tongue (Gomez-Fell et al., 2022) has been attributed to changes in fast-ice extent, as has the intermittent acceleration of the Thwaites Glacier tongue (Miles et al., 2020). However, little to no change in glacier flow velocities was reported following fast-ice breakout events at Shirase Glacier (Nakamura et al., 2010, 2022). The contrasting responses reported in these studies suggest that either seasonal or local geometric factors may ultimately control the extent to which fast ice may stabilise connected glaciers (Fraser et al., 2023).

More recently, the 2022 widespread disintegration of multi-year fast ice in the Larsen B embayment, Antarctic Peninsula, provided an opportunity to assess the impact of instantaneous fast-ice loss on multiple glaciers, aided by abundant present-day observational datasets for the region (Sun et al., 2023; Ochwat et al., 2024; Surawy-Stepney et al., 2024). Furthermore, the capacity of fast ice to buttress glaciers was brought into question due to the similarities in the response of the regions' outlet glaciers following this disintegration event and the 2002 collapse of the Larsen B Ice Shelf (Scambos et al., 2004; Rignot et al., 2004; De Rydt et al., 2015).

Following the 2002 ice-shelf collapse, Crane Glacier (Fig. 1) retreated by more than 10 km in just over 2 years (Needell and Holschuh, 2023) with ice flow across the grounding line increasing by up to 3 times over a similar period (Rignot et al., 2004). Crane's longer-term response was more complex, and following the initial rapid retreat, phases of arrest and subsequent re-advance were exhibited (Needell and Holschuh, 2023; Rott et al., 2018; Wuite et al., 2015).

Figure 1 Panel (a) shows an overview of the Larsen B embayment in polar stereographic projection displayed on a hill-shaded relief map of the REMA mosaic for Antarctica (Howat et al., 2019). Latitude and longitude lines are labelled for reference. The location of Larsen B is illustrated by the red ring on the inset outline of Antarctica (source: BedMachine v3; Morlighem, 2022). The black line outlines the extent of the model domain. The blue line denotes the approximate extent of fast ice in the months prior to its disintegration, and the red line shows the fast ice and melange extent after disintegration, with coordinates extracted from Landsat 8 imagery (24 February 2022). Panel (b) shows a close-up view of Crane Glacier captured by Landsat 9 (13 December 2021), with the black line showing the extent of the model domain.

Sea ice formed seasonally in the Larsen B embayment each year following the 2002 ice-shelf collapse, but in 2011 the sea ice became landfast before further developing into multi-year fast ice. Multi-year fast ice not only is thicker than seasonal formations of sea ice (i.e. first-year ice) (Arthur et al., 2021; Massom et al., 2010), but can also display enhanced mechanical strength through desalination over summer months (Fraser et al., 2023). This thicker, stronger ice allowed the trapping and bonding of fragments of partially calved icebergs in a pro-glacial melange (Ochwat et al., 2024; Surawy-Stepney et al., 2024). The thickness of the multi-year fast ice was estimated to be between 2.5 and 4 m across the embayment (Scambos et al., 2017), though thicknesses of magnitudes of tens to hundreds of metres were reported for melange elements close to the glacier outlets (Ochwat et al., 2024). The development of multi-year fast ice seemingly aided a decade long re-advance of Crane, but its disintegration in mid-January 2022 triggered a second period of rapid retreat (Needell and Holschuh, 2023). Ice-flow velocities were also affected with acceleration observed in the months following the disintegration, however a more significant increase in flow speeds over both grounded and floating ice was not observed until later in the year (Fig. 2e, Movie S2 in the Supplement).

Figure 2 (a–d) Crane's outlet in (a) December 2021, (b) January 2022, (c) February 2022 and (d) March 2022 (acquisition dates provided in Table 1), cropped from Landsat 8 and Landsat 9 satellite images. The fast ice was still intact in December 2021, and a degree of ice melange remained in the outlet of the fjord in the months following the fast-ice disintegration. The black line shows the modelled grounding line position, and the white line shows the digitised terminus location. The coloured dots are locations from which velocity data were extracted by ENVEO et al. (2021). Monthly averaged ice-flow velocities from 2015–2023 (ENVEO et al., 2021) are shown in panel (e). The blue and red lines represent the flow velocities at the corresponding blue and red location markers shown in the top panels with the vertical black lines corresponding to the dates of panels (a)–(d).

The similarities between the 2002 and 2022 events led to multiple studies assessing the role that the fast ice played in buttressing Larsen B's outlet glaciers, including, but not limited to, Crane Glacier (Sun et al., 2023; Ochwat et al., 2024; Surawy-Stepney et al., 2024). However, conclusions on the extent to which fast ice can buttress adjoining glaciers varied across these studies. Sun et al. (2023) argued that whilst the fast ice covered the same area as the previously existing ice shelf, the buttressing provided was low as damage could more readily propagate through the fast ice due to its weaker and thinner makeup, therefore limiting its resistive potential. This was evidenced by a lack of instantaneous acceleration at six glaciers following the disintegration and near-zero correlation between fast-ice extent and glacier velocities (Sun et al., 2023). A delayed velocity response, observed 8 months later, was attributed to the retreat of Crane's ice front rather than being directly due to the collapse of the fast ice (Sun et al., 2023). In contrast, Ochwat et al. (2024) recognised the similarities in dynamic response following the fast-ice and ice-shelf breakups, including instantaneous acceleration of ice flow, to mean there were parallels in the buttressing provided by the two different ice masses. The destabilisation and calving of floating ice tongues in different outlet glaciers in Larsen B were therefore attributed to the loss of buttressing from the fast ice (Ochwat et al., 2024). Surawy-Stepney et al. (2024) used an ice-sheet model to estimate the buttressing stresses of an idealised distribution of fast ice on Larsen B's outlet glaciers. The authors first optimised the rate factor of the glaciers to reproduce the observed velocity field without the inclusion of fast ice in the model and then added varying thicknesses of fast ice to the embayment region of the domain to calculate the resulting change in ice-flow speed. The study concluded that whilst the fast ice affected the dynamic behaviour of floating ice in different glacier outlets, this was not a result of direct buttressing in the same context as is understood for ice shelves. Instead, the impacts that fast-ice coverage had on calving behaviour and glacier flow velocities were attributed to secondary processes including ocean swell attenuation (Teder et al., 2022; Christie et al., 2022) and bonding of melange in areas downstream of glacier termini (Robel, 2017). The potential issue of this approach lies in the initialisation process: by solving the optimisation problem without including fast ice, the buttressing stresses from the fast ice, if any, have been compensated for by adjusting the rate factor of the glacier.

Here we seek to quantify for the first time the buttressing stresses that were imparted to Crane Glacier by the observed mixture of fast ice and pro-glacial melange. For this study we use the ice-sheet model Úa with a computational domain which includes the glacier as well as ambient fast ice and pro-glacial melange. We quantify buttressing using the buttressing ratio (Gudmundsson, 2013; Schoof, 2007) and further assess the change in the near-terminus stress field following the fast-ice disintegration, which may have triggered the retreat of Crane's ice front. Finally, we test the robustness of the ice-sheet model in determining the resistive stresses imparted to Crane by investigating the sensitivity of our results to changes in the prescribed input thickness of fast-ice and melange elements. This is to improve confidence in employing our methodology in scenarios where uncertainty in ice thickness measurements is of concern.

2 Datasets

In order to configure a model geometry representative of Crane Glacier prior to the disintegration of fast ice in the Larsen B embayment, we utilised a range of Reference Elevation Model of Antarctica (REMA) strip DEMs (Howat et al., 2019), timestamped between 30 September 2020 and 16 January 2022. The strips are defined at 2 m spatial resolution and have an absolute error of ±2 m in the horizontal and vertical planes. To create a continuous surface elevation profile across our model domain, each strip was added in turn, and where coordinates overlapped between strips, the most recent strip was used. As elevations at intersecting locations did not deviate beyond the associated error range, no further corrections were applied between intersecting strips. The near-terminus region and melange-filled outlet of Crane were populated by the strip dated 16 January 2022. The REMA strip DEMs are referenced to the WGS64 ellipsoid and were corrected for the geoid using values from BedMachine (Morlighem, 2022; Morlighem et al., 2020) interpolated at each nodal point of our model domain.

The bedrock elevation was defined using the Huss and Farinotti (2014) bedrock DEM, who used a mass continuity approach informed by NASA's Operation IceBridge and ground-based measurements to define a bedrock dataset across the Antarctic Peninsula at 100 m resolution. This dataset was merged with multibeam swath bathymetry data (Rebesco et al., 2014) which provided direct measurements of bedrock elevations in the vicinity of Crane's outlet and grounding line. However, due to uncertainty in the bedrock profile and recent observations suggesting that Crane's grounding line lay further downstream prior to the fast-ice disintegration (Wallis et al., 2024), we performed an additional sensitivity experiment considering the shallower BedMachine bed elevation data (Morlighem, 2022; Morlighem et al., 2020) to ensure that our conclusions would not be impacted by uncertainty in the estimated bed topography (Fig. S2).

The ice thicknesses across regions of floating ice were calculated from the flotation criterion considering the input surface elevations (Howat et al., 2019) and uniform ice and ocean densities of 917 and 1030 kg m⁻³ respectively. The grounding line is defined in our model where the calculated submerged ice thickness reaches the prescribed bedrock elevation (Huss and Farinotti, 2014).

Monthly averaged ice velocity maps at 200 m grid spacing were derived from successive Sentinel-1 IW SLC (Single Look Complex) image pairs (2014–2023) using a combination of coherent and incoherent offset tracking techniques (ENVEO et al., 2021; Nagler et al., 2021, 2015).

3 Methodology

3.1 Experiment design

We performed a series of numerical simulations to quantify the buttressing stresses exerted on Crane Glacier by the fast ice and pro-glacial melange. The simulations were performed in two steps, an overview of which is given below with further details described in the following subsections.

In the first step, the rate factor, A in Glen's flow law, was estimated through the inversion of measured velocities, allowing us to reconstruct the stress field throughout Crane and the ambient ice melange. We used satellite imagery to extract the location of the calving front prior to the fast-ice disintegration and calculated the normal resistive stress at these coordinates in order to quantify the buttressing provided by the fast ice and melange. We made the same calculations at the location of the terminus in the days immediately after the fast-ice disintegration and again at locations that the ice front had retreated to in the following months (Table 1). These terminus locations were considered in order to compare the buttressing stresses provided by the fast ice and melange with those provided by the portions of Crane's floating ice shelf which were lost to calving in the months following the fast-ice disintegration.

In the second step, we investigated the expected change in the near-terminus stress regime following the fast-ice disintegration. After having inverted for A and ensured that modelled ice velocities were in good agreement with observations, we perturbed the model in a diagnostic simulation by removing the ice melange from the domain and assessed the instantaneous change in stress regime.

Table 1 The terminus IDs referred to throughout the study represent the terminus locations extracted from satellite imagery with acquisition dates shown below. Landsat imagery was obtained courtesy of the U.S. Geological Survey.

3.2 Model set-up

We used the ice-flow model Úa (Gudmundsson, 2020), which solves the governing equations of ice dynamics using the shallow shelf approximation (SSA) (Morland, 1987; MacAyeal, 1989). The model uses a 2D vertically integrated approach which allowed us to assess the stress distribution throughout the domain in lateral and transverse directions, enabling an investigation into buttressing from the fast ice and melange where effects of lateral drag from the margins of the fjord are captured.

The model domain included the entirety of Crane, along with all in-flowing tributary regions of the glacier, and extended downstream from the outlet of Crane's fjord (Fig. 1). The locations of rock outcrops were defined using Landsat imagery, and holes in the mesh were placed at these locations where areas of thin ice and high strain rates may have caused numerical difficulties for the ice-flow model (De Rydt et al., 2015). Boundary conditions at the edges of the domain were fixed to ice-flow velocities from observational data (ENVEO et al., 2021). A zero flow condition was imposed at the interior boundaries.

The finite element mesh was refined to a minimum resolution of 100 m over areas of fast ice and melange and regions close to the terminus and grounding line positions of the glacier. To reduce computational cost, a coarser resolution was employed further upstream and varied based on flow velocities up to a maximum resolution of 2.5 km. A total of 36 268 elements with mean size of 194.2 m made up the model domain.

3.3 Inversion

Model parameters for basal slipperiness (C) and ice rheology (rate factor, A) were determined through an inversion process (MacAyeal, 1993) based on ice velocity measurements following a commonly used methodology (e.g. Hill et al., 2018; Barnes et al., 2021; Sun and Gudmundsson, 2023). Ice rheology is assumed to follow Glen's flow law (Glen, 1955), with stress exponent n=3, and basal sliding is assumed to follow Weertman's sliding law (Weertman, 1957), with stress exponent m=3.

Prior values for A and C were chosen based on assumed temperature and flow velocity respectively before iteratively updating these values to minimise the cost function $J=I+R$. Here, I represents the misfit between modelled and observed velocities, and R is a regularisation parameter for which we use the Tikhonov regularisation of the form

$$\begin{array}{}\text{(1)}& \begin{array}{rl}R& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left[{\mathit{\gamma }}_{\mathrm{s}}^{\mathrm{2}}A{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}\left({\displaystyle \frac{A}{\stackrel{\mathrm{̃}}{A}}}\right)\right)}^{\mathrm{2}}\right.\\ & +{\mathit{\gamma }}_{\mathrm{s}}^{\mathrm{2}}C{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}\left({\displaystyle \frac{C}{\stackrel{\mathrm{̃}}{C}}}\right)\right)}^{\mathrm{2}}+{\mathit{\gamma }}_{\mathrm{a}}^{\mathrm{2}}A{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}\left({\displaystyle \frac{A}{\stackrel{\mathrm{̃}}{A}}}\right)\right)}^{\mathrm{2}}\\ & \left.+{\mathit{\gamma }}_{\mathrm{a}}^{\mathrm{2}}C{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}\left({\displaystyle \frac{C}{\stackrel{\mathrm{̃}}{C}}}\right)\right)}^{\mathrm{2}}\right]\mathrm{d}\mathcal{A},\end{array}\end{array}$$

where 𝒜 is the domain area, $\stackrel{\mathrm{̃}}{A}$ and $\stackrel{\mathrm{̃}}{C}$ are prior estimates for A and C, and γ_s and γ_a are regularisation parameters penalising slope and amplitude respectively. We performed L-curve analysis (Hansen, 1992) to optimise the regularisation parameters, iterating until convergence.

The inversion was based on monthly averaged velocity measurements from November 2021. This negated the possibility of including any contamination that may be present in the January 2022 dataset, which partly accounts for time periods both prior to and following the fast-ice disintegration. Monthly averaged velocities from December 2021 were discounted due to high errors associated with the velocity measurements downstream of the terminus location; however results from sensitivity experiments show the choice of velocity dataset used in the inversion has little impact on the findings of the study (Fig. S3).

3.4 Diagnostic modelling

After having inverted for A and ensured that modelled ice velocities were in good agreement with observations, we performed a diagnostic simulation in order to assess the instantaneous change in the near-terminus stress regime following the disintegration of the fast ice. We perturbed the model by removing all ice from the model domain downstream of the December 2021 terminus location (Table 1). It is noted that for numerical reasons, removing the fast ice and melange actually involves replacing the defined thickness from the REMA strip DEMs (Howat et al., 2019) with a nominal 0.1 m thickness, which is suitably thin so as not to impact the upstream flow of ice.

Prior to running the diagnostic simulation, the boundary conditions at the domain's flow outlet were changed to a stress-free condition.

3.5 Buttressing quantification

Buttressing is commonly quantified by determining the normal resistive stress imparted by an ice shelf at the grounding line R_N and comparing this to the hypothetical ocean pressure which would be exerted at the same location in the absence of an ice shelf R_o (Schoof, 2007; Gudmundsson, 2013). Here we apply the same methodology, instead considering the vertically integrated normal resistive stress imparted by the fast-ice and melange elements at the glacier terminus:

$$\begin{array}{}\text{(2)}& {\displaystyle }{R}_{\mathrm{o}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{i}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{o}}}}\right)gh,\text{(3)}& {\displaystyle }{R}_N={\mathit{n}}^{⊺}\cdot \left(\mathbf{R}\mathit{n}\right),\text{(4)}& {\displaystyle }\mathbf{R}=\left[\begin{array}{cc}\mathrm{2}{\mathit{\tau }}_{xx}+{\mathit{\tau }}_{yy}& {\mathit{\tau }}_{xy}\\ {\mathit{\tau }}_{xy}& \mathrm{2}{\mathit{\tau }}_{yy}+{\mathit{\tau }}_{xx}\end{array}\right],\end{array}$$

where g is the gravitational constant, h is ice thickness, n is the unit vector normal to the terminus, R is the resistive stress tensor, and τ_ij are components of the deviatoric stress tensor computed in the model inversion. Values of 917 and 1030 kg m⁻³ are considered for ice density ρ_i and seawater density ρ_o respectively.

The buttressing ratio is given by Θ_N.

$$\begin{array}{}\text{(5)}& {\mathrm{\Theta }}_N={\displaystyle \frac{R_N}{R_{\mathrm{o}}}}\end{array}$$

The scenario of no buttressing is given by Θ_N=1, as in this example the resistive stress imparted by the fast ice and melange normal to the terminus would be equal to that of the vertically integrated ocean pressure at the same location. Θ_N<1 shows there are buttressing stresses resisting ice flow, with decreasing values corresponding to increasing resistive stress. Where Θ_N>1, resistive stresses are less than in the hypothetical ice-free scenario, indicating that ice is being pulled downstream. Where Θ_N<0, deviatoric stress normal to the terminus is compressive; otherwise it is tensile.

3.6 Fast-ice and melange thickness sensitivity

The inversion process provides an optimised estimate of the ice rheology over the glacier, fast-ice and melange areas of the model domain. Crucially, this process accounts for the discontinuous nature of the melange by returning a varied A (rate factor) field, derived from the measured velocities and dependent upon the input ice thicknesses by solving the conservation of momentum equations with SSA approximation:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\partial }_x\left(h\right(\mathrm{2}{\mathit{\tau }}_{xx}+{\mathit{\tau }}_{yy}\left)\right)+{\partial }_y\left(h{\mathit{\tau }}_{xy}\right)-t_{\mathrm{b}x}={\mathit{\rho }}_{\mathrm{i}}gh{\partial }_{x}s,\text{(7)}& {\displaystyle }{\partial }_y\left(h\right(\mathrm{2}{\mathit{\tau }}_{yy}+{\mathit{\tau }}_{xx}\left)\right)+{\partial }_x\left(h{\mathit{\tau }}_{xy}\right)-t_{\mathrm{b}y}={\mathit{\rho }}_{\mathrm{i}}gh{\partial }_{y}s,\end{array}$$

where h is ice thickness, t_bh is the horizontal components of the basal stress vector (equal to zero in floating regions), ρ_i is ice density, g is the gravitational constant, s is the ice surface elevation, and τ_ij are components of the deviatoric stress tensor calculated using Glen's flow law,

$$\begin{array}{}\text{(8)}& {\dot{\mathit{\epsilon }}}_{ij}=A{\mathit{\tau }}^{n-\mathrm{1}}{\mathit{\tau }}_{ij},\end{array}$$

where ${\dot{\mathit{\epsilon }}}_{ij}$ are components of the strain rate tensor, A is the rate factor, τ is the second invariant of the deviatoric stress tensor τ_ij, and n is the stress exponent.

The products of rate factor and ice thickness are found on the left-hand side of the momentum equations (Eqs. 6 and 7), with ice thickness only appearing separately from the rate factor in the body force component on the right-hand side of the equations. We therefore tested the hypothesis that uncertainty in fast-ice and melange thickness measurements is accounted for in the inversion, as the resulting ice rheology is described not solely by the rate factor but rather by the product of A⋅h.

We performed a series of sensitivity experiments by adjusting the prescribed surface elevation of fast-ice and melange elements downstream of the December 2021 terminus in increments of ±10 % up to a maximum 40 % change from the original surface elevations (Howat et al., 2019). In regions of thin ice, the minimum defined surface elevation was limited to 0.1 m.

For each adjusted ice thickness configuration, we performed a new model inversion. Using the updated rheology rate factor fields, we calculated the buttressing ratio at the terminus location in each sensitivity case. Our hypothesis would be proven if the calculated buttressing ratios were not affected by the adjustment in the input surface elevation of fast ice and melange. Results from these sensitivity experiments are presented in Sect. 4.3.

Figure 3 The buttressing ratio (Θ_N, Eq. 5) was defined at 100 m intervals along terminus positions derived from Landsat imagery representing the ice front in December 2021, January 2022, February 2022 and March 2022 (Table 1). Histograms of the buttressing ratios at these intervals along each month's terminus position are shown in panels (a)–(d) with a normally distributed probability density function overlain in black. The vertical dashed black lines show the mean buttressing ratio. Panel (e) shows the buttressing ratios in coloured bubbles along the terminus location with dashed yellow, blue, cyan and red lines showing the terminus location in December 2021, January 2022, February 2022 and March 2022 respectively, plotted on top of a surface elevation map of Crane's outlet region. The dashed black line shows the grounding line.

4 Results

4.1 Buttressing ratio

High spatial variability in the calculated buttressing ratios is seen at the December 2021 and January 2022 terminus positions (Fig. 3). Compression is seen in the central regions (red bubbles, Fig. 3), and the majority of buttressing ratios lie below 1, showing that resistive backstresses are imparted by the fast ice and melange. Buttressing ratios exceed 1 in some cases close to the edges of the fjord, indicating that ice at the margins is being pulled downstream. Mean buttressing ratios of 0.68 and 0.65 were found across the December 2021 and January 2022 terminus locations respectively with values between 0.8 and 0.9 occurring most frequently along these ice fronts.

The February 2022 and March 2022 terminus positions show a more consistent grouping of buttressing ratios with the majority of calculated Θ_N values lying between 0.6 and 0.7 (Fig. 3). Buttressing ratios do not exceed 1 at any point with mean values of 0.66 and 0.60 calculated along the February 2022 and March 2022 ice fronts respectively.

Figure 4 Panel (a) shows the principal deviatoric stress distribution in the near-terminus region of Crane determined from the inversion with fast ice and melange included in the domain. Extensional and compressive stresses are shown by the blue and red arrows respectively. A zoomed-in view of this stress distribution at the terminus is shown in panel (b). Panel (c) shows the change in the principle deviatoric stresses following the removal of fast ice and melange downstream of the December 2021 terminus, determined in the diagnostic model run. A zoomed-in view of this stress distribution at the terminus is shown in panel (d). The maximum and mean increases in the first principle stress along the terminus are 70.8 and 19.2 kPa respectively. Stress distributions in each panel are overlain on the ice surface elevation profile with the grounding line and terminus positions illustrated by the solid black and dashed yellow lines respectively.

4.2 Change in terminus stress regime

We compared Crane's stress regime from our inversion, representing the stress field in the glacier with fast ice still intact, against that from the diagnostic simulation in which fast ice and melange were removed downstream of the December 2021 terminus. The removal of melange elements from the domain led to an increase in extensional stresses in Crane's floating ice shelf running perpendicular to the ice front (Fig. 4). The magnitude of this change increased with closer proximity to the ice front where maximum and mean modelled increases in extensional stresses were found to be 70.8 and 19.2 kPa respectively. These values can be compared to the modelled basal drag values approaching the grounding line which have an average of 64.4 kPa.

4.3 Sensitivity of results to fast-ice and melange thickness

The buttressing ratios calculated across the terminus remain consistent across each sensitivity experiment with the spatial distribution of areas of high and low buttressing unchanged by the input ice thickness (Fig. 5). Greater resistive stresses are seen in the central region of the ice front with areas of negative buttressing (Θ_N>1) seen at the margins.

A 14.3 % increase in the mean value of Θ_N was calculated with ice thicknesses downstream of the terminus location reduced by 40 % (the thinnest thickness distribution assessed). In all other sensitivity experiments, the mean value of Θ_N differed by a maximum of 6.6 % from the reference configuration. No correlation between the percentage change in modelled thickness of ice melange and the deviation of mean buttressing ratio from the reference case results is observed. Uncertainty in the ice thickness profile therefore does not significantly affect our results.

Figure 5 The buttressing ratios calculated at the December 2021 terminus location with ice melange represented by varying thicknesses in the model. Updated model inversions were performed separately for each configuration. In panel (a), the fast-ice and melange thicknesses downstream of the terminus were defined using the REMA strip DEMs as per the initial intact model configuration. In panels (b)–(i), ice thicknesses downstream of the terminus were adjusted from the DEM-derived thicknesses by the ratio labelled on each panel. The solid black line shows the terminus location with buttressing ratios along this line shown by coloured bubbles at 100 m intervals.

5 Discussion

5.1 Buttressing contribution of fast ice and melange

The buttressing stresses imparted by multi-year fast ice to Larsen B's outlet glaciers have previously been assessed through the interpretation of observational data (Sun et al., 2023; Ochwat et al., 2024), as well as in an idealised modelling set-up assuming uniform thickness of fast ice across the Larsen B embayment (Surawy-Stepney et al., 2024). As no consensus as to the significance of the buttressing provided by the fast ice was reached between these studies, we sought to quantify this buttressing for the first time using a numerical modelling set-up which captured the realistic profile of pro-glacial melange in Crane's immediate outlet. With a mean buttressing ratio of 0.68 calculated at the glacier terminus (Fig. 3), our results show that this combination of fast ice and pro-glacial melange provided significant resistive backstress to Crane.

Negative buttressing ratios, which correspond to areas of compressive stress, were found in some locations across the December 2021 and January 2022 ice fronts (Fig. 3). Examination of satellite imagery (Fig. 2a) and the surface elevation profile of the terminus region (Howat et al., 2019) highlights the discontinuous nature of the damaged ice shelf and melange in the transition between glacier and fast ice. Fragments of calved icebergs were trapped by the fast ice, preventing transport away from the terminus and making the definition of an exact terminus location difficult to achieve. When we consider that the model reduces the resolution over this region of high variability in ice thickness to a 100 m grid, we anticipate that compression is found in regions of the model which are influenced by crevassing and fragmentation of the ice. These areas are unlikely to impart such high resistive stresses as the corresponding buttressing ratios suggest, which would result in the mean ratios being artificially shifted to a lower amount. If we instead consider the mode values calculated across the December 2021 and January 2022 terminus locations, we find buttressing ratios between 0.8 and 0.9 (Fig. 3). These values may therefore be more representative of the buttressing provided by the ice melange than the mean values reported above. Despite this, the mode values still support the conclusion that buttressing was imparted by the melange prior to its disintegration and that these buttressing stresses are lower than those provided by the ice shelf further upstream (Fig. 3).

The buttressing ratios presented above are subject to a degree of uncertainty due to the difficulty in defining the exact location where the transition between glacier and melange exists. By assessing the normal resistive stresses at 100 m intervals along the defined terminus location, we aimed to minimise this uncertainty by capturing the variability in ice thickness and rheology that may exist between melange elements. We have not tested the impact on the results that would be caused by changes to the defined terminus coordinates; however similar variability in the thickness and rheology between melange elements would be expected if the terminus location was moved elsewhere in the transition zone between the damaged ice shelf and melange region.

We note that further uncertainty may exist, arising from uncertainties in the datasets used in this study, as well as being caused numerically. We tested the sensitivity of results to the bed topography considered in the model geometry (Fig. S2) and also to the velocity field considered in the inversion (Fig. S3) to demonstrate the robustness of these results. We discuss the sensitivity of the results to the prescribed thickness of the fast ice and melange in Sect. 5.3.

5.2 Crane's response to the fast-ice disintegration

The quantification of buttressing from both the ice melange and the floating ice shelf (Fig. 3) provides insight into the observed response of Crane to the disintegration of Larsen B's landfast sea ice. Observational data show (1) an acceleration of ice-flow speed and (2) rapid calving of Crane's floating ice shelf.

The flow speed at the terminus increased by an average of 3.4 % in the first month following the fast-ice disintegration compared to when the fast ice was intact (Movie S2). A smaller increase of 0.5 % across the grounding line was observed over the same time period. A fluctuation between increase and decrease in flow speed across the grounding line was observed until August 2022 with average flow speeds typically within 1 % of the pre-disintegration velocities; however acceleration at the terminus region reached a maximum of 13.6 % over this period, suggesting that the mechanical support supplied to the glacier by the fast ice and melange played a greater role in the restriction of ice flow at the terminus than further upstream at the grounding line.

A total of 8 months later, in September 2022, increases in flow speed at the terminus and grounding line locations had increased to 14.10 % and 7.54 % respectively, with further increases observed in the following months (Fig. 2e, Movie S2). Sun et al. (2023) argued that such a delay in the velocity response was due to the retreat of Crane's ice front rather than being directly due to the disintegration of the fast ice. However, our results suggest that these phenomena are connected and ultimately caused by the fast-ice disintegration.

Satellite imagery showed Crane's terminus retreating by 6 km in the first month following the disintegration event (Fig. 2a–c), with iceberg calving known to be triggered by changes in the state and integrity of adjoining landfast sea ice (Christie et al., 2022; Arthur et al., 2021; Miles et al., 2017). Following the removal of the fast ice and melange from our model domain, an increase in extensional stresses was found in the vicinity of the terminus (Fig. 4), the magnitude of which mirrors that reported in the idealised model scenario of Surawy-Stepney et al. (2024). With these increases seen perpendicular to the terminus, it is likely that the loss of buttressing from the fast ice and melange played a key role in the initial calving response following the fast-ice disintegration (Benn et al., 2007; Walter et al., 2010). It is noted that calving may have occurred more readily at this time as the near-terminus area was characterised by susceptible, highly damaged ice (Fig. 2a), similar to observations reported by Amundson et al. (2010). A further assessment into the changes in stress regime following the calving of the floating ice would be required to better understand the progressive calving behaviour seen later in the year.

As Crane's floating ice shelf also provided buttressing against upstream ice flow (Fig. 3), it follows that with each subsequent calving event, further resistive backstress was lost, ultimately leading to the further acceleration of ice flow. With changes in grounding line flux being the most sensitive to the loss of buttressing close to the location of the grounding line (Mitcham et al., 2022), we attribute any perceived delayed response in ice-flow acceleration to be part of the evolving dynamic response of Crane following the disintegration of the fast ice.

5.3 Thickness of fast-ice and melange elements

We configured the thickness of the melange in our model using a high-resolution surface elevation profile (Howat et al., 2019), defined at 2 m resolution with an associated ±2 m error estimate per pixel. The average surface elevation of model elements downstream of the December 2021 terminus location is 4.57 m (41.66 m thickness), with minimum and maximum surface elevations and standard deviation of 0.022, 50.69 and 7.33 m respectively. If it is assumed that the maximum associated error of 2 m exists over the entire modelled area of melange, the upper limit to the average change in surface elevation in this area is 43.75 %.

The thinnest configuration of melange that we tested considered a 40 % reduction in surface elevation, which resulted in a 14.3 % increase in the calculated mean buttressing ratio compared to the reference configuration (Fig. 5b). Results from the other sensitivity experiments were typically within about 6 % of the reference configuration (Fig. 5), and the amount by which these sensitivities deviated from the reference case was uncorrelated with the percentage change in surface elevation. These results show that despite considerable uncertainties in the thickness of the fast ice and melange, the buttressing ratios calculated at the glacier terminus are not significantly impacted by these uncertainties. Though the rheology rate factor A is calculated in the inversion, it is the product of A⋅h that allows us to capture the stress regime throughout the melange.

Results from the sensitivity experiments support the conclusion that the melange buttressed the glacier prior to its disintegration despite variation being seen in the distribution and magnitude of the buttressing ratios calculated at each sample point across the terminus between sensitivity cases. Mean buttressing ratios at the glacier terminus were between 0.67 and 0.72 in seven out of the eight sensitivity cases assessed and below 0.8 in all sensitivity experiments.

6 Conclusions

We defined the stress regime in Crane Glacier and the surrounding ambient fast ice and melange through inverse modelling, allowing us to quantify the buttressing provided to Crane by a realistic melange profile for the first time.

Our results show that Crane was buttressed by the melange prior to the disintegration of multi-year landfast sea ice over the Larsen B embayment in January 2022. A mean buttressing ratio of Θ_N=0.68 was found at the interface between the glacier and ice melange. Further, we have shown that the stress distribution over regions of melange can be modelled effectively by considering the product of the rheology rate factor A optimised through the inversion of measured velocities and the input ice thickness h. Uncertainties in the thickness of the fast ice and melange therefore has little impact on the calculated buttressing ratios.

Observations showed a 3.4 % increase in flow speed at the ice front and the rapid loss of 6 km of Crane's floating ice shelf in the first month after the disintegration, whilst our model simulations found a perturbation in the near-terminus stress regime following the removal of melange from the model domain. An average increase in extensional stresses normal to the terminus of 19.2 kPa was found, which likely triggered the rapid calving of Crane's floating ice over a region which was already highly damaged. The consequent loss of large portions of the floating ice shelf led to the further reduction in resistive backstresses to the glacier, which was followed by the further acceleration of ice flow throughout 2022.

In order to gain a more comprehensive understanding of how significant the role of melange and fast ice may be in buttressing ice shelves and outlet glaciers, future research should focus on different regions where topography and flow characteristics differ from that of Crane. The methodology of this study can be applied to any area of interest if good-quality observational velocity fields of both the glacier and melange are available to constrain the model for accurate calculation.