A Seismological Study of the Michoacán-Colima, Mexico, Earthquake of 19
September 2022 (M_w7.6)

Introduction

In the current public perception, 19 September is the date
when large, destructive earthquakes occur in Mexico. The
Michoacan earthquake of 1985 (M_w 8.0), which caused un-
precedented deaths and damage in Mexico City, occurred on
this date. The Puebla-Morelos earthquake of 2017 (M_w 7.1),
which may have been the deadliest intraslab event in the
history of Mexico City, also occurred on the same date. So,
when on 19 September 2022 a subduction thrust earthquake
(M_w 7.6) broke the Cocos-North American plate interface
along the coast of Michoacán-Colima, there was general
consternation and disbelief. The earthquake caused severe
damage to many towns and cities in the states of Michoacán
and Colima (EERI Preliminary Virtual Reconnaissance
Report, 2022). The largest aftershock (M_w 6.7) that occurred
on 22 September caused further damage and panic. Both of
these events were felt strongly in the lake-bed zone of Mexico
City, about 450 km away. The Mexican Seismic Alert System
(SASMEX) performed well; the lead time for the arrival of
strong motion in Mexico City was about 2 minutes (https://
www.youtube.com/watch?v=NCjVeiIZADw).

The tectonic setting of the area of where the 2022 earth-
quake occurred is shown in Figure 1. In the region, the ocean-
ic Rivera (RIVE) and Cocos (COCOS) plates subduct below
Mexico which forms part of the North American (NOAM)
plate. The boundary between the RIVE and COCOS plates,
as well as the relative convergence speed between the two
plates, are controversial. Bandy et al. (1995) suggest that
the subducted RIVE-COCOS boundary lies directly beneath
the southern Colima Rift (SCR) and is parallel to it (Figure
1). The SCR extends from the city of Colima to the Middle

America Trench and forms a part of the Colima rift. CO-
COS-NOAM relative convergence rate at 17.9°N, 104.0°W
is ~ 6.0 cm/yr in the direction 32.3°N (DeMets et al., 2010).

Subduction of RIVE and COCOS plates below NOAM
gives rise to large, shallow thrust earthquakes. Large earth-
quakes that have occurred in the region since 1910 are listed
in Table 1. The aftershock areas of the events, if known, are
shown in Figure 1. The locations of the 2022 mainshock
and its largest M_w 6.7 aftershock are also given in the figure.
We note that the epicenter of the mainshock falls within the
aftershock area of the 1973 earthquake (M_w 7.6) outlined
by Reyes et al. (1979) based on seismograms recorded on a
portable network deployed in the field.

The three largest subduction thrust earthquakes in Mex-
ico since 1900 have occurred along the Michoacán-Coli-
ma-Jalisco segment of the Mexican subduction zone. The
earthquakes of 3 June 1932 (M_s 8.2) and 9 October 1995 (M_w
8.0) ruptured the RIVE-NOAM plate interface, whereas the
19 September 1985 (M_w 8.0) event broke the COCOS-NOAM
interface. The earthquakes listed in Table 1 caused damage
to towns and cities in the vicinity of their rupture areas but
two of them were also destructive to Mexico City. The earth-
quake of 7 June 1911 (M_s 7.7) destroyed the town of Ciudad
Guzmán in the state of Jalisco. It also caused considerable
damage in Mexico City (Miranda y Marron, 1911-1912). As
mentioned earlier, the 1985 Michoacán earthquake caused
unprecedented damage and deaths in Mexico City.

In this paper, we present a source study of the 2022
earthquake and its major M_w 6.7 aftershock in the context of
previous large earthquakes in the vicinity, and discuss the
characteristics of the ground motion at regional distances.

107°       -106°       -105°       -104°       -103°       -102°

Table 1. Large subduction thrust earthquakes since 1910 in the region of interest

References and notes keyed to event number in Table 1

• 1.   Location from ISC-GEM catalog; Ms from Abe (1981)

• 2.   Location from ISC-GEM catalog; aftershock area from Singh et al. (1985); M_s from Abe (1981); M_w from Wang et al.      (1982)

• 3.   Location from ISC-GEM catalog; aftershock area from Singh et al. (1985); M_s from Abe (1981); M_w from Wang et al.      (1982)

• 4.   Location from ISC-GEM catalog; M_s from Abe (1981). Location given by Kelleher et al. (1973) is: 18.85⁰N, 102.94⁰W

• 5.   Location from ISC-GEM catalog; aftershock area, Ms and M_w from Reyes et al. (1979)

• 6.   Location and aftershock area from Havskov et al. (1983); M_s and M_w from Global CMT catalog

• 7.  Location and aftershock area from UNAM Seismology Group (1986); M_s and M_w from Global CMT catalog

• 8.  Location from ISC-GEM catalog; Ms and Mw from Global CMT catalog

• 9.   Location and aftershock area from Pacheco et al. (1997); M_s and M_w from Global CMT catalog

• 10.  Location and aftershock area from Singh et al. (2003); M_s and M_w from Global CMT catalog

• 11.  Location and M_w from this study; M_s from Global CMT catalog

Our analysis is based on local and regional data as well as
teleseismic P-wave data. We also discuss the probability of
having observed three major earthquakes on the same day.

Epicentral recording

Maruata station (MMIG), located on the coast of Michoa-
cán and nearly above the hypocenter [(S-P) time 2.8 s], is

Figure 2. Acceleration, velocity, and displacement at the epicen- tral station of Maruata (MMIG) during the mainshock. Circles in the bottom frame show coseismic static displacement retrieved from GPS station TNMR collocated with MMIG (Z: +25.3 cm; EW: -3.8 cm; NS: -34.4 cm)

equipped with a broadband seismograph, an accelerograph,
and a GPS receiver. The broadband seismograms were sat-
urated on the S-wave arrival. The acceleration traces were
integrated to obtain velocity and displacement. Because
of the baseline shift in the acceleration, the velocity often
does not approach the expected zero level at the end of the
recording; instead records often show a residual velocity.
Integration of these velocity recordings without a shift
correction leads to unrealistic displacements. To correct the
shift, we selected a time, T_h after the end of the intense part
of motion and fit, in the least-square sense, a straight line to
the velocity data between T1 and the end of the record. The
line at T is then connected to time T₀ which we choose at the
P-wave arrival. These two-line segments are used to correct
the velocity record, which are then integrated to obtain the
displacement (see, Singh et al., 2020 for more details). We
followed this procedure in the integration. The traces are
shown in Figure 2. The PGA and PGV on the NS component
are 1090 gal and 28.3 cm/s, respectively.

