Pine Island glacier ice shelf melt distributed at kilometre scales

Abstract. By thinning and accelerating, West Antarctic ice streams are contributing about 10% of the observed global sea level rise. Much of this ice loss is from Pine Island Glacier, which has thinned since at least 1992, driven by changes in ocean heat transport beneath its ice shelf and retreat of the grounding line. Details of the processes driving this change, however, remain largely elusive, hampering our ability to predict the future behaviour of this and similar systems. Here, a Lagrangian methodology is developed to measure oceanic melting of such rapidly advecting ice. High-resolution satellite and airborne observations of ice surface velocity and elevation are used to quantify patterns of basal melt under the Pine Island Glacier ice shelf and the associated adjustments to ice flow. At the broad scale, melt rates of up to 100 m yr−1 occur near the grounding line, reducing to 30 m yr−1 just 20 km downstream. Between 2008 and 2011, basal melting was largely compensated by ice advection, allowing us to estimate an average loss of ice to the ocean of 87 km3 yr−1, in close agreement with 2009 oceanographically constrained estimates. At smaller scales, a network of basal channels typically 500 m to 3 km wide is sculpted by concentrated melt, with kilometre-scale anomalies reaching 50% of the broad-scale basal melt. Basal melting enlarges the channels close to the grounding line, but farther downstream melting tends to diminish them. Kilometre-scale variations in melt are a key component of the complex ice–ocean interaction beneath the ice shelf, implying that greater understanding of their effect, or very high resolution models, are required to predict the sea-level contribution of the region.


Introduction
Thinning ice shelves (Pritchard et al., 2012;Shepherd et al., 2010) and the corresponding decrease in the restraint experienced by inland ice flow (Flament and Rémy, 2012;Joughin et al., 2010;Payne et al., 2004;Pritchard et al., 2009;Zwally and Giovinetto, 2011) are recognized as major drivers of current Antarctic ice loss. This ice-shelf change is particularly pronounced in West Antarctica, where in recent decades the grounded part of Pine Island glacier (PIG) has thinned, accelerated and retreated (Rignot, 2008;Shepherd et al., 2001), in response to increased oceanic heat transport beneath its floating ice shelf and resulting feedbacks (Jacobs et al., 2011). Some efforts have been made to relate basal melt to ocean temperature  and the broad-scale iceshelf geometry (Little et al., 2009), but the detailed patterns and rates of basal melt on specific ice shelves are known only on relatively coarse scales (Payne et al., 2007). Without a thorough understanding of the processes that control the dominant scales of ice-shelf melt, future projections of changes in PIG and similar glaciers will be dependent on melting parameterizations that are poorly constrained by observations (Joughin et al., 2010;Katz and Worster, 2010).
One difficulty in studying PIG arises from the unsteady nature of the glacier. The ice geometry is evolving rapidly (years) and the ice is being advected at speeds above 4 km yr −1 whilst the ocean beneath is expected to respond at sub-annual scales to both local and remote forcing (Steig et al., 2012;Thoma et al., 2008). In addition, there are many reasons to expect that the spatial pattern of melt is complex. The floating tongue of PIG contains a series of both longitudinal and transverse channels (Bindschadler et al., 2011) (Figs. 1a, b and 2a), and on a similar ice shelf (Rignot and Steffen, 2008) basal melt was found to be strongly  (Joughin et al., 2010), while the green indicates the grounding line of the fast flowing trunk in 1996 (Rignot, 2008). The black box delineates the area covered by the airborne radar survey in (b). The black vector indicates the ice flow direction. Inset shows the location of Pine Island Glacier in West Antarctica. (b) Surface elevation from the January 2011 radar survey linearly interpolated on the SPIRIT grid. Color bar applies to concentrated along longitudinal subglacial channels (elongated in the ice flow direction). Beneath the PIG shelf basal crevasses are located above each channel apex. The formation of such crevasses may be a response to basal melting, suggesting that changes in channel-scale ice-ocean dynamics could indirectly affect the structural integrity of such ice shelves . Conversely, a recent study using a coupled ice-ocean model suggested that the presence of melt channels allows ice shelves to survive higher sub-ice ocean temperatures than they would otherwise (Gladish et al., 2012). An important step towards improving understanding of these processes, in the crucial West Antarctic setting, is the development of the ability to measure the spatial pattern of melt rate beneath PIG.

