Near-realtime source analysis of the 20 March 2012 Ometepec-Pinotepa Nacional, Mexico earthquake

Resumen

Se aplica un procedimiento de inversión de falla finita para obtener un modelo de deslizamiento sobre el plano de falla para el sismo M_W 7.4 de Ometepec-Pinotepa Nacional del 20 de marzo 2012, utilizando ondas P telesísmicas registradas por la red del Global Seismographic Network. La inversión se hace en tiempo casi-real utilizando parámetros de la fuente reportados por el Servicio Geológico de Estados Unidos (USGS) y el proyecto del Global Centroid Moment Tensor (gCMT). La orientación de la falla y el ángulo de deslizamiento se obtienen del mecanismo del gCMT, asumiendo que la falla coincide con el plano nodal con bajo ángulo de buzamiento. Las dimensiones de la falla y la duración máxima de la función fuente se definen en base a la magnitud reportada para el evento. Datos telesísmicos obtenidos de la base de datos del Continuous Waveform Buffer del USGS se utilizan en la inversión con tiempos de inicio que corresponden a los arribos de la onda P implementados en el cálculo del hipocentro. La inversión se estabiliza aplicando una transición suave de deslizamiento a lo largo de la falla y una reducción simultánea del momento sísmico. Estas restricciones se implementan utilizando un peso de suavizamiento estimado directamente del problema inverso, permitiendo así la recuperación en un solo paso del patrón de ruptura menos complicado. La inversión de los registros en desplazamiento revela una simple ruptura circular similar al área de deslizamiento determinada por el USGS utilizando ondas de cuerpo y ondas superficiales, indicando que las ondas P telesismicas pueden identificar los rasgos principales de la fuente en tiempo casi-real. Inversiones adicionales realizadas utilizando registros en velocidad identifican una fuente mas detallada que cubre una zona elíptica de 2500 km² y que se extiende hacia arriba y hacia abajo del hipocentro a lo largo de la falla. Esta zona elíptica tiene dimensiones que cubren el área definida por dos fuentes separadas determinadas por otros investigadores utilizando datos locales y ondas sísmicas globales. Los resultados indican que los registros en velocidad podrian aportar mayor información sobre los detalles de la ruptura en inversiones realizadas en tiempo casi-real utilizando ondas P telesísmicas.

Abstract

We apply a single-step, finite-fault analysis procedure to derive a coseismic slip model for the large M_W 7.4 Ometepec-Pinotepa Nacional, Mexico earthquake of 20 March 2012, using teleseismic P waveforms recorded by the Global Seismographic Network. The inversion is conducted in near-realtime using source parameters available from the USGS/NEIC and the Global Centroid Moment Tensor (gCMT) project. The fault orientation and slip angle are obtained from the gCMT mechanism assuming that the fault coincides with the shallow-dipping nodal plane. The fault dimensions and maximum rise time are based on the magnitude reported for the event. Teleseismic data from the USGS/NEIC Continuous Waveform Buffer database are used in the inversion with record start times set to the P-wave arrivals used to compute the earthquake hypocenter. The inversion is stabilized by requiring a smooth transition of slip across the fault while minimizing the seismic moment. These constraints are applied using a smoothing weight that is estimated from the inverse problem, allowing the recovery of the least-complicated rupture history in a single step. Inversion of the deconvolved, ground-displacement waveforms reveals a simple, circular rupture similar in extent to the source identified by the USGS/NEIC using body-and surface-wave data, indicating that the teleseismic P waves can provide a first-order source model for the event in near-realtime. Additional inversions conducted using velocity records identify a more-detailed rupture model characterized by an elliptical 2500 km² source region extending updip and downdip from the hypocenter. This elliptical source preserves the orientation and overall dimensions of a dual-source slip model obtained recently by other investigators using local strong motions and global seismic waveforms. The results indicate that velocity waveforms could provide additional details of the earthquake rupture in near-realtime, finite-fault inversions using teleseismic P waves.

Introduction

