An updated seabed bathymetry beneath Larsen C Ice Shelf, Antarctic Peninsula

In recent decades, rapid ice shelf disintegration along the Antarctic Peninsula has had a global impact through enhancing outlet glacier flow and hence sea level rise and the freshening of Antarctic Bottom Water. Ice shelf thinning due to basal melting results from the circulation of relatively warm water in the underlying ocean cavity. However, the effect of sub-shelf circulation on future ice shelf stability cannot be predicted accurately with computer simulations if the geometry of the ice shelf cavity is unknown. To address this deficit for Larsen C Ice Shelf, West Antarctica, we integrate new water column thickness measurements from recent seismic campaigns with existing observations. We present these new data here along with an updated bathymetry grid of the ocean cavity. Key findings include a relatively deep seabed to the southeast of the Kenyon Peninsula, along the grounding line and around the key ice shelf pinning-point of Bawden Ice Rise. In addition, we can confirm that the cavity's southern trough stretches from Mobiloil Inlet to the open ocean. These areas of deep seabed will influence ocean circulation and tidal mixing and will therefore affect the basal-melt distribution. These results will help constrain models of ice shelf cavity circulation with the aim of improving our understanding of sub-shelf processes and their potential influence on ice shelf stability.

The datasets are comprised of all the new point measurements of seabed depth. We present the new depth measurements here, as well as a compilation of previously published measurements. To demonstrate the improvements to the sub-shelf bathymetry map that these new data provide we include a gridded data product in the Supplement of this paper, derived using the additional measurements of both offshore seabed depth and the thickness of grounded ice. The underlying seismic datasets that were used to determine bed depth and ice thickness are available at https://doi.org/10.5285/315740B1-A7B9-4CF0-9521-86F046E33E9A (Brisbourne et al., 2019), https://doi.org/10.5285/5D63777D-B375-4791-918F-9A5527093298 (Booth, 2019), https://doi.org/10.5285/FFF8AFEE-4978-495E-9210-120872983A8D (Kulessa and Bevan, 2019) and https://doi.org/10.5285/147BAF64-B9AF-4A97-8091-26AEC0D3C0BB (Booth et al., 2019).

1 Introduction

The loss of Antarctic ice shelves is of global significance for two reasons. First, ice shelves provide a buttressing force – controlled by the geometry and stress regime of the ice shelf – to the glaciers or ice streams that feed them. Although loss of the floating ice shelf makes only a small direct contribution to sea level rise, the removal of buttressing results in acceleration of the tributary glaciers, enhancing their current contribution to sea level rise (Rignot et al., 2004; Scambos et al., 2004; Rott et al., 2002; Fürst et al., 2016). Secondly, basal melting of ice shelves produces cold and low-salinity water that influences Antarctic Bottom Water (AABW) formation, which in turn affects the properties of the global oceans (Jacobs, 2004).

Over recent decades, there has been a southwards progression of ice shelf loss along the eastern Antarctic Peninsula. The disintegration of the Larsen A Ice Shelf (in 1995) and the Larsen B Ice Shelf (in 2002) resulted in a step increase in flow of the grounded glaciers that formerly fed these ice shelves (e.g. Khazendar et al., 2015). This increase in glacier flow resulted in accelerated sea level rise and increased freshening of dense AABW (Jullion et al., 2013). In a number of cases, ice shelf retreat has been attributed to atmospheric warming (Vaughan and Doake, 1996; Rott et al., 1998; Skvarca et al., 1999). With the Antarctic Peninsula exhibiting one of Earth's highest rates of atmospheric warming during the late 20th century (Vaughan et al., 2003), the long-term viability of the Larsen C Ice Shelf (LCIS) is in question. However, Holland et al. (2015) demonstrated that the thinning of LCIS over the last decade is a result of both atmospheric and oceanic influence in almost equal measure. For the remaining ice shelves on the Antarctic Peninsula, the relative contribution to their future stability by basal melt from incursions of relatively warm ocean water and increased surface melting by a warmer atmosphere is still unknown.

