Comparing numerical ice-sheet model output with radio-echo sounding measurements in the Weddell Sea sector of West Antarctica

Abstract

Numerical ice-sheet models are commonly matched to surface ice velocities from InSAR measurements by modifying basal drag, allowing the flow and form of the ice sheet to be simulated. Geophysical measurements of the bed are rarely used to examine if this modification is realistic, however. Here, we examine radio-echo sounding (RES) data from the Weddell Sea sector of West Antarctica to investigate how the output from a well-established ice-sheet model compares with measurements of the basal environment. We know the Weddell Sea sector contains the Institute, Möller and Foundation ice streams, each with distinct basal characteristics: Institute Ice Stream lies partly over wet unconsolidated sediments, where basal drag is very low; Möller Ice Stream lies on relatively rough bed, where basal drag is likely larger; and Foundation Ice Stream is controlled by a deep subglacial trough with flow-aligned bedrock landforms and smooth unconsolidated sediments. In general, the ice-sheet model represents each ice-stream system well. We also find that ice velocities do not match perfectly in some locations, and that adjustment of the boundaries of low basal drag, to reflect RES evidence, should improve model performance. Our work showcases the usefulness of RES in calibrating ice-sheet model results with observations of the bed.

1. Introduction

The Weddell Sea (WS) sector of the West Antarctic Ice Sheet (WAIS) includes three major ice streams: Institute Ice Stream (IIS), Möller Ice Stream (MIS) and Foundation Ice Stream (FIS). The trunks of IIS and MIS are separated by the Bungenstock Ice Rise (BIR), with the former flowing to its grid south and west and the latter to grid east (Fig. 1a, Supplementary Fig. S1a). Jeofry and others (2018a) recently developed a revised digital elevation model (DEM) of the bed topography of the region, extending 135 km grid south of the BIR, 195 km grid east of FIS, over the Pensacola Mountains and 185 km grid west of IIS. The region is thought to be at a physical threshold of change, as the grounding lines of the IIS and MIS are perched atop a steep reverse bed slope towards the >1.8 km deep Robin Subglacial Basin (Ross and others, 2012), and ocean modelling predicts an influx of warm water that may lead to enhanced melting at the grounding line later this century (Hellmer and others, 2012). Critical to evaluating the sensitivity of the region to ocean-driven melting are reliable numerical ice-sheet models. Using the BEDMAP2 depiction of subglacial topography, a number of models have shown the region to be extremely sensitive to melting at ice stream grounding lines. For example, Wright and others (2014) used the BISICLES ice-sheet model, which includes grounding line physics (Cornford and others, 2013), to determine the response of WS ice streams under ocean warming scenarios. In this model, using today's values of ice accumulation and basal ice-shelf melting as a control experiment, the IIS lost mass within ~200 years and experienced accelerating mass loss in the subsequent 1800 years. By the end of the model run, it was losing ~8 km³ of ice annually. Similarly, MIS lost mass at a rate of 1–2 km³ per year. This behaviour suggests that the two ice streams are not necessarily in equilibrium with the imposed forcing fields. When melting was increased, the two ice streams reacted even faster by deglaciating the Robin Subglacial Basin within a few centuries (Wright and others, 2014). While this modelling exercise suggests that the IIS and MIS are highly sensitive to future ocean-driven change (see also Cornford and others, 2015), the validity of the model in terms of its set-up to mimic the present ice sheet has yet to be evaluated thoroughly.

Fig. 1. Map of the Weddell Sea sector encompassing the Institute Ice Stream, Möller Ice Stream and Foundation Ice Stream. (a) the bed elevation DEM with 500 m resolution (Jeofry and others, 2018a), with subglacial hydrological pathways (blue lines) and locations of subglacial lakes (green polygons) (Jeofry and others, 2018b); (b) the total roughness index of the IMAFI survey (Rippin and others, 2014); (c) InSAR-based ice velocity (Rignot and others, 2017); (d) BISICLES modelled ice surface speed; (e) percentage difference between measured and modelled surface ice velocity (positive values are where measurements exceed the model velocities); and (f) basal friction estimated by the ice-sheet model inversion. RES transects, displayed in Figure 2, are shown by white dashed lines. a, b and c include MODIS MOA satellite imagery (Haran and others, 2014). The black line denotes the boundary of grounded and floating ice (from MODIS MOA data - see Jeofry and others, 2018a), whereas black/white dotted lines denote the margins of Institute, Möller and Foundation Ice Streams. Note that more RES survey lines and transect compilations, referred to in the text, are provided in Supplementary Information.

