Observed meltwater-induced flexure and fracture at a doline on George VI Ice Shelf, Antarctica

Abstract

Global Navigation Satellite System (GNSS) observations and ground-based timelapse photography obtained over the record-high 2019/2020 melt season are combined to characterise the flexure and fracture behaviour of a previously formed doline on George VI Ice Shelf, Antarctica. The GNSS timeseries shows a downward vertical displacement of the doline centre with respect to the doline rim of ~60 cm in response to loading from a central meltwater lake. The GNSS data also show a tens-of-days episode of rapid-onset, exponentially decaying horizontal displacement, where the horizontal distance between the doline rim and its centre increases by ~70 cm. We interpret this event as the initiation and/or widening of a fracture, aided by stress perturbations associated with meltwater loading in the doline basin. Viscous flexure modelling indicates that the meltwater loading generates tensile surface stresses exceeding 75 kPa. This, together with our timelapse photos of circular fractures around the doline, suggests the first such documentation of meltwater-loading-induced ‘ring fracture’ formation on an ice shelf, equivalent to the fracture type proposed as part of the chain-reaction lake drainage process involved in the 2002 breakup of the Larsen B Ice Shelf.

1. Introduction

Summer meltwater ponding has been observed on several Antarctic ice shelves for many decades (David and Priestley, 1909; Stephenson and Fleming, 1940), and is widespread around the continent (e.g. Kingslake and others, 2017). Ponding is implicated in the processes of surface meltwater-induced flexure and hydrofracture (Scambos and others, 2009; Banwell and others, 2013, 2019) making it a potential source of ice-shelf instability and break-up (Robel and Banwell, 2019; Lai and others, 2020). With future surface melting across ice shelves predicted to increase (Gilbert and Kittel, 2021; Kittel and others, 2021), resulting in increased firn-air-content depletion (Dunmire and others, 2024) and ice-shelf break up leading to grounded ice-flow acceleration and sea-level rise (Scambos and others, 2004; Wuite and others, 2015), it is important to document and better understand the processes of surface meltwater-induced flexure and hydrofracture in response to surface-meltwater loading and unloading.

While surface meltwater ponding and drainage can be sensed remotely with airborne photographs (e.g. Bell and others, 2017) and satellite data (e.g. Arthur and others, 2020; Dirscherl and others, 2021; Dell and others, 2020, 2022; Banwell and others, 2023), ice-shelf deformation changes in response to meltwater-driven flexure or fracture are rarely detectable by remote sensing, despite the potential of InSAR in this regard (Li and others, 2021). A recent exception is presented by Warner and others (2021), who document, using repeat digital elevation models (DEMs) and ICESat-2 satellite laser-altimetry track data, the creation of a doline on the Amery Ice Shelf in response to the sudden (<3 d) draining of a meltwater lake by hydrofracture. This event involved a ~11 km² area of local surface lowering (averaging −24 m and with a peak value of −80 m), at the centre of a ~60 km² region that underwent uplift (reaching 36 m adjacent to the central depression) due to the ice shelf's flexural hydrostatic rebound following the drainage of meltwater into the ocean below.

The present study builds on a recent field campaign that provided the first in-situ evidence of ice-shelf flexure in response to the filling and drainage of a surface lake on the McMurdo Ice Shelf (Banwell and others, 2019). Using Global Navigation Satellite System (GNSS) observations combined with timelapse camera imagery and pressure transducer-derived lake depth data, they showed that the lake drained by overflow and channel incision, rather than by hydrofracture. The magnitude of vertical ice-shelf deflection reached ~1 m at the lake centre, declining to zero at ~500 m from the lake centre.

In the present study, we analyse in-situ observations of meltwater-driven ice-shelf flexure and fracture at a doline on north George VI Ice Shelf (Fig. 1). Using continuous observations from static GNSS stations located at the doline centre and rim crest, and photographs from a timelapse camera looking into the doline basin, we document the ice-shelf's flexure in response to meltwater ponding in the doline's centre. We additionally document what we interpret to be the initiation and subsequent horizontal widening of a ring fracture (Banwell and others, 2013) between the doline's rim and centre. This is the first real-time in-situ observation of a fracture on an ice shelf initiated by meltwater-loading-induced flexure. Additionally, through combining these observations with idealised viscous-flexure modelling of the ice shelf, our analysis shows that meltwater loading-induced stresses likely play a key role in triggering ice-shelf fracture.

Figure 1. Field area. (a) Overview of the Antarctic Peninsula showing George VI Ice Shelf within the George VI Sound between Alexander Island to the West and Palmer Land to the East (map adapted from Banwell and others (2021), their Fig. 1a). (b) Overview of the ice shelf from Antarctic REMA DEM (Howat and others, 2019). The locations of our AWS and the BAS Fossil Bluff AWS (‘FB AWS’ ) are indicated by red stars. (c) Close-up of the REMA DEM over our study site (taken from the REMA Explorer web browser). (d) WorldView-2 image (18 January 2020) of the doline with a lake in its centre. The small orange dots indicate the locations of three (of our four original) GNSS stations originally installed; data from just the two GNSS stations (GPS01 and GPS02) closest to the doline centre are used in this study (the third, GPS03, also shown, and fourth, GPS04, not shown, were flooded with meltwater). Red boxes in panels A, B and C show the regions seen in close-up view in panels B, C and D, respectively.

2. Background on ice-shelf dolines

Dolines on ice shelves are distinctive, but relatively rare, drained-lake features (Stephenson and Fleming, 1940; Reynolds, 1983; Moore, 1993; Bindschadler and others, 2002). Analogous to a sink hole in karst terrain, a doline forms when water from a lake on, or just below, the surface of the ice shelf drains through the ice shelf into the ocean below (Mellor, 1960; Mellor and McKinnon, 1960; Bindschadler and others, 2002; Banwell and MacAyeal, 2015; MacAyeal and others, 2015; Warner and others, 2021). The void created by the absent water forms a depression in the ice-shelf surface that is bowl shaped, often circular in plan view, and is sometimes filled with rough pieces of the former lake's superimposed ice lid. Dolines are thought to form rapidly, over hours to days, as was witnessed recently from satellite laser altimetry and imagery (Warner and others, 2021). Immediately following doline formation, the ice shelf responds to the missing weight of the drained meltwater by flexing upward in the centre of the doline, a process referred to as ‘hydrostatic rebound’. Due to the stiffness of the ice, there is a regional uplift surrounding the original region of mass unloading, in the form of an uplifted rim (Bindschadler and others, 2002), whose spatial extent is initially determined by the ice shelf's thickness and elastic properties (MacAyeal and Sergienko, 2013). The region of uplift is surrounded by an inverted forebulge (i.e. a surface ‘moat’), where the ice shelf surrounding the doline is flexed down to a lower elevation than the surrounding ice shelf (Banwell and others, 2013; MacAyeal and others, 2015). This inverted forebulge attracts summer meltwater, leading to the formation of lakes that have an annular or ring-like geometry, the appearance of which is commonly used to identify the presence of the doline visually (Fig. 1d; see also Fig. 1 of MacAyeal and Sergienko (2013)).

