Introduction

Nearly two decades of satellite observations have revealed that the flow of Pine Island Glacier (PIG), West Antarctica, has accelerated, along with other glaciers draining into the Amundsen Sea, causing rapid mass loss (e.g. Reference Shepherd, Wingham and MansleyShepherd and others, 2002; Reference RignotRignot and others, 2004; Reference Wingham, Wallis and ShepherdWingham and others, 2009; Reference Joughin, Smith and HollandJoughin and others, 2010). Oceanic influence has played a key role in driving the rapid mass loss. The temperature and volume of warm Circumpolar Deep Water (CDW) in Pine Island Bay increased between 1994 and 2009 (Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others, 2011). The resulting increase of sub-ice-shelf melting to several tens of meters per year has reduced ice-shelf buttressing and allowed faster flow of PIG (Reference Payne, Holland, Shepherd, Rutt, Jenkins and JoughinPayne and others, 2007). Although the current phase of thinning of PIG is suggested to have begun before 1973 when the ice was clearly in contact with a submarine ridge (Reference JenkinsJenkins and others, 2010), the grounding-line migration and flow-speed changes are poorly quantified before the satellite era. The delivery of oceanic heat to the base of the ice shelf depends on its shape and on the bathymetry of the sub-ice-shelf cavity, which are poorly known except along a few survey lines that were collected by an autonomous underwater vehicle (AUV) in 2009 (Reference JenkinsJenkins and others, 2010). Models of the future of PIG also depend on bed properties (e.g. Reference HookeHooke, 2005; Reference JoughinJoughin and others, 2009). Mixed-bed conditions, with both hard crystalline rock and sediments, have been observed for the ice-free areas in Pine Island Bay (Reference Lowe and AndersonLowe and Anderson, 2003) and inferred for the base of the main trunk of PIG through a diagnostic modeling study (Reference JoughinJoughin and others, 2009). However, the bed morphology and sediment distribution beneath the PIG ice shelf that influenced the evolution to the current behavior and imbalance remain uncharted.

Acoustic techniques (seismic or sonar) offer the most direct means of assessing depths of sediment and water, but application is difficult and expensive for the cavity beneath a heavily crevassed ice shelf such as that of PIG. The gravity technique has been widely used for estimating geological structures of various spatial scales such as subsurface rock density variations (e.g. Reference Li and OldenburgLi and Oldenburg, 1998; Reference Camacho, Fernández and GottsmannCamacho and others, 2011) and crustal thickness (e.g. Reference JordanJordan and others, 2010). Glaciological application of the gravity technique includes estimation of the depth and distribution of subglacial sediment and water bodies (e.g. Reference StudingerStudinger and others, 2001; Reference Filina, Blankenship, Thoma, Lukin, Masolov and SenFilina and others, 2008) and sub-ice-shelf bathymetry (e.g. Reference Tinto and BellTinto and Bell, 2011). Gravity data lack the resolution of seismics, but are already available for key regions of the Antarctic ice sheet (Reference Tinto and BellTinto and Bell, 2011; Reference Cochran and BellCochran and Bell, 2012) and provide a powerful tool for interpolating between seismic data, as well as allowing lower-resolution estimates even where seismic data are lacking.

In this paper, we analyze aerogravity data using a non-linear inversion approach constrained by AUV data to estimate the sub-ice-shelf bathymetry and distribution of sediment for PIG. Although the non-uniqueness of such inversions (e.g. Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005) and the lack of detailed knowledge of the in situ density structure (e.g. Reference Studinger, Bell and TikkuStudinger and others, 2004) prevent precise determinations, the data do provide useful guidance on key sub-ice-shelf features.

Methods

Data