Fig. 2. RES transects and model putputs from each of the three major ice streams in the Weddell Sea sector. (a) A–A’ (Institute Ice Stream). (b) B–B’ (Möller Ice Stream). (c) C–C’ (Foundation Ice Stream). Locations are provided in Figure 1. Methods used to calculate bed-echo reflectivity in c are given in Jeofry and others (2018b).

There is firm recognition of the importance of grounding-line processes in modulating ice-sheet change (Hellmer and others, 2012; Favier and others, 2014; Joughin and others, 2014; Rignot and others, 2014; Khazendar and others, 2016; Siegert and others, 2016). However, despite recent advances in modelling (Martin and others, 2015; Ritz and others, 2015; Thoma and others, 2015; Whitehouse and others, 2017), ice-sheet models treat the grounding zone imperfectly. Bed topography, basal conditions and subglacial processes are neither known precisely nor modelled realistically at and near this critical interface of the ice-sheet system (Jeofry and others, 2017; Jeofry and others, 2018a). Existing observations of grounding zones (i.e. Horgan and others, 2013) point to a complex transition between floating and grounded ice that needs to be understood fully and characterised well to predict how ice sheets behave in such places.

Here, we assess the BISICLES model output in the WS sector against radio-echo sounding (RES) measurements to understand: (1) whether there is geophysical support for the modelling approach, especially in its use of variations in basal traction to deliver ice velocities that match InSAR measurements of ice flow (Rignot and others, 2017); (2) whether there are inconsistencies between model depictions of bed processes and measurements; and (3) how geophysical measurements can be used to refine the ice-sheet modelling.

2. Methods

The RES data used in this study were compiled from two main sources: flights conducted by the Center for Remote Sensing of Ice Sheets (CReSIS) as a part of the NASA Operation IceBridge (OIB) missions in 2012, 2014 and 2016; and a survey of the IIS and MIS (referred to as IMAFI, which stands for Institute Möller Antarctic Funding Initiative) undertaken in 2010/2011 by the British Antarctic Survey. Details of the RES systems, acquisition and data processing are provided in Jeofry and others (2018a).

We compare the model results with RES data in two ways. First, by comparing the model fields of ice surface velocity and basal friction against published measurements of ice velocity and basal roughness, to observe regional similarities and differences. Second, by making direct comparisons of model results along RES transects, to evaluate site-specific consistencies and inconsistencies. Below we describe how basal roughness from geophysical data was calculated, and provide some details of the BISICLES numerical model and in particular its treatment of basal traction.

2.1 Total roughness

The roughness datasets we use in this paper were first published by Rippin and others (2014), who offer full details of the method. To summarise this, total roughness measurements (Fig. 1b, Supplementary Fig. S1b, Rippin and others, 2014), extracted from the IMAFI data for IIS and MIS, are essentially a depiction of bed obstacle amplitude or vertical roughness (Li and others, 2010). Roughness calculations were undertaken using the Fast Fourier Transform (FFT) approach. Prior to analysis, data gaps were identified to ensure continuous along-track data. Small data gaps (~10 missing data points = ~10 m) were replaced by a linear interpolation, while large data gaps were characterised as ‘broken’, and no FFT analysis was carried out. Additionally, the data need to be equally spaced so that data smoothing can be carried out over a constant distance. Data were therefore resampled at a regular 10 m step size. The mean of the bed elevation, Z _m(x), was subtracted from the bed elevation profile, Z _e(x), to remove topographic anomalies over a window similar in size to the FFT window, which is 2^N samples where N = 5 (32 data points = 320 m):

(1)$$Z_0\lpar x \rpar = Z_{\rm e}\lpar x \rpar -\; Z_{\rm m}(x),$$

where x is the horizontal position. A FFT was later applied to obtain the spectral power density S(k):

(2)$$S\lpar k \rpar = \displaystyle{1 \over l}\vert {{\bar{Z}}_0\lpar k \rpar } \vert ^2,$$

