# Bulletin of the Seismological Society of America

## Abstract

We investigate the portability of a spectral amplitude method in determining moment tensor solutions for mine‐related seismicity recorded by deep mining networks to tectonic earthquakes recorded by surface aftershock networks. The original methodology inverts spectral amplitudes—with polarity attached—for direct *P*, *SV*, and *SH* waves recorded at deep‐seated geophones, so surface recordings must first be cleaned from free‐surface effects to recover the incident *P* and *S* waves. For precritical incidence, the correction is easily achieved by dividing the vertical component of the *P* and *SV* waveforms by the corresponding free‐surface reflection coefficients; for postcritical incidence, a more sophisticated correction that accounts for waveform distortion introduced by the coefficient’s phase shift is needed. Correction of *SH* components is achieved through division by a factor of 2. The proposed corrections are applied to 16 earthquakes recorded at local distances (<10 km) by an aftershock network deployed in the locality of São Caetano, Pernambuco, between 15 September and 23 December 2010. Spectral amplitudes are obtained through a time‐domain technique, and polarities assigned by comparing the corrected seismic pulses with a low‐pass‐filtered delta function. Although interference with the *S*‐to‐*P* critical reflection prevented the use of many *SV* amplitude measurements, inversion of the remaining spectral amplitudes allowed the recovery of deviatoric moment tensors for most of the events. Comparison with an independent first‐motion fault‐plane mechanism developed for the area shows consistency with our moment tensor solutions. Additionally, the ported methodology allows estimation of moment magnitudes for the selected events.

## Introduction

Although intraplate seismicity is known to characterize many tectonically stable portions of the continents, how and why earthquakes occur away from tectonic plate boundaries remains an open question. The traditionally accepted view holds that intraplate earthquakes happen on old, pre‐existing faults that are favorably oriented with respect to the regional stress field (Talwani, 1999); however, perturbations of the regional stress field from either flexural stresses (e.g., Assumpção and Sacek, 2013), lithospheric potential energy variations (e.g., Nielsen *et al.*, 2014), and/or local stress concentrators such as shallow plutons (Stevenson *et al.*, 2006) may also play a role. Stress concentrators were recently introduced in Talwani (2014) as discrete structures where pockets of elevated strain rate trigger the buildup of local stresses, resulting in the rotation of the regional stress field with wavelengths of tens to hundreds of kilometers. Characterizing intraplate stress fields in detail is key to improving our understanding of the underlying causes behind the occurrence of intraplate earthquakes.

The stress field around a fault can be derived from the analysis of the P and T axes from fault‐plane solutions (e.g., Jost and Herrmann, 1989) or, more accurately, from the inversion of the slip direction and the shear stresses on fault planes (Michael, 1984, 1987). The most stable estimates of fault‐plane solutions probably arise from analysis of seismic moment tensors from regional waveform modeling. Indeed, this has been the approach of choice for relatively large intraplate earthquakes. Examples include the 21 February 2008 *M*_{w} 6.0 event in Svalvard (Junek *et al.*, 2014), or the 21 May 1997 *M*_{w} 5.8 earthquake in Jabalpur (Saikia, 2006), which were well recorded at regional distances by sparse networks. This approach has also been applied to moderate intraplate events continuously monitored by regional networks in active seismic zones (e.g., Herrmann and Ammon, 1997) and/or captured by temporary deployments (e.g., Steffen *et al.*, 2012). For small intraplate events that are well recorded at local distances, on the other hand, development of fault‐plane solutions has often relied on the analysis of first‐motion polarities (e.g., Johnson *et al.*, 2014). An excellent example is northeastern Brazil, which is characterized by swarm‐like seismic activity with maximum magnitudes of 5.0 *m*_{b} (see e.g., Ferreira *et al.*, 1998). Aftershock monitoring networks are generally deployed in the zone of activity shortly after notice of felt earthquakes is given by the population, and composite fault‐plane solutions are then developed from first‐motion polarities constrained by aftershock fault‐plane geometry (e.g., Lima Neto *et al.*, 2013, 2014; Oliveira *et al.*, 2015).

Development of full moment tensor solutions for small earthquakes is, nonetheless, quite common in deep‐mine settings. Rather than first‐motion polarities, a widely utilized technique is the inversion of spectral amplitudes for the direct *P*, *SV*, and *SH* waves with polarities attached (Trifu *et al.*, 2000). The methodology assumes the propagating medium is an infinite whole space of uniform velocity and, therefore, that straight ray paths connect sources and receivers. As no free surface is present in deep‐mine environments, a simple rotation of the seismograms into the local ray‐coordinate system suffices to obtain the *P*, *SV*, and *SH* components of the direct *P* and *S* waves. Spectral amplitudes are then measured in the frequency domain from the amplitude spectra of the *P*, *SV*, and *SH* waveforms or in the time domain through the integrals of the squared displacements and velocities (Urbancic *et al.*, 1996). The methodology was successfully utilized, for instance, by Trifu and Shumila (2002) to determine full moment tensor solutions for 21 mine‐related events with moment magnitudes (*M*_{w}) between −1.2 and 0.0, recorded by an underground network of accelerometers in a deep metal mine in Canada, and by Julià *et al.* (2009) to develop full moment tensor solutions for 76 mine tremors with moment magnitudes between 0.5 and 2.6 for a deep gold mine in South Africa.

