A geodetic study of the 2003 January 22 Tecomán, Colima, Mexico earthquake

Summary

We use coseismic displacements and aftershock information from Global Positioning System (GPS) measurements at 27 sites in western Mexico and a 12-station local seismic network to determine the characteristics of the 2003 January 22 M_w = 7.2 subduction thrust earthquake near Tecomán, Colima, Mexico. Estimates of the earthquake moment, slip direction and best-fitting slip distribution are derived by optimizing the fit to the GPS displacements for a 3-D finite element mesh that simulates the study area. The calculated moment release is 9.1 × 10¹⁹ N m (M_w = 7.2), with maximum slip of 2 m at a depth of 24 km and a maximum rupture depth of 35–40 km. The inversion indicates that coseismic rupture extended downdip from depths of 9 to 40 km along a 80 km along-strike region that is bounded by the edges of the Manzanillo Trough. The optimal solution is robust with respect to plausible changes in the subduction interface geometry and differing subsets of the data. A comparison of the cumulative post-seismic slip that can be inferred separately from earthquake aftershocks and GPS measurements within a year of the earthquake indicates that 95 per cent or more of the post-seismic deformation was aseismic. Near-term post-seismic measurements indicate that slip propagated downdip to areas of the subduction interface beneath the coastline within days following the earthquake, as also occurred after the nearby M_w = 8.0 Colima-Jalisco subduction earthquake in 1995. The similar behaviours and locations of the 1995/2003 earthquake sequence to two earthquakes in June of 1932 suggests that thrust earthquakes along the subduction interface northwest of the Manzanillo Trough may trigger earthquakes in the vicinity of the Manzanillo Trough; however, our modelling of Coulomb stress changes caused by the 1995 earthquake indicate that it induced only modest unclamping of the subduction interface in the vicinity of the Tecomán rupture. In addition, GPS measurements indicate that elastic shortening characterized areas onshore from the Tecomán rupture from mid-1997 up until the time of the rupture, consistent with progressively stronger clamping of the subduction interface during this period. This precludes any obvious triggering relationship with the 1995 earthquake. The apparent coincidence of the edge of both the 1932 and 1995/2003 rupture sequences with the edge of the Manzanillo Trough may indicate that the trough is a mechanical barrier to along-strike rupture propagation. This implies a limit to the area of potential slip and hence rupture magnitude during future large earthquakes in this region.

earthquake triggering, Middle America trench, Rivera plate, subduction thrust earthquakes

1 Introduction

As the predominant source of destructive earthquakes in much of Mexico and Central America, the Middle America subduction zone has been the target of numerous seismologic studies since the 1960s. To better understand seismic hazards in this region, seismologic studies of the epicentral parameters and rupture areas of large earthquakes along the Middle America trench (Dean & Drake 1978; Chael & Stewart 1982; Burbach et al. 1984; Singh et al. 1984, 1985; Anderson et al. 1989; Pardo & Suarez 1995), possible seismic gaps (Singh et al. 1981; Nishenko & Singh 1987), and the budget of seismic to aseismic slip (McNally & Minster 1981; Pacheco et al. 1993) have all been undertaken. Global Positioning System (GPS) measurements of crustal displacements and strain are now also being used to study the seismic cycle in southern Mexico, with particular emphasis on interseismic strain accumulation (Yoshioka et al. 2004), aseismic strain transients (Lowry et al. 2001; Kostoglodov et al. 2003) and transient post-seismic processes (Hutton et al. 2001; Márquez-Azúa et al. 2002), all of which are largely invisible to seismometers.

The most recent large earthquake to rupture the Middle America subduction interface was the 2003 January 22 Tecomán subduction earthquake offshore from Tecomán, Colima, Mexico (Fig. 1). With an estimated magnitude M_w between 7.2 and 7.6, the Tecomán earthquake caused the deaths or injuries of more than 1000 people and damaged or destroyed more than 40 000 homes and businesses, most in and near the inland city of Colima. Focal mechanisms for this earthquake indicate that it accommodated shallow underthrusting of the offshore oceanic lithosphere beneath the continental margin (Fig. 1; Ekström et al. 2004; Yagi et al. 2004). Inversions of local and teleseismic body waves suggest that the rupture initiated at a depth of 15–25 km and propagated both up and downdip (Yagi et al. 2004), with most slip concentrated in patches immediately updip and downdip from the hypocentre. Modelling of seismic waveforms (Yagi et al. 2004) also suggests that most and possibly all of the coseismic slip was limited along-strike to areas beneath the Manzanillo Trough and southern Colima graben, which are sediment filled structures that appear to accommodate trench-parallel extension (Bourgois et al. 1988; Bandy et al. 1995; Michaud et al. 2000).

Figure 1.

Upper: Location and tectonic map of study area. Grey arrows and 2-D, 1σ ellipses depict Cocos and Rivera plate velocities relative to North American plate for the past 0.78 Ma (DeMets & Wilson 1997). Open stars show epicentres for the 1995 Colima-Jalisco and 2003 Tecomán earthquakes (Pacheco et al. 1997; Yagi et al. 2004). Earthquake focal mechanisms are from the Harvard centroid-moment tensor catalogue. Solid circles show locations of GPS sites used in this study. Lower: North and east components of continuous GPS coordinate time-series for site COLI, located onshore from the 1995 and 2003 rupture zones. Non-linear post-seismic site motions represent a superposition of steady elastic strain accumulation from locked portions of the Middle America subduction zone and decaying afterslip and viscoelastic flow triggered by the two earthquakes (Márquez-Azúa et al. 2002).