We used the free-air gravity anomaly data (70s full wavelength filter) acquired during the 2009 Antarctic campaign of NASA’s Operation IceBridge (Reference Cochran and BellCochran and Bell, 2010). The mean crossover error in this dataset over the PIG area is 1.2 mgal. We also used the ice and ocean surface elevations measured by the Airborne Topographic Mapper (ATM; Reference KrabillKrabill, 2010) and ice thickness derived from the Multi-channel Coherent Radar Depth Sounder (MCoRDS; Reference AllenAllen, 2010), both of which were obtained concurrently with the gravity data. All three datasets were interpolated to a 2.5 km × 2.5 km grid covering the area shown in Figure 1b. The coverage of the gravity data is relatively sparse away from the ice shelf; hence, we focus our analyses and discussions on the area over, and close to, the ice shelf, although the inversion is performed over the entire area shown in Figure 1b to avoid edge effects. We used the seabed depths along the AUV tracks shown in red lines in Figure 1b (personal communication from A. Jenkins, 2011) as fixed constraints in the inversion. Reference JenkinsJenkins and others (2010) created a map of seabed elevations by contouring those AUV data (Reference JenkinsJenkins and others, 2010) and we used it as the starting model of the bathymetry for the inversion together with an initial sediment layer thickness of 200 m.

Fig. 1. (a) Map of the study area. Image is from Moderate Resolution Imaging Spectroradiometer (MODIS) Mosaic of Antarctica (Reference Haran, Bohlander, Scambos, Painter and FahnestockHaran and others, 2006). (b) Interpolated map of the IceBridge free-air gravity anomaly. The thin black lines show flight paths where gravity data were collected, the thick black line is the grounding line from MEaSUREs, the thick pink line is the edge of the ice shelf visually picked from a MODIS image on 9 November 2009 obtained from the US National Snow and Ice Data Center’s MODIS Antarctic Ice Shelf Image Archive (Reference Scambos, Bohlander and RaupScambos and others, 2009) and the red lines are the AUV tracks.

Gravity inversion

For our inversion of the free-air gravity anomalies for bathymetry and sediment layer thickness, our forward model followed Reference PlouffPlouff (1976). Here, the domain of interest is discretized into small three-dimensional (3-D) rectangular prisms with four layers of different densities: ice (p = 900 kgm–³), sea water (p= 1028 kg m–³), sediment (2013 kg m–³, appropriate for rock of 40% porosity with voids filled with sea water) and bedrock (2670kgm–³), as shown in Figure 2. We included a sediment layer following the inference of Reference JoughinJoughin and others (2009) for the existence of extensive weak bed under the grounded part of PIG, due most likely to deforming sediments, and the widespread occurrence of sediments in deglaciated regions (e.g. Reference Lowe and AndersonLowe and Anderson, 2003). We have no detailed knowledge of the appropriate density for the sediment; hence, we chose porosity of 40%, which leads to a density close to values previously used by other gravity studies for Antarctica (Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005; Reference Filina, Blankenship, Thoma, Lukin, Masolov and SenFilina and others, 2008).

Fig. 2. Schematic diagrams of (a) 3-D representation of four-layered source body by rectangular prisms (adapted from Reference Seber, Sandvol, Sandvol, Brindisi and BarazangiSeber and others, 2001; Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005) and (b) discretization of the domain of interest (from Reference Seber, Sandvol, Sandvol, Brindisi and BarazangiSeber and others, 2001).

The 3-D rectangular prisms used in the discretization have fixed horizontal dimension of 2.5 km x 2.5 km. The depth of each layer interface was allowed to vary, except for ice (z₁ and z₂ in Fig. 2a, available over the whole domain from the ATM and MCoRDS data) and water (z₃ in Fig. 2a) along the AUV tracks. The 3-D gravity anomaly of a prism is calculated by the formula of Reference PlouffPlouff (1976), and the total gravity anomaly at a given observation point P(_m,n) is computed as the sum of the contributions of gravity anomalies from all the prisms. The gravity anomaly for the whole domain is then given by executing the above calculation for all observation points.

The observed and calculated gravity anomalies are offset by a DC shift (Reference Tinto and BellTinto and Bell, 2011), since the method introduced above calculates the absolute gravity anomaly whereas the IceBridge gravity data are in reference to the geoid. Since we have many observations of seabed depth scattered in 2-D space, we determined this DC shift by searching for the plane that minimizes the root-mean-square (rms) difference between the observed and calculated gravity anomalies along the AUV tracks.