To improve projections of the effects of basal melt on ice shelves, knowledge of the geometry of the ocean cavity beneath is vital (Mueller et al., 2012; Jenkins et al., 2010; Grosfeld et al., 1997; Goldberg et al., 2019; Pattyn et al., 2017). Models of sub-shelf circulation are critically dependent on cavity geometry, particularly in regions where the influence of strong tides is topographically constrained (e.g. Mueller et al., 2012). Ongoing efforts to model ocean processes beneath LCIS suffer from inadequate knowledge of cavity geometry because seabed depth is poorly sampled (Brisbourne et al., 2014). Improving knowledge of cavity geometry is crucial for LCIS because the sparse existing data suggest the presence of large-scale seabed features capable of guiding ocean currents and inducing significant tidal mixing. It is impossible for computer simulations to accurately predict the future influence of the ocean on LCIS without knowledge of the geometry of such features.

Although they are labour intensive, seismic methods remain the most reliable method for determining sub-shelf cavity geometry. Airborne and ground-based radar are used extensively to map ice thickness but cannot penetrate the sub-shelf cavity. Autonomous underwater vehicles (AUV) provide another direct measurement of sub-shelf bathymetry but with limited coverage at present (e.g. Jenkins et al., 2010). Inversion for water column thickness using airborne gravity measurements is sensitive to assumptions about local density variations, such as sediment infill, and may lead to inaccurate results (Brisbourne et al., 2014). Recent studies using gravity inversion combine data from multiple methods to address these assumptions (e.g. Muto et al., 2016).

2 Location and previous work

LCIS, the largest ice shelf on the Antarctic Peninsula at around 44 000 km² (Cook and Vaughan, 2010), lies just south of the recently collapsed Larsen A and B ice shelves (Fig. 1). The geometry of LCIS's sub-shelf cavity has previously been measured in detail at specific locations only (Brisbourne et al., 2014): this campaign was designed to target locations where an existing inversion of gravity measurements indicated areas of significant control over sub-shelf circulation (Cochran and Bell, 2012). However, uncertainties associated with gravity inversions for bathymetry result in large areas of unknown geometry, specifically beneath LCIS (i) away from the western grounding line, (ii) away from the ice front and (iii) in the south.

Figure 1 Map of seismic points used in the gridded bathymetry product of this study. The approximate path of the ice shelf rift that resulted in the calving of iceberg A68 is highlighted (Jansen et al., 2015). The background is MODIS imagery (Scambos et al., 2007), predating the break-off of iceberg A68 along the rift.

We build on a number of published sources of bathymetric data with new observations from four recent field campaigns. The existing bathymetric data used in the gridding process here (Fig. 1, blue dots) are derived from a targeted seismic bathymetry survey, seismic refraction experiments and drill site measurements (Brisbourne et al., 2014; Nicholls et al., 2012). The depth to grounded ice and known offshore bathymetry of Bedmap2 is included in the gridding process (Fretwell et al., 2013). Surface elevation and ice thickness measurements at Bawden Ice Rise (BIR) are also included (Holland et al., 2015). Here, we integrate these existing data with the new measurements of seabed depth. All data are then gridded to obtain a new bathymetry map of LCIS. We recognize that users of these data are likely to create bathymetric grids using their preferred gridding method, specific to their preferred model and resolution. As such, the gridded product presented here is presented as an aid to discussion and to highlight the value of these new data.

3 Data acquisition and processing of new observations

3.1 Data acquisition

In December 2016, 14 seismic bathymetry measurements were made across LCIS, targeting areas of sparse data coverage (Fig. 1, magenta dots). The seismic source consisted of a sledgehammer with a plate stamped into the snow surface, or dug down to a shallow ice layer, to improve source consistency of the shots at that location for stacking purposes. The 24 Georod receivers (Voigt et al., 2013) were buried to 0.3 m depth, at 10 m spacing, with a 30 m offset between the shot and the first receiver. Burying sensors in this way ensures good coupling and provides protection from wind-induced noise. Georods consist of four geophone elements in a series, which improves the signal to noise ratio. We recorded 2 s records at a 0.125 ms sample interval with a 24-channel data logger. At each site, ∼20 hammer blows were recorded using a geophone adjacent to the hammer plate to initiate recording. An additional stack of 10 hammer blows was also recorded for on-site evaluation of the seismic reflection strength. To determine an accurate surface elevation a dual-frequency GPS system was deployed for the duration of the seismic acquisition at each site.

These data are supplemented by bathymetry measurements from an additional 16 seismic refraction and reflection surveys across LCIS (Fig. 1, orange, black and yellow dots). Although many of these experiments targeted depth profiles of the firn, the data are suitable for ice shelf thickness and seabed depth measurement. The acquisition procedure is similar to that described above, and therefore data quality and uncertainties are similar. Details of the acquisition parameters for each experiment are presented in Table 1 with further details in the metadata of each data archive.