In this work, we investigate the portability of the spectral amplitude method of Trifu *et al.* (2000), from developing moment tensor solutions of mine‐related seismicity recorded by deep‐seated networks to constraining focal‐plane geometry for tectonic earthquakes recorded by local aftershock networks on the Earth’s surface. The main challenge is the removal of the free‐surface reflections and pulse distortion (at postcritical distances) from the surface recordings to recover the *P*, *SV*, and *SH* components of the incident *P* and *S* waves. We show that the recovery can be successfully achieved through estimation of the modulus and phase of the free‐surface reflection coefficients and the Hilbert transform of the original recordings (Choy and Richards, 1975), although interference with the *S*‐to‐*P* critical reflection may prevent a clean recovery of the *SV* component at near‐critical distances. We then apply the proposed corrections to 16 events recorded by a surface aftershock network in the locality of São Caetano, Pernambuco, and invert the recovered spectral amplitudes to obtain the deviatoric moment tensor for each individual event. A composite fault‐plane solution was developed by Lima Neto *et al.* (2013) from the polarity of the first *P*‐wave arrivals in São Caetano aftershock recordings, with fault‐plane geometry constraints. We find that the deviatoric moment tensor solutions for our best‐constrained events are consistent with the composite fault‐plane mechanism, and that the rejection of near‐critical *SV* amplitudes during inversion is key for the successful determination of the deviatoric moment tensor. In addition, we estimate moment magnitudes for the selected events from the scalar seismic moment, which are found to be in the 0.1<*M*_{w}<1.4 magnitude range.

## Free‐Surface Corrections

One of the main differences of recording seismic waves with a surface network, with respect to deep‐mine settings, is the presence of a free surface. First of all, reflected *P* and *S* waves are generated upon impingement of an incoming *P* or *S* wavefront, so seismometers located on the free surface record both incident and reflected waves (including *P*‐to‐*S* and *S*‐to‐*P* conversions). Second, for incident *SV* waves, an *S*‐to‐*P* critical reflection is generated for *SV* waves at critical incidence angles, which may interfere with the *SV* pulse at near‐critical distances. Moreover, the *SV*‐to‐*SV* reflection coefficient becomes complex at postcritical angles, introducing a constant phase shift in the reflected *SV* wave that distorts the postcritical signal from its original pulse shape. Free‐surface effects are demonstrated in Figure 1 through synthetic seismograms computed for a half‐space of uniform seismic velocity and density, and for epicentral distances between 2 and 10 km (source depth is 4 km). Figure 1 shows that the effect on *P* and *SH* waves is a simple scaling factor, with no pulse distortion with respect to the incident *P* or *SH* wave at any epicentral distance. The effects on *SV* waves, on the other hand, are more complex. On the vertical component, *SV* distortion is small and restricted to near‐critical distances, where interference with the *S*‐to‐*P* critical reflection occurs. On the radial component, on the other hand, *SV* distortion is severe and persists at postcritical distances. In this case, both interference with a large‐amplitude *S*‐to‐*P* critical reflection (near‐critical distances) and distortion from the constant phase introduced by the free‐surface reflection coefficient seem to be at play.

To recover the incident *P*, *SV*, and *SH* amplitudes from the surface recordings, we must first rotate the seismograms from the ZNE (geographical) system into the ZRT (vertical–radial–transverse) system. As the *SH* wave is polarized parallel to the surface and perpendicular to the plane of the ray, and no *S*‐to‐*P* conversions are expected from an incident *SH* wave (Fig. 1), the transverse component naturally isolates the *SH* pulse. Recovery of the incident *SH* amplitude is then as simple as dividing the transverse component by the *SH*‐to‐*SH* reflection coefficient. *P* and *SV* amplitudes, on the other hand, are projected along both the vertical and radial components, and both experience mode conversion upon reflection. Recovery of the incident *P* and *SV* amplitudes is, therefore, more intricate.

### Recovery of the *P* Amplitude

Both *P* and *SV* waves are generated upon reflection from an incident *P* wavefront, so the recorded *P*‐wave pulse is actually the superposition of the incident *P* wave with reflected *P* and *SV* waves. For an incident *P* wave with positive polarity in the direction of propagation (left inset in Fig. 2), it can be shown that the vertical and radial components of the recorded displacements are given by (1)(e.g., Aki and Richards, 2002) and(2)in which *α* is *P* velocity, *β* is *S* velocity, *p* is the ray parameter, and represent the vertical and radial components of the total surface amplitude, respectively, and *u*^{(P)} represents the amplitude of the incident *P* wave.

The coefficients multiplying the incident *P*‐wave amplitude are always real (see Fig. 2), so recovery of the incident *P*‐wave amplitude from the seismograms can be easily achieved by either dividing the vertical component with the vertical coefficient in equation (1) or by dividing the radial component with the radial coefficient in equation (2). Figure 3a illustrates the accuracy of such recovery with the vertical synthetic seismograms displayed in Figure 1. First, *P*‐wave pulses were isolated by cutting 0.5 s before and 0.5 s after the theoretical *P*‐wave arrival time. Then, the modulus of the free‐surface reflection coefficient was computed through equation (1), assuming the same velocity and density values utilized to compute the synthetic seismograms. Finally, the *P*‐wave pulses were divided by the computed modulus of the reflection coefficients. The comparison with the predicted *P*‐wave pulses, obtained after assuming a whole space with the same velocity and density values, is excellent. Recovery from the radial component of the synthetic seismograms (not shown) performed similarly. In practice, however, we always choose the vertical component to recover the incident *P*‐wave amplitude, as the radial component may be subject to uncertainties during the rotation into the ZRT system from inaccurate source locations.