Once formed by the initial lake drainage, dolines can persist for decades, and their surface topography is observed to slowly evolve in response to secondary ice-shelf processes associated with viscous ice relaxation, ice flow, accumulation, ablation and meltwater movement and loading within and near to the doline. For example, a doline on the George VI Ice Shelf near the subsequently established British Antarctic Survey (BAS) Fossil Bluff station was initially observed in 1936 by members of the British Graham Land Expedition (Stephenson and Fleming, 1940). It was still visible in the 1980s (e.g. Reynolds, 1983), and it remains visible in satellite imagery now. For the doline that is the focus of this study, informally named ‘South Doline’ (Fig. 1, location: 71.5109 S, 67.70145 W), its estimated age since initial formation in response to a surface lake drainage event near the grounding line is a maximum of 50 years. This estimate is based on the fact that the doline must have formed on the floating ice shelf to have its characteristic uplifted rim geometry with a central depression, and, based on MEaSUREs v2 ice flow velocities (Rignot and others, 2017), that 50 years of advection has transpired between the ice shelf's eastern grounding line and the doline's present position. Our inspection of Landsat satellite optical imagery reveals the presence of the doline continuously, in close to its current geometry, since at least 1997. However, although South Doline is mature in terms of its age, its surface topography has undergone decades of evolution since its initial formation. This has occurred in response to both its westerly advection further onto the ice shelf with ice flow, resulting in South Doline's current elliptical plan geometry (Reynolds, 1983), as well as to changes in surface melt and hydrology, as focused upon in the current study.

We had two major reasons for deploying ground-based field instrumentation around South Doline. First, as ice shelves are typically relatively flat, it is difficult to predict the location and plan-form geometry of common surface meltwater lakes from season to season. This is especially the case in the spring to early summer period, when it is possible to conduct field work on ice shelves that undergo high summer melt rates, like George VI Ice Shelf (Banwell and others, 2021). However, to measure flexure in a manner that captures the key length scales associated with ice thickness and viscoelastic strength (e.g. MacAyeal and Sergienko, 2013), GNSS satellite receivers need to be located in a transect extending from a known lake centre. This attention to location was attempted with fortunate success in a field campaign on the McMurdo Ice Shelf (Banwell and others, 2019). A lake that forms in the basin of a doline (Fig. 2b) has a much more predictable geometry, and this motivated our field effort on South Doline.

Figure 2. South Doline on north George VI Ice Shelf. (a) WorldView-2 image of South Doline (18 January 2020) with instrument locations noted. Dotted lines denote the field of view of the timelapse camera. A meltwater lake (with several irregularly shaped lobes) is observed in the doline's centre. (b) Interpretation of the image in panel A showing names for various morphological features, as well as the locations of a lake-surface, height-limiting moulin (observed in timelapse imagery, indicated by a red X) and the meltwater catchment zone (i.e. the ‘doline ramp’ ) for the doline basin. (c) Photograph looking north from the location of GPS02 on the south rim of the doline taken by the field team in November 2019 after deployment of the camera. GPS01 (position indicated by the red line) and the majority of the doline basin are in the camera's field of view. (Camera viewshed shown as dotted lines in panel A.) Two points labelled ‘image stabilization points’ served as fiducial references tied to the mountainous landscape and were used to process images in a manner to remove camera wobble and pointing variations associated with wind-driven vibration of the camera mounting pole (see Methods). (d) Photograph looking South from the location of GPS01 in the doline basin taken in November 2021. GPS02 and the timelapse camera are indicated (red line) on the South rim of the doline.

The second reason for deploying our instrumentation around South Doline is that a surface lake in a doline is self-contained, as the uplifted doline rim prevents meltwater from elsewhere on the ice shelf from entering the central doline basin (Fig. 2b). Therefore, this study site provides a strict constraint on the doline's lake volume as a function of surface melt rates within a well-defined catchment, which provides an important simplification for the application of our model of meltwater-induced ice-shelf flexure.

3. Field instrumentation

3.1 Placement

Originally, we deployed four GNSS receivers in a transect out from the doline centre (following the methods of Banwell and others (2019)), an automatic weather station (AWS) and a 10 m subsurface-temperature thermistor string within the doline basin, a timelapse camera system on the doline rim, and several autonomous water-pressure sensors within and outside the doline (Figs 2a, c, d). Instrument deployment was completed during November 2019, and although annual data collection and instrument maintenance was intended, the instruments were left unattended for two years due to Covid-19. When revisited in November 2021, only two of the GNSS receivers (GPS01, near the doline centre, and GPS02, on the south side of the doline rim, Figs 2a, c, d) and the timelapse camera system (Fig. 2c) had survived and collected useful data. All other instrumentation at this site was lost or destroyed due to flooding from meltwater produced during the previous two full melt seasons, the first of which (2019/2020) was a record-high melt year for north George VI Ice Shelf (Banwell and others, 2021). In November 2019, we had also deployed an additional, AWS ~40 km north of South Doline (Fig. 1b; location: 71.14409 S, 67.69971 W), which survived the 2019/2020 melt season, and from which we use 2 m air temperature as well as in-ice temperatures (<10 m depth) in the current study.

The two GNSS systems had Trimble NetR5 receivers with Trimble Geodetic Zephyr antennas, powered by a battery with photovoltaic charging system contained in a Pelican case. Both systems collected data at 15 s intervals throughout the 2019/2020 melt season until the reduction of daylight (~1 April 2020) lowered photovoltaic charging to the point where the receivers only operated intermittently.

The timelapse camera provided an oblique view covering most of South Doline's basin and its opposing ramp and rim (Figs 2a, b), as well as the snow surface on the ramp in the foreground near the base of the camera (Fig. 2c). The camera system, commercially available through HarborTronics Inc., included a Canon model T7 (2000D) with a 18–55 mm f/3.5–5.6 lens (23–63$^\circ$ HAOV) powered by a deep-discharge absorbent glass mat lead/acid battery charged by a photovoltaic panel, and protected by a weather-proof enclosure with a glass aperture. The camera took photographs at 30 min intervals from its deployment in mid-November 2019 through the end of March 2020.