Table 1 Field acquisition parameters for new data presented here.

^* Note that due to the use of a range of acquisition geometries, the specific geometry used at every refraction experiment site is included in the data repository. BAS refraction data were published previously in Brisbourne et al. (2014) but are again presented here as they form part of the analysis of new data.

Figure 2 presents an example of a seismic gather formed of 10 hammer blows stacked during acquisition. Clear ice base and seabed arrivals, as well as multiples thereof, are observed. Where necessary, to help identify reflections, a frequency wavenumber filter was used to suppress ground roll that may mask the ice base reflection. An automatic gain control filter and semblance analysis was also used as required to identify arrivals. However, ice base and seabed reflection travel times were measured on raw seismic records, even if a filter was required to help identify arrivals. A relatively thin ice shelf will result in the ice base reflection arriving within ground roll noise (see Fig. 2). In these cases, surface multiples of the ice base reflection were used to calculate the primary two-way travel time through the ice column.

Figure 2 Example hammer and plate seismic shot gather with readily identified primary seismic reflections and multiples. The primary ice base reflection at 190 ms is masked by ground roll signal at the far offsets.

3.2 Seismic velocities in ice and water and thickness measurement

We follow the procedures outlined in Brisbourne et al. (2014) to convert from travel time to thickness. Values of seismic velocity are required to convert travel times to layer thickness or depth. A mean seismic velocity in the water column of 1445±1 m s⁻¹ was derived during conductivity–temperature–depth (CTD) measurements made beneath northern and southern LCIS by Nicholls et al. (2012).

The seismic velocity profile in the upper 100 m of the ice shelf, which includes the firn, was measured using the shallow refraction experiments presented here (see Table 1). At each of the refraction sites, a series of surface shots was recorded with increasing receiver spacing. The first arrivals were picked and converted to a velocity–depth profile using the method described by Kirchner and Bentley (1990). This method relies on a monotonic increase in velocity with depth, an assumption that is supported by observations of smoothly varying travel times. Maximum velocities at 100 m depth calculated by inversion of the refraction measurements range from 3698 to 3916 m s⁻¹. Below 100 m depth, we assume that ice density is constant and seismic velocity depends on ice temperature alone. The CTD measurements of Nicholls et al. (2012) indicate an ice base temperature of −2 ^∘C. Temperature measurements within the ice column indicate an approximately linear temperature profile with a small range (−14 ^∘C at 100 m depth to −2 ^∘C at the ice base; Nicholls, unpublished data) and therefore, using the relationship of Kohnen (1974), a small range of seismic velocities (3800–3827 m s⁻¹). Therefore, below 100 m we linearly interpolate between the velocity measured by seismic refraction at 100 m depth and an ice base velocity calculated from the temperature–velocity relationship of Kohnen (1974). Where a bathymetry measurement and seismic refraction experiment are not coincident, results from the closest seismic refraction experiment are used to determine ice thickness.

Measurement of the surface elevation allows for the estimation of ice thickness assuming freely floating ice. These estimates can guide the identification of ice base reflections in the data. The EIGEN-GL04C geoid level (Forste et al., 2008) is removed from the elevation, and an empirical relationship determined by Brisbourne et al. (2014) is used to calculate ice thickness: geoid-corrected height, $h=(\mathrm{0.113}\pm \mathrm{0.005})H+(\mathrm{5.003}\pm \mathrm{1.525})$, where H is ice column thickness in metres. This relationship accounts for firn thickness, which affects mean density. The absence of a clear ice base reflection is not necessarily a result of low-quality data. Under certain conditions, particularly in ice shelf suture zones, poorly consolidated marine ice at the base of the ice shelf may result in a weak or absent seismic reflection. At site PRHB04 (see Table 2), in the absence of a clear ice base reflection, we calculate the ice thickness and its uncertainty from the surface elevation using the empirical relationship described above.

Table 2 Location and seabed depth measurements and associated uncertainty of all new points used in this study.

3.3 Uncertainties