Below, we use measurements from 27 GPS sites (Fig. 1) and a local short-period seismic network to study several different aspects of the Tecomán earthquake and its near-term post-seismic response. We first employ GPS measurements, including coseismic offsets for 12 sites within 100 km of the rupture zone, to describe the location, distribution and magnitude of coseismic slip on the subduction interface. Modelling of the coseismic offsets is done using the finite element method (FEM) and 3-D mesh that simulates the assumed geometry of the subduction interface and first-order spatial variations in elastic properties of the crust and mantle in our study area. Aftershocks recorded by a local seismic network are used to justify the range of likely geometries for the subduction interface and provide independent evidence for the extent of the rupture zone. We compare our geodetically based solution to the seismologically derived coseismic slip distribution of Yagi et al. (2004) to demonstrate that the latter likely overestimates the downdip extent of coseismic rupture. We use continuous and quasi-continuous GPS measurements from eight regional stations to quantify the post-seismic response, including reversals in the vertical motions at two coastal sites that are consistent with the rapid onset of fault afterslip immediately downdip from the rupture zone. We also determine Coulomb stress changes that were induced by the M_w = 8.0, 1995 October 9 Colima-Jalisco earthquake, which ruptured the subduction zone immediately northwest of the Tecomán earthquake (Courboulex et al. 1997), to investigate a possible triggering relationship between the two earthquakes. Finally, we compare the Tecomán earthquake rupture zone to the rupture zones of previous large earthquakes in this region to search for patterns and differences possibly useful for furthering our understanding of the seismic cycle in this area.

2 Data

2.1 GPS observations

The Jalisco GPS network consists of 31 geodetic markers that were installed in the 1990s and have been occupied annually since early 1995. By chance, we initiated our network occupation in 2003 approximately one day before the Tecomán earthquake occurred offshore from our GPS network. Eight dual-frequency GPS receivers were deployed when the earthquake struck, half at sites that were relatively close to the rupture zone. Measurements at these eight sites give coseismic offsets that are relatively uncontaminated by post-seismic effects. We occupied an additional 18 sites within 1 week of the earthquake and one site (PORT) 5 months after the earthquake. The offsets measured at these 19 sites are thus contaminated by varying degrees of post-seismic motion, although we later describe evidence that post-seismic slip at these sites is unlikely to be more than 5–10 per cent of the coseismic motion over the time frame that we made our measurements.

All GPS code-phase measurements were analysed using GIPSY analysis software from the Jet Propulsion Laboratory (JPL). We employed a standard point-positioning analysis strategy (Zumberge et al. 1997) combined with resolution of integer phase ambiguities. Daily GPS station coordinates were first estimated in a non-fiducial reference frame (Heflin et al. 1992) employing precise satellite orbits and clocks from JPL, and were then transformed to ITRF2000 (Altamimi et al. 2002) using daily seven-parameter Helmert transformations supplied by JPL. We also estimated and removed daily and longer-period spatially correlated noise between sites using a regional-scale noise stacking technique first described by Wdowinski et al. (1997) and extended by Márquez-Azúa & DeMets (2003) to estimate and remove plate-scale, common-mode noise from campaign and continuous GPS station coordinate time-series. After removing common-mode noise, the daily site coordinate repeatabilities are 2–4, 3–5 and 8–10 mm in the north, east and down components, respectively. These are 10–40 per cent smaller than the daily coordinate repeatabilities for the uncorrected daily site coordinates. The common-mode corrections also effectively reduced longer-period noise in the GPS time-series, typically by 50–70 per cent in amplitude.

For the eight sites that were operating during the earthquake, coseismic offsets are estimated by subtracting the GPS site location derived from all pre-earthquake data from the site location determined from the first 24 hr of GPS data after the earthquake. Logistical problems at site COLI, where power was lost for 36 hr after the earthquake, and site AVAL, where a toppled antenna was not reset until 16 hr after the earthquake, preclude precise determinations of uncontaminated coseismic offsets at those locations.

The offsets for sites that were not occupied until one or more days after the earthquake and had not been occupied for as long as 2 yr before the earthquake were estimated by comparing site locations extrapolated (via linear regression) from the pre-earthquake site coordinate time-series (Figs 2–4) to the site coordinates determined from the first 24 hr of data collected after the earthquake.

Figure 2.

GPS coordinate time-series used to find the 19 quasi-coseismic offsets. North–south movement observed since 10/95. Each dot shows the relative site location for a measurement session of 12–24 h duration. North American plate motion has been removed. The linear regression used to determine pre-earthquake site motion is shown by the medium-thickness lines. Thick vertical lines represent the quasi-coseismic offsets. Horizontal medium-thickness lines after the earthquake illustrate how site position after the earthquake is obtained from the first available GPS solution. Thin lines connect observations that are not included in the calculation of quasi-coseismic displacements in order to make each site's position data visually coherent. Timescale differs before and after the earthquake.

Figure 3.

GPS time-series showing east–west site movement with North American plate motion removed. Site motion and quasi-coseismic offsets are described in Fig. 2.

Figure 4.

GPS time-series showing vertical site movements. Site motion and quasi-coseismic offsets are described in Fig. 2. Vertical scale differs from that for Figs 2–3.

Uncertainties for the eight site offsets that are based on data collected immediately before and after (or soon after) the earthquake are ±4, ±7 and ±11 mm for the north, east and vertical components based on the NED daily scatter reported above. The remaining 19 offsets include varying amounts of post-seismic movement and also depend on extrapolation to estimate the pre-earthquake site coordinates. We thus assign larger uncertainties to these ‘quasi-coseismic’ offsets: ±6, ±10 and ±16 mm for the north, east and down components. One site (PORT) that was not occupied until 5 months after the earthquake is assigned even larger uncertainties of ±9, ±15 and ±23 mm for the NED components. One highly anomalous vertical displacement at site MILN, most likely due to an erroneously measured antenna height, was downweighted to reduce its influence on the solution. A second anomalous vertical offset, that at CRIP, is likely a consequence of a shorter observation session at that site the day after the Tecomán earthquake.

Prior to estimating the coseismic offsets, all site motions were transformed into a reference frame fixed to the North American plate (Márquez-Azúa & DeMets 2003) in order to remove the contribution of steady plate movement to each site's motion before and after the earthquake. The coseismic displacements of all 27 sites are illustrated in Fig. 5 and given in Table 1.

Figure 5.

(a) Horizontal coseismic GPS site displacements for the Tecomán earthquake and their associated 2D, 1σ uncertainties. Measurements from sites occupied during the earthquake are shown in black; measurements calculated from GPS time-series (‘quasi-coseismic’) are shown by open arrows. (b) Vertical displacements. Solid and open bars indicate pure coseismic and quasi-coseismic offsets, respectively.