Methods
This study uses ice surface velocity observations and a set of high-resolution digital elevation models (DEMs) of the floating tongue of PIG surface acquired in 2008 and 2011 to quantify basal melt at short (< 10 km) to large (> 10 km) lengthscales. Analysis of the melt distribution provides some insight into the most influential glaciological and ocean processes at play. The analysis relies on the co-registration of the DEMs, by cross correlation, to allow a direct comparison of the same parcels of ice as they travel down the ice shelf (i.e. a Lagrangian comparison of the surface elevation). This method does not, as has previously been required, assume steady ice flow. Given that previous studies have shown that the surface features on PIG are, basically, reflections of the basal topography brought about by relaxation towards the hydrostatic condition , we expect that most of the surface elevation changes are due to basal melt, and we employ ice surface velocities to quantify the elevation change due to ice divergence.

Input data products
We used a pair of optical stereo-imagery, 40 m resolution DEMs produced from images acquired on 5 January (Fig. 1a) and 13 March 2008 as part of the SPIRIT project (Korona et al., 2009). Their location accuracy is better than 40 m, as verified using control points in the nearby Hudson Mountains. Obviously erroneous areas (due to high slopes or very weak surface roughness) were identified and removed by setting a maximum 20 m difference threshold between the two DEMs. Next, both DEMs were vertically adjusted by a constant offset using the mean difference with ICESat profiles acquired on 10 and 13 March 2008, which were previously corrected for tidal and accumulation effects (Pritchard et al., 2012). Horizontal variation in the vertical bias is neglected. Mean differences with ICESat are 2 ± 5.7 m and 2.9 ± 5 m (see Appendix A), confirming SPIRIT accuracy of ±6 m (Korona et al., 2009).
About three years later, in January 2011, a detailed grid of ice-penetrating radar was acquired during two airborne surveys over a section of PIG's floating tongue . The 500 m separation of transverse lines across the glacier provides a clear snapshot of the ice structure at the time (Figs. 1b and 4a). These sections were fixed using GPS with a horizontal accuracy of ±0.5 m and vertical accuracy of ±2.4 m based on cross-overs. For the purpose of comparing the surface elevation products, the airborne elevation was linearly interpolated onto the SPIRIT grid.
We also used an estimate of the ice surface velocity, from November 2008, which was derived using Advanced Land Observation System (ALOS) Phased-array L-band Synthetic-Aperture Radar (PALSAR) observations (Joughin et al., 2010). This velocity field covers most of the ice shelf with a resolution of 500 m and an error of about 10 m yr −1 , and little change in the general velocity has been observed since, indicating a hiatus in the glacier acceleration (Joughin et al., 2010).

Basic theory
Assuming (a) horizontally uniform ice/firn density and (b) vertically constant ice velocity, the depth averaged ice thickness evolution deriving from mass conservation is given by where H is the ice thickness, U is the ice velocity, M is the thickness change resulting from melt or accumulation, either surface or basal, and ∇ is the del operator. Assuming also (c) hydrostatic equilibrium and (d) negligible surface accumulation and melt, H , U and M can be replaced by surface elevation h, surface velocity u and basal-melt-induced surfaceelevation change m s in Eq. (1), respectively, giving: Assumption (c) is verified at spatial scales greater than a few ice thickness (here a few km, see Sect. 3). Although typical, (a) (also discussed in Sect. 3) and (b) are possibly less accurate, and finally, d is consistent with the few visual observations of the ice surface that suggest there is rather little snow on the fast flowing trunk of PIG's floating tongue (Bindschadler et al., 2011) and sufficiently cold surface air temperature to prevent significant surface melt (Steig et al., 2009;Trusel et al., 2012).
Basal melt can then be deduced from m s using where ρ i = 918 kg m −3 and ρ sw = 1029 kg m −3 are, respectively, ice and seawater densities.