### Recovery of the *SV* Amplitude

As for *P*‐wave incidence, both *P* and *SV* waves are generated upon reflection from an incident *SV* wavefront, so the recorded *SV* pulse is actually the superposition of the incident *SV* wave with reflected *P* and *SV* waves. For an incident *SV* wave with positive polarity toward the vertical (right inset in Fig. 2), the vertical and radial components of the recorded displacements are given by(3)(Aki and Richards, 2002) and (4)in which *α* is *P* velocity, *β* is *S* velocity, *p* is the ray parameter, and represent the vertical and radial components of the total surface amplitude, respectively, and *u*^{(S)} represents the amplitude of the incident *SV* wave.

For incident *SV* waves, there exists a critical angle *i*_{c} after which no reflected *P* wave is generated. Following Snell’s law, this critical angle is given by equation (5)which is obtained after setting the reflection angle for the converted *S*‐to‐*P* wave to 90°. The coefficients multiplying the incident *SV* amplitudes in equations (3) and (4) are real for precritical angles (*i*<*i*_{c}), but become complex for postcritical angles (*i*≥*i*_{c}; Fig. 2). Recovery of the precritical incident *SV* amplitude is therefore similar to the recovery of the incident *P* amplitude and can be achieved by either dividing the vertical component with the vertical coefficient in equation (3) or by dividing the radial component with the radial coefficient in equation (4). Recovery of the postcritical incident *SV* amplitude, on the other hand, requires more work. The phase of the coefficient introduces a constant phase into the reflected *SV* pulse at all frequencies. This means that the harmonic components making up the pulse are shifted differently in time and that the overall pulse shape will appear distorted with respect to the incident pulse (recall Fig. 1). Following Choy and Richards (1975), the distorted pulse can be written as (6)in which *f*^{R}(*t*) is the distorted pulse, *f*(*t*) is the original (undistorted) pulse, *H*[ ] denotes the Hilbert transform, and *R* and are the amplitude and phase of the coefficient, respectively. Choy and Richards (1975) also show that equation (6) can be inverted to (7)effectively providing the means to recover the incident, undistorted *SV* pulse from the surface, distorted recording.

Figure 3b illustrates the accuracy of the recovery of the *SV* incident amplitude at both pre‐ and postcritical incidence with synthetic seismograms. Similar to *P*‐wave recovery, *SV*‐wave pulses were first isolated by cutting 0.5 s before and after the theoretical *S*‐wave arrival time. Then, the modulus and phase of the free‐surface reflection coefficient was computed through equation (3), assuming the same velocity and density values utilized to compute the synthetic seismograms. Finally, the Hilbert transform of the *SV*‐wave pulses were computed and the incident *SV* pulse recovered through application of equation (7) to the synthetic seismograms. The comparison with the predicted *SV*‐wave pulses, obtained after assuming a whole space with the same velocity and density values, is excellent at both pre‐ and postcritical incidence. At near‐critical incidence (4‐km epicentral distance), the recovery is quite good, observing only a small phase shift due to interference with the *S*‐to‐*P* critical reflection in the free‐surface seismogram. Recovery from the radial component of the synthetic seismograms (not shown) displays lower performance at near‐critical distances, due to both stronger amplitudes of the *SP* critical reflection and stronger dispersion (recall Fig. 1). As with the recovery of the *P*‐wave amplitude, we always choose the vertical component to recover incident *SV*‐wave amplitudes.

### Recovery of the *SH* Amplitude

As mentioned before, only *SH* waves are generated upon reflection from an incident *SH* wavefront. The reflection coefficient is a mere factor of 2 (Aki and Richards, 2002), so recovery of the incident *SH* waveform is as simple as dividing the transverse component around the *S* wave by a factor of 2. The recovery is illustrated with synthetic seismograms in Figure 3c.

## Porting the Spectral Amplitude Method

The portability of the spectral amplitude method of Trifu *et al.* (2000) to surface aftershock networks is tested with local seismograms recorded near the locality of São Caetano, in the Brazilian state of Pernambuco. Seismic activity commenced on 11 March 2010 with an event of magnitude 2.7 *m*_{R} that was felt by the population, and continued for several months after the occurrence of the first event. A monitoring network was deployed in the area by the Laboratório Sismológico of the Universidade Federal do Rio Grande do Norte (LabSis/UFRN) between 15 September and 23 December. The network consisted of seven, three‐component, short‐period stations recording continuously at 250 samples per second. Each station was equipped with a Sercel L4A‐3D sensor, with flat response in velocity above 2 Hz, a RefTek‐130 digitizer, and a Global Positioning System receiver for timekeeping. Nominal sensor sensitivity was reported as 276.8 V/m/s by the manufacturer, although it was lowered to a value of 176.0 V/m/s in our analysis to account for the effect of the shunt resistor (9016 Ω). Digitizers were set to operate on high‐gain setting, with a nominal high‐gain value of 85 nV/count reported by the manufacturer. During this period, a total of 21 seismic events were identified and located through the aftershock network (Lima Neto *et al.*, 2013). Although magnitude was not calculated for that sequence in Lima Neto *et al.* (2013), the authors do report that seismic events in the area are typically in the 0.5<*m*_{R}<3.7 range. Source locations and origin times are listed in Table 1, whereas epicentral locations and network geometry are displayed in Figure 4.

