TY - JOUR
AU - Bordoni,, Andrea
AB - SUMMARY We consider the viscoelastic rheology of the solid Earth under the Antarctic Peninsula due to ice mass loss that commenced after the breakup of the Larsen-B ice shelf. We extend the previous analysis of nearby continuous GPS time-series to include five additional years and the additional consideration of the horizontal components of deformation. They show strong uplift from ∼2002 to 2011 followed by reduced uplift rates to 2018. Modelling the GPS-derived uplift as a viscoelastic response to ongoing regional ice unloading from a new ice model confirms earlier estimates of low upper-mantle viscosities of ∼0.3–3 × 1018 Pa s in this region but allows a wide range of elastic lithosphere thickness. The observed and modelled north coordinate component shows little nonlinear variation due to the location of ice mass change to the east of the GPS sites. However, comparison of the observed and modelled east coordinate component constrains the upper-mantle viscosity to be less than ∼9 × 1018 Pa s, consistent with the viscosity range suggested by the uplift rates alone and providing important, largely independent, confirmation of that result. Our horizontal analysis showed only marginal sensitivity to modelled lithospheric thickness. The results for the horizontal components are sensitive to the adopted plate rotation model, with the estimate based on ITRF2014 suggesting that the sum of residual plate motion and pre-2002 glacial isostatic adjustment is likely less than ∼±0.5 mm yr−1 in the east component. Satellite geodesy, Space geodetic surveys, Antarctica, Dynamics of lithosphere and mantle, Rheology: crust and lithosphere, Rheology: mantle 1 INTRODUCTION Knowledge of Antarctic glacial isostatic adjustment (GIA) is required to determine accurate estimates of ice sheet contribution to sea level change from space-gravimetry data (King et al. 2012). While most GIA models consider the solid Earth's rheology to vary only radially, a new generation of models is now being developed which considers the Earth's 3-D viscoelastic rheology reflecting information from seismic tomography (Whitehouse 2018). These seismic constraints do not, however, provide absolute constraints on mantle rheological structure. Geodetic studies of solid Earth response to recent ice mass change have begun to address this need in some areas of Antarctica, such as the Antarctic Peninsula. Glaciers in the northern Antarctic Peninsula (NAP) have been losing mass over the last two-and-a-half decades as a result of their acceleration and thinning following the retreat and collapse of several major ice shelves since 1993, such as Prince Gustav by 1993–1995 (Rott et al. 1996), Larsen-A in 1995 (Rack & Rott 2004) and Larsen-B in 2002 (Scambos et al. 2004). These represent a natural experiment that provides a rare opportunity to place absolute constraints on the region's viscoelastic rheology. The mass loss following the breakup of the Larsen-B ice shelf is well observed and was accompanied by rapid bedrock uplift that has been observed by GPS. Thomas et al. (2011) suggested that the response was due, at least in part, to the elastic deformation of the solid Earth; however, they were unable to satisfactorily reproduce the magnitude of the observed response with a purely elastic model. Building on this, Nield et al. (2014) closely reproduced the observations using an improved ice loading history combined with a viscoelastic solid earth model. They found a best-fitting earth model consisting of an upper mantle (down to 670 km) with viscosity tightly constrained to be within ∼6 × 1017 to 2 × 1018 Pa s, and an elastic lithosphere with an effective thickness largely unconstrained within the range ∼30–140 km. These results were based largely on the observed change in uplift rate over 1998.5–2013.0 at one GPS site, Palmer (PALM), located about 75 km from the rapidly changing glaciers (Fig. 1 inset). Figure 1. Open in new tabDownload slide Palmer GPS observations (grey dots) with viscoelastic predictions based on the different ice load change scaled uniformly (reduction: 0, 25, 50 and 75 per cent) from beginning 2011 for two different earth models—A: lithospheric thickness 30 km with upper-mantle viscosity 9 × 1017 Pa s (dash line) and B: lithospheric thickness 130 km with upper-mantle viscosity 7 × 1017 Pa s (solid line). The GPS and model time-series have been shifted to align over the period up to 2002. The inset shows study region and GPS sites locations PALM (or PAL2), DUPT and VNAD, corresponding to Palmer station, Duthiers Point station and Vernadsky station, respectively, and the study area (red line). The GPS and model time-series have been shifted to align over the period up to 2002. Figure 1. Open in new tabDownload slide Palmer GPS observations (grey dots) with viscoelastic predictions based on the different ice load change scaled uniformly (reduction: 0, 25, 50 and 75 per cent) from beginning 2011 for two different earth models—A: lithospheric thickness 30 km with upper-mantle viscosity 9 × 1017 Pa s (dash line) and B: lithospheric thickness 130 km with upper-mantle viscosity 7 × 1017 Pa s (solid line). The GPS and model time-series have been shifted to align over the period up to 2002. The inset shows study region and GPS sites locations PALM (or PAL2), DUPT and VNAD, corresponding to Palmer station, Duthiers Point station and Vernadsky station, respectively, and the study area (red line). The GPS and model time-series have been shifted to align over the period up to 2002. The inferred low estimates of upper-mantle viscosity (UMV) were lower than previously suggested for this region over millennial timescales (Simms et al. 2012) and are 1–2 orders of magnitude lower than suggested 500 km further south from a similar geodetic study (Zhao et al. 2017). However, a low UMV is not entirely unexpected due to the active tectonic setting of the NAP that passes from a subduction zone along the South Shetland Trench, located north of the South Shetland Islands (Barker et al. 1991; Barker & Austin 1998; Nield et al. 2014). Extending the analysis to include time-series from six additional continuous GPS sites covering the shorter period between 2009 and 2013 allowed Nield et al. (2014) to confirm the modelled low viscosity mantle and suggest a preference for higher lithospheric thickness (LT), although this analysis was limited by the assumption of a spatially uniform background rate (the rate of uplift prior to commencement of recent ice unloading). These initial experiments show the value of spatially distributed bedrock deformation data in constraining mantle rheology, and we build on these here using new data sets. The solid Earth viscoelastic deformation has both a vertical and a horizontal component. In general, the pattern observed in the vertical component is proportional to the distribution of the ice loading, while the horizontal pattern provides information regarding the origin of the ice loading (Velicogna & Wahr 2013; Barletta et al. 2018). Both components are sensitive to Earth rheology and horizontal GPS time-series are 2–3 times more precise than vertical time-series. The viscoelastic deformations of the NAP in response to current ice mass changes have not yet been analysed using the horizontal geodetic components for this region. Horizontal GPS observations are often ignored in the study of GIA because of the challenge associated with accurately removing tectonic plate motion from the observations. However, given the limited GIA observations in Antarctica, and absolute constraints on rheological structure, adding horizontal observations of GIA is highly desirable. This study focuses on understating the solid Earth rheology and temporal variation in ice loading in the NAP based on daily GPS coordinate time-series. We extend the GPS time-series by 5 yr to 2018. We consider vertical and, for the first time in the NAP, horizontal displacement fields from three continuous GPS stations. We update the ice loading history used by Nield et al. (2014) with a high-resolution ice loading history from Rott et al. (2018) as input to a viscoelastic solid earth model to compare with observations, and arrive at new constraints on NAP solid Earth rheology. 2 DATA AND METHODS 2.1 GPS We consider GPS data within the study region covering the time from 1998.5 when continuous data collection commenced to the end of 2018. Our analysis includes the long-running site PALM (Fig. 1), which commenced operating in 1998 and was the focus of the studies mentioned above. Since 2006, PAL2 has operated alongside PALM and its time-series shows almost identical site behaviour to PALM during the overlapping period. The main difference between the two series is that PALM has a short (∼2 month) data gap in 2014. Other sites considered are from the LARISSA network which was deployed from 2009 and include sites DUPT and VNAD (Fig. 1; Nield et al. 2014). GPS data were analysed using GIPSY 6.3 using approaches very similar to Wolstencroft et al. (2015). In brief, this includes adopting a precise point positioning methodology holding fixed fiducial-free NASA JPL ‘repro2’ clocks and orbit files. Time-series were then transformed into the ITRF2008 (Altamimi et al. 2011), which has its origin, at secular timescales, near to the centre of mass of the entire Earth system (Dong et al. 2003). Site time-series were converted to a local east–north–up coordinate system for further analysis. Within the GIPSY analysis, we model deformation due to solid Earth body tides and ocean tide loading (using TPXO7.2 ocean tide model) and the pole tide according to IERS2010 conventions, and use calibrated absolute antennas phase centre variations (igs08_1758.atx). Data from low elevation satellites were down-weighted in the solutions. Tropospheric zenith delays and gradients were estimated once per 5 min, using the Vienna Mapping Function 1 (Boehm et al. 2006) to map tropospheric delays between satellite elevation angles and the zenith. There are no offset evidents in the vertical time-series. The GPS site coordinates are affected by global elastic deformation over timescales of weeks and longer due to redistribution of mass in the atmosphere, oceans and hydrosphere, with the largest of these in Antarctica being associated with the atmosphere (Santamaría-Gómez & Mémin 2015). We, therefore, subtracted modelled 3-D site displacements due to atmospheric loading changes within the centre of figure frame (http://massload.net). We consider both GPS coordinate and velocity time-series. To produce continuous GPS velocities, we make use of the ensemble empirical mode decomposition (EEMD) method which is used to separate the temporally varying trend from the GPS observations (Wu & Huang 2009). The EEMD method is based on empirical mode decomposition (EMD; Huang et al. 1998), an adaptive method for studying nonlinear and non-stationary time-series. In EMD, a time-series x(t) is decomposed into a finite number of intrinsic mode functions (IMFs) Cj(t), j = 1, 2, …, n, using the local maxima or minima envelope. After any oscillations in the time-series are separated, the residual R(t), which has at most one extremum, is derived as an intrinsic trend of a non-stationary and nonlinear time-series (Wu et al. 2007). The EEMD method works on even-spaced time-series, where data gaps exist and we filled them using an extrema linear interpolation scheme. The linear interpolation does not introduce artificial extrema, and this interpolation has minimal impact on the final result. The EEMD approach is a noise-assisted EMD. In EEMD, multiple noise realizations are added to the time-series x(t), from which, an ensemble average of the corresponding IMFs are extracted to yield scale-consistent signals. By using the ensemble mean of IMFs, non-physical oscillations due to random data errors can be reduced and thus low-frequency modes, especially the trend of the time-series are more accurate (Ezer et al. 2013; Ji et al. 2014; Chen et al. 2017). To estimate the uncertainties in the determined IMFs of each time-series, we randomly pick a value of the observation within two-standard deviations at each time step to represent a sample observation, and then decompose each resulted time-series. We repeat this procedure 100 times and obtain the means of each IMF, the secular trend and their spread (uncertainty; see Supporting Information Fig. S1). In this study, we used the IMF6 (decadal/interdecadal scale) to compute the rates for the respective GPS stations as which is a timescale relevant to this study. The IMF5 (quasi-biennial/triennial scale) has a higher frequency than IMF6 and does not affect the longer-term velocities focused on here. EMD and EEMD methods have been widely used in geoscience applications (Huang & Wu 2008). Since time-series data gaps may introduce spurious low-frequency variability into the final decompositions (affecting the final estimation of the GPS velocities), we only considered vertical velocities for sites with no data gaps, limiting our analysis to PAL2, DUPT and VNAD. These sites are all no closer than ∼50–75 km from the main glaciers losing mass, and we regard them as being far-field in that elastic effects do not dominate the uplift at these distances (Nield et al. 2014). For horizontal component analysis, we consider just the time-series as the signals are smaller than the vertical by a factor of ∼3–4 and computing velocities amplifies noise. We also use the PALM time-series instead of PAL2 because of the extended observations. 2.2 Ice loading model We divide the ice history into two phases (X and Y) based on the two different sources of ice elevation change. We use the ice history described in Nield et al. (2014) as input to the viscoelastic modelling for phase X: 1995.5–2011.5 (divided into linear rates over 1995.5–2001.0; 2001.0–2006.0; 2006.0–2011.5) and a new ice loading history taken from Rott et al. (2018) for phase Y: 2011.5–2018.0. The new ice loading model is derived from the surface elevation change (SEC) data described in Rott et al. (2018). These SEC data are divided into two periods: May–June 2011 to June–July 2013 and June–July 2013 to July–August 2016 and we consider linear rates over 2011.5–2013.5 and 2013.5–2018.0. We downsample their 12 m resolution SEC data to 200 m resolution and convert elevation change into ice-mass changes using a density of 900 kg m−3 for the feeding glaciers of the respective region. For the loading computations, we convert mass change data into a set of 155 043 discs with 0.1 km2 area for input to the solid earth model, where each disc height represents an ice mass gain or loss following Nield et al. (2014). Data around the Larsen-A region was only available for 2013.5–2018.0, so we consider use the linear rate for 2013.5–2018.0 for the whole 2011.5–2018.0 period. For the earlier periods, the Nield et al. (2014) ice loading history for the Larsen-B region was used with linear rates of change over 2001.0–2006.0 and 2006.0–2011.5 based on the data availability. For the Larsen-A region, they used observations over 1995.5–2001.0 and held this linear rate constant through to 2011.5. For regions outside the Rott et al. (2018) data set, for the period 2011.5–2018.0, we use the Nield et al. (2014) rates. The spatial pattern of the ice mass changes data of Nield et al. (2014) and Rott et al. (2018) is shown in Supporting Information Figs S2 and S3. The reader is referred to Berthier et al. (2012), Nield et al. (2014) and Rott et al. (2018) for more details about the complete ice loading model used in this study. 2.3 Viscoelastic modelling 2.3.1 Elastic modelling We separate the modelling of the solid Earth response into elastic and viscous components, using the approach previously employed in this region by Nield et al. (2014). The elastic deformation time-series was computed by the elastic component of the software VE-HresV2 (Visco-Elastic High-Resolution technique for Earth deformations; Barletta et al. 2006, 2018). The load Love numbers were computed to a maximum spherical harmonic degree of 3700 based on a compressible, self-gravitating Earth using VE-CL0V3RS v1.4 (Visco-Elastic Compressible LOVe numbER Solver; Barletta et al. 2018) with a PREM-based structure (Dziewonski & Anderson 1981). We compute the Green's functions by using load Love numbers, which were spatially convolved with ice loading discs according to the methods presented in Barletta et al. (2006). The model was evaluated at semi-annual time steps over 1995.5–2018.0. While Nield et al. (2014) focused entirely on the vertical component of deformation, we also compute the horizontal deformation (Barletta et al. 2018) to allow comparison with GPS horizontal component. 2.3.2 Viscous modelling For the viscous component of the modelling, we again follow the approach of Nield et al. (2014), using the compressible viscous component of the software VE-HresV2 (Barletta et al. 2006, 2018), VE-CL0V3RS v1.4 to compute the viscoelastic time-dependent Green's functions up to maximum spherical harmonic degree 1195, which is adequate for the temporal and spatial viscous response scale considered here. We adopt a similar viscoelastic model setup to Nield et al. (2014) assuming a 4-layer viscoelastic Earth structure with a purely elastic lithosphere, and a viscoelastic upper mantle, transition zone and lower mantle with linear Maxwell rheology. This model is forced by the ice history model described above. Modelling is performed in the centre of mass of the solid-Earth system, which is different from the GPS reference frame. However, our focus here is on nonlinear motion and hence much of the already small difference due to geocentre movement is removed. We tested a total of 345 earth models in this study by varying the lithospheric thickness between 20 and 130 km and the viscosity of the upper mantle between 0.1 × 1018 and 300 × 1018 Pa s. The bottom of the upper mantle was fixed to 400 km and the modelled transition zone extends from 400 to 670 km with a viscosity of 4 × 1020 Pa s, and a lower mantle with a viscosity of 1 × 1022 Pa s. For all model configurations, we computed displacement time-series at semi-annual periods over 1995.5–2018.0. Again, we compute both horizontal and vertical deformations. We sum the elastic and viscous components to compare with the GPS time-series, note that we consider the data uncertainties for all RMS calculations. Our GPS time-series have additional deformation contribution caused by ice unloading history not considered in our loading model (e.g. due to LGM deglaciation) and we define this at the background rate. For the vertical component, we estimate this ‘background rate’ (BGR) for each viscoelastic model by deriving a trend from the difference of the observed and each model time-series over 1998.5–2002.0 (Nield et al. 2014). We adopt the same vertical BGR for all three sites used in our study prior to further analysis. We then add this background rate to our model predicted uplift rates. The horizontal components of the GPS site time-series is dominated by plate motion. To study the horizontal component due to present-day ice mass changes, we remove the plate contribution from our observations using the ITRF2014 plate motion model (Altamimi et al. 2017). The ITRF2014 plate velocity for the GPS sites with a no-net Rotation (NNR) frame from the University NAVSTAR Consortium (UNAVCO) Plate Motion Calculator platform (https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html#models). Supporting Information Table S1 gives the different GPS-site plate motions using this ITRF2014 model. There are important limitations to geodetically derived plate motion models as they are estimated using GPS data that contain uncorrected GIA (King et al. 2016) and post-seismic deformation (King & Santamaria-Gomez 2016) and hence they are subject to potential bias. Furthermore, we are unable to unambiguously remove background horizontal GIA from before 2002 (i.e. due to deglaciation following Last Glacial Maximum), for which models suggest horizontal rates in this region of around up to 0.84 mm yr−1 but of varying direction (King et al. 2016), and so we adopt zero signal and test the sensitivity to these assumptions in Section 3.2. 3 RESULTS 3.1 Vertical component results and model fit analysis Fig. 1 shows the GPS height time-series for PALM. This reveals that the increased rate of uplift, compared to before 2002 and previously shown to 2011 continues through to the end of the time-series at 2018.0. The longer time-series reveals an apparent slowing of the rate around 2010 or 2011 that was previously not clear. The ability to separate model predictions of deformation based on different rheological models depends entirely on nonlinearity in the time-series and so additional and substantial changes in the rate of mass loss and bedrock uplift provide the potential for improved understanding of NAP rheology. We first compared these observations to an extended version of the Nield et al. (2014) ice loading model, assuming no change in the rate of ice mass unloading after 2011. Combining this with the best fitting earth model of Nield et al. (2014), with an LT of 130 km and UMV of 7 × 1017 Pa s, provides a poor fit to the GPS data after ∼2011 (yellow line). However, by reducing the post-2011.0 unloading rate by 50–75 per cent from the 2006.0–2011.0 rate, we are able to fit the GPS data more closely (cyan and red lines). Nield et al. (2014) showed that alternative viscous models could also provide an adequate fit to the GPS data, notably those with low effective lithospheric thicknesses (as low as ∼20–30 km). Predictions based on a model with LT of 30 km and UMV of 9 × 1017 Pa s are shown with dashed lines in Fig. 1. With this configuration, all variants of the post-2011 ice loading significantly over-predict the amount of uplift since 2011. If the Earth rheological model is broadly correct, the region would need to have switched to mass increase since ∼2011 for this model configuration to fit the data, which was not the case (Rott et al. 2018; Shepherd et al. 2018). We then updated the ice loading history after 2011 as described in Section 2.3 and compared a wide range of model predictions (varying the LT and UMV) with the PALM time-series. We computed RMS misfit between the PALM GPS (1998.5–2018.0) and modelled the uplift time-series, with the results shown in Fig. 2(a). The results narrow the preferred earth model parameters obtained by Nield et al. (2014) to ∼50–130 km and ∼0.9–3 × 1018 Pa s. For the PALM station, the earth model with the lowest RMS misfit has a value of 65 km and 3 × 1018 Pa s (see Supporting Information Fig. S4). Figure 2. Open in new tabDownload slide (a) RMS misfit between PALM GPS and model predicted uplift. (b) RMS misfit between all stations (PAL2, DUPT and VNAD) and model predicted rates. c–e) Uplift rates from the EEMD analysis of GPS time-series with the best-fitting viscoelastic earth model. GPS observations (black line), GPS uncertainties (grey dash line) and best-fitting model (yellow line). Note that data collection commenced at PAL2 before the other sites and the y-axis differs between sites. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Figure 2. Open in new tabDownload slide (a) RMS misfit between PALM GPS and model predicted uplift. (b) RMS misfit between all stations (PAL2, DUPT and VNAD) and model predicted rates. c–e) Uplift rates from the EEMD analysis of GPS time-series with the best-fitting viscoelastic earth model. GPS observations (black line), GPS uncertainties (grey dash line) and best-fitting model (yellow line). Note that data collection commenced at PAL2 before the other sites and the y-axis differs between sites. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) To further examine the changes in the uplift rate and constrain the Earth rheology, we expand the analysis to include the three sites with continuous data in this region since 2006 (PAL2, DUPT and VNAD; Fig. 1) and consider the EEMD-filtered GPS time-series in terms of uplift rates. Studying changes in uplift rates, rather than changing uplift time-series, allows us to mitigate the effects of unknown background rates, otherwise assumed to be spatially constant. We also use the updated (phase Y) ice loading model to predict the uplift rate using viscoelastic solid earth models. We compared the GPS observations with the predictions from 345 different viscoelastic models. We computed the RMS misfit between the GPS and model rates for each station after considering the BGR calculation from the PALM site (Fig. 2b and Supporting Information Fig. S5). The analysis of the GPS station data, either individually or combined, show good agreement with UMVs of ∼0.3–3 × 1018 Pa s but show poor sensitivity to the lithospheric thickness. These findings confirm the results of Nield et al. (2014) that used a much shorter data set. The minimum RMS misfit is found for an LT of 105 km with a UMV of 9 × 1017 Pa s. However, the DUPT site shows good agreement with a slightly higher UMV (∼5 × 1018 Pa s) compared to the other stations (see Supporting Information Fig. S5b), and this might be due to its greater closeness to major ice mass loss areas (eastern Peninsula). We present the resulting GPS site velocity time-series with the best fit model rates in Figs 2(c)–(e). Time-series of the three continuous GPS stations reveal a reduction in uplift rates since ∼2011 compared to acceleration before then (Fig. 1). The step-wise velocity changes in the model reflect the periods where different time-constant ice loading rates apply. The reduction in uplift rates occurs slightly earlier at PAL2 compared with DUPT and VNAD, possibly due to differing proximity to the glaciers losing mass or a changing geographical pattern of mass loss over time. In Figs 2(c)–(e), we note a possible slight increase in uplift rates in the early part of the VNAD series although this is uncertain. The maximum observed uplift rate of ∼12.4 ± 0.2 mm yr−1 was recorded at the DUPT station between 2009 and 2010. All the GPS observations show a bias between the model and observed uplift rates inherited from the viscoelastic model before 2009 and this is particularly clear at DUPT. The reduction of the observed rates from the beginning to the end of the uplift rate time-series is ∼37, 20 and 25 per cent for PAL2, DUPT and VNAD, respectively. It suggests that ice mass loss reduced between 2011.0 and 2018.0. 3.2 Horizontal component results and model fit analysis Fig. 3 presents the horizontal displacement of the long-running PALM site after removing the ITRF2014 plate motion model. The horizontal information is a useful supplement to the vertical for areas with a more complex ice mass change pattern, such as the NAP where the mass change pattern originates from many different glaciers with different signals. The horizontal displacements in Figs 3(a) and (b) are generally directed towards the west as is expected given the major ice loss area is due east of PALM (see Supporting Information Fig. S6). As with the uplift, the east components show clear changes after 2002.0, which is indicative of rapid mass loss acceleration after the Larsen-B breakup. Figure 3. Open in new tabDownload slide PALM GPS north (a) and east (b) component time-series, after the subtraction of the ITRF2014 plate motion model. Figure 3. Open in new tabDownload slide PALM GPS north (a) and east (b) component time-series, after the subtraction of the ITRF2014 plate motion model. We compared our GPS data after removing the plate motion with the predictions from the range of numerical viscoelastic earth models and calculated the RMS misfit for each of the north and east coordinate components. There is a possible offset in the horizontal components of the Palmer time-series during 1999, so we only consider data from then onward. The RMS misfits between the observed and model-predicted north component for PALM individually and for all sites combined, are shown in Fig. 4, with the other two sites shown in Supporting Information Fig. S7. Aside from the DUPT site, we find that none of the models provide satisfactory agreement with the observations, as indicated by the 95 per cent confidence limit, this is likely due to errors in the plate motion model. The DUPT site, however, shows good agreement with almost 90 per cent of the earth models used in this study (see Supporting Information Fig. S7). Overall, the north component shows little sensitivity to the change of earth model properties, likely due to the location of our GPS stations to the west of the regions losing mass. Consequently, we focus only on the east component in the following discussion. Figure 4. Open in new tabDownload slide RMS misfit between GPS and modelled predicted north component time-series. Figure 4. Open in new tabDownload slide RMS misfit between GPS and modelled predicted north component time-series. Fig. 5 shows the RMS misfit between observed and model-predicted east components. The individual analysis PALM and DUPT show good agreement with UMV of ∼0.5–9 × 1018 Pa s, but show poor sensitivity to the effective LT. The minimum RMS misfit is found for LT 115 km and LT 25 km with UMV of 5 × 1018 Pa s for PALM and DUPT sites, respectively. The RMS fit at VNAD is poor, probably due to its proximity to the glaciers losing mass and sensitivity to errors in the loading model. This affects the combined solution presumably. The prediction based on the best-fitting earth model constrained from PALM and DUPT is plotted against the GPS east component time-series in Fig. 6. Figure 5. Open in new tabDownload slide RMS misfit between GPS and modelled predicted east component time-series. Figure 5. Open in new tabDownload slide RMS misfit between GPS and modelled predicted east component time-series. Figure 6. Open in new tabDownload slide GPS eastward observations after subtraction of plate motion from ITRF2014 and model predicted east time-series at each GPS station for the best-fitting viscoelastic earth model of PALM and DUPT present in Figs 5(a) and (b). Figure 6. Open in new tabDownload slide GPS eastward observations after subtraction of plate motion from ITRF2014 and model predicted east time-series at each GPS station for the best-fitting viscoelastic earth model of PALM and DUPT present in Figs 5(a) and (b). The studies of Klemann et al. (2008) and King et al. (2016) showed that the GPS derived plate motion models may be biased due to mismodelled or unmodelled GIA signals. King et al. (2016) demonstrate plate motion model differences up to 0.84 mm yr−1 depending on the nature of the unmodelled GIA signal. Different plate motion models also vary by around 1 mm yr−1 due to sensitivity to site selection, including post-seismic deformation effects (King & Santamaria-Gomez 2016). To explore this potential error, we recomputed the RMS misfit for the east component by considering four scenarios based on the addition or subtraction of 0.5 and 1 mm yr−1. Supporting Information Figs S8–S11 show the RMS misfit between the observations and model predictions for these four scenarios. The DUPT and VNAD show a similar pattern of RMS misfit for the four scenarios, as does the combined RMS. However, PALM shows sensitivity to the addition and subtraction of 0.5 and 1 mm yr−1. By inspecting the misfit of these four scenarios for the PALM site, we note that the RMS misfit increases if we add or subtract the values for each scenario. Moreover, the only scenario that provides reasonable fit is with +0.5 mm yr−1 (Supporting Information Fig. S8) but the results are inconsistent, with DUPT preferring lower viscosity and PALM preferring higher, and PALM results are not consistent with the finding using the vertical component base on the previous (Nield et al. 2014) and on this study. This indicates that the PALM site motion is accurately described by the ITRF2014 plate motion and that any plate motion model error or long-wavelength post-LGM GIA signal has a sum within ± 0.5 mm yr−1. 4 DISCUSSION Our continuous GPS observations show that the uplift rate decreased after 2011, while still being substantially higher than the rate before 2002, and our inferred (background) rate of uplift for the period before 1995 (∼2–3 mm yr−1 subsidence). Assuming a mantle relaxation time of a few years, this is compatible with a reduction in mass loss over this period. The reduction in GPS uplift rates since ∼2011 indicates that ice mass loss has reduced in this region since then, which agrees with the ice surface elevation and velocity change data presented by Rott et al. (2018). They found that overall ice-mass loss reduces around the NAP from mid-2011 to mid-2016. Analysis of time-series from three continuous and extended GPS sites tightly constrains the UMV in the range ∼0.3–3 × 1018 Pa s but they do not help in reducing the uncertainty on the LT which remains poorly constrained. Since 2015 all three sites show more stable rates of uplift compared with the rates during 2011–2015. The preferred model suggests a pre-1995 subsidence of ∼2–3 mm yr−1 at PALM. This contrasts with predicted uplift in this region from a range of millennial-scale post-LGM GIA models (Thomas et al. 2011). However, this region has seen centennial-scale snow accumulation increases, which could explain this rate of subsidence. Indeed, the modelling of Nield et al. (2012), for a relatively low viscosity earth model (5 × 1019 Pa s upper mantle) and reconstructed ice loading changes since 1850, suggests subsidence of ∼3.5 mm yr−1 for this region; they did not test lower viscosity models which will both increase the maximum rate of subsidence but also reduce the spatial extent of uplift. After removal of the ITRF2014 plate motion model, the horizontal GPS displacements are directed towards the west–southwest (see Fig. S6), consistent with major ice mass loss ongoing in the eastern part of the Peninsula. Our analysis of the east component provides new evidence, independent form the uplift time-series, that this region is underlain by the UMVs suggested by Nield et al. (2014) and from this study based on uplift time-series alone. The PALM and DUPT sites constrain the viscosity of the upper mantle to less than 9 × 1018 Pa s but, like the vertical component, show poor sensitivity to the effective thickness of the elastic lithosphere. Our east component constrains the viscosity range to be slightly higher than the vertical component, although substantially overlapping with it; this higher range may be due to different noise levels in the data sets or poorly sampled mass changes on the west coast of the Peninsula which has different effects on the different coordinate components. Further, our exploration of uncertainty in background GIA and plate motion indicates that the ITRF2014 model might precisely describe the motion of PALM prior to 2002 and suggests that the combination of plate model error and post-LGM GIA signal has a magnitude of ∼ ± 0.5 mm yr−1 for this region. We use a simple 1-D earth model with linear Maxwell rheology, varying just the elastic thickness of the lithosphere and viscosity of the upper mantle. We do not consider the effects of lateral heterogeneity in the solid Earth structure as our study area is small, and three are unlikely to be significant lateral variations within the area covered by the GPS sites. However, Barker et al. (2001) suggest there is a lateral variation in the solid Earth structure along the length of the Peninsula and this has been further supported by recent seismic tomographic analysis (Lloyd et al. 2019). Previous experiments in this region (e.g. Nield et al. 2014) suggested no sensitivity to modelled solid Earth properties from the transition zone downwards. We have not revisited these analyses although we note that an effect of lateral heterogeneity in Earth structure and a sensitivity to the mantle transition zone viscosity has been identified in the West Antarctic region (Barletta et al. 2018). The horizontal displacement field used to constrain the viscosity of the upper mantle used in this study can eventually be used to constrain a 3-D model of GIA in the future; such analyses will occur as more GIA models include 3-D rheological structure. 5 CONCLUSIONS With a combination of geodetic observations and viscoelastic modelling, we have improved the absolute constraint on the rheology of the upper mantle and elastic lithosphere in the Antarctic Peninsula region, assuming a linear Maxwell rheology. Such information is critical in providing robust models of Antarctic GIA that consider lateral variations in Earth properties. Our analysis extends the previously analysed time-series by 5 yr, capturing a reduction in uplift rates. We also fill a gap in the analysis of this natural experiment using horizontal displacement fields derived from GPS. Our consideration of the horizontal displacement field of GPS observations is the first such application over this region in the study of the viscoelastic response to present-day ice mass unloading. We computed model predictions based on the four significant ice mass change periods after the breakup of the Larsen B ice shelf: 2001.0–2006.0, 2006.0–2011.5, 2011.5–2013.5 and, 2013.5–2018.0 (Nield et al. 2014; Rott et al. 2018), in addition to ice loading changes from Larsen A:1995.5–2011.5 and 2011.5–2018.0 (Nield et al. 2014; Rott et al. 2018). Our continuous GPS uplift time-series indicates that the ice mass loss has reduced since ∼2011 around this region and this is supported by the ice model based on Rott et al. (2018) over 2011.5–2016.5 and then viscoelastic modelling using that ice model extended by 5 yr to 2018.0. The uplift data sets suggest the UMVs of ∼0.3–3 × 1018 Pa s for this region. A comparison of the model predicted east component with the GPS observations suggests an upper limit of the UMV of ∼<9 × 1018 Pa s for this region, consistent with that estimated from the vertical component in both this study and that of Nield et al. (2014) and this study based on analysis of the vertical component. Our research indicates that the GPS horizontal component may be a useful supplement to the vertical displacements to study viscoelastic deformation in this region, particularly with further substantial changes in mass loss rates or the spatial location of mass loss. Further, those using GPS data sets to estimate Antarctic plate rotation should be aware that horizontal GPS velocities in this region are strongly affected by viscoelastic deformation. Exploration of the potential bias in the plate motion estimation from the ITRF2014 model suggests that this plate model reasonably robustly describes the PALM site plate motion prior to 2002 and suggests that any potential plate model error or post-LGM GIA signal has a sum within ∼±0.5 mm yr−1 of zero for this region. GIA cannot be considered as time constant in this region; its time variability must be considered for studies using GRACE data. Further investigation with high spatial resolution and over extended time periods will allow more complex features of the solid Earth to be studied, both in this region and elsewhere (e.g. Zhao et al. 2017; Barletta et al. 2018). SUPPLEMENTARY INFORMATION Figure S1. Top: the IMFs of time-series in PAL2, DUPT and VNAD. The solid black line shows each mean IMF of 100 randomly constructed time-series; the shading shows the uncertainty of each IMF, note that only IMF6 is used to compute the rates for the respective GPS sites. Bottom: the time derivative of the secular trend shown in the top panel. Figure S2. The rate of observed ice mass change given in metres water equivalent per year. (a) Ice mass change rate of Larsen A and B region for the period 1995.5 to 2002.5. (b) Larsen B region ice mass change rate from 2001.0 to 2006.0. (c) Larsen B region ice mass change rate from 2006.0 to 2011.5. Blue line shows the grounding line and axes are Antarctic Polar Stereographic X, Y (km). Figure S3. The rate of observed ice mass change given in metres water equivalent per year. (a) Ice mass change rate of Larsen A and B region for the period 2013.5 to 2016.5. (b) Larsen B region ice mass change rate from 2011.5 to 2013.5. (c) Larsen B region ice mass change rate from 2013.5 to 2016.5. Blue line shows the grounding line and axes are Antarctic Polar Stereographic X, Y (km). Figure S4. Uplift of PALM time-series with the best-fitting viscoelastic earth model and background rate present in Fig. 2(a). Figure S5. RMS misfit between GPS and modelled predicted uplift rates. Figure S6. GPS observed horizontal deformation rates after subtraction of plate motion from ITRF2014. The arrows indicate the size and direction of the deformation at three sites. Figure S7. RMS misfit between GPS and modelled predicted horizontal north component time-series. Figure S8. RMS misfit between GPS and modelled predicted horizontal east component time-series with 0.5 mm yr−1 plate motion subtraction. Figure S9. RMS misfit between GPS and modelled predicted horizontal east component time-series with 0.5 mm yr−1 plate motion addition. Figure S10. RMS misfit between GPS and modelled predicted horizontal east component time-series with 1 mm yr−1 plate motion subtraction. Figure S11. RMS misfit between GPS and modelled predicted horizontal east component time-series with 1 mm yr−1 plate motion addition. Table S1. Plate motion estimation from plate model ITRF2014. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS NHS is a recipient of an Australian Government Research Training Program Scholarship and the School of Technology, Environments and Design, University of Tasmania top-up scholarship. GPS data used in this study were downloaded from public archives at UNAVCO and International GNSS Service. The ITRF2014 plate motion estimation was from the UNAVCO public archives database. We thank NASA JPL for making GIPSY available. MAK was supported by an Australian Research Council's Future Fellowship (project number FT110100207) and the Australian Research Council's Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001). This work was also supported by the Natural Environment Research Council grant number NE/L006065/1. VE-CL0V3RS is a self-funded project by VRB and AB. REFERENCES Altamimi Z. , Collilieux X., Métivier L., 2011 . ITRF2008: an improved solution of the international terrestrial reference frame , J. Geod. , 85 , 457 – 473 . 10.1007/s00190-011-0444-4 Google Scholar Crossref Search ADS WorldCat Crossref Altamimi Z. , Métivier L., Rebischung P., Rouby H., Collilieux X., 2017 . ITRF2014 plate motion model , Geophys. J. Int. , 209 ( 3 ), 1906 – 1912 . 10.1093/gji/ggx136 Google Scholar Crossref Search ADS WorldCat Crossref Barker D.H. , Austin J.A., 1998 . Rift propagation, detachment faulting, and associated magmatism in Bransfield Strait, Antarctic Peninsula , J. geophys. Res. , 103 , 24 017 – 24 043 . 10.1029/98JB01117 Google Scholar Crossref Search ADS WorldCat Crossref Barker P.F. , Dalziel I.Y., Storey B., 1991 . Tectonic development of the Scotia Arc region , in The Geology of Antarctica , pp. 215 – 248 ., ed. Tingey R.J., Clarendon Press . Barletta V.R. , Ferrari C., Diolaiuti G., Carnielli T., Sabadini R., Smiraglia C., 2006 . Glacier shrinkage and modeled uplift of the Alps , Geophys. Res. Lett. , 33 , doi:10.1029/2006GL026490. 10.1029/2006GL026490 OpenURL Placeholder Text WorldCat Crossref Barletta V.R. et al. ., 2018 . Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability , Science , 360 , 1335 – 1339 . 10.1126/science.aao1447 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Berthier E. , Scambos T.A., Shuman C.A., 2012 . Mass loss of Larsen B tributary glaciers (Antarctic Peninsula) unabated since 2002 , Geophys. Res. Lett. , 39 , doi:10.1029/2012GL051755. 10.1029/2012GL051755 OpenURL Placeholder Text WorldCat Crossref Boehm J. , Werl B., Schuh H., 2006 . Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data , J. geophys. Res. , 111 , doi:10.1029/2005JB003629. OpenURL Placeholder Text WorldCat Chen X. , Zhang X., Church J.A., Watson C.S., King M.A., Monselesan D., Legresy B., Harig C., 2017 . The increasing rate of global mean sea-level rise during 1993–2014 , Nat. Clim. Change , 7 , 492−495 . 10.1038/nclimate3325 OpenURL Placeholder Text WorldCat Crossref Dong D. , Yunck T., Heflin M., 2003 . Origin of the international terrestrial reference frame , J. geophys. Res. , 108 , doi:10.1029/2002JB002035. OpenURL Placeholder Text WorldCat Dziewonski A.M. , Anderson D.L., 1981 . Preliminary reference Earth model , Phys. Earth planet. Inter. , 25 , 297 – 356 . 10.1016/0031-9201(81)90046-7 Google Scholar Crossref Search ADS WorldCat Crossref Ezer T. , Atkinson L.P., Corlett W.B., Blanco J.L., 2013 . Gulf Stream's induced sea level rise and variability along the US mid-Atlantic coast , J. geophys. Res. , 118 , 685 – 697 . Google Scholar Crossref Search ADS WorldCat Huang N.E. , Wu Z., 2008 . A review on Hilbert-Huang transform: method and its applications to geophysical studies , Rev. Geophys. , 46 , doi:10.1029/2007RG000228. 10.1029/2007RG000228 OpenURL Placeholder Text WorldCat Crossref Huang N.E. et al. ., 1998 . The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis , Proc. R. Soc. A , 454 ( 1971 ), 903 – 995 . 10.1098/rspa.1998.0193 Google Scholar Crossref Search ADS WorldCat Crossref Ji F. , Wu Z., Huang J., Chassignet E.P., 2014 . Evolution of land surface air temperature trend , Nat. Clim. Change , 4 , 462 – 466 . 10.1038/nclimate2223 Google Scholar Crossref Search ADS WorldCat Crossref King M.A. , Bingham R.J., Moore P., Whitehouse P.L., Bentley M.J., Milne G.A., 2012 . Lower satellite-gravimetry estimates of Antarctic sea-level contribution , Nature , 491 , 586 – 589 . 10.1038/nature11621 Google Scholar Crossref Search ADS PubMed WorldCat Crossref King M.A. , Santamaría-Gómez A., 2016 . Ongoing deformation of Antarctica following recent Great Earthquakes , Geophys. Res. Lett. , 43 , 1918 – 1927 . 10.1002/2016GL067773 Google Scholar Crossref Search ADS WorldCat Crossref King M.A. , Whitehouse P.L., Van Der Wal W., 2016 . Incomplete separability of Antarctic plate rotation from glacial isostatic adjustment deformation within geodetic observations , Geophys. J. Int. , 204 , 324 – 330 . 10.1093/gji/ggv461 Google Scholar Crossref Search ADS WorldCat Crossref Klemann V. , Martinec Z., Ivins E.R., 2008 . Glacial isostasy and plate motion , J. Geodyn. , 46 , 95 – 103 . 10.1016/j.jog.2008.04.005 Google Scholar Crossref Search ADS WorldCat Crossref Lloyd A. et al. ., 2019 . Seismic structure of the antarctic upper mantle based on adjoint tomography , J. geophys. Res. , 125 , doi:10.1029/2019JB017823. OpenURL Placeholder Text WorldCat Nield G.A. et al. ., 2012 . Increased ice loading in the Antarctic Peninsula since the 1850s and its effect on glacial isostatic adjustment , Geophys. Res. Lett. , 39 , 17 . OpenURL Placeholder Text WorldCat Nield G.A. et al. ., 2014 . Rapid bedrock uplift in the Antarctic Peninsula explained by viscoelastic response to recent ice unloading , Earth planet. Sci. Lett. , 397 , 32 – 41 . 10.1016/j.epsl.2014.04.019 Google Scholar Crossref Search ADS WorldCat Crossref Rack W. , Rott H., 2004 . Pattern of retreat and disintegration of the Larsen B ice shelf, Antarctic Peninsula , Ann. Glaciol. , 39 , 505 – 510 . 10.3189/172756404781814005 Google Scholar Crossref Search ADS WorldCat Crossref Rott H. , Skvarca P., Nagler T., 1996 . Rapid collapse of northern Larsen ice shelf, Antarctica , Science , 271 , 788 – 792 . 10.1126/science.271.5250.788 Google Scholar Crossref Search ADS WorldCat Crossref Rott H. et al. ., 2018 . Changing pattern of ice flow and mass balance for glaciers discharging into the Larsen A and B embayments, Antarctic Peninsula, 2011 to 2016 , Cryosphere , 12 , 1273 , doi:10.5194/tc-2017-259. 10.5194/tc-12-1273-2018 Google Scholar Crossref Search ADS WorldCat Crossref Santamaría-Gómez A. , Mémin A., 2015 . Geodetic secular velocity errors due to interannual surface loading deformation , Geophys. J. Int. , 202 , 763 – 767 . 10.1093/gji/ggv190 Google Scholar Crossref Search ADS WorldCat Crossref Scambos T.A. , Bohlander J., Shuman C.U., Skvarca P., 2004 . Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica , Geophys. Res. Lett. , 31 , doi:10.1029/2004GL020670. 10.1029/2004GL020670 OpenURL Placeholder Text WorldCat Crossref Shepherd A. et al. ., 2018 . Mass balance of the Antarctic Ice Sheet from 1992 to 2017 , Nature , 556 , 219 – 222 . Google Scholar Crossref Search ADS PubMed WorldCat Simms A.R. , Ivins E.R., Dewitt R., Kouremenos P., Simkins L.M., 2012 . Timing of the most recent Neoglacial advance and retreat in the South Shetland Islands, Antarctic Peninsula: insights from raised beaches and Holocene uplift rates , Quat. Sci. Rev. , 47 , 41 – 55 . 10.1016/j.quascirev.2012.05.013 Google Scholar Crossref Search ADS WorldCat Crossref Thomas I.D. et al. ., 2011 . Widespread low rates of Antarctic glacial isostatic adjustment revealed by GPS observations , Geophys. Res. Lett. , 38 , doi:10.1029/2011GL049277. 10.1029/2011GL049277 OpenURL Placeholder Text WorldCat Crossref Velicogna I. , Wahr J., 2013 . Time-variable gravity observations of ice sheet mass balance: precision and limitations of the GRACE satellite data , Geophys. Res. Lett. , 40 , 3055 – 3063 . 10.1002/grl.50527 Google Scholar Crossref Search ADS WorldCat Crossref Whitehouse P.L. , 2018 . Glacial isostatic adjustment modelling: historical perspectives, recent advances, and future directions , Earth Surf. Dyn. , 6 , 401 – 429 . 10.5194/esurf-6-401-2018 Google Scholar Crossref Search ADS WorldCat Crossref Wolstencroft M. et al. ., 2015 . Uplift rates from a new high-density GPS network in Palmer Land indicate significant late Holocene ice loss in the southwestern Weddell Sea , Geophys. J. Int. , 203 , 737 – 754 . 10.1093/gji/ggv327 Google Scholar Crossref Search ADS WorldCat Crossref Wu Z. , Huang N.E., 2009 . Ensemble empirical mode decomposition: a noise-assisted data analysis method , Adv. Adapt. Data Anal. , 1 , 1 – 41 . 10.1142/S1793536909000047 Google Scholar Crossref Search ADS WorldCat Crossref Wu Z. , Huang N.E., Long S.R., Peng C.-K., 2007 . On the trend, detrending, and variability of nonlinear and nonstationary time series , Proc. Natl. Acad. Sci. USA , 104 , 14889 – 14894 . 10.1073/pnas.0701020104 Google Scholar Crossref Search ADS WorldCat Crossref Zhao C. , King M.A., Watson C.S., Barletta V.R., Bordoni A., Dell M., Whitehouse P.L., 2017 . Rapid ice unloading in the Fleming Glacier region, southern Antarctic Peninsula, and its effect on bedrock uplift rates , Earth planet. Sci. Lett. , 473 , 164 – 176 . 10.1016/j.epsl.2017.06.002 Google Scholar Crossref Search ADS WorldCat Crossref © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Reduced ice mass loss and three-dimensional viscoelastic deformation in northern Antarctic Peninsula inferred from GPS
JO - Geophysical Journal International
DO - 10.1093/gji/ggaa229
DA - 2020-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/reduced-ice-mass-loss-and-three-dimensional-viscoelastic-deformation-tqDn16lIKk
SP - 1013
VL - 222
IS - 2
DP - DeepDyve
ER -