We performed the inversion using very fast simulated annealing (VFSA), which was employed in modeling the water depth and the sediment layer thickness of Vostok Subglacial Lake by Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others (2005) and Reference Filina, Blankenship, Thoma, Lukin, Masolov and SenFilina and others (2008). Details of VFSA are provided by Reference Sen and StoffaSen and Stoffa (1995) and Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others (2005). VFSA is preferred here despite being computationally slower than some other inversion methods such as the Fourier transform-based algorithm (Reference ParkerParker, 1973; Oldenberg, 1974; our large number of model parameters requires a large number of iterations to find the optimal solution). However, simulated annealing (SA) is especially powerful when the relation between model and data is nonlinear (Reference Sambridge and MosegaardSambridge and Mosegaard, 2002), as in our case, and also makes it relatively easy to constrain the inversion with a priori information such as the AUV tracks (Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005). Also, the VFSA version of SA rapidly tests thousands of models, which allows us to derive the posterior probability density (PPD) for Bayesian statistics (Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005). We ran the VFSA algorithm five times to test 5000 models. The mean of the PPD is our most likely model and we report uncertainties as the 95% confidence interval (approximately two standard deviations).

We solve the inverse problem by minimizing the misfit function:

(1)

where g is the free-air gravity anomaly and subscripts obs and cal denote the observed and calculated gravity anomaly respectively. The normalization in Eqn (1) prevents large local anomalies from dominating the inversion through the tendency for larger anomalies to have larger misfits (Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others, 2005). The general technique described by Reference Roy, Sen, Blankenship, Stoffa and RichterRoy and others (2005) includes the option of fitting spatial gradients in gravity anomalies as well as the anomalies themselves, which is useful for very sharp gradients such as those associated with faults, but we follow guidance for relatively smooth variations and consider only the anomalies. Sub-glacial topography derived from aerogravity data shows that PIG flows in a graben-like valley (Reference JordanJordan and others, 2010), which may be better resolved by incorporating the spatial gradient in gravity anomalies. However, it is not clear at this point if such a subglacial valley extends to the area of the PIG ice shelf, and without further knowledge we prefer relatively simple results (smoothly varying bathymetry). We evaluated E for the grid locations within the target area shown in Figure 1b that are within 5 km (two gridcells) of the nearest point along the aerogravity flight-lines.

To speed up the forward calculation of the gravity anomaly, we calculated the gravity effect of only those prisms within a 20 km radius of a particular point. We found after several exploratory computations that the contribution of the blocks beyond 20 km can be neglected because it is no larger than ∼0.12 mgal versus the ∼1.2 mgal error in the IceBridge free-air anomaly data. We further reduced the computational time for convergence to the optimal model by using a smoothing filter when perturbing the model, employing a 5 :1 weighting of a prism with the eight nearest neighbors in a layer.

Because of the inherent non-uniqueness of the gravity inversion (e.g. Reference Jacoby and SmildeJacoby and Smilde, 2009), it is best to introduce as much a priori information as possible to constrain the problem. In addition to the AUV bathymetry data, we used the grounding-line position from the NASA Making Earth Science Data Records for Use in Research Environments (MEaSUREs; Reference Rignot, Mouginot and ScheuchlRignot and others, 2011). We also fixed the sediment layer thickness to zero for the Hudson Mountains nunataks, which are eroded volcanoes (Reference Rowley, Laudon, La Prade, LeMasurier and LeMasurierRowley and others, 1990).

Resutls and Discussion

Figure 3a shows the modeled bathymetry (the depth below sea level of sea water where it occurs and of ice where no sea water is present), with the residual between the observed and calculated free-air gravity anomaly shown in Figure 3b. Water column thickness (measured between the ice-shelf base and sediment or from sea level to sediment in front of the ice shelf) and its 95% confidence interval are shown in Figure 3c and d, respectively, with sediment thickness and its 95% confidence interval in Figure 3e and f, respectively. Although not all parameters are normally distributed, we found that the majority of PPDs can be approximated by the Gaussian distribution, motivating our use of the 2σ values here.