The GPS receiver at MMIG had stopped working 20 days
before the mainshock due to a problem with the solar panel.
The station was reestablished 4 days after the event. Succes-
sive measurements show post-seismic creep. Correcting for
the lost time series before the 2022 earthquake by extrapola-
tion and for the post-seismic creep, the estimated coseismic
static NS, EW, and vertical, Z displacements from GPS are
-34.4 cm, -3.8 cm, and +25.3 cm, respectively. These values
are marked in the bottom frame of Figure 2, which shows
the displacement seismograms. Not surprisingly, the static
displacement from GPS differs from that estimated from
integration. The MMIG traces are reminiscent of the epi-
central recording at Caleta de Campo (CALE) during the 19
September 1985, Michoacán earthquake (M_w 8.0) (Anderson
et al., 1986) with some differences: PGA at CALE during
the 1985 earthquake was much smaller (141 cm/s²; NS and
EW), PGV was about the same (24.7 cm/s; NS), and PGD
was greater (78 cm; NS).

Basic source parameters of the mainshock and the
major aftershock

Since 2014, the Servicio Sismológico Nacional (SSN, Mex-
ican National Seismological Service) routinely calculates
and publishes M_w through W-phase inversion (Kanamori
and Rivera, 2008) using an algorithm modified by Hayes
et al. (2011) and revised by Duputel et al. (2012). For M >
5.2 earthquakes, the algorithm automatically gets triggered
10 minutes after the origin time and uses broadband data of
the SSN stations (Pérez-Campos et al., 2019). It starts with
the preliminary, automatically obtained, SSN location and
magnitude, and looks for the best half duration and then the
best location. For the 19 September 2022 mainshock and its
major aftershock of 22 September we revised the routine
near-realtime W-phase solution by checking and, if required,
updating the response files and eliminating data with obvious
problems. The revised solutions of the mainshock and the
major aftershock are listed in Tables 2 and 3, respectively.
The tables also give the source parameters reported by the
United States Geological Survey (USGS) and the Global
Centroid Moment Tensor (GCMT) project.

There are some differences in the focal mechanism and
seismic moment (M₀) given by the three sources. For ex-
ample, M0 of the mainshock estimated in this study and by
the USGS are nearly the same, 2.7x10²⁰ N-m (M_w 7.55) but
the value listed in the GCMT catalog is 1.7 times greater.
Henceforth, we shall take M₀ of the mainshock and the ma-
jor aftershock as 2.7x10²⁰ N-m (M_w 7.6) and 1.6x10¹⁹ N-m
(M_w 6.7), respectively. We note that, with respect to the SSN
epicenter, the USGS epicenter is shifted by 44 km towards
N53°E for the mainshock and 29 km towards N36°E for the
aftershock. A consistent NE shift of the epicenters of Mex-

Table 2. Source parameters of the 19 September 2022, Michoacán-Colima earthquake

Table 3. Source parameters of the major aftershock of 22 September 2022

* Depth fixed.

+ A grid search was performed for the depth and the centroid location.

^x Global CMT and USGS/NEIC source parameters last accessed on 06/12/2022.

ican subduction zone earthquakes reported by international
agencies has been documented earlier (Singh and Lermo,
1985; Hjorleifsdóttir et al., 2016).

Aftershock distribution

Aftershocks that occurred in the first 30 days (805 events with
coda-wave magnitude Mc > 3.5) are shown in Figure 3. We
determined CMT solutions of seven significant aftershocks
in addition to the major aftershock (Table 4). Focal mech-
anisms of the mainshock and the eight aftershocks (thrust:
five; normal: two; strike slip: one) are displayed in Figure 3.

Several features of the aftershocks are worth noting in
Figure 3. They overlap the elliptical 1973 aftershock area
outlined by Reyes et al. (1979). Relatively few aftershocks
occurred within the large coseismic slip area of the 2022
earthquake (see next section). Relative lack of aftershocks
over the areas of large slip has been reported for many
earthquakes (see Das and Henry, 2003 for a review). Most

Table 4. Source parameters of seven additional, significant aftershocks

of the aftershocks of the 2022 earthquakes were concentrated
to the south of the mainshock epicenter and in the SE part
of the 1973 aftershock area. A similar concentration was
observed in the aftershock distribution in 1973 which led
Reyes et al. (1979) to suggest that the rupture initiated to
the SE and propagated to the NW.

Finite fault model of the mainshock and the major
aftershock

We determined slip models for the earthquakes of 19 and 22
September 2022 using the rapid finite-fault inversion meth-
odology described by Mendoza and Martínez-Lopez (2022).