### Recovering *P*, *SV*, and *SH* Pulse Shapes

Recovery of the incident *P*, *SV*, and *SH* pulse shapes on the recordings of the 2010 São Caetano aftershocks was achieved through the free‐surface corrections outlined in the Free‐Surface Corrections section. First, the original seismograms were rotated from the ZNE system into the great‐circle path to obtain the ZRT components. The same event locations and origin times as determined by Lima Neto *et al.* (2013)—with reported confidence bounds of less than 0.5 km in both epicentral distance and depth—were assumed (see Table 1). The rotated waveforms were then demeaned, detrended, and tapered, with the instrument response then removed through deconvolution to obtain the displacement seismograms in the 2–100‐Hz frequency band. Second, *P* and *S* travel‐time arrivals were manually picked and the *P* and *S* waveforms windowed from the displacement seismograms. We used windows between −0.5 s and one‐tenth of the *S*–*P* travel time around the *P*‐wave arrival (vertical component), and between −0.5 s and twice the *S*–*P* travel time around the *S*‐wave arrival (vertical and transverse components). Finally, *P* and *SV* pulse shapes for the incident *P* and *S* waves were recovered through application of equation (7) to the vertical component of the windowed seismograms. The *SH* pulse shape was obtained after dividing the windowed transverse component by a factor of 2.

The coefficients and phases in equation (7) were calculated from equations (1) and (3), after assuming a *P* velocity of 5.9 km/s, an *S* velocity of 3.45 km/s, and a density of 2658 kg/m^{3}. *P*‐ and *S*‐velocity values were taken from the location model of Lima Neto *et al.* (2013), whereas density was determined with an empirical Birch‐type relationship (Berteussen, 1977). The ray parameter for the incident *P* and *S* waves was obtained after assuming a straight ray path between source and receiver, and calculating the incidence angle according to geometry. Examples of rotated displacement waveforms and recovered *P*‐, *SV*‐, and *SH*‐wave pulses are displayed in Figure 5 for the recordings of event 08 (Table 1) at station SOPS.

### Measuring Spectral Amplitudes

Spectral amplitudes were measured on the recovered *P*, *SV*, and *SH* waveforms in the time domain through the squared integrals of displacement and velocity. Following Urbancic *et al.* (1996), the spectral amplitude is given by (8)in which *u* is the spectral amplitude, *S*_{D2}=∫*D*^{2}(*t*)*dt* with *D* being displacement, and *S*_{V2}=∫*V*^{2}(*t*)*dt* with *V* being velocity. Velocity seismograms were obtained after differentiation of the instrument‐corrected displacement seismograms, and windowed around the *P*‐ and *S*‐wave arrival times as described before for the displacement seismograms. Before calculating the integrals, the seismograms were further filtered between 10 and 30 Hz. To assess the accuracy of the measurements, we compared the measured spectral amplitudes with the amplitude spectra of the recovered *P*‐, *SV*‐, and *SH*‐wave pulses. This is illustrated in Figure 5b, where the amplitude spectra for the recovered incident waves for event 08 at station SOPS are superimposed to the time‐domain estimates from equation (8). Figure 5b demonstrates that observed and measured spectral levels match well within the 10–30-Hz frequency band. We experimented with a wide range of frequency bandwidths and found that the selected frequency band gave the most satisfying agreement between time‐ and frequency-domain estimates.

The method of Trifu *et al.* (2000) is based on the inversion of *P*, *SV*, and *SH* spectral amplitudes with polarity attached, so polarities need to be assigned to all spectral amplitudes. To achieve this goal, we followed the approach of Julià *et al.* (2009). First, we obtained a positive, unit‐amplitude pulse by convolving a delta function with the nominal instrument response. Then, both the unit‐amplitude signal and the time‐domain *P*, *SV*, or *SH* pulses were low‐pass filtered below 10 Hz to obtain a simple pulse shape. Visual correlation or anticorrelation with the low‐pass‐filtered delta function would then reveal a positive or negative polarity, respectively. The process is illustrated in Figure 5c, where amplitude‐scaled pulses are superimposed to the recovered *P*, *SV*, and *SH* amplitudes in the time domain. The comparison reveals the *P* and *SH* components have positive polarity, whereas the polarity for the *SV* component is negative.

Instrument responses for the seismic sensor and digitizer were developed from the natural frequency, damping, sensitivity, and gain reported by the manufacturers. As mentioned earlier, the nominal sensor sensitivity was corrected to account for the effect of the shunt resistor. Additional corrections included flipping the polarity of the waveform recordings, which turns out to be reversed in the Sercel L4A‐3D seismometers (J. M. Ferreira, personal comm., 2014). Finally, we also detected that stations within the aftershock network might have operated at different gains. The issue was easy to correct, as RefTek‐130 digitizers may only be set at either high‐ or low‐gain values (which differ by a factor of 32). Plotting the seismograms on the same vertical scale allows for a quick identification of the waveforms that operated at a low‐gain setting, and multiplying those waveforms by 32 equalizes all waveforms to the same high‐gain setting.

### Inverting Spectral Amplitudes