Fig. 3. (a) Modeled bathymetry. The white dashed line indicates the cross-section profile shown in Figure 4, and the black dashed circle indicates the submarine ridge. (b) Residual between the observed and calculated gravity anomaly. (c) Water column thickness. (d) 2σ uncertainty for bathymetry and water column thickness. (e) Sediment layer thickness. (f) 2σ uncertainty for the sediment layer thickness.

The mean absolute residual for the target area is 5.6 mgal, just over four times the error in the gravity data. However, Figure 3b shows that the fit is much better than this in almost all of the target area. Large residuals in excess of +20 mgal (shown in dark red) are widespread outside our target area.

Of greater interest are localized residuals of between -15 and -20 mgal in a few places near the grounding line. We suggest that these arise from the difference between the grounding line used in our inversion and the grounding line at the time of the gravity acquisition. We constrained the grounding line in our model by using the MEaSURES differential interferometric synthetic aperture radar (InSAR) data collected from 1996 to 1999, but used the IceBridge gravity and ice-thickness data from 2009, and the grounding line likely retreated between these measurement campaigns as shown by Reference Joughin, Smith and HollandJoughin and others (2010). The ice ungrounded from the submarine ridge in the 1970s (Reference JenkinsJenkins and others, 2010). Subsequently, the ice thinned at a rate of 5.5 m a^-1 between 1992 and 2001 (Reference Shepherd, Wingham and RignotShepherd and others, 2004), accompanied by grounding-line retreat of 1.2 km a^-1 along the main trunk between 1992 and 1996 (Reference RignotRignot, 1998). Thinning has continued at a higher rate in excess of -7 m a^-1 between 2003 and 2008, as observed by Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others (2012), so the grounding line likely has continued to retreat. By forcing the model to include ice grounded to the 1990s position using 2009 ice thickness, we have essentially forced the model to put sediment rather than water into the region crossed by the retreating grounding line. Indeed, the inversion returns anomalously thick sediment in those regions of negative residual gravity anomaly near the grounding line; replacing some of this sediment with lower-density water would greatly reduce the residuals. We are not promoting aerogravity as an optimal means to map grounding-line retreat; however, the ability of the IceBridge data and the inversion to detect the effects of the migrating grounding line may be of interest to readers who are not familiar with the sensitivity and fidelity of the technique.

One additional large negative residual occurs away from the grounding line near the Hudson Mountains (coordinate -1625, -270; -21 mgal), close to large positive residuals (near coordinate -1632, -277; up to +23 mgal). In this region of steep topography and complex volcanic geology (Reference Rowley, Laudon, La Prade, LeMasurier and LeMasurierRowley and others, 1990), even the flight-line spacing of 5 km longitudinal and 10 km latitudinal to the flow direction of PIG may be too sparse. Furthermore, the choices we made for an optimal inversion beneath the ice shelf (such as 2.5 km x 2.5 km prisms, a uniform bedrock density and no faults) are likely to be suboptimal here.

Figure 3d and e show that uncertainties in the thickness of the water column and sediment layer are generally in the range ±150-180 m, but as large as ±220 m for the latter. As expected, uncertainties are larger for areas away from the AUV track. Figure 3a shows that the AUV tracks detected the southern part of a deep basin near the modern grounding line of PIG, deepening up-glacier. A second deep basin mapped by the AUV near the modern ice-shelf front also has a second lobe, to the south. Uncertainties in water column thickness (Fig. 3d) are sufficiently small that these basins are unlikely to be artifacts. The grounding line is strongly topographically controlled, extending seaward around bathymetric highs.

The previously unmapped deep basin near the modern grounding line (coordinate -1605, -280) has similar bathymetry (∼700m deep; Fig. 3a) and water column thickness (∼400-500m thick; Fig. 3c) to the basin detected by the AUV (coordinate -1585, -285). It is possible that the characteristics of the sea water in this newly mapped basin are similar to those measured by AUV in the basin to the south, where the water temperature was found to be ∼3.5°C above the in situ freezing point (Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others, 2011). If so, then the recent ice thinning and grounding-line retreat caused by the access of warm CDW, as suggested by Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others (2011), likely occurred across a large part of the ∼30km wide main trunk of PIG.

