Pine Island Glacier (Antarctica) velocities from Landsat7 images between 2001 and 2011: FFT-based image correlation for images with data gaps

Abstract

Remote sensing of ice motion by tracking displacement of surface features is a valuable tool in glaciology. Efficient image feature-tracking programs, such as IMCORR, based on fast Fourier transform methods can produce misleading correlations if there are data gaps in either or both of the reference and search images. This is particularly problematic if the data gaps are regular in character, such as for Landsat7 images collected after the failure of the Scan Line Corrector (SLC-off images). We demonstrate that this situation can be alleviated by filling the data gaps with suitably chosen random data. We modified IMCORR to achieve this automatically (source code is included), but generic image-processing software could be used to modify inputs for other correlation packages. We test our method using images of Pine Island Ice Shelf, Antarctica, and document the acceleration of the velocity field for the floating extension of Pine Island Glacier over the decade 2001–11. We also combine our velocities with recent NASA Operation IceBridge ice thickness data from CReSIS to estimate the basal melt rates.

1. Introduction

The estimation of ice-sheet, ice-shelf and glacier surface velocities from optical imagery relies on calculation of the displacement of features between pairs of co-registered images (Bindschadler and Scambos, 1991). Automated methods for evaluating the displacement between subregions of the image pairs rely on methods for calculating the optimal correlation between numerous potential optimal subregion image pairs. Automated techniques can recognize the displacement of coherent and persistent patterns as well as distinctive individual features. Various methods have been applied to this problem, such as spatial domain correlation (e.g. Ahn and Howat, 2011), least-squares matching of digital terrain models (e.g. Kaufmann and Ladstädter, 2003), orientation correlation (e.g. Haug and others, 2010) and frequency domain correlation (e.g. Scambos and others, 1992). Frequency domain methods offer substantial computation efficiency gains compared to spatial domain methods, especially when fast Fourier transform (FFT) methods are employed (Haug and others, 2010).

The Landsat7 satellite was launched in April 1999, and as of January 2013 continues to provide valuable data from the Enhanced Thematic Mapper Plus (ETM+) instrument including 15 m spatial resolution panchromatic band imagery. Landsat7 imagery is of particular value for displacement tracking due to its relatively high spatial resolution, long history and associated archive of many images. The large area covered by an individual image makes matching efficient and enables co-registration using reference points even with extensive regions of moving ice. These advantages have led to a large body of work and expertise involving image correlation using Landsat7 images.

On 31 May 2003 the Scan Line Corrector (SLC) in the ETM+ instrument, which compensated for the forward motion of the satellite, failed. This resulted in the duplication of imaging of some subregions of each scene and the failure to acquire images of other portions. The data gaps are aligned normal to the direction of satellite motion (i.e. cross-track), and become more pronounced towards the edges of the image, since the satellite travels furthest in the period before the scan returns to each cross-track edge. Only a strip through the centre of the image is essentially complete, while the gap at the extreme edges is approximately 390–450 m, corresponding to one full scan line. Typically 22% of a full scene is lost. Figure 1 shows an example of the combined data-gap masks from a pair of Landsat7 images of Pine Island Glacier and Ice Shelf, Antarctica (location shown in inset). In the sub-images we use for testing our modifications to IMCORR in Section 3, the data gaps in the masks have a maximum width of 360 m, with the combined mask gaps spanning up to 540 m in Figure 1.

Fig. 1. Superposition of the Landsat7 ETM+ SLC-off data-gap masks from Landsat images of the Pine Island Ice Shelf, Antarctica (location shown in inset), for 8 December 2010 and 11 February 2011 that were applied to the 13 January 2001 and 16 November 2002 images to assess the impact of SLC-off mask data gaps on visible feature correlation tracking. The 1999 grounding line from Rignot and others (2011a) is shown in grey, and the 2009 grounding line, after Joughin and others (2010), is shown in black.

These data gaps have limited the use of these so-called SLC-off images in surface displacement studies, although two previous studies have developed strategies to allow the use of Landsat7 SLC-off images to determine surface feature displacements. Heid and Kääb (2012) include a review of correlation techniques for Landsat7 SLC-off images. Haug and others (2010) use frequency correlation applied to normalized complex images constructed from spatial derivatives of the image intensity, termed orientation correlation by Fitch and others (2002). The complex fields are set to zero in regions where the image is uniform and has vanishing derivatives. This reduces the sensitivity to data gaps as only their edges produce any signal. Ahn and Howat (2011) overcome the data-gap problem by performing the correlation in the spatial domain, and deliberately excluding any data-gap points from the relevant correlation sums. The current method, described in Section 2, takes a different approach, that retains the computational advantages of operating in the frequency domain and utilizing FFT algorithms, by adapting existing well-tested software for image intensity correlations.

