Interactive comment on “ In-situ GPS records of surface mass balance and ocean-induced basal melt for Pine Island Glacier , Antarctica ”

This paper describes a GPS-based technique to make estimates of basal mass balance (BMB) and surface mass balance estimates for Pine Island Glacier ice shelf (PIGIS). The inclusion of GPS-derived snow surface height relative to the GPS antenna is a particularly interesting addition to standard ice-shelf GPS analyses, and the results are tied in well to other utilized datasets including the stereo-imagery of the DEM and the surface snow/firn state from a RACMO-based firn density model (FDM). I have the advantage of being late in my review, so I don’t need to make the detailed


Introduction
The widespread availability of precise Global Positioning System (GPS) measurements has revolutionized the study of ice dynamics and glacier mass balance (e.g., Gao and Liu, 2001).Continuously operating dual-frequency GPS receivers provide high-frequency (1 Hz or less), highly accurate (< 1-3 cm) measurements of position, which can be used to derive surface velocity and elevation change.For applications involving ice dynamics, these measurements 30 offer important constraints for the mass continuity equation, which equates surface elevation change with surface mass balance, basal mass balance, and ice flux divergence.
Surface mass balance (SMB) processes include precipitation, sublimation, wind redistribution of surface snow, and melt water runoff.Regional climate models forced by reanalysis output (e.g., RACMO, MAR) now provide daily to The Cryosphere Discuss., doi :10.5194/tc-2016-288, 2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.
monthly estimates of Antarctic SMB on a relatively coarse grid (~5.5 to 27 km).In-situ SMB measurements are, however, still essential for model calibration and validation.Traditionally, SMB can be measured using stake networks, automated weather stations, near-surface radar surveys, and firn/ice cores, all of which require substantial field operations in remote locations.These measurements also tend to bias model calibration towards locations that are either accessible or of high scientific interest for other reasons.Reanalysis of glacier mass balance records indicate 5 that these measurement biases can significantly affect results, often resulting in overestimates of cumulative balance, as dynamic areas are less well sampled (Andreassen et al., 2016).Antarctic firn/ice core records indicate that SMB variability over most of Antarctica during the last 800 years was statistically insignificant, but that accumulation increased more than 10% for high-accumulation coastal regions (e.g., the Amundsen Sea Embayment) since the 1960s (Frezzotti et al., 2013).These areas are precisely those that that are most poorly sampled for SMB reanalysis 10 calibration and validation using traditional methods (i.e., firn/ice cores).
Accurate knowledge of firn compaction and its spatiotemporal variability is essential for interpreting observed surface elevation change in remote sensing data (e.g., satellite altimetry), and for partitioning this change into components related to ice dynamics and SMB (e.g., Shepherd et al., 2012;Wouters et al., 2015).Depth-dependent compaction rates can be estimated from a number of different methods, including vertical strain measurements (Arthern et al., 15 2010;Hamilton and Whillans, 1998), borehole optical stratigraphy (Hawley and Waddington, 2011), repeat phasesensitive radio-echo sounding (pRES) measurements (e.g., Jenkins et al., 2006) and ice-penetrating radar observations of internal layers over time (e.g., Medley et al., 2014Medley et al., , 2015)).In the absence of in-situ measurements, dynamic firn models forced by modelled SMB can provide estimates of compaction rates throughout the firn column, which can be integrated to obtain estimates for surface elevation change over time (e.g., Ligtenberg et al., 2011).20 Basal mass balance (BMB) beneath ice shelves (i.e., bottom melting, accretion) is driven by complex ice-ocean interaction.State-of-the-art ice-shelf cavity ocean circulation models offer some insight into sub-shelf ice-ocean interaction, but these models lack validation, as in-situ hydrographic observations are limited, especially within the sub-shelf cavity and the ice-ocean boundary layer.Some direct measurements are available from autonomous submersibles (e.g., Dutrieux et al., 2014) and instrumentation packages deployed through ice-shelf boreholes (e.g., 25 Stanton et al., 2013), but available data are limited to short time periods and small spatial extents.Precise measurements of surface elevation change from remote sensing observations (e.g., laser altimetry, digital elevation models (DEMs)) can also be used to infer BMB (e.g., Dutrieux et al., 2013;Pritchard et al., 2012;Shean, 2016), but time intervals between repeat observations are typically several months to years.
Here, we use temporally dense in-situ GPS records from the Pine Island Glacier ice shelf to constrain local SMB, firn 30 compaction, flux divergence, and BMB.We use changes in observed GPS surface elevation and antenna height to validate SMB and firn model output, and then use firn model results to estimate the time-variable downward velocity of the GPS poles.Flux divergence is estimated from observed strain rates between GPS stations.These constraints are then used to isolate elevation change due to local BMB.This approach yields temporally dense estimates of basal melt rates at spatially sparse GPS locations, which are combined with high-resolution DEMs from the same time period to 35 provide spatial context for the observed elevation change.These complementary results for the PIG ice shelf provide new information about the time-variable magnitude and spatial distribution of basal melting, offering indirect The Cryosphere Discuss., doi :10.5194/tc-2016-288, 2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.
observations of ice-ocean interaction and BMB sensitivity to ocean heat content variability, with implications for other rapidly evolving "warm-cavity" Antarctic ice shelves.
The main shelf has complex surface topography, including km-scale surface ridges and troughs that correspond to basal keels and channels, respectively (Bindschadler et al., 2011;Vaughan et al., 2012).A series of longitudinal (along-flow) ridges and troughs are present along the shelf centerline, with transverse (across-flow) ridges and troughs along the lateral margins (Figure 1).Local basal melt rates vary considerably across these features (Dutrieux et al., 20 2013;Shean, 2016).
Hydrographic observations seaward of the PIG calving front in Pine Island Bay suggest that basal melting declined by ~50% between 2010 and 2012 (Dutrieux et al., 2014).Long-term 2009-2015 mooring records seaward of the southern calving front (Figure 1) show a significant decrease in ocean temperature (~1-1.5°C)over ~450-770 m depths from late 2011 to early 2012, and then again from mid-2012 to early 2013 (Christianson et al., 2016).These 25 observations show that the ocean heat content at the PIG ice-shelf front varies considerably over monthly to interannual timescales.