Both GNSS receiver antennas and the camera were mounted on aluminium poles that were drilled at least 2.5 m into the ice-shelf surface in November 2019 (Figs 2c, d). Thus, the camera and antennas were fixed within the material frame of reference of the surface ice. These three poles remained stable, and did not tilt or sink into the surface ice throughout the entire 2019/2020 melt season. This stability was determined from various considerations including (1) the pole for GPS01 was constantly recorded in photographs and remained vertical, (2) the camera field of view did not change appreciably through the 2019/2020 melt season (except for temporary vibrations due to wind), (3) although not in the field of view of the camera, the pole on which GPS02 was mounted was located in the same area as the camera and would have experienced the same minor effects of surface ablation as the pole for the camera, and (4) as we will show, the individual GNSS station data demonstrate that the record of relative horizontal and vertical displacements between GPS01 and GPS02 could not be explained as a result of tilt of the pole on which GPS02 was mounted.

3.2 GNSS processing

Following methods used by Banwell and others (2019), GNSS data were collected at a 15 s resolution and processed to yield positions every 5 min using TRACK (Chen, 1998), the kinematic mode of the GAMIT (v.10.71) software (Herring and others, 2018). For a base station located on bedrock, we used FOS1, the GNSS station operated near to the BAS Fossil Bluff Research Station (Koulali and others, 2022), ~30 km from South Doline. However, the FOS1 base station was not operational until 15 December 2019, thus two to three weeks of data recorded by GPS01 and GPS02 prior to that time were not used in the analysis here. We used default values for TRACK processing parameters, which included specification of minimum satellite elevation (15$^\circ$); horizontal and vertical noise tolerance (1  cm); the satellite constellation(s) and number of included satellites (33, which was the functioning GPS constellation operated by the US during the period the GNSS receivers were operating); and the signals to process (carrier signals L1 and L2). For specification of GPS constellation orbits and clock errors, we used the Massachusetts Institute of Technology compilation of final satellite precise ephemerides (.sp3 files) provided by the US National Geodetic Service. TRACK outputs reported standard deviation of horizontal and vertical positions on the 5 min time-series of ~2 cm and ~5 cm, respectively. Vertical-elevation time series were provided by TRACK relative to the geoid, which in the region of the George VI Ice Shelf was taken to be 7.86 m.

To restrict attention to only long-period monotonic vertical displacements associated with meltwater-load ice-shelf flexure, the average tidal signal from the two GNSS stations was subtracted from the 5 min raw elevation data from TRACK. The tidal displacement, which is common to both GNSS stations, was estimated by fitting 34 tidal constituents ranging in period from 3.1 h to 27.5 d (including compound tides, tidal harmonics, semidiurnal, diurnal and long-period tides) to the vertical-elevation time series using least squares following the method of Pawlowicz and others (2002) (see also, Banwell and others, 2019). Following de-tiding, these GNSS data were further corrected by subtracting the inverse barometer effect (IBE) deduced from the barometric pressure observed at the BAS Fossil Bluff AWS (Fig. 1b, location: 71.329 S, 68.267 W, 66 m a.s.l.) using a sensitivity of 97.27 hPa m⁻¹.

The 5 min timeseries of horizontal and vertical position were cleaned by visual inspection of outliers and obvious short-term errors associated with multipath signal interference. Segments of both the vertical and horizontal time-series of <30 min were discarded where the vertical elevation was >10 cm different from the trend associated with daily tidal variation. The data used to determine long-term (multiple day) vertical elevation and horizontal separation changes between the doline centre and rim were subsequently smoothed using a 5 d running mean (i.e. data shown in Figs 6d, 7c).

3.3 Timelapse imagery processing

Initial processing of camera photographs was undertaken to remove mild effects of random wind-induced mounting-pole vibration and swaying. These movements were corrected by re-registering each photo by hand to a fixed reference frame based on surrounding mountain features using the MatLab routines cpselect and fitgeotrans. Two registration points (Fig. 2c) easily recognised in the distant mountain landscape were used to perform the re-registration in a manner that corrected for camera wobble due to wind vibration.

Secondary processing of camera photographs was done to obtain a quantitative timeseries of the total meltwater lake area through the 2019/2020 melt season. Using photographs taken at local noon on each day from mid-December 2019 to late April 2020, the images were orthorectified to yield a view of the doline basin as if photographed from directly overhead. This was accomplished by using an affine transformation in the MatLab fitgeotrans set up to match the WorldView-2 satellite image (1.84 m pixel resolution) taken on 18 January 2020 (Fig. 2a). To orthorectify the images, between 10 and 15 registration points were identified in each photograph (depending on visibility conditions) using the features in the doline basin, such as the prominent edges of the water-covered areas and surrounding streams. A sample of this orthorectified photography is shown in Figures 3a,b.

Figure 3. Image orthorectification of meltwater areas in the doline basin. (a) Oblique timelapse photo from local noon on 17 January 2020. (b) Orthorectified image of the doline basin created from the oblique photo in panel A. (c) Digitised meltwater lake areas (red outline).

The meltwater lake features visible in the orthorectified images were identified subjectively by eye, and their boundaries were digitised by hand. An example of this ad hoc digitization is shown in Figure 3c. The areas of the polygons enclosing the area covered by meltwater were then determined using the MatLab function polyarea and summed (if more than one lake boundary was needed) to produce an estimate of the total lake area. Not all orthorectified images were analysed in this manner, because digitisation by eye and hand was often found to be too difficult as a result of poor visibility or sun glint in the image. Another difficulty was that interpretation of meltwater features was at times inconsistent, particularly when there was superimposed ice during the period from early February 2020. It is thus difficult to estimate the uncertainty of lake area derived from the imagery. Repeated, but independent classifications of the same image led to differences in lake area of about 5%, and day-to-day variations in lake area during periods when superimposed ice was present was about 10%. This informal error analysis suggests that lake area derived from the imagery is accurate to ~15%.

4. Results: timelapse camera imagery

Samples of the timelapse photos are shown in Figure 4, and a timelapse video composed from the processed imagery is included in Supplementary Materials (Video S1). The timeseries of lake area data calculated from our sequence of orthorectified images is shown in Figure 5b.

Figure 4. Timelapse camera photos taken from the South rim of the doline (panels A–C taken at 12:00 local time, and corrected for camera wobble and vibration using fiducial points indicated in Fig. 2c. Panel D taken at 16:00 local time, although not corrected for camera wobble and vibration). The red line in all panels denotes the location (just above the centre of the line) of GPS01. The sequence of images shows the doline (a) prior to the onset of melting (27 November 2019), (b) when the lake had its maximum observed area (9 January 2020), (c) after a moulin visible at the base of the ramp on the North side of the doline basin has limited the lake height (photo is from 26 January 2020), and (d) at the end of the melt season (17 March 2020). Insets in panels B and C display cropped close-ups of the area with the moulin. Ring-fractures (i.e. ‘crevasses’ ) which align with the curve of the doline rim are indicated in panel C. In panels A, C and D, the foreground in the photograph is dry firn or snow. In panel B, the foreground is slush. A timelapse video of all noon-time photos (1 image per day) is provided as Video S1.

