Incomplete separability of Antarctic plate rotation from glacial isostatic adjustment deformation within geodetic observations

S U M M A R Y Geodetic measurements of Antarctic solid Earth deformation include signals from plate rotation and glacial isostatic adjustment (GIA). Through simulation, we investigate the degree to which these signals are separable within horizontal GPS site velocities that commonly define plate rotation estimates and that promise new constraints on models of GIA. Using a suite of GIA model predictions that incorporate both 1-D and 3-D Earth rheologies, we show that, given the present location of GPS sites within East Antarctica, unmodelled or mismodelled GIA signal within GPS velocities produces biased estimates of plate rotation. When biased plate rotation is removed from the GPS velocities, errors as large as 0.8 mm yr−1 are introduced; a value commonly larger than the predicted GIA signal magnitude. In the absence of reliable forward models of plate rotation or GIA then Antarctic geodetic velocities cannot totally and unambiguously constrain either process, especially GIA.


I N T RO D U C T I O N
Models of Antarctic glacial isostatic adjustment (GIA) show large differences despite recent developments (Whitehouse et al. 2012b;Ivins et al. 2013;Argus et al. 2014).Errors in such models have substantial impact on GRACE-based estimates of ice mass contribution to sea level (King et al. 2012).GPS vertical velocities are now one of the primary sources of information that can be used to assess the reliability of forward models of GIA or empirical solutions of GIA-induced surface deformation, especially in Antarctica where other data are sparse (Whitehouse et al. 2012b;Ivins et al. 2013;Argus et al. 2014;Gunter et al. 2014).Historically, the more precise horizontal GPS component has not been widely used, primarily due to the sensitivity of GIA model predictions of horizontal deformation to commonly unmodelled lateral (3-D) variations in Earth rheology (Kaufmann et al. 2005;Latychev et al. 2005;King et al. 2010).Models that consider these variations are now being developed (e.g.van der Wal et al. 2015) but a second issue that must be addressed before GPS estimates of horizontal deformation can be compared with model predictions is how to determine and remove the component of observed velocity that is due to plate motion.
One approach is to solve for global plate motions using only GPS sites that are known, or presumed, not to be strongly affected by GIA (Sella et al. 2007;Altamimi et al. 2012).The result can be used to correct observations at all GPS sites for the underlying effects of plate motion, and the residual velocities should reflect deformation due to GIA.This approach assumes that the GIA signal, when sampled at the sites used to estimate plate motion, does not include a component that resembles rigid-body rotation.
However, the possibility exists that horizontal motion due to GIA, sampled at discrete GPS locations, does contain a rotation-like signal (Klemann et al. 2008).If this is not accounted for then estimates of plate motion will be biased to some extent and comparisons between residual GPS velocities and GIA model predictions will be flawed because the rigid-body component of GIA would have been removed.A second approach has therefore sometimes been adopted (Klemann et al. 2008;Argus & Peltier 2010) whereby GPS site velocities are corrected for GIA prior to solving for plate rotation.This approach attempts to eliminate the potential bias due to unmodelled GIA that is present in the first approach.However, different bias may be introduced if the GIA model used to correct the GPS velocities is inaccurate; given that the GPS data are being used to validate models of GIA the argument becomes somewhat circular.
While most sites globally are subject to horizontal GIA signals >∼0.5 mm yr −1 (Klemann et al. 2008), the problem of separating tectonic and GIA signals is particularly acute in regions where the GIA signal is significant across large portions of the exposed continental crust, to which GPS antennas are attached, notably North America and especially Antarctica.In this study we therefore investigate the degree to which current methods for removing plate rotation in areas affected by GIA introduce bias into the inferred GIA velocity field because of the imperfectly known rotation-like component of GIA being absorbed into the plate rotation estimate.In particular, we use simulations to test the degree to which we can separate Antarctic GIA and plate rotation.We base these simulations on GIA model output and the existing network of geodetic sites, thereby testing the accuracy to which we can estimate both plate rotation and horizontal velocities associated with GIA.

