GPS constraints on deformation in northern Central America from 1999 to 2017, Part 1 – Time-dependent modelling of large regional earthquakes and their post-seismic effects

SUMMARY

We use continuous and campaign measurements from 215 GPS sites in northern Central America and southern Mexico to estimate coseismic and afterslip solutions for the 2009 M_w = 7.3 Swan Islands fault strike-slip earthquake and the 2012 M_w = 7.3 El Salvador and M_w = 7.4 Guatemala thrust-faulting earthquakes on the Middle America trench. Our simultaneous, time-dependent inversion of more than 350 000 daily GPS site positions gives the first jointly consistent estimates of the coseismic slips for all three earthquakes, their combined time-dependent post-seismic effects and secular station velocities corrected for both the coseismic and post-seismic deformation. Our geodetic slip solutions for all three earthquakes agree with previous estimates that were derived via static coseismic-offset modelling. Our time-dependent model, which attributes all transient post-seismic deformation to earthquake afterslip, fits nearly all of the continuous GPS site position time-series within their several-millimetre position noise. Afterslip moments for the three earthquakes range from 35 to 140 per cent of the geodetic coseismic moments, with the largest afterslip estimated for the 2012 El Salvador earthquake along the weakly coupled El Salvador trench segment. Forward modelling of viscoelastic deformation triggered by all three earthquakes for a range of assumed mantle and lower crustal viscosities suggests that it accounts for under 20 per cent of the observed post-seismic deformation and possibly under 10 per cent. Our results thus point to afterslip as the primary and perhaps dominant mode of post-seismic deformation for these three earthquakes. Forward modelling of post-seismic deformation associated with the larger M_w = 7.6 September 2012 Costa Rica thrust earthquake suggests that afterslip, viscoelastic flow, or some combination thereof was responsible for a significant change in motion observed at a GPS site on San Andres Island in the Caribbean Sea more than 500 km from all four earthquakes. The measurable effects of the 2009 and 2012 earthquakes on the motions of GPS sites in nearly all of northern Central America underline the importance of time-dependent calibrations for transient, earthquake-related effects for studies of steady-state deformation processes.

Seismic cycle, Transient deformation, Space geodetic surveys, Satellite geodesy

1 INTRODUCTION

GPS measurements in northern Central America and southern Mexico, including Guatemala, Honduras, El Salvador, and the Mexican state of Chiapas (Fig. 1), began between 1999 and 2003, with broad goals of better understanding plate motions and seismic hazards in this complexly deforming region. Positioned at the western end of the Caribbean plate, northern Central America is located near the intersection of three major fault zones, consisting of the Motagua-Polochic fault zone, which accommodates ≈20 mm yr⁻¹ of left-lateral strike-slip movement between the Caribbean and North America plates (Fig. 1), the Middle America trench, where the Cocos plate subducts northeastwards at 70–80 mm yr⁻¹ (DeMets et al. 2010), and the Central America volcanic arc faults, which accommodate 10–15 mm yr⁻¹ of northwestward translation of the Central America forearc sliver towards a poorly understood continental triple junction in southern Guatemala (e.g. DeMets 2001; Authemayou et al. 2011).

Figure 1.

Seismotectonics and geography of the study area. Dashed lines denote subducting slab contours. Red stars show epicentres for major Central American earthquakes since 2009 and the 1976 Guatemala earthquake. Associated earthquake focal mechanisms are the global centroid-moment tensors (Ekström et al. 2012). Black and blue vectors show Cocos plate velocities relative to the North America and Caribbean plates, respectively (DeMets et al. 2010).

Efforts to better understand these faults are strongly motivated by their long histories of destructive earthquakes, including the M_w = 7.5 1816 Polochic fault earthquake (White 1985), the M_w = 7.5 1976 Motagua Fault earthquake (Plafker 1976), and numerous destructive M6–6.5 20th-century tectonic earthquakes along the volcanic arc (White & Harlow 1993). In particular, an improved understanding of the regional seismic hazards requires the following: (1) Well-determined estimates of the angular velocities that best describe the motions of the regional plates and blocks. (2) The degree of interseismic locking along faults. (3) The depths at which interseismic locking and afterslip occur. (4) Strain-rate tensors to quantify rates and directions of distributed deformation across the region’s numerous lesser, but still-hazardous faults (e.g. Caceres et al. 2005). (5) Modelling of Coulomb stress changes to quantify possible earthquake triggering relationships (e.g. Martinez-Diaz et al. 2004; Graham et al. 2012).

The results described below are the outcome of an long-term international effort to achieve these goals. In this first part of our two-part analysis, we use time-dependent modelling of GPS measurements between 1999 and 2017 at more than 200 stations in northern Central America to estimate geodetic coseismic slip and post-seismic afterslip solutions for the 2009 May 28 M_w = 7.3 Swan Islands fault earthquake, the 2012 August 27 M_w = 7.3 El Salvador earthquake, and the 07 November 2012 M_w = 7.4 Champerico (Guatemala) earthquake (Fig. 1). Using a static-deformation approximation and subsets of the data used herein, previous authors have estimated geodetic coseismic slip solutions for each of these earthquakes (Graham et al. 2012; Ellis et al. 2015; Geirsson et al. 2015) and afterslip solutions for two of the three earthquakes (Ellis et al. 2015; Geirsson et al. 2015). Here, we simultaneously invert all available GPS data from the region to jointly estimate coseismic and afterslip solutions for all three earthquakes, thereby assuring that their summed elastic effects correctly describe the evolving space-time pattern of GPS station positions in northern Central America between 1999 and the present. Viscoelastic deformation may also have been triggered by these earthquakes or the M_w = 7.6 05 September 2012 Costa Rica earthquake several hundred kilometres south of our study area; for completeness, we evaluate whether this is consistent with transient post-seismic deformation recorded in the region since 2009.

GPS site velocities that are corrected in a consistent manner for the coseismic and post-seismic effects of all three earthquakes are a key outcome of our analysis. In Part 2 of our analysis, we invert these newly determined GPS velocities to estimate angular velocities for plates and blocks in our study area, interseismic locking along several major faults, and strain-rate tensors for areas where distributed deformation occurs (Ellis et al. ‘GPS constraints on deformation in northern Central America from 1999 to 2017 , Part 2: Block rotations, fault slip rates and locking, and distributed deformation’, manuscript in prep. 2018).