Pine Island Glacier (PIG), its floating continuation which passes through the southern portion of the Pine Island Ice Shelf (PIIS) into Pine Island Bay, and the wider Amundsen Sea sector of the West Antarctic ice sheet are much studied (e.g. Rignot, 1998; Shepherd and others, 2001; Rignot and others, 2002; Rignot, 2008; Jenkins and others, 2010; Jacobs and others, 2012; Lee and others, 2012), due in part to the substantial acceleration of the PIG system since the period 1974–87 (Joughin and others, 2003), an acceleration that has been interspersed with periods of relatively stable velocities (Joughin and others, 2003; Rignot, 2008; Lee and others, 2012).

The PIG acceleration has been driven by enhanced Circumpolar Deep Water (CDW) intrusions onto the continental shelf and into the ocean cavity beneath the ice shelf (Shepherd and others, 2004), causing thinning of the buttressing ice shelf and retreat of the grounding line (Rignot, 1998, 2002; Jenkins and others, 2010; Joughin and others, 2010). The speed-up of PIG has resulted in drawdown of the surrounding grounded ice sheet (Shepherd and others, 2001, 2004; Pritchard and others, 2009; Scott and others, 2009; Wingham and others, 2009) and a significant loss of ice from the region as observed from analysis of the Gravity Recovery and Climate Assessment (GRACE) satellite gravity data (see, e.g., Lee and others, 2012).

Pritchard and others (2012) present elevation changes over Antarctic ice shelves, as measured by Ice, Cloud and land Elevation Satellite (ICESat) laser altimetry, and, after applying modelled corrections for elevation change due to non-steady accumulation and firn densification, estimated that the ice thickness of the southern part of the PIIS decreased at 6.8 m a⁻¹ between 2003 and 2008, which they attribute to increased basal melting.

Basal melt underneath ice shelves is driven by heat exchange with the underlying ocean, and recognition that it is a major loss mechanism for Antarctic ice shelves and major outlet glaciers is relatively recent (e.g. Rignot and Jacobs, 2002). Remote sensing and oceanographic measurements both indicate that averaged melt rates under the PIIS are large, with much higher values near the PIG grounding line where the ice shelf is thicker, and that the rates have been increasing over recent years. The increases appear to be associated with greater circulation of CDW in the evolving sub-ice-shelf cavity, as well as small increases in ocean temperatures. For example, Rignot (1998) estimated melt rates under the floating part of PIG averaging 50 ± 10 m a⁻¹ over the first 20 km from the grounding line, and a mean rate between the grounding line and the ice-shelf calving front of 24 ± 4 m a⁻¹, during the period 1992–96 when the grounding line was retreating at 1.2 ± 0.3 km a⁻¹. Using oceanographic measurements, Jacobs and others (2011) estimated that meltwater production for the main (southern) component of the PIIS was 85 ± 6 km³ a⁻¹ in 2009, an increase of ∼50% since 1994 (56 ± 7 km³ a⁻¹of melting then, applying the same analyses to sparser observations). They inferred corresponding melt rates of 33 and 22 m a⁻¹, if all the losses were attributed to the main fast-flowing PIG stream (assuming an area of 2500 km²), although they pointed out that basal melting patterns are expected to be far from uniform. They also noted that small increases in ocean temperatures over that period were probably insufficient to explain the increased melting and that changes in the sub-ice cavity geometry and circulation were likely to be significant. In a study of variability of basal melt beneath the PIIS, Bindschadler and others (2011) also estimated ‘non-seasonal’ basal melt rates for four profiles along the floating extension of PIG. For ice-shelf thicknesses and velocities corresponding to the 2004–05 epoch they estimated an average basal melt rate of 19 m a⁻¹ over the flowband from the PIG grounding line to the calving front, from which they estimated a loss rate of 44.2 km³ a⁻¹ for the PIG stream. For the zone nearer the grounding line they found considerably higher melt rates locally, with two profiles indicating 78.2 and 92.7 m a⁻¹ over the first 18 km downstream of the grounding line. Bindschadler and others (2011) also interpreted a sequence of basal voids in one of their profiles as reflecting a seasonal process of additional melting near the grounding line. They correlated this with ocean modelling results for summer ocean heat content in the Amundsen Sea region by Thoma and others (2008), suggesting it indicated an influence of interannual variations in ocean forcing.

Rignot (2008) presented estimates of ice velocities in PIG and the PIIS from a variety of epochs, and argued that the ice flux at the PIG grounding line increased from 77 ± 2 Gt a⁻¹ in 1996, to 100 ± 5 Gt a⁻¹ in 2006. He suggested that scaling the measured 2006 fluxes in proportion to general changes in grounding line ice velocities provides an estimate of changes in ice discharge from PIG into the southern part of the PIIS, and estimated the flow in 2007 would be 107 ± 15 Gt a⁻¹. Clearly this line of reasoning assumes that changes in grounding line ice thickness can be neglected.

Previous studies of the PIG velocity field, both on the ice shelf and the upstream grounded glacier, include the use of both radar (Rignot and others, 2002; Joughin and others, 2003; Rabus and Lang, 2003; Rignot, 2008) and visible wavelength (Joughin and others, 2003; Lee and others, 2012) imagery, and ground-based GPS measurements (Scott and others, 2009).