PIG GPS sites
Several long-term GPS stations were installed on the PIG shelf as part of a larger investigation involving ice-sheet and ice-shelf dynamics (Bindschadler et al., 2011;Stanton et al., 2013).During the early part of this effort, two GPS 30 stations continuously collected data from January 2008 to January 2010 on the southern PIG ice shelf (PIG2) and the fast-flowing, grounded ice upstream of the grounding line (PIG1) (Figure 1).In addition, a ~2x2 km array of five stations (SOW1-4, BOAR, Figure 2 High-resolution optical imagery and DEM data (section 3.6) over the 2012-2014 sites show that SOW1, BOAR, and SOW3 were oriented approximately along a flowline within a longitudinal surface trough (Figure 2) that overlies a 5 longitudinal basal channel.Transverse ice-penetrating radar profiles across this channel provide ice thickness estimates of ~450-460 m near the apex and ~540 m over adjacent keels (Stanton et al., 2013).Figure 2 shows estimated hydrostatic ice thickness for longitudinal and transverse profiles across the GPS array.
A borehole was drilled through the ice shelf, and an instrument package with an upward-facing ice-bottom altimeter (acoustic ranger) was deployed beneath the shelf from January to February 2012.Measurements from this bottom 10 altimeter and complementary pRES experiments provided basal melt rate estimates of ~14-25 m/yr within the longitudinal channel (Christianson et al., 2016;Stanton et al., 2013).
The 2012-2014 GPS array was located near several transverse surface depressions (Figure 2).Local surface slopes were ~0.6-0.9°within the largest of these depressions, immediately downstream of SOW3 and SOW4.A notable linear depression located approximately 1 km upstream of SOW1 (arrow in Figure 2) opened as a rift in ~2014 (R1 in 15 Jeong et al. (2016)), and was subsequently the site of a large iceberg calving event that occurred around July 2015.
The fortuitous placement of the 2012-2014 GPS array near these features complicates interpretation of GPS records, but also provides new constraints on the spatiotemporal evolution of longitudinal strain and rift formation for the PIG shelf.

GPS position/velocity processing
As described in Christianson et al. (2016), GPS data were processed using differential-carrier-phase positioning relative to bedrock GPS sites (Backer Island [BACK,  Daily-averaged positions of these base stations were calculated using GAMIT and stabilized relative to a fixed circum-25 Antarctic reference frame using a Kalman filter (GLOBK, (Herring et al., 2015)).Geodetic GPS positions relative to WGS84 ellipsoid were calculated every 30 seconds; we analyzed a subset of these positions sampled at 10-minute intervals.All position data were converted to a local Cartesian horizontal coordinate system with orthometric heights above the EGM2008 geoid (Pavlis et al., 2012).Positions with uncertainty >8 cm were removed.The BOAR record was curtailed on April 29, 2013 (1.31 years duration), when an abrupt ~2.0 m elevation decrease and corresponding 30 horizontal offset occurred, suggesting that the pole fell over.
Horizontal velocities for each GPS station were computed from daily mean positions.Relative distances between stations were used to calculate strain, with linear fits to estimate strain rates.

GPS corrections
We estimated vertical tidal displacement for all GPS positions on the PIG ice shelf using the CATS2008A model (Padman et al., 2002).We used mean sea level pressure values from the 0.75°-grid-cell ERA-Interim reanalysis products (Dee et al., 2011) to correct for vertical displacement due to the inverse barometer effect (e.g., Padman et al., 2003).To do this, we removed the 2002-2016 median (985.21hPa) from 6-hour sea level pressure and scaled the 5 residuals by ~1 cm/hPa.
Figure 3 shows that tidal amplitudes for the GPS sites range from approximately -0.9 to +1.3 m and IBE amplitudes range from -0.3 to +0.3 m.These signals were removed from the GPS antenna phase center elevations, and residual high-frequency noise was removed with a low-pass filter (1.5-day cutoff), yielding smoothed time series for further analysis (Figure 3). 10 The constant 3.66 m pole length and 5.32 cm antenna phase center offset was removed from filtered absolute antenna elevation (h a ) to estimate corresponding absolute pole-base elevation (h pb ) relative to the EGM2008 geoid.Figure 4 shows a schematic of this geometry with additional observables as described below.

GPS antenna height and surface elevation
We computed mean daily antenna height above the surface (z a ) from L1 C/A code multipath surface reflections using 15 the GPS interferometric reflectometry methodology outlined in Larson et al. (2015) (Figure 4).This method takes advantage of the fact that the interference between the direct and reflected GPS signals produces characteristic frequencies in signal-to-noise ratio data recorded by the GPS receiver; these frequencies are directly related to the distance between the GPS antenna phase center and the reflecting surface.Geodetic antennas are designed to suppress multipath, so these interference patterns are best resolved at low GPS satellite elevation angles.20 Antenna height solutions were calculated for elevation angles of 5-25°, which sample the surface within a radial extent of ~5-50 m.Local surface slopes at each site are negligible, eliminating the need for an azimuthal correction (e.g., Larson and Nievinski, 2013).
Antenna height (z a ) was subtracted from filtered absolute antenna elevation (h a ) to obtain daily records of absolute surface elevation h (i.e., firn-air interface) for all sites.The absolute antenna elevation accuracy was ~1 cm and daily 25 antenna height precision was ~1 cm (Larson et al., 2015), with resulting surface elevation accuracy of ~1-2 cm.These surface elevation values are directly comparable with satellite/airborne laser altimetry data and stereo DEM products.
The multipath antenna height above surface (z a ) was also used to estimate pole-base depth below the surface (z pb ) at each site.
Continuous antenna height time series were generated for all seven PIG GPS sites.The SOW3 record was curtailed 30 on August 22, 2013, when antenna heights decreased below the ~0.5 m minimum threshold (Nievinski, 2013).
closest to the 2012-2014 GPS array (-75.07°N,-100.80°E, Figure 1).The values for adjacent grid cells are 0.74 mwe/yr near the grounding line of the main shelf (-75.15°N,-99.88°E) and 0.84 mwe/yr over the south shelf (-75.30°N,-101.14°E),providing some information on large-scale spatial variability.These values are consistent with SMB estimates of ~0.5-1.0 mwe/yr derived from CReSIS Snow Radar data collected upstream of the PIG grounding line (Medley et al., 2014(Medley et al., , 2015) ) and SMB estimates of 0.99 and 1.06 mwe/yr for stake measurements near 2006-2008 5 GPS sites over the upstream PIG trunk (Scott et al., 2009).
We also analysed 2011-2015 temperature data with 3-hour interval from the Evans Knoll (-74.85°N,-100.40°E, Figure 1) automated weather station (AWS) (Lazzara et al., 2012), located at an elevation of ~178 m (height above EGM2008 geoid) on a bedrock outcrop approximately 40 km north of the 2012-2014 GPS array (Figure 1).The AWS temperature values were scaled to the surface elevation of the GPS array (~66 m above EGM2008 geoid) assuming a 10 dry lapse rate of 9.8° C/km.
To provide historical context, we extracted 2-m air temperature from 1979-2015 for grid cells over the PIG shelf in 0.75°-resolution ERA-Interim reanalysis products (Dee et al., 2011) with 6-hour interval.The ERA-Interim temperature data display a bias of -2.81°C (median of offsets) compared to the scaled 2011-2015 AWS data, which is consistent with previous evaluations (Jones et al., 2016).This bias was removed from the ERA-Interim temperature 15 record, with residual offsets displaying root-mean-square error (RMSE) of 3.26 °C and normalized median absolute deviation (NMAD) of 2.77 °C.

Firn densification model
Model SMB output from RACMO2.3 was used to force the semi-empirical 1-D IMAU-FDM dynamic firn model (Ligtenberg et al., 2011) with 3-hour timesteps.The IMAU-FDM output is available at 2-day intervals.Velocities 20 (v ice ) across the firn-ice transition (defined as the layer with 917 kg/m 3 density) were assumed to be in equilibrium with average 1979-2015 SMB (a).Vertical velocity components for surface accumulation, surface sublimation, surface snow drift erosion/deposition, surface melt, dry firn compaction (v fc ), and a vertical mass difference buoyancy correction were computed for the 2008-2010 and 2012-2014 periods.These estimates were combined to provide time series of simulated surface (ℎ) and pole-base elevations (ℎ #$ ) for each GPS station.25

Basal melt rate
Mass conservation for a column with ice-equivalent thickness H (after removing a thickness correction d that accounts for total air content in the firn column) relates Eulerian thickness change (dH/dt) with dynamic thinning/thickening due to local flux divergence (∇ • , positive for extension), surface mass balance  (meters ice equivalent), and basal mass balance  (meters ice equivalent, defined as positive for melt): 30 This approach assumes that the total firn air content (d » 12 m for the PIG shelf (Shean, 2016)) remains constant for the period dt, which precludes the need to account for thickness change due to firn compaction, while allowing for The Cryosphere Discuss., doi :10.5194/tc-2016-288, 2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.
variable ice-equivalent SMB input .This assumption is supported by the limited variability (~1-3% or +/-0.3 m) in modelled IMAU-FDM total firn air content for the three PIG shelf grid cells during the relevant ~2-year study periods.
The material derivative definition relates Eulerian (fixed reference frame) and Lagrangian (reference frame moving with the ice column) thickness change: Rearranging Equation 2and substituting into Equation 1, we obtain the mass conservation equation for Lagrangian 5 thickness change: For a floating ice shelf in hydrostatic equilibrium, we can estimate ice-equivalent freeboard surface elevation ℎ A = (ℎ − ), where h is measured surface elevation and d is firn air content.We can then estimate total ice-equivalent thickness H: Assuming a constant density for ocean water ( E =1026 kg/m 3 ) and ice ( F = 917 kg/m 3 ), we can substitute Equation 410 into Equation 3, and rearrange to estimate basal melt rate from observed Lagrangian elevation change: Equation 5 depends on ice-equivalent surface elevation change.We now consider elevation change for the GPS pole base within the firn column of this simplified ice shelf.
The short-term effects of surface processes (e.g., accumulation, melting) are largely absent below the upper few meters of the firn column.Thus, the elevation change rate of the GPS pole base (Dh pb /Dt) is much more sensitive to 15 compaction within the underlying firn.The downward velocity of the pole base due to firn compaction (v fc ) varies as a function of pole-base depth z pb within the firn column.Values for v fc can be estimated by integrating firn model compaction rates from the firn-ice transition to ℎ #$ at each timestep.If SMB for the time period dt () is approximately equal to the long-term average SMB (a), then: which can be substituted into Equation 5to estimate basal melt rate from observed pole-base elevation change: 20 We neglect the slight reduction in total ice-equivalent thickness at the pole base (h pb ) vs. surface (h), as pole-base depth is negligible compared to total ice thickness ( #$ ≪ H).For the local flux divergence term, we use local freeboard thickness (as in Equation 7) or H from radar measurements (Stanton et al., 2013).Uncertainty is estimated as ~0.15 m/yr for downward firn-compaction velocity (Ligtenberg et al., 2011) If modelled pole-base velocities are correct, then basal melt rate estimates from Dh f /Dt (Equation 5) and Dh pb /Dt (Equation 7) should be similar.We also note that inferred basal melt rates are ~9 times more sensitive to pole-base elevation change (Dh pb /Dt) and local v fc than SMB () for a floating ice shelf.

Velocity 5
Figure 5A shows the velocities of the PIG1 and PIG2 stations.On the floating ice at PIG2, surface velocity increased from ~355 m/yr to ~380 m/yr between 2008 and 2010 as the GPS migrated downstream, with increased speedup beginning in late 2008.Velocities at PIG1 increased at a relatively steady rate from ~420 m/yr to ~460 m/yr as the station moved toward the fast-flowing PIG trunk (Figure 1).
In general, the stations displayed similar relative velocity evolution, with several abrupt >10-20 cm/day velocity changes (Figure 5).

Strain rate
For this analysis, we assume that SOW1, BOAR, and SOW3 GPS stations were oriented approximately along a 20 flowline, and that the observed displacements represent longitudinal strain.Figure 6A shows that the observed cumulative displacement between SOW1 and SOW3 (initial distance 2073.0 m) was ~6.7 m (~3.4 m/yr), with extensional strain rates of ~0.0017 yr -1 .Observed strain rates between SOW1 and BOAR (~2.1 m/yr over 1045.2 m, 0.0020 yr -1 ) were greater than those between BOAR and SOW3 (~1.5 m/yr over 1029.1 m, 0.0014 yr -1 ).Transverse strain rates were relatively low, with compression (-0.0004 yr -1 ) between BOAR and SOW2, and extension (0.0007 25 yr -1 ) between BOAR and SOW4.Subtle changes in strain rates between SOW1 and SOW3 were observed from 2012-2014 (Figure 6C).In general, increased (decreased) extensional strain rates occurred between SOW1 and SOW3 following an increase (decrease) in absolute GPS array velocity.These longitudinal strain rates correspond to shelf thinning rates of ~0.5-0.9 m/yr assuming no lateral spreading, with limited expected surface Dh/Dt (<0.07-0.13m/yr) for measured ice thickness of ~450-460 m (Stanton et al., 2013) 30 or for thickness estimated from local freeboard (~430-500 m).

Downslope flow
In addition to elevation change associated with strain rates and corresponding local flux divergence, some component of observed Dh/Dt may be related to deformation due to local gradients in the driving stress.Over grounded ice, this component of Dh/Dt is typically dominated by advection over basal topography, and the vertical component of the corresponding motion (V 0 ) can be estimated using observed horizontal GPS displacement and local surface gradients 5 from an independent DEM (e.g., Larson et al., 2015).
While this approach is irrelevant for a freely-floating ice shelf in hydrostatic equilibrium, we attempt to estimate an upper bound for the vertical component of Dh/Dt due to diffusion of local topography.To do this, we again consider observed relative horizontal displacements within the 2012-2014 GPS array.Local surface slopes near SOW1, SOW2, and BOAR are negligible (<0.2°), so we assume V 0 = 0 for these stations.If all of the observed ~3.4 m/yr relative 10 displacement between SOW1 and SOW3 was attributed to flow down ~0.6° local surface slopes at SOW3, then the associated V 0 magnitude would only be ~0.03m/yr, which is negligible compared to the observed 5.2 m/yr Dh/Dt.
Considering that much of the observed relative displacement appears to be related to longitudinal extension (Section 3.2), we neglect any potential Dh/Dt contribution related to V 0 .

GPS antenna and surface elevation change 15
Figure 7 shows antenna elevation (h a ) change relative to the initial absolute antenna elevation (h a0 ) at each station.These results show negative, highly linear (R 2 0.98-1.00)Dh a /Dt (Table 1), with rates of -1.6 to -2.1 m/yr at SOW1, SOW2, and BOAR, and higher rates of -5.2 m/yr and -3.8 m/yr at SOW3 and SOW4, respectively.The PIG1 rates over grounded ice are -7.6 m/yr with apparent concave-downward curvature.This is consistent with V 0 expected from local surface topography and dynamic thinning over the PIG trunk associated with velocity increases in 2006-2008 20 GPS observations (Scott et al., 2009) and satellite records (Joughin et al., 2010;Mouginot et al., 2014).
Figure 9C shows corresponding surface elevation (h) change relative to initial surface elevation (h 0 ) at each station.
The 2008-2010 surface elevation change at PIG2 is negligible (-0.13 m/yr).By contrast, surface elevations decreased significantly at all 2012-2014 sites, with rates of -0.9 to -1.3 m/yr for SOW1, SOW2 and BOAR, and rates of -4.1 m/yr and -3.0 m/yr at SOW3 and SOW4, respectively.25 Residuals about these linear fits (Figure 9D+E) are small (Dh a /Dt RMSE of 0.095 m and Dh/Dt RMSE of 0.143 m for sites on the PIG shelf), with some seasonal to annual variability observed at all sites.We also note apparent abrupt (~days-weeks) elevation changes that occurred across all stations in the 2012-2014 array (e.g., June 2012).

Antenna height evolution
Assuming that the pole base remains fixed within its original firn layer (effectively behaving as an isochron tracer in 30 the firn column), any observed decrease in antenna height above the surface (z a in Figure 4, Section 2.3) can be attributed to surface accumulation (e.g, snowfall, deposition of snow by wind).Conversely, an increase in antenna height can be attributed to surface ablation (e.g., melt, sublimation, removal of snow by wind) and compaction of snow/firn above the pole base.For convenience, we define the reflector height (z rh ) to track cumulative accumulation and ablation relative to the initial surface (h 0 ).35 At both PIG1 and PIG2, there are periods of relatively rapid reflector height increase (e.g., from May-August 2008), 5 followed by a steady decrease (e.g., August 2008-February 2009).This is consistent with periods of increased snow accumulation followed by several months of ongoing firn compaction with limited snowfall.The 2012-2014 records show similar periods of abrupt increase and steady decrease, with more limited duration.
All records show an abrupt reflector height decrease (~0.2-0.3 m) between December 2012 and January 2013, which is consistent with surface melting and/or enhanced firn compaction rates in the upper few meters of the firn column 10 (see Section 3.7).

High-resolution DEMs
In addition to the GPS elevation data, we generated WorldView/GeoEye stereo DEMs (Shean et al., 2016) with 32-m posting over the PIG shelf (Shean, 2016) to provide spatial context for the GPS time series.A total of 7 WorldView DEMs intersected the GPS positions between 2012-2014.We sampled DEM surface elevation at corresponding GPS 15 positions and compared with GPS-derived surface elevation (as described in Section 2.3) where possible.Table 2 shows statistics for the sampled DEM elevation compared with GPS surface elevation at each site (Figure 8).
In general, we observe good agreement between the two datasets, with RMSE of 0.72 m and NMAD of 0.57 m for the full sample (n=25).The DEMs display a slight bias (+0.43 m) relative to the GPS surface elevation.We note that the January 14, 2012 WorldView DEM is anomalously high, which biases DEM Dh/Dt trends with limited sample count 20 (Table 2).
High-resolution Lagrangian Dh/Dt elevation-change maps (see methodology in Shean, 2016) were computed for the 2012-2014 GPS sites by forward propagating 32-m DEM pixels from two initial DEM products (February 2, 2012 and October 23, 2012) using interpolated, time-variable surface velocity maps from Joughin et al. (2010) and Christianson et al. (2016).Lagrangian Dh/Dt maps were generated for all valid DEM combinations.Products with 25 time interval between 0.5-2.5 years were aggregated, and median Dh/Dt values were assigned to initial DEM pixel locations.
Figure 10 shows the composite Lagrangian DEM Dh/Dt maps, which provide spatial context for the GPS Dh/Dt estimates.Little or no elevation change was observed over longitudinal ridges, while areas within and near transverse depressions experienced enhanced thinning (Figure 10).This thinning was concentrated on the upstream side of the transverse depressions.The Dh/Dt products relative to the October 23, 2012 DEM (Figure 10D) also show the pattern 35 of thinning associated with the rift that opened upstream of SOW1 in ~2014 (Jeong et al., 2016).

Surface mass balance
The 1979-2015 RACMO average SMB over the central PIG shelf is ~0.9 mwe/yr.Monthly SMB climatology shows low accumulation rates of ~0.01-0.04mwe/month over the PIG shelf during the austral summer (November to February), and high accumulation rates of ~0.08-0.1 mwe/month during austral winter (March to October) (Figure 9F).Daily SMB products show periods of days to weeks with increased accumulation (e.g., March 2013) that can be 5 correlated with abrupt increases in surface reflector height.By contrast, the reflector height records show a steady decrease due to ongoing near-surface firn compaction during extended periods with little/no accumulation.
The ~3-4 week period between December 24, 2012 and January 17, 2013 was relatively warm, with scaled AWS atmospheric temperatures of ~1-5°C for most days (Figure 9F).Surface elevations decreased by ~0.2-0.3 m across the entire GPS array during this period (Figure 9A), which is consistent with significant surface melting.The daily 10 RACMO SMB data also show two accumulation events during the last week of December 2012 (Figure 9E), which likely involve rain on snow.Some component of the observed surface elevation decrease could be associated with enhanced near-surface melting and firn compaction during these events.No corresponding short-term changes were recorded by the pole-base GPS elevations during the ~3-4 week period, suggesting that the processes responsible for the observed surface changes were limited to the upper ~1.

Firn model results
Estimated downward pole-base velocities due to firn compaction (v fc ) were all ~0.70-0.75m/yr (Figure 7, Table 1) 20 from 2008-2010 and 2012-2014, despite the range of initial pole depths.A slight decrease in the compaction rate occurred over time, but the curves appear linear at these depths (Figure 7).These values are consistent with observed Dz rh /Dt and Dh/Dt during periods with little or no surface accumulation Figure 9C shows that the IMAU-FDM simulated surface elevation (ℎ) ranges from -0.1 to +0.4 m from 2008-2010 and -0.2 to +0.2 m from 2012-2014.The observed ℎ/ trend is +0.17 m/yr from 2008-2010, with no significant 25 trend from 2012-2014.The magnitude and timing of the modelled surface elevation variability is consistent with the detrended observed surface elevation change (Figure 9D).The observed Dh/Dt trends (-1 to -4 m/yr) (Figure 9C), however, cannot be explained by modelled elevation change due to only SMB and firn.

Basal melt rates
We computed basal melt rates from surface Dh/Dt and pole-base Dhpb/Dt using Equations 5 and 7, respectively.The 30 resulting values range from ~2-4 m/yr at PIG2 to ~40-43 m/yr at SOW3, with good agreement between the rates obtained by surface and pole-base elevation change (Table 1).
The 2012-2014 melt rate estimates show significant spatial variability.The three upstream stations (SOW1, SOW2 and BOAR) show lower melt rates of ~9-14 m/yr, while the downstream stations near the transverse depression (SOW3 and SOW4) have higher rates of ~29-43 m/yr for the same time period.35 The Cryosphere Discuss., doi:10.5194/tc-2016-288,2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.

Assumptions
The results in the Section 3 relied on several simplifying assumptions.We now offer further discussion of these assumptions and their potential influence on our results.

SMB spatial variability
We used modelled SMB from a single RACMO2.3grid cell to drive the IMAU-FDM dynamic firn model, and applied 5 the result to all GPS stations.We expect these parameters to vary spatially (e.g., Medley et al., 2015) due to local environmental conditions (e.g.PIG2 is >400 m higher than SOW1-4 stations on the shelf) and local surface topography (e.g., km-scale ridges/troughs), which will affect near-surface winds and snow redistribution.
In addition, the IMAU-FDM values do not account for horizontal advection of the 1-D firn column through spatiallyvariable RACMO fields (accumulation, surface temperature, etc.) over time.The GPS sites over the PIG shelf are 10 moving northwestward at ~4 km/yr (Figures 1 and 5), which is nearly double the observed PIG shelf velocities from the mid-1970s (Mouginot et al., 2014).Thus, the local firn columns beneath the GPS sites likely experienced variable SMB input over their ~50-100 km horizontal path during the corresponding 1979-2015 time period.This suggests that the true firn column thickness and compaction rates may differ from the IMAU-FDM estimates.For this reason, we use a firn air content estimate (d » 12 m) derived from ice-penetrating radar two-way travel time and altimetry surface 15 elevation measurements (see appendix A of Shean, 2016).
We also expect local SMB and firn column thickness variability within the array.The greater z rh surface reflector heights (a proxy for surface accumulation) at SOW3 indicate that enhanced local accumulation occurred within the transverse depression (Figure 2), likely due to wind-blown snow.However, we note that all of the remaining 2012-2014 stations display similar z rh reflector height evolution, as do the two 2008-2010 stations (Figure 9A).Future 20 RACMO and IMAU-FDM products with improved resolution should provide additional constraints on spatial variability.

Pole-base settling
Some of the reflector height increase and/or observed negative Dh/Dt could be related to settling of the poles within the firn.We assume that the poles froze in place shortly after installation, and the contact area (~1200 cm 2 for a ~1-25 meter-long cylinder with ~3.8 cm diameter) with surrounding snow/firn should be sufficient to counter the downward gravitational force.Thus, we expect that the Dh/Dt recorded by the GPS antenna pole represents rates at the base of the pole, rather than rates within an overlying layer.
A related consideration involves heating of the exposed pole during summer, which might lead to decoupling from the surrounding snow/firn and thus, additional penetration of the pole base within the firn.However, we do not see 30 any indication of such settling from December 2012 to January 2013, when surface elevations decreased by ~20-30 cm and pole base elevations showed little change.The lack of pole-base elevation change suggests that surface meltwater did not percolate below ~1-2 m depth into the snow/firn.

Strain rate length scales
As discussed in Section 3.2, the expected surface Dh/Dt from local flux divergence is <0.1 m/yr, assuming that the observed ~1.5 m/yr relative displacement is evenly distributed across the ~1 km distance between GPS stations (BOAR and SOW3).Even if this strain is concentrated over a shorter distance (e.g., ~200 m), this contribution only increases to <0.5 m/yr, which is still significantly less than the observed ~3-4 m/yr Dh/Dt signals.Interestingly, 5 despite similar extensional strain rate estimates for SOW1-BOAR (0.0020 yr -1 ) and BOAR-SOW3 (0.0014 yr -1 ), we note a large difference in local Dh/Dt values at these stations (Figures 9 and 10).This supports the assumption that local flux divergence for these sites is a minor component of the observed Dh/Dt, and consequently, that the observed Dh/Dt is primarily caused by basal melt.

Long-term SMB and firn compaction
The observed reflector height increase (z rh ) of ~0.7-1.0 m/yr at all GPS sites (Figure 9A) is consistent with observed long-term SMB estimates (a) and downward near-surface velocity estimates due to firn compaction (v fc ).If we assume a near-surface mean snow/firn density of 0.5 kg/m 3 , the ~0.7-0.9 mwe/yr RACMO SMB estimates correspond to a surface elevation increase ~1.4-1.8 m/yr.Removing the expected ~0.7-0.75 m/yr surface elevation decrease 15 associated with firn compaction, we arrive at expected relative reflector height increases of ~0.7-1.0 m/yr, which matches observed Dz rh /Dt values (Figure 9A+B).
The fact that we observe similar trends for surface Dh/Dt and pole-base Dh pb /Dt (Table 1) also supports the assumption

Residual Dh/Dt variability
The surface (Figure 9D) and pole-base (Figure 9E) residuals appear to be uncorrelated.This suggests that seasonal surface processes (e.g., accumulation influencing near-surface compaction rates) are not responsible for driving the observed pole-base variability.We considered several possible sources for the observed sub-annual variability in the surface and pole-base records, including ocean (e.g., currents, sea surface height), atmospheric (e.g, pressure, 30 temperature), and dynamic processes (e.g, resistive stress from sea ice and/or mélange in shear margins).
Unfortunately, we were unable to definitively determine the cause(s) for these variations in the ~2 year GPS records.
The Cryosphere Discuss., doi :10.5194/tc-2016-288, 2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.Some of the short-term (days-weeks) variability (e.g.June 2012) observed across all five 2012-2014 stations (Figure 9D+E) could be related to insufficient or incorrect IBE correction.The magnitude and timing of these systematic anomalies, however, suggests that they are likely related to grounding/ungrounding events (Joughin et al., 2016).

Strain rate history, rifting, and grounding evolution
The lateral shear across the GPS array is related to increased longitudinal extension closer to the PIG centerline, likely 5 due to locally enhanced ductile deformation (i.e., "necking" (Bassis and Ma, 2015)) across transverse depressions, and/or expansion of basal/surface crevasses and rifts.The SOW3 station, which lies within a large transverse depression (Figure 2), displays a slight acceleration in negative elevation change (Figure 9E), potentially due to increased local extension within the depression.
The fact that the rift immediately upstream of SOW1 ultimately became the site of the 2015 calving event (Jeong et 10 al., 2016) suggests that the velocities observed at the GPS array are not necessarily representative of the velocities near the PIG grounding line.An upstream regrounding event would slow ice upstream of the GPS array, initially resulting in increased extensional strain rates across the transverse rifts/depressions, followed by a velocity decrease at the GPS array.Conversely, an upstream ungrounding event would initially lead to decreased extensional strain rates across the transverse rifts/depressions, followed by an increase in observed GPS velocities . 15 We suggest that an upstream regrounding event in ~June 2012 (Joughin et al., 2016) could be responsible for increased strain rates across the GPS array (Figure 6).Similarly, an ungrounding event in ~April 2013 followed by a grounding event in ~November 2013 could explain the decrease and subsequent increase in strain rates.
There is an abrupt ~10-20 cm pole-base elevation decrease at both SOW3 and SOW4 in late 2013, near the end of the records (Figure 9E).No surface elevation information is available at SOW3 due to missing antenna height data, but a 20 corresponding surface elevation decrease is observed at SOW4 (Figure 9D).These elevation decreases do not appear to be related to site servicing.Rather, these observations are consistent with relatively abrupt local extension within the transverse depression at SOW3 and SOW4, but not the upstream sites.The timing of this event corresponds with observed lengthening of the large rift (R1) upstream of SOW1 (Jeong et al., 2016), supporting the hypothesis that relatively rapid, localized extension occurred across the transverse depressions and rifts during this period.25

Comparison with in-situ basal melt rate observations
The GPS basal melt rate estimates (~9-14 m/yr for SOW1, SOW2 and BOAR sites) appear consistent with those from bottom altimeter (~14.7 m/yr from January-February 2012) and pRES (~15-25 m/yr) measurements of Stanton et al. (2013).These measurements provide some validation for the GPS results, but a direct comparison may be imprudent.
The Stanton et al. (2013) measurements were taken at a site approximately 1.34 km upstream of SOW1 (K.Riverman, 30 personal communication, 2016), within the same longitudinal channel, but much closer to the site of the R1 rift (Figure 2A).Given the complex and evolving local surface topography, we also expect local variation in basal topography and associated melt rates during the 2012-2014 period.Furthermore, the altimeter sampled a ~5 cm diameter spot with unknown upstream/downstream orientation, approximately 30-40 cm from the edge of the 20 cm borehole.Aside from local melt variability expected due to turbulent flow near the altimeter pole or borehole edge, the altimeter 35 The Cryosphere Discuss., doi :10.5194/tc-2016-288, 2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.provided a relatively small spatial sample compared to the GPS results, which are sensitive to changes in a column of ice with much larger spatial extent (10s to 100s of meters).

Basal melt rate spatial variability
The GPS records at SOW1, SOW2 and BOAR show similar Dh/Dt rates and residuals, which is consistent with their apparent orientation on the same "block" between transverse rifts/depressions.They are also located over the same 5 set of longitudinal basal channels, and should be exposed to similar sub-shelf circulation.
The DEM Dh/Dt maps show enhanced surface elevation change rates, and thus higher basal melt rates, on the north side (left side in Figure 10) of longitudinal basal channels and on the upstream side of transverse depressions.The latter is consistent with increased Dh/Dt observed at the SOW3 and SOW4 sites, located on the upstream side of a transverse depression.10 This relationship is potentially related to enhanced buoyant flow over increased basal slopes (e.g., Jenkins, 2011) and/or turbulence as water within the upstream meltwater channel first enters the "cavern" of the transverse depression.
We also suggest that the transverse depressions may serve as conduits connecting flow between the longitudinal channels, potentially leading to increased circulation velocity and higher melt rates within the transverse depressions.

Implications for melt rate sensitivity to ocean variability 15
The 2012-2014 GPS data reveal subtle (~2-4%) velocity changes that appear correlated with observed variations in ocean temperature records from moorings (Figure 1) in Pine Island Bay (Christianson et al., 2016).Our analysis supports the alternative Christianson et al. (2016) hypothesis, that these velocity variations are primarily related to upstream grounding evolution (Joughin et al., 2016), and extension across a series of transverse depressions (Figure 6).20 Both the GPS pole-base Dh pb /Dt and surface Dh/Dt fits are highly linear, with no significant variation in inferred basal melt rates for this 2-year time period.If sub-shelf melt rates beneath the GPS array were directly related to ocean heat content beyond the shelf front in Pine Island Bay, and had decreased by ~50% as suggested by Dutrieux et al. (2014), a significant change in observed Dh/Dt would be expected during this period.The lack of any significant Dh/Dt deviation suggests that melt rates at these sites are not significantly affected by observed ocean temperature variability.25 This finding suggests that either: 1) these sites are not representative of melt rates for the greater shelf (e.g., those near the grounding line), 2) the oceanographic measurements near the PIG ice front (Figure 1) are not representative of water circulating beneath these ice-shelf sites, or 3) local melt rates are less sensitive to the observed oceanographic changes than previously assumed.

Future work 30
High-resolution velocity maps derived from sub-meter imagery could potentially constrain local velocity divergence and length scales for observed strain between GPS receivers.In addition, seismic data from stations deployed near the GPS array and regional sites could help constrain the timing and location of rift propagation and grounding/ungrounding events.
It may be possible to further constrain firn-compaction rates, and thus long-term SMB, using relative layer thicknesses observed in CReSIS snow radar measurements (e.g., Medley et al., 2015) or in-situ pRES observations (e.g, Jenkins et al., 2006).However, airborne radar data over the PIG shelf suffer from clutter due to km-scale topography and crevasses, while the available intermittent pRES records (Stanton et al., 2013)  The GPS antenna height records document a relative surface increase of ~0.7-1.0 m/yr, which is consistent with modelled SMB of ~0.7-0.9 mwe/yr and downward firn compaction velocities of ~0.7-0.75 m/yr.An abrupt ~0.2-0.3 m surface elevation decrease, likely due to surface melt, is observed during a period of warmer atmospheric 15 temperatures from December 2012 to January 2013.Observed longitudinal strain rates for the 2012-2014 GPS array are ~0.001-0.002yr -1 , with negligible associated surface elevation change.
Observed surface Dh/Dt is highly linear for the PIG shelf sites, with trends of -1 to -4 m/yr and residual RMSE of ~0.1 m for the 2012-2014 sites.Similar estimates with reduced residual variability are obtained after removing simulated downward GPS pole-base velocity due to firn compaction from GPS antenna elevation records.Estimated 20 basal melt rates are ~10 to ~40 m/yr near the center of the fast-flowing PIG shelf, and ~2-4 m/yr for the southern shelf.These melt rates are similar to those derived from complementary in-situ instrument records (Stanton et al., 2013) and high-resolution stereo DEMs (Shean, 2016).
Both GPS and DEM records show higher melt rates within and near transverse surface depressions and rifts associated with longitudinal extension.Basal melt rates for the 2012-2014 period show limited temporal variability, despite 25 significant changes in ocean heat content at the ice front and likely in the ice-shelf cavity.Residual elevation change variability is likely related to upstream grounding/ungrounding events and variable strain rates across transverse depressions/rifts.

Site Time period
Days h a0 (m)              1).E) Detrended GPS pole-base elevation (h pb ) after removing v fc (see Figure 7), with arbitrary y-axis offset.Legend lists linear fit, which can be directly compared with Dh/Dt fits in D (see Table 1).Note reduced residual magnitude and dampened ~seasonal signal compared to D. Unlike surface records, no significant change is observed from Dec. 2012 to Jan.
, ~0.1 m/yr for elevation change associated with local flux divergence, and ~5 kg/m 3 for the density of ice and ocean water.25 The Cryosphere Discuss., doi:10.5194/tc-2016-288,2017 Manuscript under review for journal The Cryosphere Published: 3 January 2017 c Author(s) 2017.CC-BY 3.0 License.
1979-2015 ERA-Interim 2-m air temperatures over the PIG shelf shows many previous warm periods with greater magnitude and duration than the 2012/2013 period.However, the 2012/2013 warm period stands out in the past decade, with the most recent comparable periods during 2005/2006 and 2006/2007.
that 2008-2010 and 2012-2014 SMB () is consistent with long-term 1979-2015 SMB (a) and firn-compaction rates.The limited variability in surface elevation at PIG2 (Figure 9C) suggests that the observed 2008-2010 SMB over the 20 South PIG shelf was approximately in balance with ongoing firn compaction and basal melt.The 2012-2014 sites, however, show significant surface elevation decrease that greatly exceeds simulated IMAU-FDM surface elevation (ℎ) variability (+/-0.2m) due to SMB and firn compaction.No significant Dℎ/Dt trend is expected during this period, suggesting that SMB and firn compaction are consistent with long-term values, and the large observed Dh/Dt values must be attributed to basal melt.25 likely lack the sensitivity to detect small changes in firn layer thickness during the ~3-week observation period.These limitations highlight the value of long-5 term GPS records to constrain surface evolution where observations are sparse and model results are poorly constrained.Expanding the scope of our study to include the full archive of GPS data for the Antarctic and Greenland ice sheets would offer a valuable new dataset for ice-sheet SMB and model validation.6 Summary and conclusions We analyzed 2008-2010 and 2012-2014 GPS records for PIG.These data provide validation for surface elevation 10 and Dh/Dt derived from high-resolution WorldView stereo DEM records, with sampled DEM RMSE of ~0.72 m and NMAD of ~0.57m.
, reflector height change Dz rh /Dt, antenna elevation change Dh a /Dt (equal to pole-base elevation change Dh pb /Dt), modelled downward pole-base vertical velocity due to firn compaction v fc , surface elevation change Dh/Dt, and 5 corrected pole-base Dh pb /Dt after removing v fc .Corresponding ice-equivalent basal melt rates  computed for both the surface Dh/Dt and corrected pole-base Dh pb /Dt (equations 5 and 7, respectively).* Note: PIG1 values over grounded ice do not include correction to remove expected Dh/Dt due to advection along local surface slope (V 0 ).

