Ice shelf basal melt rates from a high-resolution DEM record for Pine Island Glacier , Antarctica

Ocean-induced basal melting is directly and indirectly responsible for much of the Amundsen Sea Embayment ice 10 loss in recent decades, but the total magnitude and spatiotemporal evolution of this melt is poorly constrained. To address this problem, we generated a record of high-resolution Digital Elevation Models (DEMs) for Pine Island Glacier (PIG) using commercial sub-meter satellite stereo imagery and integrated additional 2002–2015 DEM/altimetry data. We implemented a Lagrangian elevation change (Dh/Dt) framework to estimate ice shelf basal melt rates at 32–256-m resolution. We describe this methodology and consider basal melt rates and elevation change over the PIG shelf and lower catchment from 2008–2015. 15 We document the evolution of Eulerian elevation change (dh/dt) and upstream propagation of thinning signals following the end of rapid grounding line retreat around 2010. Mean full-shelf basal melt rates for the 2008–2015 period were ~82–93 Gt/yr, with ~200–250 m/yr basal melt rates within large channels near the grounding line, ~10–30 m/yr over the main shelf, and ~0– 10 m/yr over the North and South shelves, with the notable exception of a small area with rates of ~50–100 m/yr near the grounding line of a fast-flowing tributary on the South shelf. The observed basal melt rates show excellent agreement with, 20 and provide context for, in situ basal melt rate observations. We also document the relative melt rates for km-scale basal channels and keels at different locations on the shelf and consider implications for ocean circulation and heat content. These methods and results offer new indirect observations of ice-ocean interaction and constraints on the processes driving sub-shelf melting beneath vulnerable ice shelves in West Antarctica.

A detailed understanding of the processes (e.g., ocean forcing, marine ice sheet instability) responsible for these observed changes, and their relative importance over time, is critical for future projections of PIG dynamics, mass loss, and contributions to global sea-level rise.