The rapid derivation of the extended properties of large earthquakes is an important component of earthquake early-alerting programs and is necessary for effective post-earthquake response. Preliminary earthquake source models are useful, for example, for estimating the ground motions expected in the epicentral region. Several international seismological agencies currently calculate preliminary fault-slip models for large earthquakes following their occurrence. The U.S. Geological Survey National Earthquake Information Center (USGS/NEIC), for example, provides fault-slip information following significant earthquakes using body-and surface-wave data recorded worldwide. Such timely analyses are not currently performed by researchers in Latin America.

The timely recovery of coseismic slip models requires the rapid acquisition and processing of the recorded seismic waveform data followed by the implementation of a waveform-inversion methodology usually based on a simplified parameterization of the earthquake source. Different data types, including regional records and teleseismic body- and surface-wave data, have been examined to rapidly infer the extended properties of large earthquakes using various inversion schemes (e.g., Mendoza, 1996; Ji et al., 2002; Dreger et al., 2005, Ammon et al., 2006; Hayes, 2011; Mendoza et al., 2011). Studies using teleseismic P waves (e.g., Mendoza et al., 2011; Mendoza and Hartzell, 2013) indicate that these data alone can be used to rapidly identify preliminary, first-order rupture models. These studies employ the kinematic, finite-fault inversion scheme of Hartzell and Heaton (1983) to model the earthquake source, with the coseismic slip computed in each of a fixed number of subfaults along a simple planar fault. The Hartzell and Heaton (1983) formulation generally requires running multiple inversions using varying amounts of stabilization to identify the simplest solution that reproduces the observed records. For near-realtime applications, Mendoza (1996) suggested running the inversion only once by fixing the seismic moment to a prescribed value. This approach was found to provide a preliminary image of the coseismic slip pattern for large earthquakes in the Mexico subduction zone, although the models exhibit significant fluctuations in slip due to the lack of spatial smoothing along the fault (Mendoza et al., 2011). More recently, Mendoza and Hartzell (2013) examined the multiple-inversion process of Hartzell and Heaton (1983) and found a relationship between the amplitudes of the observed teleseismic P-wave records and the amount of stabilization required to obtain the simplest rupture model. The relation allows an application of the methodology in a single step, using both moment-minimization and spatial-smoothing constraints estimated from the inverse problem. Mendoza and Hartzell (2013) used this single-step approach to derive first-order slip models for several different-size events, including the 2011 M_W 7.1 East Turkey, the 2011 M_W 9.0 Japan, and the 2012 M_W 8.6 Northern Sumatra earthquakes using P-wave displacement records.

In this study, we use the single-step finite-fault inversion scheme of Mendoza and Hartzell (2013) to analyze the large M_W 7.4 interplate thrust earthquake that occurred in the Guerrero-Oaxaca border region, Mexico, on 20 March 2012 and examine the results to explore the development of an automated inversion process that could be implemented in near-realtime following the occurrence of large subduction earthquakes along the Middle America subdduction zone. The earthquake resulted from shallow subduction of the Cocos plate beneath North America and occurred in a section of the plate boundary that has experienced more than ten M_W > 7 events in the last century (UNAM Seismology Group, 2013). The earthquake caused extensive damage in the towns of Ometepec, Guerrero and Pinotepa Nacional, Oaxaca, and was also strongly felt in Mexico City (Juarez Garcia et al., 2012). Here, we present the results of a near-realtime application of the inversion procedure conducted the day of the event using teleseismic ground-displacement P waveforms, consistent with the data processing used by Mendoza and Hartzell (2013). We also perform additional inversions of the teleseismic dataset using velocity records to examine the extent to which principal features of the rupture history can be derived. We find that the velocity records provide more detail of the coseismic slip compared to the displacement waveforms, due to the higher frequency content. Thus, it would be useful to include these data in near-realtime analyses of the teleseismic P waves to derive provisional, first-order source models for large M_W > 7 interplate subduction earthquakes.

Inversion method