Alternative Eulerian framework
Previous studies of this and other ice shelves have employed an Eulerian framework, where the Lagrangian elevation change (the first term in Eq. 2, Dh/Dt), is decomposed into two terms: the Eulerian elevation change (dh/dt) and the surface height advection (u · ∇h), following The Eulerian elevation change (dh/dt) is measured by satellite altimeters and combined with precipitation estimates to infer ice shelf thinning rates (Pritchard et al., 2012;Shepherd et al., 2001Shepherd et al., , 2010Wingham et al., 2009). Alternatively, the height advection (u · ∇h) and surface divergence (h∇ · u) are calculated together as the ice flux divergence, which gives the melt rate if the ice is assumed to be in steady state (e.g. Rignot and Steffen, 2008). If changes in flux divergence are neglected, thinning rates (deduced from dh/dt) can then be associated with changes in melting (e.g. . In the case of steady ice geometry, or over a very short time period, the Eulerian and Lagrangian framework are likely to give similar results. But for Pine Island, where channels are being formed and deformed over the course of the observational period (3 yr here), sparsely sampled Eulerian estimates are likely to suffer from aliasing. To illustrate the benefits of the Lagrangian methodology presented here, both the Eulerian elevation change and the height advection will be computed.

Lagrangian elevation-change-based melting calculation
A process of normalised two-dimensional cross correlation between SPIRIT and 5-by-5 km mosaic boxes from the airborne radar survey is used to match the DEMs and account for horizontal advection of channelised ice patterns. Basically, each 5 km × 5 km parcel of ice from one DEM is compared with all possible parcels of the second DEM, and a match is found where the correlation between the two is maximum. The Lagrangian elevation change (Dh/Dt) is computed by difference, with an error of ±6.5 m over 3 yr, or ±2.2 m yr −1 . Horizontal displacement or ice velocity u averaged between 2008 and 2011 is retrieved as a by product.
Equations (2) and (3) allow for the derivation of basal melt in the unsteady framework of the PIG's floating tongue at an effective resolution of around 1 km (Fig. 1d). The calculation was made by comparing the airborne radar survey DEM with the two SPIRIT DEMs, and the results shown are averaged between the two estimates. A similar, direct comparison of the two SPIRIT DEMs provided an estimate of the error associated with the SPIRIT products. We found that the latter can lead to 2.5 m yr −1 error in Lagrangian elevation change, but the spatial patterns of the error are not significantly correlated with the patterns presented here, and so do not modify our conclusions.
The ALOS PALSAR velocity product provides a means to evaluate the velocity field we obtained from the coregistration of the DEMs, and indeed the two velocity fields agree within their respective margin of uncertainties, leaving no doubt that the DEM matching exercise was successful. The second term of Eq. (2), the elevation change due to icevelocity divergence, can be interchangeably calculated using this and the previously described velocity field.

Large-scale budget
Finally, we will find that melting can be decomposed in two dominant scales: large scale (> 10 km, labelled with a bar) and small scale (< 10 km, labelled with a prime), such that Small-scale anomalies were computed by subtracting the moving 16 km 2 square window smoothed field (large-scale field) from the original field. By definition, so that substituting Eqs. (5) and (4) in Eq.
(2) and applying the large-scale smoothing to the resultant gives dh/dt+ū· ∇h +h (∇ ·ū)+u · (∇h )+h ∇ · u =m s . (7) The first three terms on the left side of Eq. (7) are, respectively, the large-scale Eulerian elevation change, the largescale surface height advection, and the large-scale flow divergence. The following two terms represent the averaged impact of small-scale advection and small-scale flow divergence on the large-scale surface expression of melting (right side of Eq. 7).

Overall melting patterns
The DEMs we have used fully capture the complex surface topography of PIG's floating tongue ( Fig. 1a and b), and the ice flow divergence at high resolution ( Fig. 1c) allows us to derive the pattern of sub-shelf melt (Fig. 1d) from our Lagrangian methodology. At large scales (> 10 km), rapid surface lowering (Fig. 1a) in the area downstream from the grounding line and in the shear margins of the fast flowing trunk indicates a relatively high melt rate. Although it does not cover the entire grounding line area, our map ( Fig. 1d) shows basal melt rates of 100 m yr −1 near the grounding line, generally reducing to 30 m yr −1 around 20 km downstream of an isolated area of grounding identified in 2009 (Joughin et al., 2010). Peak melt rates are expected in these areas where the ice comes into contact with the warmest waters . This pattern of melt roughly agrees with the results of a 2-dimensional plume model that used a smoothed basal ice geometry (Payne et al., 2007), which also suggested an area-averaged melt rate for the entire ice shelf of ∼ 27 m yr −1 . Assuming similar spatial variability and extrapolating into areas that remain unsampled here, i.e. inferring stronger basal melt near the grounding line and weaker melt close to the ice front, then our results are consistent with averaged melt values of 33 m yr −1 deduced from 2009 oceanographic constraints (Jacobs et al., 2011). Ice flow divergence ( Fig. 1c) accounts for less than 20 % of the ice-shelf thinning near the grounding line.

Large-scale balance
A more detailed view of the surface elevation balance at large scale can be obtained by using Eq. (7), computing each term individually in the Eulerian framework. Figure 2 shows the main terms of the balance, in which the averaged impact of small-scale advection and small-scale divergence on the large-scale surface expression of melting are negligible. A clear result is that between 2008 and 2011, during the period when a hiatus in the glacier acceleration has been reported (Joughin et al., 2010), the system appears to be roughly in steady state at large scale, and an apparent balance is found between the large-scale height advection (Fig. 2b) and the large-scale surface expression of melting (Fig. 2d). In other words, to a first approximation large-scale basal melting seaward of the grounding line imparts a form to the ice base as it passes by.
Assuming this balance is exact, and using the broader coverage of the large-scale height advection term, basal melting can be estimated over the entire fast flowing trunk of  the ice shelf (using Eq. 3). Basal melting is thus found to reach over 200 m yr −1 near the grounding line but averages to 37 m yr −1 over the sampled area ( Fig. 2) giving an integrated melt-induced steady state mass loss of 87 km 3 yr −1 to the ocean. Although the extreme numbers found very close to the grounding line should be treated with caution, since the ice there may not be freely floating , this estimate is practically identical to the 85 ± 6 km 3 yr −1 determined in 2009 from oceanographic constraints (Jacobs et al., 2011). These results are thus consistent with the conclusion that the ice shelf was nearly in steady state between 2008 and 2011, albeit in association with continued strong thinning upstream of the grounding line. It is worth noting, however, that our error estimate for large-scale elevation change is 2.2 m yr −1 , equivalent to large-scale thickness change of about 20 m yr −1 , well above the typically quoted 5 m yr −1 (Pritchard et al., 2012;Shepherd et al., 2004;Wingham et al., 2009) for Pine Island ice shelf. We therefore cannot draw a firm conclusion about the system steadiness.

Channel-scale melting patterns
We may also invert this quasi-steady state concept by considering only the small-scale, channel-scale anomalies remaining after the large-scale means are removed (Fig. 3). Overlaid on the large-scale ice surface topography is the reflection of a series of 3-4 km-wide, 100-200 m-deep sinuous longitudinal basal channels aligned with the ice flow on the central part of the ice shelf (but see Appendix B) and 1-2 kmwide, 50 m-deep transverse channels on its flanks (Figs. 1a, b and 3a). These features have been discussed in detail elsewhere (Bindschadler et al., 2011;Vaughan et al., 2012). Another salient feature is the confluence of two 3 km-wide channels at the southernmost part of the ice front (Fig. 3a), where concentrated meltwater outflows have been observed (Jacobs et al., 2011).
The scale and nature of this basal channel system and their correlation with basal melt water outflows (Hellmer et al., 1998;Mankoff et al., 2012;Payne et al., 2007)  to speculation about their role in structuring ice-ocean interactions. In particular, based on observations and idealized modelling, the longitudinal  and transverse (Bindschadler et al., 2011) channels on the floating part of PIG have been assumed to be created by basal melt. Furthermore, recent ice-ocean simulations indicate that basal channels can be initiated by irregularities in the ice thickness at the grounding line and subsequently enlarged by basal melt, with ice dynamics playing only a minor role (Gladish et al., 2012).
Our results support the hypothesis of melt channel formation by showing that channel features appear near the grounding line, quickly reach maximum surface expression, and slowly vanish downstream towards the ice front (Fig. 3a). Interpretations based on an implicit assumption of steady state, however, should be treated with caution on this rapidly changing ice shelf. Indeed, the 2009 grounding line position shows a clear retreat from 1996, but with a grounded area 15 km downstream, where longitudinal channel signatures are deepest (Fig. 3a). In 2008, at least, channels upstream of the newly grounded area were weak and clear relationships between evolving grounding line, basal melt and channels are not easy to establish.
Just downstream of the isolated area of grounding identified in 2009, the spatial scales of channels are evident in our calculated basal melt rate (Fig. 3b). The channel-scale melt variability decreases from 40 m yr −1 near the grounding line to about 15 m yr −1 downstream. In other words, largescale melt is modulated at the channel-scale by 40 % near the grounding line, and by close to 50 % downstream. Alternatively, the amplitude of channel-scale melt reaches 80 and 100 % of the total. Moreover, channel-scale melt patterns are strongly anti-correlated (correlated) with the surface (basal) topography near the grounding line, but correlated 10 km downstream in the central part of the ice shelf where longitudinal features dominate, and weakly correlated closer to the ice margins where transverse features appear (Fig. 3c). These results therefore support the interpretation that where the ice goes afloat and the large-scale melt rates are high, basal melt carves the channels. As the channels are advected downstream they reach an area where basal melt is lower and the opposite effect occurs, i.e. melting is larger on keels between the subglacial channels and the subglacial topography is therefore smoothed.

Ice adjustment to melting
While depth-averaged ice deformation is not anticipated to play a significant role in reaction to the carving of basal channels (Casassa and Whillans, 1994;Gladish et al., 2012), vertically sheared deformation within the ice is probably important. Indeed, observations made near the grounding line of Rutford Ice Stream have shown that as the ice within (between) channels sags (rises) toward hydrostatic compensation, the associated bending of the ice shelf sets up transverse  compression (extension) in the near-surface layers and corresponding extension (compression) near the ice base (Jenkins et al., 2006). A recent modelling study showed that the bending stresses caused by channelisation on PIG are sufficient to explain the observations of basal crevasses at the apex of basal channels and surface crevasses above basal keels . Both symptoms are consistently observed here, lending support to this concept. Firstly, the radar data show that the surface features do not fully reflect the basal channels (see below). Secondly, extracting the channelscale anomalies in ice velocity (Fig. 4) shows that the pattern of ice deformation (compression transverse to the surface depressions and extension transverse to keels) is evident in the ice-surface velocity field.

Relevance of channel-scale melting
Given the amplitude of the channel-scale melting signature, it is clear that kilometre-scale processes crucially affect the overall ice shelf-ocean interactions. These spatial scales are also likely to affect observations and methodologies used for investigating ice shelf dynamics. Figure 5 presents all the terms from Eqs. (2) and (4): the Eulerian elevation change (Fig. 5a), the height advection (Fig. 5b), and the ice flow divergence (Fig. 5c) that should all sum to give the surface expression of basal melt (Fig. 5d). Ice divergence is found to be playing a minor role a few kilometres away from the grounding line, where the Lagrangian elevation change (not shown) roughly mirrors the surface expression of basal melting. The Eulerian elevation change is not negligible at the channel scale (Fig. 5a), confirming that channel-ocean interactions are evolving in time as channels are advected with the ice flow. The height advection term presented in Fig. 5b is much larger than any other, violating the expected balance between the four presented terms. This large and noisy pattern derives entirely from the advection of channel features by the largescale ice velocity (not shown). In addition to aliasing issues, we surmise that the computation of the surface elevation gradient increases the level of noise initially contained in the SPOT5 DEM. Regardless, this apparent discrepancy actually reveals a strength of the Lagrangian method presented here, where the sum of the Eulerian elevation and height advection ( Fig. 5a and b) is directly computed as a single quantity and the noise thereby removed by means of methodology.
Another point demonstrated by these diagnostics is the relative importance of the channel advection signal in the observed surface elevation change (Fig. 5a). This signal will contaminate altimeter-elevation-change-based estimates of ice shelf thinning, but the level of noise it introduces largely depends on the instrumental footprint (which can vary from a few hundred metres to a few kilometres), while the sampling size determines whether or not such noise can be reduced to yield meaningful estimates.

Discussion
Although it is clear from previous observations ) that sub-glacial channels grow substantially in magnitude in the region where PIG goes afloat, it is not possible to determine if the location and separation of these channels is controlled by small-scale topography imparted as the ice crosses the grounding line. Another contributing factor at the grounding line is the drainage of subglacial melt (1.7 km 3 yr −1 , Joughin et al., 2009) from beneath the grounded glacier that provides a source of buoyancy, initiating ocean plumes that drive the melt (Jenkins, 2011). It is also possible that part of the observed surface signal results from sudden un-grounding from a rapidly evolving channelized subglacial terrain.
The methodology developed in this study assumes a uniform ice density-depth profile and approximate hydrostatic balance of the ice. Such assumptions are necessary to derive basal melt from surface elevation change. The first assumption could be invalidated locally by the presence of surface crevasses or variation in the amount of air in the firn. The second is expected to be valid at spatial scales greater than a few ice thicknesses (here a few km). Channel scales (1-3 km) are therefore in a range where bridging stresses in the ice are able to compensate for small geometric irregularities arising from shear-processes within the ice or, indeed, basal melt. In an attempt to shed light on such effects we present two diagnostics here.
The apparent air thickness in the firn can be computed by assuming a geoid model and by making the hydrostatic approximation (Holland et al., 2011). This results in an apparent air thickness pattern varying from 0 to 20 m that is very closely related to the channel-scale geometry ( Fig. 6c and  d). In areas of relatively strong basal melt near the grounding line (see Fig. 1d), and along longitudinal channels, we find that apparent air thickness anomalies correlate positively with basal elevation anomaly (Fig. 6d). The firn within surface channels apparently contains about 20 m more air than between channels. Such a large amplitude is difficult to explain as a result of surface snow accumulation and temperature variations associated with the ridge/valley structure of the channels, even before the presence of surface crevasses between the channels  acting to reduce bulk firn density there is taken into account. Conversely, in some areas of transverse channels, the surface between the channels apparently contains more air than within the Figure 6 a) Ice base elevation from airborne radar observations assuming 10 m air thickness in the firn, and interpolated on the SPIRIT grid, overlaid with 0 m contour of the corresponding small-scale surface elevation anomaly. Same contour applies to all panels. b) Departure from hydrostatic equilibrium, or difference between base elevation in a and base elevation deduced from the surface assuming flotation and constant ice density of 918 kg m -3 (i.e. no air). c) Air thickness (Holland et al., 2011) necessary for hydrostatic equilibrium assuming air density of 2 kg m -3 and constant geoid anomaly from WGS84 of -26 m. d) Correlation over moving 4 km 2 square boxes between departure from hydrostatic equilibrium (c) and small-scale surface elevation anomaly (black contours). Note that a very similar result applies for the correlation between air thickness and small-scale surface elevation anomaly.  (Holland et al., 2011) necessary for hydrostatic equilibrium assuming air density of 2 kg m −3 and constant geoid anomaly from WGS84 of −26 m. (d) Correlation over moving 4 km 2 square boxes between departure from hydrostatic equilibrium (c) and small-scale surface elevation anomaly (black contours). Note that a very similar result applies for the correlation between air thickness and small-scale surface elevation anomaly. channels (Fig. 6d). This is consistent with the presence of crevasses between the channels, but the contradictory results for the longitudinal channels leads us to believe that although variations in bulk firn density probably play a minor role in modifying the ice surface/base relationship, they should not modify our conclusions significantly.
To test the assumption of hydrostatic equilibrium, we assume a constant 10 m air thickness in the firn and compute the difference between the predicted draft of free floating ice and the ice shelf draft measured by the airborne radar survey ( Fig. 6a and b). In areas of relatively high basal melt near the grounding line and over longitudinal channels, the difference indicates that basal channels are only partially reflected at the surface ( Fig. 6b and d). This can be explained by the fact that bridging stresses within the ice allow for recently carved basal channels to be not fully represented at the surface, implying a bending signal of relaxation in the ice and the observed transverse surface convergence and divergence in the ice surface velocity field (Fig. 4). This also suggests that where channels are being carved (near the grounding line, see Fig. 3), basal changes are only partially reflected at the surface, and our surface elevation-change-based results for basal melt assuming floatation are contaminated by the signature of ice relaxation. However, the area dominated by lower melt and transverse channels shows the opposite: overrepresentation of the basal geometry at the surface. The latter is consistent with the removal of channels in this region by rapid melting of the keels and slower melting at the crests, with relaxation to free floatation again operating on a longer timescale. We surmise that, at this smaller channel scale, the exact relationship between the ice surface and base results from a complex interplay between non-uniformity of the firn density and bridging stresses. Regardless, our methodology should provide a good approximation to channel-scale basal melt pattern, albeit slightly smoothed by the longer timescale ice relaxation processes. However, sub-channel or sub-kilometre-scale anomalies are unlikely to be meaningful.