Geographic setting
The fast-flowing portion of the PIG ice shelf ("Main shelf", Figure 2D) is ~25 km wide and nearly 100 km long, with ice thickness of ~1-1.5 km near the main grounding line, and ~300-400 m near the calving front.Surface velocities over the main shelf are currently ~4 km/yr (~11 m/day) (Figure 2C).Two additional ice shelves are located to the northeast ("North shelf") and southwest ("South shelf') of the main shelf (Figure 2D), separated by ~2-4-km-wide shear margins.In general, surface velocity and ice thickness are relatively small (<100-500 m/yr and <300-700 m, respectively) over the North and South shelves, except for a fast-flowing tributary of the South shelf with velocity of ~1 km/yr and thickness of ~1 km near the grounding line (Figure 2).Total shelf area in recent decades varied between ~5500-6000 km 2 , due to changes in the grounding line and calving front positions.The full PIG catchment (Figure 2A) covers ~1.8-2.0x10 5 km 2 with annual surface mass balance (SMB) estimates of ~68+/-6 Gt/yr (Medley et al., 2014).
The surface of the PIG shelf is characterized by a series of longitudinal (approximately along-flow) ridges/troughs near the centerline and transverse (cross-flow) ridges/troughs toward the lateral margins that correspond to basal keels/channels (Vaughan et al., 2012) (Figure 2D).
The sub-shelf bathymetry shows a large transverse seabed ridge with relief of ~400 m above the adjacent seafloor (Figure S1).
This ridge has been the site of intermittent grounding since the mid-1940s (Smith et al., 2016), and it affects circulation within the cavity, effectively blocking some of the deep, warm CDW from entering the inner cavity (De Rydt et al., 2014;Dutrieux et al., 2014b).We further subdivide the main shelf into "inner" and "outer" regions relative to the transverse seabed ridge (Figure S1).The "ice plain" (e.g., Thomas et al., 2004) mentioned throughout the text describes a region over the inner shelf with relatively smooth, gently sloping bed (Figure S1).The lightly grounded "ice plain" was the site of significant grounding line retreat from ~1990s to ~2008, with average rates of ~1 km/yr (Park et al., 2013;Rignot et al., 2014).Our DEM record begins near the end of this progressive retreat, when the "ice plain" region was afloat except for a few isolated grounded spots (Joughin et al., 2016).

Oceanographic setting
Westerly surface winds near the continental shelf edge drive northward Ekman transport of surface water away from the continent.This draws deep, relatively warm CDW onto the continental shelf where it flows toward Pine Island Bay along two broad bathymetric troughs carved by previous glacial advances (e.g., Jakobsson et al., 2012;Kirshner et al., 2012).
The circulation pathway beneath the PIG ice shelf is less certain, but should generally be clockwise in nature, with modified CDW inflow at depth along the north side of the outer cavity, and outflow of relatively fresh meltwater along the south side of the outer cavity (Dutrieux et al., 2014b).Deep, inflowing water that encounters the large transverse seabed ridge is likely diverted to the south, flowing alongside the ridge within the outer cavity and moving toward the South cavity.Water at intermediate depth is expected to overtop the seabed ridge, creating a sharp density front and a northward jet at the ridge crest.
Eventually, these waters continue down local bathymetric slopes within the inner cavity toward the grounding line.Once in the inner cavity, the dense, modified CDW reaches the grounding line (Jenkins et al., 2010), with expected cyclonic (clockwise) circulation along the main shelf grounding line, and fresh, buoyant meltwater outflow along the centerline and south side of the shelf closing the circulation loop.The temporal evolution of this general circulation pattern, and exchange between the inner, outer, and South shelf ocean cavities depends on a number of factors, including cavity geometry defined by the evolving ice shelf base and grounding line.
Past studies offer a general picture of PIG basal melt rate spatial distribution, with relatively high rates (>100 m/yr) near the main shelf grounding line and lower rates over the outer shelf (Bindschadler et al., 2011;Dutrieux et al., 2013;Payne et al., 2007).Little is known, however, about basal melt rate temporal variability.Bindschadler et al. (2011) concluded that transverse channels/keels formed annually near the grounding line due to seasonal variability in available ocean heat content (Thoma et al., 2008;Webber et al., 2017), while simulations by Sergienko (2013) showed that similar features may be a spontaneous byproduct of the coupled ice-shelf-plume system with constant ocean heat content.

Data and methods
We present high-resolution surface elevation observations to investigate the spatial and temporal evolution of PIG.The following sections describe data sources and relevant processing methodology.

Elevation data
We use elevation data from a number of sources, including stereogrammetric DEMs and satellite and airborne altimeters.

WorldView/GeoEye stereo DEMs
We generated DEMs from high-resolution commercial stereo satellite imagery (DigitalGlobe WorldView-1, WorldView-2, WorldView-3, and GeoEye-1) using the NASA Ames Stereo Pipeline (ASP, (Beyer et al., 2018)) and methodology described in Shean et al. (2016).A total of ~3000 along-track stereopairs from October 2010 to May 2015 were processed for the West Antarctic coastline (Figure 1).For this study, we focus our analysis on a ~260x240 km region with dense WorldView/GeoEye DEM coverage covering the PIG shelf and lower trunk (Figure 1C).
The Level-1B (L1B) images were orthorectified using a smoothed version of the Bedmap2 surface DEM (Fretwell et al., 2013) before stereo correlation.Processing settings for ASP (see Shean et al., 2016) include "seed-mode 3" (sparse_disp utility) to initialize the correlation, a 2-level correlation pyramid limit, a correlation timeout of 360 seconds, parabolic sub-pixel refinement, and filtering of isolated disparity map clusters with area <1024 pixels.
We generated additional "cross-track" or "coincident mono" DEMs from pairs of independent mono images with geometry suitable for stereo reconstruction.We identified candidate pairs in the DigitalGlobe image archive based on the criteria in Table 1, and generated 24 DEMs from images acquired between October 2011 and January 2012.Some of these cross-track pairs were acquired on the same orbit, while others were acquired on different orbits, sometimes by different spacecraft.Final time offsets between the images ranged from 0.007 to 1.6 days.
The cross-track DEMs potentially have increased error due to residual horizontal displacement errors (i.e., errors due to ice flow between image acquisitions), non-ideal stereo geometry (e.g., smaller convergence angles) and the fact that some errors in ephemeris data for the two images are independent (as opposed to highly correlated errors for along-track pairs).In practice, these issues can result in increased DEM vertical/horizontal bias and increased relative error (e.g., more "tilt").Despite potentially increased error, we include these cross-track DEMs in our analysis to fill critical gaps in coverage near the PIG grounding line, and to increase overall DEM sample size for the 2011/2012 season.As described in Section 2.2, these errors are mitigated through subsequent DEM co-registration and correction.

SPIRIT DEMs
We incorporated all six available SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies,  (Bouillon et al., 2006;Korona et al., 2009), which we improve substantially using control points as described in section 2.2.1.2.
We used the DEM V1 products (generated with correlation parameters tuned for gentle slopes), applied the corresponding "CC" mask to preserve correlation scores between 50-100% (masking most interpolated areas), reprojected to a standard polar stereographic projection (EPSG:3031), and removed the EGM96 geoid offset to obtain elevations relative to the WGS84 ellipsoid.We filtered the resulting products to remove isolated pixels, mask elevations <20 m above sea level, and remove any pixels with >30 m absolute elevation difference from the per-pixel median of 2010-2015 WorldView/GeoEye DEMs, effectively removing spurious values associated with clouds in the original imagery.

Satellite and airborne altimetry
The NASA Operation IceBridge (OIB) mission has collected airborne altimetry data over PIG during annual campaigns from 2009/2010 to 2014/2015, except for the 2013/2014 season.Most campaigns occurred during October-November, with data acquisition flights for a particular site typically occurring over ~1-3 days.We assembled all available NASA Airborne Topographic Mapper (ATM, (Krabill et al., 2002;Martin et al., 2012)) and Land, Vegetation, and Ice Sensor (LVIS, (Blair et al., 1999;Hofton et al., 2008)) airborne lidar data for use in our study area.A total of 25 ATM flights and 7 LVIS flights We processed all altimetry data as described previously (Shean et al., 2016), and produced gridded 32-m and 256-m DEMs with sparse coverage for each campaign using the ASP point2dem utility.This utility assigns the output value for each grid cell by computing the weighted mean of all points within a 1-grid-cell-width radius.
We also include available 2003-2009 NASA ICESat-1 Geoscience Laser Altimeter System (GLAS, (Schutz et al., 2005;Zwally et al., 2002)) satellite altimetry data.These data were clustered by ~33-day campaign and gridded as described above, providing 18 additional sparse DEMs.While caution must be exercised during interpretation of these sparse data over rough surfaces or steep slopes, we include them in our analysis to extend the observational record to 2003.

DEM co-registration and correction
The following sections describe a cascading co-registration and correction workflow used to improve both absolute and relative DEM accuracy over the PIG study area.

Co-registration with altimetry
Where possible, a point-to-point iterative closest point (ICP) algorithm (Shean et al., 2016) was used to co-register DEMs to filtered altimetry data from the sources described in Section 2.1.3.The altimetry data were queried for each DEM extent and the returned points were limited to "static" (e.g., nunataks) and "dynamic" (e.g., slow-moving ice with limited slope/roughness) control surfaces.Any points over floating portions of the PIG shelf were excluded.The remaining points were further filtered using a maximum expected displacement (product of measured surface velocity magnitude and time offset between the point and DEM timestamp) threshold of 10 m and maximum DEM-point time offset of 1 year.All control points were assumed to have vertical accuracy of ~0.1 m.

WorldView/GeoEye DEM co-registration
The majority of the WorldView/GeoEye DEMs had 10 6 -10 8 filtered points available for co-registration, with DEM-altimetry time offsets of only a few months.The ICP co-registration provided translation corrections for 368 of 575 DEMs over the PIG catchment.Figure 3 and Table 2 show a significant improvement in multiple quality metrics following co-registration.
Uncorrected DEMs had an initial mean vertical bias of +3.1 m above the altimetry data (Figure 3), and we applied a -3.1 m vertical correction to the remaining 207 DEMs that lacked adequate control data.

SPIRIT DEM co-registration
The filtered SPIRIT DEMs were co-registered with the ICP routine described in Section 2.2.1, and the results are shown in

Elevation correction for ocean and atmospheric variability
After DEM co-registration, we corrected all elevation data (including altimetry) over the floating portions of the PIG shelf to remove the effects of ocean tides, atmospheric pressure (Inverse Barometer Effect (IBE), (e.g., Padman et al., 2003)) and mean dynamic topography.
We computed tidal amplitude ∆ℎ # using the CATS2008A inverse barotropic tide model (an updated version of the model described in Padman et al. (2002)).The inverse barometer effect magnitude ∆ℎ %&' was computed from 6-hour interval ERA- Interim mean sea level pressure reanalysis data (Dee et al., 2011).We removed the 2002-2016 median pressure (985.21hPa), and scaled residuals by ~1 cm/hPa to obtain the approximate inverse barometer correction.Tidal amplitude for DEM timestamps ranged from -0.75 to +1.04 m ( = 0.33 m), while the inverse barometer effect amplitude ranged from -0.3 to 0.35 m ( = 0.11 m) (Figure S3).These high-frequency (hourly-daily) corrections show good agreement with observed surface elevation records from GPS receivers on the PIG shelf (Shean et al., 2017).
The mean dynamic topography (∆ℎ /01 ) correction removes residual offsets between the geoid and mean sea level due to ocean circulation.Estimates for mean dynamic topography near ASE are approximately -1.2 m (Andersen and Knudsen, 2009).
Corrected ice surface elevation above sea level is calculated as: Where ℎ 3 is measured elevation above the WGS84 ellipsoid and ∆ℎ 5 is the EGM2008 geoid offset (Pavlis et al., 2012) (approximately -27.6 to -24.4 m across PIG shelf).We defined the coefficient  to increase linearly with distance l beyond the grounding line: which provides a smooth transition from floating to grounded ice for these corrections.For this study, the grounding line (Figure 2) was defined with a single composite polygon derived from DInSAR (Joughin et al., 2016;Rignot et al., 2014) and high-resolution DEM data, with an approximate timestamp of 2011.
After correction using Equation (1), surface elevation from airborne altimetry approaches 0 m above sea level over open water.

WorldView/GeoEye DEM "tilt" correction
As identified in Shean et al. (2016), a subset (~5-10%) of the WorldView/GeoEye DEMs appear to have a slight along-track and/or cross-track "tilt" of ~1-3 m over the ~111 km strip length, likely due to small errors in spacecraft attitude metadata.For most of these "tilted" DEMs, the available control point spatial distribution is insufficient to constrain a rigid-body ICP rotation.
Initial attempts using bootstrapping and least-squares minimization of offsets between adjacent, overlapping DEMs to solve for a "tilt correction" failed due to overfitting and the propagation of larger errors near some DEM edges.To correct these problematic DEMs, we developed an optimization approach that simultaneously solved for interannual dh/dt and planar corrections to remove individual DEM tilt.In principle, this is similar to the SERAC method used for altimetry over the Greenland ice sheet (Csatho et al., 2014;Schenk and Csatho, 2012).

2017
).Thus, while the dynamic response to earlier rapid grounding line retreat and speedup continues to propagate upstream across the PIG catchment, we expect relatively limited variability in elevation change rates during this period.
We manually masked the main shelf and fast-flowing trunk within ~30 km of the grounding line, and then used the criteria listed in Table 3 to identify "dynamic control surfaces" (Figure 4) over grounded ice with limited linear trend (dh/dt) and limited residual variance about this trend.Over these surfaces, the elevation at any particular DEM pixel i (with spatial coordinates xi and yi) at time j is given by: ℎ @,A = ( @  A +  @ ) + ( A  @ +  A  @ +  A ) where ai and bi represent the slope and offset of a linear model fit to elevation values at pixel i, and coefficients cj, dj and ej define a planar correction for all i within a DEM at time j.
We solved for these coefficients using least-squares minimization with regularization and a smoothness constraint designed to penalize large spatial gradients.Elevation values from filtered, gridded altimetry data were included in the solution with increased weight.Stereo DEMs with <40 km along-track length were limited to a vertical offset correction (ej), with no tilt correction (cj = dj = 0).Tilt tolerances were increased for cross-track DEMs (Section 2.1.1)and vertical offset tolerances were increased for input DEMs that were not initially co-registered using ICP.Tilt tolerance was limited in the DEM cross-track direction, as most of the observed tilt was in the DEM along-track direction.Figure 4 and Table 4 summarize the results of these corrections, with considerable improvement in all metrics.