Figure 5. (a) Mean daily air temperature at 2 m above the surface (red line) from the AWS located 40 km northwest of South Doline (see Fig. 1b for ‘AWS’ location), and cumulative positive-degree days (PDD; black dots) calculated from these same air temperature data. (b) Area of surface ponding in the doline basin through the 2019/2020 melt season (blue dots/line with black error bars), calculated from orthorectified timelapse photos. As a point of reference, we estimate the total area of the doline basin floor to be ~6 × 10⁴ m².

The timelapse imagery record shows a cold, snow-covered surface from mid-November to about 3 December 2019, at which point melting was seen on the snow surface in the camera foreground, and small surface meltwater ponds began to replace the snow surface within parts of the doline basin (Figs 4a, b). Meltwater continued to accumulate in the basin until about 9 January 2020 (Figs 4b, 5b) at which time a moulin opened up at the base of the doline ramp on the north side of the lake basin (Fig. 4c), which subsequently limited the height to which the meltwater could rise within the basin (see insets in Figs 4b, c). The moulin appears to have acted in a manner analogous to a ‘safety drain’ in a bathtub that prevents over filling. After 9 January, lake area declined rapidly for 2–3 d, then fairly steadily until late March, with minor fluctuations in the rate of decline, including a period of slight enlargement in mid to late February. The lake area continued to decrease even after the water level dropped below the observed surface elevation of the moulin (Fig. 5), which may be due to vertical incision of an outflow stream flowing into the moulin.

From around late December 2019, ring fractures (visible in Fig. 3a on 17 January 2020, and in Fig. 4c on 26 January 2020) are visible in the timelapse imagery around the lake basin (Video S1) as well on the north side of the doline ramp in WorldView-2 imagery (e.g. Fig. 2a, dated 18 January 2020). The lake surface never reaches the level of these specific observed ring fractures, so they do not contribute to the drainage of the lake. However, we hypothesize that these ring fractures open as a result of stress perturbation associated with increased meltwater loading in the doline's central basin.

It is notable that although the total lake area decreased from 9 January 2020 (Fig. 5b), mean daily air temperatures at our AWS continued to fluctuate around 0 °C until early March 2020, resulting in increasing cumulative positive-degree days (PDDs) through the melt season (Fig. 5a). It is likely, therefore, that meltwater continued to be produced, and to flow from the doline ramp into the doline basin, until mid-March 2020.

5. Results: GNSS derived displacements

5.1 Vertical motion

The vertical-elevation data for the doline rim (GPS02) and basin (GPS01) between 15 December 2019 and 1 April 2020 are displayed in Figure 6. The raw (5 min) elevation data relative to the geoid shows that the overall mean elevation of the rim is ~25 m higher than that of the basin, and reveals the typical mixed diurnal and semidiurnal tidal variation of the order ~1.5 m, but also longer-term cycles (Figs 6a, b).

Figure 6. Vertical elevation relative to the local WGS84 geoid (7.86 m) of the GNSS antennae sited on (a) the doline rim (GPS02) and (b) in the doline basin (GPS01) (when displaying GNSS position data, we use the terms ‘rim’ and ‘basin’ to signify the GPS02 and GPS01 stations). Segments of the timeseries in panels A and B where data appear to be missing correspond with times when either the basin, the rim or both GNSS stations were judged to be outliers or were in gross disagreement with the local tidal trend. (c) Vertical elevation of the doline rim (blue) and basin (red) corrected for tide and inverse barometer effect, and plotted relative to zero metres elevation. Solid black lines represent 5 d running means of the relative elevation. Scatter of the blue and red points (the 5 -min data timeseries) around the 5 d running-mean curves is ~±10 cm, which we interpret as 1σ errors for the vertical elevation data. (d) Elevation difference between the basin and rim as a function of time (i.e. difference between the red and blue timeseries shown in panel C). Negative values of relative elevation indicate that the basin has dropped relative to the rim.

The long-term tide- and IBE-corrected residual elevation (relative to the mean first day elevation for each GNSS station) displays a slow increase of elevation for both stations from about the beginning of the timeseries (15 December 2019) to about mid-February 2020 (Fig. 6c). This general increase in elevation of ~0.8 cm per day is consistent with the ice shelf moving upward as a result of surface ablation at a rate of ~0.9 cm per day. The GNSS station antennae are fixed in the reference frame of the ice; so this means that when 0.9 cm of ice is ablated from the surface, the remaining ice must move upward by ${\rho _{{\rm i}}}/{\rho _{{\rm sw}}}\times 0.9 \ {\rm cm} = 0.8$ cm to maintain hydrostatic equilibrium, where ρ _i and ρ _sw are the densities of ice and salt water, respectively (Table 1).

Table 1. Parameters used in model

After about mid-February 2020, the long-term residual vertical elevations of both the rim and basin drop until late March, but the elevation of the basin drops at a greater rate than that of the rim (Fig. 6c). This observation is also apparent in Figure 6d, which shows the change in the basin elevation relative to the rim (calculated by subtracting the relative elevations of the two GNSS plotted in Fig. 6c). Over the 3.5 months of data collection, we observe an overall sinking of the basin relative to the rim of about 60 cm; first at a gentle rate with some fluctuations until mid-February, then at a faster rate until late March (Fig. 6d). We suggest that the overall sinking of the doline basin relative to the doline rim is due to surface ablation (i.e. mass removal) on the doline's ramp and corresponding meltwater accumulation (mass loading) in the doline's basin. Furthermore, we suggest that variations in the rate of differential vertical displacement may be due to temporal variations in the rate of surface ablation on the doline ramp, changes in the rate at which meltwater flows down the ramp to the basin, and/or the effects of ocean currents beneath the ice shelf. We also hypothesize that short-term (order days) reversals in the overall sinking of the basin relative to the rim observed in early January and early February may represent basin unloading as water leaves the lake basin via a moulin or fractures to the ocean below.