Conclusions
Previous work has indicated high melt rates near the grounding line of PIG (Payne et al., 2007;Rignot and Jacobs, 2002) and the presence of basal channels in the ice (Bindschadler et al., 2011;Mankoff et al., 2012;Vaughan et al., 2012). Our observations show that the pattern of melting on PIG ice shelf is highly complex. Within 10 km of the grounding line, the melt rate is at least 100 m yr −1 . Only 20 km downstream this reduces to 30 m yr 1 . Between 2008 and 2011, basal melting was largely compensated by ice advection, allowing us to estimate an average loss of ice to the ocean of 87 km 3 yr −1 , in close agreement with 2009 oceanographically constrained estimates. Close to the grounding line, melting is concentrated in the basal channels and carves out those channels at 80 m yr −1 . Further downstream, on the central part of the ice shelf dominated by kilometres-wide longitudinal structures, melting on the keels is 30 m yr −1 faster than in the channels, which explains the gradual loss of channels in the downstream part of the ice shelf and the inversion of the surface elevation anomalies relative to free floatation. Our method does not give significant results for transverse or smaller (< 1 km) structures, but large spatial variations in melting likely occur over such smaller scales as well (Stanton et al., 2013).
The gradual regime shift in channel melt could be explained by the initial formation, near the grounding line, of buoyant meltwater plumes rising up the ice base and most efficiently entraining heat to the channel crests, and a decrease in the heat entrainment efficiency downstream as the slope weakens, the ice base shallows and the warm water source gets further away. At some stage, in this scenario, the plumes within the channels deliver less heat to the ice shelf than the warmer deeper waters bathing the channel keels. But more complex scenarios are not excluded.
With the advent of ice surface DEMs of even higher resolution (few metres) taken at regular time intervals, we can expect that the methodology developed here will reveal unforeseen details about the distribution of surface elevation changes and by inference of basal melt where the underlying assumptions are valid, thereby increasing our understanding of atmosphere-ice-ocean interaction dynamics and their temporal and spatial variability.
Our observations of the area close to the grounding line therefore indicate melt rates that are 80 % higher in channels than on neighbouring keels, and point to high spatial variability in the melt rates across the ice shelf, indicating strong modulation of ice-ocean interactions at kilometre scales. This implies that in situ observations need to be interpreted within their contextual position relative to the channels. Possibly the most important implication of this work concerns the modelling of sub-ice shelf cavities. Accurately represent-ing sub-kilometre scales using conventional ocean models is challenging even for dedicated regional studies, and will remain impossible for global coupled climate models for some time to come. One approach to solving this problem is to use unstructured computational meshes to focus the model resolution on features of interest, such as these channels (Kimura et al., 2013;Timmermann et al., 2012). A more conventional alternative would be to parameterize their effect on the larger scales that models are able to resolve. For either of these approaches to be successful, an essential pre-requisite is a detailed observational understanding of the channels, for which the present study provides a significant advance.