The method is based on the finite-fault waveform inversion scheme developed by Hartzell and Heaton (1983) and modified by Mendoza and Hartzell (2013) to identify the earthquake rupture history from the recorded teleseismic P waveforms in a single step. The method uses multiple consecutive slip intervals of finite duration to identify a long rise time, if required by the observations. In the analysis, a fault plane of prescribed orientation, dimensions, and depth is placed in the crustal structure of the source region and divided into a given number of subfaults. Fault dimensions are chosen to generously encompass the expected extent of rupture based on the moment magnitude of the earthquake. For the 2012 Ometepec-Pinotepa Nacional earthquake, the fault plane corresponds to the gently-dipping nodal plane in the Global Centroid Moment Tensor (gCMT) source mechanism (strike, dip, and rake of 295º/13º/91º), assumed to represent faulting along the Cocos-North America plate boundary. This gCMT mechanism was available from the USGS website the day of the event. The fault length and width were set to 140 km and 70 km, respectively, with the fault divided into 20 subfaults along the strike and 10 subfaults down the dip. The hypocenter was placed at the center of the fault at the 17-km depth computed by the USGS/NEIC.

Synthetic P-wave seismograms are then computed for each subfault at all the recording sites by summing the responses of a set of point sources distributed uniformly across each subfault, taking into account the time delays for a rupture front propagating at a prescribed rupture speed. The point-source responses are obtained using generalized ray theory (Helmberger and Harkrider, 1978) and a near-source crustal structure that corresponds to the upper 35 km of the AK135 velocity model (Kennett et al., 1995). The rupture velocity is fixed at 2.5 km/s, the value suggested by Mendoza and Hartzell (2013) corresponding to about 70% of the average crustal shear-wave speed in the velocity model. A boxcar source-time function is used to calculate the point-source responses, with the boxcar duration depending on the size of the event.

The synthetic records are placed end-to-end to form a matrix A of subfault synthetics. The observed records are similarly placed end-toend for all stations to form the data vector b. Together, these form an overdetermined system of linear equations Ax=b that is solved for the solution vector x whose elements correspond to the subfault slips required to reproduce the observations. Multiple slip intervals are accommodated by adding columns to the A matrix that are constructed by successively lagging the subfault synthetics by the width of the boxcar source-time function, resulting in a specified number of time windows that discretize the rise time on the fault. The inversion thus recovers the slip required of each subfault in each of the specified time windows. For the 2012 Ometepec-Pinotepa Nacional earthquake, we used five time windows of 1s duration, allowing up to 5s for the rise time at any point on the fault. This flexibility allows for the rise-time variations observed by Somerville et al. (1999) for earthquakes of this size.

In the Hartzell and Heaton (1983) formulation, the constraint equations λFx = 0 are added to the linear system to stabilize the inversion, where the smoothing weight λ controls the tradeoff between applying the constraints and fitting the observations. The constraints take two forms, with the linear system to solve given by

where F₁ corresponds to the difference in slip between adjacent subfaults, thus requiring a smooth transition of slip from subfault to subfault. F₂ is the identity matrix, effectively reducing the length of the x vector and the total seismic moment. Cd^-1 is a data covariance matrix that normalizes the data record at each station to its maximum amplitude. The coefficient matrix Cd^-1A thus contains the observed-amplitude information following the data normalization process (Mendoza and Hartzell, 2013).

The linear system corresponds to a discrete ill-posed problem, which can be solved using Tikhonov-type regularization through an iterative L-curve analysis (Hansen, 1998). This iterative process recovers the simplest solution by identifying the optimum smoothing weight λ_s=λ₁=λ₂ that quantifies the balance between fitting the data and meeting the constraints (Mendoza and Hartzell, 2013). An independent estimate of this optimum smoothing weight, however, can be obtained from the average (|a|_avg) of the absolute values of the elements of the coefficient matrix using the relation λ_{s =} 90 |a|_avg (Mendoza and Hartzell, 2013). Thus, the amount of stabilization to apply can be estimated directly from the inverse problem and then used to invert the observed data, allowing the recovery of an earthquake source model in a single step. Note that the seismic moment could be fixed to a prescribed value, instead of minimizing the moment in the F₂ constraint. However, this would require a reevaluation of the relationship between the smoothing weight and the elements of the coefficient matrix to identify the proper smoothing value to use in the inversion.

Near-realtime source analysis