Table 1.

Observed and modelled coseismic GPS site offsets.

2.2 Seismic observations

We used seismograms from a 12-station, short-period seismic network (RESCO) operated by the University of Colima (Fig. 6) to locate hypocentres from first arrivals of P and S waves to learn more about the limits of the 2003 rupture along the subduction interface (Fig. 6b), the geometry of the subduction interface beneath our study area (Fig. 6a), and the amount of post-seismic deformation that was accommodated by seismic processes (Fig. 6c). Aftershocks recorded within the first 30 d after the earthquake (shown by red hypocentres in Fig. 6) are largely confined to offshore areas beneath the Manzanillo Trough (Fig. 6b). If the aftershock epicentres are reliable, they suggest that the Tecomán earthquake largely ruptured the subduction interface beneath the Manzanillo Trough, in accord with the seismologically constrained coseismic slip distribution of Yagi et al. (2004).

Figure 6.

(a) Trench-normal cross-section of aftershocks detected by RESCO seismic network. The hypocentre for the main shock is from Yagi et al. (2004). Three different model subduction interfaces are shown. Our preferred geometry (bold line), adopted from Hutton et al. (2001), implies a depth of 20.4 km for the Tecomán earthquake. End-member models with shallower and more steeply dipping slabs (dashed lines) imply respective hypocentral depths of 15 and 25 km. (b) Map view of RESCO aftershocks (coloured circles) and seismometer locations (yellow triangles). (c) Histogram showing number of earthquakes per day, with best fit based on modified Omori's law shown by the red line (discussed in text).

Surprisingly, nearly all of the aftershocks are located at depths shallower than the subduction interface that we adopted in our mesh (Fig. 6a). The unrelocated hypocentres are thus less useful for defining the location and dip of the subduction interface than we hoped. Several explanations for the unexpectedly shallow aftershock depths are possible. The most likely in our view is that the aftershock depths are systematically in error. Alternatively, the subduction interface may be significantly shallower beneath the Manzanillo trough and southern Colima graben than is assumed in our preferred mesh, or the Tecomán earthquake may have triggered numerous smaller earthquakes within the overlying plate. Testing these alternative hypotheses will require relocating the aftershocks, which is beyond the scope of this analysis.

3 Methods

Our principal objective is to estimate the distribution of coseismic slip along the subduction interface beneath and offshore from our study area. Rather than employing elastic half-space modelling and assuming homogeneous elastic properties throughout the half-space, we instead opted to use finite element modelling so as to better simulate the geometry of the subduction interface and the likely large-scale lateral and vertical variations in the elastic properties of the crust and mantle. Below, we describe the principal attributes of the finite element mesh. Details of the inverse procedure we employ to estimate a best-fitting distribution of coseismic slip are given in Appendix A.

3.1 Finite element mesh

We use the commercially available finite element package ABAQUS to model the region. Our 3-D finite element mesh (Fig. 7) employs spherical coordinates, a spherical earth geometry, and is designed to approximate the geometrical characteristics of the northern Middle America subduction zone. The mesh includes 6 km-thick oceanic crust and continental crust divided into a 16 km-thick elastic upper crust, a 24 km-thick viscoelastic lower crust, and viscoelastic mantle beneath both to a depth of 250 km. Variations in the rigidity of the upper crust seem likely given the existence of the sediment-filled southern Colima graben within our GPS field area; however, we are unaware of any tomographic studies of the upper crust in this region that could be used to prescribe variations in the elastic properties of the upper crust within our finite element mesh. We thus model the upper crust as homogeneous. The model configuration and material properties (Table 2) share many similarities to those employed by Masterlark (2003).

Figure 7.

Exploded view of finite element mesh employed for modelling. Values of Young's modulus E are shown. Coseismic deformation is modelled as a static step in the FEM, so the entire model behaves elastically. Other material properties are given in Table 2.

Table 2.

Material properties assumed for finite element model.

The subduction interface is represented as a surface of contact node pairs (Fig. 8) that replicate the interface geometry derived by Hutton et al. (2001) from GPS and seismologic data (Pardo & Suarez 1993, 1995). Node spacing is 10 km along strike and varies downdip. The subduction interface we use agrees well with aftershock locations for the Tecomán earthquake (Fig. 8), particularly above depths of 30–35 km where most slip is focused. The influence of alternative geometries for the subduction interface on our modelling results is described in a later section.

Figure 8.

Cross-section of finite element mesh. RESCO aftershock locations that delineate the Benioff zone correspond well with the assumed location of the subduction interface above depths of ∼40 km, coinciding with the maximum extent of seismogenic rupture typical of the Middle America subduction zone (Suárez & Sánchez 1996).

4 Results

4.1 Best-fitting coseismic slip distribution

We determined the best coseismic slip distribution by systematically inverting the 27 GPS offsets to identify the smoothing weight ε that minimizes χ²_ν. During early iterations of the data, we allowed slip to occur along an area of the subduction interface substantially larger than indicated by the distribution of aftershocks. Later iterations eliminated large areas of zero fault slip from the problem domain, thereby improving both the computational efficiency and resolution of our model in areas where slip occurred. Both the fit to the GPS offsets and best-fitting distribution of coseismic slip along the subduction interface (Fig. 9a) remain approximately constant for a broad range of smoothing weights, indicating that our best-fitting model is robust with respect to the values we adopt for the smoothing weight.

Figure 9.

(a) Best-fitting coseismic slip distribution and observed (solid arrows) and predicted (red arrows) horizontal coseismic site displacements. The grey rectangular border defines the region of allowed fault slip. (b) Residuals (observed minus predicted displacements) for the best-fitting model. (c) Observed (black) and predicted (red) vertical displacements during the Tecomán earthquake. (d) Predicted (red) and observed (black) vertical displacements for sites inland from the rupture area.

The best-fitting slip distribution (Fig. 9a) is an approximately elliptical region located offshore from the southern Colima graben and focused on the RESCO earthquake epicentre. Most of the slip along the subduction interface is concentrated between depths of 9 km and 35–40 km, with maximum slip of 2 m at a depth of 24 km. For the shear moduli in Table 2, the calculated moment for this slip distribution is 9.1 × 10¹⁹ N m, corresponding to M_w = 7.2 (Table 3). Slip does not appear to extend beyond the edges of the Manzanillo Trough, suggesting that the faults that bound this structure played a role in limiting the rupture.