Output elevation data
We prepared a resampled "stack" of all co-registered, corrected DEMs over the PIG shelf using a common 256-m grid.
Additional stacks with increased grid resolution (64-m and 32-m, respectively) were prepared over high-priority areas (i.e., inner shelf and GPS validation sites (Shean et al., 2017)).

Post-correction DEM accuracy
As discussed in Shean et al. (2016), the uncorrected vertical/horizontal accuracy for the along-track stereo DEMs is <5.0 m CE90/LE90.After systematic artifact removal and co-registration, vertical accuracy can be less than <0.2-0.4 m for surfaces with <10° slope.For the PIG shelf, we conservatively estimate the final DEM accuracy to be ~1 m after co-registration and least-squares "tilt" correction.We initially expect increased uncertainty for 2013/2014 DEMs due to reduced availability of OIB altimetry data during this season.This uncertainty, however, was reduced after the least-squares correction, which leveraged altimetry data and corrected WorldView/GeoEye DEMs from adjacent years.
Several factors can reduce the effectiveness of DEM co-registration with altimetry.The primary problems for PIG include sparse control data with limited variation in surface slope and aspect, and longer offsets between the DEM and altimetry timestamps (~1-12 months).Over these timescales, surface processes (e.g., accumulation/ablation, wind redistribution of

Annual surface elevation composites and mosaics
We generated weighted-average composites using the ASP dem_mosaic utility for all available elevation data in a given year (September-April, but typically October-March), with a nominal January 1 timestamp (Figure 5).For each output pixel, the weighted averaging algorithm assigns greater weight to input pixels from spatially continuous sources (e.g.DEMs with few data gaps) and penalizes isolated pixels or clusters of pixels (see ASP documentation for details).The resulting composites appear seamless, but can include smoothing artifacts due to variable temporal sampling of input elevation data, especially for features that advect in the along-flow direction.
Adjacent WorldView/GeoEye stereo DEMs are often acquired weeks or months apart during a particular season due to clouds and/or competition for resources.Even after DEM co-registration and correction, this sampling can introduce horizontal and vertical feature offsets between adjacent DEMs in fast-flowing regions.Generally, this temporal sampling is not a problem for smaller targets covered by a single WorldView/GeoEye DEM footprint (e.g., Greenland outlet glacier termini).Larger targets like the PIG shelf, however, require >10 WorldView/GeoEye DEMs for complete coverage, and more sophisticated mosaicking approaches are necessary to preserve local features.
To obtain full-shelf coverage while also preserving timestamps and relative elevation values within individual input DEMs, mosaics without averaging or blending were generated for the ~October-March period each year.We used a "reverse" ordering scheme for input DEM timestamps, so that the last DEM from each season was mosaicked on top.Finally, we generated WorldView/GeoEye DEM mosaics when complete shelf-wide coverage was available over a relatively short time span (e.g., October-December 2012, Figure 2D).In such cases, input DEM products were manually selected and ordered to minimize feature offsets.

Surface velocity
Surface velocity data constrain horizontal ice-shelf advection rates and aid interpretation of observed elevation change.In an effort to generate self-consistent velocity and DEM products, we estimated velocity using feature tracking with normalized cross-correlation of two DEMs, similar to the approach described by Dutrieux et al. (2013).However, this approach is susceptible to spurious correlations and data gaps over flat, featureless areas, especially for low-resolution inputs (e.g., 40-m SPIRIT DEMs).This technique also fails for longer time intervals (>2 years), as surface processes, deformation and rotation due to velocity gradients, and spatially variable basal melt decreased coherence.For these reasons, we used an independent set of gridded velocity products, which enabled reconstruction of particle paths for arbitrary elevation data, including sparse altimetry.
We derived spatially and temporally continuous velocity fields for the full PIG shelf using piecewise linear interpolation via 3-D (x,y,t) Delaunay Triangulation.Linear barycentric interpolation was then used to extract spatially continuous velocity grids with 512-m resolution for a regular time interval of 122 days from January 1, 2008 to June 1, 2016.The interpolated velocity products were smoothed in the time dimension with a 610-day, 2 nd -order Savitzy-Golay filter, and then in the spatial dimension with a 2.5-km rolling median filter to mitigate artifacts in the input mosaics.To obtain velocity fields with increased spatiotemporal sampling, we performed secondary interpolation with a high-resolution timestep (e.g., 5-20 day) and increased spatial sampling (e.g., 32-256 m), with a final Gaussian smoothing filter (~0.17-km sigma) applied in the spatial dimension to reduce any residual interpolation artifacts.The basal melt rate calculations described in Section 2.8.2 required estimates of the velocity divergence, which we calculated from these interpolated, smoothed velocity products for each high-resolution time step using a central-difference approach.
We produced a new combined bed dataset using aerogravity/Autosub data, existing open-water bathymetry, and all available quality-controlled CReSIS MCoRDS and British Antarctic Survey (BAS) Polarimetric Airborne Survey Instrument (PASIN) ice thickness measurements collected over grounded ice.We fit a smooth surface to these data using an inversion procedure that preferentially minimizes bed curvature in the along-flow direction, while matching the bed elevation at data points to within the estimated data errors (see methods of Medley et al., 2014).While some local "peaks" over the longitudinal seabed ridge beneath the PIG shelf may be biased high, this bed appears most consistent with observed recent grounding line evolution (Joughin et al., 2016).

2.7
Surface mass balance (SMB) The Regional Atmospheric Climate Model (RACMO) v2.3 (Ettema et al., 2009;Lenaerts et al., 2012;Van Meijgaard et al., 2008;Van Wessem et al., 2014) provides continent-wide estimates of surface mass balance on a 27-km grid.To estimate SMB over the PIG shelf, we used monthly average SMB products available through December 2013, and repeated the observed 2013-2014 SMB signal for calculations spanning 2014-2015.The RACMO SMB data were interpolated over the study area to generate gridded products with the same extent and spatial sampling as the DEM and velocity products.

Elevation change
We consider elevation change for PIG using both Eulerian dh/dt (fixed reference grid) and Lagrangian Dh/Dt (grid moving with the surface) descriptions.These two approaches are complementary and provide distinct information over grounded and floating ice.

Theory
Assuming incompressibility, constant ice density, and column-average velocity u, the Eulerian description of mass conservation for a column of ice with ice-equivalent thickness H (after removing firn-air content d) can be expressed as: where ̇ is surface accumulation rate (meters ice equivalent for time interval dt) and  ̇ is basal melt rate (meters ice equivalent, defined as positive for melt).
The flux divergence term, ∇ • (), can be expanded as: where ∇ •  is the velocity divergence (positive for extension) and ∇H is the thickness gradient.
The relationship between Lagrangian (denoted by material derivative operator 0 0# ) and Eulerian thickness change is provided by the material derivative definition: Equations 4, 5, and 6 can be combined to obtain Lagrangian thickness change for the column: Over grounded ice, we assume that the bed elevation remains constant, and can substitute Eulerian surface elevation change dh/dt for Eulerian thickness change dH/dt.This assumption does not hold for floating ice.If we assume hydrostatic equilibrium, however, we can estimate freeboard ice thickness from observed surface elevation.We remove firn-air content d from observed surface elevation h to obtain ice-equivalent freeboard surface elevation, and then compute ice-equivalent freeboard thickness:

Eulerian long-term interannual trend
To characterize long-term (~5-10 year) elevation change over the PIG shelf during and after the period of rapid groundingline retreat, we computed interannual per-pixel trends for the 2003-2010 and 2010-2015 periods.These trends were determined using a linear fit to surface elevation for each grid cell with 3 or more observations, with >6 valid samples for most cells.No smoothness constraint was imposed -all fits were computed independently, although adjacent elevation values are highly correlated.

Basal melt rate
Both Eulerian and Lagrangian frameworks can be used to estimate the basal melt rate.The Lagrangian description tracks elevation change for the same column of ice over time, eliminating potential aliasing due to advection of high-frequency surface gradients (i.e., ice-shelf surface ridges and troughs).If velocity divergence and surface accumulation are known, Equation 9 can be rearranged to solve for the component of observed elevation change due to basal melt:

Lagrangian Dh/Dt basal melt rate implementation
Past studies of Lagrangian Dh/Dt basal melt rates used in situ observations (e.g., Jenkins et al., 2006), a single pair of gridded DEM observations (e.g., Dutrieux et al., 2013), or a series of sparse altimetry data (e.g., Moholdt et al., 2014).The approach presented here uses hundreds of independent DEM observations with variable spatial coverage over an 8-year time period.
The observed cumulative particle Dh/Dt is then used to estimate evolving surface elevation h at each time step n along the particle path (assuming the Dh/Dt rate is constant), and local velocity divergence (∇ • ) values are sampled at each time step n along the particle path from the continuous velocity products described in Section 2.5.The corresponding local ℎ(∇ • ) is then integrated over the full path: This approach should accurately capture time-variable thinning/thickening due to local velocity divergence experienced along each path, rather than sampling velocity divergence from single, fixed velocity grid.We also sampled time-variable SMB grids at each time step, but in practice, the spatiotemporal variability for the monthly 27-km products is limited along the ~8 km particle paths, and we use a time-averaged estimate for ̇ extracted at the particle path midpoint.Finally, we substitute the cumulative particle Dh/Dt and local ℎ(∇ • ) into Equation 10, which provides an integrated basal melt rate estimate for a single particle from a single pair of DEMs.

Basal melt rate path distribution
We consider two end members for the spatiotemporal distribution of ice shelf basal melt rates.End member #1 assumes a horizontally and vertically variable melt rate field that does not vary over time, so that features with variable draft (i.e., keels and channels) melt at different rates as they advect through this field.End member #2 assumes that melt rate spatial variability is highly correlated with local shelf thickness gradients (and associated basal slope), so that local melt rates advect with features on the shelf (e.g., once formed, a transverse basal channel will continue to melt rapidly as it advects downstream).In reality, basal melt rates are likely sensitive to some combination of these two end-member scenarios.
The methodology described in Section 2.9.1 provides basal melt rate estimates for each particle in a Lagrangian reference frame.For subsequent analysis on a regular grid, we must remap these observations into a common, global Eulerian reference frame.This step is complicated by the fact that the long time intervals between DEM observations (~2 years) and high advection rates (~4 km/yr) on the main PIG shelf result in particle path lengths (~8 km) that greatly exceed the input DEM grid cell size and the desired melt rate grid cell size (256 m).To address this issue and to evaluate the two basal melt rate endmember scenarios, we developed two approaches to work with observed Lagrangian Dh/Dt basal melt rates in an Eulerian reference frame: "along-flow distribution" and "initial-pixel" (Figure 6).

Along-flow distribution products
The "along-flow distribution" approach partitions observed particle basal melt rates (Section 2.9.1) evenly across each path, and computes statistics for each cell in a fixed Eulerian grid using all paths that pass through that cell (Figure 6).This approach potentially provides a more realistic map of the melt rate "field" (end member #1), but it effectively smooths basal melt rate estimates in the along-flow direction, especially for longer path lengths.This leads to reduced resolving power for basal melt rate spatial variability (end member #2), especially for features with transverse orientation (e.g.transverse channels/keels, rifts).
The path history of all valid particles for a particular DEMi-DEMj combination is reduced to identify a unique set of occupied grid cells in the global Eulerian reference frame.For each particle path, basal melt rate  ̇ is calculated as described in Section 2.9.1 and these values are distributed evenly along encountered cells.This procedure yields a spatially variable particle count within each cell in the global Eulerian coordinate system; only one particle will pass through a cell on the upstream edge of the domain, while ~10-100 unique particles could pass through a cell near the center of the domain over the full Dt interval.
We then compute the median and NMAD for each cell (Figure 6).This approach reduces noise and provides metrics to evaluate variance and uncertainty in derived basal melt rates.

Initial-pixel products
The "initial-pixel" approach assigns particle basal melt rate values to the corresponding path origins in DEMi, so the resulting basal melt rates grids have the same spatial extent as DEMi.This approach is relatively straightforward, and was used in earlier work (e.g., Dutrieux et al., 2013).It preserves the relative spatial distribution of basal melt rates across individual features in DEMi (e.g., channels and keels) to evaluate end member scenario #2, but doesn't resolve where along the ~8 km particle path that melt actually occurred.
For a given DEMi-DEMj combination, the initial-pixel approach assigns particle basal melt rate values to DEMi pixel locations.
For each initial DEMi, we then create stacks of available DEMi-DEMj initial-pixel basal melt rate products and compute a perpixel stack median map.In other words, basal melt rates calculated from each valid downstream DEMj are assigned to the initial DEMi pixel locations, and median values for each DEMi pixel are computed assuming no temporal variability in basal melt rates for all valid ti-tj intervals.

Path distribution considerations
Under melt-rate end-member #1, the initial-pixel approach will introduce a negative bias for a fixed basal melt rate field with relatively large negative spatial gradient (e.g., 200 m/yr to 100 m/yr over 8 km in the inner cavity), as the mean path basal melt rate (150 m/yr) will be assigned to the initial-pixel locations (where rates are locally 200 m/yr).We experimented with an approach using path mid-point locations rather than initial-pixel locations, but this resulted in large gaps near the grounding line and prevented direct comparison of basal melt rates with the original DEMi elevations.Under melt-rate end-member #2, the initial-pixel approach provides more realistic basal melt rate magnitude and spatial distribution than the along-flow distribution approach.The difference between the two approaches will be negligible for areas of the PIG shelf with low surface velocity (<250 m/yr).We present results for both approaches and evaluate the two melt-rate end-member scenarios.

Shelf-wide basal melt rate composites
In the above sections, we described basal melt rate calculations for a single DEMi-DEMj combination with sufficient overlap and a ti-tj time interval that falls within the chosen Dt range (~2 years), which represents only one of many potential valid DEMi-DEMj combinations that can be formed from the full set of DEMs over the PIG shelf.
For a given DEMi, after we calculate basal melt rates using the first viable DEMj, the DEMi particles are further propagated and the process is repeated for all other viable DEMj until the ti-tj time interval exceeds the maximum Dt interval.The entire process is then repeated for all possible DEMi.The individual DEMi-DEMj basal melt rate products can have relatively high uncertainty and/or limited spatial extent, so we created annual melt-rate mosaics and composites to reduce noise and increase total spatial coverage.We used different methodology for the "along-flow distribution" and "initial-pixel" approaches, as described below.

Along-flow distribution composites
We generated weighted-average basal melt rate composites from individual "along-flow distribution" basal melt rate products.
This approach provides basal melt rate grids centered on January 1 for the ~2-year interval products.For each grid cell in the output mosaics, the weighted-average approach favors pixels near the center of input products with larger areal coverage.Perpixel standard deviation is also calculated for each ~2-year basal melt rate composite, providing maps that capture the spatial distribution of basal melt rate uncertainty (and any true basal melt rate temporal variability during the ~2-year period).The annual composites were then used to generate a mean basal melt rate composite for the full 2008-2015 period.

Initial-pixel mosaics
The per-pixel stack median products described in Section 2.9.2.2 provide high-resolution maps of local basal melt rates, but they are limited to the DEMi spatial extent.To overcome this limitation, we generated mosaics of the stack median products using a reverse time ordering scheme, so basal melt rate estimates for the most recent DEMi timestamp were mosaicked on top.This approach preserves the local basal melt rate distribution within each stack median product, while providing coverage over as much of the shelf as possible, with limited time offset between spatially adjacent observations.These products can be directly compared with surface elevation (and corresponding freeboard thickness estimates) from the "reverse" ordering DEM mosaics described in Section 2.4.

Uncertainty and sources of error
Surface elevation uncertainty over the PIG shelf includes errors due to the geoid model (0.1-0.4 m), mean dynamic topography (0.2 m), and tide/IBE correction (0.1 m).For simplicity, we assume a constant firn-air content of 12 m with uncertainty of 2 m to account for spatial and temporal variability (see Appendix of Shean, 2016).We used a depth-averaged density for ice and underlying ocean water of 917 +/-5 kg/m 3 and 1026 +/-1 kg/m 3 , respectively, and assume that these densities are constant in both space and time.We assume uncorrelated errors of 1 m for surface elevation, 50 m for bed elevation, 30 m/yr for velocity (for ~37.5° look angle and +/-0.5 m tide) (Joughin, 2002) and 28% for SMB (Depoorter et al., 2013).
Our conversion from surface elevation to ice thickness assumes that the ice shelf is in hydrostatic equilibrium (Shean et al., 2017).We use a consistent methodology and the same assumptions of hydrostatic equilibrium for the full 2008-2015 study time period, which increases confidence in observed temporal change.We do not update the grounding line mask for basal melt rate calculations, and some of the persistent high and low basal melt rate values <1-2 km downstream of the grounding line may be related to evolving grounding line position and insufficient masking over grounded ice (Joughin et al., 2016;Milillo et al., 2017).Transient regrounding of keels will yield increased surface elevations and larger apparent freeboard thickness values.This may also lead to localized ice deformation and non-zero vertical strain rates that are inconsistent with the assumption that surface velocity equals the column-average velocity.
Uncertainty for elevation change and basal melt rate products depends on the time interval.For example, assuming that errors are uncorrelated, a 1-m absolute error in surface elevation should result in ~1.4 m root-mean-square error in elevation change.
This elevation change uncertainty should remain constant, so integrating observations over longer periods will result in greater signal-to-noise for annual elevation change rates (e.g., ~1.4 m/yr uncertainty for a 1-year interval or ~0.7 m/yr for a 2-year interval, assuming constant rates).This estimate does not, however, include slope-dependent vertical error due to cumulative horizontal displacement error, which will increase for longer time intervals.It is challenging to quantify this Dh/Dt uncertainty contribution in a forward sense, as multiple sources (e.g., cumulative displacement error from velocities, DEM co-registration, DEM resampling) can lead to slope-and aspect-dependent errors.Basal melt rate products can also include artifacts over shear margins and near the ice front due to anomalously large Dh/Dt values (+/-20-40 m) from advection of near-vertical surface gradients (e.g., ice front, icebergs, rifts) and errors in velocity divergence.
The stacking and averaging approaches described in Section 2.9.3 should reduce many of these errors, but this improvement is difficult to capture with formal error estimates.The initial-pixel stack per-pixel NMAD (Section 2.9.2.1) and along-flow per-pixel standard deviation (Section 2.9.3.1)metrics can provide maps of uncertainty, but these estimates will also include any true basal melt rate temporal variability during the observation period.We also note significantly increased 2010-2015 thinning rates (~3-4 m/yr) over upstream ice stream shear margins within ~60 km of the grounding line, especially the north shear margin.A series of curvilinear elevation anomaly "bands" with orientation approximately transverse to flow is apparent over the catchment ~40-100 km upstream of the grounding line (Figure 7D).

Basal melt rate spatial distribution
Figure 8 shows mean 2-year Lagrangian Dh/Dt basal melt rate products for the 2008-2015 period.Full ice-shelf basal melt rates were ~82 Gt/yr for "initial-pixel" and ~93 Gt/yr for "along-flow distribution" composite 2-year products.
In general, basal melt rates are >150-200 m/yr near the main shelf grounding line, with highest rates of >250 m along the north side of the grounding line (Figure S4).Basal melt rates are generally ~50-100 m/yr over the main shelf inner cavity, where ice thickness exceeds ~600-700 m, and ~10-30 m/yr over most of the outer shelf, where ice thickness is ~300-500 m.We observe considerable anisotropy, with longitudinal spatial correlation over lengths scales of ~20 km and significant transverse ~km-scale variability.This is true for both the "initial-pixel" and "along-flow distribution" products (Figure 8), suggesting that this anisotropy is not a result of smoothing in the along-flow direction.The northern third of the outer main shelf displays ~3-4 longitudinal features with elevated basal melt rates of ~30-40 m/yr (red arrow in Figure 8).Upstream of these features, a broad (~10 km wide x 20 km long) region of low-relief transverse ridges/troughs displays reduced basal melt rates of ~5-10 m/yr (green arrow).
Basal melt rates are ~0-10 m/yr over the South shelf and ~0-5 m/yr over the North shelf (Figure 8).High basal melt rates of ~60-90 m/yr are observed near the relatively deep (~900 m) grounding line of the fast-flowing (~0.7-1.0 km/yr) South shelf tributary.Elevated basal melt rates of ~20-50 m/yr are also observed within large channels on the south shelf (blue arrow in Figure 8).Integrated basal melt rates over the North and South shelves are ~5 and ~10 Gt/yr, respectively.

Channel-scale melt distribution
We used the "initial-pixel" basal melt rate mosaics to evaluate observed basal melt rates for basal channels and keels on the main shelf.We applied a high-pass filter (1.5 km sigma Gaussian) to annual "reverse" order DEM mosaics (Section 2.4), and defined masks for channels and keels using filtered elevations less than -1 m and greater than +1 m, respectively (Figure 9).These masks were applied to corresponding 2-year "initial-pixel" basal melt rate products, and weighted-average composites were generated from all available years to document the spatial distribution of main shelf channel and keel melt rates for the full 2008-2015 period.The value at any given pixel in the channel (keel) composite is derived from melt rates for several advecting channels (keels) that intersected that pixel over time, providing a sample of background melt rates for channel (keel) features at different locations in the cavity.
The highest basal melt rates are associated with longitudinal surface ridges (basal keels) within ~3-4 km of the grounding line.
In the inner cavity (between ~4-15 km from the grounding line), high basal melt rates (>100 m/yr) are associated with both longitudinal surface troughs (basal channels) and surface ridges (basal keels).Several persistent channels display high basal melt rates throughout the 2008-2015 record, but there is more apparent temporal variability associated with deep keels due to grounding and ungrounding.
Over the mid to outer shelf, we observe relatively high basal melt rates over keels (~20-40 m/yr) and limited basal melt rates in transverse channels (~0-10 m/yr).We note that both channels and keels display higher basal melt rates over the northern portion of the outer shelf (red arrow in Figure 8).Higher basal melt rates of ~10-20 m/yr are observed over ~50-70 km-long longitudinal keels near the shelf centerline, while ~0 m/yr basal melt rates are observed within adjacent longitudinal channels.
One prominent longitudinal keel displays basal melt rates of ~30-40 m/yr (black arrow in Figure 8).

Long-term elevation change
Grounding line retreat and speedup through 2010, combined with inherent marine ice sheet instability, are primarily responsible for the strong thinning observed upstream of the grounding line at PIG (Joughin et al., 2010) (Figure 7).Our observations show that this thinning decreased after 2010, which is consistent with results from model simulations documenting the inland migration of the associated speedup (Joughin et al., 2010).The end of rapid grounding line retreat and the regrounding of deep keels on the seabed ridge (Christianson et al., 2016;Joughin et al., 2016) likely contributed to decreased thinning rates immediately upstream of the grounding line after 2010.The continued thinning over upstream shear margins (Figure 7) can also be attributed to this evolution, as sustained thinning rates of >5-10 m/yr over the main trunk prior to 2010 (Flament and Rémy, 2012;Joughin et al., 2010;Wingham et al., 2009) led to an increase in surface slopes and transverse driving stress across the shear margins.

Basal melt rate spatial distribution
Our results show a ~11 Gt difference between the shelf-wide along-flow distribution and initial-pixel basal melt rate estimates, with most of this difference over the inner cavity.This discrepancy is likely related to steep spatial gradients in the "fixed" melt rate field (end member #1), which we would expect to introduce a negative bias in the initial-pixel basal melt rate estimates, as described in Section 2.9.2.Thus, the along-flow distribution melt-rate of ~93 Gt/yr estimate is likely a better shelf-wide estimate.The along-flow distribution and initial-pixel basal melt rates are comparable on the outer shelf and slowmoving areas of the North and South shelves, with both offering good "resolution" of basal melt rates for individual features (e.g., channels and keels).The spatial distribution of high basal melt rates near the grounding line is likely a function of modern (post-2006) cavity geometry (Figure S1) and sub-shelf circulation.Mass conserving bed reconstruction for the 1990s configuration revealed a large longitudinal seabed ridge (~4 km wide x 30 km long) near the centerline of the inner cavity (Rignot et al., 2014).The highest basal melt rates of >200-250 m/yr are observed on the north side of this longitudinal seabed ridge, where warm, salty water circulating at depth through the inner cavity first reaches the grounding line (e.g., Dutrieux et al., 2014b).
We note that the enhanced ~30-40 m/yr basal melt rates over the northern portion of the outer shelf (red arrow in Figure 8) are located just downstream of the transverse seabed ridge.Both the Autosub observations and ocean GCM simulations show increased ocean current velocity and enhanced variability due to cold water intrusion near this location (Dutrieux et al., 2014b), suggesting that this local high in basal melt rates could be related to local circulation patterns and/or upwelling.This location is also one of the expected pathways for warm CDW inflow into the inner cavity (e.g., St-Laurent et al., 2015), and we suggest that as this water flows over the transverse seabed ridge, it could lead to enhanced turbulence, vertical heat transport towards the ice base, and associated basal melting.

Channel-scale melt distribution
Our results are generally consistent with past work (e.g., Dutrieux et al., 2013) suggesting that increased melt occurs within deep basal channels closer to the grounding line, and on basal keels over the outer shelf (Figure 9).Inner-cavity channels/keels have much higher relief than outer-shelf channels/keels, so we might expect higher basal melt rates due to faster plume-driven flow along inner-cavity channels.However, we also observe high basal melt rates over deep keels in the inner cavity, suggesting that high heat content and local circulation can dominate melting at these depths.
Our results demonstrate the potential for high-resolution Lagrangian Dh/Dt analyses of channel-scale features on ice shelves, even with known methodological limitations (Section 2.10 and discussion in Dutrieux et al., 2013;Shean et al., 2017).Keels on the mid to outer PIG ice shelf typically reach water depths up to ~400-450 m, while channels are typically ~300-350 m.
These features should intersect the observed thermocline, with temperature gradients of over 1.0°C possible between ~300 and ~450 m depth (Dutrieux et al., 2014b).Our results are consistent with the hypothesis that enhanced melting of outer shelf keels is related to their exposure to warmer water at depth, with reduced plume-driven flow in the channels due to limited ice thickness gradients.We note that the transverse surface ridges and troughs on the south side of the main shelf display greater relief than those along the north side of the shelf, with correspondingly higher basal melt rates over the deeper keels.Based on these observations, we suggest that statistical analyses of keel melt rates over time could provide new information about the spatiotemporal evolution of the thermocline in the outer cavity.

Comparison with past basal melt-rate assessments
The local basal melt rates observed within the deep inner cavity (>200 m/yr) are significantly higher than past estimates of ~100 m/yr (Bindschadler et al., 2011;Dutrieux et al., 2013;Payne et al., 2007).They are more consistent with the Payne et al. (2007) estimates, which show high rates of ~200-300 m/yr near the mid-1990s grounding line.Furthermore, these studies relied on interpolation of sparse ICESat-1 tracks to estimate spatially continuous Eulerian dH/dt for the entire PIG shelf (e.g., -5.32+/-0.3m/yr (Rignot et al., 2013)).The ICESat-1 GLAS laser spot was ~30-70 m in diameter with ~170 m along-track spacing and ~20 km cross-track spacing between repeat tracks over PIG (e.g., Figure 3 of Pritchard et al. (2009)).Limited measurements were available to constrain local slopes sampled by repeat ICESat-1 tracks over the PIG shelf, and aliasing of advecting km-scale surface ridges and troughs can lead to significant errors in thinning rates inferred from smoothed ICESat-1 repeat tracks (e.g. Figure 13 of Sergienko ( 2013)), especially after converting inferred elevation change to freeboard thickness change.While this may not be relevant for flat, smooth ice shelves like the Ross and Ronne-Filchner ice shelves (e.g., Moholdt et al., 2014), this issue complicates analysis of the sparse ICESat-1 dh/dt measurements over the relatively rough PIG ice shelf, and previous uncertainty estimates for shelf-wide basal melt rates based on ICESat-1 observations are likely too low.Thus, while basal melt rates may have been higher from ~2003-2008, we cannot rule out the possibility that no long-term change occurred between ~2003-2008 to ~2008-2015.
We note that observations with ~20-km spatial resolution (e.g., ICESat-1 or radar altimetry (e.g., Paolo et al., 2015)) can capture the long-term temporal evolution of Eulerian elevation change and basal melt for the full PIG ice shelf, but they cannot directly capture changes associated with dynamic ice-ocean processes that operate on shorter spatial scales.The high-resolution DEM record and methodology presented here allows for both full-shelf basal melt rate estimates and analysis of the detailed spatiotemporal evolution of km-scale features that are coupled to sub-shelf circulation and local basal melting.As the highresolution DEM record continues to grow, future analyses for PIG and other Antarctic ice shelves will provide new insight into the underlying processes controlling ice-ocean interaction, with implications for future ice sheet stability.

Summary/conclusions
We developed a method to correct and integrate high-resolution DEM observations with satellite/airborne altimetry data and           Equations 10 and 12, the cumulative basal melt rate along the A12 path is estimated.This procedure is repeated for "particle" B and all other "particles" in DEM1 that intersect DEM2.For the "along-flow distribution" approach, the cumulative basal melt rate for path A12 is assigned to each Eulerian grid cell along path A12, including grid cell C.This assignment is repeated for path B12 and all other paths for DEM1-DEM2 particles, so that many basal melt rate values will be assigned to grid cell C. The median basal melt rate is calculated from all paths intersecting C.This median value at C (and all other grid cells with nonzero path count) are used to populate the along-flow distribution basal melt rate map for DEM1-DEM2.This process is repeated for DEM1-DEM3 and all other valid downstream DEM1-DEMj combinations for the specified ~2-year time period.The same process is then repeated for all initial DEMi, and shelf-wide composites are generated as described in Section 2.9.3.1.For the "initial-pixel" approach, the cumulative basal melt rate for path A12 is assigned to cell A1.This process is repeated for the basal melt rate along path A13 and all other valid downstream DEMj to estimate "initial-pixel" stack median basal melt rate for A1, and all other pixels in DEM1.This "initial-pixel" stack median process is repeated for all valid DEMi, and these products are combined to create shelf-wide composites as described in Section 2.9.3.2.

(
Korona et al., 2009)) 40-m posting DEMs that covered some portion of the PIG shelf between January 5, 2008 and January 18, 2010.Unlike the sub-meter WorldView/GeoEye imagery, the ~5 m GSD SPOT-5 images are unable to resolve meter-scale ice sheet texture, and stereo image correlation often fails for flat, effectively featureless surfaces, leading to gaps in the output DEM.The km-scale ridges/troughs, ~100-1000 m wind-sculpted surface features, and rifts on the main PIG shelf, however, provide adequate texture for successful correlation.Compared to the WorldView DEMs, the SPIRIT DEMs include increased noise and additional artifacts, but cover a much larger area (~120 km swath width).Elevation values in the SPIRIT DEMs are represented as integers, with horizontal and vertical accuracy estimates of <10 m crossed the study area between 2009-2015, with data collection for each flight typically lasting <4 hours.The high-altitude LVIS surveys on October 20, 2009 and October 10, 2011 cover a significant portion of the main shelf, while other LVIS/ATM flights generally consist of a few sparse flightlines distributed across the shelf.

Figure S2 .
Figure S2.In addition to the filtered airborne data, a large sample of near-contemporaneous ICESat-1 GLAS data were available for co-registration of the 2008-2009 SPIRIT DEMs.After co-registration we estimate that the lower-resolution SPIRIT DEM products have 3-4 m or better absolute vertical accuracy (1-sigma).One of the DEMs (January 3, 2009) had large residual offsets between control point and DEM elevation, and we performed a secondary round of vertical bias correction (-3.1 m) to minimize offsets between this DEM and a 2010-2015 WorldView/GeoEye DEM per-pixel median elevation composite over flat, smooth surfaces near the main shelf.

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.snow) can potentially lead to surface elevation changes of ~1 m, and advection of small-scale surface features can lead to horizontal co-registration errors.We used a network of five 2012-2014 GPS sites on the outer shelf (Shean et al., 2017) as independent check points for WorldView DEMs.Corrected DEM elevations show good agreement (~0.72 m root mean squared error [RMSE] and ~0.57normalized median absolute deviation [NMAD]) with cm-accuracy surface elevations derived from GPS interferometric reflectometry (GPS-IR) antenna height records at each site.Unfortunately, no valid SPIRIT DEM pixels were available near the 2008-2010 GPS sites.

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.density for sea water ( V ) and ice ( @ ).Equation 8 can then be substituted into Equation 7 to obtain a mass-conservation expression for Lagrangian surface elevation change for the floating ice, dropping the constant d from the material derivative term: This set of DEMs provides thousands of combinations for Lagrangian Dh/Dt basal melt rate computation, with the flexibility to vary the time interval Dt.Most of the PIG shelf DEM data were acquired seasonally from ~October-March, so we computed interannual Dh/Dt for time intervals of ~1 and ~2 years.Longer time intervals decrease spatial resolution, as the observed Dh/Dt values are integrated across a longer path, but they provide improved signal-to-noise ratio for Dh/Dt, and we use the ~2 year products for further analysis.To calculate basal melt rate, we track each pixel in an earlier DEM acquired at time ti (DEMi) to its corresponding downstream location where it intersects a later DEM acquired at time tj (DEMj).Since our velocity fields vary over time (Section 2.5), an The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.appropriate time step Dt for this tracking is automatically determined based on the grid cell size and maximum velocities (e.g., ~10-20 days for 256-m grid).For each time step n { | 0 <  ≤ /∆,  ∈ ℤ b }, all valid pixels ("particles") from DEMi are propagated along flow paths (Figure 6) computed from the time-variable velocity fields.This propagation yields updated DEMi particle positions at time ( @ + ∆).For those particles whose paths intersect DEMj, we calculate the observed Lagrangian elevation change rate as:

For
our chosen Dt of ~2 years, a total of 117 unique DEMi with initial ti timestamps spanning 2008-2013 and sufficient DEMj intersection area were available over the PIG shelf.Each DEMi formed ~2-40 valid DEMi-DEMj combinations, yielding a final set of >1000 independently generated DEMi-DEMj basal melt rate products.

Figure 7
Figure 7 shows long-term Eulerian elevation change (dh/dt) for the full study area.From 2003-2010, thinning rates <30 km upstream of the grounding line were ~5-10 m/yr, while those farther upstream over the catchment were ~1 m/yr.From 2010-2015, thinning rates near the grounding line decreased to ~0-1 m/yr, with increased thinning of ~1-2 m/yr over the catchment.
These features are related to dense series of arcuate surface crevasses (e.g., Scott et al., 2010) that display elevation change due to advection during the 2010-2015 period.Individual DEMs show elevation differences of ~0.5 m between these crevasse bands and inter-band surfaces.Over the PIG shelf, we observe 2010-2015 dh/dt signals with spatial scales of ~10-15 km that are unrelated to advection of km-scale surface features (Figure 7D).We observe ~1-2 m/yr thickening downstream of grounding line on the north side of the inner main shelf and ~1 m/yr thinning over the south side of the outer main shelf.The South PIG shelf shows <1 m/yr thinning from 2010-2015, with ~3 m/yr thinning over upstream ice within ~10 km of the grounding line.The North shelf shows little elevation change with <0.5-1 m/yr upstream thinning.

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.Our full-shelf mean basal melt rates for 2008-2015 (~82-93 Gt/yr) are less than past estimates of ~95+/-14(Depoorter et al., 2013) and 101+/-8 Gt/yr(Rignot et al., 2013) for the 2003-2008 period.While our estimates fall within reported uncertainty, it is possible that this difference could represent a true decrease in shelf-wide melt from2003-2008 to 2008-2015, which    would be consistent with estimates from oceanographic observations that indicate a decrease from ~100 Gt/yr to ~73 Gt/yr between 2007 and 2009(Dutrieux et al., 2014b).However, this apparent decrease may be at least partially attributable to methodological differences (e.g., shelf area, flux-gate placement).The previous studies also mixed observations from different time intervals during a highly dynamic period in PIG's recent history, with dh/dt from 2003-2008 ICESat-1 data, velocities from an InSAR mosaic with approximate time stamp of2007-2008 (Rignot et al., 2011), and average SMB from 1979-2010.
surface velocity data to estimate Eulerian dh/dt, Lagrangian Dh/Dt, and ice shelf basal melt rates.Mean 2008-2015 basal melt rates for the full PIG ice shelf were ~82-93 Gt/yr.Local basal melt rates were ~200-250 m/yr near the grounding line, ~10-30 m/yr over the outer main shelf, and ~0-10 m/yr over the North and South shelves, with notable exception of ~50-100 m/yr near the grounding line of a fast-flowing tributary on the South shelf.The Lagrangian Dh/Dt basal melt rates show excellent The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.agreementwith, and provide spatial/temporal context for, in situ basal melt rate observations.Basal melt rates vary substantially across ~km-scale ice shelf thickness variations, with greater melting observed in basal channels and over deep keels near the grounding line and relatively shallow keels over the outer shelf.The methods and general results presented here provide a foundation for further analysis of the detailed spatiotemporal evolution of basal melt rates and connections with ocean observations for the PIG shelf during the 2008-2015 period.

Figure 1 :
Figure 1: Cumulative and annual DEM mosaics for West Antarctica.A) Cumulative mosaic of ~3000 WorldView/GeoEye stereo DEMs from 2010-2015, overlaid on Bedmap2 shaded relief map.Black outline shows location of Figure 2. B) DEM mosaic adjusted to EGM2008 geoid and stretched to show surface elevation over floating ice shelves.C) Total count of DEMs for the 2010-2015 time period.Note dense coverage over PIG.D) Annual DEM mosaics with same color scale as in A. Total cumulative 2-m DEM coverage 920