Once spectral amplitudes have been measured on the free‐surface‐corrected seismograms and polarities have been attached, inversion for the seismic moment tensor closely follows the scheme of Trifu *et al.* (2000). As amplitudes are associated with the flat, low‐frequency portion of the seismic spectrum, it is assumed that ray propagation occurs along straight ray paths within a half‐space of uniform velocity and density. We assumed the same parameters as for the free‐surface corrections (*V*_{P}=5.9 km/s, *V*_{S}=3.45 km/s, *ρ*=2658 kg/m^{3}). Under such assumptions, the forward problem can be formulated as (9)in which is the vector of *P*, *SV*, or *SH* spectral amplitudes, *c*^{(i)}=1/(4*pv*^{3}*R*) with *v* either the *P*‐ or *S*‐wave velocity, *R* is the hypocentral distance, *ρ* is the density, **M** is the seismic moment tensor, and **F**^{(i)} is an excitation matrix that depends on the takeoff angle and event azimuth and combines with **M** to express the radiation pattern for *P*, *SV*, or *SH* waves. The colon in equation (9) represents the double dot product between **F** and **M**. Equation (9) defines a system of linear equations that, considering the symmetry of the moment tensor, can be expressed as (10)in which *n* is the number of data points (spectral amplitudes), *m*_{1}=*M*_{11}, *m*_{2}=*M*_{12}=*M*_{21}, *M*_{3}=*M*_{22}, *m*_{4}=*M*_{13}=*M*_{31}, *m*_{5}=*M*_{23}=*M*_{32}, and *m*_{6}=*M*_{33} and in which *f*_{1}=*F*_{11}, *f*_{2}=*F*_{12}+*F*_{21}, *f*_{3}=*F*_{22}, *f*_{4}=*F*_{13}+*F*_{31}, *f*_{5}=*F*_{23}+*F*_{32}, and *f*_{6}=*F*_{33}. The system is inverted with a zero‐trace constraint *m*_{1}+*m*_{3}+*m*_{6}=0 (deviatoric moment tensor) after a singular value decomposition (SVD) of the matrix (see e.g., Menke, 1984). As we need all six components of the moment tensor, no truncation of the SVD spectrum is applied, which is equivalent to minimizing the difference between observed and predicted spectral amplitudes in a purely least‐squares sense. More details on the inversion methodology can be found in Trifu *et al.* (2000).

As equation (10) represents a linear system of equations, confidence bounds can be easily assessed after propagating uncertainties in spectral amplitudes through (11)(see e.g., Taylor, 1997), in which (*σ*_{m})_{i} is the uncertainty of the *i*th moment tensor component, (*σ*_{d})_{k} is the uncertainty of the *k*th spectral amplitude measurement, and *a*_{ik} is the inverse matrix associated with the solution of equation (10). Uncertainties for the spectral amplitudes (*σ*_{d})_{k} can be estimated as a given percentage of the measured amplitude values (Fletcher and McGarr, 2005).

To illustrate the inversion method, Figure 6 displays the inversion results for event 08 (see Table 1). The inverted moment tensor has the following components: *m*_{xx}=0.32±0.03×10^{11} N·m, *m*_{xy}=0.06±0.07×10^{11} N·m, *m*_{xz}=−1.38±0.11×10^{11} N·m, *m*_{yy}=−0.13±0.02×10^{11} N·m, *m*_{yz}=0.08±0.04×10^{11} N·m, *m*_{zz}=−0.19±0.03×10^{11} N·m (*x* is north, *y* is east, and *z* is down). Confidence bounds were obtained after assuming that uncertainties in spectral amplitude are 10% of the measured amplitude. The total scalar moment is *M*_{0}=1.47±0.28×10^{11} N·m. The figure reveals that the moment tensor solution predicts the polarity and amplitude of the spectral amplitudes utilized in the inversion quite satisfactorily. The correlation coefficient of observed versus predicted amplitudes in the frequency domain is close to 100%, and the comparison between predicted and observed amplitudes in the time domain is visually satisfying. The only exceptions are the *P*‐wave amplitudes for stations SOSB and SOVJ, which are clearly underpredicted by the moment tensor solution, one of them (SOSB) with opposite polarity. Also, as expected, the moment tensor solution seldom predicts amplitudes that were not utilized during the inversion.

More insight can be gained from the eigenvalues (*σ*_{i}) and eigenvectors (**e**_{i}) associated with the moment tensor. Realize that the eigenvalues reveal the tensional (positive) or compressional (negative) character of the deviatoric strains at the source, whereas the eigenvectors provide their direction. For this moment tensor, we obtain that *σ*_{1}=(1.47±0.14)×10^{11} N·m, *σ*_{2}=(−0.12±0.01)×10^{11} N·m, and *σ*_{3}=(−1.34±0.07)×10^{11} N·m and that **e**_{1}=(0.769±0.005,−0.007±0.052,0.640±0.007), **e**_{2}=(−0.058±0.221,−0.997±0.014,0.059±0.021), **e**_{3}=(−0.637±0.009,0.083±0.006,0.766±0.008). Further decomposition of the deviatoric moment tensor into its double‐couple (DC) and compensated linear vector dipole contributions (Knopoff and Randall, 1970), reveals a predominantly DC component of up to 83%. Assuming that local strain eigenvectors from our moment tensor solutions are close to the local stress eigenvectors, we infer that the intermediate stress is horizontal and pointing in the east–west direction, whereas the maximum (positive) and minimum (negative) deviatoric stresses lie at about 45° in the north‐vertical plane. Accordingly, the best DC for this moment tensor solution displays near‐vertical and near‐horizontal focal and auxiliary planes. The strike , dip *δ*, and rake *λ* for the best‐DC planes are , *δ*=85°, *λ*=−93°, and , *δ*=6°, *λ*=−57° which, considering uncertainties reported for the moment tensor components, fall within ranges of 86°±2°, 85°±1°, and −94°±4° and 294°±29°, 7°±3°, and −62°±31° for the fault angles of the fault and auxiliary planes, respectively. The focal mechanism plots for the deviatoric and best‐DC solutions are displayed in Figure 6c. A near‐vertical fault‐plane orientation striking in the east–west direction was also obtained from the composite first‐motion polarity analysis of Lima Neto *et al.* (2013).