Figure 3 :
Figure 3: A) Original absolute GPS antenna elevation (light gray), after tide correction (mid-gray), and after tide+IBE correction (black) for SOW4.Red line shows smoothed time series and yellow dashed line is linear fit (-3.76 m/yr).Sampled DEM elevations (cyan) show surface elevation, which is offset from GPS antenna elevation by antenna height (see Figures 4 and 8).B) High-frequency (<1.5 days) component of GPS record and CATS2008A tide model prediction, showing excellent

Figure 4 :
Figure 4: Schematic of GPS station geometry.Absolute surface elevation (h, dark blue line) is computed by removing multipath antenna height above surface (z a ) from absolute antenna elevation (h a , black line).With known pole length, we can calculate pole-base depth below the surface (z pb ) and absolute pole-base elevation (h pb , red line).At time t 1 (right panel), ongoing firn compaction resulted in decreased antenna and pole-base elevation, while new snow accumulation resulted in

Figure 5 :
Figure 5: Station velocities derived from daily mean positions for A) 2008-2010 GPS sites, and B) 2012-2014 GPS sites.Note offset between SOW2 and SOW4, indicative of lateral shear across the ~2 km wide array, with greater extension near the center of the PIG shelf.

Figure 6 : 5 The
Figure 6: A) Observed cumulative displacement between SOW1 and other 2012-2014 stations.Legend lists initial distances.B) Observed cumulative strain, with best fit strain rate listed in legend.C) Smoothed residuals from linear fit, highlighting subtle variations and inflections in strain rates.D) Daily GPS velocity.Note timing of abrupt absolute velocity changes observed at all sites and inflection points in observed strain across the array.

Figure 7 : 5 Figure 8 :
Figure 7: Filtered (solid line) and original (transparent) GPS antenna elevation (h a ) relative to initial absolute antenna elevation (first item in legend).Legend includes linear Dh a /Dt fit to filtered antenna elevations.Dashed black lines show IMAU-FDM estimated downward velocity due to firn compaction (v fc ) at the pole base.

Figure 9 :
Figure 9: Plots of observed and modelled values for GPS sites (see Figure 4 for reference).A) Reflector height (z rh ) relative to initial surface (h 0 ).Positive values indicate surface height increase relative to GPS antenna.Legend indicates initial and monthly RACMO2.3SMB.G) Observed AWS and ERA-Interim 2-m T records for PIG shelf, with above-zero AWS T plotted in red.Note period of extended warm T in Dec. 2012 to Jan. 2013, which corresponds to ~0.2-0.3 m surface elevation decrease.

Figure 10 :
Figure 10: WorldView DEMs and composite median Lagrangian Dh/Dt products generated using A-B) initial DEM from