where $\bar{Z}_0\lpar k \rpar$ is the FFT of Z ₀(x) and l is the length of the profile. Total roughness index, ξt was later defined as the integral of the FFT power spectra (Siegert and others, 2004; Bingham and Siegert, 2009; Li and others, 2010) in each moving window:

(3)$$\xi t = \mathop \int \limits_0^{\infty ^{}} S\lpar k \rpar {\rm \; d}k.$$

In Rippin and others (2014), the roughness measurements were interpolated using the ‘Topo to Raster’ function in ArcGIS on a 1 km grid and smoothed over a circular area with a radius of 10 km (Hutchinson, (1988)) (Fig. 1b). The data window in which roughness calculations are made is between 10 and 320 m, meaning that roughness wavelengths within this frame are recorded well. However, under the interpolation procedure, it is inevitable that there is some loss of fidelity in the roughness. In our comparison between RES data and model output, we used the along-track roughness calculations.

The results reveal a significant zone of very low roughness in an area encompassing part of the IIS trunk and the BIR; there is no discernible change in bed roughness across the shear margin (Fig. 2a). Siegert and others (2016) concluded that this low roughness patch comprises unconsolidated weak sediment, and that the shear margin is controlled by thermo-dynamics and hydrology; where the bed of IIS is wet and the BIR is frozen. In contrast, the shear margin between BIR and MIS is marked by a change in bed roughness, which likely represents a border between the sediment and more undulating terrain, probably in relation to the Pagano Shear Zone that aligns parallel to the transition between fast and slow ice flow (Jordan and others, 2013).

2.2 BISICLES ice-sheet model and basal traction

The BISICLES ice-sheet model (Cornford and others, 2015) is a finite-volume model that applies a 2D force balance approximation to the solution of the Stokes free-surface problem, with the addition of a vertically integrated stress component (Schoof and Hindmarsh, 2010):

(4)$$\nabla. \lsqb {\emptyset h\bar{\mu} \lpar {2{\dot{\rm \varepsilon}} + 2{\rm tr}\lpar {{\dot{\rm \varepsilon}}} \rpar {\bf I}} \rpar } \rsqb + {\bf \tau} ^{\rm b} = \rho _{\rm i}gh\nabla _{\rm s},$$

where $\emptyset h\bar{\mu}$ is the vertically integrated effective viscosity, $\emptyset$ is a stiffening or ‘damage’ factor, h is the ice thickness and $\bar{\mu}$ is the vertically averaged viscosity. ${ \dot{\rm \varepsilon}}$ is the horizontal rate-of-strain tensor, ${\bf I}$ is the identity tensor, ${\bf \tau} ^{\rm b}$ is the basal traction, ρ _i is the density of ice, g is gravitational attraction and s is the ice surface. In the model, basal sliding resistance is governed by a linear viscous friction law, hence basal traction is defined as:

