Postseismic Deformation in the Northern Antarctic Peninsula Following the 2003 and 2013 Scotia Sea Earthquakes

Large earthquakes in the vicinity of Antarctica have the potential to cause postseismic viscoelastic deformation affecting measurements of displacement that are used to constrain models of glacial isostatic adjustment (GIA). In November 2013, a Mw 7.7 strike‐slip earthquake occurred in the Scotia Sea, 650 km from the Antarctic Peninsula. GPS time series from the northern Peninsula show a change in rate after this event, indicating a far‐field postseismic deformation signal is present. In this study, we use a finite element model with a suite of 1D and 3D Earth structures to investigate the extent of postseismic deformation in the Antarctic Peninsula. Model output is compared with GPS time series to place constraints on the Earth structure in this region. The preferred Earth structure has a thin lithosphere combined with a Burgers rheology with steady‐state viscosity of 4 × 1018 Pa s and transient viscosity one order of magnitude lower. Our study shows that including 3D Earth structure does not improve the fit. Using the best fitting Earth structure, we run a forward model of the nearby 2003 Mw 7.6 strike‐slip earthquake and combine the predictions for both earthquakes. We show that postseismic deformation is widespread across the northern Peninsula with rates of horizontal deformation up to 1.65 mm/yr for the period 2015–2020, a signal that persists for decades. These results suggest that much of Antarctica may be deforming due to recent postseismic deformation and this signal needs to be accounted for when using GPS observations to constrain geophysical models.