### Deviatoric Moment Tensor Solutions

The procedure outlined in the Measuring Spectral Amplitudes and Inverting Spectral Amplitudes sections was applied to the 21 events listed in Table 1 to determine individual deviatoric moment tensor solutions for each of them. We report results on 16 of them, corresponding to those for which enough (minimum of 5) acceptable amplitude measurements could be obtained. By “acceptable”, we mean that polarities could be visually attached after low‐pass filtering below 10 Hz. The results are summarized in Table 2, which lists individual moment tensor components, correlation coefficient between observed and predicted spectral amplitudes, normalized root mean square (nrms), DC percentage (Knopoff and Randall, 1970), and predicted polarity ratios (number of correctly predicted polarities over number of observations) for *P*, *SV*, and *SH* amplitudes. We define the nrms as (12)in which *u*^{(obs)} means observed amplitude and *u*^{(pred)} is predicted amplitude. The nrms will vary between 0 (perfect fit) and ∞ (bad fit). Very few *SV* amplitudes were actually utilized during the development of deviatoric moment tensor solutions. Most of the stations were located at near‐critical distances from the source, so that *SV* amplitudes could have interference with the *SP* critical reflection (recall Fig. 1). Only *SV* amplitudes for stations well above or well under the critical distance were considered, considerably reducing the number of available observations for this amplitude type.

Given the expected tectonic origin of the aftershocks under analysis (Lima Neto *et al.*, 2013), one anticipates solutions with large DC percentages (defined here as 80% or more). A closer inspection of Table 2 reveals that, in general, such solutions are restricted to those constrained by a large number of amplitudes (events 01, 04, and 16). Moreover, for those solutions, *SH* polarities are always predicted correctly, whereas the few incorrectly predicted polarities are mostly of the *P*‐wave type. We suspect this is due to the smaller size of *P*‐wave amplitudes, which gives them less weight during minimization of the root mean square (rms) norm. Correlation coefficients are above 70% for all inversions, and above 80% for those with a large DC percentage; nrms values are always under 0.6, but we do not see particularly small values for events with large DC percentages. Realize that recovery of the focal mechanism depends mainly on a good azimuthal coverage of the radiation pattern. Large nrms values just indicate good prediction of the measured amplitudes, which is easier to achieve when few measurements are available (so that azimuthal coverage is poor). Two notable exceptions to the general pattern outlined above include event 03, which is constrained by 12 amplitudes, and yet, has a relatively low DC percentage of 65%, and event 18, which displays a DC percentage of 95% while being constrained by just six amplitudes. Event 03 has somewhat smaller spectral amplitudes, on the order of 10^{−11} m·s (i.e., when compared to amplitudes of the order of 10^{−10}–10^{−9} m·s for large DC events), which makes them more likely to be contaminated by noise. Event 18 is constrained by just *SV* and *SH* amplitudes from up to four different stations, which might or might not be providing sufficient azimuthal coverage. The feasibility tests described in the Feasibility Analysis section will bring more light into this issue.

## Discussion

### Feasibility Analysis

The deviatoric moment tensor solutions reported in Table 2 were obtained through quite an inhomogeneous set of spectral amplitude measurements. Moreover, although there seems to be a quite direct relationship between number of spectral amplitudes utilized in the inversion and DC percentage of the moment tensor solutions, the relationship is far from infallible. To investigate this issue further, we conducted a couple of feasibility tests on the 16 moment tensor inversions summarized in Table 2. The first test consisted of creating a synthetic set of spectral amplitudes—contaminated up to a given percentage of random noise—for a known moment tensor solution and with the same distribution of stations and amplitude types as observed for a given event. As the known moment tensor solution, we chose that from the composite fault‐plane solution of Lima Neto *et al.* (2013), with an arbitrary scalar seismic moment of 10^{14} N·m. The synthetic amplitudes were then inverted and the moment tensor solution compared to the true moment tensor. The result of such a test is displayed in Figure 7a, where the nrms (as in equation 12) between inverted and true moment tensor components is plotted against event number for noise levels of 10% and 20%. It is apparent that events 15, 18, 20, and 21 are the most unstable, whereas events 01 and 08 are among the most stable. Recall that event 18 had a high DC percentage while being constrained by just six amplitudes (Table 2), whereas events 01 and 08 had high DC percentages and were constrained by up to 11 amplitude measurements.

The second test consisted of replicating the same exercise, but using the inverted moment tensor solution in Table 2 as the true solution for each event, instead of the composite focal mechanism of Lima Neto *et al.* (2013). The results of this second test are displayed in Figure 7b. Interestingly, events 15, 20, and 21 are now quite stable, whereas event 18 remains unstable. Events 01 and 08, nonetheless, remain stable under this new test. We suspect the change in stability reflects the variability in azimuthal coverage of the radiation pattern associated with each true focal mechanism. Realize that both compressional and dilational quadrants must be sampled to robustly constrain a focal mechanism and that, given a fixed set of stations, sampling will vary with focal mechanism.