Simple modifications to an existing, widely used software package (IMCORR; Scambos and others, 1992) are presented here, which extend the analysis of displacement of surface features between image pairs to include the use of images with data gaps. The method is then applied to a sequence of Landsat7 images of the PIIS, where later images have large data gaps due to the SLC failure.

The PIG system was chosen as a demonstration site for this new technique since it is an area of considerable interest, as shown by numerous previous studies discussed above. The Landsat7 image archive offers a long sequence of images of the PIG region, all from the same instrument, and correlation techniques that cope with the SLC-off data gaps offer the capacity to provide additional information about this important and dynamic region.

2. Method

The presence of regular and hence highly correlatable data-gap artefacts in images creates difficulties for automated correlation techniques, and this is especially the case for Landsat7 SLC-off images. High-efficiency spectral or frequency domain array correlation methods using FFT algorithms require regularly spaced data without gaps (Swan, 1982). To apply these methods to imagery with data gaps that can be defined and masked (e.g. the SLC-off data voids), we fill the data-gap regions of both the respective reference and search sub-images with random data drawn from the intensity histogram of the corresponding sub-image (calculated excluding the data-gap value of 0). The added random data should have no consistent correlations with the other image, and so should not influence correlations from surrounding valid data.

While this infilling could be undertaken as a separate image pre-processing step, for convenience we implement it as part of the IMCORR (Scambos and others, 1992) software package by modifying the IMCORR source code (see Appendix for details).

For the purpose of validating the randomized data infilling technique, we apply it to Landsat7 images of the PIIS that pre-date the SLC-off problem to generate a reference measurement, and compare the results when SLC-off masks from corresponding later images are applied to the pair of images. Specifically 97.5 × 97.5 km sub-images from the Landsat7 ETM+ images from path 233 row 113 (see inset in Fig. 1 for location) acquired on 13 January 2001 and 16 November 2002 were used. These sub-images contain the main trackable features of interest. These images, and those used later to analyse acceleration of PIG and the PIIS are ‘Systematic Terrain Correction’ (Level 1 GT) products in standard Polar Stereographic projection (reference latitude 71° S) from the US Geological Survey.

The band 8 images were contrast-enhanced to strengthen correlations. They were then high-pass filtered by subtracting a smoothed image that was generated by convolution with a Gaussian kernel using a characteristic scale of 177 m for 50% suppression. This filtering removes smooth, longer-scale variations in intensity that likely reflect the underlying bedrock as expressed in ice surface topography and are typically stationary rather than reflective of ice motion.

To allow high-pass filtering of the SLC-off images, without the data gaps causing undesirable artefacts extending into the valid portions of the image, the data-gap mask for the image was extracted and the data gaps filled using nearest- neighbour interpolation. The image was then high-pass filtered as above, and finally the data-gap mask was reapplied. This was only to stop the data-gap regions from influencing the high-pass filtering, as the data-gap regions were reset to zero and subsequently filled with suitable random data by the modified version of IMCORR.

For the validation tests, the SLC-off data-gap masks from the corresponding regions of 8 December 2010 and 11 February 2011 Landsat7 images were superimposed on the 13 January 2001 and 16 November 2002 images respectively. Visible feature displacements between the 2001 and 2002 images were then calculated using the modified version of IMCORR (for both the reference unmasked and the masked image pairs) and the unmodified version of IMCORR for the masked image pairs. The parameters used in all these IMCORR analyses were: a reference window size of 64 × 64 pixels (corresponding to a 960 × 960 m window at the 15 m resolution of the band 8 Landsat7 image), a search window of 512 × 512 pixels to cover a large range of displacements, and displacement offsets of 70 and −160 pixels in the image easting and northing directions, respectively. This choice of reference window size means that it always contains some potentially valid data in spite of the SLC-off mask. We take a 50% overlap between reference windows, corresponding to a 32-pixel stride, to generate a displacement field on a 480 × 480 m grid.

The uncertainty in displacement measurements from IMCORR is 0.1 pixels (Scambos and others, 1992) corresponding to 1.5 m for Landsat ETM+ band 8 imagery, although the influence of orthorectification and image registration errors can lead to considerably larger systematic errors. Orthorectification, the task of unfolding the influence of surface topography in converting the image to a map plane, relies on elevation models, and errors in these models can artificially displace features, including potential no-motion registration points. This is a generic issue in feature tracking and one that should have limited influence on the present study, since the area we study has relatively low surface relief. We deliberately selected Landsat image pairs separated by large time intervals. While this reduces the ability to track shorter-lived features, it does result in large displacements (typical values are ∼100 pixels in easting and ∼280 in northing directions). Accordingly, even pixel-level errors in the displacements in the present study would typically lead to small velocity errors, at the 1% level or below.

3. Validation using images of the pine island ice shelf