GIA Modelling
The GIA model accounts for both lateral and radial variations in Earth rheology (van der Wal et al. 2015).The surface load history is determined by solving the sea level equation (Farrell & Clark 1976).This is forced by the W12 Antarctic deglaciation model (Whitehouse et al. 2012a), which is incorporated into the ICE-5G model (Peltier 2004) in order to define the global ice-loading history during the last glacial cycle.
Global solid Earth deformation is calculated using the finite element software package ABAQUS, where elements have a 2 • × 2 • resolution at the surface and viscoelastic parameters are defined for a series of layers with boundaries at 35,70,120,230,400,670,1170 and 2890 km depth. 2 • × 2 • resolution is a coarse resolution to represent details in the coastline and ice extent in, for example, the Antarctic Peninsula.However, the important features in the velocity field around the former ice domes are adequately represented by this model resolution.The model is incompressible, but we include an experiment with Poisson's ratio smaller than 0.5 which simulates 'material compressibility' but neglects the effect of compressibility on buoyancy and self-gravitation (Klemann et al. 2003;Wang et al. 2008).Since these effects are of second order (Klemann et al. 2003), we do not expect them to have a large influence on the variability in horizontal velocities.
The model does not consider deformation due to 'rotational feedback' (Mitrovica et al. 2001), but these effects are long wavelength and would be near-identical across all models as the same ice history and the same lower mantle viscosity are adopted in all our model variants.Therefore this feedback is not expected to have a strong effect on the variability between different models.
In most experiments the lithosphere is represented by assuming purely elastic deformation down to a depth of 120 km, similar to the approach of Latychev et al. (2005).However, similar to this earlier study, we also include one experiment in which we account for the presence of plate boundary zones.In this additional experiment, plate boundary elements in the uppermost layer (0-35 km) are assumed to have rheological properties identical to the 'weakest' (least viscous) element in the next layer down (35-70 km).And in order to capture variations in lithospheric thickness and the broadening of plate boundary zones at depth, rheological properties in this second layer and the layer below (70-120 km) are assigned using seismic velocity perturbations, as described in the next paragraph.
In all experiments, viscoelastic deformation below 120 km is assumed to occur via two modes of deformation: diffusion and dislocation creep.Further details on the parametrization of these modes of creep are given in van der Wal et al. (2010Wal et al. ( , 2013)), but we note here that the main factor contributing to lateral variations in rheology is mantle temperature.In this study we use velocity perturbations from five different global seismic models (see Table 1) to determine vertical and horizontal temperature perturbations using the depth-dependent temperature derivative of seismic wave velocities given in Karato (2008).The scenario in which rheological Table 1.Global seismic velocity models used to define lateral and vertical variations in mantle rheology.All experiments use the W12 ice-loading model (Whitehouse et al. 2012a), and the experiment name (see the main text) also indicates the grain size and water content adopted, for example, W12_SL_dry_4 mm represents the model that uses the Schaeffer & Lebedev (2013)  properties only vary in the radial direction ('1-D' models) is also investigated; in this case a standard mantle geotherm (Turcotte & Schubert 2002) is used to define temperature within the constitutive equation (Hirth & Kohlstedt 2003).For each global seismic model, velocity perturbations are calculated relative to the reference Earth model associated with that model, and we assume that 100 per cent of the seismic perturbations are due to temperature perturbations; hence our predictions represent an upper bound on the influence of lateral variations in rheology on GIA.In addition to defining the temperature distribution within the mantle, the governing equations for diffusion and dislocation creep (Hirth & Kohlstedt 2003) require information about average grain size, water content, and melt content.In our experiments we assume no melt content, the average grain size is taken to be either 4 or 10 mm, and water content is varied between a somewhat wet (1000 ppm H 2 O) and a fully dry state.Details of the full set of experiments carried out are listed in Table 1.
In total we produced 26 sets of modelled horizontal velocities due to GIA.Velocities from each model, interpolated at existing Antarctic GPS site locations, are shown in Supporting Information Figs S1 and S2.Across this suite of Earth models, both the magnitude and the direction of the horizontal velocity prediction at any location can vary substantially, with significant differences being exhibited between models that adopt the same underlying global seismic model but different grain size and water content.These differences are typically larger than the detection limits of modern GPS, at least in the absence of confounding signals.
The 1-D models produce similar results to the models that include lateral Earth structure, but differences are apparent in the Antarctic Peninsula and along the Transantarctic Mountains (TAM), suggesting that the effects of low mantle viscosity (e.g.Antarctic Peninsula) and strong viscosity gradients (e.g.perpendicular to the TAM) both warrant further investigation.

Plate rotation estimation from simulated data
To examine the degree to which plate rotation estimates are affected by GIA we simulate an 'observed' horizontal velocity field from the output of each GIA model run.We do not introduce any source of actual plate rotation.For each 'observed' horizontal velocity field we estimate a plate rotation parametrized as an Euler pole location and rate of rotation.Once the plate rotation bias is estimated we evaluate it in terms of GPS velocity at all GPS locations, including those not used in the plate rotation estimation.
Our simulations consider two different strategies which differ in terms of whether a different GIA model is first removed from the simulated observations or not, as described below.