Table 3.

Comparative earthquake sizes and depths.

The average slip direction of the best-fitting slip distribution is N28°E, about 1° clockwise from the slip direction derived by Yagi et al. (2004) from their inversion of local seismic and teleseismic data. Given the good agreement between these two estimates, we conclude that slip direction is well determined. Our own and Yagi et al. (2004) best estimates for the slip direction lie 9°–10° counter-clockwise from the dip direction, indicating that the earthquake accommodated oblique downdip slip. Our estimates are also rotated 7–10° counter-clockwise from the predicted convergence directions for the Rivera-North America (N36°E±2°) and Cocos-North America (N33°E±2°) plate pairs (DeMets & Wilson 1997). The Tecomán earthquake thus does not help to resolve whether lithosphere subducting beneath the Manzanillo Trough moves with the Cocos or Rivera plates.

Weighted rms misfits for the best-fitting slip distribution are 7, 8 and 22 mm in the north, east and vertical components, respectively. χ²_ν for the best solution is 2.13, which indicates that the data are misfit by 1.5 times their estimated uncertainties. The predicted patterns of horizontal and vertical coseismic offsets capture the first-order features of the surface deformation defined by our observations (Fig. 9b–d and Table 1). These include the rapid decrease in the magnitudes of the site displacements with distance from the earthquake epicentre and the rapid variation in displacement directions at coastal locations near the epicentre. Most notably, the model does a relatively good job of fitting the coseismic offset for site SJDL, whose coseismic motion points ∼90° counter-clockwise from the coseismic offsets exhibited at other nearby sites (Figs 5 and 9a).

The best-fitting slip distribution predicts surface displacements that systematically misfit the observations in two areas. The model predicts coseismic vertical motions substantially different than observed at all three sites closest to the rupture zone (MIRA, CRIP and SJDL in Fig. 9d). For example, the model predicts 6 mm of uplift at site MIRA, whereas 24 mm of subsidence was measured. That the vertical motions at these sites are fit poorly despite their high formal data importance (Table 1 and Fig. 10) suggests that one or more of our modelling assumptions could be wrong. For example, if sympathetic slip along one or more crustal faults was triggered by the earthquake, our assumption that all slip was focused along the subduction interface would be incorrect. Alternatively, lateral variations in material properties not accounted for in our model could influence the model predictions. The misfit could also result from a suboptimal finite element mesh geometry for the subduction interface, although this source of misfit seems unlikely (see Section 4.3.1).

Figure 10.

Graphical depiction of model resolution m_r and data importance n_r. Model resolution is shown in a surface projection of the model area; slip is better resolved in darker regions. Data importance are expressed by the vertical bars; the length of each bar corresponds to the sum of the E, N and D data importance for each GPS site. The sum of all data importance is 18.7.

The second coherent pattern in the model residuals consists of east-directed residual offsets (Fig. 9b) for inland sites on the Jalisco block (AUTA, AYUT, GUAC, LIM2, SEBA and TAPA, see Fig. 5). The low formal data importance for the offsets at these sites (Fig. 10) indicate that they contribute a negligible amount of information to the best-fitting slip solution and influence it minimally. Possible reasons for this misfit are discussed and tested in Section 5.1.

4.2 Model resolution and data importance: constraints on downdip rupture limit

Fig. 10 shows the model resolution m_r and relative data importance n_r/ǀn_rǀ, which are defined in Appendix A. As expected, sites located closer to the rupture zone contribute the most information about the location and magnitude of slip. The data from sites AVAL, COOB and COLI, though more distant from the rupture zone, also contribute significant information to the model because of their low uncertainties. Overall, 95 per cent of the model information is contributed by displacements from 10 of the 27 sites (Table 1). All ten of these are located within 100 km of the earthquake epicentre.

The model resolution is nearly zero in the shallowest portions of the subduction interface (Fig. 10), as is typical for GPS networks in subduction zone settings. Any coseismic slip that occurred updip from depths of ∼5 km is thus not detectable with our network. The model resolutions are significantly better from depths of 10–60 km. The apparent absence of slip below depths of ∼35–40 km in our best slip distribution is thus a well-resolved feature of our model, indicating that the GPS sites sample the surface deformation in a manner that imposes strong constraints on the downdip limit of coseismic slip.

Simple elastic models of thrust earthquakes predict that the transition or hinge line that separates areas of coseismic uplift from areas of coseismic subsidence should occur approximately at the location on the surface that lies above the downdip limit of coseismic rupture. That all of the coastal sites onshore from the Tecomán rupture zone either subsided or exhibited insignificant vertical motion during the earthquake (see vertical motions of sites MIRA, CRIP, NOVI, SJDL and UCOL in Fig. 9d) implies that the rupture extended down the subduction interface no farther inland than the Pacific coast and more likely somewhere offshore. For the interface geometry that we employ (Fig. 8), this implies an approximate maximum rupture depth of 30–35 km, approximately the same as is implied by our inverse modelling. This result is robust with respect to inversions of various subsets of our data.

4.3 Model sensitivity

4.3.1 Influence of subduction interface geometry

The imperfectly known geometry of the subduction interface adds uncertainty to our modelling results, particularly since the geometry we use (Figs 6 and 8) was optimized to fit geodetic observations from the 1995 October 9 Colima-Jalisco earthquake, which ruptured the subduction interface northwest of the Manzanillo Trough (Hutton et al. 2001). Since it is computationally impractical to explore numerous alternative mesh geometries, we instead examined the fits and slip distributions for two geometries (indicated by dashed lines in Fig. 6) that bracket the 15–25 km range of hypocentral depths for which seismograms from the Tecomán earthquake are well fit by synthetic seismograms (Yagi et al. 2004).