### Comparison to Composite Fault‐Plane Solution

A comparison of our most robust deviatoric moment tensor solutions (see the Feasibility Analysis section) with the composite fault‐plane solution of Lima Neto *et al.* (2013) is given in Table 3. The table lists the strike, slip, and rake of the fault and auxiliary planes—along with the ranges associated with uncertainty in the inverted moment tensors—for the best DC component of our deviatoric moment tensor solutions, as well as the strike, dip, and rake angles for the fault and auxiliary planes of the composite fault‐plane solution. To aid in the comparison, the *P*‐wave radiation patterns (or focal mechanism plot) for each of the listed fault planes are superimposed in the geological map presented in Figure 4. The composite fault‐plane solution suggests normal faulting dominated the São Caetano area during the seismic sequence under study; except for event 13, which indicates reverse faulting, and event 16, which indicates strike‐slip faulting, our deviatoric moment tensor solutions also display predominant normal faulting. If we focus on the large DC‐percentage events (01, 04, 08, and 16), which are likely to be the most robustly constrained, we see the majority of our moment tensor solutions are consistent with fault planes in the east–west and northeast–southwest directions with a mixture of normal and strike‐slip mechanisms. This distribution of focal mechanisms is consistent with known orientation of shear zones in northeast Brazil and with the known state of stress in the region (e.g., Bezerra *et al.*, 2011); a few focal mechanisms, nonetheless, differ substantially from that of Lima Neto *et al.* (2013).

To assess whether the variety of focal mechanisms exhibited in Table 3 is reliable, we computed synthetic seismograms for some of our large DC deviatoric moment tensor solutions and the composite fault‐plane solution of Lima Neto *et al.* (2013). Let us consider first event 08, for which our fault mechanism predicts a radiation pattern that is nearly opposite to that predicted by the composite first‐motion solution (see Table 3). Comparison between observed and synthetic seismograms (Fig. 8a) reveals that, in general, our fault‐plane solution predicts size and polarity for *P*, *SV*, and *SH* waveforms better than the composite solution. Moreover, as expected from the *P*‐wave radiation patterns in Table 3, the composite fault‐plane solution has reversed *P*‐wave polarities when compared with our moment tensor predictions, or to observations. Better agreement in size and polarity is also observed when comparing observed and synthetic waveforms for events 01 and 04 (not shown), although it is less apparent for event 04 due to the similarity in the predicted radiation patterns. Regarding event 16, for which our spectral amplitude inversion predicts a strike‐slip mechanism, the comparison is displayed in Figure 8b. The *SV* and *SH* radiation pattern is clearly better predicted by our inverted focal mechanism. Prediction of the *P*‐wave radiation pattern is harder to assess due to the small *P*‐wave amplitudes at low frequencies; although polarities seem to be predicted by both mechanisms with similar success, our fault‐plane solution does a better job at predicting amplitude size.

The comparison with synthetic seismograms displayed in Figure 8 demonstrates the composite, fault‐plane solutions cannot satisfactorily explain the radiation patterns associated with events 01, 04, 08, and 16. Incidentally, none of those events were utilized to develop the composite solution of Lima Neto *et al.* (2013), as only events with very tight values in rms misfit, horizontal error location, and azimuthal gap were considered in that work. Nonetheless, Lima Neto *et al.* (2013) showed that those events lie along a well‐defined plane, suggesting that they indeed happened within the same fault. The composite fault‐plane solution may just be one more focal mechanism for the region, and may not necessarily be representative for the entire seismic zone. The variation in fault‐plane orientation and slip suggested in Table 3 for the São Caetano sequence of 2010, therefore, demonstrates that the seismic zone is not as simple as a single fault plane, but is a complex fault system where a variety of fault planes might be interacting.

### Dependency on Earth Model Parameters

