Since the earlier times of documented seismological observations, it was noticed that an earthquake (usually a large one) was followed by a sequence of many smaller earthquakes, originating in the epicentral region; the first, larger, earthquake is called the mainshock, or main shock or main event, and the following, smaller, earthquakes are called aftershocks. These sequences, and their spatial and temporal distributions, depend on the characteristics of the mainshock and on the physical properties and the state of stress, strain, temperature, etc., of the region of occurrence (Kisslinger, 1996).
The observation that in many aftershock sequences the magnitude of the largest aftershock is about less than that of the mainshock is known as Båth’s law (Helmstetter & Sornette, 2003; Richter, 1958); it implies an energy ratio between largest aftershock and mainshock. An earthquake following some mainshock but not small enough to be considered an aftershock by the Båth’s law criterion is often considered to be a mainshock in its own right and to constitute, together with the mainshock, a multiple event or multiplet. A large event followed by an even larger one is demoted from mainshock to foreshock; a probabilistic calculation, based on the observation that aftershocks follow the Gutenberg-Richter relation (Gutenberg & Richter, 1954; Kisslinger & Jones, 1991; Shcherbakov et al., 2005), indicates that for the seismicity following a given mainshock has a ~ 6.3% probability of including a “daughter” event larger than the mainshock (Helmstetter et al., 2003; Holliday et al, 2008).
Thus, the simple definition of aftershock as an earthquake occurring after a mainshock and in its epicentral region, although implying some causal relation with the mainshock, is partly semantic and largely circumstantial. Indeed, smallish earthquakes that constitute the background seismicity occur all the time in a seismic region in the absence of large events, and continue occurring whether or not a large earthquake occurs, so that not all earthquakes occurring in the region after a mainshock are necessarily aftershocks.
If aftershocks are a result of the occurrence of the mainshock, then they should be related in a physical way with its rupture process. The aftershock-producing mechanism is not yet known, but it is conceivably related with adjustments to the post-mainshock stress field (Lay & Wallace, 1995) possibly through viscolelastic processes or through fluid flow (Nur & Booker, 1972); whichever the actual process, aftershocks should be related with the rupture plane of the mainshock. Kisslinger (1996) qualitatively defines three kinds of aftershocks: Class 1 is those occurring on the ruptured area of the fault plane or on a thin band around it. Class 2 is events that occur on the same fault but outside of the co-seismic ruptured area. Class 3 is events occurring elsewhere, on faults other than the one ruptured by the mainshock; these events, whether in the same region or not, will not be considered here as aftershocks, but rather will be classified as triggered earthquakes.
where t is time measured from that of the mainshock, and k, c, and p are positive parameters which vary with the lithologic, tectonic, and other conditions of the study area. Commonly, p ranges from 0.9 to 1.8, with most instances between 1.1 and 1.4 (Utsu, 1961); values which do not show dependence on the magnitude of the mainshock (Utsu, 1962).
Aftershocks occurring within 24 to 48 hours after a large earthquake locate mostly over the co-seismic rupture area, and provide a good estimation of it (Lay & Wallace, 1995), which indicates that seismicity at the time is mostly class 1; over longer times the aftershock area increases (Felzer & Brodsky, 2006; Helmstetter & Sornette, 2002; Mogi, 1968; Tajima y Kanamori, 1985; Valdés et al., 1982), including at the edges class 2 events.
Among other reasons why aftershock identification is important, we can mention the following few examples. Aftershocks can give important information about the rupture area; also, from estimations of co-seismic slip on the fault plane by inversion of seismic waves, several authors have found that aftershocks are scarce in areas of maximum slip and concentrate around their edges (Dreger et al., 1994; Engdahl et al, 1989; Hauksson et al., 1994; Mendoza & Hartzell, 1988), so that aftershocks give information about the rupture process of the mainshock. Aftershocks can also yield information about the properties of the epicentral region (Knopoff et al., 1982; Kisslinger, 1996; Kisslinger & Hasegawa, 1991; Figueroa, 2009), and about possible triggering mechanisms (Roquemore & Simila, 1994).
Since aftershocks can be large enough to contribute to the damage (particularly after structures have been debilitated by a mainshock), it is important to estimate the hazard associated with aftershock activity (Felzer et al., 2003; Reasenberg & Jones, 1989).
It has been proposed that, since aftershock activity depends among other factors on the stress status of the region, there is research on whether some characteristics of the aftershock activity from intermediate-sized earthquakes can be useful as precursory data for large earthquake hazard estimation (e.g. Jones, 1994; Keilis-Borok et al., 1980; Wyss, 1986).
Finally, for some studies concerned with large earthquakes, aftershocks can be considered as noise, and have to be eliminated from the catalogs (e.g. Gardner. & Knopoff, 1974; Habermann & Wyss, 1984).
In order to use or eliminate aftershocks it is first necessary to identify them. Many methods have been used, ranging from visual inspection (Molchan & Dmitrieva, 1992) to sophisticated numerical techniques. A common method identifies aftershocks as those shocks locating within temporal and spatial windows having lengths which usually depend on the magnitude of the mainshock (Gardner & Knopoff, 1974; Keilis-Borok et al., 1980; Knopoff et al., 1982). More sophisticated methods identify aftershocks as belonging to a spatial cluster, consisting of events within a given distance of at least one other event belonging to the cluster, which includes the mainshock (Davis & Frohlich 1991a,b; Frohlich & Davis, 1985). A variation of the window method considers events from larger to smaller magnitudes with the size of the spatial windows a function of the magnitude, and density of events (Prozorov & Dziewonski, 1982; Prozorov, 1986). Other methods include recognizing some statistical property (e.g. Kagan & Knopoff 1976, 1982; Vere-Jones & Davies 1966) or interpreting the relations between events according to some statistical chain or branching model (Molchan & Dmitrieva, 1992; Reasenberg, 1985).
Our method includes some of the above mentioned techniques used to discard events which cannot be aftershocks, and then proceeds to identify aftershocks based on the physical model of a rupture plane and on recognized statistical relationships. An early unsophisticated application of the rupture plane model, which proved that this principle of aftershock identification was feasible, was part of an unpublished MSc. thesis (Granados, 2000).
2. The method
We work with seismic catalogs containing occurrence time (days), hypocentral x (East), y (North), and z (up), and, optionally, horizontal and vertical location uncertainties and, all these in kilometers. If location uncertainties are not included in the catalog, optional horizontal and vertical uncertainties are assigned equally to all events.
Any events occurring before the event with the largest magnitude, the mainshock, are eliminated. All spatial coordinates are then referred to those of the mainshock.
A rough time cutoff, eliminates events occurring after more than an optional cutoff time (default is), because after a few years it is difficult to distinguish aftershock activity from background seismicity.
The extent of the aftershock area depends on the energy, i.e. on the magnitude, of the mainshock (Utsu, & Seki, 1955), as does the rupture area. A first rough spatial discrimination, based on an average of the empirical magnitude M vs. source-length r (km) relationships of Wells and Coppersmith (1994):
or (Kagan, 2004):
eliminates events farther away from the hypocenter than 1.5 times the length estimated by (2).
Next, a spatial clustering analysis, where events separated by no more than a given critical distance r, of the order of hundreds of meters to a few kilometers, depending on the spatial coverage of the catalog, are considered to be related, is used to eliminate events which do not relate to the mainshock or to other possible aftershocks.
The parameters in the modified Omori’s law (1) are not known a priori; they are estimated from the statistics of the aftershocks (Davis & Frohlich, 1991b; Guo & Ogata, 1997; Ogata, 1983), which we do not yet have. However, this relation tells us that for long enough times after the occurrence of the mainshock the number of aftershocks decreases until seismic activity in the epicentral region returns to its background level (Ogata & Shimazaki, 1984). When aftershock occurrence shows gaps comparable to those characterizing the background seismicity, we can consider that the aftershock activity is, if not ended, at least scarce enough to be comparable to the background activity and can no longer be distinguished from it. The critical gap length depends on the region and the magnitude threshold of the observations; we use a default critical gap length. The gap is measured as the average of a given number of inter-event times (default). Events occurring after a critical gap are discarded from the possible aftershocks.
Weights, based on relation (1), are optionally assigned to the remaining events, using typical values for c (default) and p (default) and by setting, so that those events which have a large likelihood of being aftershocks have weights ~1, while later events have smaller weights:
Next, plane fitting is carried out iteratively; at each iteration, a plane that passes through the mainshock hypocenter is fitted to all remaining events, through a genetic scheme described below, and fit outliers (events too far away from the plane) are discarded. Iteration continues until the goodness-of-fit criterion is met (successful fit) or until a preset maximum number of iterations is attained (unsuccessful fit).
The, c, and p values can be refined using the final results of a first, tentative aftershock determination, to do a second one.
The genetic plane search for the plane, characterized by its azimuth and dip, which minimizes the L1 norm of the perpendicular distances from the fault plane to all events, is as follows. An initial set of equispaced and values, covering the acceptable ranges, is built and the fit to the aftershock candidates is evaluated for each pair of values.
To estimate the error of fit for each candidate plane, for each azimuth and dip pair, event coordinates are rotated as:
and in the coordinate system the equation of the plane is, so that the L1 fit error for the j’th parameter pair is easily computed as
where is the number of remaining aftershock candidates.
The parameter pairs corresponding to the best fits are chosen as the parents of the next generation; they are sorted by the absolute value of their respective errors, so that parent number 1 is the best fit and parent number has the largest error. The standard deviations of the parent’s parameters, and, are evaluated.
Next, children are constructed with parameters which are averages of those of each pair of parents; other children are constructed for each parent j by randomly modifying one or the other of the parameter values, according to
where N designates the normal distribution, , , and is a minimum allowable value that ensures significant variations.
Errors are computed for the children and the best fits among the whole population, parents plus children, are chosen as the parents for the next generation. The process is repeated until the goodness of fit criterion is met (and the process ends) or until a preset number of generations is attained.
For the current iteration the best fit is for the plane corresponding to parent number one, and; using these values event coordinates are rotated as in (4) and the standard deviation, , of the distances is evaluated. Those events with
where is a damping factor (default is 1.25) and is the location uncertainty, adjusted for the plane dip as
are eliminated as outliers.
This method has been implemented as a Matlab program, aftplane.m, which allows the user to interactively set most parameters, 3D plots the input hypocenters, identifies the aftershocks of one mainshock, optionally plots the aftershocks and/or independent shocks (in different colors), and optionally outputs the corresponding catalogs. A definite advantage of using Matlab for this algorithm is that both data and trial models are handled as matrices, so that rotating, sorting, and identifying values is done more efficiently and with less lines of code than would be possible in other programming languages like FORTRAN or BASIC.
A variation of this method is used as a function by program cleancat.m which identifies and eliminates the aftershocks of all mainshocks in a catalog. For each event not previously identified as an aftershock, aftershocks are identified following the same steps described above, except that, after fitting the plane a search is made for events with magnitudes, which will not be considered as aftershocks, and if any are found, the aftershock list is cut so as to exclude all events beginning with the first of these non-aftershocks. Thus, for mainshocks occurring at times, instead of the total Omori number of aftershocks (Utsu, 1970; Ogata, 1983)
where is the Heavyside function, we are actually evaluating
is the boxcar function, and is the time of the succeeding event occurring along the same fault plane and cluster or the time when the gap criterion is met (if no such event or gap occur). An event at may “inherit” some of the aftershocks from the one at, but they will be identified as aftershocks and be duly eliminated. Båth’s law indicates that should be 1.2, but this criterion results, for most catalogs, in too many mainshocks; our default value is thus, which means that we accept as aftershocks those with energy ratio, but this value can be easily changed by the user. The largest events in the catalog are plotted vs. time above the initial and resulting event densities, to illustrate which aftershocks were eliminated. The program iterates the whole process, as many times as needed, until no more aftershocks are found
In cleancat, parameters are not set interactively, but can be easily adjusted in a list of adjustable parameters at the beginning of the code.
The program optionally outputs a catalog excluding identified aftershocks.
We will now show some examples of the application of the method. The aftplane program will be used to identify fault planes and aftershocks from three mainshock-aftershock sequences from two different parts of the world featuring different tectonic environments. The cleancat program will be used to clean the catalog of a fault system.
3.1. Aftplane: transcurrent regime, Joshua Tree and Landers earthquakes
In 1992, two large earthquakes, the Joshua Tree, April 23 6.2, 33.9°N, 116.3°W and the Landers, June 28 7.3, 34.2°N, 116.5°W, occurred close together on a line previously unrecognized as a potential throughgoing seismogenic fault (Nur et al., 1993). We chose these events as illustration because, although both events have mainly strike-slip mechanisms, they have slightly different strikes and dips, so that we wanted to test whether the method could identify these small differences.
Figure 1 shows the location of the study area in souther California, USA, and its recent seismicity; the faults ruptured during the Joshua Tree and Landers earthquakes are located within the red diamond.
For both events we used maximum allowable horizontal and vertical location uncertainties of 0.2 km and0.5 km, respectively; a priory cutoff times; maximum distance for spatial clustering km; maximum permissible time gap estimated as the average of inter-event times; Omori weighting parameters days and; aftershock magnitude criterion; fit criterion maximum error km.
3.2. Aftplane: Joshua Tree earthquake
The catalog for the Joshuea tree earthquake contained 5075 events spanning 66.30 days (~0.182 yr). Figure 2 (left) shows the 3497 remaining events after the first rough elimination by acceptable uncertainties and by an estimated expected fault length of ~7.97 km corresponding to a critical distance of 11.96 km. Figure 2 (right) shows as blue circles the 3379 shocks identified as clustering with the main event.
The main Joshua Tree shock and the identified 1094 aftershocks are shown in figure 3, both in a plan view (left) which clearly shows the resulting 171.6º faultplane azimuth, and a cross section along the fault plane azimuth (right) which shows the resulting 86.6º faultplane dip. The values found by aftplane agree extremely well with those estimated by Velasco et al. (1994) of strike 171º, dip 89°. Figure 4 shows a cross section parallel to the fault plane, illustrating aftershock concentrations.
3.2.1. Aftplane: Landers earthquake
The catalog for the Landers earthquake contained 49,605 events spanning 4,932.52 days (~13.514 yr). Figure 5 (left) shows the 17,553 remaining events after the first rough elimination by acceptable uncertainties and by an estimated expected fault length of ~76.78 km corresponding to a critical distance of 115.17 km. Figure 5 (right) shows as blue circles the 12,834 shocks identified as clustering with the main event.
The main Landers shock and the identified 3,225 aftershocks are shown in figure 6, in a cross section seen along the determined 340.6º fault plane azimuth, which shows the resulting 70.1º faultplane dip. The values found by aftplane agree extremely well with those estimated by Velasco et al. (1994) of strike 341°º, dip 70°. Figure 7 shows a cross section parallel to the fault plane, illustrating aftershock concentrations.
3.3. Aftplane: subduction regime, Armería (Tecomán, Colima) earthquake
The Armería 7.4 earthquake, also known as the Tecomán or Colima 2003 earthquake, occurred on January 22, 2003, on the subduction zone along the boundary of the Northamerican and Pacific plates in central-western Mexico.
Figure 8 shows the location of the study area, the mainshock epicenter (red star) and the subsequent seismicity recorded and located by the Colima Seismic Network (RESCO).
Nuñez et al (2004) and Yagi et al (2003) estimated a fault plane with a 300° strike and a quite shallow 20° dip, which agrees with the 20° to 30° dip of the subduction zone determined by Andrews et al (2010).
For aftplane we used the RESCO catalog with the same parameter values mentioned above, except for horizontal and vertical location uncertainties of 0.075 km and0.50 km, and maximum distance for spatial clustering km;
The catalog for the Armería earthquake contained 11,475 events spanning 1,529.9 days (~4.192 yr). Figure 9 (left) shows the 10,275 remaining events after the first rough elimination by acceptable uncertainties and by an estimated expected fault length of ~92.73 km corresponding to a critical distance of 139.09 km. Figure 9 (right) shows as blue circles the 7,109 shocks identified as clustering with the main event.
The main Armería shock and the identified 460 aftershocks are shown in figure 10, in a cross section seen along the determined 86.2º fault plane azimuth, which shows the resulting 33.2º faultplane dip. The values found by aftplane agree extremely well with the above mentioned estimates strike 300° and 20° to 30° (Nuñez et al., 2004; Yagi et al., 2004; Andrews et al., 2010). Figure 11 shows a cross section parallel to the fault plane, illustrating aftershock concentrations.
3.4. Cleancat: whole Joshua Tree-Landers fault system
To illustrate the use of program cleancat we chose the catalog covering the whole Joshua Tree-Landers fault system (Figure 1), because this is a system with many close-lying, subparallel, faults, which gives scope to the iterative aftershock recognition scheme of the program.
The parameters used are 0.2 km and0.5 km; a priory cutoff times; maximum distance for spatial clustering km; maximum permissible time gap estimated as the average of inter-event times; Omori weighting parameters days and; aftershock magnitude criterion; fit criterion maximum error km.
Figure 11 shows all events in the catalog (black crosses), and identified aftershocks as yellow circles. Total processing consisted of 10 iterations which identified and eliminated 11,665, 4,212, 1,702, 86, 94, 30, 80, 18, 49, and 1 aftershocks, respectively, for a total of 17,937 aftershocks.
Figure 12 shows the occurrence times and magnitudes of the largest events in the catalog (top), and below them (middle) is plotted the ocurrence density (per day) vs. time (blue line); peaks in the occurrence rate after large shocks aftershocks are clearly seen. The bottom panel of figure 12 shows the occurrence density after aftershock elimination using (red line); although densities are much lower than before filtering, the peaks coinciding with the occurrence of the Joshua Tree and Landers earthquakes indicate that many aftershocks are not being identified.
Use of (Fig.13) effectively diminishes the troublesome above mentioned peaks to background seismicity level; the seismicity rate peak arounddays is associated with distributed seismicity with events about the same size. Thus, we see that Båth’s law is not appropriate for the mainshock-aftershock relationships in the seismicity of the Joshua Tree-Landers fault system; a small (0.1 unit) change in the magnitude criterion can make a large difference in the aftershock recognition capability.
We present a simple method for identification and/or elimination of aftershocks, based on the generally accepted assumption that aftershocks are related to the fault rupture of the mainshock. The method has been tried on various catalogs with good results and, when aftershocks are numerous enough, good estimates of rupture planes that agree very well with those reported in the literature.
A variation of the method used for eliminating all aftershocks from a seismicity catalog (“catalog cleaning”) uses, iteratively, a variation of the same principle. Using the seismicity occurrence time rate as illustration and criterion of the effectiveness of the method, indicates that the required difference between mainshock and aftershocks is a key parameter for correct aftershock identification, and that (Båth’s law) may be too strict for some geographic areas and/or seismo-tectonic settings.
Many thanks to José Frez, Juan García A., and María Luisa Argote for useful criticism and comments. We are grateful to RESCO and Gabriel Reyes, and to the SCEC for the use of their catalogs. We also thank Mr Davor Vidic of Intech for the kind invitation to participate in the present book.