2 GPS DATA

2.1 Data and processing methods

To minimize possible biases in our results from different strategies for processing or post-processing GPS data, we strove to analyse all the data using identical methods. We compiled daily GPS RINEX files from 74 continuous and 141 campaign GPS sites, including all sites in Belize, El Salvador, Guatemala, Honduras, the Mexican state of Chiapas, the western portion of the Caribbean Sea, and some stations in Nicaragua (Fig. 2; Supporting Information Table S1). The campaign observations include some or all data from DeMets (2000, 2004, 2007a,b, 2008a,b, 2009, 2011a,b,c), DeMets & Tikoff (2015a,b,c,d,e,f,g,h), Dixon (2001, 2003, 2004, 2010), Franco et al. (2012), LaFemina (2013a,b,c,d,e), Lyon-Caen et al.(2006), Newman (2010), Schwartz & Dixon (2000), and Stalleret al. (2016). The earliest GPS data from our study area are from 1999 (Supporting Information Table S1), although observations at several sites on the Caribbean plate or nearby areas of the North America plate extend back to 1993 (Supporting Information Table S1). Seventy per cent of the GPS sites in our study area became operational before the May 2009 Swan Islands earthquake, the earliest earthquake modelled in this study (Fig.2).

Figure 2.

Locations of GPS stations used in this study. Light green patches show the rupture areas of three major earthquakes since 2009 (see the text). Focal mechanisms are from Fig. 1 and coincide with the earthquake epicentres. EQ; earthquake.