Our inversion of the 27 coseismic offsets that employs the steeper subduction interface yields a rupture area nearly identical to that for the best-fitting slip distribution (Fig. 11a), a modestly larger (17 per cent) seismic moment, and a weighted misfit that is indistinguishable from that for our best-fitting solution. Neither the slip distribution nor the misfits for our best-fitting model change measurably if we instead use a steeper subduction interface geometry. In contrast, our inversion of the data that employs a finite element mesh with a shallower subduction interface yields a model misfit that is ∼40 per cent larger than for our best-fitting model. The shallower geometry thus yields a coseismic slip model that fits the data more poorly. The data misfits for our best-fitting model (Fig. 9), therefore, are not reduced significantly by varying the geometry of the subduction interface within plausible bounds suggested by the available seismic data.

Figure 11.

Slip distributions and predictions for finite element meshes with (a) deeper and (b) shallower subduction interfaces. Both models predict coastal uplift immediately onshore from the rupture, in disagreement with the observations. Observed and predicted coseismic movements are shown with black and red colour, respectively.

4.3.2 Influence of quasi-coseismic data

Nineteen of the 27 GPS stations where we measured coseismic displacements were not occupied until days to weeks after the earthquake. This renders our measurements of their ‘coseismic’ offsets susceptible to any transient post-seismic deformation that might have occurred. The continuous GPS coordinate time-series at site COLI (Fig. 1) clearly shows that Tecomán earthquake triggered transient post-seismic deformation, as does the GPS time-series for the site closest to the rupture zone (MIRA). After the earthquake, both COLI and MIRA moved southwestward towards the rupture zone at rates that decayed rapidly with time (Figs 1–3). At COLI, the post-seismic movement within one day of the earthquake equalled ∼5 per cent of the coseismic offset at this site and increased thereafter to 10–15 per cent 1 week after the earthquake and 20–25 per cent 1 month after the earthquake. Given the directions and sense of the measured post-seismic and coseismic motions at both sites (i.e. southwestward), any estimate of the coseismic offset that is based on measurements first made days to weeks after the earthquake would yield an apparent coseismic offset that is greater than the actual coseismic offset.

The observations at COLI and MIRA indicate that the horizontal offsets we measured at the 19 sites that we were unable to reoccupy until days to weeks after the earthquake are systematically larger than the true coseismic offsets at these sites, most likely by amounts that range from 5 to 15 per cent depending on when a given site was first occupied after the earthquake. We therefore, reduced all 19 of these quasi-coseismic offsets by amounts of 5–15 per cent, depending on when a given site was first occupied after the earthquake. We then re-inverted the 19 adjusted offsets and eight offsets determined for the sites that were operating during the earthquake to derive new optimal-fit models for the coseismic slip distribution. This yielded a modified slip model (not shown) that closely resembled our best-fitting coseismic slip solution. The modified model failed to improve significantly the fit relative to that of our best-fitting slip distribution and at some sites, yielded worse fits.

The major characteristics of our best-fitting coseismic slip distribution (Fig. 9), are therefore, robust with respect to likely biases in our data from uncorrected post-seismic motion. Moreover, the maximum upward bias in the seismic moment that we estimate from our inversion of the pure and quasi-coseismic offsets listed in Table 1 is likely to be no greater than 10 per cent and is more likely less than 5 per cent given our strategy of downweighting all of the quasi-coseismic offsets to reduce their influence on the best-fitting solution.

4.4 Comparison to seismological solution

Yagi et al. (2004) invert local strong-motion observations and teleseismic data to estimate the spatial distribution of coseismic slip during the Tecomán earthquake. Their rupture model is dominated by two slip patches connected by a narrow saddle of lower magnitude slip where the earthquake nucleated (Fig. 12). Maximum slip in their model is 3.4 m, with an overall earthquake moment of 2.3 × 10²⁰ Nm (M_w = 7.5). Assuming that the cumulative slip within the first few days of the earthquake was dominated by the coseismic rupture, as is suggested by the modest near-term post-seismic transient at site COLI (Fig. 1), the geodetically based and seismologically derived slip solutions should be similar in their locations, magnitudes and slip directions.

Figure 12.

Predictions of coseismic surface displacements at GPS sites using the seismologically based slip distribution from Yagi et al. (2004) and FEM modified to mimic elastic properties employed by Yagi et al. (2004). (a) Observed (black) and predicted (red) horizontal coseismic site displacements. The grey rectangular border defines the region of allowed fault slip. (b) Residuals (observed minus predicted displacements) for the best-fitting model. (c) Observed (black) and predicted (red) vertical displacements during the Tecomán earthquake. (d) Predicted (red) and observed (black) vertical displacements for sites inland from the rupture area.

We tested whether the seismic slip distribution is compatible with the measured GPS offsets by sampling the seismic slip distribution into our finite element mesh and using it to predict 3-D surface displacements. To ensure that the predicted displacements accurately reflect those implied by Yagi et al. (2004), we also modified our mesh properties to incorporate shear moduli derived from their shear velocity model (Table 4). Since rocks usually have larger shear moduli under dynamic (seismic) conditions than under static (geodetic) conditions (Hooper et al. 2002), we divide the shear modulus obtained from V_s by 1.1, following the relationship between dynamic and static moduli presented for appropriate rock types in Lama & Vutukuri (1978). Relative to our observed coseismic offsets, the rms misfits of the surface displacements predicted by the Yagi et al. (2004) model are 13, 16 and 47 mm in the north, east and vertical components (Fig. 12), respectively. These are 84, 107 and 118 per cent larger than the misfits for our best-fitting slip distribution, respectively. Seismic slip in the Yagi et al. (2004) model, as sampled into our mesh, extends tens of kilometers inland from the coast to depths of nearly 50 km. This moves the hinge line that separates regions of coseismic uplift and subsidence far inland relative to our best geodetic model (Figs 9d and 12d). The deep slip causes large misfits (>10σ) to the observed vertical motions of the coastal sites. Given these large misfits to the observed surface displacements, the seismologic evidence for coseismic slip as deep as 50 km may be an artefact of the seismological data or may reflect a fundamental limitation in using the available seismic data to establish the downdip limit to rupture.

Table 4.

Shear moduli calculated from stated V_s and density from Yagi et al. (2004).