The SLC-off data-gap masks correspond to 12.2% and 12.8% of the subimages, and, due to slight overlap in these masks, the combined mask accounts for 22.4% of the sub-image pairs (Fig. 1). While the presence of the SLC-off data-gap mask clearly results in fewer valid IMCORR correlation matches (see Fig. 2) our modification to IMCORR for filling the data gaps with suitable random data clearly mitigates this degradation to a significant extent. We found that adjusting the threshold of the IMCORR parameter ‘correlation strength’ to (subjectively) eliminate most spurious correlation matches, while maintaining the vast majority of the valid correlation matches in the spatially coherent regions, suggests a threshold in correlation strength of around 6–7 is suitable. The correlation strength parameter in IMCORR is the sum of the magnitude of the peak height of the correlation map above the mean background and the difference between the largest and second largest distinct peaks, (both normalized by the standard deviation of background values), and an additional contribution proportional to the number of other distinct ‘large’ peaks in the correlation map.

Fig. 2 Number of valid IMCORR correlation matches over the PIIS as a function of correlation strength threshold, for a pre-SLC-off image pair (13 January 2001 and 16 November 2002) as a reference case (dots) and with artificially imposed Landsat7 SLC-off data-gap masks. Processing of the masked images was done with the standard IMCORR (squares) and with the data gaps filled with random data (crosses).

Using the unmodified IMCORR code, the significant decrease in the number of valid correlation matches for SLC-off masked images is obvious, whereas the modified code operating on the masked images performs as well as could be expected given the smaller area of valid data available. There are ∼7.7% fewer valid correlations using the masked images, which have a 22.4% reduction in valid data. The 50% overlap of reference windows in both spatial directions accounts for the smaller reduction in the number of valid correlations. Furthermore, the comparison of the spatial distribution and coherence of good correlations is even more telling (see Fig. 3). The unmodified IMCORR code produces virtually no valid matches in the mask-affected area (Fig. 3c) while the modifications result in very little data loss in the mask-affected areas (Fig. 3b).

The modified IMCORR produces very little bias when comparing the displacements calculated for the unmasked and masked image pairs. In particular the average difference in the easting and northing components of displacement between the unmasked and masked images decreases with increasing correlation strength (as expected). For correlation strengths over 6 (which culls most of the spurious and spatially non-coherent correlation matches) the average difference is less than a pixel and is below 1/4 of a pixel for correlation strengths of 9 and above. In contrast, the unaltered IMCORR analysis shows overall biases that are consistently higher (for each correlation strength threshold), typically by a factor of 2–3. This occurs even though there are very few matches in the mask-affected area (see Fig. 3c), with a larger proportion of results accordingly coming from the unmasked central area of the image where displacements are identical.

Fig. 3. Spatial coverage of the direction of measured displacements (angles measured anticlockwise relative to positive easting direction) for: (a) unmasked test image pair, (b) masked image pair with suitable random number filling and (c) masked image pair with original IMCORR code. The background image is the first ETM+ band 8 image of the pair (13 January 2001). The 1999 grounding line from Rignot and others (2011a) is shown in white, and the 2009 grounding line, after Joughin and others (2010), is shown in black

The distribution of the differences in both x and y displacements between the reference image and the masked image is highly leptokurtic (with an excess kurtosis of >175 for IMCORR correlation strengths above a threshold of 6), indicating that most of the distribution variance is due to relatively few extreme deviations. Accordingly, we assess the performance of our modifications to IMCORR using the median and interquartile range (i.e. the range that contains 50% of the data). For thresholds above 6, the interquartile ranges for x and y displacements are 0.06 and 0.08 pixels respectively. Median differences are zero, reflecting the presence of the regions of the images that are completely unaffected by the SLC-off data gaps.

For the flow direction data shown in Figure 3a and b, the mean difference in directions decreases from 0.23° for a correlation strength threshold of 6, to 0.08° for a threshold of 13.5. The median difference is zero, and the interquartile range is 0.01° for threshold values greater than 6.

These tests indicate that the simple modification works satisfactorily, extending the usefulness of IMCORR into the region where SLC-off data gaps appear.

4. Pine island glacier and ice shelf: dynamics and change

Velocity time series

To examine the acceleration of ice flow in the PIIS discussed earlier, the modified IMCORR correlation code was applied to a sequence of Landsat7 images of the PIIS spanning the period 2001–11, with and without the SLC-off data (see Table 1 for details).

Table 1. Dates of the Landsat7 images of Pine Island Glacier used in this study

During this period there was a continued retreat of the grounding line of PIG, which we indicate in our figures by plotting the 1999 grounding line, determined by interferometric synthetic aperture radar (InSAR) from Rignot and others (2011a), and the 2009 position determined from Terra-SAR-X by Joughin and others (2010). While the location where the PIG first begins to float retreats ∼20 km over this period, Joughin and others (2010) show a central isolated ‘lightly grounded’ region persisting at, or even advanced from, the 1999 position.

Rather than making subjective choices of correlation strength thresholds, we used an IMCORR results flag to select good displacement matches. To reduce variability in the calculated displacements, the results were binned into 1 × 1 km bins based on least trimmed squares (Rousseeuw and Hubert, 1997) requiring at least three data points in each bin. The least trimmed squares method is very robust to outliers. To further reduce noise, the resulting dataset has nearest-neighbour culling applied, so that a data point was only accepted if four of the nine data points in a moving 3 × 3 window had both displacement components within 150 m. Finally the displacements were converted to velocities by subtracting the displacement offsets for nearby stationary points and dividing the result by the time difference between the images.