As mentioned previously, our timelapse photo record shows a maximum total lake area in the doline basin on around 9 January 2020, with steady lake-area decrease thereafter (Fig. 5b and Video S1). We suggest that the observed continued downward vertical displacement (Fig. 6d), following the maximum observed lake area, is due to continued viscous-flexure of the doline basin in response to the loading from ponded meltwater still in it. However, it is interesting that the rate at which the central basin continues to lower stays relatively steady, despite the observed decrease in total lake area. We offer four possible explanations for this: (1) preferential draining by the moulin (mentioned above) of a relatively shallow region of the lake, which was not fully connected to other deeper regions of the lake, meaning that the total volume of meltwater in the lake may have remained constant, or even increased, after 9 January, despite the observed decrease in the total lake area; (2) filling of a meltwater feature not visible at the surface, such as a buried lake (Dunmire and others, 2020) connected to the shrinking surface lake; (3) downwarping of the thin ice-shelf under the doline basin floor due to the weight of the surface (and potentially subsurface) meltwater, due to continued melt and runoff from the doline ramp into the lake basin until early to mid-March (Fig. 4a); and (4) enhanced lake-bottom ablation due to the lower albedo of water relative to that of ice (Tedesco and others, 2012).

Here we expand upon the four possible explanations listed above that may explain the discrepancy between the lake area and vertical elevation change timeseries. For possible explanation (1), unfortunately we do not have lake volume observations as the two pressure sensors that we deployed in the lake basin were irretrievable after our missed field season due to Covid-19, and cloud-free, high-spatial resolution optical satellite imagery of the doline basin is lacking. For possible explanation (2), we envisage that such a subsurface meltwater feature could form either in an open void within impermeable ice (Dunmire and others, 2020), and/or within a firn aquifer (Montgomery and others, 2020), whereby an increasing volume of meltwater input through the melt season would exploit the surrounding unsaturated firn. For possible explanation (3), flexure-induced lake deepening would control the relative vertical displacement of the doline basin and the foot of the doline ramp where the lake-surface height limiting moulin is located. The estimated ice shelf thickness in the doline basin (40 m) is roughly 1/6 the thickness of the surrounding ice shelf, hence it is likely for the vertical displacement between the doline rim and basin to be concentrated by flexure within the basin. If so, then the lake volume would continue to grow after mid-February as the lake bottom sinks relative to the lake-height limiting moulin on the flank of the doline ramp. This sinking would thus be some sizable fraction (e.g. 3/4) of the observed overall vertical displacement (~60 cm) between the doline rim and basin. For the final possible explanation (4), which we think will be least significant, we note that enhanced lake-bottom ablation does not change the mass loading of the lake, as the ice that melts to form meltwater remains in the lake. However, lake-bottom ablation does make room (e.g. 1/10 extra volume of the melted ice) for more meltwater from elsewhere in the doline basin (i.e. the doline ramp) to flow into the central lake, because the volume of meltwater produced at the lake bottom occupies a smaller volume than the volume of the original ice that melted (as ice has a lower density than water).

5.2 Horizontal motion

The horizontal distance, as a function of time through the 2019/2020 melt season, between the doline rim and basin is shown in Figure 7. The mean horizontal velocity of the basin and rim GNSS stations over the observation period is 11.7 cm per day. Between 15 December 2019 and 25 January 2020, the initial 329.75 m distance between the two stations increased monotonically with time by ~5 cm. Starting on 25 January 2020, the horizontal separation suddenly and rapidly increased by ~30 cm in just 4 d (Fig. 7a, black box), including ~10 cm of change on 28 January 2020. This rapid increase in horizontal separation is a result of the GNSS antenna on the rim moving in a southwesterly direction relative to the GNSS antenna in the basin. Thus, the overall trajectory of the rim station departs from the straight-line trajectory of the basin station (Fig. 7b). We interpret this event as the initiation and/or opening of a fracture between the doline rim and basin (such as those fractures observed in our timelapse imagery on the doline ramp; Fig. 3c and Video S1). On 31 January 2020, the rate of horizontal separation increases between the two GNSS slowed, however it still caused the separation distance to increase another 40 cm over a period of ~55 d, ending in late March, when the separation rate between the two GNSS returned to a similar value to that recorded at the start of our observations.

Figure 7. (a) Horizontal trajectories of the doline basin (GPS01) and rim (GPS02) GNSS antennae (red and blue, respectively) plotted on a WGS84/Antarctic Polar Stereographic projection (axes are labelled in metres of the projected coordinates relative to the starting location of the basin antenna). To compare the trajectories, the rim trajectory was displaced to overlay the basin trajectory for observations on 16 December 2019. A curious event begins on 25 January 2020 (indicated by the black box), where the rim trajectory departs from the basin trajectory by increasing its distance in the southwesterly direction by ~70 cm over ~60 d, with ~35 cm displacement in the first ~4 d. We interpret this to be the initiation and/or opening of a fracture. The true azimuth from the rim antenna to the basin antenna is ~5 degrees East of North. (b) Close up of the two trajectories corresponding to the black box shown in panel A. Trajectory positions are coloured by the relative height of the tide measured by the vertical elevation data shown in Figures 6a,b. The displacement begins on 25 January and proceeds through ~29 January with a curious pattern. Periods of low tide (blue) correspond to when the rim trajectory displaces at approximately a right angle to the trend of both the basin and rim trajectories. Periods of high tide (red) correspond to when the two trajectories are parallel in the direction of the basin's undisturbed trajectory. (c) Change in the horizontal distance between the doline basin and rim GNSS antennae as a function of time. The solid black line represents a 5 d running mean and semi-transparent black dots show error in the horizontal data.

A notable and curious feature of the increase in horizontal separation between 25 and 29 January 2020 is that the trajectory of the doline rim GNSS (GPS02) ‘zig-zags’ relative to the local tidal elevation (as deduced from our tidal analysis described in the Methods) (Fig. 7b). Prior to the sudden increase in the doline rim and centre's horizontal separation starting on 25 January, the rim's trajectory is parallel to that of the basin (GPS01), which is when the tide elevation is high (red colour in Fig. 7b). But, when the tide elevation is low (blue colour in Fig. 7b), the rim's trajectory takes an abrupt right-angle turn and proceeds in a direction that is perpendicular to the basin's long-term trajectory. As mentioned above, we interpret this separation as a result of fracture initiation and/or opening between the rim and basin. That this fracture opens only during one phase of the tide (when the tide is low) implies that tidal currents beneath the ice shelf may contribute to the driving of fracture initiation and/or opening, in addition to the meltwater loading in the doline basin.

Comparison of the timeseries of horizontal separation between the rim and basin (Fig. 7c) with vertical elevation difference between the rim and basin (Fig. 6d) shows that there is no clear corresponding event in the vertical elevation timeseries beginning on 25 January 2020. The lack of a clear correspondence between the horizontal and vertical components of relative motion confirms that the sudden horizontal movement starting on 25 January was not due to a disturbance (e.g. tilt) of either of the two GNSS antenna poles. Thus, we interpret the horizontal separation change occurring after 25 January to be associated with an impulsive onset of ice-shelf horizontal strain between the doline basin and rim, i.e. indicative of the initiation and/or opening of a fracture.

6. An exploratory model of doline-related ice-shelf movement