As is summarized in Table 3, the geodetically derived earthquake moment of 9.1 × 10¹⁹ N m (M_w = 7.2) is smaller by a factor of 2.2–2.5 than the seismologically derived moments of 2.3 × 10²⁰ N m (M_w = 7.5) reported by Yagi et al. (2004) and by Ekström et al. (2004)(2.0 × 10²⁰ N m, M_w = 7.5). Since the elastic properties assumed by these studies differ, a more useful comparison is the seismic potency (Ben-Menahem & Singh 1981), which is defined as , where μ is the shear modulus and is the slip integrated over the rupture area. The potency of our best-fitting slip distribution is 5.1 × 10⁹ m³, relatively close to the potency of 4.7 × 10⁹ m³ for the slip distribution of Yagi et al. (2004). We conclude that the earthquake magnitude is well determined. Ekström et al. (2004) do not provide potency or shear modulus estimates.

4.5 Post-seismic response

4.5.1 Seismic versus aseismic slip

We fit the available RESCO aftershocks using the modified Omori's law (Fig. 6c) and found that aftershock occurrence followed a typical power-law decay, with p = 1.045. Like past earthquakes on the northern Middle America subduction zone, the total number of aftershocks was relatively small (Singh et al. 2003). Based on the local magnitudes of the aftershocks, an upper bound on the cumulative moment release of aftershocks that occurred within 6 months of the main shock is ∼2 × 10¹⁸ N m, equal to 2–3 per cent of the moment release of the main shock. Assuming that these aftershocks released slip over an area approximately equal in size to the area of coseismic rupture (∼2500 km²), the integrated elastic response at the surface from these aftershocks is only a few percent of the coseismic displacements.

The nearly complete GPS coordinate time-series for site COLI (Figs 1 and 14) offers a useful independent estimate of the magnitude of post-seismic slip along the subduction interface after the Tecomán earthquake. Site COLI moved 127 ± 8 mm towards S34°W during the Tecomán earthquake. During the first year after the earthquake, COLI moved an additional 50 mm towards S27°W, the same direction within the uncertainties as the coseismic motion and equal to 40 per cent of the coseismic offset. Post-seismic motion at the nearby continuous site COOB (Fig. 13) was also directed towards S30°W and exhibited a ratio of post-seismic to coseismic movement of ∼50 per cent 1 yr after the earthquake. The 40–50 per cent ratio of post-seismic-to-coseismic motion observed at the two sites exceeds by more than an order-of-magnitude the small elastic displacement that is predicted to have occurred at these sites based on the cumulative moment release for the RESCO aftershocks for the same period.

Figure 14.

Best-fitting slip distributions from the 1995 and 2003 earthquakes from inversions of the 1995 and 2003 GPS site displacements using identical FEMs.

Figure 13.

North motions of continuously operating GPS sites in the study area and post-seismic response to the 1995 and 2003 earthquakes. All site motions are corrected for the steady motion of the North American plate. Motions before October 1995 at sites COLI and INEG represent steady interseismic strain accumulation, in contrast to transient motion from 1995 October to 2003 January caused by a combination of transient afterslip downdip from the 1995 earthquake rupture zone, transient viscoelastic flow triggered by the earthquake and strain accumulation from relocked, seismogenic areas of the subduction zone (Márquez-Azúa et al. 2002).

The aftershock and post-seismic GPS observations, therefore, indicate that ∼95 per cent of the post-seismic movement during the year after the earthquake was aseismic. In addition, the nearly identical coseismic and post-seismic directions of motion at COLI suggest that afterslip was focused along or directly downdip from the coseismic rupture zone (see below). In contrast, the directions of coseismic and post-seismic site movements at site COLI during and after the 1995 Colima-Jalisco earthquake differed by 44° (Márquez-Azúa et al. 2002) and by tens of degrees at other GPS sites in western Mexico (Hutton et al. 2001), consistent with focusing of afterslip triggered by that earthquake along areas of the subduction interface that differed significantly from its coseismic rupture area.

4.5.2 Evidence for rapid downdip migration of post-seismic slip

Coastal sites CRIP and MIRA both experienced sustained uplift in the days following the Tecomán earthquake (Fig. 4), in contrast to the downward motion of both sites during the earthquake. The reversal in the sense of vertical motion at these two sites repeats identical behaviour at site CRIP during and after the 1995 Colima-Jalisco earthquake (Melbourne et al. 2002), which was interpreted by Melbourne et al. as evidence that significant earthquake afterslip occurred along the subduction interface directly downdip from the offshore coseismic rupture zone. We concur with this interpretation. Our observations in 1995 and again in 2003 show that coseismic subsidence at coastal sites onshore from a subduction zone can be rapidly offset by post-seismic uplift, raising concerns about the usefulness of estimates of coseismic vertical motion that are based on measurements first made days to months after an earthquake.

5 Discussion

5.1 Model misfits for Jalisco block sites

The spatially coherent residual motions for sites on the Jalisco block are a robust feature of our best-fitting (Fig. 9b) and alternative slip solutions (e.g. Figs 11–12). Possible explanations include the following: (1) the offsets for these sites contain intersite correlated errors, (2) one or more modelling assumptions regarding coseismic slip or the elastic properties of the crust are incorrect or (3) an unrecognized tectonic effect causes the misfit. We briefly discuss each below.

Correlated errors in our estimated offsets at inland sites could be introduced by the linear-motion assumption we employed for estimating coseismic offsets at those sites. None of the six inland sites with systematically misfit coseismic offsets had been occupied for two full years before the earthquake (Figs 2–4). Their estimated quasi-coseismic offsets are therefore sensitive to any errors in the linear-motion approximation that we applied to estimate their locations just before the earthquake. Hutton et al. (2001) and Márquez-Azúa et al. (2002) describe evidence for significant transient motion at these and other sites following the 1995 Colima-Jalisco earthquake—although most transient motion appears to have decayed away by mid-1997, the viscoelastic response could require a decade or longer to reach insignificant levels and thus cannot be discounted as the possible underlying cause for our inability to fit the coseismic motions of some inland sites.