We applied spatially uniform displacement offsets for collocation for each image pair, based on a cluster of displacement measurements associated with a group of prominent features (nunataks and ice rises) on the north side of the PIIS core flow. For the 2001–02 image pair used in our test case, the uniform displacements were 0.036 pixels in easting and −0.97 pixels in northing, calculated by trimmed least squares. The standard deviations of the displacements used for this example were 1.8 and 1.3 pixels in easting and northing respectively. The corresponding kurtosis values for these displacement distributions were greater than 5 and 8 respectively, showing that to a large extent these standard deviations are due to a small number of outliers among the selected reference points. Similar offsets were generated using these features for each of the image pairs.

The velocity distribution from the 24 January 2005 and 14 January 2007 image pair is shown in Figure 4. On the more northerly side of the ice shelf, flow is nearly constant along streamlines, with speeds increasing towards the centre of the flow, while the southerly portion of our velocity map shows a region of steadily increasing speed in the half of the shelf nearer the calving front.

Fig. 4. Velocity magnitudes for the PIIS (colours) and representative direction vectors, derived from 24 January 2005 and 14 January 2007 image pair. Also shown are the −1595000 easting (vertical white line) and the −304000 northing (horizontal white line) lines for the velocity profiles in Figure 5. The background image is the first ETM+ band 8 image of the pair (24 January 2005).The 1999 grounding line from Rignot and others (2011a) is shown in white, and the 2009 grounding line, after Joughin and others (2010), is shown in black.

The time series of the magnitude of velocities along the −1595000 m easting and −304000 m northing profiles (locations marked in Fig. 4) are shown in Figure 5. In agreement with earlier work (e.g. Joughin and others, 2003; Rignot, 2008; Lee and others, 2012) we see an acceleration in ice flow in the PIIS and near the grounding line of PIG in Figure 5. The magnitudes of ice velocities from the MEaSURES compilation (Rignot and others, 2011b), interpolated onto our two profiles, are also shown for comparison.

Fig. 5. Pine Island Ice Shelf speed history, showing the acceleration of ice flow along (a) −1595000 easting profile and (b) −304000 northing profile, over the five epochs of feature-tracked velocities of the present study. The corresponding ice flow speeds from the MEaSURES compilation (Rignot and others, 2011b) are shown for comparison. The intersections of the profiles are indicated by vertical dashed lines. The locations of the profiles are shown in Figure 4.

Lee and others (2012) have applied the spatial domain correlation method for processing SLC-off images of Ahn and Howat (2011) to a sequence of eight pairs of Landsat7 ETM+ images of PIG/PIIS spanning the period 2003–10, with much smaller time intervals (16–34 days) separating the image pairs, compared to our selections. While the results presented in figure 4 of Lee and others (2012) are for ice velocities along an unspecified profile (presumably an approximate flowline), they show essentially the same behaviour as portrayed in Figure 5, of a main ice-shelf core velocity of ∼3000 m a⁻¹ before mid-2003, increasing to just over 3500 m a⁻¹ by the end of 2006, increasing further to ∼4000 m a⁻¹ by 2008–09, with a slight additional increase by January 2010. Our latest velocities, using images from 2009 and 2011 (centred on June 2010), show the maximum speed at 4200 m a⁻¹ (Figs 5 and 6b), in good agreement with Lee and others (2012).

Fig. 6. Pine Island Ice Shelf ice speed and strain-rate components from feature tracking. Left-hand column from 13 January 2001 to 16 November 2002 and right-hand column 3 November 2009 to 10 February 2011. (a, b) Speed (note the different velocity scales); (c, d) strain rate along the local flow direction; (e, f) strain rate normal to the local flow direction; and (g, h) magnitude of the shear strain rate in the same local reference frame. The convention for the normal strain rates, along and perpendicular to flow, is that positive values correspond to extending flow (and tensile deviatoric stresses). The background images for the left- and right-hand columns are the ETM+ band 8 images for 13 January 2001 and 3 November 2009 respectively. The 1999 grounding line from Rignot and others (2011a) is shown in white, and the 2009 grounding line, after Joughin and others (2010), is shown in black.

In Figure 6a and b we show the speed of PIG and the PIIS in 2001 and 2010 respectively. The increase in speed over the decade is so marked that a shift in the colour scale by 1100 m a⁻¹ is required. The range selected to show the structure in the velocity field remains fixed, and the comparison highlights the extent to which the spatial pattern of speeds remains the same. In both cases, the fastest flow occurs in the southern half of the stream.