Errors in picking reflections, seismic velocities and seabed topography all contribute to uncertainties in ice and water column thickness calculations. Picking of ice and seabed reflections was repeated three times for each site in order to quantify the error, indicating a maximum picking error in the seismic reflections of 0.5 ms. Ice column thickness is determined by the difference in seabed and ice base arrival times and therefore has an uncertainty of 1.0 ms. We assume a conservative estimate of the uncertainties in seismic velocity in the ice column of 30 m s⁻¹ (Kirchner and Bentley, 1990; Rosier et al., 2018). The presence of marine ice in suture zones (Kulessa et al., 2019) or significant warm refrozen ice within the firn column (Hubbard et al., 2016; Ashmore et al., 2017) may result in seismic velocities that deviate from the standard model and introduce greater uncertainty in measured velocities. However, a previous study highlighted the consistency between seismically derived ice thickness measurements and those from surface elevation measurements (Brisbourne et al., 2014). Importantly, where an ice base reflection can be identified the calculated thickness of the water column is independent of the ice velocity–depth profile used.

Ice base and seabed topography can introduce additional uncertainty to thickness measurements (Nost, 2004). Calculations of ice and water column thickness from travel times assumes that reflectors are planar and horizontal. Such a geometry results in a characteristic curvature, or move-out, of travel times with increasing receiver offset. Assuming an isotropic seismic velocity structure, any deviation from standard move-out is indicative of dip at the reflecting interface. Brisbourne et al. (2014) used observed deviations from standard move-out to demonstrate that topography across LCIS causes a maximum error in seabed depth of <10 m. However, a full assessment of the error introduced by bed topography would require multiple measurements at each site at different angles across the slope and this uncertainty is therefore not included here.

We calculate seabed depth by removing ice and water-cavity thickness from surface elevation data. Uncertainties in dual-frequency GPS elevation measurements are up to ±40 mm. Where direct surface elevation measurements are not available (dual-frequency GPS measurements not made; see Table 1), the REMA surface digital elevation model (DEM) of 2017 at 8 m resolution is used (Howat et al., 2019), resulting in an absolute elevation uncertainty in these areas of ±2 m. No tidal correction is made to surface elevations, resulting in additional uncertainty of ±2 m (Brisbourne et al., 2014).

Based on the above uncertainty sources but excluding the unknown seabed slope, we calculate uncertainties in seabed depths at each site and present these with the seabed depths in Table 1. Uncertainty at site PRHB04, where no ice base reflection was observed, has been calculated using the range of ice thickness values indicated by the empirical surface elevation relationship of Brisbourne et al. (2014) and the resultant uncertainty in water column thickness.

4 Bathymetry gridding and results

We include a 1 km horizontal resolution bathymetry grid of LCIS's cavity in the Supplement. To produce this grid we interpolated between all existing and new seabed depth measurements. We augmented these data with the grounding line position, grounded-ice bed depths and offshore bathymetry derived from Bedmap2 (Fretwell et al., 2013). The measured bed geometry of Bawden Ice Rise was derived from Holland et al. (2015) (see Table A1). We use a natural neighbour interpolation implemented in the grid data function of MATLAB (Release 2019a), which is well suited to a dataset with an uneven distribution of data points (Sibson, 1981). Importantly, the fit to these points does not “overshoot”, which would result in interpolated values that are higher or lower than known values. A weakness of this method is that where seabed topography changes rapidly with respect to data coverage the seabed may not be well constrained. This results in a discrepancy where the interpolated seabed depth is shallower than the ice draft as reported in Bedmap2. Therefore, in the gridded product we deepen the seabed where this discrepancy occurs and assign a seabed depth to be equivalent to the Bedmap2 ice draft plus a minimum water column thickness of 10 m. This ensures that all interpolated seabed depths are consistent with the Bedmap2 ice thickness.

Figure 3 presents a map of seabed elevation in the LCIS region, resulting from gridding of all available data, as described above, with a minimum cavity thickness of 10 m when compared to Bedmap2 ice draft. The analysis and application of such grids is of course dependent on the limitations of data coverage and the gridding method used. Our method does not superimpose any additional constraints on the resultant bathymetry, as highlighted by the bullseye deepening around the single data point to the south of the Kenyon Peninsula that in reality may form a linear trough. However, the sparse data coverage in that area precludes any definitive knowledge of bed geometry and any gridded product will require careful interpretation.

Figure 3 Updated bathymetry map of Larsen C Ice Shelf with large-scale features highlighted. For clarity, elevations above 0 m are unscaled, and we label only the −400 and −600 m contours. The brown line represents the Bedmap2 grounding line (GL) (Fretwell et al., 2013). The orange line represents the ice front on 20 December 2018, highlighting the new ice front following the calving of iceberg A68 determined from MODIS imagery (Vermote and Wolfe, 2015). MOI – Mobiloil Inlet. The black arrows highlight likely oceanographic conduits, as discussed in the text.