The optimum model shows significantly nonzero sediment thickness everywhere except over some of the volcanic peaks of the Hudson Mountains that were assumed to have no sediments. The sediment layer is ∼700 m thick along the AUV track near the grounding line, with similar or thicker sediment seaward of the ice-shelf front. Thinner sediment is indicated further up-glacier beneath PIG and in association with the bathymetric ridge under the ice shelf. Recall that sediment near the grounding line is likely thinner than shown due to the effects of the retreating grounding line, as discussed above.

The sediment layer is well resolved along the long AUV track (Fig. 4a). The PPD of the model is shown with the 95% confidence interval as the uncertainty bounds; the optimal model of Figure 3 falls in the middle of the dark band of the PPD. The PPD is also shown for the one gridpoint between the end of the AUV track and the MEaSURES grounding line, at 67.5 km from the seaward end of the AUV track where the bathymetry had to be obtained from the inversion. The neighboring constraints led to a skewed PPD. The sediment layer thickness at the seaward end of the profile is 5 3 8 ± 138 m and thins towards the seaward slope of the ridge. The sediment layer is thinnest (193 ± 120 m) 22.5 km from the seaward end of the profile. The thickness then increases towards the grounding line, reaching 479 ± 143 m on the ridge crest and a maximum of 7 3 4 ± 189m at ∼10 km in front of the MEaSURES grounding line.

Fig. 4. (a) Cross section of the model along the long AUV line (indicated by white dashed line in Fig. 3a) and the PPD function. The 95% confidence interval of the PPD is indicated by dashed lines. The red dashed line is the depth of the sediment layer manually perturbed to decrease the residual between the observed and calculated gravity anomaly. (b) Observed and calculated free-air gravity anomalies (solid lines) and the gravity anomaly calculated with the manually perturbed sediment layer depth indicated by the red dashed line in (a).

The gravity anomaly modeled on the long AUV track has residuals of up to ± 7 mgal (Fig. 4b). We were able to reduce the residuals further by manually perturbing the sediment layer thickness proportional to the residual, decreasing by ∼150m at ∼27km from the start of the profile where the residual was -7 mgal and increasing by 100-200 m between 65 and 75 km points where the residual was +7 mgal. This shifted the site of thinnest sediment ∼5 km upslope toward the grounding line and the thickest sediment ∼10 km inland without notably changing the main features of the sediment distribution.

We note that with the data and constraints available, it is difficult with the gravity technique to separate recently deposited tills from more lithified, older sedimentary deposits. Thus, our results are not necessarily in disagreement with the suggestion by Reference JenkinsJenkins and others (2010), based on the bedforms of the sea floor imaged by the multi-beam acoustic sensor of an AUV, that the sediment layer beneath the PIG ice shelf is thinnest or is absent at the crest of the ridge, is thicker on the seaward slope from deposition of sediments scoured off the ridge crest, and is thin or absent again near the edge of the ice shelf. We simply note that if Reference JenkinsJenkins and others (2010) are correct, then the distribution of the most recent sediments differs from that of older sediments, perhaps suggesting that erosion beneath advancing ice and deposition by retreating ice have both been concentrated in the lee of the topographic high.

Conclusion

Gravity, laser and radar data from NASA’s Operation IceBridge, constrained by AUV measurements and grounding-line positions from InSAR, allow inversion for the thicknesses of the water column and sediment layer beneath the PIG ice shelf and surrounding areas. Relatively deep water is found adjacent to, but distinct from, basins mapped by the AUV. The newly mapped deep basin near the modern grounding line likely indicates that the recent ice thinning and grounding-line retreat caused by intrusion of warm CDW occurred across most of the width of the main trunk of PIG. Significant thicknesses of sediment are resolved everywhere except over a few volcanic peaks of the Hudson Mountains that were assumed to have no sediments. In light of the generally greater erodibility and smoothness of sediments and sedimentary rocks relative to crystalline bedrock, it is likely that over recent glacial–interglacial cycles PIG was able to generate lubricating till, much like Whillans Ice Stream on the Siple Coast (Reference Alley, Blankenship, Rooney and BentleyAlley and others, 1989), and was not restrained by rough bedrock.