The moment tensor solutions reported in Table 2 were developed after matching theoretical spectral amplitudes with observed spectral amplitudes measured at a number of seismic stations. Theoretical amplitudes were computed assuming a half‐space with *P* velocity of 5.9 km/s, *S* velocity of 3.55 km/s (*V*_{P}/*V*_{S} ratio of 1.71), and density of 2658 kg/m^{3}, consistent with the study of Lima Neto *et al.* (2013). These Earth model parameters, however, are also subject to uncertainty. In this section, we assess the significance of the Earth model by exploring the dependence of the inverted moment tensor components on Earth model parameters.

The Earth model is used at two different stages during the processing and inversion of spectral amplitudes: first, during the correction of the dataset for free‐surface effects; and second, during the inversion of the spectral amplitudes. To assess the effect on the correction of the free‐surface effects, we computed free‐surface coefficients (equations 1–4) for a range of velocity models, with *P* velocity ranging between 5.0 and 6.5 km/s and *V*_{P}/*V*_{S} ratios ranging between 1.65 and 1.75. Density varied with *V*_{P} according to the empirical relationship of Berteussen (1977). The resulting coefficients did not change significantly with *V*_{P}, and variations with *V*_{P}/*V*_{S} ratio were within ±0.1 of the central value. Variations were largest near grazing incidence for *P* waves, and largest near‐critical incidence for incident *SV* waves. Also, for a given incidence angle, variation was larger for the radial coefficients (equations 2 and 4) than for the vertical coefficients (equations 1 and 3). Away from grazing and near‐critical incidence, recovered *P* and *SV* amplitudes were seldom sensitive to variations in the free‐surfaces coefficients.

To assess the sensitivity of the inverted moment tensor components to variations in the Earth model parameters, we inverted the spectral amplitudes for the events listed in Table 2 with the same range of Earth models mentioned above. Changes in individual moment tensor components could deviate significantly from the original value, but that seemed to affect the overall radiation pattern just minimally; changes seem to be mostly absorbed by the scalar seismic moment. This behavior is illustrated in Figure 9, where the *P*‐wave radiation patterns for the full‐deviatoric and best‐DC moment tensor solutions for the suite of inversions corresponding to event 08 are displayed. The largest change in the moment tensor is observed for the *xz* component, which varies between 0.71×10^{11} N·m (*V*_{P}=5.0 km/s, *V*_{P}/*V*_{S}=1.75) and 2.16×10^{11} N·m (*V*_{P}=6.5 km/s, *V*_{P}/*V*_{S}=1.65). This translates into a progressive increase of the scalar moment tensor from 0.74×10^{11} N·m to 2.34×10^{11} N·m. Yet, fault angles offer moderate variations of less than 10° with no perceptible changes in the orientation of the fault planes.

### Moment Magnitudes

In addition to fault‐plane orientation and direction of slip, the spectral magnitude method allows estimation of moment magnitude from the inverted moment tensors. Moment magnitudes (*M*_{w}) are determined from the scalar seismic moment according to (13)(Hanks and Kanamori, 1979), in which *M*_{0} is the scalar moment tensor in dyn·cm. Using the scalar moments from our deviatoric moment tensor solutions (Table 3), we obtain moment magnitudes in the 0.1<*M*_{w}<1.4 range, which contrasts with the 0.5<*m*_{R}<3.3 range reported in Lima Neto *et al.* (2013) for seismic activity in the region. The *m*_{R} scale was developed by Assumpção (1983) from maximum *P*‐wavetrain amplitudes recorded in the 1–10 Hz frequency band by seismographs in the 200–1500 km epicentral distance, and was calibrated with teleseismic events in the 4<*m*_{b}<5 magnitude range. It was postulated to be a good proxy of the *m*_{b} scale down to magnitude 2. Our magnitude estimates, on the other hand, are based on moment tensor solutions that predict seismic amplitudes, and provide a better physical base for assessing the magnitude range for seismic activity in the São Caetano area.

Finally, we note that our spectral amplitudes were not corrected for attenuation, which means our reported magnitude estimates might be somewhat underestimated. Assuming quality factors of *Q*_{P}=600 and *Q*_{S}=350, as observed in other Precambrian regions worldwide, we estimate an amplitude reduction at 10 Hz and 5 km epicentral distance of just 5% for *P* waves and 12% for *S* waves. Indeed, the amplitude spectra displayed in Figure 5 show a quite flat plateau below the corner frequency, suggesting minimal amplitude loss by attenuation. We must then conclude that our inferred magnitude range of 0.1<*M*_{w}<1.4 is accurate.

## Conclusions

Summarizing, we ported the spectral amplitude method of Trifu *et al.* (2000) for mine‐related events recorded in deep mines to tectonic events recorded in surface aftershock networks. Portability relies on recovery of the incident *P*, *SV*, and *SH* amplitudes after correction of the free‐surface effects using *a priori* reflection coefficients, and on avoiding *SV* amplitude measurements at near‐critical distances. The ported methodology has been applied to surface recordings of a sequence of events monitored by a local aftershock network in the locality of São Caetano, Pernambuco. The majority of the inverted (deviatoric) moment tensor solutions has relatively large DC contents (above 60%), and associated fault‐plane mechanisms suggest stresses have been released through both normal and strike‐slip faulting. Scalar seismic moments revealed moment magnitudes in the 0.1–1.4 *M*_{w} range for the seismic events.

## Data and Resources

Plots and basic data processing were done with the Seismic Analysis Code, v.101.4 (http://ds.iris.edu/files/sac-manual/, last accessed May 2016). Maps were made using the Generic Mapping Tools, v.4.5.8 (http://gmt.soest.hawaii.edu, last accessed May 2016; Wessel and Smith, 1998). Radiation patterns and fault angles were obtained through the moment tensor plotting and decomposition (MoPaD) tool of Krieger and Heimann (2012). Moment tensor inversion codes were developed by J. Julià.

## Acknowledgments

We thank Karen Assatourians, an anonymous reviewer, and Associate Editor Cezar Trifu for the time and effort dedicated to carefully reading the original article. Their many comments and suggestions greatly improved the presentation and clarity of this article. This work was supported by the Brazilian Centro Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Grant Number 456791/2013‐2). Data were collected under the Instituto Nacional de Ciência e Tecnologia em Estudos Tectônicos (INCT‐ET of the CNPq, Grant Number 57.3713/2008‐01). S. L. E. F. d. was supported by a two‐year scholarship from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) to complete his M.Sc. degree at Universidade Federal do Rio Grande do Norte (UFRN). J. J. thanks CNPq for his research fellowship (CNPq, Process Number 304421/2015‐4).

- Manuscript received 24 May 2016.