10.1029/2023JB026685
2 of 19 importance of large, far-field earthquakes to Antarctic deformation will depend to some extent on the rheological properties of the Antarctic lithosphere and upper mantle, notably when considering sustained postseismic deformation that can persist for decades and longer (e.g., Huang et al., 2020;Khazaradze et al., 2002;Klotz et al., 2006).
The continent of Antarctica is made up of two geologically distinct regions with a thick cratonic lithosphere and high viscosity mantle in the East, and a lower viscosity active rift system in the West (Morelli & Danesi, 2004;Wiens et al., 2023).Analysis of gravity data has shown the effective elastic thickness of the lithosphere to be thinner in West than in East Antarctica (Pappa & Ebbing, 2023;Paxman, 2023).Recent seismic tomographic models (Lloyd et al., 2020;O'Donnell et al., 2019) have also highlighted strong lateral variations in the upper mantle structure across Antarctica with notable slower than average shear waves speeds for the Antarctic Peninsula, indicating lower viscosity there relative to the rest of Antarctica.Studies of GIA have also inferred a low upper mantle viscosity for the Antarctic Peninsula by using GPS to constrain modeled viscoelastic uplift in response to decadal scale ice-mass loss (Nield et al., 2014;Samrat et al., 2020Samrat et al., , 2021;;Zhao et al., 2017).The way in which this low viscosity zone is affected by far-field external forcing is not yet clear.
Non-negligible postseismic deformation in Antarctica not only depends on rheological properties but also requires earthquake magnitude to be sufficiently large.For very large earthquakes (greater than magnitude 9) postseismic deformation can be observed over 1,000 km from the source (Shao et al., 2016) and over very long time scales, deformation from the magnitude 9.5 earthquake that occurred in Chile in 1960 was observed 40 years later (Khazaradze et al., 2002).Postseismic deformation could also be mis-interpreted as inter-seismic deformation when geodetic records are short or information on historic earthquakes is limited (Trubienko et al., 2013).
The study of postseismic deformation in Antarctica lags far behind other continents such as Asia and North America mainly due to the, until now, limited geodetic stations in Antarctica and the low frequency of large earthquakes surrounding it.To date, the only study published on Antarctic postseismic deformation (King & Santamaria-Gomez, 2016) demonstrates that Antarctica has been undergoing postseismic viscoelastic deformation since the 1998 magnitude 8.2 Balleny Islands earthquake, at the time the largest intraplate earthquake recorded globally (Figure 1).This earthquake occurred 600 km from the Antarctic coastline (Figure 1) but was a large enough magnitude that it caused measured horizontal deformation on the continent of up to 4.8 mm/yr (King & Santamaria-Gomez, 2016).The modeling presented by King and Santamaria-Gomez (2016) suggests this deformation is ongoing and widespread across the continent.For a full summary on the topic in an Antarctic context, the reader is referred to van der Wal et al. (2023).(Bird, 2003).Focal mechanisms of earthquakes mentioned in the text shown in black and all post-1976 earthquakes with M w ≥ 7.0 shown in red (from the Global CMT Catalog (Ekström et al., 2012)).

10.1029/2023JB026685
3 of 19 In this study we investigate postseismic deformation in the Antarctic Peninsula in response to two large earthquakes which occurred approximately 650 km north-east of the tip of the Peninsula (Figure 1) at the Scotia Ridge-the left-lateral strike-slip transform boundary between the Scotia Plate and the Antarctic Plate.These large strike-slip earthquakes occurred on 4 August 2003 (M w 7.6) and 17 November 2013 (M w 7.7) just to the east and north of the South Orkney Islands, respectively.Several continuous GPS stations were in operation in the northern Antarctic Peninsula at the time of the 2013 earthquake and close examination of these time series reveals a change in displacement rate following the earthquake.This is the first time GPS data from the northern Peninsula has been used to investigate postseismic deformation.We use the GPS-observed change in rate to place constraints on the regional Earth structure by testing the fit of the modeled postseismic deformation against the GPS time series, and compare our findings to the regional Earth structure inferred from GIA studies.Due to the lack of GPS stations in operation at the time of the 2003 earthquake it is not possible to do the same for this earthquake.Instead, we use the inferred range of Earth structures constrained from the 2013 earthquake to forward model the postseismic deformation in response to the 2003 earthquake to estimate the total potential magnitude and spatial pattern of postseismic deformation across the northern Antarctic Peninsula.

GPS Data
There are 18 GPS stations located on the Antarctic Peninsula and South Orkney Islands that lie within 1,200 km of the 2003 and 2013 Scotia Sea earthquakes (locations and details given in Table 1, Table S1 in Supporting Information S1 and subset shown in Figure 2).The stations have been installed at various times from 1995 onwards as part of several different projects, and at the time of writing, around half remain operational.Time series records from some of the stations cease prior to the 2003 or 2013 earthquakes or contain data gaps at key times precluding them from being included in the analysis in this study.We exclude any sites more than 1,200 km from the earthquakes as it is unlikely that deformation observed at these distances would be dominated by postseismic deformation and, on the southern Antarctic Peninsula, observations are more sensitive to local effects of changing ice loads.
The time series are used in two ways.First, where data exists at the time of the 2013 earthquake, the coseismic offset is estimated from the processed time series and used to refine coseismic fault slip that is input to the model (Section 2.2 and 3.2).Second, the change in displacements after the 2013 earthquake is estimated from the time series (Section 2.3) and used to constrain a forward model of postseismic deformation where the input Earth structure and properties are varied, and the predicted deformation is tested against the detrended time series.

GPS Processing
The GPS data were processed and uncertainties estimated with the GAMIT/GLOBK 10.71 software package (Herring et al., 2016) using the technique described by Koulali et al. (2022).We used daily GPS raw phase observations to estimate a loosely constrained solution of station coordinates, the zenith delay of the atmosphere at each station, and orbital and Earth orientation parameters (EOP).This solution was then used as quasi-observations in a Kalman filter to estimate a consistent set of coordinates and velocities.We examined all the position time series for outliers and offsets or "jumps."We aligned our solution to the ITRF2014 reference frame using the six-parameter transformations (three-rotations, three-translations) of daily coordinates based on a globally-distributed set of ITRF sites.Our GPS solution includes data from the POLENET-ANET and UKANET networks as well as the core IGS stations for the period between 1995 and 2020.

Coseismic Offset
Coseismic offsets were estimated from the processed GPS time series at the time of the 2013 earthquake for those stations that were operational during the event.We manually inspected all the GPS time series for coseismic offsets associated with the 2013 event.Postseismic deformation in GPS time series is often modeled using a logarithmic function (ln(1 + Δt/τ)) within a time series trajectory model.This requires an empirical determination of the relaxation constant τ, which depends on the size of the earthquake, and the choice of this can affect the determined offsets.To avoid this trade-off, we opted to estimate the coseismic offsets using a window of 108 days covering the rupture date using a model with a linear rate and a step function (Montillet et al., 2015).The 108-day window is a compromise in minimizing biases caused by seasonal and sub-seasonal variability, postseismic deformation and time-correlated noise.Offset uncertainties were determined assuming a First-Order Gauss-Markov Extrapolation (FOGMEX) algorithm (Floyd & Herring, 2020).The magnitude of the GPS-estimated coseismic displacement and associated 95% uncertainty at each GPS location is shown in Figure 2 (blue arrows).We were unable to estimate the coseismic offset from the 2003 earthquake as very few GPS stations were operational at this time and, for those that were, there are not enough observations to give a robust estimate.

Postseismic Displacement
To isolate the transient postseismic deformation that is not attributable to the long-term linear rate nor the seasonal deformation, we fitted time series using a trajectory model using the following expression: where A is a constant, B is the secular velocity, C and D are sine and cosine terms to represent the annual cycle, E is the coseismic offset, λ is amplitude of the logarithmic term, τ is the logarithmic constant parameter and Θ(t eq ) is a Heavy-side step function at the time of the earthquake (t eq ).In this model we consider two logarithmic terms corresponding to the 2003 and 2013 events.The tau parameter was set to 490 days for the 2013 earthquake.Then, we correct each time series for the background linear rate (B) and the seasonal terms (C and D).We used the linear least squares inversion method to determine the six constants, including the postseismic amplitude, using a fixed relaxation parameter tau of 490 days for the 2013 earthquake.Fixing the relaxation parameter enables the estimation of the amplitude using a linear least-squares approach.Unlike in Section 2.2, where we only needed to estimate the coseismic offsets to validate the slip model, here we considered fitting the entire observed time series to capture the observed relaxation's time dependence and therefore isolate the postseismic signal.
For GPS sites that recorded observations before and after the 2013 earthquake, the time series were detrended to remove the secular trend and reveal the change in deformation.The length of the time series prior to the earthquake clearly impacts our ability to robustly estimate this change in deformation, with shorter time series giving larger uncertainties on the long-term trend.For each site, the length of the time series prior to the 2013 earthquake is given in Table 1.By removing the long-term trend from the time series, we aim to isolate as much as possible the displacement that is due to postseismic deformation, however we acknowledge that there may be other signals that remain uncorrected from the time series such as deformation related to ice-mass change (see Section 2.4).

Deformation Due To Ice-Mass Change
Recent ice-mass loss in the northern Antarctic Peninsula caused by the disintegration of the Larsen A (in 1995) and Larsen B (in 2002) ice shelves and subsequent acceleration of glaciers (Rignot et al., 2004;Rott et al., 1996;Scambos et al., 2004) is causing ongoing viscoelastic deformation of the Earth on a decadal timescale (Nield et al., 2014;Samrat et al., 2020).The viscoelastic deformation related to ice loss is predominantly vertical, but there is also some smaller horizontal motion as well.This signal may not be fully removed by detrending the time series if the deformation is non-linear through time since the detrending process assumes that a change in deformation rate since 2013 is solely due to postseismic deformation.
The recent study by Samrat et al. (2020) used ice surface elevation change data as ice loading input to a high resolution 1D GIA model to predict viscoelastic deformation in the Northern Peninsula between 1995 and 2018.They tested model-predictions from a suite of 345 Earth structures against vertical and horizontal GPS time series from from locations to constrain a lithospheric thickness and upper mantle viscosity.Importantly, the authors showed that the uplift rate at PALM has reduced since around 2011 indicating a reduction in ice-mass loss.Therefore, there may be some bias when using the detrended GPS to constrain the model of postseismic deformation by not accounting for this signal.
In this paper, we present results of best fitting models when comparing model output to GPS time series.To explore the potential bias, we make an additional correction to the GPS time series to try and account for ice-loss related viscoelastic deformation.We present results both with and without this correction.
The correction is made by calculating the model-predicted ice-related viscoelastic deformation rate at each GPS location before and after the 2013 earthquake from the model predictions presented by Samrat et al. (2020).The change in deformation rate, representing non-linearity in deformation related to ice-mass change, is then removed from the detrended GPS time series for the period following the earthquake.A change in rate of zero, for example, would indicate that ice-loss related deformation over this time period was linear and therefore no additional correction is made.
We test several variations of this correction based on different "best-fitting" viscoelastic models from Samrat et al. (2020).We disregard the north component in (3.) and (4.) as Samrat et al. (2020) showed that there was little sensitivity to changing Earth properties, likely due to the position of the GPS sites in relation to the ice-mass loss locations.Results are shown in Section 4 using no correction to the GPS time series as well as corrections 1-4 detailed above.

Fit to the Data
As detailed in the previous section we make model predictions of postseismic deformation from the 2013 earthquake using a range of different Earth structures.To place bounds on the Earth structures that best fit the data we calculate the WRMS error (Equation 2) between the model predictions and GPS observations as follows: where GPS i and MODEL i are the displacements (in mm) at a given time, and w i is the weighting derived from the GPS uncertainty ( ).We use the time series of displacements rather than rates because estimating a rate would introduce further uncertainties and rates may change through time.Before calculating the WRMS error however, we must first make some decisions about which GPS data to use, and the following considerations are necessary: 1. which direction of motion to use: north, east, and/or vertical; 2. whether or not to correct GPS observations for ice loss-related deformation, and if so, which correction to use (see Section 2.4); 3. which GPS stations to use.
For (1.), we use the north and east components of motion only.This is because the earthquakes are both strike-slip earthquakes meaning the amount of vertical displacement will be minimal.Furthermore, the vertical GPS time series are much more affected by ice-loss related deformation and would therefore be more reliant on making an accurate correction for this.
For consideration of point (2.), we show results both with and without making a correction for the ongoing deformation related to recent ice-mass loss.We discuss the effects of this correction in Section 5.
When considering (3.), we apply some selection criteria to derive a subset of GPS data with which to test model fit.There are 13 GPS sites with data both before and after the 2013 earthquake (locations shown in Figure 2 and the time series are shown in Figures S1 to S13 in Supporting Information S1).We exclude all sites within 100 km of the earthquake (BORC) as the finite-element model is not able to accurately reproduce displacements very close to the earthquake (Nield et al., 2022).We exclude sites that have less than 1 year of data before the earthquake (MBIO and SPGT) as estimating a deformation rate before the earthquake for the purposes of detrending the time series will be unreliable.We also exclude ROBN because it shows an anomalous signal.We further exclude the time series PALM and PAL2 from Palmer station as comparison to nearby sites suggests it may be responding in part to some local source of deformation, and one of the time series from O'Higgins so as not to over-rely on this geographical location (OHI2).This provides a set of GPS time series (CAPF, FONP, VNAD, HUGO, DUPT, OHI3) to test the data.

Finite Element Model
The model used in this study is a global, spherical, finite element model of postseismic deformation (Nield et al., 2022), built using the software package Abaqus (version 2018, Hibbitt et al., 2016).The model is briefly described in this section, and full details of the model, including the results of benchmarking tests, can be found in the study by Nield et al. (2022).
The main purpose of the model is to estimate far-field postseismic deformation, that is, more than approximately 300 km from the source of the earthquake.The model consists of a global, spherical, low-resolution mesh, with a high-resolution region of interest embedded within it.The low-resolution mesh has elements that are up to 500 km in plan view, and the high-resolution mesh has elements that are 10 km at the fault plane, increasing to 50 km at the edge of this region.The depth resolution of the elements for each region is given in Table 2.The two meshes are tied together so that they cannot move independently of each other, and that displacement and stress are continuous through the boundary.A fault plane is defined by creating a cut within the high-resolution mesh and coseismic slip is prescribed at the node pairs that lie on either side of fault plane using kinematic constraint equations (Masterlark, 2003;Nield et al., 2022).Due to difficulties in creating a valid mesh in Abaqus, the geometry of the fault is restricted to a single plane, and multiple planes of different strike and dip cannot yet be accommodated.However, far-field postseismic deformation is less sensitive to simplifications made to the fault geometry than near-field deformation (Tregoning et al., 2013), and simplifying fault plane geometry, while retaining earthquake moment, has negligible effect on far-field estimates (Takeuchi & Fialko, 2013).Although the geometry is restricted to a single fault plane, slip and rake can be varied, and a different value specified at each node on the fault plane if required.A separate model is constructed for each of the 2003 and 2013 earthquakes since the geometry of the slip planes is different.
Following the application of coseismic slip in the first (elastic) step of the model run, the model was run for seven 1-year timesteps to allow viscoelastic deformation, during which time the fault is fixed so that there is no further relative displacement between the node pairs.We use 1-year time steps due to the heavy computational cost of the model.Viscoelastic deformation is governed by the stress resulting from the imposed fault slip and the specified Earth properties.The Earth models used in this study are described in Section 3.3.
Boundary conditions are applied to the model center to ensure it is fixed in space.Following Wu (2004), we apply elastic foundations at each material layer boundary where a difference in density occurs across it (including the surface).The model is also capable of self-gravitation, whereby the change in gravitational potential due to vertical displacement is applied iteratively as an additional load.However, due to the increased processing time we chose not to implement this feature since the effects are negligible for strike-slip earthquakes (Nield et al., 2022).
Since the focus of the study is far-field, the model does not include afterslip (e.g., Pollitz et al., 2021) or poroelastic rebound (e.g., Barbot & Fialko, 2010;Peña et al., 2022) as these effects are considered to be small at distances greater than ∼300 km (Peña et al., 2019).The outputs from the model are coseismic displacement immediately following the application of the slip (i.e., the earthquake), and postseismic displacement at each time step after the earthquake.

Fault Properties
Information about earthquake properties such as fault geometry and slip properties are required inputs to the model.For the 2013 Scotia Sea earthquake we use the finite-fault inversion by Ye et al. (2014) as a basis for our input.The slip model they infer consists of three separate planes each with a different strike and dip (their Figure 7).Because our model is limited to a single fault plane, we take the weighted average of the strike and dip (weighted by the length of each plane) and the total length and width which we round to the nearest 10 km so it can be accommodated by our mesh.We then divide the fault plane into three sections, the length of each matching those deduced by Ye et al. (2014), and assign each section a single slip and rake using an average of their model.
The Ye et al. (2014) solution used only seismic data for the inversion, so we modified the slip and rake to minimize the residuals between the model-predicted coseismic displacement and calculated geodetic GPS offsets.This was necessary to compensate for the simplifications made to the fault geometry.We verify that our adaptation of this solution for our finite element model is suitable by calculating the seismic moment (9.5 × 10 20 Nm) and comparing to the value published by Ye et al. (2014) (1.1 × 10 21 Nm).A small difference of 13.5% in the seismic moment would not significantly affect our far-field results (e.g., King & Santamaria-Gomez, 2016).The GPS and model-predicted coseismic displacement for the 2013 earthquake is shown in Figure 2a and resulting fault plane is shown in Figure 2c.A good match between GPS-observed and model-predicted coseismic displacement (modeled vectors within the GPS uncertainty ellipsoids) is obtained at most sites apart from PALM/PAL2 and HUGO which may be responding to some unaccounted-for local loading or unloading causing greater northwards motion.
For the 2003 earthquake, we are not able to use GPS estimated offsets to tune slip properties as there are not enough observations from the time of this earthquake to estimate coseismic offset.Instead, we take the fault geometry, strike, dip, and rake from the initial model of Ye et al. ( 2014

Earth Models and Rheology
The finite element model offers flexibility in terms of the Earth structure and rheological models that can be implemented.Maxwell and Burgers rheology are implemented through the use of a Prony time series directly in Abaqus (e.g., Nield et al., 2022;Takeuchi & Fialko, 2013), whereas more complex 3D or power-law rheologies can be implemented through a user-defined subroutine which provides properties on an element-by-element basis (e.g., Freed et al., 2012;van der Wal et al., 2015).
In this study, we aim to explore a wide parameter space to place constraints on the possible Earth properties in the northern Antarctic Peninsula for earthquake related deformation over timescales of years.Due to the computation time of the model, it is impossible to test every Earth structure and rheological model.However, we select a set of candidate Earth structures and vary parameters within each model and suggest a best fit model from those tested.For all models, we keep the same layer structure in terms of density and elastic parameters derived from the Preliminary Reference Earth Model (PREM, Dziewonski & Anderson, 1981) (Table 2), but we vary the viscosity assigned to each layer.Upper layers that are modeled as the lithosphere are set as a purely elastic material in Abaqus.The lower mantle is too deep to influence the results since the postseismic response to the earthquake occurs at shallow depths and hence the lower mantle retains the same properties for all 1D and 3D models, with a linear Maxwell viscosity of 5 × 10 21 Pa s.The Earth models detailed here relate to the parameter space tested for the 2013 earthquake only.The 2003 earthquake is modeled with the resulting "best fit" Earth model.
We begin by testing simple 1D Earth structures with Maxwell and Burgers rheology, both of which have been extensively used to model postseismic deformation (Pollitz et al., 2000;Takeuchi & Fialko, 2013;Wang et al., 2012).We vary the thickness of the lithosphere, the asthenosphere and the upper mantle, as well as varying the viscosity of the asthenosphere and upper mantle, testing a total of 267 different configurations.The upper mantle always has a viscosity the same or higher than that of the asthenosphere.For the Burgers models, the asthenosphere and upper mantle are assigned a transient and steady-state viscosity, where the transient viscosity is one order of magnitude less than the steady state, as commonly used in postseismic deformation modeling (e.g., Hearn et al., 2009;King & Santamaria-Gomez, 2016;Trubienko et al., 2013).We keep the ratio of transient to steady state viscosity fixed for all Burgers models to prevent the parameter space becoming prohibitively large.
We test a total of 164 Earth models with Burgers rheology.Elastic parameters remain the same for transient and steady-state materials.All parameters and ranges tested are given in Table 2.
In addition to 1D Earth models we also test two different 3D models to see whether they can give a better fit to the GPS-observed postseismic deformation.
There is some evidence from seismic and geophysical studies of lateral heterogeneity in Earth structure in the Northern Antarctic Peninsula, and in particular, the presence of the subducted former Phoenix Plate below (Lloyd et al., 2020;Parera-Portell et al., 2021;Yegorova et al., 2011).Lloyd et al. (2020) attribute fast anomalies (i.e., probable higher viscosity) at 300-700 km depth beneath the Peninsula to being part of the subducted slab.The first 3D (linear Maxwell) model simulates this structure with a 50 km thick high viscosity slab underlying a low viscosity wedge (Figure 3).A suite of models (72 in total) is run with two wedge viscosities, and an asthenosphere of varying thickness and viscosity.
The second 3D model uses a composite power-law rheology for the uppermost 400 km of the Earth.The main aspects of the 3D model implementation are described below and a full description of the methods can be found in van der Wal et al. (2015) and van der Wal et al. (2013).Instead of directly prescribing mantle viscosity, as is the case for the Maxwell and Burgers rheologies, creep properties are calculated for each element and fed into the model on an element-by-element basis.The creep strain components in Abaqus are then computed as follows: where q is the von Mises stress, t is time, n is the stress exponent (here we use n = 3.5), and B diff and B disl are the creep parameters for diffusion and dislocation creep, respectively.The material behavior becomes non-linear (i.e., power-law creep) when n is greater than 1.The creep parameters are calculated according to laboratory-derived flow laws for olivine (Hirth & Kohlstedt, 2003): In Equation 4, values for A and α (constants), p (grain size exponent), r (water fugacity exponent), E (activation energy), and V (activation volume) are taken from Hirth and Kohlstedt (2003), and listed in Table S2 in Supporting Information S1.R is the gas constant and following van der Wal et al. ( 2015) we set the melt fraction (φ) to zero and use a pressure (P) gradient of 0.033 GPa/ km with depth.The three remaining parameters are temperature (T), grain size (d) and water content (fH 2 O).For temperature, variations with respect to a background temperature (the geotherms of Turcotte and Schubert (2002)) are obtained from the temperature derivative of seismic wave velocities detailed in Karato (2008), where the seismic wave velocities are taken from the seismic velocity perturbations of the ANT-20 seismic model (Lloyd et al., 2020).Full details of this method can be found in the study by van der Wal et al. (2013).For grain size and water content, we try six combinations of values, using grain sizes 1, 4, 10 mm, and completely wet/dry, following van der Wal et al. (2015).
As the flow laws derived by Hirth and Kohlstedt (2003) are only valid to a depth of approximately 400 km, below this we use a 1D Maxwell linear rheology with viscosity of 1 × 10 20 Pa s and a lower mantle viscosity of 5 × 10 21 Pa s.

2013 Earthquake
Studying the postseismic deformation following the 2013 earthquake presents an opportunity to place bounds on the Earth structure and properties in the northern Antarctic Peninsula by comparing model predictions with GPS observations.Results are presented in this section for the four different Earth models tested-1D Maxwell (Section 4.1.1),1D Burgers (Section 4.1.2),3D slab model (Section 4.1.3),and 3D composite power-law rheology (Section 4.1.3).

1D Maxwell
Figure 4 shows the viscosity profiles for the suite of 267 1D Maxwell models with the full range of parameters tested shown in gray, and the "best fit" models, that is, those with WRMS errors (see Equation 2) in the 95% confidence interval shown in red.We show the fit to the subset of GPS time series in the north and east directions only.Figure 4a shows the range of best fit models where no ice-related correction has been made to the GPS observations and Figures 4b-4e show the equivalent where ice-mass loss related deformation has been accounted for using the 4 different models described in Section 2.4.
The Earth structure is more tightly constrained when ice related deformation has been corrected for, that is, the range of best fitting parameters is narrower, although the minimum WRMS value is slightly lower when this correction is not made.The ice-related correction that gives the best fit is show in Figure 4e, and the resulting best fitting Earth models have a thin lithosphere (20 km) with an underlying asthenosphere of viscosity 3-4 × 10 18 Pa s.Below 210 km depth a high viscosity of 1 × 10 20 Pa s is required for the rest of the upper mantle.

1D Burgers
Figure 5 shows the viscosity profiles for the suite of 164 1D Burgers models tested in gray with best-fitting models in blue/green.As for the Maxwell model, the Earth structure is more tightly constrained when using GPS time series that have been corrected for the effects of ice-mass change deformation.The minimum WRMS is slightly lower for the best fit Burgers model than the Maxwell model indicating a slightly better fit to the GPS is achieved.The best fitting model to GPS time series that have been corrected for the effects of ice-mass change (Figure 5e) has an asthenosphere layer between 20 km depth and 140 km depth with transient viscosity 4 × 10 17 Pa s and steady state viscosity 4 × 10 18 Pa s.Below this, the upper mantle has a viscosity of 5 × 10 20 Pa s.
As mentioned in Section 3.1, the finite element model has been run at 1-year time steps for the full suite of Earth structures tested.This is due to the large computational time of each model run.However, by using 1-year time steps we may be missing much of the response from the transient viscosity.We set the transient viscosity to be one order of magnitude lower than the steady-state viscosity and the range of values tested is 1 × 10 17 Pa s-1 × 10 18 Pa s.Due to this low viscosity we expect that transient deformation is occurring over the weeks and months immediately following the earthquake and therefore it may be missed by calculating output at 1-year intervals.
To investigate this, we re-compute postseismic deformation with the best fitting Burgers model at a higher temporal resolution (with 20 steps per year, equivalent to approximately 2.5 weeks) for the first 2 years.Figure 6 shows the GPS time series at OHI3 in the east and north directions, along with the model-predicted deformation for the best fitting model at the two temporal resolutions.When calculated at shorter time steps, the model computes more of the transient deformation.At approximately 1 year after the earthquake the two models converge, demonstrating that the transient signal has diminished within this time.This demonstrates that when using Burgers rheology, the finite element model needs to be run at a suitable temporal resolution to capture the initial transient deformation.
The WRMS for the best-fit model increases marginally from 1.68 to 1.69 mm when we use higher temporal resolution possibly because no clear transient signal is apparent in the GPS time series.

3D Earth Models
We tested two different 3D models to see whether introducing lateral heterogeneity in Earth structure could provide a better fit to the GPS data.The best fitting simple 3D subducted slab structure has an asthenosphere viscosity of 4 × 10 18 Pa s and low viscosity wedge above the slab of 9 × 10 17 Pa s giving a minimum WRMS of 1.81 mm.The best model with composite power-law rheology uses a grain size of 1 mm and zero water content to give a minimum WRMS of 2.31 mm.These WRMS results are both higher than the best fitting Burgers model, suggesting that whilst 3D Earth structure may be important, many iterations of the model would be needed to fully investigate all potential Earth structures or input parameters to find a better fit, if one exists.
Figure 7 shows the model-predicted postseismic displacement time series at GPS locations used in the WRMS calculation for each best fit Earth model (Maxwell, Burgers-high temporal resolution, 3D slab, 3D composite powerlaw rheology) along with the GPS time series uncorrected (in cyan) and corrected (in gray) for changes in ice-related deformation (using correction  number 4 in Section 2.4).Prior to the earthquake (dashed line in Figure 7) these GPS time series overlap, with the correction only affecting positions after the earthquake.There is little difference in model-predicted displacement over this time scale between the 1D models and the 3D slab model.The 3D composite power-law model shows larger differences at some sites, particularly in the east direction.

2003 Earthquake and Combined Displacement
Since very few of the GPS sites (OHI2, OHI3, FREI, PALM) were operational at the time of the 2003 earthquake, it would not be a rigorous test to attempt to constrain Earth properties from modeling this earthquake.Instead, we use the best fitting Earth model with Burgers rheology from the analysis of the 2013 earthquake and run a forward model to predict the postseismic deformation across the Antarctic Peninsula from the 2003 event.
Figure 8a shows horizontal postseismic deformation rates for GPS locations on the Peninsula for the period 2005 to 2010.At all stations the deformation rate is less than 0.5 mm/yr.Figure 8b shows model-predicted postseismic deformation rates for the period 2015 to 2020 for the 2003 earthquake (in blue) and the 2013 earthquake (in green).There is more than twice the postseismic deformation from the larger magnitude 2013 earthquake (maximum 1.3 mm/yr) than the 2003 earthquake (maximum 0.3 mm/yr).Figure 8c shows the combined displacement for both earthquakes for 2015 to 2020, peaking at 1.65 mm/yr at O'Higgins, and Figure 8d shows the combined displacement 50 years after the 2013 earthquake.Whilst the magnitude of postseismic deformation diminishes with time, there is still 1.2 mm/yr deformation decades later.The model predicts that the direction of motion will change with time toward the east.

Earth Models and Rheology
In this study we have used GPS observations before and after the 2013 Scotia Sea earthquake to constrain Earth properties in the Antarctic Peninsula.We have selected a range of candidate Earth structures and rheologies and tested a wide parameter space.However, due to the large number of parameters it is clearly not possible to test all combinations with a computationally expensive model.The best-fitting models presented here therefore represent the statistical best fit to the data from the parameter space tested.Other combinations of parameters may produce a good fit that lies just outside the confidence interval, with trade-offs, for example, between layer thickness and viscosities.
The best fitting model (WRMS of 1.68 mm, Figure 5e) has a Burgers rheology with an asthenosphere transient viscosity of 4 × 10 17 Pa s and steady-state one order of magnitude higher, a ratio that is fixed in our model.A Maxwell model also provides a good fit to the data with WRMS only slightly higher (1.79 mm, Figure 4e) and the same steady state viscosity (4 × 10 18 Pa s), indicating that the "steady state" is the dominating signal.When we run a higher temporal resolution model to enhance the transient signal, we do not find a better fit to the data, and in fact the transient signal is undetectable by eye in any of the GPS time series.A different ratio of transient to steady-state viscosity may be required to improve the fit, however, we expect only a marginal improvement in results given the limited number of observations to test against.Transient deformation following other large earthquakes has been observed in far-field GPS time series (e.g., in North America (Freed et al., 2007)), however, in the Northern Peninsula this rapid deformation may be obscured by noise or other processes such as seasonal signals.
Initial tests with 3D Earth models do not improve the fit to the GPS.The best fit simple 3D slab model has an asthenosphere viscosity of 4 × 10 18 Pa s and low viscosity wedge above the slab of 9 × 10 17 Pa s.Because the asthenosphere viscosity is the same as it is for the best fitting 1D Earth models, this demonstrates that a 1D model is a good approximation for this small area when testing against observations from 6 GPS sites.Only an increase in the number or geographical distribution of GPS stations could enable a 3D Earth model to be constrained.The composite power-law rheology combined with high resolution seismic model was the worst fitting model tested, likely due to the large number of unknown input parameters to this flow-law (e.g., water content, grain size) and our results do not suggest that such a model setup could not produce equal or better fit to the data.This method of 3D modeling for postseismic deformation still has potential if enough variations of input parameters can be tested, or new data become available to inform of likely values for some of the inputs.This should be considered in future studies.mantle viscosity of 1 × 10 18 Pa s, using data that pre-dates the 2013 earthquake, although their preferred lithosphere was much thicker (120 km) than that found in our study (20 km).Samrat et al. (2020) also found a thick lithosphere (105 km) but a lower value of upper mantle viscosity was preferred (9 × 10 17 Pa s), although their GPS data included the post-2013 period and they did not correct for postseismic deformation.
Our study has investigated a wide parameter space by including an asthenosphere of varying thickness and viscosity as well as different rheological models.The asthenosphere thickness for our best fitting Maxwell and Burgers model is 190 and 140 km, respectively, thinner than the equivalent layer in the GIA studies.This layer also has slightly higher steady-state viscosity than the GIA studies.However, because the GIA studies mentioned above did not consider Burgers rheology, it may mean that the upper mantle viscosities they constrained are some average of the transient and steady-state values suggested by our study.Nevertheless, given the model uncertainties and the sparsity of GPS observations in the Antarctic Peninsula the Earth structure constrained by the GIA studies and our postseismic study are quite similar.

Sensitivity to Choice of Data
The choice of which data to use when calculating the WRMS clearly impacts the resulting best fitting models.
Here we discuss the impact of these choices (detailed in Section 2.5).
We use GPS observations in the north and east direction only since the strike-slip earthquakes produce minimal vertical displacement and the vertical is more dependent on correcting for ice-related deformation.Overall, the results do not change when using all three directions of motion, but we get a poorer fit to the data with the minimum WRMS increases for all Earth models tested.The best fitting Earth model remains the same, but the WRMS increases from 1.68 to 2.21 mm.
The second choice was which ice-related correction to make to the GPS data prior to the WRMS calculation.We tested four different corrections along with not making this correction at all.As described in Section 4.1, applying an ice-related deformation correction to the GPS time series results in a more tightly constrained Earth model than when we don't use a correction, although the minimum WRMS error is slightly higher.The four different ice-related corrections are based on the work by Samrat et al. (2020) and can be divided into two sub-groups based on upper mantle viscosity of 9 × 10 17 Pa s and 5 × 10 18 Pa s.The lower WRMS, and therefore better fitting models to the GPS data, are obtained with an ice-related correction based on mantle viscosity of 5 × 10 18 Pa s (Figures 4d and 4e and Figures 5d and 5e).A much poorer fit is achieved when using the correction based on the published best fit model from Samrat et al. (2020) with upper mantle viscosity of 9 × 10 17 Pa s (Figures 4b  and 5b).Samrat et al. (2020) constrained this best fit model by comparing ice-related deformation to GPS time series that were uncorrected for postseismic deformation, using east and vertical deformation between 2001 and 2018 (i.e., includes post-2013 earthquake related deformation).This means that there will be bias in these results which, when used as a correction in our study, carries this bias through and results in a poorer fit.By adjusting the Earth model used in the ice correction, we get a better fit.
Using this method to model either postseismic deformation or ice loading deformation and constraining an Earth model using GPS observations creates a circular problem.Because the GPS time series first need correction for the other process, which is currently a model-based estimate also been constrained using GPS time series, there is always a bias in the results.To further advance our understanding of these deformations, either a joint inversion is required or a model that is capable of simulating both processes occurring concurrently.

Postseismic Deformation
Our modeling predicts a horizontal postseismic deformation of up to 1.65 mm/yr at O'Higgins, the northern-most GPS site on the Peninsula, for the period 2015-2020.The vertical postseismic deformation is much smaller, predicted to be −0.35mm/yr.The horizontal deformation persists decades after the events, reducing to 1.2 mm/ yr 50 years later (Figure 8d).As the deformation diminishes with time, the direction of motion changes toward the east.Since sites in the Antarctic Peninsula are contaminated by horizontal postseismic deformation they should not be used for studies of plate rotation, confirming the site selection used in the ITRF2014 plate motion model (Altamimi et al., 2017).
In their study of postglacial rebound, Argus et al. (2011) identified residual horizontal motion at O'Higgins of 2.3 mm/yr south east relative to the Antarctic Plate.The authors were unable to explain the cause of this residual motion and highlight that it is the opposite direction as that expected from ice-mass loss.The data used in that study pre-dates the 2013 earthquake, but we suggest that postseismic deformation from the 2003 Scotia Sea earthquake could be contributing 0.3 mm/yr east motion to this residual.This suggestion is further supported by the study of Dietrich et al. (2004) who report a lower horizontal residual deformation rate relative to rigid plate motion at O'Higgins of 1-2 mm/yr, using GPS data collected prior to the 2003 earthquake.
We have shown that large (magnitude 7.6 and 7.7) earthquakes occurring at a distance of 650 km from the tip of the Peninsula result in widespread deformation of over 1 mm/yr which persists for decades.There may be more large earthquakes, or smaller but closer earthquakes, that have occurred in the last 50 years that are contributing to the present-day deformation rates on the Peninsula, and this is worth further investigation.

Conclusions
Our modeling estimates the northern Antarctic Peninsula is deforming with a displacement rate up to 1.65 mm/ yr in the horizontal direction during 2015-2020 following the two large strike-slip earthquakes that occurred in the Scotia Sea in 2003 and 2013.By comparison, model-predicted vertical postseismic deformation is small, with a displacement rate peaking at −0.35 mm/yr for the same time period.Comparing model-predicted deformation from a suite of possible Earth structures to GPS time series suggests that a Burgers rheology with an asthenosphere transient viscosity of 4 × 10 17 Pa s and steady-state one order of magnitude higher best explains the observations.
A previous study (King & Santamaria-Gomez, 2016) also found postseismic deformation in East Antarctica following a large earthquake that occurred ∼600 km from the coastline.Significant postseismic mantle flow and associated movement of mass may cause continent-wide changes in the gravity field which may have an effect on GRACE gravity data.Furthermore, geodetic measurements such as GPS are routinely used to constrain models of GIA which in turn are used to correct GRACE gravity data to estimate present-day ice-mass change and hence sea-level change in Antarctica.The postseismic signal has not previously been accounted for meaning that there are potential biases in results.More work needs to be done to further investigate the extent and magnitude of postseismic deformation across the continent.
A question that remains however, is how to accurately isolate deformation from specific processes in the GPS records before using them to constrain Earth properties, particularly when processes are occurring on the same timescales.The best-fitting Earth models after such analysis may produce different viscosity profiles and inconsistent rheologies if signals have not been correctly accounted for.This problem is non-trivial.It may be possible to separate the various spatial-temporal processes using principal or independent component analysis, or other data separation approaches such as those based on machine learning (Dong et al., 2006;Gao et al., 2022;Liu et al., 2018;Mandler et al., 2021).A viscoelastic model that includes both earthquake-related slip and a changing ice load would also be a step forward in modeling Antarctic deformation in a consistent way.

Figure 1 .
Figure1.The Antarctic Plate and associated tectonic plate boundaries(Bird, 2003).Focal mechanisms of earthquakes mentioned in the text shown in black and all post-1976 earthquakes with M w ≥ 7.0 shown in red (from the Global CMT Catalog(Ekström et al., 2012)).

Figure 2 .
Figure 2. (a) Coseismic offsets estimated from GPS (blue) and model (red) for the 2013 earthquake.Error ellipses are 2-sigma and shown at the same scale as the vectors.Panel (b) slip and rake of the 2003 earthquake, each point represents a node on the slip plane.Model inputs given in top left, longitude and latitude refer to the fault plane center point on the surface.Panel (c) as for b but for the 2013 earthquake.Focal mechanisms are from the Global CMT Catalog(Ekström et al., 2012).
1. Published best fit model: lithospheric thickness 105 km and upper mantle viscosity 9 × 10 17 Pa s. 2. Best fit model using GPS time series only up until the time of the 2013 earthquake (to avoid any potential contamination from postseismic deformation): lithospheric thickness of 110 km and upper mantle viscosity 9 × 10 17 Pa s. 3. Best fit model using the PALM east time series only up until the time of the 2013 earthquake: lithospheric thickness of 55 km and upper mantle viscosity 5 × 10 18 Pa s. 4. Best fit model using the DUPT east time series only up until the time of the 2013 earthquake: lithospheric thickness of 75 km and upper mantle viscosity 5 × 10 18 Pa s.
Density, Young's Modulus and Poisson's Ratio derived from PREM.Shaded regions depict the thickness range of each layer tested where EL = elastic lithosphere, A = asthenosphere, UM = upper mantle, LM = lower mantle, labels indicate the range of viscosities tested.
) (their Figure6), rounding the geometry to accommodate 10 km elements.We then split the fault plane into three sections each with a different slip magnitude averaged from the final model ofYe et al. (2014) (their Figure7).Again, we verify the resulting seismic moment (3.5 × 10 20 Nm) is close to the value reported byYe et al. (2014) (4.0 × 10 20 Nm).The model values are shown in Figure 2b.

Figure 3 .
Figure 3. (a) Plan view of 3D model containing subducted slab structure with cross section X-X' shown in (b).Earthquake locations shown by black stars in (a).

Figure 4 .
Figure 4. Viscosity profiles for 1D Maxwell rheological models for the upper most 400 km of the model.Gray lines show the full range of values tested and red lines indicate models with WRMS within the 95% confidence interval.Minimum WRMS given in each panel (in mm).WRMS calculated using: (a) GPS data not corrected for the effects of ice loading; (b)-(e) GPS data corrected for the effects of ice loading using a different model (1-4 in Section 2.4).

Figure 6 .
Figure 6.GPS (a) east and (b) north time series (with ice load correction) at OHI3 (gray dots) with model-predicted postseismic displacement using Burgers rheology at 1-year time steps (blue line) and ∼2.5 week time steps (maroon dashed line).Minimum WRMS (in mm) given for each model in panel (b).

Figure 5 .
Figure 5.As for Figure 4 but for 1D Burgers rheological models where the green (steady-state viscosity) and blue (transient viscosity) lines indicate the models with WRMS within the 95% confidence interval.Minimum WRMS given in each panel (in mm).

Figure 7 .
Figure 7. Time series of GPS observations with pre-earthquake trend removed to reveal postseismic deformation (cyan) and corrected for the effects of ice-loss related deformation (gray) in east (a), north (b), and vertical (c) directions.Best fit modeled postseismic deformation for Maxwell model (red lines), Burgers model (blue dashed lines), 3D slab Maxwell model (orange lines) and 3D power-law model (green dashed lines).** indicates that the GPS time series were used in the WRMS calculation.Gray dashed vertical line shows the time of the 2013 earthquake.

TheFigure 8 .
Figure 8.(a) modeled horizontal displacement rates for 2005-2010 for the 2003 earthquake.Panel (b) modeled horizontal postseismic displacement rates for the period 2015-2020 for the 2003 earthquake (blue), the 2013 earthquake (green); (c) combined modeled horizontal displacement rates for the 2003 and 2013 eq for 2015 to 2020; (d) as for (c) for 2060 to 2065.Earthquake locations shown by black stars.

Table 1
Details of the GPS Sites Located Within 1,200 km of the 2013 Earthquake

Table 2
Earth Model Layers, Properties and Mesh Resolution for the Low-and High-Resolution Parts Respectively (High Resolution Part Stops at the Lower Mantle Boundary)