No correction for GIA (strategy 1)
In this first strategy, the simulated observations are not corrected for GIA in any way prior to estimating plate rotation, a strategy adopted in several previous global (Kreemer et al. 2014) or regional (Lidberg et al. 2007;Sella et al. 2007) studies.Velocities are extracted from the GIA model runs at the locations of GPS sites commonly used to define the Antarctic Plate; we adopt the sites used by Argus et al. (2014) (green circles in Fig. 1).We discuss variants to this site selection below.
For this set of velocities we then estimate the location and rotation rate of the Euler pole using standard techniques (Goudarzi et al. 2014).The Euler pole parameters are sensitive to the uncertainties assigned to each plate-defining site and we adopt the site-specific uncertainties of Argus et al. (2014), assuming that their north and east velocity uncertainties equally contribute to their tabulated speed uncertainty.The plate rotation realization is dominated by the longrunning, high-precision sites in coastal East Antarctica.
We obtain one set of Euler pole parameters per GIA model, and use these to predict plate-rotation velocities at all GPS site locations (pink lines in Fig. 1a), including those not used to define the rigid plate.Bearing in mind that, by design, our simulated velocity fields do not include any component of plate motion, any non-zero platerotation velocity is the degree to which GPS velocities, after removal of plate rotation, would be in error because part of the GIA signal is removed.

Correction for GIA using model predictions (strategy 2)
Some previous studies have recognized the potential importance of accounting for GIA when estimating plate rotation, and have attempted to remove the effects of GIA using model predictions (of poorly known accuracy) before estimating Euler pole parameters (e.g.Argus & Peltier 2010).In this second strategy we adopt output from one GIA model as the simulated 'observed' velocity field as before, and output from a different GIA model as an analyst's chosen GIA model which is subtracted from the observations prior to solving for the Euler pole parameters.In other words, differences in pairs of GIA model predictions are used in the Euler pole computations.