(5)$${\bf \tau} ^{\rm b} = \left\{ {\matrix{ {-C{\bi u}} & {{\rm if\;} \displaystyle{{\rho_{\rm i}} \over {\rho_{\rm w}}}h \gt \! -b} \cr 0 & {{\rm otherwise}} \cr}} \right.,$$

where C(x, y) is a basal friction coefficient, u is the horizontal ice velocity, ρ _w is the density of water and b(x, y) is the bed elevation. Inverse modelling of basal friction is ‘solved’ by optimising the fit between InSAR ice velocity (Fig. 1c, Supplementary Fig. S1c) and modelled ice speed (Fig. 1d, Supplementary Fig. S1d) to find both C and $\emptyset$ (Wright and others, 2014). A gradient-based optimisation method (non-linear conjugate gradient method) is used to quantify the difference between satellite-derived and modelled ice velocities. The basic optimisation is, however, underdetermined since both C and $\emptyset$ are unknown and sensitive to large fluctuations in response to noise in the observations. The percentage difference between measured and modelled ice surface velocities is provided in Figure 1e and Supplementary Figure S1e. The bed friction required for this solution is provided in Figure 1f and Supplementary Figure S1f. It should be noted that the model output used in this paper has a grid resolution of 2 km.

3. Model vs RES comparison

To our knowledge, aside from Siegert and others’ (2005) broad assessment of East Antarctic dynamics, there have been no analyses of how continental-scale numerical ice-sheet model results compare with RES-derived observations of basal roughness, although we note that Ashmore and others (2014) compared radar bed returned power to basal drag, Kyrke-Smith and others (2017) compared seismic-derived physical properties to basal drag, and Smith and others (2013) looked at a range of geophysical data to one another, including a basal drag inversion. The acquisition of RES data from IIS, MIS and FIS provides a good opportunity to assess if and how the ice-sheet model's depiction of basal conditions, necessary to optimise the calculation of surface ice velocity, agrees with the measurements of the ice–bed interface.

3.1 Institute Ice Stream and Bungenstock Ice Rise

To investigate the model's performance against basal measurements, five RES transects were compared with ‘synthetic’ model-derived transects (compiled from geo-rectified gridded model output) across the IIS and BIR (Fig. 1, Supplementary Fig. S1; Profiles A–A’, D–D’, E–E’ [0 to ~140 km], F–F’ and G–G’).

In general, the modelled ice speed accurately characterises the pattern of InSAR ice velocity as recorded in the RES transects, in terms of the magnitude and the gradients between slow and fast flow (orange and pink lines in Fig. 2; Profile A–A’, Supplementary Fig. S2; Profiles D–D’ to G–G’). There is also a good qualitative correlation between both the modelled and InSAR ice velocity across the BIR (Fig. 2; Profile A–A’, Supplementary Fig. S2; Profiles D–D’ to G–G’), and at the west of the IIS trunk (Supplementary Fig. S2; Profiles D–D’ and E–E’).

However, the modelled ice speed is noticeably lower than the InSAR ice velocity across the IIS trunk and towards the IIS-BIR shear margin (Figs 1e and 2, Supplementary Figs S1e and 2): Profiles A–A’ at ~80 km where the difference is ~100 km a⁻¹ (50% less); D-D’ at ~70 km where the modelled ice speed is up to ~200 m a⁻¹ (50%) slower than the InSAR-derived velocity; E–E’ at ~100 km where the difference is 100 m a⁻¹ (50% less); F-F’ at ~80 km where the difference is ~80 m a⁻¹ or (40% less); and G–G’ at ~60 km where the difference is ~40 m a⁻¹ (30% less). Conversely, the modelled ice speed is higher relative to the InSAR-derived ice velocity between 0 and ~60 km in Profiles D–D’ (up to 60 m a⁻¹, or a 20% increase).

The highest modelled velocities in the trunk of IIS correspond to a well-demarcated region over which the modelled basal traction is very low (Fig. 1f, Supplementary Fig. S1f). This low-drag region matches extremely well with the well-marked area of low roughness basal sediments (Siegert and others, 2016). Details of this match can be observed further in the transects at the western side of IIS (Supplementary Fig. S2; Profiles D–D’ and E–E’), where the slow modelled ice velocity due to the high basal friction/drag corresponds well with relatively higher total roughness. Hence, there is strong geophysical evidence for the association between enhanced ice flow, low basal drag and the presence of weak water-saturated sediment in and around the IIS trunk.

However, this relationship does not hold over the BIR (Fig. 2; Profiles A–A’ from 0 to 60 km, Supplementary Fig. S2; Profiles D–D’ from 80 km, F–F’ from 0 to 50 km; and G–G’ from 0 to 45 km), where both ice velocities and basal roughness are low (due to the extension of the sediment drape observed beneath the fastest component of the IIS). The model forces low velocities over BIR by introducing a region of high basal friction. The RES evidence suggests that the bed of BIR is frozen, as opposed to melting beneath the IIS trunk, according to reflected RES power (Siegert and others, 2016). Hence, the shear margin between IIS and BIR is a hydrological one, with otherwise identical basal material. Thus, the model's explanation of the IIS/BIR transition due to low/high bed friction simplifies the glaciological reality and limits understanding of how ice-sheet thermo-dynamics affect the ice sheet here.

3.2. Möller Ice Stream and Bungenstock Ice Rise

Six RES transects were analysed across the MIS and BIR to investigate how well the ice-sheet model is performing here (Fig. 1 and Supplementary Fig. S1; Profiles B–B’, E–E’ [~140 km onwards] and H–H’ to K–K’). InSAR ice velocity data across the MIS show that the ice stream is slower relative to both IIS and FIS. The model is able to capture the generalities of this ice stream feature very well (Figs 1c–e and Supplementary Figs S1c–e) and, unlike for IIS, the magnitude and position of modelled ice flow show a very good correlation with the InSAR ice velocity in all transects. The relatively slow satellite-derived ice speed is likely to be due to the apparent rough bed topography beneath the MIS trunk. This roughness can be seen in RES profiles (Fig. 2; Profile B–B’, Supplementary Fig. S2; E–E’ and H–H’ to K–K’). The ice-sheet model is able to account for lower ice velocities within the MIS trunk by increasing the basal friction. The association between the area of measured high bed roughness and modelled enhanced basal friction is strong, indicating geophysical support for the modelling approach to the MIS.

The MIS-BIR shear margin is dealt with in the model by a sharp transition between moderate and very high basal friction. As discussed above, the bed of BIR likely comprises extremely smooth frozen unconsolidated sediments. In addition, Jordan and others (2013) showed that the western margin of the MIS is aligned with the Pagano Shear Zone, which separates Jurassic intrusions from Cambrian/Permian metasediments. As the calculation of basal water flow indicates, the subglacial hydrology (after Le Brocq and others, 2009; Fig. 1b, Supplementary Fig. S1b) is greatly influenced by the Pagano Shear Zone and by subglacial thermo-dynamics. Hence, the geophysical explanation for the MIS-BIR border is a combination of the transition between frozen and wet beds, and a geological transition between smooth frozen sediment beneath the ice rise and rough bedrock beneath the ice stream. While the ice-sheet model is capable of matching measured ice velocities in MIS, the glaciological complexity revealed by geophysical analysis is not accounted for in the model.

3.3. Foundation Ice Stream

The ice-sheet model was run using Bedmap2 (Fretwell and others, 2013), which did not include much data for FIS. Hence, with new RES data now available (Jeofry and others, 2018a), FIS provides an excellent test for the ice-sheet model to examine whether the approach of optimising surface velocities holds in regions where bed topography is not depicted well.

Five RES transects were analysed across the FIS trunk (Fig. 1 and Supplementary Fig. S1; Profiles C–C’ and L–L’ to O–O’). While there are some consistencies between the satellite-derived and modelled ice velocities across the central region of the FIS trunk, the modelled ice speed is shown to be higher and lower than the InSAR ice velocity at the grounding zone and upstream FIS/Academy Glacier, respectively (Figs 1c and d and Supplementary Figs S1c and d). The modelled ice velocity shows reasonable correlation with the InSAR measurements across relatively rough subglacial terrain in all transects, suggesting the model's basal friction parameter has some geophysical support here. However, inconsistencies between the satellite-derived and modelled ice velocity occur across and upstream of the main trunk of FIS (Fig. 1e and Supplementary Fig. S1e, Fig. 2; Profile C–C’, Supplementary Fig. S2; Profiles L–L’, N–N’ and O–O’). In transect C–C’, the model generates an ice velocity that is ~100 m a⁻¹ less than the InSAR value (a reduction of ~25%), whereas in profile L–L’, the model velocity is ~100 m a⁻¹ higher than the InSAR value in the trunk of the ice sheet (15% increase), and over 200 m a⁻¹ higher (25%) at the shear margin (e.g. Supplementary Fig. S2; Profile L–L’ at between 20 and 40 km). Other differences of similar magnitude can be seen in profiles C–C’ and O–O’, where the InSAR velocity is ~20% higher than in the model.

The satellite-derived and modelled ice velocity is relatively low across the Pensacola Mountains and a large subglacial hill (from the Patuxent Range) at the boundary between the upstream FIS/Patuxent Ice Stream and Academy Glacier (Fig. 1e and Supplementary Fig. S1e; Fig. 2; Profile C–C’, Supplementary Fig. S2; Profiles L–L’ to O–O’). The absence of high qualitative reflectivity values across these regions, commonly used as a basic proxy for basal water (e.g. Siegert and others, 2016), indicates they are likely frozen at the bed. Thus, the model's introduction of high basal friction here can be justified by a combination of thermo-dynamics and rough terrain. The highest modelled velocities are produced as a consequence of low basal friction across the main trunk and upstream FIS. RES data show this region corresponds with a flat-based and deep subglacial trough (Fig. 2; Profile C–C’, Supplementary Fig. S2; Profiles L–L’ to O–O’). Hence, although the trough is not accounted for in Bedmap2, the model deals with this problem by reducing bed friction in regions where the trough is now known to exist.

While it is easier to diagnose the relationship between bed elevation and the modelled basal friction/drag across the main trunk and upstream FIS (Joughin and others, 2006), it is trickier across the Academy Glacier (Supplementary Fig. S2; Profiles N–N’ and O–O’). Bed elevation beneath the Academy Glacier is characterised by deep rough-bedded low elevation topography (Paxman and others, 2019). High modelled basal friction corresponds well with the position of the troughs despite their absence in Bedmap2. High RES reflectivity beneath Academy Glacier (Carter and others, 2007; Jordan and others, 2018), and scores of ‘active’ subglacial lakes (Wright and Siegert, 2012), suggest that there is a lot of water present at the ice–bed interface (Fig. 1a and Supplementary Fig. S1a). Hence, while the model depicts Academy Glacier reasonably well, it fails to account for known glaciological processes acting at the bed.

Given that basal topography is critical to the flow and form of FIS and its tributaries, and because the ice-sheet model was run with a bed that did not accurately represent this topography, it is remarkable to observe that the model accounts for FIS reasonably well through its optimisation of ice flow through adjustment to bed friction. Low friction is generally associated with the position of deep troughs and, while the exact location is often offset by a few kilometres, there does appear to be geophysical support for the approach. However, processes driving ice flow in FIS are simplified by the model, and thus changes to them (such as to hydrological conditions) are unaccounted for.

4. Summary and conclusions

This paper evaluates how the set-up of the BISICLES ice-sheet model (see Cornford and others, 2015) performs in realistically characterizing the flow and form of ice in the WS sector of the WAIS, with respect to gross ice flow patterns and basal conditions. The model was run using the Bedmap2 DEM, which incorporated crude in-field derived IMAFI ice thickness measurements, and without new CReSIS RES data in key areas (Jeofry and others, 2018a). Compared with the most recent DEM of the region Bedmap2 is known to be less realistic in characterizing topography, especially across FIS, which is influenced heavily by a ~2 km deep trough that extends upstream and connects with two distinct tributaries. Rather than modelling basal processes to determine enhanced ice flow, the model adjusts the amount of friction at the bed in order to match its ice flow with InSAR measurements.

Generally, the model characterises the complexities of enhanced ice flow across the WS sector well. However, the model's approach to simulating enhanced ice flow, and shear margins between slow and fast flow, solely by adjustment of basal drag is unsophisticated and misses important glaciological processes associated with subglacial topography, geology and ice-sheet thermo-dynamics. For example, RES demonstrates that IIS, MIS and FIS each have different glaciological regimes that are unaccounted for by simply adjusting bed friction. Hence, while the model describes the flow of ice reasonably well, its use of bed friction paramaterisation to achieve this overlooks processes that may be subject to change, either dynamically (such as the influence of active subglacial lakes causing changes to basal water pressure) or through modification to ice-sheet systems when the ice sheet responds to external drivers (such as thermo-dynamics).

Re-running the model across the new bed DEM (Jeofry and others, 2018a) is clearly worthwhile. It may also be beneficial to match the areas of enhanced flow, and its borders, against geophysical evidence in order to evaluate where reduced bed friction is most relevant (such as over wet-based sediments) and where other processes unaccounted for in the model (such as hydrology) may be influential. In doing this, the model may be better able to determine not only where the ice sheet is most sensitive to external drivers but why change may occur. The adaptive mesh function in the BISICLES model may also need to be investigated with respect to RES data, as the resolution of the model output will vary across the ice sheet. Our work highlights the benefit of examining RES along-track data to understand the function of ice-sheet models relative to the ice sheets they are aiming to mimic. Until now, surprisingly little model-data comparison has been undertaken, which suggests a great deal of future work is possible and necessary.