Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading

The full-text may be used and/or reproduced, and given to third parties in any format or medium, without prior permission or charge, for personal research or study, educational, or not-for-pro t purposes provided that: • a full bibliographic reference is made to the original source • a link is made to the metadata record in DRO • the full-text is not changed in any way The full-text must not be sold in any format or medium without the formal permission of the copyright holders. Please consult the full DRO policy for further details.


Introduction
Rapid changes in climate in the Antarctic Peninsula (AP) over the past 50 years have led to the retreat and eventual collapse of several major ice shelves (Fig. 1), such as Prince Gustav by 1993-1995(Rott et al., 1996, Larsen A in 1995 (Rott et al., 1996), and Larsen B in 2002 (Rack and Rott, 2004) (see Cook and Vaughan (2010) for a complete summary). In response to ice shelf collapse, tributary glaciers have exhibited acceleration and thinning (e.g., De Angelis and Skvarca, 2003;Rignot et al., 2004;Scambos et al., 2004) and this dynamic ice loss induces a solid Earth response which may be observed in GPS records.
The study of Thomas et al. (2011) identified markedly-increased uplift in GPS coordinate time series from the Northern Antarctic Peninsula (NAP) that they associated with ice unloading related to the breakup of Larsen B Ice Shelf in 2002. This uplift was best captured in the near-continuous cGPS record at Palmer station which exhibited an increase in uplift rate from 0.1 mm/yr prior to 2002.2, to 8.8 mm/yr thereafter. Thomas et al. (2011) suggested that the effect was due to the elastic response of the solid Earth but they were not able to satisfactorily reproduce the increased uplift rates with an elastic model, which they suggested was at least partly due to the weakly defined magnitude and spatial pattern of icemass loss in their model.
The NAP lies in a complex tectonic setting which passes from active subduction along the South Shetland Trench, located north   of the South Shetland Islands, to passive margin west of 65 • W at the intersection of the Hero Fracture Zone with the Shetland Platform (see Fig. S1 in the supplementary material). The Bransfield Basin separates the South Shetland Islands from the NAP and has active volcanism along a mid-axial region, suggesting a back arc tectonic setting (Barker and Austin, 1998;Barker et al., 1991). The conversion from active subduction to passive margin, and hence mantle conditions more likely reflective of low viscosity, took place relatively recently at ∼ 4 Ma before present along the AP margin just south of Hero Fracture Zone . One of the GPS stations used in this study (ROBI, see Section 2.1) lies along a presumed incipient rift axis as expressed by the active volcanic chain of the Seal Nunataks (González-Ferrán, 1983). However, there is little known in relation to the mantle or crustal configuration beneath the Seal Nunataks region. The reader is referred to Barker (1982) and Larter and Barker (1991) for a full tectonic history of the region. Due to the active tectonic setting of the region, the mantle is likely to have a relatively low viscosity compared with other locations undergoing deformation in response to changes in ice-mass e.g. East Antarctica or Fennoscandia. Using a combination of inferred ice history, GPS and GRACE data, Ivins et al. (2011) suggested this region has a relatively thin lithosphere (20-45 km) and a low viscosity mantle (3-10 × 10 19 Pa s). Due to the low viscosity nature of the upper mantle, the Earth's viscous response to ice-mass change in the AP is much more rapid than in other regions of Antarctica, and post-1995 unloading events may hence be contributing to the observed uplift in the NAP through viscoelastic rebound. Likewise, assuming a Maxwell rheology, there may be very little, to no, residual response of the NAP to unloading events associated with recession from the Last Glacial Maximum. In this study we use cGPS data from the NAP to constrain a local model of solid Earth response to a high resolution present-day mass loss field (Scambos et al., in review). By comparing the modelled elastic uplift time series with the observed cGPS uplift time series from Palmer Station we show that elastic uplift alone cannot reproduce the cGPS observations. A viscoelastic model is then used to predict uplift based on a range of Earth models, and results are compared with the Palmer record to obtain the range of best fitting models. Finally, the Earth model is further refined using six shorter cGPS time series from the NAP (see Table 1). Fig. 1 shows the locations of available cGPS sites in the NAP. Of these, we used the seven sites closest to the region of ice-mass change (see Table 1). Palmer is a long term station with an excess of 15 years of data, and the remaining six sites were installed in 2009-2010 as part of the LARISSA project (LARsen Ice Shelf System, Antarctica, 2013) (http://www.hamilton.edu/expeditions/ larissa). We did not use the record from O'Higgins (a compilation of three records from two adjacent stations, OHIG, OHI2, OHI3; labelled OHI2 in Fig. 1) as a constraint as it lies far from the region of largest mass loss and as such may be affected by potential lateral heterogeneity in Earth structure. Spring Point (Bevis et al., 2009) (SPPT in Fig. 1) was also excluded due to the lack of data at this site; however we compare our results with both of these records in the supplementary material.

GPS
The cGPS data from 1998 through to the end of 2012 were processed using a Precise Point Positioning strategy (Zumberge et al., 1997) within GIPSY-OASIS v6.1.2. Homogeneously reprocessed (as of 2012) satellite clocks and fiducial-free orbits were fixed to values provided by the Jet Propulsion Laboratory. The receiver clocks, tropospheric zenith wet delay, tropospheric gradients and station coordinates were estimated in standard ways (e.g. Thomas et al., 2011). For the troposphere, we adopted the Vienna Mapping Function v1 (Boehm et al., 2006) and we assumed hydrostatic zenith delays based on the European Centre for Medium-Range Weather Forecasts. Higher order ionospheric effects (Petrie et al., 2011) were not corrected. Solid Earth tides were corrected according to the IERS 2010 conventions (Petit and Luzum, 2010) and ocean tide loading displacements were corrected based on the TPXO7.2 ocean tide model (Egbert et al., 2009), which has been shown to perform very well in this region , using the SPOTL software (Agnew, 1996(Agnew, , 1997 in a centre of mass of the whole Earth system (CM) reference frame (Fu et al., 2012). We fixed carrier phase ambiguities to integers where possible (Bertiger et al., 2010). Observations were weighted according to their elevationdependent scatter from a preliminary, but otherwise identical, solution.
Fiducial-free daily site coordinates were transformed into the ITRF2008 (Altamimi et al., 2011) using a 6-parameter transformation and then corrected for changes in atmospheric pressure loading using global 2.5 × 2.5 degree grids of loading, in a centre of figure of the Earth (CF) reference frame, derived from the National Center for Environmental Protection (NCEP) reanalysis surface pressure dataset (van Dam, 2010). The correction was applied by removing a daily average displacement with respect to a mean reference value. Large outliers from the DUPT time series were manually identified and removed, and then a median filter was applied to all time series. We only consider the height time series in this paper.
Velocities and realistic uncertainties were estimated using the CATS software (Williams, 2008), along with annual and semiannual harmonics. cGPS time series contain time-correlated noise which inflates the true velocity uncertainties relative to the formal errors obtained from a conventional linear regression. We consider a, now common, white noise plus flicker noise model using the CATS software from which we determine velocity uncertainties. We scale these to 2-sigma for subsequent use. Velocities and 2-sigma uncertainties for each cGPS site are given in Table 1. Below we compare model output both to the height time series and to vertical velocities derived from the time series. For consistency, all model-predicted uplift rates were estimated over the same time period as the cGPS observed uplift rate.

Ice-mass loss
The input ice load model is based on a dataset of elevation change derived from Digital Elevation Model (DEM) differencing and ICESat data covering the NAP region north of 66 • S. The time span of the data varies for different sub-regions. For the Larsen B tributary glaciers data are available for two time periods, 2001(Shuman et al., 2011(Berthier et al., 2012. Comparing the two time-periods reveals differences in spatial patterns of elevation change (Figs. 1b, c) but the overall estimated mass loss during these two periods differs by less than 10% (Berthier et al., 2012). For the areas north and west of Larsen B, the data span the period 2003. In all cases we take the rate of ice-mass change as being constant throughout the respective data periods; we discuss extrapolations to other time periods later.
The data were converted to a set of 17 846 loading discs with areas between 0.9 and 1.1 km 2 for input to the model, where the height of each disc represents a mass loss or gain, using an ice density of 900 kg/m 3 to convert to equivalent water height following Berthier et al. (2012). Data gaps over large glaciers were infilled using an inverse distance weighting algorithm (inpaint_nans within Matlab). Discs with very small mass change in the range ±0.5 m weq /yr have a negligible effect on deformation at sites tens to hundreds of km distance from the source of loading and were discounted from the ice load model to speed computation time. This was tested using the best fitting Earth model and resulted in no more than ±0.2 mm/yr differences in uplift rates at the cGPS sites. The resulting ice-mass change model is shown in Fig. 1a with the two periods of mass change for Larsen B glaciers, 2001-2006 and 2006-2011, shown separately in Figs. 1b and 1c, respectively. The uncertainty on the elevation change dataset is ±1 m/yr (2-sigma) (Scambos et al., in review), and this error bound was used to create upper and lower limits on our input ice load model.
The pattern of mass loss for glaciers feeding the former Prince Gustav and Larsen A ice shelves prior to 2001 is not known. However, there is evidence that these glaciers reached their current velocities within a few years of the breakup in 1995 (Rott et al., 2008), and this velocity has been maintained 15 years after ice shelf collapse (Rott et al., 2011). Observations of the glaciers feeding the more southern former Wordie Ice Shelf (Wendt et al., 2010), which disintegrated in a series of events between 1966 and 1989, also suggest that high rates of mass loss are sustained over decades. We therefore assume that the observed mass changes for these northerly regions are representative of ice-mass loss from 1995 to the present-day. We do not model any ice history before 1995, but instead take account of any ongoing deformation related to ice-mass changes before this time by estimating a linear background rate from the Palmer cGPS observations (see Section 3.1). In detail, we assumed that mass loss in a region began at the half year mark after the collapse of the corresponding ice shelf (i.e. 1995.5 for Prince Gustav and Larsen A glaciers, and 2002.5 for Larsen B glaciers), and continued at the same rate to present-day (2013.0). For the Larsen B glaciers this is clearly justified, as Scambos et al. (2004) show glacier acceleration and thinning commenced a few months after ice shelf collapse. We assumed that the Larsen B tributary glaciers were not losing significant mass before 2002.5 and set these discs to zero change between 1995.5 and 2002.5 accordingly. Any elevation changes that occurred away from former ice shelf regions were assumed to be part of a multi-decadal trend and associated mass changes were applied for the full modelling period. These were generally small and have little effect on our modelling. Although widespread glacier retreat is seen on the western Peninsula (Cook et al., 2005;Pritchard and Vaughan, 2007), thinning appears to be generally limited to a small section at the front of the glaciers and, importantly for this study, the pattern is changing relatively slowly with time (Kunz et al., 2012). Because the Larsen A and Prince Gustav glaciers are distant from our cGPS sites (150-300 km from Palmer), and the uplift signals due to mass changes in these areas and on the western Peninsula are linear, errors in the above mentioned assumptions have only a small effect on the conclusions drawn from our modelling based on the non-linearity in the Palmer record.

Elastic modelling
The elastic uplift was computed with the elastic component of the software VE-HresV2 output (Visco-Elastic High Resolution technique for Earth deformations) (Barletta et al., 2006;Barletta and Bordoni, in preparation), which is based on a compressible, spherically symmetric, self-gravitating Earth. Green's functions were spatially convolved with the ice loading discs according to the methods presented in Barletta et al. (2006). Load Love numbers, based on the PREM Earth structure (Dziewonski and Anderson, 1981), were computed in the centre of mass reference frame up to a maximum spherical harmonic degree of 3700 using VE-CL0V3RS v1.4 (Visco-Elastic Compressible LOVe numbER Solver) (Barletta and Bordoni, in preparation) and the degree-1 Love number was computed as described by Spada et al. (2011). By using the assumption that the elastic Love numbers become asymptotic after the maximum degree, the software implementing the High Resolution technique allows us to capture the loading concentrated on glaciers a few km wide (Barletta et al., 2006). In this way, the resolution is limited only by the resolution of the input loading discs.
A time series of modelled elastic uplift was computed at the location of Palmer and compared with the cGPS observations (Fig. 2). As the cGPS observations are recorded relative to an arbitrary reference height and the model output is relative to zero uplift at the start of the modelled time period, the cGPS observations have been offset to correct for this based on their pre-2002 mean. To account for the effects of centennial-or millennial-scale glacial isostatic adjustment (GIA) in the cGPS record we also estimated a 'background' vertical rate by subtracting the modelled elastic uplift rate from the cGPS uplift rate over the linear part of the Palmer record (1998)(1999)(2000)(2001)(2002). This gives the uplift rate due to any ice-mass changes prior to the start of our ice loading model, assuming an elastic-only response to post-1995 events. This rate was then included in our model-predicted time series so that model output could be directly compared with cGPS observations.

Viscoelastic modelling
The viscous uplift of the Earth in response to the ice-mass loss was computed using the compressible viscous component of the software VE-HresV2 output, which uses VE-CL0V3RS v1.4 to compute the elementary viscoelastic time-dependent Green's functions (convolved with Heaviside function) up to degree 1195, and assumes that at higher degrees they do not change with time so the combined Green's function is negligible. This was verified when choosing the maximum degree so that the results do not suffer from effects of truncation and, as for the elastic modelling, we were able to capture the response from glaciers a few km wide. We limit our study to a Maxwell rheological model. It is worth noting that models with alternative and more complex rheologies may also sufficiently explain the observations, however at present the dataset is too sparse to resolve or require them; we return to this point in the discussion.
We adopt a four-layer viscosity structure consisting of an elastic lithosphere, and a viscoelastic upper mantle, transition zone and lower mantle, as shown in Table 2. The density structure of the model consists of 31 finer layers with densities from the PREM Earth structure (Dziewonski and Anderson, 1981). We define a simple four-layer viscosity model as the limited data do not allow a more complex model to be resolved. This is discussed further in Section 5.
We searched for the range of plausible best-fit Earth models by varying the lithospheric thickness between 20 km and 160 km, and the upper mantle viscosity between 1 × 10 17 and 1 × 10 20 Pa s. Given that Simms et al. (2012) suggest a value of 1-2 × 10 18 Pa s for the South Shetland Islands, which lie closer to the active subduction zone than our study region, and typical values of man-  (Auriac et al., 2013;Lange et al., 2014;Sato et al., 2011), mantle viscosities below 1 × 10 17 Pa s are not thought to be physically realistic for this region of the Earth. All other parameters remained fixed. Below the upper mantle layer is a transition zone between 400 and 670 km depth with a fixed viscosity of 4 × 10 20 Pa s, as suggested by Sato et al. (2011) in their study of the Earth's response to ice-mass change in Alaska; and below this, a lower mantle with a fixed viscosity of 1 × 10 22 Pa s.
Sensitivity to different mantle layer thicknesses and a more complex Earth structure will be discussed later in Section 5.

cGPS constraints
The uplift time series output from the viscous model were added to the modelled elastic uplift and the background rate, which was recalculated as described in Section 3.1, this time by subtracting the modelled viscoelastic uplift rate from the cGPS uplift rate between 1998 and 2002. The resulting uplift time series for each Earth model in the parameter space was then compared, first of all, with the Palmer cGPS observations only. In order to determine the range of Earth models consistent with our data, the RMS misfit between the modelled uplift and the cGPS uplift was calculated and is shown in Fig. 3a.
In an attempt to place further constraints on the range of wellfitting Earth models, we repeated the viscoelastic modelling to calculate uplift at the six LARISSA cGPS locations (Fig. 1) for the full range of Earth models. By assuming that any lateral changes in Earth structure are minimal over the distance spanned by the cGPS stations (a maximum of 300 km), all sites can be used as constraints on a 1-D Earth model. The RMS misfit was computed again by comparing the model-predicted uplift (viscoelastic + background) with the cGPS observed uplift at all seven stations.
When computing the modelled time series at the LARISSA stations, which were not occupied prior to the Larsen B break-up, we assumed that the background rate previously calculated for Palmer was representative of the whole region. That is, we assumed a spatially constant background rate across all seven cGPS sites; this is supported by the closeness of fit of the initial Palmer-constrained model to most of the LARISSA sites (residual uplift rates in Table 3). Our assumption implies that the sites would have been uplifting at lower rates prior to 2002 and the time series would be non-linear, similar to that observed at Palmer. The implications of assuming a spatially-constant background rate are discussed later in Section 5.1. At present, we do not attempt to include geologic constraints on the background uplift rate, such as from marine limits and deglacial chronologies, as most sites (but not all) lack evidence suitable for long-term estimates of glacial isostatic adjustment.

Elastic modelling
The Palmer cGPS record displays significant non-linearity after 2002; however, the results of the elastic modelling show that even within the uncertainty bounds of the ice-mass change data (±1 m/yr), these changes in rate cannot be explained by elastic  −0.1 uplift only. In fact, more than five times the amount of observed mass loss (i.e. five times the mass loss shown in Fig. 1, applied to each disc) would be required to reproduce the magnitude of the observed uplift (modelled uplift shown by the green line in Fig. 2). This is not plausible and so we reject the possibility of missing ice unloading in our model, as the missing mass would not only need to be large, or be very close to Palmer, but also sustained from 2002 to present. Such a large signal would require a major ice shelf collapse or substantial glacier mass-loss adjacent to Palmer and neither scenario exists. Consequently, we conclude that less than half of the rapid increase in uplift at Palmer can be accounted for by elastic rebound, and examine if additional viscous uplift may help explain the remaining cGPS uplift signal.

Viscoelastic modelling constrained by PALM
The RMS misfit between the modelled uplift and the Palmer cGPS uplift is shown in Fig. 3a. The best fitting Earth models, lying within the 95% confidence limit of observational residuals, have a lithospheric thickness in the range 20-160 km and an upper mantle viscosity in the range 1 × 10 17 -2 × 10 18 Pa s. There is some trade-off between the two parameters, with thicker lithosphere models accompanying a lower viscosity mantle and vice versa. The Earth model with the lowest RMS misfit (4.67 mm) has values of 130 km and 7 × 10 17 Pa s. Computing the RMS again with a shortened time series ending in 2011 to coincide with the ice-mass change data, results in a best fitting model with a lithospheric thickness of 20 km. There is a possible offset in the Palmer time series during 1999 and only using data after this time (and recomputing the background rate appropriately) results in a best fitting model with a lithospheric thickness of 30 km. This highlights that the lithospheric thickness is poorly constrained within our model, although the upper mantle viscosity is robustly found to be less than 2 × 10 18 Pa s in all cases. For the Earth model with lowest RMS misfit to the Palmer time series, we compared the model-predicted velocities at all cGPS sites with observed velocities (Table 3). The model over-predicts uplift at CAPF by 2.8 mm/yr (compared with a 2-sigma uncertainty of 2.9 mm/yr) and under-predicts uplift at DUPT by 2.4 mm/yr, which is the only significant residual, but the misfit here is only ∼23% of the modelled uplift. The model performs well at the other four LARISSA sites with misfits in uplift rate < 2 mm/yr and within the 2-sigma observational error.

Viscoelastic modelling constrained with all cGPS records
Constraining the Earth model using uplift data from only one cGPS location results in an upper mantle viscosity that is relatively well constrained, but with a broad range of acceptable values of lithospheric thickness. Fig. 3b shows the RMS misfit between modelled uplift and the cGPS observed uplift for all sites. When using all the cGPS data to constrain the Earth models, the range of lithospheric thickness for the best fitting models narrows to 100-140 km and the acceptable range of upper mantle viscosity narrows to 6 × 10 17 -2 × 10 18 Pa s, as indicated by the 95% confidence limit. The Earth model with the lowest RMS misfit (4.38 mm) has values of 120 km and 1 × 10 18 Pa s.
The model-predicted time series for the best fitting Earth models in Figs. 3a and 3b are plotted against the cGPS time series in Fig. 4. For comparison, predicted time series using the "VM2" Earth model, the viscosity structure which accompanies the global ICE-5G GIA model (Peltier, 2004), are also plotted, along with time series calculated using an Earth model within the ranges suggested by Ivins et al. (2011) Ivins et al. (2011) model, and only 1 mm for VM2, compared with 123-130 mm for our best fitting models. In the supplementary material we compare the model-predicted uplift with GPS records at two further locations: OHI2 and SPPT (see Fig. 1 for locations), and this is discussed in Section 5.
The spatial distribution of model-predicted uplift rates (estimated over the same time period as the cGPS observed uplift rate, i.e. 2009.0-2013.0) for the elastic and viscous components are shown in Figs. 5a and 5b, respectively, the latter based on the best-fitting Earth model from Fig. 3b, as constrained by all seven cGPS sites. Fig. 5c shows the sum of panels a and b and represents the viscoelastic uplift rates including the uniform background rate, with cGPS uplift rates over-plotted (as per Table 3). The observed cGPS uplift rates are well reproduced by the model.

Earth model
Using the ice-mass change dataset and observations from seven cGPS stations we have been able to constrain a range of Earth models for the NAP. We first used Palmer station only to constrain the Earth model, as the background uplift rate due to long-term glacial isostatic adjustment is well constrained by the pre-2002 data, available only at this site. The addition of the six LARISSA cGPS sites narrows the overall range of best fitting Earth models to a lithospheric thickness between 100 km and 140 km and upper mantle viscosity in the range 6 × 10 17 -2 × 10 18 Pa s. Using a uniform background rate for the whole region may introduce some error in these results if the signal is in fact spatially variable, as suggested by Nield et al. (2012). Nonetheless, using the LARISSA cGPS data provides some verification for the inferences from the Palmer dataset. To test the sensitivity to a spatially constant background rate, we computed a new background rate based on that estimated from Palmer, but scaled at the other cGPS locations according to the spatial pattern of the W12a Antarctic GIA model . Computing the RMS again for all cGPS sites does not change the best fitting model and reduces the minimum RMS misfit by only 0.01 mm, giving further confidence in our results. Furthermore, comparing the best fitting model-predicted uplift with campaign GPS observations between 1997 and 2013 at the location of Spring Point, which were not included as constraints on the model, shows a qualitatively good fit and strengthens our conclusions (Fig. S3 in the supplementary material).
The Earth model with minimum RMS misfit in both cases (Figs. 3a and 3b) has a thick lithosphere (120-130 km) and a low upper mantle viscosity (7 × 10 17 -1 × 10 18 Pa s). However, it is important to note that the solution is non-unique and within the 95% confidence limit the RMS misfit varies by less than 1 mm, so any model within this limit can provide a reasonable fit to the data. The combination of thick lithosphere and low upper mantle viscosity is somewhat unexpected, as low viscosity regions are commonly accompanied by a thinner lithosphere (e.g. Auriac et al., 2013;Lange et al., 2014;Simms et al., 2012). However, as described in Section 4.2, computing the RMS misfit of Palmer using only pre-2011 data gives a best fitting lithospheric thickness of 20 km, which highlights our poor sensitivity to lithospheric thickness. Similarly, when discounting data before 1999.5 which may be affected by a possible offset in the Palmer time series, the best fitting lithospheric thickness is reduced to 30 km. In contrast, we robustly determine an upper mantle viscosity of less than 2 × 10 18 Pa s (the upper limit of the 95% confidence interval) is required to fit the available data. A low viscosity upper mantle is consistent with the back-arc setting and evidence of recent volcanism in the region.
Our range of Earth models is different to those determined by Ivins et al. (2011) for a larger region encompassing ours. Fig. 4 reveals that the Ivins et al. (2011) Earth model, when combined with post-1995 ice unloading, cannot explain the rapid uplift at Palmer since 2002. There are a number of possible reasons for this. First, Ivins et al. (2011) considered a single Earth model for a region approximately three times larger than ours, and hence their model is an average for this wider region. Second, the ice unloading sce-narios considered by Ivins et al. (2011) are based on relatively few observations and their Earth model may be compensating for ice load errors in poorly constrained regions. Third, Ivins et al. (2011) were limited in their ability to consider the non-linearity in the Palmer record, as their analysis required them to combine it with the GRACE time series which started after 2002, and therefore they treated it as a single linear rate over the post-2002 data period. Finally, it needs to be verified if our Earth model could satisfactorily fit the observations of the kind used by Ivins et al. (2011); if it cannot, then this may be an indication that use of higherorder constitutive theories that exhibit non-linear creep response functions (see Ivins and Sammis (1996), Fig. 7 therein), and whose material parameters may be independently tested by both laboratory and geophysical observation, must be considered.
The peak uplift predicted by our best fit Earth model is 47 mm/yr located in the Hektoria/Green glacier basin (Fig. 5c), corresponding to the geographical location of the largest mass loss (Berthier et al., 2012). The peak uplift is dominated by the elastic signal and has a small spatial extent, diminishing to 30 mm/yr or less within ∼ 30 km. Our nearest cGPS site is located at Foyn Point (FONP), ∼ 40 km away, where the observed uplift rate of 14.9± 2.7 mm/yr agrees well with the 16.4 mm/yr predicted by the model. Our rates may differ from the true uplift for parts of the model domain if the long-term background uplift rate is substantially different to the spatially constant term we have adopted; however, the closeness of the fit between the LARISSA cGPS data and our model predictions suggests our first approximation is reasonable.

Lateral variations in Earth structure
In using a 1-D symmetric Earth model we do not consider the effects of lateral heterogeneity in Earth structure. As our study region is small, there are unlikely to be large variations within the area covered by the cGPS stations. However, the long-term tectonic history of the region suggests that there may be a gradient in Earth structure along the length of the Peninsula (Barker, 2001;Larter et al., 1997). This is supported by the recent study of Simms et al. (2012) who predict a thinner lithosphere (15 km) for the South Shetland Islands, which lie 100 km off the northern tip of the AP. Due to the likely lateral variations in Earth structure, we did not include any cGPS data far from the site of largest ice loss as constraints (e.g. O'Higgins which lies 300 km to the north, OHI2 in Fig. 1), although we compare the model-predicted uplift to these cGPS observations in Fig. S2 of the supplementary material. Fig. S2 shows that the linear uplift recorded at O'Higgins cannot be explained by our model, both in terms of the magnitude of uplift and linearity of the time series. This is likely due to a combination of increased uncertainty in ice unloading near to O'Higgins over 1995-2001, the different Earth structure and our assumption of a spatially constant background rate.
The Earth models inferred here show that the NAP cannot successfully be modelled as part of a continent-wide GIA model, as the Earth structure is too different from that traditionally used for the rest of Antarctica (e.g., Whitehouse et al., 2012). Our work has important implications for forthcoming GIA models which incorporate 3-D Earth structure, and it identifies a location where upper Earth structure may be constrained by observations.

Sensitivity to a complex Earth structure
We found that the cGPS observations can be explained reasonably well by a simple four-layer viscous model, in which only the lithospheric thickness and upper mantle viscosity parameters were varied. The depth over which the load induces mantle flow was tested by decreasing the depth of the modelled upper mantletransition zone boundary to 350 km. The RMS of the two best fit models increased by 6-11% suggesting that mantle flow occurs to at least this depth. Increasing or decreasing the transition zone viscosity by an order of magnitude made less than 1% difference to the RMS for each best fit model. Our results are, therefore, not sensitive to changes in Earth model parameters below 400 km depth due to the small spatial extent of the load and observations. In view of this we consider the implications of an Earth model with a more complex structure in the top 400 km.
Several studies that have used a more complex Earth structure include a low viscosity zone (LVZ) directly beneath the elastic lithosphere (e.g., Sato et al., 2011). We performed some sensitivity tests to investigate whether such a model could provide a better fit to the data. Two thicknesses of the LVZ were tested (100 km and 200 km), along with several different ratios of LVZ viscosity to upper mantle viscosity (1:5, 1:10, 1:20) (six experiments in total). The results showed that whilst the single best fitting Earth model changed for each experiment, the overall range within the 95% confidence limit remained broadly the same. Furthermore, the minimum RMS misfit did not improve by more than 1%, demonstrating that a significantly better fit to the data could not be achieved with a more complex Earth model.
A more spatially extensive cGPS network might in the future enable a more complex Earth structure and lateral variations to be resolved and this network expansion is currently underway as an extension of the LARISSA project, with a permanent site now installed at Spring Point.

Sensitivity to ice model uncertainties
As described in Section 2, we modelled ice-mass change over a longer time period than is covered by the elevation change data by linearly extrapolating a constant rate of ice loss backwards and forwards in time from a few months after ice shelf break up to the present day, with no ramp up of mass change included. Studies have suggested that glaciers in the NAP that have accelerated following an ice shelf collapse remain at elevated speeds for decades. Rott et al. (2008) showed that Drygalski glacier, feeding the former Larsen A Ice Shelf, did not accelerate between 1999 and 2007, and Rott et al. (2011) state that the increased velocity of the Larsen A and Prince Gustav glaciers has so far been maintained 15 years after ice shelf collapse. The uncertainties in our ice model therefore relate to how quickly Larsen B glaciers accelerated to reach the 2002-2006 rates (Fig. 1b), and the acceleration history of Larsen A and Prince Gustav glaciers between 1995 and 1999.
To investigate this we simulated the acceleration of the glaciers in our ice loading model by linearly increasing the rate of mass loss from 0 m weq /yr at the time of ice shelf collapse, to full rates (as shown in Fig. 1) one year later for Larsen B and Larsen A/Prince Gustav glaciers separately. For Larsen B glaciers this ramp up scenario of ice-mass loss improved the RMS misfit by less than 5%, and similarly for Larsen A and Prince Gustav glaciers, confirming that the effect of errors in our ice loading assumptions is small.

Elastic effects of surface mass balance anomalies
Whilst the model output closely matches the cGPS time series overall, there are localized misfits. This is likely due to local time variable mass changes which are not included in our ice loading model. We investigated the possibility that these fluctuations are caused by an elastic response to variations in surface mass balance (SMB) over shorter periods of time than the DEMs allow us to resolve. SMB is dominated by iceberg calving, some melt runoff, and accumulation. Using the RACMO2.1/ANT27 dataset of SMB up to 2011.0 (Lenaerts et al., 2012), we removed the long term trend to obtain anomalies to the mean at each grid point, and converted them to a set of loading discs. The elastic response to the SMB anomalies was calculated at Palmer, and superimposed onto the time series for the best fit viscoelastic model from Fig. 3a. We perform this calculation with some caution due to the low resolution of the SMB model (27 km) compared with the small valley glaciers that dominate much of this region (Trusel et al., 2013). Nevertheless, Fig. 6 shows that the additional elastic response improves the fit between the modelled and observed time series and explains most of the short-term variations, which likely reflect seasonal signals and multi-year perturbations in SMB. The RMS misfit, calculated over the shorter time period of the SMB model, reduces marginally from 4.74 mm to 4.56 mm. It is not known how the effect of SMB anomalies could improve the fit at the other cGPS sites as the RACMO2.1/ANT27 model output is presently only available up until the end of 2010, providing minimal overlap with the LARISSA time series. A viscous response to SMB load changes is feasible but was not considered, as due to the fluctuating nature of the SMB loads it is likely to be small.

Conclusions
• Non-linear cGPS observed uplift from the NAP cannot be explained by the elastic response of the Earth to ice-mass loss events relating to Larsen A, Prince Gustav, and Larsen B ice shelf break-up.
• A linear Maxwell viscoelastic model can produce a good fit to the Palmer cGPS record, which constrains the upper mantle viscosity to less than 2 × 10 18 Pa s, but the lithospheric thickness remains poorly constrained.
• Shorter time series from the six cGPS stations of the LARISSA network verify this result, finding a best fitting Earth model with an upper mantle viscosity of 6 × 10 17 -2 × 10 18 Pa s and narrowing the lithospheric thickness to 100-140 km, although the analysis is limited by the assumption of a spatially uniform background rate.
• A more complex Earth structure, in terms of vertical stratification, could explain the observed data equally well, but additional cGPS stations are required to resolve this structure further.
• Reconciliation of present-day rates of bedrock uplift with average rates over longer timescales may require rheological models that exhibit a non-linear response.
• In regions of low upper mantle viscosity and on-going icemass change, cGPS data cannot be correctly interpreted without considering the viscoelastic response to present-day icemass change.