The field observations presented above open an interesting window into how an ice shelf responds to the filling and draining of a meltwater lake within the central basin of a pre-existing doline. Two significant questions arise from our observational timeseries of lake-area, and vertical and horizontal ice-shelf motion:

1. 1. To what extent are our observations of vertical ice-shelf motion in the doline centre and on the doline rim caused by the movement of meltwater into the doline basin from the doline ramp, and the subsequent loading of that meltwater in the doline centre?

2. 2. Are the viscous stresses associated with meltwater loading in the doline sufficient to cause fracture initiation and/or opening?

We explore these questions using a simple, idealised numerical model of viscous ice-shelf flexure described below and in the Supplemental Material. We emphasize that the model is not used to predict or assimilate data, but rather as an exploratory tool to investigate how various observed aspects of the doline's geometry, its central meltwater lake and the state of the surrounding ice shelf have a bearing on the two questions we pose above.

6.1 Model domain

For the purpose of exploring ice-shelf flexure in response to meltwater loading, we adopt a simple geometrical representation of the doline (Fig. 8). The idealized doline is taken to be an axisymmetric circular feature (Fig. 8a) where ice deformation is a function of the radial and vertical coordinates, r and z, respectively. For initial conditions, all ice columns in the domain are assumed to be in local hydrostatic equilibrium, with surface and base elevations of (1 − ρ _i/ρ _sw)H(r) and ( − ρ _i/ρ _sw)H(r), respectively, where ρ _i is the density of ice, ρ _sw is the density of salt water (Table 1), and H(r) is the ice thickness varying as a function of r (Figs 8b, c). Further detail is provided in the Supplementary Material (’Additional Model Details’).

Figure 8. Idealised, axisymmetric model domain used to simulate the doline. (a) Plan view; (b) and (c) cross-section views (full domain and close-up of basin and rim, respectively).

Our assumption of hydrostatic equilibrium as an initial condition means we ignore any processes that could have created disequilibrium conditions. Such processes might have included the initial doline formation and/or subsequent episodes of meltwater loading and unloading in its central basin. However, the flexural readjustment to such processes would likely have been completed within a few months of such perturbations. Another process leading to hydrostatic disequilibrium might be the doline's evolving topography in response to its advection with ice flow. However, as we show below, the stresses involved are small, and perturbations to hydrostatic equilibrium will be minimal as the ice shelf can readjust within days to weeks. Thus, our assumption of hydrostatic equilibrium at the beginning of the melt season (when our observations begin) is reasonable, and, in any case, we do not have the necessary observations to constrain any other ice-shelf thickness distribution.

Given our assumption of hydrostatic equilibrium, and as the mean observed surface elevations of GPS01 and GPS02 averaged over the 2019/2020 melt season (Figs 6b, a) were ~6 m and ~31 m, respectively, which agree with the surface elevations according to the Reference Elevation Model of Antarctica (REMA) tile dated 15 December 2020 (Howat and others, 2019), the corresponding ice thicknesses H of the doline basin and rim are calculated as ~40 and ~240 m, respectively. In our model domain, an annulus surrounding the doline basin (130 < r < 330 m) is taken to represent the doline ramp, and ice thickness linearly increases from 40 m at r = 130 m to 240 m at r = 330 m, which is defined to be the position of the doline rim (Figs 2b, 8b). Beyond the doline rim, ice thickness in the model domain is taken to be uniform at 240 m. We ignore any changes in ice-shelf thickness due to snow accumulation or basal melting through the melt season. We also do not include the surface topographic or ice thickness effects of the meltwater-filled moat surrounding the doline (at about 1 km distance from the doline centre, Fig. 1d) in our model simulations, as sensitivity tests (not shown) showed that this had no substantial effect on model results.

We identify the centre of the doline at r = 0 to represent where GPS01 is located, and take the doline rim at r = 330  m to represent where GPS02 is located. The horizontal dimension of the doline is chosen to be consistent with the horizontal separation of GPS01 and GPS02 and the approximate size of the central basin determined from satellite imagery. Beyond the doline rim, the numerical domain includes an annular area that extends out to r = 2  km. This 2 km scale is chosen from the appearance of meltwater lake features in satellite imagery, which, beyond this 2 km distance, no longer appear to be influenced (in terms of their surface topography and shape) by the presence of the doline (e.g. Fig. 1d). This scale is also consistent with the characteristic length scale for viscoelastic flexure of an ice shelf of a thickness similar to that surrounding South Doline (e.g. MacAyeal and others, 2021).

6.2 Model dynamics and rheology

The model is constructed using COMSOL, a commercial finite-element package, that solves for the ice deformation in the idealised axisymmetric model domain using the Stokes Equations for incompressible laminar flow (e.g. Durand and others, 2009; Christmann and others, 2019; Mosbeux and others, 2020) and an effective viscosity, ν _eff, based on Glen's flow law (e.g. MacAyeal, 1989):

(1)$$\nu_{\rm{eff}} = {B( z) \over 2} \dot{e}_{II}^{( {1}/{n}) -1}$$

where $\dot {e}_{II}$, n and B(z) are, respectively, the second invariant of the strain-rate tensor, the flow-law exponent (which we take to be 3) and a flow-law rate constant appropriate for the ice temperature-depth profile,

(2)$$B\;( z) = B_o \exp{\left(\displaystyle {Q\over 3RT( z) } \right)}$$

where B _o, Q, R and T(z) are, respectively, a flow-law rate constant, the activation energy for ice deformation, the universal gas constant and the temperature (in degrees Kelvin) as a function of z (Table 1). For the temperature-depth profile, we take a linear function that varies from −8 °C at the surface (corresponding to the annual average surface ice temperature observed with a 10 m thermistor string deployed between November 2019 and November 2021 in the ice at our AWS; Fig. 1b) to −2 °C at the base (the freezing point of ocean water at pressure).