We processed all of the GPS code-phase data with release 6.3 of the GIPSY software suite from the Jet Propulsion Laboratory. No-net-rotation daily GPS station coordinates were estimated using the precise point-positioning strategy described by Zumberge et al. (1997). Our processing methodology includes constraints on a priori tropospheric hydrostatic and wet delays from Vienna Mapping Function parameters (http://ggosatm.hg.tuwien.ac.at), elevation-dependent and azimuthally dependent GPS and satellite antenna phase centre corrections from IGS08 ANTEX files (available via ftp from sideshow.jpl.nasa.gov), and FES2004 corrections for ocean tidal loading (http://holt.oso.chalmers.se). Phase ambiguities were resolved using GIPSY’s single-station ambiguity resolution feature (Bertiger et al. 2010). Daily no-net-rotation station location estimates were transformed to IGS08, which conforms to ITRF2008 (Altamimi et al. 2011), using daily seven-parameter Helmert transformations from the Jet Propulsion Lab. We estimated daily correlated noise between stations from the coordinate time-series of linearly moving continuous stations outside the study area (Marquez-Azua & DeMets 2003). Corrections of the raw daily GPS site positions for this common-mode noise reduced the daily scatter and amplitude of the longer-period noise in the GPS time-series by 20–50 per cent. All GPS coordinate time-series were also corrected for equipment-related offsets and other discontinuities not related to earthquakes. Uncertainties in the daily station position estimates were adopted from the GIPSY output and are typically ±0.6 mm in longitude, ±0.5 mm in latitude, and ±2.5 mm in elevation.

One of our continuous GPS sites, VMIG in El Salvador, is impacted by multipath noise due to a partially blocked antenna(Supporting Information Fig. S1). To mitigate this noise, we subdivided the station’s non-linear position time-series into five approximately linear segments. We then applied a filter based on a linear Theil Sen regression (Sen 1968; Blewitt et al. 2009) to identify and remove daily position outliers for each segment of ±5 mm in the horizontal components and ±10 mm in the vertical.

3 METHODS: TIME-DEPENDENT MODELLING

Fault afterslip and viscoelastic relaxation are the primary causes of transient post-seismic deformation following large earthquakes, although their relative contributions are hard to determine given limited information about crust and mantle rheologies and the location, magnitude, and temporal characteristics of afterslip (Hu et al. 2004; Suito & Freymueller 2009; Hu & Wang 2012; Wang et al. 2012; Kogan et al. 2013; Sun & Wang 2015). In our analysis, we model the two processes separately to estimate a maximum bound for afterslip and to find approximate upper and lower bounds for viscoelastic deformation. None of the observed GPS time-series are corrected for the modelled viscoelastic deformation prior to our TDEFNODE inversion. All transient post-seismic deformation is thus attributable to fault afterslip estimated in our model inversions.

3.1 Inverse modelling with TDEFNODE

TDEFNODE, which calculates time-dependent or static deformation in a homogeneous elastic half-space based on equations from Okada (1985), can concurrently estimate coseismic slip solutions, afterslip solutions, afterslip decay rates, interseismic fault locking, linear station velocities, plate/block angular velocities, and strain-rate tensors from campaign and continuous GPS position time-series and other types of geodetic, seismologic, and plate kinematic data (McCaffrey 2009). For the inversion described in this paper, coseismic and afterslip solutions, afterslip decay rates, linear station velocities, and seasonal periodic signals are estimated. At a given GPS station, the time-dependent movement d(t) due to afterslip triggered by an earthquake at time t_eq is modelled in TDEFNODE as d(t)∼ = ∼A*log₁₀(1∼ + ∼(t − t_eq)/t_c), where A is an amplitude term related to a site’s elastic response to earthquake afterslip and t_c is the logarithmic decay constant.

Faults in the TDEFNODE elastic half-space are defined by sets of nodes that approximate the fault trace and dip. Green’s functions that quantify the 3-D surface elastic response to unit slip at each fault node are calculated and gathered to form the basis for the time-dependent inversion. Simulated annealing and grid search iterations are employed to minimize the reduced chi-square statistic. For our inversion, slip on the Swan Islands Fault is estimated as uniform slip on six rectangular patches that span the 2009 earthquake rupture zone (Section 4.1). In contrast, slip values for the Middle America Trench are estimated at each fault node, with smoothing applied to avoid implausible node-to-node variations in slip values (Sections 4.2 and 4.3).

3.2 Fault meshes and slip resolution: Swan Islands fault and Middle America trench

For the Swan Islands fault, we adopt the fault geometry of Graham et al. (2012), who subdivided the fault into six ≈50-km-long rectangular patches (Fig. 3) using the fault trace defined by a Sea Marc II seafloor survey (Rosencrantz & Mann 1991; Rogers & Mann 2007) and assuming that the fault is vertical to a depth of 15 km. Beginning with a 2009 GPS network geometry similar to ours, Graham et al. (2012) use a series of checkerboard tests to determine how well slip could be resolved along the six rectangular patches they used to represent the Swan Islands fault. An inversion of noisy synthetic data recovered the starting slip values for five of the six slip patches that comprise the fault to within 10-20 per cent, failing only for the easternmost fault patch farthest from the GPS network.

Figure 3.

(a) TDEFNODE geodetic slip solution for the 2009 Swan Islands earthquake. (b) Graham et al. (2012) geodetic slip solution. (c) TDEFNODE earthquake afterslip solution through 2017.0. (d) Site motions estimated with our best coseismic and afterslip solutions are shown in panels (a) and (c). The white and black stars are earthquake centroids from the Global CMT catalogue (Dziewonski et al. 1981; Ekström et al. 2012) and United States Geological Survey, respectively.

We approximated the Middle America subduction interface using the Slab 1.0 geometry of Hayes et al. (2012). Fault node spacings are ≈20 km along-strike and 10-km down-dip. From a standard checkerboard test, Ellis et al. (2015) find that inversions of synthetic data created from slip patches at various depths off the coast of Guatemala near the Champerico earthquake recover the starting models relatively well at depths of 10 km to 60 km (details are given in the Ellis et al. (2015) supplementary materials). Checkerboard tests using the sparser GPS network configuration from El Salvador, Nicaragua, and Honduras that existed during the 2012 El Salvador earthquake (Supporting Information Fig. S2) indicate that a shallow slip patch offshore from Nicaragua and the Gulf of Fonseca can be recovered (Supporting Information Fig. S2B), but that assumed deep slip (Supporting Information Fig. S2C) is more difficult to recover (Supporting Information Fig. S2D). Care is thus required in interpreting the coseismic and post-seismic results for the El Salvador earthquake described below.

Given the limited ability of the land-based GPS stations to resolve offshore slip during all three earthquakes modelled herein, we applied smoothing when estimating their coseismic and afterslip solutions. For modelling the 2009 Swan Islands earthquake, Graham et al. (2012) found that subdividing the Swan Islands fault into along-strike patches smaller than ≈50 km did not significantly improve the model fit, whereas using larger sub-fault dimensions degraded the fit. We thus adopt their optimized fault characterization and do not smooth the independent slip values estimated in our inversion for the six rectangular patches. For modelling of the 2012 El Salvador and Champerico subduction zone thrust earthquakes, we used a spread penalty smoothing parameter to penalize large slip values at distances progressively farther from the slip centroid (McCaffrey 2009). For the latter two earthquakes, we also tested but rejected simpler 2-D Gaussian coseismic or afterslip slip sources because they significantly degraded the fits to the GPS time-series relative to our preferred solutions.

3.3 Viscoelastic deformation: forward modelling with VISCO-1D

In order to predict post-seismic viscoelastic deformation, we use VISCO-1D software (version 3), which solves for the viscoelastic deformation of a spherically layered Earth given a user-defined viscosity structure and a linear Maxwell or bi-viscous rheology for each layer (Pollitz 1997). To explore the likely range of short-term and long-term viscoelastic responses, we used three different viscosity structures appropriate for continental crust (Supporting Information Fig. S3). We estimate an upper limit for the long-term viscoelastic deformation using the Hearn et al. (2013) M1 Earth structure (left panel in Supporting Information Fig. S3), which Hearn et al. used to find an upper bound for viscoelastic deformation triggered by historic earthquakes in California. We estimate an upper bound for shorter-term viscoelastic deformation (i.e. during the first few years after an earthquake) by decreasing the mantle viscosity assumed in the M1 Earth model from 3 × 10¹⁸ to only 5 × 10¹⁷ Pa s, as estimated by Hu & Wang (2012) from their analysis of short-term post-seismic deformation for the 2004 Sumatra earthquake. We refer to this Earth structure, which is illustrated in the right panel in Supporting Information Fig. S3, as the maximum-response model. Finally, we approximated a lower limit for viscoelastic deformation by increasing the M1 mantle viscosity from 3 × 10¹⁸ to 1 × 10¹⁹ Pa s. We refer to this Earth structure as our minimum-response model (centre panel in Supporting Information Fig. S3).

4 RESULTS: COSEISMIC AND POST-SEISMIC DEFORMATION

For the 215 GPS sites in our study area, 357 627 total observations constrain our best-fitting model, consisting of the northing, easting, and vertical daily position estimates at the sites. The best-fitting model, consisting of the locations and magnitudes of coseismic slip for all three earthquakes, afterslip decay constants, locations, and amplitudes for each earthquake, and a linear velocity for each GPS site, is described by 170 adjustable parameters.

Reduced chi-square for the best-fitting solution is 14.0 of which the misfit between the model and data contributes 70 per cent and the smoothing penalties contribute the remainder. The portion attributable solely to the data misfit implies that the misfits are, on average, 2.8 times larger than the assigned data uncertainties, which are typically ±0.5–0.6 millimetres in the daily horizontal components of a station position and ±2.5 mm in the vertical component. At the 74 continuous GPS stations, the weighted root-mean-square (WRMS) misfits of our best-fitting model are ≈2 mm for the horizontal position components and 7 mm in the vertical. At the 141 campaign GPS sites, the WRMS misfits are ≈3 mm for the two horizontal components and 8 mm in the vertical. These misfits are close to the 1–3 mm scatter that is typical of the daily positions of well-behaved continuous GPS sites. We conclude that our best-fitting model is consistent with the GPS observations. Results for each of the three earthquakes, including their predicted viscoelastic deformation, are presented next. Supporting Information Table S2 lists the 3-D coseismic offsets at all 215 sites for all three earthquakes modelled herein. In order to facilitate future analyses with GPS data from northern Central America, the Supporting Information also include an archive with the observed and modelled daily movements for each of the 215 stations that are included in the analysis. The daily position corrections in these files can be used to correct geodetic position time-series for the coseismic and post-seismic effects for all three earthquakes modelled herein.

4.1 2009 Swan Islands earthquake

The 2009 May 28 M_w = 7.3 earthquake on the Swan Islands fault was the largest earthquake on the North America–Caribbean plate boundary since the 1976 M_w = 7.5 Motagua fault earthquake in Guatemala (Plafker 1976; Fig. 1). Inversions by Graham et al. (2012) of GPS station offsets and seismic waveforms for this earthquake give geodetic (Fig. 3b) and seismic slip solutions both with up to 1 m of slip along the western ≈250 km of the Swan Islands fault. Afterslip associated with this earthquake has not been estimated previously.

Fig. 3 compares our new coseismic slip solution (Fig. 3a) to the Graham et al. (2012) geodetic slip solution (Fig. 3b). Our new estimate uses the same fault geometry as did Graham et al. (2012), consisting of six rectangular patches with assumed uniform slip on each patch. Our coseismic slip solution is constrained to have slip amplitudes on each fault patch within 30 per cent of the Graham et al. (2012) solution, which is consistent with seismic waveform constraints. Of the five continuous GPS sites that were operating during the earthquake, the position time-series of the closest site, ROA0 on Roatan Island (Fig. 3d), is well fit by our model (Supporting Information Fig. S4). The position time-series for the other four continuous sites, which are located ≈400 km south of the rupture zone, are also well fit, with residual positions of only a few millimetres relative to the positions estimated with our best-fitting TDEFNODE model (Supporting Information Fig. S5).

For an assumed shear modulus of 40 GPa, the best-fitting geodetic moment from our inversion is 1.2 × 10²⁰ N m (M_w = 7.35), close to the 1.3 × 10²⁰ N m global centroid-moment-tensor (gCMT) catalogue estimate (Dziewonski et al. 1981; Ekström et al. 2012). Maximum slip on the fault was ≈900 mm, close to the ≈1000 mm maximum estimated by Graham et al. (2012).

Most of the information about afterslip triggered by the Swan Islands earthquake comes from the continuous GPS site ROA0 due to its proximity to the earthquake (Fig. 3d). The ROA0 position time-series (Fig. 4 and Supporting Information Fig. S4) shows clear evidence for transient post-seismic motion, particularly in the east position component, which records motion nearly parallel to the N70°E-striking fault. For the best-fitting decay time constant of 30 d, the afterslip estimated with our time-dependent TDEFNODE inversion fits the curvature of the east component of motion within the ±2 mm daily noise (Fig. 4 and Supporting Information Figs S4 and S5). The post-seismic deformation is thus well described by a simple logarithmically decaying afterslip model. In the following section, we consider whether some deformation may also be viscoelastic.

Figure 4.

North (a), east (b) and vertical (c) daily positions (grey circles) for GPS site ROA0 reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2 to emphasize the transient site motion due to earthquake afterslip and viscoelastic rebound. White circles with red outlines are 30-d averages of observed daily positions. Black curve shows TDEFNODE afterslip model. Coloured curves show viscoelastic deformation solely attributable to the 2009 Swan Islands earthquake for all three rheological models described in Section 3.2 and shown in Supporting Information Fig. S2. Dashed lines denote earthquake times. Red area in inset shows the 2009 earthquake rupture zone.

If all of the post-seismic deformation were attributable to fault afterslip, our model implies that up to 480 mm of afterslip occurred by 2017.0 (Fig. 3c). As of 2017.0, the cumulative moment released by afterslip was 3.55 × 10¹⁹ N m, ≈35 per cent of the geodetic earthquake moment. Nearly all of the afterslip was concentrated along the eastern two-thirds of the rupture area (Fig. 3c), with only minimal afterslip at the fault’s western end where it transitions onshore to the seismically hazardous Motaguafault.

A comparison of the post-seismic deformation at sites ROA0 (Fig. 4) and SNJE (Fig. 5d), which are located respective distances of 10-20 km and ≈400 km from the 2009 earthquake rupture zone, shows that the elastic effects of fault afterslip diminish rapidly with distance from the rupture zone, as expected. Despite this rapid drop-off, the afterslip significantly altered the velocities of GPS sites at distances hundreds of kilometres south of the fault for years after the earthquake. Along the north coast of Honduras, ≈50–100 km south of the rupture zone, motion attributable to afterslip was 20 mm by 7.6 yr after the earthquake (red arrows in Fig. 3d), implying a 2.5 mm yr⁻¹ average change in post-seismic versus pre-earthquake site velocities in this region. At site SNJE ≈400 km from the earthquake rupture zone, an observed, 2–3 mm yr⁻¹ northward acceleration of the station motion beginning at the time of the 2009 earthquake persisted until at least 2012 (Fig. 5). At GPS site SAN0 700 km south of the rupture zone, no obvious change in the site velocity occurred after the May 2009 earthquake. The post-seismic effects of the 2009 earthquake were thus small enough to be ignored at locations in Nicaragua and farther south.

Figure 5.

North (a), east (b) and vertical (c) daily positions (grey circles) for GPS site SNJE reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2. The coloured curves show viscoelastic deformations that are individually attributable to the 2009 Swan Islands earthquake and 2012 El Salvador earthquake for all three rheological models described in Section 3.2 and shown in Supporting Information Fig. S2. Red area in the inset shows the 2009 earthquake rupture zone. See Fig. 4 caption for further information.

Fig. 4 shows the viscoelastic deformation predicted through early 2017 at site ROA0 for our Swan Islands earthquake slip solution for all three Earth models described in Section 3.2 and shown in Supporting Information Fig. S3. The maximum-response model predictions (red curve in Fig. 4) overestimate the post-seismic deformation at ROA0 by only ≈10 per cent by 2017.0, but exceed by a factor of three to four the observed post-seismic deformation at site SNJE ≈400 km from the earthquake (Fig. 5). The M1 model of Hearn et al. (2013) correctly predicts the post-seismic deformation measured at SNJE (Fig. 5) but greatly underestimates thedeformation measured at site ROA0 (Fig. 4). None of the three viscoelastic models we tested match all the observations. If the good fit of the M1 model to the deformation measured at SNJE indicates that it better approximates the rheological structure of northern Central America than does the maximum-response model, then its poor fit at ROA0 implies that fault afterslip was responsible for most (≈70 per cent) of the post-seismicdeformation.

4.2 El Salvador earthquake

Slip solutions for the 2012 August 27 El Salvador thrust earthquake have been estimated from seismic waveforms (Ye et al. 2013) and geodetic data (Geirsson et al. 2015). The Ye et al. inversion of 83 broad-band P waves indicates that coseismic slip was concentrated above depths of 30 km and had peak amplitude of 1.2 m. The Geirsson et al. inversions of 22 GPS site offsets also suggest that slip was shallow (above 20 km depth) and generally averaged ≈1 m (Fig. 6b).

Figure 6.

(a) TDEFNODE geodetic slip solution for the 2012 El Salvador earthquake. Subduction depth contours are spaced every 20 km. Black circles indicate the fault nodes where slip is estimated. (b) Geodetic earthquake slip solution from fig. 6a of Geirsson et al. (2015). (c) TDEFNODE earthquake afterslip solution through 2017.0. (d) Site motions predicted by the coseismic and afterslip solutions shown in (a) and (c). White and black stars show the earthquake centroids from the Global CMT catalogue (Dziewonski et al. 1981; Ekström et al. 2012) and the United States Geological Survey, respectively. Named GPS sites are continuous stations.

Coseismic slip (and afterslip) for the 2012 El Salvador earthquake were estimated in our TDEFNODE inversion via smoothing of the slip values at independent nodes that represent the Middle America subduction interface. Our new coseismic slip solution (Fig. 6a) agrees well in location and amplitude with previous estimates (Ye et al. 2013; Geirsson et al. 2015) and has a geodetic moment of 1.3 × 10∼²⁰ N m (M_w = 7.3), the same as the gCMT moment. The coseismic offsets, which point nearly due south (blue arrows in Fig. 6d), are well fit by our time-dependent model (Supporting Information Fig. S6). Our new solution confirms that the earthquake was an unusually shallow, low-slip event, as concluded by Ye et al. (2013) and Geirsson et al. (2015) (Fig. 6b).

Our inversion indicates that most afterslip was concentrated west-northwest of and 10-km downdip from the coseismic slip (compare Figs 6 a and c). The migration of afterslip to a location WNW of the coseismic slip is consistent with a ≈10° clockwise rotation of the post-seismic GPS station directions (red arrows in Fig. 6d) relative to the more southward-pointing coseismic directions (blue arrows in Fig. 6d). By 2017.0, 4.5 yr after the earthquake, the cumulative afterslip moment was 1.8 × 10∼²⁰ N m, 1.4 times larger than the geodetic earthquake moment. The GPS displacements attributable solely to afterslip were 2-10 times larger than the coseismic offsets (red versus blue arrows in Fig. 6d).

Our new afterslip solution fits most of the GPS station position time-series within their 1–2 mm daily noise (Figs 5 and 7 and Supporting Information Figs S6 and S7). The largest misfit occurs at site AIES in central El Salvador, where the station latitudes (north component) are systematically misfit by ≈5 mm for 1 yr after the earthquake (Fig. 8 a and Supporting Information Figs S6 and S7). In an effort to reduce this misfit, we explored numerous alternative combinations of the afterslip decay time constant, slip rake, and different degrees of smoothing, but none reduced or eliminated the misfit. We also explored whether the fit was improved if we relaxed our assumption that afterslip was stationary and instead modelled post-seismic slip with a source that migrated slowly westward along the subduction interface. We did not however find any migrating-source models that significantly improved the fit.

Figure 7.

North (a), east (b) and vertical (c) daily positions (grey circles) for GPS site SSIA reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2. Coloured curves show viscoelastic deformation solely attributable to the 2012 El Salvador earthquake for all three rheological models described in Section 3.2 and shown in Supporting Information Fig. S2. Red area in the inset map shows the 2012 El Salvador earthquake rupture zone. See Fig. 4 caption for further information.

Figure 8.

North (a), east (b) and vertical (c) daily positions (grey circles) for GPS site AIES reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2. See Fig. 7 caption for further information.

The persistent misfit to the north component of motion at AIES also occurs, though to a lesser degree, at other continuous GPS sites in southern El Salvador. In particular, a comparison of the observed motions at sites VMIG, SNJE and SSIA during the first 1–2 months after the 2012 El Salvador earthquake (shown by the shaded area in Supporting Information Fig. S7) shows that all three moved more rapidly southward during this period than predicted by our TDEFNODE afterslip model (horizontal line in Supporting Information Fig. S7). The more rapid trenchward motions at all four sites during this period might be due to viscoelastic deformation given that the predicted viscoelastic response at all these sites includes a southward component (e.g. Figs 5 a and 8a).

The afterslip decay constant estimated from our time-dependent inversion is 620 d, much longer than our 30-day estimates for the Swan Islands earthquake (Section 4.1) and Champerico earthquake (Section 4.3). A sensitivity test of the fit to the value of the decay time constant indicates that decay times as short as 300 d do not significantly degrade the fits or alter the slip solutions for this earthquake, but that shorter decay times lead to unacceptably large misfits. Afterslip triggered by the 2012 El Salvador thrust earthquake thus appears to have lasted much longer and had higher amplitude than either of the other two earthquakes that we studied.

Figs 5 and 7 compare post-seismic deformation recorded at sites SNJE and SSIA to the viscoelastic deformation calculated with VISCO-1D for the geodetic earthquake slip solution shown in Fig. 6(a). None of the three rheological models we considered (Supporting Information Fig. S3) match the observed deformation (Figs 5 and 7). The maximum-response model predicts motion of only ≈10 per cent of that observed. Our results concur with those presented by Geirsson et al. (2015), who conclude that relatively little post-seismic deformation can be attributed to viscoelastic rebound.

4.3 Champerico, Guatemala earthquake

The 2012 Champerico earthquake ruptured the Middle America trench offshore from southern Guatemala, where locking transitions from strong near Chiapas to weak offshore El Salvador (Lyon-Caen et al. 2006; Correa-Mora et al. 2009; LaFemina et al. 2009; Franco et al. 2012). Inversions by Ellis et al. (2015) of coseismic offsets and afterslip amplitudes at 19 continuous GPS sites in Guatemala, El Salvador, and southern Mexico indicate that coseismic slip of up to 2 m was concentrated in a relatively compact region above 30 km depth, immediately updip from the afterslip. Their coseismic slip estimate agrees well with an independent solution from an inversion of regional and teleseismic P waves (Ye et al. 2013).

Fig. 9 shows our new 2012 Champerico earthquake slip solution, which was derived via smoothing of the slip values estimated at nodes that approximate the Middle America subduction interface. Slip of 0.5–1.2 m was concentrated at depths of 10–30 km, close to the global CMT earthquake centroid (white star in Fig. 9a) and in good agreement with the Ellis et al. (2015) slip solution shown in Fig. 9(b). The estimated geodetic moment, 1.45 × 10²⁰ N m, is close to the 1.3 × 10²⁰ N m moment estimated by Ellis et al. (2015) and 1.33–1.5 × 10²⁰ N m seismological estimates (Ekström et al. 2012; Ye et al. 2013). The coseismic offsets recorded in the position time-series of nearby continuous GPS stations are well fit by the new TDEFNODE model (Supporting Information Fig. S8) and agree well with the static offsets estimated and inverted by Ellis et al. (2015).

Figure 9.

(a) TDEFNODE geodetic slip estimate for the 2012 Champerico earthquake. (b) Geodetic earthquake slip solution from Ellis et al. (2015). (c) TDEFNODE earthquake afterslip estimate through 2017.0. (d) Site motions predicted by the coseismic and afterslip solutions in (a) and (c). White and black stars show the earthquake centroid from the Global CMT catalogue (Dziewonski et al. 1981; Ekström et al. 2012) and the United States Geological Survey, respectively. Named GPS sites are continuous stations. Subduction depth contours are spaced every 20 km. Black circles indicate the fault nodes where slip is estimated.

Our new afterslip solution suggests that afterslip was concentrated at depths of 40–60 km (Fig. 9b), ≈20 km downdip from the region of coseismic slip. A ≈30° clockwise rotation of the post-seismic station directions compared to the coseismic offset directions (compare red to blue arrows in Fig. 9d) supports the downdip migration of the afterslip relative to the coseismic slip. By 2017.0, 4.2 yr after the earthquake, the cumulative moment for the afterslip was 1 × 10²⁰ N m, ≈70 per cent of the geodetic earthquake moment. The best-fitting decay constant is 30 d, implying that ≈25 per cent of the first year’s deformation due to afterslip accrued in the first month after the earthquake.

The TDEFNODE model fits the position time-series for nearby continuous GPS sites within their several mm daily noise (Fig. 10 and Supporting Information Figs S7 and S8). The largest misfit occurs at site COAT directly inland from the earthquake, where measurements for the first few months after the earthquake are systematically misfit in the east component by 2–5 mm (Fig. 10 b and Supporting Information Fig. S7b). We experimented with different degrees of smoothing and decay time constants to reduce the misfit, but were unable to improve the fit without concurrently degrading the fits at other continuous sites.

Figure 10.

North (a), east (b) and vertical (c) daily positions (grey circles) for GPS site COAT in southern Guatemala reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2. Coloured curves show the viscoelastic deformation that is attributable to the 2012 Champerico earthquake for all three rheological models described in Section 3.2 and shown in Supporting Information Fig. S2. Red area in the inset map shows the 2012 Champerico earthquake rupture zone. See Fig. 4 caption for further information.

Fig. 10 compares post-seismic deformation measured at site COAT (located directly onshore from the Champerico earthquake) to the viscoelastic deformation predicted by VISCO-1D and our new coseismic slip solution. The maximum-response model underestimates the observed deformation by a factor of two or more. The M1 and minimum-response estimates, which may approximate the viscoelastic deformation more accurately than the maximum-response model (Section 4.1), predict less than 10 per cent of the observed deformation. Afterslip thus appears to have been the primary cause of post-seismic deformation for this earthquake.

4.4 Far-field effects of the 2012 Nicoya Peninsula earthquake

On 2012 September 9, 9 d after the El Salvador earthquake, an M_w = 7.6 megathrust earthquake ruptured the Middle America trench at the Nicoya Peninsula in Costa Rica (Protti et al. 2013). Although we excluded this earthquake and the relevant GPS data from our TDEFNODE inversion because it falls outside our study area, we evaluated its far-field effects, which may impact some stations used in our analysis. Supporting Information Fig. S3 shows coseismic and afterslip solutions determined by Protti et al. (2013) and Malservisi et al. (2015) from inversions of the coseismic offsets and post-seismic deformation that were recorded by the dense Nicoya Peninsula GPS network and nearby sites. The maximum fault slip during the earthquake was more than 4 m and afterslip was as high as 300–500 mm on significant portions of the subduction interface by 70 d after the earthquake (Supporting Information Fig. S3b). The cumulative afterslip moment 70 d after the earthquake was ≈2 per cent of the geodetic earthquake moment.

Using the Protti et al. (2013) geodetic slip solution (Supporting Information Fig. S9a), we predicted the viscoelastic response for the 2012 Nicoya earthquake for all three Earth structures shown in Supporting Information Fig. S3. At locations on the Nicoya Peninsula directly onshore from the rupture zone, the minimum- and maximum-response models predict cumulative viscoelastic displacements as large as 200 mm by 2017.0 (not shown). The predicted deformation exceeds that for the other three earthquakes considered above (Fig. 11), as expected given that the Nicoya earthquake was the largest of the four earthquakes by a factor of two or more.

Figure 11.

Cumulative viscoelastic deformation integrated over the period 2009 May 28 to 2017.0 as predicted with VISCO-1D software and the maximum-response rheological model shown in Supporting Information Fig. S2 and described in the text. Four earthquake slip solutions were used as input to VISCO-1D, as follows: (1) the 2009 May 28 M_w = 7.3 Swan Islands strike-slip earthquake, (2) the 2012 August 27 M_w = 7.3 El Salvador earthquake, (3) the 2012 September 5 M_w = 7.6 Nicoya earthquake (Costa Rica), (4) the 2012 November 11 M_w = 7.4 Champerico (Guatemala) earthquake (Protti et al. 2013). Their respective coseismic slip solutions are shown in Figs 3(a), 6(a) and 9(a), and Supporting Information Fig. S9. The tick marks in (a) show the directions of the horizontal viscoelastic displacements. Small symbols in (b) show the GPS site locations from Fig. 2.

Relevant to our study, we compared post-seismic deformation measured at the continuous GPS station SAN0, 530 km northeast of the Nicoya earthquake rupture zone to that predicted by our viscoelastic models (Figs 11 and 12, and Supporting Information Fig. S10). The observations (grey and red/white circles in Fig. 12a) show that SAN0 accelerated southwards after the earthquake relative to its motion during the 5 yr before the 2012 Nicoya earthquake. Overall, the observations indicate that the site’s transient motion was ≈4 mm southward and ≈1 mm westward by early 2015, ≈2 yr after the earthquake (Fig. 12). For comparison, displacements predicted by the minimum- and maximum-response viscoelastic models range from 2 to 15–20 mm through 2017.0, respectively (Fig. 12; also compare Fig. 11 and Supporting Information Fig. S10). The maximum-response model, shown the red curves in Fig. 12, overestimates the observed southward displacement by a factor of four or more. The M1 or minimum-response viscoelastic models better approximate the observed deformation(Section 4.1).

Figure 12.

Time-dependent viscoelastic deformation at GPS site SAN0 (located in the previous figure) for the four earthquakes listed in the previous figure and all three rheological models shown in Supporting Information Fig. S2. The daily north (a) and east (b) station positions, shown by the grey circles, are reduced by the long-term site velocity and corrected for coseismic offsets given in Supporting Information Table S2. Coloured curves show the viscoelastic deformation individually attributable to the 2009 Swan Islands earthquake and 2012 El Salvador, Champerico, and Costa Rica subduction zones earthquakes. See Fig. 4 caption for further description.

We also compared the post-seismic deformation observed at SAN0 to the elastic deformation predicted by the 70-d afterslip solution of Malservisi et al. (2015) (Supporting Information Fig. S3). Using Malservisi et al.’s subduction fault geometry and a simple elastic-half space model, their afterslip solution predicts a ≈1–2 mm southwest-directed elastic displacement at SAN0 by 70 d after the earthquake. Best estimates of the net displacement at SAN0 from a linear regression of the daily site positions during the same 70-day period give 1 mm of southward and 0.9 mm of westward movement after correcting for the motion of the Caribbean plate. The observations are thus consistent with afterslip as the source of all the deformation recorded at SAN0 during the first 70 d after the Nicoya earthquake.

5 RESULTS: USE OF INTERSEISMIC VELOCITY FIELD TO EVALUATE SOLUTION ROBUSTNESS

A major outcome of our time-dependent TDEFNODE model is a series of daily positions for 219 GPS sites that are systematically corrected for the coseismic and afterslip effects of the 2009 Swan Islands and 2012 El Salvador and Champericoearthquakes.

We evaluated the robustness of the newly estimated station velocities and effectiveness of our time-dependent modelling via a comparison of the RMS misfits to a sample of the GPS station time-series before and after our time-dependent modelling. For example, a linear regression of the original, uncorrected position time-series for site COAT in Guatemala (Supporting Information Fig. S8) gives RMS misfits of 8.1 and 13.5 mm in the east and north components, respectively. In contrast, a linear regression of the COAT position time-series corrected for the coseismic and post-seismic effects of all three earthquakes gives RMS values of 1.9 mm in both components (Supporting Information Fig. S7), ≈80 per cent smaller than the original misfits. Similar 50–92 per cent improvements in fit are observed at nine other continuous GPS sites in our study area, all of which exhibited large transient deformations in their original position time-series. The TDEFNODE RMS misfits at these nine sites range from 1.2 to 3.1 mm (Supporting Information Fig. S7), consistent with the 1.9 mm RMS misfit atCOAT.

As a basis of comparison, RMS misfits for linear regressions of the position time-series for 910 continuous GPS stations in the North America plate interior range from 0.8 to 7 mm in the horizontal components and average 1.6 mm in both the north and east components. The 1.2–3.1 mm TDEFNODE RMS misfits are thus similar to the misfits for the mostly well-behaved plate interior GPS sites, where transient effects associated with earthquakes are negligible. On this basis, we conclude that transient deformation in the position time-series for sites in our study area is minimized effectively by our time-dependent model.

6 DISCUSSION AND CONCLUSIONS

6.1 Earthquake comparisons

Of the three earthquakes we modelled with TDEFNODE, the 2009 Swan Islands and 2012 Champerico earthquakes occurred at normal seismogenic depths and had afterslip-to-earthquake moment ratios (35 per cent and 70 per cent) typical of other large strike-slip (Freed et al. 2006) and subduction-thrust earthquakes (Lin et al. 2013). Post-seismic observations for both earthquakes are consistent with our assumptions of logarithmicallydecaying afterslip at a fixed location on the fault. Our modelling suggests that most afterslip for the 2009 Swan Islands earthquake occurred on the eastern end of the rupture zone rather than its western end (Fig. 3c), where large afterslip could have altered the stresses acting on the seismically hazardous Motaguafault.

Slip during the 2012 El Salvador earthquake averaged ≈1 m and was shallower than 20 km (Fig. 6a), consistent with previous results (Ye et al. 2013, Geirsson et al. 2015). The earthquake triggered an additional ≈1 m of nearby afterslip (Fig. 6c) with a moment equivalent to 140 per cent of the geodetic earthquake moment. Lin et al. (2013) report that afterslip following most subduction thrust earthquakes has an equivalent moment release less than 50 per cent of the coseismic moment. The 2012 El Salvador earthquake thus triggered an unusually large amount of afterslip. Given that geodetic measurements onshore from the El Salvador trench segment suggest that locking at normal seismogenic depths is either weak or zero (Correa-Mora et al. 2009), we infer that little or none of the plate convergence is accommodated by thrust earthquakes below depths of ≈25 km offshore El Salvador. Our modelling suggests that afterslip plays an important role in accommodating the plate convergence above depths of ≈25 km (Fig. 6c).

Finally, our good model fits suggest that afterslip was the process responsible for most post-seismic deformation after the 2009 Swan Islands and 2012 El Salvador and Champerico earthquakes. This finding is consistent with modelling evidence that post-seismic deformation associated with earthquakes with M_w <7.5 can be well approximated via an elastic afterslip model (Sun & Wang 2015).

6.2 Implications of far-field post-seismic deformation for western Caribbean tectonic studies

Our observations and modelling indicate that coseismic and post-seismic deformation associated with the four Central America earthquakes described above are resolvable at distances of 500 km or more from their rupture zones, encompassing most locations in Central America. By implication, estimates of plate/block rotations and interseismic fault locking in the earthquake-prone, western portion of the Caribbean plate that ignore the transient effects of regional earthquakes are likely to be biased in unpredictable ways. The good fit of our time-dependent model (Section 5 suggests that most of the transient deformation present in our original GPS position time-series was successfully removed, even given our simplifying assumption of fault afterslip as the source of all the post-seismic deformation.

Future realizations of the regional velocity field may require corrections for viscoelastic deformation given that such deformation decays more slowly (decades to centuries) than does afterslip. Fig. 11 and Supporting Information Fig. S10, which approximate the maximum and minimum summed viscoelastic responses to the 2009 Swan Islands earthquake and 2012 El Salvador, Nicoya, and Champerico earthquakes, illustrate some of the challenges and pitfalls associated with any such corrections. The magnitude and direction of a viscoelastic response predicted for a given location can vary significantly depending on the depths, viscosities, and rheological behaviours that are assigned to the lower crust and mantle for a given Earth structure. Uncertainties in the coseismic slip solutions for recent and historic earthquakes (i.e. the 1976 M_w = 7.5 Motagua fault earthquake) also cause uncertainties in the predicted viscoelastic response, particularly at sites close to earthquake rupture zones, where viscoelastic deformation is sensitive to the assumed location and slip distribution of the modelled rupture.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Figure S1: Position time-series for GPS site VMIG, which has abnormally high multipath noise. Site motion is relative to the Caribbean plate. Daily site positions are indicated by the sage-coloured circles. Red/white circles show 30-d average positions. Light-grey circles show the daily common-mode noise that has been removed from the daily positions. The red line denotes the slope that best fits the daily positions before the August 2012 El Salvador earthquake. Coseismic offsets caused by the 2009 Swan Islands earthquake and the 2012 El Salvador and Champerico (Guatemala) subduction thrust earthquakes have been removed for clarity.

Figure S2: Checkerboard tests for the Middle America trench segment that ruptured during the 2012 August 27 El Salvador earthquake, whose rupture area is defined by the red-dashed enclosed areas in panels (A) and (C). Panels (A) and (C) show starting models and their predicted (synthetic) GPS velocities for continuous (red) and campaign (black) sites. Panels (B) and (D) show slip solutions recovered from inversions of the synthetic GPS velocities and residual GPS velocities from the best-fitting solutions.

Figure S3: Diagrams of layer viscosities and depths for the three Earth structure models used for our viscoelastic modelling.

Figure S4: Horizontal fits of TDEFNODE time-dependent model at site ROA0. Shaded circles show changes in the east, north, and vertical components of the daily station positions after removing a best-fitting slope. Dashed lines denote times of the 2009 Swan Islands and 2012 El Salvador and Champerico earthquakes.

Figure S5: Residual daily GPS site positions (coloured circles) for our TDEFNODE time-dependent model. The residual daily site position is defined as the modelled position subtracted from the observed position. Dashed line denotes time of the 2009 Swan Islands earthquakes. The locations for sites AIES, SNJE, SSIA and VMIG are shown in Fig. 5(d) of the main document. Vertical signals associated with this strike-slip earthquake are small and not shown.

Figure S6: TDEFNODE fits (red and blue lines) to daily north, east, and vertical station positions reduced by their best-fitting slopes. The sites selected are continuous sites in El Salvador near the 2012 El Salvador earthquake. The locations for all four sites are shown in Fig. 6 d in the main document. Dashed lines denote times of the 2009 Swan Islands earthquake.

Figure S7: Residual daily GPS site positions (red and blue circles), defined here as the positions estimated with the TDEFNODE model subtracted from the observed positions. The shaded area shows the 72-d period between the 2012 El Salvador and Champerico earthquakes, whose times are indicated by the dashed lines. Locations for sites BARI, CHPO, COTZ, COAT and MTP1 are variously shown in Figs 9(a) and (d). RMS values for the TDEFNODE misfits include station positions for all time spanned by the data for each station, as specified in Table S1. The fits to the vertical station positions are not shown.

Figure S8: TDEFNODE fits (red and blue lines) to daily north, east, and vertical station positions reduced by their best-fitting slopes. The sites selected are continuous sites in Guatemala near the 2012 Champerico (Guatemala) earthquake. The locations for all five sites are shown in Fig. 8 of the main document. Dashed lines denote times of the 2009 Swan Islands earthquake and 2012 El Salvador and Champerico earthquakes.

Figure S9: (a) The geodetic coseismic slip solution of Protti et al. (2013) for the September 2012 M = 7.6 Nicoya earthquake. (b) The geodetic afterslip solution for the same earthquake, representing the cumulative afterslip for the first 70 d after the earthquake (Malservisi et al. 2015).

Figure S10: Cumulative viscoelastic deformation integrated over the period 2009 May 28 to 2017.0 as predicted with VISCO-1D software and the minimum-response rheological model shown in Fig. S2 and described in the text.

Table S1: GPS site information.

Table S2: Estimated coseismic offsets.

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

AE was supported by a UNAVCO Graduate COCONet Fellowship. We thank Roland Burgmann, Laura Wallace and the associate editor for constructive suggestions. This work was funded by the NSF grant EAR-1144418 (DeMets) and also benefited from the National Science Foundation support for UNAVCO under grant EAR-1042906. We thank Central American institutions that provided technical, data, and logistical support, including Universidad Nacional Autonoma de Honduras, Universidad Nacional Autonoma de Mexico, Universidad Mariano Galvez (Guatemala City), Universidad San Carlos (Guatemala City), Ministerios de Medio Ambiente y Recursos Naturales El Salvador, Instituto Geografico Nacional of El Salvador and Guatemala, INSIVUMEH of Guatemala, and CONRED of Guatemala. Figures were prepared using Generic Mapping Tools software (Wessel & Smith 1991).