The epochs of the velocities in the MEaSURES compilation are not completely clear from Rignot and others (2011b) but may come largely from speckle-tracking on Japanese Advanced Land Observing Satellite (ALOS) Phased Array-type L-band SAR (PALSAR) imagery from 2006 or 2007. PIIS velocities for 2006 using PALSAR in Rignot (2008) are broadly consistent with our 2005–07 velocities. Accordingly, in Figure 7 we present a comparison of our 2005–07 velocities (centred on 2006) with the MEaSURES dataset. We show the percentage differences between our velocities for 2006 and the MEaSURES velocities, presenting the differences in speed, and in velocity components resolved into the easting and northing directions. The differences in speed in Figure 7a indicate agreement is typically better than 5%, except for a zone near the northern end of the ice calving front. There is a second well-defined boundary between positive and negative differences running down the shelf, with the largest excess for our speeds lying along the southern boundary of our velocities. Component differences are also typically a few per cent. The speed difference is understandably dominated by the northing component of velocity, as shown in Figure 7c, given the general direction of flow in the PIIS. A possible explanation for these systematic differences is that velocities from different sources or epochs may have been used in assembling the comprehensive MEaSURES coverage.

Fig. 7. Differences between ice velocities from the present study for the 24 January 2005 to 14 January 2007 image pair and MEaSURES compilation (Rignot and others, 2011b), as a percentage of the MEaSURES values. Comparison of (a) magnitudes, (b) easting components and (c) northing components. The background image is the ETM+ band 8 image for 24 January 2005, where the SLC-off data loss is evident. The 1999 grounding line from Rignot and others (2011 a) is shown in white, and the 2009 grounding line, after Joughin and others (2010), is shown in black.

Dynamics

The strain-rate components from the 2001–02 and 2009–11 image pairs are shown in Figure 6c–h. These are presented in local reference frames aligned with the ice flow directions. The magnitudes of the shear strain rates (Fig. 6g and h) show a very fast-flowing central core in the PIIS from the PIG inflow, which is clearly affected by lateral drag. Since we do not capture velocities for the whole of the ice shelf it is difficult to comment on the width of these shear margins, and coverage is uneven, but in 2010 (Fig. 6h) the transverse dimensions of the central fast flow clearly increase seawards of the isolated region shown as grounded by Joughin and others (2010). The shear rates along the southern margin of our velocities mostly exceed 0.02 a⁻¹, reflecting the contrast between the PIG inflow and the rest of the ice shelf.

The underlying satellite image shows bands of heavy surface crevassing, indicative of possible shearing flow, flanking the regions where we have retrieved ice displacements, particularly near the grounding line. It is likely that features in these regions change too rapidly to be tracked over intervals of a year or more, and that image pairs spanning shorter intervals could have more success capturing the motion here, although at some cost in accuracy. Another indication of the restraining forces at work in Pine Island Bay is given by the longitudinal strain-rate components (the spatial rate of speed change along the local flow direction) shown in Figure 6c and d. These do not show the large extending deformation (positive strain rate) that would be expected of an unconfined ice shelf, but even display zones of along-flow compression, particularly on the northern side extending for ∼20 km downstream of the grounding line in 2001–02 (Fig. 6c). The more southerly portion of the 2009–11 flow shows a small, steady, positive longitudinal strain rate (Fig. 6d) once clear of the grounding zone, as expected from the developing flowband of increasingly fast-flowing ice seen in Figure 6b. Bindschadler and others (2011) comment that the longitudinal strain rate along one of their 2004–05 profiles from the grounding line to the calving front averaged 0.005 a⁻¹, largely consistent with these measurements. The colour scales for the longitudinal and transverse strain rates are asymmetric, to concentrate on the preponderance of large extending (tensile) deformations (positive strain rates), particularly for the transverse case.

The transverse strain rates (Fig. 6e and f) show large positive values in the first 25–30 km of the PIG stream. In the 2001–02 case this diverging flow occurs over the grounded ice as well, perhaps due to the extending central projection of the grounding zone. In the 2009–11 case, the diverging flow seems more closely associated with crossing the initial grounding line, although this may be a matter of the more complete coverage. It is not clear whether there is a significant influence by the isolated grounded region mapped by Joughin and others (2010). The region of higher transverse strain rates extends a few kilometres further downstream in 2009–11 (Fig. 6f). In both cases there is a transition, aligned essentially orthogonal to the flow, to much lower transverse strain rates and then to regions, well within the PIIS where the transverse strain rates become negative, indicating converging flow, particularly near the margins of our measurements. This may reflect the merging of additional ice inflow in the PIIS on the southern margin of our measurements, and the proximity of the grounding line on the northern side.

In the final half of the flow across the PIIS there are also some indications of possible bands of alternating compressive (converging) and extensional (diverging) motion parallel to the flow (particularly in the 2001 case of Figure 6e) at the transverse scale of a few kilometres. Vaughan and others (2012) pointed out that there are numerous sub-glacial channels in the PIIS (many of them parallel to the flow), which they attribute to basal melting processes. They have proposed that the resulting disturbance of local hydrostatic equilibrium generates stresses that would locally distort the flow while creating arrays of surface crevasses and a basal crevasse where the ice surface and base respectively are in tension. It is an interesting suggestion of possible finer-scale structure in ice-shelf flow dynamics. With the transverse strain rates in this region varying mainly within the range ±0.005 a⁻¹, which corresponds to our estimate of uncertainty for strain rates, further examination of the underlying displacement fields and investigation of correlation with surface features, ice thickness, etc., would be required to decide whether this is physically significant.

