# Seismological Research Letters

- © 2008 by the Seismological Society of America

## INTRODUCTION

A shallow moderate earthquake of MN (Nuttli magnitude) 4.1 occurred near Sudbury, Ontario, at 07:22:55 UTC (02:22:55 local time) on 29 November 2006. The event was located 6 km from Lively, Ontario, and was induced by mining activity at Inco's Creighton mine. Inco's in-mine monitoring system, and its subsequent visual inspection, places the event at 46.478 N 81.203 W, at a depth of 2.35 km. This event offers us a rare close-up view of a moderate earthquake, because it was recorded at close distances by five three-component broad-band POLARIS (Portable Observatories for Lithospheric Analysis and Research Investigating Seismicity; http://www.polarisnet.ca) seismographic stations that are installed above and within the Creighton mine. Two stations are located on the surface, near the epicenter, while three subsurface stations are at depths of 1 to 2 km. The closest station to the event is within the Sudbury Neutrino Observatory (SNO) lab, which is 2 km beneath the surface and about 300 m above the event. The event caused no damage to the SNO lab and only minor damage to the mine workings; there were no injuries. The event was widely felt in the Sudbury, Ontario, region, as expected based on instrumental ShakeMaps (http://www.shakemap.carleton.ca) and confirmed by Internet intensity reports (http://www.earthquakescanada.ca). The mainshock was followed by two significant aftershocks at 2:36 (MN 2.0) and 2:38 (MN 3.1) local time and by many smaller aftershocks. In addition to the near-source records, this event was well-recorded regionally at distances from 20 km to hundreds of kilometers by the seismographic stations of the Canadian National Seismographic Network (CNSN) and POLARIS. Figure 1 shows the location of the Lively event in relation to regional seismographic stations; the nearest regional station is SUNO, 22 km away. In this paper we provide a preliminary examination of the recorded ground motions from this event and discuss their implications for our understanding of source and attenuation properties in eastern North America (ENA).

## TECTONIC SETTING AND EARTHQUAKE SEQUENCE

### Tectonic Setting

One of the world's deepest underground mines (> 2 km), the Creighton mine has produced more than 160 metric tons of nickel-copper ore since opening in 1901. The mine is located along the southern periphery of the Sudbury Structure, interpreted by most researchers as a relict giant impact (*e.g.*, Riller 2005; Grieve and Therriault 2000). The structure consists of a differentiated melt sheet (Sudbury Igneous Complex) overlain by impact breccias and post-impact sedimentary rocks (Pye *et al.* 1984). The structure, together with older Precambrian rocks of the Canadian Shield, was deformed by post-impact orogenic processes to produce a large-scale asymmetric synform (Milkereit *et al.* 1992). Like most other ore deposits in the Sudbury region, the Creighton Mine is situated near the base of the impact melt sheet, near its contact with the footwall. The SNO laboratory, where the closest seismic instruments are located, is situated within a massive igneous unit (norite) within the Sudbury Igneous Complex.

### The Lively Sequence

The Lively earthquake sequence consisted of the MN 4.1 mainshock and a sequence of aftershocks. Accurate characterization of aftershock statistics is important for hazard calculations, particularly in a mining environment where they may determine protocols for re-entry into a mine following a significant mining-related earthquake. In this case, two moderate aftershocks occurred 13 and 15 minutes after the mainshock. Although these larger aftershocks were observed at regional distances, most other aftershocks were too small to be observed at regional distances and were recorded only by local instruments. The close proximity of the deep underground seismometer at SNO enabled the measurement of a rich record of aftershock activity at the MN > 0 level that would otherwise be largely undetected. By examining the raw seismic data from the deep SNO site, we measured the arrival time and amplitude of 70 aftershocks that occurred within a 12-hour period after the mainshock.

Figure 2 shows the rate of observed aftershocks during the initial 12-hour period following the mainshock, measured using a moving 1-hour time window. We fit the observed trend using a modified form of Omori's law (Shcherbakov *et al.* 2004), (1) where *t* is time since the mainshock, *M* is magnitude threshold (here estimated to be MN ∼ 0), τ(*M*) is the reciprocal of the initial rate of aftershocks, *c*(*M*) is the characteristic decay time, and *p*(*M*) is a parameter that controls the rate of decay of aftershocks. We used a grid search method to find the values of τ, *c*, and *p* that minimize the least-squares misfit with respect to the observed aftershock rate. The resulting curve (figure 2) fits the overall aftershock trend reasonably well but fails to reproduce several conspicuous peaks in aftershock activity that occurred 2-8 hours after the mainshock. In order to fit this aftershock decay pattern, a more complex model is required, such as the epidemic-type aftershock sequences model (ETAS), which includes triggering events within a cascading sequence (Ogata 1988, 1999). We remark that, of the three parameters needed for the generalized Omori representation of aftershock rate, only τ is well-constrained by the data considered here. Large tradeoffs between *c* and *p* mean that other combinations of these parameters provide nearly equivalent fits to the data. Thus, although more work is needed to characterize the aftershock model more accurately, these preliminary results illustrate the potential insights into aftershock processes that may be afforded by the near-field aftershock recordings obtained here.

## GROUND-MOTION DATA FROM THE LIVELY EARTHQUAKE

### Instrumental Data

Figure 1 shows the locations of regional stations in relation to the Lively earthquake. Regional stations typically sample at 40 to 100 Hz. The five SNO stations are located beneath the symbol marking the event location; the sample rate for these local stations is 200 samples/s. The two local stations on the surface (SSNO and 11SNO) are sited on hard rock, while the three subsurface stations (DSNO, LSNO, and 46SNO) are within underground openings, all in hard rock, that are > 4 m in diameter (both directions) and approximately 3 m in height. Table 1 shows the distances from the MN 4.1 event to SUNO and the five SNO stations. Figure 3 shows plots of example seismograms. Note that the very short duration of the records at SNO (about 1 s) limits the analysis of the local seismograms to frequencies greater than ∼ 2 Hz. Furthermore, the absolute timing of the underground stations is not accurate, due to the lack of availability of a GPS time signal underground.

Regional and local data were processed using procedures described by Atkinson (2004). Briefly, the window of strongest shaking (shear window, including direct, reflected, and refracted phases) was selected for each record, and a 5% taper was applied at each end of the window. The Fourier spectrum of acceleration was determined, correcting for instrument response. (Fourier displacement spectra were also calculated.) The spectra were smoothed and tabulated in increments of 0.1 log frequency units, for log frequencies of 0 to 1 (*e.g.*, 1 Hz to 10 Hz). Instrument-corrected acceleration time series were used to compute response spectral amplitudes (5% damped pseudo-acceleration, PSA). In computing the instrument-corrected time series, the records were baseline corrected and Butterworth filtered (4th order). Filter frequencies for the regional records are 0.1 Hz to 50 Hz, while the local records were filtered from 2 Hz to 90 Hz, due to their higher sampling rate and short duration.

### Felt Intensity Data

The Lively event was widely felt in the Sudbury region, with the maximum reported individual estimate of modified Mercalli intensity (MMI) being VI (6). One hundred twenty-four observations were reported to the Geological Survey of Canada's community Internet intensity program (http://www.earthquakescanada.ca), which is a minor modification of the U.S. Geological Survey's “Did You Feel It?” program (modifications allow users to report location by postal code rather than ZIP Code, for example). Figure 4 shows the intensity observations. These were binned by distance to establish stable mean trends, as discussed by Atkinson and Wald (2007); the mean binned intensities near the source are approximately IV (4). The mean binned intensities were used to infer MMI-based estimates of PGA and PGV by inverting the following equations of Atkinson and Kaka (2007) to solve for PGA and PGV: (2a) (2b)

## SOURCE AND ATTENUATION PROCESSES

### Analysis of Fourier Spectra

The processed spectral data are used to evaluate the overall source and attenuation processes for the Lively event. We can correct the Fourier acceleration spectral amplitudes at each station back to the reference distance of 1 km, using the trilinear spectral ENA attenuation model of Atkinson (2004). This provides an apparent source spectrum from each station, given by: (3a) (3b) (3c) where *r* is hypocentral distance, *A _{obs}*(

*f*) is the observed acceleration spectrum at the station, and

*c*

_{4}(

*f*) is the coefficient of anelastic attenuation from Atkinson (2004). (Note: The anelastic attenuation coefficient

*c*

_{4}is of positive sign in this formulation, increasing from 0 at

*f*= 0.2 Hz to 0.00271 at

*f*= 20 Hz. The geometric spreading coefficient is 1.3 in the first 70 km.) The values of log

*A*(

_{src}*f*) are averaged over all stations to obtain the mean apparent source spectrum and its standard deviation. Only stations sited on hard rock, on the surface, were used to obtain the apparent source spectrum; the subsurface recordings will be compared to this spectrum later.

Figure 5 shows the apparent source spectrum of acceleration (*r* = 1 km) for the event, for both the horizontal and vertical components. The displacement spectrum (obtained by dividing the acceleration spectrum by (2π*f*)^{2}) is also shown for the vertical component. The vertical component is interpreted to be representative of the unamplified horizontal component (Atkinson 2004; Beresnev and Atkinson 1997; Siddiqqi and Atkinson 2002), as the amplification on the vertical component is considered negligible for hard-rock sites. The horizontal component shows minor amplification relative to the vertical, due to the impedance effects of the crustal shear wave velocity profile; velocities at depth (3 km) are likely about 3.5 km/s, while those near the surface are typically 2 to 3 km/s (Beresnev and Atkinson 1997); this would produce an amplification of about a factor of 1.3, according to the quarter-wavelength approximation (Boore and Joyner 1997). The observed spectrum follows the simple omega-squared model (Brune 1970) well. From the vertical-component Fourier displacement (FD), we infer a long-period spectral level of 0.02 cm/s (logFD = -1.7). The seismic moment (*M*_{0}) can be estimated from Boore (2003) as (4) where (Note: The constants in the numerator represent radiation pattern, partition onto two horizontal components, and free surface amplification, respectively.) We assume ρ = 2.8 g/cm^{3} and β = 3.5 km/s. The moment magnitude (**M** = 2/3log*M*_{0} - 10.7; Hanks and Kanamori 1979) obtained from the displacement spectrum is then **M** = 3.8 (*M*_{0} = 4.76 × 10^{21} dyne-cm). This is in good agreement with the moment tensor solution for this event, presented in the next section. Under the Brune (1970) model, the acceleration spectra can be written (Boore 2003) as (5) where *f*_{0} is the corner frequency, which is related to the stress drop by (Boore 2003): (6) for Δσ in bars, *M*_{0} in dyne-cm, and β in km/s. Thus we can use the high-frequency level of the acceleration spectrum, along with the inferred moment, to estimate stress drop. Based on the vertical component, we would estimate a high-frequency level of approximately 0.8 log units. The horizontal component high-frequency level is about 1.1 log units, but it has likely been amplified by about a factor of 1.3 (0.1 log units) by the impedance gradient. Overall then, our best estimate for the high-frequency level is 0.9 log units, resulting in an estimated corner frequency of 2.86 Hz and a stress drop of 22 bars using equations 3a, 3b, 3c and 4. In comparison, studies of regional ENA earthquakes of small-to-moderate magnitude, in which the stress drop is calculated using the same procedure as that described above, typically suggest stress drops of 100 to 200 bars for events of **M** > 4 (Atkinson 2004). However, these same studies also reveal a decreasing trend in stress drop value with decreasing magnitude, with typical values for an event of **M** 3.8 being in the range of 30 to 80 bars. Thus the value of 22 bars obtained for the Lively event is a relatively low stress drop, but not anomalous for an event of this magnitude.

An important issue in evaluating the robustness of the apparent source spectra is the goodness-of-fit of the assumed attenuation model to the trends actually shown by the data. Figure 6 shows plots of the inferred source levels (by equation 1) for each station as a function of distance, for two selected frequencies. There are no trends of the inferred source levels with distance, which corroborates the assumed attenuation model. In particular, the source levels inferred from the closest stations are not significantly different than those inferred from stations hundreds of kilometers away.

### Moment Tensor Analysis

A moment magnitude of 3.8 is inferred for this event based on the long-period level of the apparent source spectra at local to regional distances. A more robust measure of the seismic moment can be obtained from a moment tensor inversion. Figure 7 shows the moment tensor inversion, based on modeling of the waveforms at five stations. A standard eastern Canada velocity model (Brune and Dorman 1963, with crust thickened to 40 km) was used at all stations. At each station, a frequency range was selected for the modeling; the signal to noise ratio was checked, and the inversion code used only those frequencies that have a signal to noise ratio > 2. Each component was assigned a weight in the inversion based on its judged quality. For example, VLDQ is weighted 1, 1, 0, which means that the vertical and radial are given full weight and the tangential is not used; the tangential component was discounted because it looks noisy, although the actual fit to this component is similar to that achieved for the other components. In displaying the results, the misfit is calculated even for the components that weren't used in the inversion, which tends to make the misfit appear large. If the misfit is recalculated using only the components that were used in the inversion, it is reduced from 0.69 to 0.36, which is considered a reasonably good fit. The moment magnitude of 3.7 derived from the moment tensor inversion is in good agreement with the estimate of 3.8 obtained from the regional seismographic data.

On the focal mechanism plot of figure 7, the dashed lines are for the best fit double couple, and the shading is for the complete moment tensor solution. The CLVD (compensated linear vector dipole) component, which gives an indication of the non-double-couple component of the source, is a bit high. The high CLVD could be an indication that there was a non-double-couple component to the rupture, possibly related to the mining environment. However, caution should be taken not to overinterpret it. The CLVD could also result from noise in the data, not unrealistic given that we are modeling long periods for a relatively small earthquake. The best depth from the moment tensor inversion is about 7 km, according to the misfit plot as a function of depth on figure 7 (including all components). However, the depth is not well-constrained. Moment tensor inversions are often not good at resolving very shallow depths because of the longer periods modeled (Ekström and Dziewonski 1985). In this case, the depth of the event is known from the mine monitoring and inspection. The depth can be more accurately estimated from seismic records based on the focal depth modeling method discussed in the next section. Note that if the depth is forced to be shallower in the moment tensor inversion, the best fit occurs for a 3-km depth, and the calculated thrust component will be higher (figure 7 inset). The non-double-couple component also increases to 49%, which would normally be an indication of a poorer fit but for reasons previously discussed should be interpreted with caution in this case. Comparing the misfit on a station by station basis, by forcing the solution to a shallower depth there is no change (no change being defined as less than 0.01 difference in either direction) for OTT, KAPO, and BUKO; a slight improvement (0.03) for 3 km at SADO; and a stronger preference (0.21) for 7 km at VLDQ. A comparison of individual component mis-fits reveals a similar distribution, with no change for VLDQ-R, OTT- all components, SADO-Z, KAPO-R, and BUKO-R. A 3-km depth improves the fit at SADO-R and T, BUKO-Z, and KAPO-T with a mean improvement of 0.09. The 3-km depth deteriorates the fit at VLDQ-Z and T, BUKO-T, and KAPO-Z with a mean difference of 0.22. In summary, it appears that while the 7-km depth does provide a better overall fit to the long-period waveforms, the well-constrained shallower depth determined from the local data does not lead to a significantly worse fit. We also note that the thrust mechanism obtained using the shallower depth is more similar to the focal mechanisms of other earthquakes in western Quebec and eastern Ontario than is the predominantly strike-slip mechanism obtained with the depth of 7 km.

### Focal Depth Analysis from Regional Waveforms

The depth of the Lively event is well-constrained at 2.35 km by very local seismic records in the mine and has been confirmed by inspection of the focal area by Inco. It is nevertheless of interest to see how well the depth can be inferred from regional seismographic information, as this bears on the reliability of methods of estimating depth for the more usual case where we have few or no near-source records. A particularly promising method of estimating focal depths in regions of sparse instrumental coverage is the regional depth phase modeling approach described by Ma and Atkinson (2006). This method is based on modeling the regional depth phases *sPg, sPmP*, and *sPn*, with their reference phases *Pg, PmP,* and *Pn*, respectively. These phases are often clearly visible on regional seismograms. The method is based on the calculation of synthetics using the reflectivity method (Randall 1994) with a default focal mechanism and a default crustal model (*e.g.*, Mereu *et al.* 1986). By comparing the synthetics with the observations for selected station distances over a reasonable range of focal depths, the depth that provides the best agreement (in relative arrival times) between the synthetics and observations can be identified. For this event, depth information comes from the *sPg* phase. We selected vertical displacement waveforms recorded at distances from 162 to 296 km that clearly show this phase. We simulated synthetic waveforms in a suitable range of depths, in order to choose the depth at which the synthetic best matches the selected waveform. We used four seismic stations, EEO, BUKO, SADO, and CRLO to calculate the focal depth. The best fit for the synthetic with the observed data is at a depth range of 2.5-3.0 km, just slightly larger than the known depth of 2.35 km. Figure 8 illustrates the fit for one of the stations (EEO). A well-developed Rg phase on many stations (including EEO) also supports the shallow depth.

### Effects of Near-field Source Terms

In an elastic medium, the *n*th component of ground velocity due to a point-dislocation source may be written as (7) where *r* is distance; α, β, and ρ are the *P*- and *S*-wave velocity and density of the medium; *M _{pq}* is the moment tensor;

*s*(

*t*) is the source time function; and

*a*,

_{pq}*b*,

_{pq}*c*,

_{pq}*d*, and

_{pq}*e*are geometrical terms (see Aki and Richards 1980, 79 for details). The first three terms in this expression represent near-field terms that decay as 1/

_{pq}*r*

^{4}and 1/

*r*

^{2}respectively. These contribute significantly to the total ground motion only within a few wavelengths of the source and are generally neglected in the creation of synthetic seismograms. The last two terms decay as 1/

*r*and represent far-field terms that are normally considered for modeling ground motion.

In the case of near-field records of the Lively earthquake, ground-motion observations are available at a fraction of the *P* and *S* wavelengths. For the dominant frequency produced by this earthquake (∼ 1.25 Hz) at the closest station (DSNO), near-field terms are similar in amplitude to the far-field terms. Figure 9(A) shows the unfiltered vertical-component recording of ground velocity at the DSNO station caused by the mainshock, after removal of the mean and normalization by its maximum value. To approximate the long-wavelength part of the source time function (figure 9B), we used (8) with a period *T* = 0.8 s. Based on the moment tensor obtained from regional waveform modeling, we computed the near-field ground velocity without (figure 9C) and with (figure 9D) inclusion of the near-field terms. Since the seismograph is close to the focal depth, even the vertical record is dominated by shear-wave signals. The far-field shear-wave produces a waveform with the shape of the time derivative of the source time function given above. This waveform does not provide a satisfactory fit to the observed ground motion. However, interference between this and the near-field terms results in initial negative motion followed by a longer-duration positive signal (figure 9D), consistent with the observed data (figure 9A).

This example illustrates the significance of near-field terms in the ground-motion record at stations very close to the earthquake source. More detailed studies using inversion techniques may permit us to constrain the source time function more precisely.

## RESPONSE SPECTRAL AMPLITUDES

The response spectral amplitudes, showing how a single-degree-of-freedom oscillator would respond to the ground shaking, are plotted as a function of distance on figure 10. The horizontal component of pseudo-acceleration (PSA) for 5% damping is shown in comparison to various prediction equations that are used in seismic hazard and ShakeMap applications. Estimated PGA and PGV based on binned MMI observations (using the equations of Atkinson and Kaka 2007 [equations 2a, 2b] to convert MMI to PGA and PGV) are also plotted and agree well with actual instrumental data.

The response spectral amplitudes are overpredicted by the Atkinson and Boore (2006) ground-motion relations for hard rock in ENA. This is not surprising, as the event had a stress drop of only 22 bars, while the Atkinson and Boore (2006; hereinafter referred to as AB06) relations are based on a median stress drop of 140 bars. The AB06 relations contain an adjustment factor that can be used to convert the equations to equivalent values for other stress drops. The predictions for 20 bars are shown on figure 10 and agree quite well with the observations. The AB06 relations are based on an empirically calibrated theoretical model of ground motions and are mainly applicable to larger earthquakes (**M** > 4). By contrast, Kaka and Atkinson (2005) developed an entirely empirical ground-motion prediction equation for use in Ontario ShakeMap applications that is applicable to small-to-moderate events (**M** < 5). The Kaka and Atkinson (2005) relation is also shown on figure 10. It matches the observational data from the Lively earthquake well. Thus the ground-motion amplitudes that were experienced in this event are consistent with our expectations for an event of this size, at distances from near-source to hundreds of kilometers.

## SUBSURFACE OBSERVATIONS

A unique aspect of this dataset is the availability of observations from below the surface, within very hard rock, only a few hundred meters from the hypocenter. Figure 11 shows a comparison of horizontal-component observations on the surface to those at depth, while figure 12 provides the same plot for the vertical components. Note that the acceleration spectra are nearly constant with frequency over the range available (3-40 Hz), as the corner frequency of the event is about 2.8 Hz (equation 6). In these plots, the surface observations (from SSNO and 11SNO, which are very near the epicenter, and from SUNO, 22 km away) are corrected to a hypocentral distance of 1 km, using equation 1. The average source spectrum at 1 km from all stations is also shown. The subsurface observations are at hypocentral distances of 0.4 to 1.3 km. They are *not* corrected to a distance of 1 km. The reason is that we are uncertain as to the appropriate distance correction in this case. The geometric spreading factor of 1.3 found by Atkinson (2004) applies to observations on the surface, from hypocentral distances of about 10 km and greater. It results from the travel of the waves through the crustal velocity structure. At close distances, below the surface, it would be reasonable to suppose that the theoretical geometric spreading factor of 1.0 for body waves might apply—although for the closest stations the near-field terms of equation 7 may also need to be considered. However, it can be observed from figures 11 and 12 that the closest underground station, DSNO (0.4 km), has lower amplitudes than the other underground stations, LSNO and 46SNO, both of which are at distances close to 1 km from the hypocenter. Thus any distance correction to normalize the records to a reference distance of 1 km would make little change to LSNO or 46SNO while reducing the amplitude at DSNO.

Overall, from the comparison of surface and subsurface records, we can conclude that 1) the distance-corrected surface records are in agreement with our expectations based on the apparent source spectrum obtained from distant stations; and 2) the subsurface records, which are at a distance close to the reference distance of 1 km for the surface records, have amplitudes that are lower than those at the surface by a factor in the range from 3 to 10. This may have important implications for ground motions experienced in underground waste repositories and may also explain the lack of significant damage to the Creighton mine workings.