Figure 4 presents the difference between the grid corrected for Bedmap2 ice draft and the original interpolated grid, highlighting the areas of the cavity where the data coverage and interpolation method are least reliable. In general, this discrepancy occurs where topography is changing rapidly with respect to data coverage, such as along the grounding line. To avoid issues such as this, other methods of interpolation may be preferred. For example, knowledge of past ice flow may be used to prescribe channels around isolated data points, or onshore slopes may be continued into the cavity. Similarly, the prescribed minimum cavity thickness of 10 m may be adapted to the resolution and type of oceanographic model used. This gridded product is therefore not viewed as a definitive bathymetry for use by the oceanographic community but is used here to highlight the value of these new data. No matter what method is applied, there are intrinsic weaknesses of gridding with interpolation in regions of sparse data coverage, and, by their nature, uncertainties cannot be quantified readily.

Figure 4 Difference between the grid corrected for the Bedmap2 ice draft mismatch and the original natural neighbour interpolation grid in metres, highlighting areas where the cavity has been forced to 10 m.

5 Data availability

The underlying seismic datasets that were used to determine bed depth and ice thickness are available at https://doi.org/10.5285/315740B1-A7B9-4CF0-9521-86F046E33E9A (Brisbourne et al., 2019), https://doi.org/10.5285/5D63777D-B375-4791-918F-9A5527093298 (Booth, 2019), https://doi.org/10.5285/FFF8AFEE-4978-495E-9210-120872983A8D (Kulessa and Bevan, 2019) and https://doi.org/10.5285/147BAF64-B9AF-4A97-8091-26AEC0D3C0BB (Booth et al., 2019).

6 Discussion and significance of the dataset

A number of key features that will influence tidal and oceanic circulation through the sub-shelf cavity and thus affect basal melt rates and meltwater circulation are apparent: (1) a relatively deep seabed surrounds Bawden Ice Rise, a key pinning point of LCIS, where LCIS is closest to floatation (Holland et al., 2015). The bathymetry of this area therefore plays a key role in the ice shelf's future stability. A deep seabed here may alter the strong tidal currents that are thought to induce melt in this region (Mueller et al., 2012) and also reduce the likelihood of regrounding following any slight thickening. (2) The southern trough, to the north of the Kenyon Peninsula, extends from Mobiloil Inlet to the ice front. Nicholls et al. (2012) highlighted this deepening in southern LCIS as a potential conduit for High Salinity Shelf Water (HSSW) that may access the deeper ice at the grounding line, providing vigorous melting. Similarly, the updated bathymetry also confirms that the Jason Trough in the north also continues through to the open ocean, to the north of Bawden Ice Rise. (3) The sub-shelf cavity to the southeast of the Kenyon Peninsula is relatively deep. Again, Nicholls et al. (2012) highlight this location as potentially important to the supply of HSSW that sustains melt at the grounding line. (4) All additional point measurements confirm that the sampled sub-shelf cavity is particularly deep close to the grounding line between Mobiloil Inlet and the Cole Peninsula. Sub-shelf circulation models highlight that the grounding line, where shelf ice is thickest and therefore deepest provides a key site for basal melt (Mueller et al., 2012). However, the inconsistency between the interpolated grid and Bedmap2 bathymetry at the grounding line (Fig. 4) highlights the remaining shortcomings of the bathymetric data close to the grounding line.

These data provide a valuable product for the study of ice–ocean interaction beneath LCIS. The interaction of these newly determined cavity features with the sub-shelf circulation pattern requires detailed oceanic modelling to ascertain their importance. The updated bathymetry is a prerequisite to estimating the contribution of sub-shelf melt to thinning of the ice shelf and the contribution of that melt to the global ocean system. Previous studies have highlighted the importance of accurate bathymetry beneath LCIS but until now have lacked information on the major troughs delineated by these data. The provision of the spot measurements will allow users to re-grid using other algorithms if required and allow for rapid assimilation of any new data points that become available in the future. This is not a definitive dataset and additional data points that address gaps in the current coverage will always be of value to reduce uncertainty where interpolation has been necessary. As the resolution of ocean models improves, the requirement for greater certainty with regards to small-scale features will also increase.

Appendix A

Table A1 Previously published seabed depth measurements included in the gridding process (Holland et al., 2015; Brisbourne et al., 2014).