Mass budget and basal melting

A gridded dataset of recent ice thicknesses of the PIG and PIIS has been derived (Gogineni, 2012) from radio-echo sounding (RES) data gathered as part of NASA’s Operation IceBridge in 2009 and 2010. In combination with ice velocities, this allows the calculation of ice fluxes for the PIIS, and provides the opportunity to disentangle the roles of ice flow dynamics (i.e. strain thinning) and basal melting in determining ice-shelf thickness and ice-shelf thickness change.

Mass conservation relates the divergence of the horizontal ice flux to the local rate of ice thickness change plus basal melting less surface accumulation, and this extends by the divergence theorem to total ice flux through the perimeter of a region and area averages of the other terms.

We constructed ice fluxes by multiplying our 2009–11 velocity data with the L3 gridded RES ice thickness data (Gogineni, 2012), interpolated onto the same 1 × 1 km grid as the velocities. These are nearly simultaneous datasets, which is important since this portion of the PIIS is reportedly decreasing in thickness at 6.8 m a⁻¹ (supplementary table 1 of Pritchard and others, 2012). The resulting ice flux distribution is shown in Figure 8a, together with streamlines constructed from the ice velocity data, and a set of ‘gates’ a–d that each span the same set of streamlines.

Fig. 8. (a) Pine Island Ice Shelf annual average ice flux density (colour shading) and streamlines (black lines) from the numerical integration of the velocity field for the period 2009–11. (b) Along-flow derivative of the flux density in (a). Also shown are the streamlines and ice flux gates a–d. The background image is the ETM+ band 8 image for 3 November 2009. The 1999 grounding line from Rignot and others (2011 a) is shown in white, and the 2009 grounding line, after Joughin and others (2010), is shown in black.

From consideration of the crossover analysis of individual RES flight lines by Gogineni (2012) and the actual ice thicknesses we estimate that the uncertainty propagated into ice fluxes is <10%. The ice thickness data are not corrected for the density of a firn layer, but since the dominant mass loss is basal melting, this bias in total ice fluxes does not seriously affect ice loss calculations for the PIIS, provided the firn layer is reasonably uniform.

The extent of loss of ice is evident, with the majority of the net loss occurring near the grounding line. To demonstrate this, Figure 8b shows the along-flow derivative of the flux field of Figure 8a. Specifically, the integrated ice fluxes through flux gates a–d in Figure 8a are 75, 52, 41 and 36 km³ a⁻¹ respectively. These fluxes come from integrating the ice flux component normal to the gate, and involve some minor linear interpolations to fill missing values for gates c and d. The ice flux gates were selected to cover the region where a stream of ice flow could be tracked, delineated by two streamlines of our velocity and flux fields. Regrettably this does not span the full inflow of PIG. The width of the selected band was a compromise between capturing as much breadth of flow as possible and starting with gate a as far upstream as possible. With the retreat of the central part of the PIG grounding line >20 km since the 1990s (Joughin and others, 2010), our gate a lies quite close to the 2009 grounding line.

The loss of 39 km³ a⁻¹ between gates a and d, occurring over a region of 1005 km², corresponds to an area-averaged loss rate of 39 m a⁻¹. This is almost 60% greater than the corresponding overall loss rate for the entire PIG stream in 1996 of 24 ± 4 m a⁻¹ calculated by Rignot (1998). To put our 39 m a⁻¹ rate of ice loss in context, compare it with the estimated −6.8 m a⁻¹ recent rate of change of ice-shelf thickness for the PIIS (Pritchard and others, 2012) and a contribution from surface accumulation in the PIIS region of 0.7–0.9 m a⁻¹ (Lenaerts and others, 2012). This highlights that significant decreases in ice-shelf thickness may involve relatively small changes in the delivery of ocean heat. In the situation where the ice-shelf thickness is decreasing even while accumulating ice at the surface, the basal melt rate must exceed the loss estimated from the ice fluxes to close the mass budget, so our estimate of basal melt rate averaged over the band between gates a and d is 47 m a⁻¹. This value is higher than the estimates of Bindschadler and others (2011) (19 m a⁻¹, for 2005–06) from four profiles across the PIIS, and the recent integrated oceanographic estimate of Jacobs and others (2011) (33 m a⁻¹, for 2009).

The decrease in integrated ice flux between gates a and b is 23 km³ a⁻¹ over an area of 305 km², corresponding to a local average loss rate of 76 m a⁻¹ (a few per cent higher if one excludes the relevant portion of the 30 km² of the suggested lightly grounded remnant of the ‘ice plain’). Allowing for ice-shelf thickness decreasing at around 7 m a⁻¹ and surface accumulation (<1 m a⁻¹) yields an estimated basal melt rate of 84 m a⁻¹ for this zone. Our first two gates are separated by ∼17 km, so assuming that gate a is close to the 2010 grounding line, it is of interest to compare this localized melt rate with previous studies. Our melt rate is in reasonable agreement with estimates mentioned earlier of 78.2 and 92.7 m a⁻¹ (for similar proximity to the grounding line) from two of the profiles studied by Bindschadler and others (2011) for the 2005–06 epoch. Rignot (1998) estimated an average melt rate from ice flux losses over the first 20 km from the PIG grounding line of 50 ± 10 m a⁻¹ for 1996. Our corresponding estimate of 76 m a⁻¹ is ∼50% higher. Of course, the relevant region has itself also migrated ∼20 km landwards. From Figure 8 the region of maximal ice loss appears even more localized.