Our model only accounts for viscous deformation; it does not account for the elastic component of ice shelf deformation, e.g. as is accounted for in the purely elastic ice-shelf model of Warner and others (2021). The Maxwell timescale (the effective viscosity scale divided by Young's modulus), which estimates the timescale over which elastic effects are significant, is short (~1 d) compared to the time period of the melt season (~100 d). This means that, in our application, viscous deformation will dominate over elastic deformation. We note that a treatment of viscoelastic deformation of similar idealised geometries is possible using a thin-plate approximation (MacAyeal and others, 2021), however the comparability of the scales of the doline basin radius and the ice thickness suggests that a thin-plate approximation is not appropriate; so a full-Stokes stress balance equation is applied in this study.

6.3 Model boundary conditions

Boundary conditions applied to the stress balance equation consist of specification of pressures on the upper and lower surfaces of the idealised model geometry and velocities on the outer extent of the model domain. At the centre of the model domain, the conditions enforcing axisymetry are specified (radial velocity is zero and vertical velocity has no r gradient at r = 0). The upper surface pressure is assumed to be zero for both the doline ramp and the surrounding ice shelf. For the doline basin, the pressure was specified to be ρ _w g d where d is water depth. At the base of the model geometry, the hydrostatic pressure of sea water is specified. The depth of the base of the ice geometry is allowed to vary with time so as to accommodate flexure. At the outer vertical boundary of the model domain (r = 2 km) the radial ice velocity u _r is specified either to be zero or a value of 12.5 m a⁻¹ (see description of model experiments below). The non-zero radial ice velocity value is chosen empirically to reduce convergence of ice flow into the doline centre, and is consistent with a large-scale radial strain rate of 2.0 × 10⁻¹⁰ s⁻¹. The vertical ice velocity at r = 2 km is allowed to vary according to a free slip condition (zero r-derivative). At the centre of the doline, the radial velocity is zero and the vertical velocity is allowed to vary according to free slip (zero r-derivative). These boundary conditions, and other details of the numerical implementation, are described in greater length in the Supplementary Material (Additional Model Details).

Following from our two motivating questions presented above, we conducted three model experiments in which forcing conditions were varied systematically to see the effects of ice-shelf strain rates and stresses (Table 2). The first two experiments are used to determine how the velocity boundary condition at the outer vertical edge of the model domain (r = 2 km) should be specified. In experiment 1, the radial velocity u _r at r = 2 km was specified to be zero. In experiment 2, a radial velocity (at r = 2 km) was specified as 12.5 m a⁻¹, which acted as an imposed strain rate that eliminated the strong ice-flow-convergence into the doline basin that was observed in experiment 1. Experiment 2 was therefore taken as a ‘control experiment’, with which to compare experiment 3. Experiment 3 used the same radial velocity boundary condition at r = 2 km as that used in experiment 2, but was also forced with ablation, meltwater movement, and hence lake filling.

Table 2. Numerical experiments

For experiment 3, the effects of melting, meltwater movement and lake filling were treated as follows. The lake filling the basin of the model domain was assumed to be filled with meltwater coming from just the annular ramp region of the model domain. Melting elsewhere in the model domain was assumed to be zero. The relationship between the ablation rate on the ramp $\dot A( t)$ and the meltwater-lake depth d(t) was taken to be:

(3)$$\displaystyle {\rm d}( t) = \displaystyle {\rho_{\rm i}\over \rho_{\rm w}}\left({( r_{\rm {rim}}^2 - r_{{\rm basin}}^2) \over r_{{\rm basin}}^2} \right)\int_0^t \dot{A}( t') {\rm d}t'$$

where t is time. For simplicity, we took $\dot A( t) = \dot A$ to be constant in time. We note that in reality, the lake's water depth will also increase due to enhanced lake-bottom ablation, which occurs because water has a lower albedo than bare ice or snow (Tedesco and others, 2012). However, as lake-bottom ablation results in no mass change within the lake, we do not account for this process in our model. Another simplification is that we do not account for any change in basin floor geometry, i.e. its downwarping due to being loaded with meltwater. To treat the loss of ice mass caused by ablation and runoff on the surface of the ramp, a negative (upward) load was applied to the upper surface of the ramp equal to $\rho _{\rm i} g \dot A \cdot t$. The lateral flow of meltwater from the ramp to the lake was assumed to be instantaneous.

6.4 Model results

Results of the three 100 d experiments (experiments 1–3 in Table 2) are shown in Figure 9a, where the modelled vertical displacements (red lines) are compared with the observed vertical elevation displacement (black lines). For experiment 1, where both regional strain rate ($\dot {e}_{{\rm rr}}$) and ablation rate and lake depth ($\dot {A}$ and d(t)) are zero, the change in vertical elevation difference and horizontal distance are both of opposite sign compared to the observed values. This indicates that, in the absence of either an externally imposed radial velocity (i.e. producing an externally imposed strain rate) or a meltwater-loading process, the doline geometry tends to decay with time by convergent ice-shelf flow into the doline basin, thereby increasing the ice thickness of the doline basin and raising its surface to be more in accord with the elevation of the ice shelf surrounding the doline. The ice shelf, in other words, deforms to shrink the doline basin radius and reduce the elevation difference between the doline basin and rim.

Figure 9. (a) Comparison of the observed vertical elevation difference (basin minus rim elevations, initiating at 0 cm) over a 100 d period against model experiments 1–3 (red lines labelled with numbers). Black line shows the 5 d running mean of vertical elevation difference (also plotted in Fig. 6d). (b) Maximum deviatoric tensile stress in the radial direction for the top 10 m of ice in the doline basin as a function of time for experiment 3 (red line; right y-axis), and observed horizontal separation change (black line; left y-axis, horizontal distance of separation as a function of time minus the initial horizontal distance of separation) as plotted in Figure 7c. The correspondence of the inferred fracture event on 25 January with the tensile stress curve suggests that a value of about 75 kPa is sufficient to have induced fracture.

For experiment 2, the addition of an externally imposed radial velocity reduces the convergent flow of the ice into the doline basin. This is because this velocity boundary condition creates a radial, azimuthal and vertical strain-rate field in the ice shelf that cancels out the convergence of ice flow into the doline centre that is seen in experiment 1. However, ablation is not specified in experiment 2, and therefore there is only minimal change in both the ice thickness in the doline basin and in the surface elevation of the doline basin relative to the rim. As shown in Figure 9a, the vertical elevation change for experiment 2 is roughly zero over 100 d. However, for the first 15 d, the elevation change of experiment 2 agrees with the slightly negative slope of the observed elevation change. This suggests that the divergent strain rate externally imposed on the doline by the surrounding ice flow is controlling the doline geometry prior to the onset of substantial melting and meltwater-loading-induced flexure.

For experiment 3, the combination of the externally imposed radial velocity (and hence strain rate) of experiment 2, alongside the addition of ablation and lake filling (this experiment) produces the best agreement between modelled and observed vertical elevation change (Fig. 9a).

The results of experiment 3 were examined further to assess whether stresses induced by the combination of the externally-imposed radial velocity and meltwater-loading are sufficient to induce fracture of the ice (either to partial or full thickness) such as implied by our observation of increased horizontal displacement between our two GNSS for the 4 d after 25 January 2020 (Fig. 7c). Figure 9b compares the horizontal separation change between the rim and basin GNSS antennae and the maximum deviatoric tensile stress in the radial direction for the top 10 m of ice in the basin, as simulated by experiment 3. The initial stresses at day 0 of experiment 3 are slightly below 60 kPa. They rise to over 100 kPa at the end of the 100 d simulation. At the time at which the fracture occurs (about day 40 after the initial date of the observations and the model initial condition), the maximum tensile stress is about 75 kPa. We note, however, that our simplistic viscous modelling approach is not explicitly able to simulate brittle fracture initiation.

We can use our model results (described above) to answer the two questions that motivated the modelling effort. The first question, whether ablation and lake filling are sufficient to explain the elevation change of the doline basin relative to the rim, is answered affirmatively by comparison of experiments 2 and 3. Experiment 3, which produces a downward movement of the basin relative to the rim roughly consistent with observations, simulates the deformation associated with of meltwater loading in the central lake. Experiment 2, taken as a control experiment where the effects of ablation, meltwater movement and lake filling were taken to be absent, did not produce downward movement of the basin relative to the rim in a manner consistent with observations.

The second question, whether tensile stresses caused by the combination of the imposed radial velocity at r = 2 km and meltwater-induced flexure are sufficiently high to initiate fracture, is answered affirmatively by the comparison shown in Figure 9b. Maximum tensile stress in the top 10 m of the basin (i.e. where we assume fractures initiate) reaches a threshold value of ~75 kPa after ~40 d. Taking the initial condition of the model run to be equivalent to the state of the doline at the beginning of the observation period (15 December) implies that the apparent fracture event on 25 January corresponds to the tensile stress reaching a threshold value of ~75 kPa (Fig. 9b), which is a reasonable stress value to initiate ice-shelf fracture (e.g. Albrecht and Levermann, 2012; Banwell and others, 2013). We note that the background tensile stress in the ice shelf near its upper surface is slightly negative (indicating compressive stress) prior to the time when the lake starts to fill. This is because the tendency for ice-flow convergence into the centre of the model domain is not entirely eliminated by the specification of a radial velocity at the outer boundary (r = 2 km) in experiments 2 and 3.

7. Conclusions and future directions

The field data presented in this study consist of records from two GNSS stations and one timelapse camera from South Doline on the north George VI Ice Shelf, Antarctic Peninsula (Fig. 1), during the record high 2019/2020 melt season (Banwell and others, 2021). The doline was chosen as a study site due to its natural constraints on the central lake's geometry, its catchment area, and hence its water depth as a function of the mean surface ablation rate in the doline.

Our timelapse photos documented a meltwater lake forming in the centre of the doline, beginning in early December 2019 and reaching its maximum area around 9 January 2020 (Fig. 5b and Video S1). The surface water level (and hence lake area) appears to have been limited by the formation of a moulin close to the bottom of the doline ramp, which we observed in our timelapse imagery (Fig. 4 and Video S1). We suggest that this moulin acted like a safety drain in a bathtub, limiting the lake's surface water-level, and thus potentially the loading potential of the surface lake on the surrounding ice shelf.

Between mid-December 2019 and late March 2020, our GNSS measurements captured changes in vertical motion (Fig. 6d) and horizontal displacement between the doline rim and its centre (Fig. 7c). The vertical-displacement timeseries showed that the doline centre dropped relative to the rim, which we infer was a result of surface ablation (i.e. unloading) on the doline ramp and meltwater loading in the doline basin. Notably, this downwards motion of the doline centre was observed to continue until late March, which is long after the observed peak in lake area (9 January 2020). This observation suggests continued viscous-flexure of the doline basin in response to ponded meltwater that remained in it. To explain this, we suggested four possibilities: (1) preferential draining via the moulin (mentioned above) of a relatively shallow region of the lake, which was not fully connected to other deeper regions of the lake, meaning that the total volume of meltwater in the lake may have remained constant, or even increased, after 9 January, despite the observed decrease in the total lake area; (2) draining of the surface lake water into an expanding buried lake or firn aquifer below; (3) downwarping of the thin ice-shelf under the doline basin floor due to the meltwater loading in the doline; and (4) enhanced localised lake-bottom ablation.

The horizontal GNSS displacement timeseries between the basin rim and centre revealed a sudden and short-lived (4 d) acceleration and change in trajectory, beginning on 25 January 2020 followed by a gradual return to the background ice velocity. This event was interpreted as a fracture initiation and/or opening in response to meltwater-loading in the central lake basin. We also observed fractures such as this on the doline ramp in our timelapse imagery (Fig. 4c and Video S1). We also noted that this fracture opened at low tide, implying that tidal currents beneath the ice shelf may help to drive fracture initiation and/or opening, in addition to the meltwater loading in the doline basin. Further investigation of this phenomenon in a future study is warranted.

Our idealised modelling addressed two specific questions arising from the field observations that we explored with three complementary model experiments (Fig. 9). Results are as follows. First, the comparison between experiments 1 and 2 showed that a large-scale extending (positive) strain rate between the doline basin and rim is needed to counteract convergent ice flow into the basin that would otherwise cause the doline to fill. Second, the comparison between experiments 2 and 3 showed that the long-term relative vertical motion observed in the GNSS data is consistent with ablation (unloading) on the doline ramp area and meltwater movement into the doline basin (loading). Third, the stress distribution in experiment 3 showed that the stresses associated with flexure due to meltwater loading are of a magnitude that is sufficient (~75 kPa (Albrecht and Levermann, 2012)) to initiate fracture on the doline ramp, which we also inferred from our horizontal displacement timeseries (Fig. 7), and observed in our timelapse imagery (Fig. 4c and Video S1).

Comparisons between our field observations and model experiments suggest several avenues for future research. Regarding field observations, either at South Doline or another doline, it would be useful to have measurements of ice-shelf thickness and vertical strain rate changes (e.g. from airborne or ground-based radar, or ApRES) across the doline during the summer melt season, as well as a higher spatial-resolution dataset of vertical and horizontal surface motion (from a denser array of GNSS stations) across the ice shelf. Seismometer arrays would be helpful to detect better the timing of fracture initiation and/or opening around the doline ramp, as well as to provide evidence for the possibility of hydrofracturing to the ocean beneath. Finally, pressure transducers would be useful to measure lake depths at various locations in a doline, in order to estimate lake volumes (as well as areas, as in this study). An increase in observational field data would facilitate more sophisticated modelling to be undertaken, e.g. using a viscoelastic approach, and including the representation of non-circular lakes and non-uniform ice thicknesses.

While future research is needed, our observations represent the first in-situ documentation of the initiation and/or opening of a ‘ring fracture’ in response to the filling and loading of a meltwater lake on an ice shelf. This phenomenon is at the root of the ‘chain reaction’ mechanism that was proposed to have helped trigger the drainage of thousands of lakes on the Larsen B Ice Shelf, and hence may have contributed to that ice shelf's rapid break up in 2002 (Banwell and others, 2013).