Calibration of SPIRIT DEMs
A number of ICESat track lines are available to calibrate the SPIRIT DEMs (Fig. A1a). However, due to large biases induced by the advection of small-scale features (e.g. channels), a direct comparison using all ICESat measurements would be incorrect. We therefore select ICESat observations taken within 5 days of the SPIRIT DEMs, thus reducing the available observations to 2 lines, both partly covering grounded or slowly moving parts of the ice stream. ICESat observations taken on the same day as the second SPIRIT DEM (13 March 2008) covered an area near the southern end of the ice front ( Fig. A1a and A1b). A direct comparison there gives biases of 3.1 ± 6.1 m and 2.5 ± 5.9 m for the SPIRIT DEMs taken in 5 January and 13 March 2008, respectively. Using the other ICESat line near the grounding line ( Fig. A1a and A1c) slightly reduces the variability, and is therefore used to estimate final calibration.
Tidal amplitude is expected to be of the order of 1 m in the Amundsen Sea in general (McMillan et al., 2011). This alone could explain the mean difference between the two SPIRIT DEMs. Given the generally unknown error in the tidal phase and amplitude that could occur in tidal models due to poorlyknown topography over the continental shelf and under the ice shelf, and the relatively large signal determined by our method (10 to > 30 m Lagrangian elevation change in about 3 yr), we believe additional vertical adjustment of the DEMs to correct for tidal effects is unnecessary.

