Glaciological setting of the Queen Mary and Knox coasts , East Antarctica , over the past 60 years , and implied dynamic stability of the Shackleton system

The discovery of the deepest subglacial trough beneath the Denman Glacier, combined with high rates of basal melt at the grounding line, have caused significant concern over its vulnerability to retreat. Recent attention has therefore been focusing on understanding the governing dynamic controls, although knowledge of the wider regional context and timescales 15 over which the future responses may occur remains poor. Here we consider the whole Shackleton system, comprising of the Shackleton ice shelf, Denman Glacier and adjacent Scott, Northcliffe, Roscoe and Apfel glaciers, about which almost nothing is known. We widen the context of previously observed dynamic changes in the Denman Glacier into the wider region of the Queen Mary and Knox coasts; with a multi-decadal timeframe and an improved biannual temporal frequency of observations in the last seven years (2014-21). We integrate new satellite observations of ice structure, changes in ice front position and ice20 flow velocities to investigate changes in the system. We furthermore use the BISICLES ice sheet model to assess the sensitivity and simulate the response times of the Queen Mary and Knox coasts to hypothetical disintegration of its floating ice areas, in response to coupled ocean and atmospheric forcing. Over the 60-year period of observation, the Queen Mary and Knox coasts do not appear to have changed significantly and higher frequency observations have not revealed any significant annual or sub-annual variations in ice flow. A previously observed increase in the ice flow speed of the Denman Glacier has not continued 25 beyond 2008, and we cannot identify any related change in the surface structure of the system since then. We do, however, observe more significant change in the Scott Glacier, with an acceleration in ice flow associated with calving and progressing from the ice front along the floating tongue since early 2020. No changes in surface structure or ice flow speed are observed closer to the grounded ice. Our upper limit numerical simulations for a 400-year period are consistent with noticeable grounding line retreat in the Denman Glacier in the next two centuries if all floating ice were lost, before stabilising again in 30 the third century from now. This equates to around 6 cm of sea level rise, a small contribution when compared to other areas of East Antarctica expected to change over the same time frame. It is clear that current knowledge is insufficient to explain the https://doi.org/10.5194/tc-2021-265 Preprint. Discussion started: 6 October 2021 c © Author(s) 2021. CC BY 4.0 License.