We applied the single-step inversion algorithm of Mendoza and Hartzell (2013) to the teleseismic P waves recorded for the 2012 Ometepec-Pinotepa Nacional earthquake to recover a first-order image of the rupture history the day of the event. The inversion uses the broadband, vertical (BHZ) waveforms recorded by stations of the Global Seismographic Network (GSN) with record start times corresponding to the P arrival times used by the USGS/NEIC to locate the events. The data were obtained from the USGS/NEIC Continuous Waveform Buffer (CWB) repository that are made available to the scientific community through dedicated client servers shortly after they are recorded at the various worldwide stations (ftp://hazards.cr.usgs.gov/CWBQuery/CWBQuery.doc). The data-retrieval procedure uses the phase-arrival times posted online by the USGS/NEIC following the earthquake (http://earthquake.usgs.gov/earthquakes) to request waveforms for stations located 25-95 degrees from the epicenter. Records from the available stations were then deconvolved to ground displacement, interpolated to a 1s time step, passband-filtered between 5 and 100s, and windowed to 50s record lengths. For the 2012 earthquake, a total of 26 teleseismic P-wave records were recovered from the CWB database (Figure 1 and Table 1) for use in the inversion. This data retrieval and waveform processing took less than 5 minutes on a 2-GHz Intel Core 2 Duo iMac desktop computer.

Figure 1.

GSN stations (triangles) used to study the 20 March 2012 Ometepec-Pinotepa Nacional earthquake plotted relative to the USGS/NEIC epicenter (star).

Table 1.

Stations used in teleseismic inversion.

Station	Distance (°)	Azimuth (°)
PASC	25.0	-42.1
BORG	69.5	26.5
ESK	79.3	35.4
FFC	38.1	-3.6
KDAK	57.3	-31.2
NNA	35.4	142.1
COLA	58.5	-22.4
COR	35.0	-32.0
HRV	34.4	36.0
KBS	77.5	10.7
KIP	56.5	-75.1
KONO	84.4	29.1
KOWA	89.8	75.0
OTAV	25.4	128.0
PAB	82.5	51.1
POHA	54.3	-77.4
RCBR	65.5	104.5
RSSD	27.8	-9.0
SAML	43.0	123.8
SDV	27.9	102.6
SFJD	58.9	19.8
SJG	30.6	82.5
SSPA	29.7	32.3
ECSD	27.0	2.5
YMN	31.7	8.6
JFWS	27.0	12.9

The slip model recovered for the earthquake using the deconvolved, teleseismic P-wave displacement records is shown in Figure 2. This result was obtained within 1 minute following the processing of the observed records. The model shows a single, circular source centered at the earthquake hypocenter with a radius of about 30 km and a peak slip of 384 cm. The observed records are fit very well (Figure 3), although there are some misalignments at several stations (e.g., BORG, KDAK, COR) that require small shifts to the USGS/NEIC P arrival times. Runs conducted after adjusting these record start times, however, yield slip models that are very similar to that shown in Figure 2.

Figure 2.

Slip model obtained for the 2012 Ometepec-Pinotepa Nacional earthquake from the inversion of the GSN teleseismic P-wave ground displacements filtered at 5-100s. View is from the top of the fault with distance measured from the hypocenter (star) placed at the 17-km depth obtained by the USGS/NEIC.

Figure 3.

Fits between observed (solid) and predicted (dashed) teleseismic P-wave displacements corresponding to the slip model of Figure 2. Waveforms are normalized to their maximum amplitudes, and numbers to the right indicate the ratio of synthetic-to-observed amplitudes corresponding to a total moment of 1.5 × 10²⁷ dynecm. This total moment, however, includes non-source related contributions from beyond the circular source area shown in Figure 2.

The moment-rate function (Figure 4) shows a major 15s pulse with a seismic moment of 1.3 × 10²⁷ dyne-cm (M_W 7.4) that corresponds to fault slip in the circular region near the hypocenter. There is also a second, smaller contribution to the total moment about 5s later due to slip beyond 50 km from the hypocenter. This second, small pulse has a minimal effect on the predicted teleseismic P waveforms and is attributed to the mapping of unmodeled, propagation effects onto the fault. That is, the observed records are fit just as well when slip from this outer portion of the fault is excluded from the source model. The circular, coseismic source region accounts for about 85 percent of the total calculated seismic moment of 1.5 × 10²⁷ dyne-cm.

Figure 4.

Moment-rate function corresponding to the slip model shown in Figure 2. Moment release in the first 15s is due to coseismic slip within the circular source area of the slip model and corresponds to a seismic moment of 1.3 × 10²⁷ dyne-cm. The second small pulse observed beyond 20s is attributed to the mapping of non-source effects onto the fault.

The slip distribution obtained in near-realtime for the 2012 Ometepec-Pinotepa Nacional earthquake is very similar to the finite-fault model derived by the USGS/NEIC (http://earthquake.usgs.gov/earthquakes/eqinthenews/2012usc0008m6h/) from the analysis of body- and surface-wave data using the method of Ji et al. (2002). The USGS/NEIC model shows a single 30-km x 35-km slip area at the hypocenter with a peak slip of ~4.5m and a seismic moment of 1.5 × 10²⁷ dyne-cm. A recent study of the event by UNAM Seismology Group (2013), however, suggests that two separate sources may have contributed to the seismic-wave radiation. They use teleseismic body and surface waves and local strong-motion data to identify two areas of slip (up to 4m peak) updip and downdip of the hypocenter. This source detail was obtained assuming an arbitrary 2s time shift in the alignment of the observed body-wave records (UNAM Seismology Group, 2013). Independent of whether this 2s shift is necessary, rapid finite-fault analyses such as the one described in this paper rely on initial P arrival measurements and would not be able to discern rupture patterns that depend on the identification of weak precursors on the teleseismic records. Nonetheless, it may be possible to derive more source detail than that shown in Figure 2 by using higher-frequency records that could potentially resolve smaller source dimensions.

Analysis using instrumental records

With this in mind, we perform a second finite-fault analysis of the 2012 Ometepec-Pinotepa Nacional earthquake using the instrumental BHZ waveforms recorded at the 26 GSN teleseismic stations with record start times set to the USGS/NEIC P arrival times. These records correspond to velocity and are given in counts/second in the CWB database. In this case, instrument responses are added to the subfault synthetics rather than deconvolved from the observed records, and both observations and synthetics are interpolated to a common time step of 0.25s prior to inversion. The resulting slip model (Figure 5) shows a more elliptical source with slip extending updip and downdip from the hypocenter.

Figure 5.

Finite-fault model obtained for the 2012 Ometepec-Pinotepa Nacional from the inversion of the broadband, teleseismic GSN instrumental P waveforms. Distances and symbols are the same as in Figure 2.

The direction of elongation is consistent with the N-S slip pattern obtained by UNAM Seismology Group (2013). The peak slip is 413 cm, and the total seismic moment is 1.9 × 10²⁷ dyne-cm. This total moment, however, includes contributions from outside the elliptical source area, more than 45 km from the hypocenter. For a rupture velocity of 2.5 km/s, this outer slip would result from an attempt to fit details of the observed records beyond 18s. Record amplitudes beyond this time, however, are relatively small (Figure 6) and correspond either to unmodeled propagation effects or to minor unresolved features of the source. Some of the record start times are again observed to require time shifts, although runs conducted after making these adjustments do not appreciably modify the inferred rupture pattern. The additional detail provided by the instrumental BHZ records is due to the higher frequency content of the velocity waveforms, compared to the ground-displacement records. That is, similar results are obtained using deconvolved, ground-velocity records in the inversion (e.g., records deconvolved within frequency bands of 0.01667-1 Hz and 0.08-1 Hz).

Figure 6.

Fits between observed, instrumental teleseismic P waveforms (solid) and those predicted (dashed) for the slip model shown in Figure 5. Records are normalized to their maximum amplitudes, and numbers to the right indicate the ratio of synthetic-to-observed amplitudes corresponding to a total calculated moment of 1.9 × 10²⁷ dyne-cm. This total moment, however, includes non-source related contributions due to attempts to fit fine details of the waveforms beyond 18s into the records.

The moment-rate function obtained using the instrumental velocity records (Figure 7) shows a 16s pulse with a seismic moment of 1.4 × 10²⁷ dyne-cm that corresponds to the elliptical source area of Figure 5. This source duration and seismic moment are consistent with those identified for the circular slip region in the near-realtime inversion of the deconvolved displacement records. The lower-amplitude pulse observed beyond 16 seconds in the moment-rate function corresponds to slip mapped onto the outer portions of the fault due to the fitting of unmodeled effects in the observed waveforms.

Figure 7.

Moment-rate function corresponding to the slip model of Figure 5. The moment release in the first 16s is 1.4 × 10²⁷ dyne-cm and corresponds to coseismic slip within the elliptical source region observed in the inferred rupture model. The second pulse observed beyond 16s is attributed to the mapping of non-source effects onto the fault.

Conclusions and discussion

We derived a first-order rupture model for the 2012 Ometepec-Pinotepa Nacional, Mexico earthquake in near-realtime using 26 broadband, GSN teleseismic P-wave records obtained from the USGS/NEIC Continuous Waveform Buffer (CWB) database. The inversion uses the deconvolved, ground-displacement waveforms to identify a relatively simple, circular 2800-km² source near the hypocenter that is comparable in both slip amplitude and rupture extent to the model identified by the USGS/NEIC using body- and surface-wave data. This result indicates that the single-step finite-fault inversion procedure can provide provisional rupture models in near-realtime for large, shallow interplate earthquakes along the Mexico coast using teleseismic P waves. For the Ometepec-Pinotepa Nacional earthquake, the analysis took less than 6 minutes to complete, including the retrieval and processing of the teleseismic P waves and the plotting of the results. For stations located at teleseismic distances, P waveforms would be recorded and stored in the CWB database within 15 minutes following the event, so that a near-realtime inversion could be performed within the first hour provided a source mechanism is available. For large events, moment-tensor mechanisms are generally computed by the USGS/NEIC within this time frame, and these mechanisms could be used in the near-realtime analysis. Alternatively, predetermined fault geometries for subduction events in specific regions could be used to accelerate the process.

We also conducted a second finite-fault analysis of the Ometepec-Pinotepa Nacional earthquake using the instrumental velocity waveforms recorded at the 26 GSN teleseismic stations. The resulting slip pattern shows a more elliptical source that covers an area of about 2500 km² and extends updip and downdip from the hypocenter. Slip is oriented SW-NE, consistent with the detailed, dual- source model obtained by UNAM Seismology Group (2013) using local strong-motion data. The velocity records thus provide additional details of the source compared to the ground-displacement waveforms due to their higher frequency content. It would be useful then to include velocity records in near-realtime applications using teleseismic P waves to recover the earthquake rupture history. These could be either deconvolved ground-velocity records or the instrumental waveforms themselves. An advantage of using the instrumental records is that, since the instrument responses are added to the synthetics rather than deconvolved from the observed records, potentially unstable computational steps are avoided in the preparation of the waveforms for inversion. Some long-period filtering, however, may be required for events greater than M_W 8 to adequately account for the finiteness of the source. For example, in their analysis of the 2011 M_W 9.3 Japan and 2012 M_W 8.6 Sumatra earthquakes, Mendoza and Hartzell (2013) low-pass filtered the teleseismic P waveforms at 200s. The number of time windows and the duration of the source-time function would also have to be greater to allow for longer rise times for these larger events. For near-realtime applications, the rise-time parameters used by Mendoza and Hartzell (2013) to model different-size events could probably be adopted. Mendoza and Hartzell (2013) used these parameters to derive distributions of coseismic fault slip that compare well with models derived using more sophisticated techniques and more complete data sets.

Acknowledgments

H. Benz and D. Ketchum facilitated access to the USGS Continuous Waveform Buffer database. The use of broadband instrumental records was prompted following collaborative work conducted with S. Hartzell. Comments provided by two anonymous reviewers greatly improved the manuscript. Partial support for this work was provided by PAPIIT Project IN104013.