Longitudinal channel misalignment with the ice flow
One expected feature of basal melt in channels is a concentration of melt to their left side (in the Southern Hemisphere) as the rising plume of meltwater is deflected by the Coriolis force (Gladish et al., 2012). This effect should preferentially melt the left-hand walls of the channels. Our observations lack sufficient spatial resolution and accuracy to resolve this Overlaid white-circled-coloured dots show 2008 along-track ICESat satellite observed surface elevation. Magenta-circled dots indicate the track-lines selected in b and c that were measured within 5 days of the SPIRIT DEMs, and covered more slowly evolving or grounded part of the glacier and part of the Hudson mountains. b) The black line shows along-track surface elevation measured by ICESat in March 13th, 2008 near the southern part of the glacier calving front. Red and green lines shows surface elevation along the same line interpolated from SPIRIT DEMs taken in January 5th and March 18th 2008, respectively, before calibration is applied. c) Same as b, but for the easternmost selected ICESat line measured in March 10th, 2008 and covering grounded ice. process directly, but we do observe that longitudinal channels tend to divert to the left from ice streamlines (Fig. B1), perhaps then, as a simple result of Coriolis-driven asymmetry in the oceanic melt relative to the ice advection.

Figure B1
Small scale surface elevation anomaly (colour) is overlaid by selected coloured streamlines of the ice surface velocity for 6 of the available observed years. All streamlines start at the same point each year but end wherever ice velocity observation is lacking. The choice of starting point is optimised to consistently show as many years and cover as much ground as possible.  Fig. B1. Small-scale surface elevation anomaly (colour) is overlaid by selected coloured streamlines of the ice surface velocity for 6 of the available observed years. All streamlines start at the same point each year but end wherever ice velocity observation is lacking. The choice of starting point is optimised to consistently show as many years and cover as much ground as possible.
While Lagrangian pathways for particles travelling with the ice flow have been derived assuming the ice velocity field to be constant, it is well known that the flow speed of Pine Island Glacier ice shelf has increased over recent decades. The acceleration impairs our comparison between the direction of longitudinal channels and pure advection within the ice, but our selections of streamlines are deduced from a series of ice surface velocity observations (Joughin et al., 2010) from 1996 to 2008, which covers the expected history of the ice parcels currently in PIG. Although variability in the streamlines suggests that any single streamline does not exactly represent a pathway for Lagrangian particles travelling with the time-varying ice flow, all available streamlines tend to divert to the right of the main longitudinal channels in the centre of the ice shelf. We therefore expect that our conclusions would not be significantly modified if full-Lagrangian trajectories were computed using an (unavailable) time-varying field of ice velocity.