Localization of strain along faults in the upper crust either during the Tecomán earthquake or in the years preceding it offers an alternative possible explanation for the east-directed residual misfits at inland sites. For example, the Tamazula Fault described by Núñez-Cornú & Sánchez-Mora (1999) ruptured in a M_w = 5.3 earthquake on 2000 March 7 (Pacheco et al. 2003), thereby indicating that some strain is localized along this fault. We tested whether sympathetic slip triggered by the Tecomán earthquake along a fault coinciding with the postulated location of the Tamazula fault would allow for a better fit to the coseismic offsets for sites in the Jalisco block. To do so, we added to our mesh a slip surface that dips southeast 60° and extends inland from the northwest wall of the Manzanillo Trough, coinciding with the approximate location of the Tamazula fault. We then re-inverted the GPS offsets while allowing for downdip thrusting along the subduction interface and additional dip-slip and strike-slip coseismic motion along the hypothetical crustal fault. The crustal fault shielded the region to its northwest (i.e. the Jalisco block) from significant elastic deformation associated with a subduction thrust earthquake to the southeast. None of the aforementioned data misfits were significantly reduced for this more complex model.

Finally, lateral inhomogeneities in crustal rocks may be responsible for significant regional differences in their elastic responses (e.g. a possible rheological contrast between the volcanic rocks of the Jalisco block and volcanic and sedimentary fill of the Colima Graben). We did not attempt to test whether possible lateral variations in elastic properties (e.g. Poisson's ratio) could substantially improve the fit to our data, primarily because of the underdetermined and non-unique nature of such an effort.

5.2 Comparative analysis of the 1995 and 2003 earthquakes

Given their proximity in time and space, the M_w = 8.0 1995 Colima-Jalisco and 2003 Tecomán earthquakes offer a rare opportunity to contrast earthquakes in nearly identical tectonic settings. To ensure consistent treatment of the 1995 and 2003 geodetic data, including a thrust-only slip requirement for both earthquakes, we re-inverted coseismic offsets reported by Hutton et al. (2001) for the 1995 Colima-Jalisco earthquake using our finite element mesh to create the data kernel G. The resulting best-fitting coseismic slip distribution (Fig. 14a) strongly resembles those of Melbourne et al. (1997) and Hutton et al. (2001) and fits the data equally well (not shown).

The 1995 and 2003 earthquakes ruptured distinctly different parts of the subduction interface (Fig. 14) and overlap by only ∼15 km at the northwest edge of the Manzanillo Trough. Due to the sparser station coverage in 1995, it is unclear whether the overlap suggested by the geodetic data is significant. We suspect however that it is because aftershocks associated with the 1995 and 2003 earthquakes also overlap in the same area (Singh et al. 2003). Except for isolated areas of higher slip, the cumulative downdip slip for both events (Fig. 14c) is consistently 2–3.5 m.

5.3 Triggering and Coulomb stress calculations: the 1932, 1995, and 2003 earthquakes

Much or possibly all of the subduction interface northwest of the Manzanillo Trough ruptured during the 1932 June 3 M_w = 7.9 and 1932 June 18 M_w = 7.8 earthquakes (Singh et al. 1985; Pacheco et al. 1997), after which the subduction interface in this region remained nearly aseismic until the 1995 October 9 M_w = 8.0 Jalisco earthquake, which initiated near the northwest edge of the Manzanillo Trough and propagated northwest for ∼150 km (Courboulex et al. 1997; Escobedo et al. 1998; Hutton et al. 2001). The pair of earthquakes in 1932 share attributes with the 1995 and 2003 earthquakes that suggest a possible triggering relationship for earthquakes along this part of the trench. Singh et al. (1985) approximate the rupture areas of the June 3–18 earthquakes from isoseismal lines defined from 1932 damage reported by regional newspapers (Fig. 15). Like the 1995 and 2003 earthquakes, the 1932 June 3 earthquake ruptured the northwest portion of the Rivera-North America subduction interface, after which the 1932 June 18 earthquake appears to have ruptured the subduction interface to the southeast (Singh et al. 1985).

Figure 15.

Isoseismals from the 1932 June 3 and 1932 June 18 earthquakes (Singh et al. 1985) and rupture areas from the 1973 Colima, 1995 Colima-Jalisco and 2003 Tecomán earthquakes. The isoseismals for the 1932 earthquakes lie approximately onshore from the respective rupture areas of the 1995 and 2003 earthquakes.

Too little is known about the 1932 earthquakes to investigate whether the Coulomb stress changes induced by the 1932 June 3 earthquake significantly unclamped the subduction interface in the vicinity of the eventual 1932 June 18 aftershock. As a proxy, we instead investigated a possible triggering relationship between the 1995 and 2003 earthquakes by using our FEM-derived slip model (Fig. 14) to calculate the Coulomb stress change on the subduction interface induced by the 1995 earthquake. The Coulomb stress change on a fault is defined as ΔCFS = Δτ + μΔσ_n, where τ is the shear stress in the direction of slip and σ_n is the normal stress across the fault (unclamping of the fault is represented by positive Coulomb stress). Following Lin & Stein (2004), we use μ = 0.4 for the coefficient of friction.

Fig. 16 shows ΔCFS along the northern Middle America subduction zone following the 1995 Colima-Jalisco earthquake. Similar to results reported by Lin & Stein (2004) for other large subduction thrust earthquakes, a reduction in the stress component normal to the subduction interface occurs in areas surrounding the coseismic rupture zone. Fault-normal unclamping at the southeast end of the 1995 rupture at the Manzanillo Trough was small, with ΔCFS ≈ 0.1 MPa at the epicentre of the eventual 2003 Tecomán earthquake. Unclamping is also predicted for areas located downdip from the 1995 rupture zone, possibly facilitating afterslip known to occur along those areas of the subduction interface (Hutton et al. 2001; Márquez-Azúa et al. 2002).

Figure 16.

Coulomb stress change on the subduction interface resulting from the 1995 M_w = 8.0 Jalisco earthquake. Slip is encouraged in regions where the fault is unclamped (yellow to red). The epicentre of the 2003 Tecomán earthquake (star) is located in an area of modest unclamping.

Although our Coulomb stress calculations indicate that the 1995 earthquake reduced the fault-normal stress in the eventual rupture zone of the Tecomán earthquake, the change was relatively small and failed to trigger a near-term thrust earthquake, unlike in 1932. By 1997, the combined effects of fault afterslip and viscoelastic rebound had released an equivalent seismic moment that was equal to more than half of the coseismic moment of the 1995 rupture (Hutton et al. 2001), further unclamping the subduction interface beneath the Manzanillo Trough. The additional post-seismic unclamping however also failed to trigger a near-term rupture. By mid-1997, GPS sites COLI, MANZ and CRIP, which are located onshore from the 1995 and 2003 rupture zones, were once again moving northeastward towards the plate interior, as they had prior to the 1995 earthquake (Fig. 1 and Fig. 13) (Márquez-Azúa et al. 2002). This motion indicates that by mid-1997, surface deformation at these sites was dominated by interseismic strain resulting from relocking of the subduction interface. Clamping of the subduction interface was, therefore, increasing in magnitude for more than 5 yr before the 2003 Tecomán earthquake.

Historical evidence clearly indicates that large earthquakes along the northern Middle America trench adjacent to the Manzanillo Trough do not always trigger ruptures of the subduction interface beneath the trough. For example, the 1973 January 30 M_w = 7.5 Colima earthquake, which ruptured a 100-km-long stretch of the Cocos-North America subduction interface immediately southeast of the Manzanillo Trough (Fig. 16), did not trigger a subsequent rupture in the Manzanillo Trough.

5.4 Implications for seismic hazard and the earthquake cycle

That no reliably recorded earthquake has ruptured completely across the Manzanillo Trough and that the Tecomán earthquake was largely limited to areas beneath the Manzanillo Trough may indicate that mechanical discontinuities within either the upper North American plate or subducting Rivera or Cocos plates constitute a barrier to along-strike propagation of large subduction thrust earthquakes in this region. Our Coulomb stress modelling indicates that thrust earthquakes that occur immediately beyond the edges of the Manzanillo Trough modestly unclamp the subduction interface beneath the Manzanillo Trough and thus could trigger a near-term thrust earthquake, as may have occurred in 1932. Whether or not such a triggering relationship applies, the presumed mechanical isolation of the subduction interface beneath the Manzanillo Trough limits the likely rupture extent and hence magnitudes of earthquakes that originate near or beneath the trough and earthquakes that rupture the subduction interface northwest of the trough. Assuming a uniform 3 m rupture of the seismogenic zone beneath the Manzanillo Trough, the maximum implied moment magnitude is M_w = 7.6 for an assumed maximum rupture area of ∼5000 km² (70 km along-strike by 70 km downdip) and shear modulus of 20–30 GPa. Similarly, a hypothetical 4 m rupture of the entire ∼200-km-long seismogenic zone northwest of the Manzanillo Trough to a downdip distance of 80 km (Hutton et al. 2001) would have an equivalent moment magnitude of M_w = 8.0–8.1, representing a likely maximum-size earthquake for the northern end of the Middle America subduction zone.

Acknowledgments

Numerous individuals assisted us in our field efforts. We thank Bill Unger of the University of Wisconsin Department of Geology and Geophysics for the long hours he has dedicated to this project for the past decade. We also thank Bill and Ana Douglass of Ajijic for their hospitality and field assistance, the Universidad de Guadalajara for access to a field vehicle and GPS equipment, Antonio Hernandez of INEGI for access to data from their GPS sites and Proteccion Civil of the state of Jalisco for logistical support. Support from NSF grants EAR-9909377 and EAR-9909321 made this work possible.

Appendix

Appendix A: Inverse Modelling Technique Description

The forward solution for surface displacements associated with slip along a subduction interface (or other surface at depth) is given by the usual system of linear equations, G · m = d. In the present work, G is the data kernel, a matrix dimensioned 3n × 2m, where n is the number of GPS sites and m is the number of node pairs on the subduction interface. Vector m contains 2m elements representing dip-slip and strike-slip displacement at each of m node pairs, and d is the 3n-element data vector containing each of the three components of motion observed at each of n GPS sites.

The data kernel G is constructed by using our 3-D finite element model to predict the components of surface displacement that occur at each of the n sites in response to unit slip at each element in m. Slip on the subduction interface is represented by dislocations between pairs of initially collocated nodes (Masterlark 2003). The dislocation is an antiparallel displacement of the nodes in the pair. Local coordinate transformations are defined in downdip, along-strike and fault-normal directions.

The best distribution of slip is defined by minimizing reduced chi-square χ²_ν, the weighted least-squares misfit χ² divided by the number of degrees of freedom for the model, while accounting for data uncertainties and enforcing both model smoothness and a thrust-only slip constraint. We use a first-order Tikhonov regularization method (Menke 1984; Aster et al. 2004) to estimate the model parameters as follows:

(A1)

The initial estimate m₀ is determined using truncated singular value decomposition and is modified during subsequent iterations by the second term on the right side of eq. (A1). Matrix W_d is a diagonal matrix containing the data weights, which are defined as the inverse square of the data uncertainties. Matrix W_m is a first-order smoothing matrix constructed to penalize large differences in slip between neighboring nodes. Its influence is controlled by the smoothing factor ε.

The thrust-only slip constraint is accomplished using the non-negative least- squares method (Lawson & Hanson 1974). Estimates of m are iteratively refined by following the negative gradient of the weighted misfit prediction error while simultaneously restricting all elements of m to be positive. Slip is further constrained to vanish at nodes that define the edge of the slip surface.

The inversion output also includes the model resolution vector m_r and data importance vector n_r, which are respective diagonals of the model resolution matrix M_r = G^−g · G and data resolution matrix N_r = G · G^−g. The generalized inverse G^−g is defined as

(A2)

Elements of m_r describe how well their corresponding elements in m are resolved, with values that range from 0 (no resolution) to 1 (perfect resolution). Elements of n_r measure the relative amount of information given by their corresponding elements in the data vector d to parameters in the model vector m. Practically speaking, geodetic data that are three-dimensionally closer to the dislocation surface carry higher importance. Similarly, areas on the dislocation surface that are closer to GPS sites are better resolved. The sum of data importance is a measure of the number of model parameters the inversion can resolve.