R E S U LT S
Considering the simulations where no attempt is made to correct for GIA (strategy 1) we find that the effect of the modelled GIA signal on the Euler pole parameters is non-zero.An example of the bias introduced to the Euler pole parameters by the GIA signal, in terms of site velocities, is shown in Fig. 1(a) where we arbitrarily adopt the simulated observations given by GIA model W12_1D_dry_4 mm (see caption to Table 1).Here, the brown lines show the original horizontal velocity field predicted by the W12_1D_dry_4 mm experiment, that is our hypothetical 'observed' velocity field, the pink lines show the velocity bias due to biased plate rotation, while the green lines show the velocity field after the subtraction of the plate rotation.With the model adopted for this example, the effect of removing the biased plate rotation is to alter both the magnitude and direction of the simulated site velocities.Velocities at the East Antarctic reference sites are significantly reduced through the estimation and subtraction of plate rotation; this is expected because the Euler pole estimation attempts to minimize the rotational signal at these sites.In West Antarctica, the plate-rotation-corrected velocities (green) shift in direction by ∼10 • -40 • and are generally increased in magnitude, compared to the unmodified velocities (brown).
Next we consider the simulations where a correction for GIA is applied to the simulated observations prior to solving for plate rotation (strategy 2).A bias will be introduced by this method if the GIA correction that is applied does not accurately represent 'true' GIA.In Fig. 1(c), the line colours are as for Fig. 1(a), with the 'observed' field defined by output from the W12_1D_dry_4 mm experiment.Prior to estimating plate rotation, a correction for GIA is made using output from an arbitrarily chosen alternative experiment; here we chose the W12_SAW_dry_10 mm experiment (Supporting Information Fig. S1).The site velocity biases using this pair of models (i.e. the biased plate rotation velocities, pink lines, Fig. 1c) show a different spatial pattern of bias to that in Fig. 1(a) but with a similar magnitude.Indeed, the sense of the biased plate rotation is opposite, highlighting that subtracting a GIA model of unknown accuracy from GPS velocities does not guarantee a more accurate estimate of plate rotation than not subtracting a model at all (Fig. 1a).
Comparing the final velocities (green) in Figs 1(a) and (c) reveals velocities that are different by as much as 90 • in direction in the Antarctic Peninsula and magnitudes that vary routinely by 50 per cent.This result suggests that using GPS velocities which are affected by GIA, or residual GIA, to estimate and remove plate rotation may lead to erroneous conclusions being drawn about the validity of GIA model predictions.
To examine a more complete set of possibilities we extended the results shown in Fig. 1 to consider site biases based on all our GIA model predictions.The velocity biases shown in Fig. 2 are based on all models (Fig. 2a) or all models differenced from the W12_1D_dry_4 mm model (Fig. 2b).
The site velocity biases vary between models, with maximum biases per model, or model pair, lying in the range 0.11 to 0.84 mm yr −1 , considering both strategies.Across all models the mean maximum biases are 0.46 and 0.36 mm yr −1 for strategies 1 and 2, respectively.That is, one can reasonably expect that GPS horizontal velocities are biased, after the removal of a GPS-determined plate rotation model, at these levels.Given that the mean predicted horizontal speed across all our GIA models is 0.46 mm yr −1 , GIA-biased plate rotation represents a substantial potential source of error.Indeed, the maximum bias in Fig. 2 is larger than 90 per cent of all modelled speeds across all our GIA models.
Examining the results from our models with 1-D Earth structure shows site velocity biases that are within the range of those associated with 3-D Earth models.Models that include plate boundaries and lithospheric thickness variations did not perform systematically differently to models that do not include these effects.Interestingly, the direction of the biases in Fig. 2 can be broadly partitioned into GIA models with dry or wet rheology, as shown in orange or blue, respectively.Differences arise as a result of wet rheologies being generally weaker and reflecting more recent and more local changes in ice loading than models with dry rheologies.In the model scenarios considered here, dry rheologies tend to predict southerly motion of East Antarctica whereas wet rheologies predict more northerly motion (compare Supporting Information Figs S1 and S2).
Adopting a different set of reference stations for the Euler pole estimation changed the spatial pattern of the biases but did not appreciably mitigate their magnitude.For instance, adopting just longrunning GPS International GNSS Service (IGS; Dow et al. 2009) sites across all of Antarctica (Altamimi et al. 2012), changed the maximum bias to between 0.54 mm yr −1 (strategy 1) and 0.60 mm yr −1 (strategy 2).We note that very recent ice load changes are known to produce large site motions at IGS sites (O'Higgins, Palmer and Rothera) in the Antarctic Peninsula (see, e.g.Thomas et al. 2011;Nield et al. 2014); such recent load changes are not considered in W12 and hence the test that includes these sites for Euler pole estimation may be biased by this omission.These motions, plus the observation of other localized motion at O'Higgins (Argus et al. 2011), leads us to recommend that sites in this region are not used to define plate rotations.

D I S C U S S I O N
Our simulations suggest that imperfectly known GIA has the potential to bias geodetically derived Euler pole estimates for the Antarctic Plate; our estimates suggest that the magnitude of bias spans 5 • in longitude, 0.7 • in latitude and 0.007 • Ma −1 in rotation rate.This bias is solely due to unmodelled or mismodelled GIA, something that cannot presently be removed from geodetic observations with confidence.These biases are quite small relative to the total plate rotation signal (Figs 1b and d) but are often large relative to horizontal GIA (Figs 1a and c).
Our results also have important implications for studies attempting to use horizontal GPS velocities to validate or tune models of GIA.In the absence of an error-free and independent model of plate rotation, this signal cannot be removed from GPS velocities with complete confidence.For instance, it is not clear that using horizontal GPS velocities to validate ICE-6G_C in Antarctica (Argus et al. 2014) is a robust test because such a test is predicated on the rotation-like component of horizontal GIA, as expressed at select East Antarctic GPS sites, being negligibly small.Argus et al. (2014) found the land surrounding the Ronne Ice Shelf to be moving radially outward from an inferred former ice centre.They show this pattern to be in agreement with predictions derived using a model of GIA that just considers radially varying Earth properties, but the comparison relies on the assumption that unmodelled horizontal GIA does not substantially affect estimated plate motion, the topic of this paper.
Our tests suggest that the application of GPS horizontal velocities to the problem of validating or tuning Antarctic GIA models must be done with great caution.One robust approach is to compare GPS velocities with GIA model velocities after both their rotational components are removed, but only when the GIA rotational component is computed using the same sites as those used to remove the rotational component from the GPS velocities.The only other truly unambiguous approaches to this problem are to either compare modelled and observed strain rates, as discussed by Marotta et al. (2004), which are invariant to rotation, or to consider small regions over which variations in the expression of plate rotation errors are small.
We focus on just one ice history in this study, combined with a large number of Earth rheologies, both 1-D and 3-D.Adopting a different ice history would only substantially affect our results if the predicted horizontal velocities contained no apparent rigidbody rotation.Similarly, several aspects of our modelling approach could be improved; the coarse spatial resolution of the finite element model, the neglect of rotational feedback, the simplistic treatment of compressibility, and the assumption that temperature variations dominate seismic velocity perturbations.However, addressing these factors is again unlikely to affect our main conclusions.Ultimately it is not known what rotational component exists in the real GIA field of Antarctica and hence caution must be exercised in the interpretation of horizontal GPS velocities to validate or constrain GIA models.

C O N C L U S I O N S
We use a suite of three-dimensional and one-dimensional GIA models to predict horizontal velocities of Antarctica, and show that these vary substantially in magnitude and direction according to the assumed mantle grain size and water content.Using these plausible horizontal GIA velocities as simulated observations, we show that GIA can introduce small biases in geodetic estimates of rigid Antarctic Plate rotation.When such biased plate rotation estimates are subtracted from observed GPS velocities they are propagated across all of Antarctica.These biases routinely reach 0.48 mm yr −1 and, in the worst case scenarios, 0.84 mm yr −1 in West Antarctica and the Antarctic Peninsula.This maximum bias is larger than 90 per cent of the predicted GIA speeds across all models, and hence this is a non-negligible potential error.In the absence of accurate models of either GIA or plate rotation then GPS measurements of horizontal velocity cannot be confidently used to validate models of GIA, unless the rotational component of the modelled GIA is removed using the same sites as used in the plate rotation estimate.Considering modelled and observed strain rates, or comparing velocities over small regions, are other robust ways to undertake this comparison.
While the GIA signal is typically small relative to plate motion, its incomplete treatment in estimating plate rotations may introduce small but possibly important biases in these estimates when horizontal deformation due to GIA extends over a large portion of the landmass and contains a rotation-like component given the geodetic station distribution.This occurs notably in North America, Eurasia and Antarctica.While studies have recognized the potential impact of GIA in these regions and attempted to mitigate it in various ways (e.g.Lidberg et al. 2007;Argus & Peltier 2010;Blewitt et al. 2013;Kierulf et al. 2014), more extensive simulations that consider a full suite of ice histories and 3-D Earth rheologies are required to more fully quantify the potential impact of unmodelled or mismodelled GIA.

A C K N O W L E D G E M E N T S
The research was supported by an NERC Independent Research Fellowship to PLW (grant number NE/K009958/1) and a University of Tasmania Visiting Scholarship.MAK is a recipient of an Australian Research Council Future Fellowship (project number FT110100207).We thank Mohammad Ali Goudarzi for making his plate rotation estimation code available.Model output as shown in Supporting Information Figs S1 and S2 is available from PLW on request.Donald Argus and Jerry Mitrovica provided helpful reviews that improved the clarity of the manuscript.

Figure 1 .
Figure 1.Illustrations of the effect of estimating and removing plate rotations from simulated Antarctic GPS data containing (a) unmodelled or (c) mismodelled horizontal deformation due to GIA.Reference frame sites used in Euler pole estimates are shown as green dots (Argus et al. 2014) and a spatially distributed subset of other sites is shown as red dots.Site velocities shown are those generated using the GIA model (brown), those due to an evaluation of biased Euler pole parameters (pink), and the difference of the two (green).Panels (b) and (d) show observed GPS velocities and 2σ uncertainties in black at selected sites (Argus et al. 2014; D. Argus, personal communication, 2014) alongside the biases repeated from panels (a) and (c) (green) at a different scale.Plate rotations in panels (a) and (b) were estimated using simulated 'observations' from W12_1D_dry_4 mm, while in panels (c) and (d) the simulated 'observations' from W12_1D_dry_4 mm were corrected for GIA using W12_SAW_dry_10 mm, meaning their difference is propagated into the plate rotation estimates.For clarity, not all velocities are shown at the reference frame sites.

Figure 2 .
Figure 2. Site velocity biases using the full suite of GIA models (equivalent to the pink lines in Fig. 1).Bias distribution for six representative sites is shown.Line colours are for models with dry (orange) and wet (blue) rheologies.Panels (a) and (b) show the biases where GIA is unmodelled or mismodelled, respectively.
seismic model, 4 mm grain size and zero water content.Each model listed in Table1is combined with either a wet or a dry rheology, and 4 mm or 10 mm grain size.