Structure and feature mapping from optical and SAR imagery
Surface structures and features of the Queen Mary and Knox coasts were mapped from satellite imagery using standard GIS techniques (following Glasser et al. (2009)). Structural features have been mapped every 6 months from February 2015 to 100 February 2021 using freely available datasets from Landsat 8 OLI (all downloaded from https://earthexplorer.usgs.gov/) and Sentinel 2A and 2B (all downloaded from https://scihub.copernicus.eu/dhus/#/home), where cloud cover is <15%, in combination with Sentinel 1A and 1B GRD (all downloaded from https://scihub.copernicus.eu/dhus/#/home) to improve spatial and temporal coverage. To include multi-decadal changes in the extent and structure of the whole system we used several different datasets from three time periods, only choosing the datasets that covered the entire area of interest. These  (Scambos et al., 2007). Datasets were registered to the Sentinel-2 imagery as required.

110
Mapped features included (where visible); the ice-shelf or floating-glacier edges, rifts, crevasses and crevasse traces and longitudinal surface features (following Glasser et al. (2009)). Interpretation of the optical imagery was performed using multiple band combinations to provide natural colour (Landsat-1 MSS bands 7-4-3, Landsat-5 TM bands 5-2-1 and Sentinel-2 bands 4-3-2) and for all imagery standard enhancement procedures (contrast stretching and histogram equalisation) were used to improve the contrast across features. The spatial resolution of the data sets varies from 10 m to 150 m and is thus a 115 limitation on the minimum size and accuracy of the features mapped in each data set.

Feature tracking from Sentinel-1
Glacier surface velocities were derived using feature tracking between pairs of synthetic aperture radar (SAR) images acquired by the Sentinel-1 satellite (all downloaded from https://scihub.copernicus.eu/dhus/#/home). Following commonly adopted methods feature tracking uses cross-correlation to find the displacement of surface features between pairs of images, which 120 are then converted to velocities using the time delay between those images (Luckman et al., 2007). We used image patch sizes of ~ 1 km in ground range, sampled at ~ 100 m in range and azimuth. Where the time-delay between images is sufficiently short, and surface change is minimized (for instance by very cold temperatures), trackable surface features include fine-scale coherent phase patterns (speckle) and the quality of velocity map is maximized. We applied feature tracking to many image pairs and selected the best velocity map in terms of minimum noise and maximum coverage of high-quality matches. 125 Uncertainties in velocity magnitude are around 0.2 m day -1 (Benn et al., 2019). This approach allowed us to optimize the quality of the surface strain map derived from the surface velocity. The image pair chosen was 18 th to 24 th November 2017 (6day repeat between Sentinel-1B and Sentinel-1A).

BISICLES ice flow model 130
We used the BISICLES ice flow model to investigate the response of the Queen Mary and Knox coasts to hypothetical sustained disintegration of floating ice, while surface mass balance remains at present day values. A similar investigation carried out for the whole of Antarctica (Martin et al., 2019) was unable to comment on this region because the Bedmap2 bedrock (Fretwell et al., 2013) did not resolve the Denman trough. Here, we address that short coming by using ice thickness and bedrock elevation data from BedMachine v2 , which includes a deep trough beneath Denman 135 Glacier. The model simulates the flow of ice numerically, employing finite-volume discretizations of an ice thickness transport 140 and a two-dimensional vertically integrated stress balance equation, to determine the ice thickness ℎ and the velocity field . In Eqn. (1) and (2), is the rate of ice accumulation and the upper 145 surface and is the melt rate applied to the base of floating ice. ́ is the horizontal rate of strain tensor, is the ice surface elevation, is the ice density, and is the gravitational acceleration. The vertically-integrated effective 150 viscosity ℎ ̅ is expressed in terms of the ice thickness, a temperature and rate-strain dependent part ̅ A given by Glen's flow law, and a 'stiffness factor' . Temperature data are from Pattyn ( 2010). The basal friction is calculated according to a regularized Coulomb law (Joughin et al., 2019;Zoet and Iverson, 2020) In this case, | |= 300 ma -1 . In common with previous work (Cornford et al., 2015;Martin et al., 2019), the fields ( , ) and & ( , ) are found by minimizing the initial mismatch between modelled speed | | and the observed speed | |, Where & is set to 1 where speed data are available and 0 elsewhere, and is the whole domain. is minimized using the nonlinear conjugate gradient method, as in all previous realistic BISICLES applications, while the observed speed is that of Rignot et. al. ( 2011).

165
The model domain is shown in Fig. 2  The single future simulation follows Martin et al. (2019). Surface accumulation is held at present day values, while an unrealistically large melt rate, up to 1000 ma -1 in places, is applied Under this forcing, the existing ice shelves thin over the course of a few years to a thickness of 100 m, and any regions of newly formed ice shelf following grounding line retreat thin 175 in the same fashion. There is no control simulation as we are not attempting to evaluate the variation in response to plausible forcing but to quantify the upper bound of the ice stream response to ice shelf ablation given the basic physics of the model.

Ice front positions 180
Between 2015 and 2021 the floating ice front of the Shackleton Ice Shelf advanced steadily over a distance of ~ 0.3 km year -1 . Calving occurred from the western side of the ice front in 1991, and a portion of the calved ice has since remained grounded just offshore of it (Fig. 3a). Between 1962 and 2021 the shelf's central front advanced a total of 18 km with no obvious change in rate of advance (Fig 3a). The Denman Glacier's ice front advanced at a rate of 1.8 km year -1 between 2015 and 2021, with a uniform pattern of advance and no seasonal variability in rate observed (Fig. 3b). An iceberg from a calving event on the 185 Denman Glacier, hypothesized to have occurred in the late 1940s (> 70 km in length; Miles et al., 2021) is visible in 1962, roughly 100 km offshore the ice front (Fig. 3a). The Denman ice front position retreated in 1984 due to another major calving event (54 km in length; Miles et al., 2021). By 1991 the floating margin was still located 10-15 km south of the 1962 position but has since advanced ~ 61 km (Fig. 3a). Remarkably, the floating ice front of Scott Glacier has experienced more variability than that of Denman or Shackleton (Fig. 3a). Between 2015 and 2019, the front advanced at a steady rate of ~ 0.75 km year -1 190 ( Fig. 4a). Since early 2020, small scale calving has caused the ice front of the eastern half of the glacier to retreat ~ 5 km further south of the 2015 front (Fig 4b-d). The western side of the Scott Glacier ice front is in a similar position in 2021 to that of 1962 but the whole ice front was ~ 10 km further south in 2009 (Fig 3a, 4c).

Ice structure 195
Two major rift systems dominate on the Shackleton Ice Shelf, both of which extend westwards from its eastern lateral margin (labelled '1' and '2' in Fig. 1b). This margin is bordered to the east by a region of heavily fractured ice, ~ 2,300 km 2 in size ( Fig. 1b), held in place by fast ice and ~ 150 m thinner than the adjacent ice shelf body (Fretwell et al., 2013). System 1 is a maximum of ~ 15 km wide at the eastern margin and extends over 40 km into the ice shelf, narrowing and eventually terminating at a spatially extensive suture zone that originates in the leeside cavity of Masson Island ( Fig. 1b & 3a). A 200 subsidiary rift branches off to the north and connects with the ice shelf front (Fig. 3a). The geometry of system 1 has not changed significantly since 1962, although its width increased by ~ 5.3 km between 1962 and 1991 ( Fig. 3a) and in 1962, there was no clear connection between the infant subsidiary rift and a front-parallel rift, visible by 1991 (Fig. 3a). System 2 is a maximum of ~ 5 km wide at the eastern margin and extends into the shelf for 16.5 km, before branching into two rifts that trend in opposing directions, ~ 14 km and ~ 21 km in length respectively ( Fig. 1b & 3a). System 2 changed more substantially 205 than system 1, branching towards the grounding line and lengthening by 3.8 km between 2015 and 2021. In 1962 the rift is only visible as a crack, opening to a rift 2.3 km wide by 1991 and at the eastern edge, 4.5 km wide by 2017 ( Fig. 1b & 3a).
Between 1991 and 2015 its southwestern branch increased in length from ~ 10 km to ~ 16 km and in width by ~ 1 km at the ice margin ( Fig. 3a & 3a). The northern crack increased in length from ~ 13 km to ~ 16 km over the same time-period. Both systems advected with ice flow towards the ice front between 2015 and 2021, with no significant changes in geometry (Fig.  210 3a).
The surface of the Denman glacier is heavily featured with a combination of crevasses, flow lines and channel-like features.
A number of small rifts (< 6km long) are evident along the western margin, separated by fast ice from Shackleton Ice Shelf and there is no identifiable change on the length or position of these rifts relative to the ice front between 2015 and 2021 (Fig  215   3b).
The floating portion of Scott Glacier is dominated by a series of rifts striking perpendicular to the flow direction (Fig. 6) through to 2009 (Fig. 5b). A rift on the western side of Scott Glacier, initiating close to the margin with Chugunov Island, now connects with an opening rift at the ice front, detaching a portion of the front of Scott, > 24 km in length (Fig. 5c). There is little observable change in Roscoe Glacier with the exception of a rift opening in the vicinity of the margin with Shackleton Ice Shelf (Fig 6). In 2021 the rift is 15 km long and 2 km wide at the widest point and extends to within 3.2 km of the ice front, 225 a significant increase in dimensions of 5.3 km and 0.35 km, observable in 2015 when the feature terminated 8.5 km from the ice front (Fig 6). (abutting Apfel Glacier and Taylor Islands) and west (abutting Denman Glacier) of Scott Glacier (Fig. 7). The eastern margin appears as a series of small rifts and crevasses, largely perpendicular to ice flow (Fig 7a-b). Over the 7-year period there has been lengthening of the features into the ice to both the east and the west of the eastern margin, as well as opening of existing features (Fig. 7). The western shear margin of Scott Glacier is more clearly defined and has been widening into Denman Glacier in the vicinity of Chugunov Island (Miles et al., 2021). In 2015 this margin is relatively straight, in line with the 235 floating margin of Denman and ~ 1.3 km wide. Remarkably, the shear margin appears to bulge progressively into Denman Glacier and is double the width by 2021 as compared with 2015 (Fig 7c-d).

Ice flow speed
Mean ice flow speed derived at biannual temporal frequency from Sentinel-1 data between 2017 and 2021 varies across the 240 Queen Mary and Knox coasts. It ranges from ~0.2 m day -1 in the area between the grounding line and Masson Island on the Shackleton Ice Shelf to ~5 m day -1 on the floating tongue of Denman Glacier (Fig. 8a). Roscoe Glacier shows speeds of 1-2 m day -1 , higher than the surrounding ice shelf and the floating portion of Scott Glacier reaches 2-3 m day -1 . Recent changes in ice speed, derived by differencing the mean speed across the whole system between 2020 and 2019, are confined to the lower ~ 60 km of the floating tongue of Scott Glacier (Fig. 8b). Increases of > 0.2 m day -1 occur across the outer 25 km, decreasing 245 to ~ 0.08 m day -1 close to the large rift adjacent to the Taylor Islands (Fig. 8b). No acceleration (or deceleration) is observed elsewhere in the system (Fig. 8b).
Ice flow speed extracted from Sentinel-1 provided a timeseries between 2017-2021, extended back to 2002 using Measures and ITS_LIVE in locations where it is available (Rignot et al., 2017), highlighting variability through time across the system 250 (Figs 9 and 10). Point locations on Shackleton Ice Shelf vary between ~ 0.2 m day -1 at the grounding line and ~ 1 m day -1 towards the front of the floating ice, with no consistent temporal trends (Fig. 9a). Denman Glacier exhibits higher speeds, from < 2 m/day upstream of the grounding line to ~ 5 m day -1 on the floating tongue but speeds remain constant through time at each point location (Fig. 9b). Scott Glacier has a similar spatial pattern with speeds increasing from ~ 1.2 m day -1 at the grounding line to > 3 m day -1 close to the floating ice front (Fig. 9c). Speeds ~ 10 km either side of the grounding line show no change through time however, the outer 30 km of the floating ice tongue show significant acceleration from the beginning of 2020 through to May 2021. Over the 17-month period ice speeds increase ~ 30 % to 2.5 m day -1 30 km from the ice front and ~ 40 % to 3.2 m day -1 close to the front (Fig. 9c). Roscoe Glacier has similar spatial patterns to both Shackleton and Denman, with slower speeds of ~ 0.4 m day -1 at the grounding line, increasing to ~1.2 m day -1 close to the floating ice front and no significant change in speed through time (Fig. 9d). 260 The magnitude of the principal strain rate, derived from mean velocity maps, highlights the shear margins of the Denman Glacier and those of Scott and Roscoe Glaciers (Fig. 10), as well as pinning points of the Shackleton Ice Shelf (Fig. 10).
Pinning points have previously been identified at the front of the Roscoe -Shackleton shear margin (labelled (a) in Fig. 10) and upstream of rift 2 on the Shackleton Ice Shelf ((b) in Fig. 10) (Fürst et al., 2015). There is evidence of two additional 265 pinning points at the front of the Denman-Scott shear margin ((c) in Fig. 10) and at the ice margin of Shackleton Ice Shelf ((d) in Fig. 10). The latter coincides with a rise in the ocean floor (Arndt et al., 2013). Basal traction over the region varies by around one order of magnitude from ~10 kPa to ~100 kPa, with stripes of soft and hard bed appearing both in the glacier trunks and outside. The stiffness factor appears similar to other BICISLES optimizations, with weak shear margins apparent in both ice stream and ice shelf regions. The future simulation shows the response of the system to hypothetical rapid, immediate, and sustained disintegration of the floating ice. A notable and rapid retreat of Denman 275 Glacier occurs after 2150, but otherwise the overall dynamic response is modest (Fig. 12). Between 2010 and 2110, the loss of essentially all floating ice in the Denman Glacier shelf leads to around 20 km of grounding line retreat of the Denman and neighbouring Scott glaciers (Fig. 12). At the same time, flow across the grounding line discharges 80-100 Gt year -1 ice volume above flotation (Fig. 13). In the following century, the grounding line in the Denman trough retreats by around 100 km over the region of lowest elevation (Fig. 12), while discharge increases to more than 200 Gt year -1 volume above flotation in 2150 280 (Fig. 13). This period of rapid change is only temporary on a centennial scale with a cumulative sea level rise contribution of ~ 6 cm, and by the following century ice grounding line retreat has all but ceased as the bedrock shoals (Fig.12) and discharge returns to ~ 100 Gt year -1 (Fig. 13), similar to the present day.

Discussion
Over the ~ 60-year period of observation, the Queen Mary and Knox coasts have not changed significantly. More frequent 285 satellite observations in recent years have not revealed any distinct temporal patterns of change in the system. With prominent suture zones and pinning points as likely agents of stability (Kulessa et al., 2014(Kulessa et al., , 2019, the front of the Shackleton ice shelf is only slowly advancing and little change in the geometry of the main surface features occur. The flow speed of the ice shelf is regulated by the extensive suture zone downstream of Masson Island (Fig 1b, 3a) that arrests the two large rift systems, and by pinning points at the front of the Roscoe-Shackleton shear margin, on the western side of the ice shelf and to the inland side 290 of the large rifts (Fig. 10). An increase in ice flow speed was observed just upstream of the Denman Glacier groundling line between 1972-4 and 1989and, to a lesser extent, through to 2008(Miles et al., 2021). Variability in ice flow speed then became insignificant through to 2016-17 (Miles et al., 2021, their Fig. 3c), a pattern that has continued since (Fig.   8b, 9b) and accordingly we cannot identify any related change in the structure of the system (Fig. 3). Scott Glacier has received less attention until now being thinner and slower than Denman, with an overall decrease in velocity observed between 1972-4 295 and 2016-7 (Miles et al., 2021). However, this part of the system is where we observe more significant change since early 2020 (Fig 5,7,8,9). Ice flow acceleration is progressing from the calving front along the floating tongue of Scott Glacier and is particularly pronounced within the frontal ~ 60 km (Fig. 8b) where small-scale calving has been observed (Fig. 5). No changes in structure or ice flow speed are observed up flow of the large rift to the west of the Taylor Islands, and the acceleration does not currently appear to have any connection to the grounded ice (Fig. 9b). Surface meltwater features, reported to be 300 frequent around the Scott and Apfel grounding lines do not appear to be increasing in area or frequency between 2000 and 2020 (Arthur et al., 2020) and are unlikely to be contributing to the changes observed on Scott Glacier.
The upper limit scenario of forcing in our BISICLES model runs suggests that noticeable grounding line retreat occurs in the Denman Glacier over the simulated 400-year time period, briefly doubling the discharge of mass above flotation to 200 Gt 305 year -1 around the year 2150, before stabilising again with discharge returning to around 100 Gt year -1 (Fig. 13). With a cumulative contribution of ~ 6 cm to sea level rise the increased discharge is not insignificant, but small compared to possible contributions from other areas of East Antarctica over the same time frame (Martin et al., 2019). Our newly discovered and any previous reported changes in the Denman and Scott Glaciers are much smaller than those considered here as an upper limit of forcing in BISICLES. Any real future dynamic changes and contributions to sea level rise will therefore very likely be less 310 than those illustrated here.
One possible interpretation of the model output is that Queen Mary and Knox coasts are relatively stable and insensitive to reasonable forcing in next 400 years, and that ice loss from the Shackleton system poses no imminent threat to the Aurora and Wilkes subglacial basins. This interpretation is reinforced in that implied rates of basal melting of the Denman ice tongue, 315 albeit high on a continental scale, are likely much lower than those required to precipitate full disintegration of the Denman Glacier floating tongue on the short timescales simulated here. However, our BISICLES simulation may suffer from poor boundary constraints due to unknown or poorly known subglacial substrate, basal hydrology, geothermal heat flux, ice properties, oceanographic conditions, and bathymetry.

320
On the nearby Totten Glacier inferences from combined geophysical exploration and numerical modelling are consistent with areas of high basal melt rate coinciding with significant grounding line retreat, possibly linked to channelized subglacial meltwater discharge (Dow et al., 2020). By analogy the deep trough beneath the Denman Glacier is also likely to favour vigorous channelization of the subglacial meltwater system close to the grounding line. Basal water and ice sheet behaviour are also affected by geothermal heat flow. Regional values of ~ 55-85 mW m -2 are interpreted from magnetic (Martos et al., 325 2017) and seismic  data, similar to the average global value of ~ 67 mW m -2 for continental regions (Lucazeau, 2019). Multivariate analysis of Antarctic and global geophysical and geological datasets is consistent with elevated geothermal heat flow, > 70-80 mW m −2 , west of the Denman region near the Gaussberg Volcano (Wilhelm II Coast) and in the Knox interior (Stål et al., 2021). These high heat flow anomalies remain poorly resolved but are likely driven by the geological evolution including volcanism, neotectonics and variation in crustal heat production, thermal conductivity and topography 330 (Stål et al., 2021). The geometry and infill of the subglacial Knox Sedimentary Basin (Maritati et al., 2016) likely imparts important controls on the Denman region heat flow distribution; basement rift flanks could be channelized regions of elevated heat (Willcocks et al., 2021), and groundwater could also be an important mechanism promoting intra-basin advection of heat (Siegert et al., 2018;Kulessa et al., 2019). A series of very large 200-300 km wide deep-seated granitic bodies have been inferred within the interior of Wilkes-Queen Mary lands, 600-700 km to the south (Aitken et al., 2014); if these subglacial 335 batholiths are strongly heat producing, then they could contribute to hot spots in geothermal heat flow as has been modelled along the Prydz Bay coast (Carson et al., 2014).

Conclusions
We conclude that over the 60-year period of observation, the Queen Mary and Knox coasts do not appear to have changed 340 significantly and higher frequency observations have not revealed any significant annual or sub-annual variations in ice flow.
The velocity changes on the Denman Glacier recently described (Miles et al., 2021;Rignot et al., 2019) appeared to be shortlived events focused on the glacier itself. We do observe more significant change in the Scott Glacier, with an acceleration in ice flow likely triggered by calving and progressing from the ice front along the floating tongue since early 2020. These shortterm changes in the flow speed and structure of Scott glacier have not yet had any noticeable impact on ice dynamics close to 345 the grounding line. Nonetheless, the BISICLES ice flow model provides us with initial insights into the extent to which reasonable upper limit changes in oceanographic or surface forcings can precipitate major grounding line retreat on centennial timescales. In the absence of better understanding of future changes in the oceanographic and atmospheric conditions in the Queen Mary and Knox coasts we consider the imminent break-up of the floating area of the whole Shackleton system as a reasonable upper limit of forcing, although real changes to this area will likely be rather less dramatic. Given the potential 350 vulnerability of the system to accelerating retreat better data recording the glaciological, oceanographic, and geological conditions in the Queen Mary and Knox coasts are urgently required to improve the certainty of numerical model predictions.
Current knowledge is insufficient to explain the observed spatial and temporal changes in the dynamic behaviour of the grounded and floating sections in the Shackleton system. With access to these remote coastal regions a major challenge, coordinated internationally collaborative efforts are required to quantify how much the Queen Mary and Knox coastal region 355 is likely contribute to sea level rise in the coming centuries.