There is also general agreement between our relatively localized melt rate increases and the ∼50% increase in basal meltwater production between 1996 and 2009 reported by Jacobs and others (2011).

We estimated the ice flux through gate a in 2010 at 75 km³ a⁻¹ which is less than the usual estimates of the total ice flow from PIG into the PIIS (see, e.g., Rignot, 2008), reflecting ice flow not captured by our gate a. We conclude this section with some more speculative estimates of recent PIIS melting. First we estimate the total input at the PIG grounding line in 2010 following the approach of Rignot (2008) and rescaling the 2006 flux estimate using the increase in ice velocities. Velocity increases between 2006 and 2010 appear quite spatially uniform (see Fig. 5), so we estimated the total flux from PIG in 2010 to be 125 ± 10 km³ a⁻¹, a factor of 1.17 higher than the 2006 value from Rignot (2008). For simplicity we used the ratio of peak speeds on the profiles in Figure 5b. We deducted 3 km³ a⁻¹ to reflect 4 years of decreasing grounded ice thickness upstream of the PIG grounding line at 6 m a⁻¹ (Pritchard and others, 2009).

We have also estimated the ice flux out of the PIIS in 2010 at 49 ± 5 km³ a⁻¹ by integrating our ice fluxes across the entire span of valid flux values in Figure 8, as close to the calving front as practicable. Combined with the scaled estimate of input from PIG this provides an approximate estimate for the imbalance of ice transport fluxes into and out of the PIIS as 76 ± 11 km³ a⁻¹ involving the fast-moving PIG stream. This is ∼50% more than the corresponding loss measured by Rignot (1998) for 1996, the similarity to the increase in melt rates indicating that our scaling argument is not too far astray. The mass flux imbalance computed using our flux gates (39 km³ a⁻¹) apparently captures approximately half the total ice flux imbalance for the flow from PIG through the PIIS.

Information about the other inflows to this part of the PIIS seems scarce. A recent modelling sensitivity study (Larour and others, 2012) suggests two significant glaciers contributing ∼9 km³ a⁻¹ (using 1996 ice velocities).

If we assume the recent 6.8 m a⁻¹ decreasing ice-shelf thickness (Pritchard and others, 2012) occurs uniformly over an ice-shelf area of ∼3000 km² (Schodlok and others, 2012) the additional basal melting would be 20 km³ a⁻¹. Aggregating these contributions we can estimate net basal melting from the PIIS in 2010 at 105 km³ a⁻¹, though with considerable uncertainty, say ±20 km³ a⁻¹.

These admittedly speculative glaciological estimates of the contributions to recent basal melting are in reasonable agreement with the 2009 oceanographic estimate for the PIIS of 85 ± 5 km³ a⁻¹ (Jacobs and others, 2011). On firmer ground are the indications of change, with ∼50% increases in basal melting seen in both the glaciological and oceanographic estimates between 1996 and 2010.

5. Conclusions

We have shown that a slight modification to the highly efficient and widely used IMCORR image correlation software (Scambos and others, 1992) permits automated displacement detection on images with data gaps, particularly relevant to Landsat7 images acquired after the SLC hardware failure in 2003. While the focus here has been on the application to Landsat7 SLC-off images, the method is general and applies to any image with identifiable data gaps, whether inherent in the image or intentionally applied to mask out unwanted correlations on irrelevant features such as clouds.

We have demonstrated the utility and skill of the method by testing it with artificial masking of Landsat ETM+ images, and then constructing a time history of the velocities in the Pine Island Ice Shelf for a sequence of pre- and post-SLC-off Landsat7 images, and showing the progressive acceleration of the flow in the floating portion of the Pine Island Glacier system. The resulting velocity datasets are freely available from the Australian Antarctic Data Centre (http://data.aad.gov.au/).

We made some comparisons with the MEaSURES velocity dataset of Rignot and others (2011b) and commented on the dynamics revealed by the flow. The freely available gridded RES ice thickness data (Gogineni, 2012) from NASA’s Operation IceBridge have enabled us to look briefly at the recent pattern of net ice loss from the flow in the ice shelf, and interpret it in terms of high rates of basal melt. We infer melt rates near the PIG grounding line in 2010 that are in reasonable agreement with estimates from Bindschadler and others (2011) for 2005–06. These represent a ∼50% increase over the corresponding melt rates in 1996 (Rignot, 1998). This is consistent with similar increases in total meltwater production at the PIIS over the same interval by Jacobs and others (2011).