## CONCLUSIONS

A mining-induced earthquake of MN 4.1 at a depth of 2.35 km, at Lively, Ontario, near Sudbury, offers a rare close-up view of a moderate event in an intraplate setting. The mainshock had a moment magnitude of 3.7 to 3.8 and a stress drop near 20 bars. The close proximity of deep underground sensors enabled recordings of numerous aftershocks and rarely seen near-field components of the seismic wavefield. Overall, ground motions were as expected, according to the prediction equations of Kaka and Atkinson (2005) and Atkinson and Boore (2006) for an event of this magnitude and stress drop, at epicentral distances ranging from 0 to hundreds of kilometers. Subsurface ground motions, recorded at distances from 0.4 to 1 km above the hypocenter, were a factor of 3 to 10 lower than the apparent source spectra.

An interesting finding is that when regional waveforms are used in a moment tensor inversion, the resulting best-fit depth (7 km) exceeds the known depth (2.35 km) of the event. However, a careful examination suggests that the degree of fit to the observed waveforms is not degraded significantly if the moment tensor solution is forced to a shallower depth. Regional waveform modeling is an alternative method for obtaining focal depths for moderate events. Regional waveform modeling provides a well-constrained estimate of focal depth for the Lively event in the range of 2.5 to 3 km, in good agreement with its known depth.

## Acknowledgments

We gratefully acknowledge the cooperation and invaluable assistance to our monitoring efforts given to us by the Sudbury Neutrino Observatory and Inco. We thank Eliza Richardson for her constructive review comments. Financial support for this study was provided by the Natural Sciences and Engineering Research Council of Canada and by Ontario Power Generation.