Figure 2 :
Figure 2: Context for the PIG catchment: A) Surface elevation, WorldView/GeoEye DEM mosaic over Bedmap2 DEM.Black outline shows location of PIG ice shelf in panel D. B) Bed topography from anisotropic interpolation.Note deep channels beneath main 925

Figure 3 :
Figure 3: PIG WorldView/GeoEye DEM co-registration results.A-C) ICP translation vector components required to co-register each DEM with filtered altimetry data.D) Median DEM error (DEM -GCP) with 16-84% spread, before (red) and after (blue) coregistration (see Shean et al. (2016) for details).The 2011/2012 cross-track stereo DEMs display larger errors before co-registration.After co-registration, bias is removed and residual error spread is typically <0.5-1 m, as summarized in Table2.935

Figure 4 :
Figure 4: Statistics for 2010-2015 WorldView/GeoEye DEMs and available 2009-2015 ATM/LVIS altimetry data over the PIG study area.Top row shows maps of elevation standard deviation, second row shows linear trend, and third row shows standard deviation of residuals from linear regression at every pixel.Left column ("Original") shows values for original DEM products before 940

Figure 5 :
Figure 5: Annual DEM composites using all available elevation data.Primary DEM sources are SPIRIT (top row), and WorldView/GeoEye (middle and bottom rows).

Figure 6 :
Figure 6: Illustration of Lagrangian Dh/Dt calculation and basal melt rate distribution on a Eulerian grid (light gray).Three DEMs (medium gray) acquired at times t1, t2, and t3 are resampled on this grid, with the same "features" A and B indicated as colored pixels.The position history for "particle" A is estimated using the velocity products described in Section 2.5, with paths indicated by dotted lines.Lagrangian Dh/Dt for A is calculated as (hA2-hA1)/(t2-t1).At each timestep along the path from A1 to A2 (A12), we estimate h (from observed Dh/Dt), velocity divergence (from observed velocity time series), and the local flux divergence.Using The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-209Manuscript under review for journal The Cryosphere Discussion started: 16 October 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 9 :
Figure 9: Relationship between km-scale ridge/trough features and initial-pixel basal melt rates for main shelf.Top row shows 2008-980