The method automatically assigns fault parameters based
on the earthquake size and derives a coseismic slip model
using teleseismic P waveforms obtained in near-realtime
from the Incorporated Research Institutes for Seismology
Data Management Center (https://ds.iris.edu/).

For the 19 September M_w 7.6 earthquake, we used the
hypocenter and moment-tensor source mechanism reported
by the USGS following the event (Table 2; earthquake.usgs.
gov/earthquakes/search/). The slip model for the shallow,
northeast-dipping plane shows two separate zones of high
slip: one downdip of the hypocenter with a peak slip of 1 m
and a second zone about 40 km to the northwest with a maxi-

mum slip of 1.3 m (Figure 4). This result was obtained within
three hours of the occurrence of the event. The rapid P-wave
inversion methodology was also applied following the M_w
6.7 aftershock of 22 September. We used fault dimensions
of 80 km by 80 km, the minimum size allowed in the rapid
P-waveinversion procedure designed to analyze earthquakes
of magnitude M_w 7 or greater (Mendoza and Martinez-Lo-
pez, 2022). For this event, we used the epicenter calculated
by the SSN (Table 3; http://www2.ssn.unam.mx:8080/
sismos-fuertes/) and the focal depth obtained by the USGS.
The distribution of coseismic slip for the shallow-dipping
plane (Figure 5) shows a single 20 km by 20 km rupture area
with a peak of 1.1 m extending primarily downdip from the
hypocenter. Although the results obtained for both events are
preliminary, they provide a general overview of the locations
of high slip and the possible direction of coseismic rupture.

Both inversions use five 1 s time windows to parameterize
the slip duration on the fault.

On 7 October 2022, USGS updated its previously pub-
lished finite fault model of the mainshock (http://earth-
quake.usgs.gov/earthquakes/eventpage/us7000i9bw/
finite-fault). The new fault model is based on the analysis
of a more extensive dataset: 41 teleseismic P waves, 23 te-
leseismic SH waves and 55 long-period surface waves, and
observations from 7 high-rate GNSS stations and 11 static
GNSS sites. The model also uses the hypocenter reported
by the SSN (Table 2) to correct for the location bias. Figure
3 reproduces the USGS finite fault model. As this model
is based on a more extensive dataset, we shall use it in our
further analysis. In this model, M₀ and maximum slip (D_max)
are 2.73 • 10²⁰ Nm and 3.2 m, respectively. Following Ye et
al. (2016) and Lay et al. (2016), we ignore subfaults with
slip D < 0.15D_max as the low slip areas are likely to be poorly
resolved. The trimmed area, A, enclosing D > 0.15D_max, is
3600 km². M₀ released over this area is 2.52 • 10²⁰ Nm and the
average slip, <D>, is 1.48 m. The relation Ao_s = (7n^3/2/16)
(M₀ /A^3/2), where Ao_s is the static stress drop (Kanamori and
Anderson, 1975), yields An_s of 3.7 MPa.

Moment-scaled radiated seismic energy, REEF, and
number of aftershocks

Radiated seismic energy, E_R, for the mainshock, from te-
leseismic data, is 3.44 ± 0.13 • 10¹⁵ J (Me = 7.46). In the
estimation of E_R, we followed the methodology of Boatwright
and Choy (1986), and included a stronger attenuation correc-
tion for subduction earthquakes discussed by Pérez-Campos

and Beroza (2001) and Pérez-Campos et al. (2003). Follow-
ing Boore and Joyner (1997) we applied a correction for
generic hard site. E_R estimation shows a strong azimuthal
dependence that can also be appreciated from the moment
rate spectrum (MRS) at each station (Figure 6). The larger
values are obtained at stations to the north, while the smaller
once occur to the south. We build the source spectrum by
patching, at low frequencies (< 0.2 Hz), the moment rate
function obtained from the source time function reported
by the USGS (http://earthquake.usgs.gov/earthquakes/
eventpage/us7000i9bw/finite-fault), and, at high frequencies
(> 0.2 Hz), the source spectrum obtained from teleseismic
data. The resulting MRS fits the theoretical spectrum from
the Brune source model (Brune, 1970) with a stress drop
of 1 MPa (Figure 6c). The moment-scaled radiated energy,
E_r / M₀, is 1.27 • 10^-5, a value similar to those reported for
other large Mexican thrust earthquakes, which range between
1.0 and 3.3 • 10^-5 with the exception of earthquakes whose
rupture areas extend up to the trench, e.g., Colima-Jalisco
earthquake of 9 October 1995 (E_R / M₀ = 5.6 • 10^-6) (Table
5). E_r / M₀ of the 2022 mainshock is close to the world-wide
average of ~ 1 • 10^-5 (e.g., Ye et al., 2016a).

For the major M_w 6.7 aftershock of 22 September, E_R is
1.33 ± 0. 06 • 10¹⁴ J (M_e = 6.51) so that Er / M_o = 8.31 • 10⁶;
in this case, E_R and source spectrum at individual station
do not show any azimuthal dependence (Figure 6b). MRS
of the earthquake is well fit by a Brune source model with
a stress drop of 0.5 MPa (Figure 6d).

We computed radiated energy enhancement factor,
REEF, for the 2022 mainshock. REEF, a measure of rupture
complexity, recently introduced by Ye et al. (2018). It is
the ratio of measured radiated energy, ER, to the calculated
minimum energy for a source of the same M0 and duration,
E_r/ E_R-min. A smaller REEF value corresponds to a simpler
source and vice versa. The duration, T, of the moment rate
function (MRF) of the 2022 earthquake from the USGS
finite-fault modeling is 32 s. E_R-min, corresponding to M₀ =
2.73 • 10²⁰ Nm and T of 32 s, is 4.1 • 10¹⁴ J (Equation 1 of Ye
et al., 2018), which gives a relatively low REEF value of
8.5. REEF values are consistently low for southern Mexico
to Middle America subduction thrust earthquakes (Table
5; Ye et al., 2018), reflecting the simplicity of the MRF of
the earthquakes along this segment of the subduction zone.

The relatively small number of m_b > 5 aftershocks is also
a characteristic of large Mexican subduction thrust earth-
quakes (Singh and Suárez, 1988). For the 2022 earthquake
there were four aftershocks with m_b > 5 in 30-day period.
Table 5 gives the number of aftershocks, N, in a 30-day pe-
riod with m_b > 5 and log (N / Ne), where Ne is the expected
number of aftershocks derived from regression analysis of

Table 5. Moment-scaled radiated seismic energy, REEF, and number of m_b > 5 aftershocks in one-month period of large Mexican subduc-
tion thrust earthquakes (Modified from Iglesias et al., 2022)

*REEF: Radiated energy enhancement factor (Ye et al., 2018)
N count includes mainshock as one event

+logNe = M_w - 6.34 (Singh and Suárez, 1988)

world-wide data: log Ne = M_w - 6.34 (Singh and Suárez,
1988). Log (N/Ne) is negative for six earthquakes including
the 2022 earthquake and close to zero for the remaining
three. Thus, along the Mexican subduction zone both low
REEF and relative lack of aftershocks prevail. Similarly to
Iglesias et al. (2022), we envision a plate interface that is
relatively smooth, containing discrete, compact asperities.
Asperities rupture smoothly, generating relatively simple
moment rate functions and low values of REEF. As the
rupture area and adjacent plate interface is also smooth
and homogeneous, there is a relative lack of aftershocks
at m_b > 5 level.

Comparison with earthquakes of 30 January 1973
(M_s 7.5, Mw 7.6) and 15 April 1941 (Ms 7.7)

From the aftershock locations, and the relative locations of
the main shock and aftershocks, Reyes et al. (1979) suggested
that the rupture during the 1973 earthquake began to the SE,
near the region of high aftershock activity, and propagated
to the NW. For the 2022 earthquake, the unilateral rupture
propagation to the NW is, of course, well established. In as
much as the aftershock areas of the 2022 and 1973 earth-
quakes overlap (Figure 3), and their magnitudes are similar
(Table 1), it is possible that the two events broke roughly the
same area, had similar gross source characteristics perhaps
even with similar source directivity.

We note, however, that the finite fault model of the
1973 earthquake constructed by Santoyo et al. (2006) using
teleseismic P waves does not show a NW directivity. This
may be due to poorer quality and limited quantity of data
(8 stations) used in the inversion for the 1973 earthquake.
Even with far more data of better quality (20 stations) for
the 2022 earthquake, the inversion of teleseismic P waves
yields a solution that is only a rough approximation of the
one obtained by the USGS finite fault modelling based on a
more extensive dataset (compare Figures 4 and 3).

For the 1973 earthquake, Reyes et al. (1979) noted that
M0 increased by a factor of about 2 as the period increased
from 100 to 300 s. They attributed this increase to possible
slow slip before or after the main slip or to unknown errors
in the estimation of M₀ at lower periods. For the 2022 earth-
quake, we computed M₀ from W-phase CMT inversion of the
regional broadband seismograms with different band-pass
filters and found negligible change in M₀ (M_w) with period
(Table 6). Thus, either the source processes of the 1973 and
2022 earthquakes differed or else the dependence of M₀ on
period for the 1973 earthquake was due to unknown errors.

Much less is known about the 1941 earthquake. Kelleher
et al. (1973) relocated the mainshock and two of its after-
shocks. This area roughly coincides with the aftershock areas
of the 1973 and 2022 earthquakes.

To test whether the 1941, 1973, and 2022 earthquakes
ruptured roughly the same area, we compared their Galitzin
seismograms (Z component) at DBN. We note that repeating
events have the same rupture area and slip and give rise to
identical seismograms.

The 1941 and 1973 analog records were vectorized and
the time series was sampled at an evenly time interval using
TIITBA-GUI (Corona-Fernández and Santoyo, 2022). Gal-
itzin record of the 2022 earthquake was synthesized from
broadband DBN seismogram as the operation of the Galitzin
seismograph was discontinued in December 1994 (Dost
and Haak, 2006). We first note that the three events have
complex P waves that bear some resemblance (Figure 7).

In Figure 8, the seismograms of 2022 and 1973 are com-
pared over three different time windows. The waveforms
are clearly not identical. Our tests, however, show that the
surface waves on the Galitzin seismograms at DBN of events
along the Mexican subduction zone which are 20 to 30 km
apart greatly differ from each other (Singh et al., 2022). In
as much as the character of the surface waves from the 1973
and 2022 earthquakes are similar (bottom frame, Figure 8),
we surmise that the rupture areas of the two events were
less than 30 km apart. From the similarity of the aftershock
areas, the waveforms at DBN, and the magnitudes of the
2022 and 1973 earthquakes, we conclude that they were
quasi-repeated events. In other words, these two events may
have ruptured roughly the same area. If so, the return period
was 50 years. We recall that the finite fault modeling yields
an average slip of about 1.48 m for the 2022 earthquake.
As the plate convergence rate is 6.0 cm/yr (DeMets et al.,
2010), this gives a coupling ratio of 0.49.

The Galitzin seismograms of 2022 and 1941, shown in
Figure 9, exhibit little resemblance. The difference is marked
in the character of the surface waves (bottom frame). Since
the waveform of the 1941 earthquake differs significantly
from those of the 2022 and 1973 events, it most likely did
not rupture the same area as the other two.

Directivity and azimuthal dependence of ground
motion

As discussed above, a source directivity towards NW during
the mainshock is clearly seen in the results of inversion of slip
on the fault as well as in plots of MRS and E_R as a function of
azimuth. A downdip directivity is also visible, albeit weakly,
in the slip inversion of the M_w 6.7 aftershock. In this section,
we examine, in detail, the effect of the source directivity on
the ground motion at regional distances. The stations whose
recordings are used in the analysis are shown in Figure 10.

• (a) Visual examination of the recordings

Figure 11a compares mainshock waveforms at stations CJIG
(azimuth 0 = 308⁰) and ZIIG (0 = 109⁰). The stations are
located at nearly the same distance but in opposite directions
(Figure 10). The shorter duration and higher amplitude of
the intense part of the motion at CJIG compared with ZIIG
strongly suggests a rupture propagation towards the NW. The
waveforms during the aftershock at the same two stations

are shown in Figure 11b. The accelerations, in this case,
are higher at ZIIG (which may be due to site effect) than
at CJIG, while velocities and displacements are about the
same. These waveforms do not support along strike direc-
tivity during the major aftershock; rupture propagation to
the east is certainly viable.

• (b) Spectral ratios of the mainshock to the aftershock
  ground motions

Under the assumption that the mainshock and the after-
shock are collocated and have similar focal mechanisms,
the spectral ratio of the ground motion at a given station
provides the ratio of their moment rate spectrum, MRS. In
the absence of directivity, the MRS is expected to be inde-
pendent of azimuth. Figures 12a,b,c,d illustrate the spectral
ratios at selected stations, each frame comprising stations in
a range of azimuth with respect to the mainshock directivity.
Frame (a): rupture propagating towards the stations; frame

• (b) : station perpendicular to the rupture propagation; frame

• (c) and (d): rupture propagating away from the station. The
  spectral ratios were computed for each of the three compo-
  nents of the ground motion. The figures also show the geo-
  metric mean of the ratios in each frame. For reference, the
  theoretical spectral ratio corresponding to Brune or source
  model (Brune, 1970) with constant stress drop of 3 MPa is
  included in each frame.

A strong dependence of the ratios on azimuth is immedi-
ately obvious. With respect to the theoretical spectral ratio,
the observed ratios are higher in frame (a), about the same
in frame (b), but lower in both frames (c) and (d). Directivity
towards NW during the mainshock and an absence of ESE
directivity during the aftershock are consistent with the
observations.

(c) PGA and PGV ratios of mainshock to aftershock
The directivity effect should also be reflected in the azimuthal
dependence of PGA and PGV ratios of the mainshock to the
aftershock. Horizontal and vertical PGA ratios, plotted in
Figure 13 a, are a strong function of station azimuth p but not
of distance R. Here, horizontal PGA = [(An² + Ae²) / 2]^1/2,

_i_i_i_i_i_i__i_i_i_i_i_i_-0.5-¹-¹-¹-¹-¹-¹-■—-0.5-¹-¹-¹-¹-¹-*—
0      40     80     120         0      40     80     120               ^{0      40     80     120         0      40     80     120}

Time, s                     Time, s                          Time, s                     Time, s

where A_N and A_E are maximum accelerations on NS and EW
components. The ratios rapidly decrease from about 12 to 3
at stations in the azimuthal range 300° <   < 360°. These

stations are in the forward direction for the mainshock and,
possibly, in the backward direction for the aftershock. The
ratio slowly decreases from about 3 to 1 in the range 0° <
< 115°. Stations in this azimuthal range are in the backward
direction for the mainshock and, probably, in the forward
direction for the aftershock. Again, the ratios in the figure
are in agreement with the directivity of the two earthquakes.
The effect of the directivity on the PGA ratios is better ap-
preciated by comparing them with the horizontal PGA ratio
of 2.5 expected from the ground motion prediction equation
(GMPE) for Mexican subduction thrust earthquakes of M_w 7.6
and M_w 6.7 (Arroyo et al. 2010).

PGV ratios, shown in Figure 13b, follow the same trend
as the PGA ratios. However, the maximum PGV ratios in
the azimuthal range 300° <  < 360° exceed 20.

• (d) PGA and response spectral amplitudes as function
  of azimuth and distance

PGA and Sa (T=2 s) for the mainshock and the aftershock are
plotted in Figure 14 as a function of the closest distance from
the fault surface, Rrup. Only stations with Rrup < 600km are
included in the figure. The stations are grouped in 3 bins as
a function of their azimuth: bin 1: 330⁰ < y < 30⁰; bin 2: 30⁰
< < 90⁰; bin 3 : 90⁰ < < 120⁰. All data, except one, are
contributed by stations at Rrup > 120 km. Superimposed on
the data are the predicted curves from the GMPE of Arroyo
et al. (2010). We note that: (i) In general, PGA values are
above the predicted curves for both events irrespective of
the bin. (ii) Sa (T=2 s) values for the aftershock in all bins
are greater than the predicted curve. The values are smaller
than predicted in bin 3 for the mainshock, consistent with
its NW source directivity.

Ground motion in the Valley of Mexico

Since there was a difference of 0.9 in the magnitude of the
mainshock and the major, M_w 6.7 aftershock, it was surprising
that they were felt with nearly equal intensity in the Valley
of Mexico. At CU, a hill-zone reference site in Mexico City,
the PGA on the NS, EW, and Z components during the
mainshock and the aftershock were (5.5, 4.5, 2.9 gal) and
(6.3, 4.2, 2.5 gal), respectively. Was the source directivity
the cause of the similarity of the PGA s?

A site-specific GMPE for CU from subduction thrust
earthquakes has been recently developed by Arroyo et
al. (2022). Figure 15 compares the observed Sa with the
predicted ones for M_w 7.6 and M_w 6.7 earthquakes. As ex-
pected, the observed Sa curves are similar. Predicted Sa for
an M_w 7.6 earthquake, on the other hand, is significantly
higher than the observed one. The converse is true for the

M_w 6.7 aftershock: the observed Sa is much greater than
expected. The source directivity away from CU during
the mainshock explains the smaller Sa. Greater Sa during
the aftershock may be attributed to rupture towards CU.

To further appreciate the role played by directivity,
we used the CU recording of the M_w 6.7 aftershock as
an empirical Green's function (EGF) and synthesized
ground motion from a target M_w 7.6 earthquake. A meth-
od developed by Ordaz et al. (1995) was followed in the
synthesis. The stress drop, Ao, was assumed to be the
same for both events and taken as 3 MPa. The median of
Sa simulations for the postulated M_w 7.6 event as well as
the observed Sa during the mainshock are shown in Figure
16. We find that an M_w 7.6 earthquake with directivity
similar to the aftershock would have produced Sa at CU
about 2.5 times greater than the observed one.

Probability of having observed three major earth-
quakes on the same day

Let us estimate the probability of having observed what we
observed: three major earthquakes (M > 7) that occur on
exactly the same day in the last 120 years. Let's start with
basic data and an assumption:

(1) In central Mexico, an average of 0.46 earthquakes
occur with M > 7 per year; that is, on average one
every 2 years, more or less.

(2) We assume that, in time, earthquakes occur as a
Poisson process; this is relevant in order to know the
probability distribution of the number of events that
we would observe in any given year

It is not difficult to calculate the probability that, in the
span of 120 years, we would have observed 2 or 3 earth-
quakes occurring on the same day. On September 19, let's say,
but the probability would be the same if we chose another
date. Although the calculation is not difficult, it is easier to
calculate the probabilities by simulation.

Using a sample of 10 million possible realizations from
those 120 years, we obtained that the probability of hav-
ing observed 2 events on September 19 is 0.0103 and the
probability of having observed 3 is 0.0005154.They seem
improbable events. But, in reality, this is not the probability
that interests us. Deep down, it strikes us that we have had 3
major earthquakes on the same date, not specifically on 19
September. Indeed, if we had observed 3 major earthquakes
on, say, 24 May, we would be just as surprised.

So, the probability we are interested in is the probabil-
ity of having observed 3 major earthquakes on the same
date, not necessarily September 19. This is more difficult
to calculate with combinatorial analysis (there is a closed
formula, but complicated to apply), although just as easy to
calculate with simulations. We obtain that the probability of
having observed 2 earthquakes on the same date, whatever
it may be, is 0.98 and the probability of having observed 3
earthquakes is 0.18.

This is amazing. According to this analysis, it is extremely
likely to have two large earthquakes on the same date if we
look at 120 years at a rate of 0.46 earthquakes/year. And on
the other hand, the probability of observing 3 is low, but
not astronomically low. If the probabilities are that big, both
events should have already happened. Well, yes: before 19
September 2022 there were already 6 pairs of large events
that had occurred on the same dates and another triad of
events that occurred on 7 June (7 June 1911, M_s 7.7; 7 June
1982, 06:52, Ms 6.9; 7 June 1982, 10.59, Ms 7.0). We just
didn't remember.

Why choose an observation period of 120 years? We
chose 120 years because, on the one hand, it is the period

(1900-2022) in which we consider that the earthquake cat-
alog is complete for M > 7. On the other hand, we chose
it because it seems that we would be equally surprised if
the first event of the sequence of 3 on the same date had
occurred in 1910, let's say, and not in 1985; but maybe we
wouldn't be so surprised, so we calculated the probabilities
for other lapses (Table 7). We confirm that, in reality, what
we observed was not so improbable.

Discussion and conclusion

There is evidence suggesting that the 2022 earthquake (Ms
7.6, M_w 7.6) is a quasi-repeat of the 1973 event (M_s 7.5,
M_w 7.6): their aftershock areas approximately coincide, the
Galitzin seismograms of the two events at DBN are reason-
ably similar, and the magnitudes are the same. Curiously, the
aftershocks of both earthquakes were also concentrated at
the SE end of the rupture area. This distribution of the 1973
aftershocks led Reyes et al. (1979) to suggest that the rupture
during 1973 propagated towards the NW. This directivity is
certainly true for the 2022 earthquake. However, finite fault
modelling of the 1973 earthquake by Santoyo et al. (2006),
using teleseismic P waves recorded at 8 stations, does not
show the NW directivity. Also, an increase in the seismic
moment by a factor of about 2 for the 1973 earthquake as
the period increased from 100 to 300 s, noted by Reyes et al.
(1979), is entirely absent from the 2022 earthquake (Table
6). These differences could be a consequence of increase
in the quality and quantity of data and improvement in the
analysis technique since 1973. It is also possible that the
differences are real and the details of the rupture process
of the two events differed even if their source areas were
roughly the same.

Reyes et al. (1979) suggested that the 1973 earthquake
may have been a repeat of the 1941 event (M_s 7.8). Galitzin
seismogram of the 1941 earthquake at DBN, however, bears
little resemblance with those of the 1973 and 2022 events
(Figures 8 and 9) which suggests that the source region of
the 1941 earthquake was different from those of the other
two events.

A unilateral rupture propagation, along the strike towards
the NW, during the 2022 mainshock is a robust feature of
the finite-fault models. Azimuthal variation of moment rate
spectrum and radiated seismic energy estimated from tele-
seismic P waves also support the NW directivity. According
to the USGS finite fault model, the rupture area over which
the slip is greater than 15% of its maximum value (3.2 m) is
3600 km² (90 km x 40 km). The average slip over this area is
1.48 m, which yields a static stress drop of 3.8 MPa.

If we accept that the 1973 and 2022 earthquakes ruptured
the same area, then the recurrence period is 50 years. For a
plate convergence rate of 6.0 cm/yr and perfect coupling, the
accumulated slip deficit in 50 years would have been 3.0 m.

If we ignore post-seismic slip, then the estimated coupling
ratio on this segment of the plate boundary, corresponding
to a coseismic slip of 1.48 m, is 0.49. This estimate agrees
surprisingly well with the GPS-derived coupling ratio for
this segment (Cosenza-Muralles et al., 2022a) and is slightly
smaller than the coupling of about 0.6 estimated from InSAR
and GNSS data (Maubant et al., 2022). The post-seismic slip,
however, may not be negligible. Similar or larger seismic
moments than the coseismic moments were released in the
post-seismic slip following the earthquakes of 2003 Tecomán
(M_w 7.5) and 1995 Colima-Jalisco (M_w 8.0) (Cosenza-Mu-
ralles et al., 2022b). The areas of post-seismic slip of these
two earthquakes partly overlap their rupture areas and, in both
cases, extend further downdip. These earthquakes, however,
occurred on the RIVE-COCOS plate boundary (Figure 1).
Characteristics of post-seismic slip on the COCOS-NOAM
plate interface, where the 2022 earthquake occurred (Figures
1 and 3), might be different, an issue that future studies will,
no doubt, address.

Our analysis of ground motions at regional distances
confirm the mainshock directivity to the NW. In our study,
we focused on the ratio of ground motions during the main-
shock and the major M_w 6.7 aftershock, thus minimizing the
site effect. These results can be interpreted by a strong NW
directivity during the mainshock and an ENE or negligible
directivity during the M_w 6.7 aftershock. Because of the
directivity, the ground motions in the Valley of Mexico
during the 2022 mainshock and the M_w 6.7 aftershock were
about the same in spite of 0.9 difference in their magnitudes.

It is well known that the source directivity has a profound
effect on the azimuthal variation of ground motion and,
hence, in the damage distribution (e.g., Somerville et al.,
1997; Koketsu et al., 2016). Directivity has been reported
even during small earthquakes (Boatwright, 2007; Calderoni
et al., 2013; Seo et al., 2022). Strong directivity was reported
during two moderate earthquakes in the Guerrero seismic
gap (8 May 2014, M_w 6.5; 11 May 2014, M_w 6.1) (Singh et
al., 2019). The recent Acapulco earthquake of 8 September
2021 (M_w 7.0) had a strong NE directivity (Iglesias et al.,
2022). Directivity, almost certainly, played a major role in
causing damage to Mexico City during the1985, Michoacán
earthquake (e.g., Anderson et al., 1986; Kanamori et al.,
1993). Similar to the 2022 event, the great Colima-Jalisco
earthquake of 1995 (M_w 8.0) had a NW directivity (e.g.,
Courboulex et al., 1997; Hjorleifsdóttir et al., 2018). Mi-
randa y Marron (1911-1912) mentions that the 7 June 1911
earthquake (M_s 7.7), whose location is poorly known but was
in the Michoacán - Colima region, was very strongly felt
in Mexico City, causing considerable damage and leaving
40 persons dead. The intensity of the earthquake in the city
was much stronger than for earthquakes of similar magnitude
that occurred along the coast of Guerrero between 1907 and
1911. Isoseismic contours of the 1911 earthquake are elon-
gated towards the east. Eastward directivity towards Mexico
City provides a plausible explanation for the intensity with
which it was felt in the city. The earthquakes of 2022 and
others events mentioned above once again bring into focus
the importance of source directivity in the recorded and
simulated ground motion in Mexico.

Finally, we find that observing three major earthquakes
(M > 7) on the same day in central Mexico is not so im-
probable.

Acknowledgments

Data used in this study were provided by the Servicio Sis-
mológico Nacional (SSN, Mexican National Seismological
Service), Red Acelerográfica del Instituto de Ingeniería
(IING), Universidad Nacional Autónoma de México (UNAM,
National Autonomous University of Mexico), and Centro de
Instrumentación y Registros Sísmicos (CIRES). We thank
personnel of SSN, IING, and CIRES for station maintenance,
data acquisition and distribution. In the estimation of the
radiated seismic energy and the source spectra, the data pro-
vided by the following network were used: CU (Caribbean
Network, doi: 10.7914/SN/CU), G (French Global Network
of Seismological Broadband Stations, GEOSCOPE, doi:
10.18715/geoscope.g), II (Global Seismograph Network
- IRIS/IDA, doi: 10.7914/SN/II), IU (Global Seismograph
Network, GSN - IRIS/USGS, doi: 10.7914/SN/IU), NL
(Netherlands Seismic and Acoustic Network, doi: 10.21944/
e970fd34-23b9-3411-b366-e4f72877d2c5). The data was
accessed through the IRIS DMC. The research was partially
supported by UNAM, PAPIIT project IN108221 (S.K.S.).
X.P.-C. had a sabbatical fellowship from DGAPA-UNAM
and thanks the Seismological Laboratory at Caltech for
partial funding for her sabbatical.

References

Abe K. (1981). Magnitude of large shallow earthquakes from 1904 to
1980. Physics of the Earth and Planetary Interiors, 27, 72- 92.

Anderson J. G., Bodin P., Brune J.N., Prince J., Singh S.K., Quaas R., Oñate
M. (1986). Strong ground motion from the Michoacan, Mexico, earth-
quake. Science, 233, 1043-1049. DOI:10.1126/science.233.4768.1043.
Arroyo D., Garcia D., Ordaz M., Mora M.A., Singh S.K. (2010). Strong
ground-motion relations for Mexican interplate earthquakes. Journal
of Seismology, 14, 769-785. doi:10.1007/s10950-010-9200-0

Arroyo D., Ordaz M., Singh S.K. (2022). Prediction of Fourier ampli-
tude spectrum of ground motion in Mexico City from interplate thrust
earthquakes, in preparation

Bandy W., Mortera C., Urrutia J., Hilde, T.W.C. (1995). The subducted
Rivera-Cocos Plate boundary: Where is it, what is it, and what is its rela-
tionship to the Colima rift?. Geophysical Research Letters, 22, 3075-3078.
Boatwright J. (2007). The persistence of directivity in small earthquakes.
Bulletin of the Seismological Society of America, 97, 1850-1861.

Boatwright J., Choy G.L. (1986). Teleseismic estimates of the energy
radiated by shallow earthquakes. Journal of Geophysical Research.
91, no. B2, 2095-2112.

Boore D. M., Joyner W. B. (1997). Site amplifications for generic rock
sites. . Bulletin of the Seismological Society of America, 87, 327-341.
Brune J. N. (1970). Tectonic stress and the spectra of seismic shear waves
from earthquakes. Journal of Geophysical Research, 75, 4997-5009.

Calderoni G., Rovelli A., Singh S.K. (2013). Stress drop and source
scaling of the 2009 April L’Aquila earthquakes. Geophysical Journal
International, 192, 260-274, doi: 10.1093/gji/ggs011.

Corona-Fernández, R.D., Santoyo, M.A. (2022). Re-examination of the
1928 Parral, Mexico earthquake (M6.3) using a new multiplatform
graphical vectorization and correction software for legacy seismic data.
Geoscience Data Journal. https://doi.org/10.1002/gdj3.159

Cosenza-Muralles B., DeMets C., Marquez-Azúa B., Sánchez O., Stock
J., Cabral-Cano E., McCaffrey R. (2022). Coseismic and postseismic
deformation for the 1995 Colima-Jalisco and 2003 Tecoman thrust
earthquakes, Mexico subduction zone, from modelling of GPS data.
Geophysical Journal International, 228, 2137-2173, https://doi.
org/10.1093/gji/ggab435

Cosenza-Muralles B., DeMets C., Marquez-Azúa B., Sánchez O., Stock J.,
Cabral-Cano E., McCaffrey R. (2022). GPS-derived interseismic fault
locking along the Jalisco-Colima segment of the Mexico subduction
zone. Geophysical Journal International, 228, 2174-2197, https://doi.
org/10.1093/gji/ggab436

Courboulex F., Singh S.K., Pacheco J.F., Ammon C. (1997). The 1995
Colima-Jalisco, Mexico, earthquake (M_w 8): A study of the rupture
process. Geophysical Research Letters, 34, 1019-1022.

Das S., Henry C. (2003). Spatial relation between main earthquake slip
and its aftershock distribution. Review of Geophysics, 41 (3), 1013,
doi: 10.1029/2003RG000119.

DeMets C., Gordon R.G., Argus D.F. (2010). Geologically current plate
motions. Geophysical Journal International, 181, 1-80.

Dost B., Haak H.W. (2006). Comparing waveforms by digitization and
simulation of waveforms for four Parkfield earthquakes observed in
station DBN, The Netherlands. Bulletin of the Seismological Society
of America, 96, S50-S55, doi: 10.1785/0120050813.

Duputel Z., Rivera L., Kanamori H., Hayes G. (2012). W-phase fast source
inversion for moderate to large earthquakes (1990 - 2010). Geophysical
Journal International, 189(2), 1125-1147.

EERI Preliminary Virtual Reconnaissance Report. (2022). Aquila, Mi-
choacán, Mexico September 19, 2022, M_w 7.6 earthquake. Release date
30 September, 2022, 109 pp.

Hayes G.P., Rivera L., Kanamori H. (2009). Source inversion of the
W phase: real-time implementation and extension to low magnitude.
Seismological Research Letters, 80(5), 817-822.

Hjorleifsdóttir V., Singh S.K., Husker A. (2016). Differences in epicentral
location of Mexican earthquakes between local and global catalogs:
An update, Geofísica Internacional, 55, 79-93.

Hjorleifsdóttir V., Sánchez-Reyes H. S., Ruiz-Angulo A., Ramirez-Herrera
M. T., Castillo-Aja R., Singh S. K., Ji C. (2018). Was the 9 October

1995 M_w 8 Jalisco, Mexico, earthquake a near-trench event? Jour-
nal of Gephysical Research: Solid Earth, 123, 8907-8925. https://
doiorg/10.1029/2017JB014899

Iglesias, A., Singh, S. K., Castro-Artola, O., Pérez-Campos, X., Co-
rona-Fernández, R. D., Santoyo, M. A., Espindola V. H., Arroyo D.,
and Franco, S. I. (2022). A Source Study of the M w 7.0 Acapulco,
Mexico, Earthquake of 8 September 2021. Seismological Research
Letters, 93(6), 3205-3218.

Kanamori H. Anderson D. L. (1975). Theoretical basis of some empiri-
cal relations in seismology. Bulletin of the Seismological Society of
America, 65, 1073- 1095.

Kanamori H., Jennings P. C., Singh S. K., Astiz L. (1993). Estimation of
strong ground motions in Mexico City expected for large earthquakes
in the Guerrero seismic gap. Bulletin of the Seismological Society of
America, 93, 811- 829.

Kanamori H., Rivera L. (2008) Source inversion of W phase: speeding
tsunami warning. Geophysical Journal International. 175, 222-238.

Kelleher J. A., Sykes L.R., Oliver J. (1973). Possible criteria for predicting
earthquake locations and their applications to major plate boundaries
of the Pacific and Caribbean. Journal of Geophysical Research, 78,
2547-2585.

Koketsu K., Miyake H., Guo Y., Kobayashi H., Masuda T., Davuluri
S., Bhattarai M., Adhikari L. B., Sapkota S. N. (2016). Widespread
ground motion distribution caused by rupture directivity during the

2015 Gorkha, Nepal earthquake. Scientific Reports, 6, 28536. http://
doi.org/10.1038/srep28536

Lay T., Ye L., Koper K.D., Kanamori H. (2017). Assessment of teleseis-
mically-determined source parameters for the April 25, 2015 M_w 7.9

Gorkha, Nepal earthquake and the May 12, 2015 M_w 7.2 aftershock. Tec-
tonophysics 714, 4-20, http://dx.doi.org/10.1016/j.tecto.2016.05.023

Maubant L., Radiguet M., Pathier E., Doin M. P., Cotte N., Kazachkina
E., Kostoglodov V. (2022). Interseismic coupling along the Mexican
subduction zone seen by InSAR and GNSS. Earth and Planetary Sci-
ence Letters 586, 117534.

Miranda y Marron M. (1911- 1912). El temblor de 7 de junio de 1911.
Société Scientifique "Antonio Alzate," -Mémoires, 32, 27-66.

Mendoza C., Martinez-Lopez M.R. (2022). Rapid finite-fault analysis
of large Mexico earthquakes using teleseismic P waves. Journal of
Seismology, 26, 333-342. https://doi.org/10.1007/s10950-022-1083-y
Ordaz M., Arboleda J., Singh S. K. (1995). A scheme of random sum-
mation of an empirical Green's function to estimate ground motions
from future large earthquakes. Bulletin of the Seismological Society of
America, 85:1635-1647

Pacheco J., Singh S. K., Dominguez J., Hurtado A., Quintanar L., Jimé-
nez Z., Yamamoto J., Gutiérrez C., Santoyo M., Bandy W., Guzmán
M., and Kostoglodov V. (1997). The October 9, 1995 Colima-Jalisco,
Mexico earthquake (M_w 8): An aftershock study and a comparison of
this earthquake with those of 1932. Geophysical Research Letters 24,
2223-2226. https://doi.org/10.1029/97GL02070

Pérez-Campos X., Beroza G.C. (2001). An apparent mechanism depen-
dence of radiated seismic energy. Journal of Geophysical Research,
106, no. B6, 11,127-11,136.

Pérez-Campos X., Singh S.K., Beroza G.C. (2003). Reconciling teleseismic
and regional estimates of seismic energy. Bulletin of the Seismological
Society of America, 93, 2123-2130.

Pérez-Campos X., Espíndola V.H., Pérez J., Estrada J.A., Cárdenas Mon-
roy C., Zanolli B., Bello D., González-López A., González Ávila D.,
Maldonado R., Montoya-Quintanar E., Vite R., Martínez L.D., Tan
Y., Rodríguez Rasilla I., Vela Rosas M. Á., Cruz J.L., Cárdenas A.,
Navarro Estrada F., Hurtado A., Mendoza Carvajal A.J. (2019). Servicio
Sismológico Nacional, México. Bulletin of the International Seismo-
logical Centre, 53(II), 29-40, https://doi.org/10.31905/SZ7RYBTM
Reyes A., Brune J.N., Lomnitz C. (1979). Source mechanism and after-
shock study of the Colima, Mexico earthquake of January 30, 1973.
Bulletin of the Seismological Society of America, 69, 1819-1840.

Santoyo M. A., Mikumo T., Quintanar L. (2006). Faulting process and
coseismic stress change during the 30 January, 1973, Colima, Mexico
interplate earthquake (M_w=7.6). Geofísica Internacional, 45,3, 163-178.
Seo M.-S., Kim W.-Y., Kim Y. (2022). Rupture directivity of the 2021
ML 2.2 Gwangyang, Korea, microearthquake: Toward resolving high-
resolution rupture process of a small earthquake, The Seismic Record
2(4), 227-236, doi: 10.1785/0320220030.

Singh, S.K., Pacheco, J.F., Alcántara, L., Reyes, G., Ordaz, M., Iglesias,
A., Alcocer, S.M., Gutierrez, C. Valdés, C., Kostoglodov, V., Reyes, C.,
Mikumo, T., Quaas, R., Anderson, J.G. (2003). A Preliminary Report
on the Tecomán, Mexico Earthquake of 22 January 2003 (M_w 7.4) and
its Effects. Seismological Research Letters, 74, 279-289.

Singh S. K., Plata-Martínez R., Pérez-Campos X., Espíndola V.H., Iglesias
A., Arroyo D. (2019). Evidence of directivity during the earthquakes of

8 and 10 May 2014 (M_w 6.5, 6.1) in the Guerrero, Mexico seismic gap
and some implications. Journal of Seismology, 23, 683-697.

Singh S. K., Pérez-Campos X., Ordaz M., Iglesias A., Kostoglodov V.
(2020). Scaling of peak ground displacement with seismic moment
above the Mexican subduction thrust. Seismological Research Letters,
91, 956-966, doi: 10.1785/0220190155.

Singh, S.K., Corona-Fernández, R.D., Santoyo, M., Iglesias, A. (2022).
Repeating large earthquakes along the Mexican subduction zone, in
preparation.

Singh S. K., Ponce L., Nishenko S.E. (1985). The great Jalisco, Mexico,
earthquakes of 1932: Subduction of the Rivera Plate. Bulletin of the
Seismological Society of America, 75, 1301-1313.

Singh S.K., Lermo J. (1985). Mislocations of Mexican earthquakes as
reported in international bulletins. Geofísica Internacional, 24, 333-351.
Singh S.K., Suárez G. (1988). Regional variation in the number of af-
tershocks (mb > 5) of large, subduction-zone earthquakes (M_w > 7.0).
Bulletin of the Seismological Society of America, 78, 230-242.

Somerville P. G., Smith N. F., Graves R. W., Abrahamson N. A. (1997).
Modification of empirical strong ground motion attenuation relations
to include the amplitude and duration effects of rupture directivity.
Seismological Research Letters 68, 199-222.

UNAM Seismology Group (1986). The September 1985 Michoacan
earthquakes: Aftershock distribution and history of rupture. Geophysical
Research Letters, 13, 573-576.

Wang S.-C., McNally K.C., Geller R.J. (1982). Seismic strain release
along the Middle America Trench, Mexico. Geophysical Research
Letters, 9, 182-185.

Ye L., Kanamori H., Lay T. (2018). Global variations of large megathrust
earthquake rupture characteristics, Science Advances 4, 1-8. https://
doi.org/10.1126/sciadv.aao4915

Ye L., Lay T., Kanamori H., Rivera L. (2016). Rupture characteristics
of major and great (M_w 7.0) megathrust earthquakes from 1990-2015:
I. Source parameter scaling relationships. Journal of Geophysical
Research, Solid Earth 121, 826-844. http://dx.doi.org/10.1002/
2015JB012426