S-wave attenuation in northeastern Sonora, Mexico, near the faults that ruptured during the earthquake of 3 May 1887 Mw 7.5

We used a new data set of relocated earthquakes recorded by the Seismic Network of Northeastern Sonora, Mexico (RESNES) to characterize the attenuation of S-waves in the fault zone of the 1887 Sonora earthquake (Mw 7.5). We determined spectral attenuation functions for hypocentral distances (r) between 10 and 140 km using a nonparametric approach and found that in this fault zone the spectral amplitudes decay slower with distance at low frequencies (f < 4 Hz) compared to those reported in previous studies in the region using more distant recordings. The attenuation functions obtained for 23 frequencies (0.4 ≤ f ≤ 63.1 Hz) permit us estimating the average quality factor QS = (141 ± 1.1 )f(0.74 ± 0.04) and a geometrical spreading term G(r) = 1/r0.21. The values of Q estimated for S-wave paths traveling along the fault system that rupture during the 1887 event, in the north–south direction, are considerably lower than the average Q estimated using source-station paths from multiple stations and directions. These results indicate that near the fault zone S waves attenuate considerably more than at regional scale, particularly at low frequencies. This may be the result of strong scattering near the faults due to the fractured upper crust and higher intrinsic attenuation due to stress concentration near the faults. Electronic supplementary material The online version of this article (doi:10.1186/2193-1801-3-747) contains supplementary material, which is available to authorized users.


Introduction
On May 3, 1887, a mayor seismic event (M w = 7.5) took place on northeastern Sonora (Mexico), destroying the town on Bavispe and its surroundings (Aguilera, 1888). This event, as well as most northern Mexico, is located south of the Basin and Range province (Suter and Contreras, 2002), and generated the longest recorded normal fault surface rupture (101.8 km) in historic time (Suter, 2006).
The Red Sísmica del Noreste de Sonora (RESNES) was installed in 2002 to monitor the seismicity related with the Basin and Range normal faults of northeastern Sonora (Castro et al., 2002;Romero et al., 2004). The seismic instruments of the RESNES array consist of Kinemetrics digital recorders (model K2) with internal Episensors that record the three components of ground acceleration, sampled at 200 samples per second, and an external shortperiod seismometer (model L4C) to record the vertical ground velocity. All stations are autonomous and have a GPS built-in timing system (Castro et al., 2010).
Since the occurrence of the major event in 1887, several studies have been made, including contemporary field studies (Goodfellow, 1888;Aguilera, 1888), studies of intensity and attenuation (DuBois and Smith, 1980;Sbar and DuBois, 1984;Bakun, 2006;Castro et al., 2008;Castro et al., 2009), regional seismotectonics (Suter and Contreras, 2002), geomorphology (Bull and Pearthree, 1988;Pearthree et al., 1990) and microseimicity (Natali and Sbar 1982). Condori (2006) and Castro et al. (2008) studied the spectral amplitude decay of body waves with hypocentral distance and proposed local and regional attenuation curves to characterize the attenuation near the fault zone and at distances beyond 100 km from the center of the network.
An important element to evaluate seismic hazard is the attenuation, particularly to estimate the intensity of ground-motion at different epicentral distances. Thus, the quality factor Q and the geometrical spreading are key parameters for ground-motion predictions and for seismic hazard analysis. We use in this article a new and more complete data set to study in more detail the seismic attenuation near the fault zone. This new data set is composed by earthquakes relocated by Castro et al. (2010) near the faults that rupture during the 1887 event. We also compare results from previous studies of spectral attenuation in the Sonora region with our new results and with results from other studies within the Basin and Range province.

Data
We selected 50 events, 398 horizontal (north-south and east-west components) and 199 vertical acceleration records recorded by RESNES from 2003 to 2007 that were located near the fault zone area. Figure 1 shows the location of the earthquakes used (black circles) and the main faults of the region. The fault segments shown with thicker lines are the segments that ruptured during the 1887 events, the Pitáycachi (P), Teras (T) and Otates (O) faults. The stations of the RESNES array (triangles) are distributed around these faults. These events were relocated by Castro et al. (2010), have local magnitudes between 0.5 and 3.5, hypocentral distances between 10 and 140 km and focal depths of less than 40 km. Table 1 lists the source coordinates and magnitudes of the events used and Table 2 the station coordinates. Figure 2a and b show the distribution of the magnitudes with epicentral distance and the histogram of the number of records per station, respectively. Most events have magnitudes between 1 and 2 ( Figure 2) and were recorded by stations NAC and OAX, located in the center of the array.
The records were baseline corrected, subtracting the mean, to remove long-period biases. We chose time windows starting before the first S-wave arrival and containing 80% of the S-wave energy to calculate the Fourier acceleration spectra of the north-south (NS) and east-west (EW) Figure 1 Tectonic map showing location of stations and earthquakes used. Black circles represent the events used for the attenuation study and triangles the recording stations. The stars represent centers of population. P (Pitáycachi), T (Teras), and O (Otates) are the faults that ruptured in the 1887 earthquake. Modified from Castro et al. (2010). components. We verify visually that the windows contain the strong-ground motions including the peak acceleration. The first and last 5% of the window are cosine tapered and the spectral amplitudes smoothed averaging within a variable frequency band of ±25% of 23 predefined central frequencies between 0.4 and 63.1 Hz. The spectral amplitudes were smoothed using a variable frequency band of +/− 25% over the 23 predefined central frequencies. For each spectral record, we selected the signal frequency band above noise level by visual inspection of the spectrum. Figure 3 shows examples of signal and noise spectra from stations ELO and MOR located on bedrock and conglomerates, respectively. This figure compares noise and signal spectral levels recorded at both horizontal components from events recorded at epicentral distances less than 20 km (left frames) and at 138.9 km (right frames). The signal spectra (solid lines) correspond to the 80% of the energy of the Swaves and the noise spectra (dashed lines) to 6-seconds time windows starting before the first P-wave arrivals. Examples of the acceleration spectra calculated for two events with magnitudes M L 3.5 and 1.9 are also shown in Figures 4 and 5, respectively. The signal to noise ratio for these stations is above one in the frequency band 1.0 -63.1 Hz for local events but increases for events with lager magnitude (M > 1.7). For regional events the signal to noise ratio of

Attenuation functions using a nonparametric method
We explore the dependence of the spectral amplitudes with hypocentral distance by considering that for a fixed frequency, the spectral amplitudes can be represented by empirical determined functions obtained following a model that separates source size from attenuation effects (e.g. Castro et al., 1990;Anderson and Lei, 1994;Castro et al., 1996;Castro et al., 2008). This technique was first proposed by Brillinger and Preisler (1984) for the analysis of attenuation relations of peak accelerations. The spectral amplitudes can be modeled with the following equation: Where U i (f, r) is a datum from event i recorded at hypocentral distance r at a frequency f, S i (f ) is a scalar that depends on the size of the ith event at a frequency f. A(f, r) is the empirically determined attenuation function that describes the distance decay trend. We assume that A(f, r) Figure 2 Data set used. a) Distribution of earthquake magnitude with epicentral distance. b) Histogram of ecents recorded by station. Figure 3 Examples of signal and noise spectra from stations ELO (located on bedrock) and MOR (located on conglomerate). Left panels are spectra from local events (r < 20 km) and right panels spectra from regional events (r = 138.9 km). Continuous lines represent the 80% energy S-wave spectra (both horizontal components) and discontinuous lines are the noise spectra of a 6-seconds time window choose previous to the first P-wave arrival.
implicity contains the effect of both the geometrical spreading and the quality factor Q, but we do not limit their behavior to a predetermined parametrical function, instead, we constrain A(f, r) to decay smoothly with distance. The basis of the smooth decaying restriction lies on the principle that the inelastic properties in the crust tend to decrease the amplitudes gradually with distance and that undulations may be related to other factors such as site and wave propagation effects that would reflect on the residuals when solving equation (1) (Castro et al., 1990;Castro et al., 2008).
The main assumptions of the model described by equation (1) are that A(f, 0) = 1.0, because at r = 0 there is not attenuation and the spectral amplitudes are fully governed by the source term S i (f ). The model assumes that the shape of the attenuation function is the same for all the sources, for a given frequency, regardless of the size of the event (Castro et al., 1990;Castro et al., 2008). Thus, the source factor S i (f ) shifts upwards or downwards the attenuation function depending on the magnitude of the event without modifying its shape. One major advantage of this last assumption is that the observed amplitudes of events recorded at different distances at a given frequency complement each other and permits to define the attenuation function at a wider distance range.
To solve equation (1), we took logarithms at both sides of the equation and formed a set of linear equations for a given frequency of the form: Where u ij = logU i (f, r) is a datum from earthquake i at distance j, s i = logS i (f ), and a j = logA(f, r) is the value of the attenuation function at distance j. We discretized the distance range of the observations using a constant interval. Equation (2) represents a system of equations that is solved by a constrained least-squares inversion. Equation (2) can be expressed on matrix form as: 1 0 0 : : 0 1 0 : : : : : : : : : : : : w 2 : : : : : : : : a j s 1 : : Where w 1 and w 2 are weighting factors provided as input to constrain a 1 = 0 at r = 0, and to weight the second derivate, for smoothing purposes, respectively. The weighting  factor w 1 forces the attenuation function to cross the origin. We tested diverse values of the weighting factors and discretized the hypocentral distances at different intervals in order to achieve the expected monotonic and smooth decay of the attenuation function without losing the characteristic trend of the data.
The resulting model parameters a j define the shape of A(f,r) for each frequency considered. A more detailed description of this method can be found in Castro et al., (1990 and1996). Figure 6 shows examples of attenuation functions found for 10 of the 23 frequencies (0.40-63.1 Hz) analyzed and the respective spectral amplitude data used to calculate them. The amplitudes in Figure 6 were scaled according with the value of s i obtained from the inversion. Figure 7 shows spectral amplitudes from event 9 (M L 3.5) and the corresponding attenuation function (dashed line) scaled with the value of s 9 obtained. This figure illustrates how good the attenuation function obtained for event 9 fits the observed spectral amplitudes.

The quality factor Q using a homogeneous attenuation model
The attenuation functions A(f,r) can be used to analyze various sources of S-wave attenuation and to estimate the quality factor Q. With this purpose, we model these empirically determined functions as the product of two different attenuation mechanisms (Castro et al., 1999): Where 1/r b represents the effect of geometric spreading, N is a normalization distance and the exponential function accounts for the amplitude decrease due to total Q (intrinsic plus scattering). The normalization distance N is a reference distance that is chosen as close to the source as the data permits. Note that (r -N) in equation (4) is the S-wave path distance length where Q takes effect to normalize the spectral amplitudes at the reference distance. The parameter υ is the average S-wave velocity (3.4 km/sec), based on the crustal structure reported by Harder and Keller (2000) for this region. This model has a crustal thickness of 35 km and consists of three layers: the uppermost layer is one km thick and has a S-wave velocity of 2.86 km/s; the second layer has a thickness of 21 km and a velocity of 3.3 km/s; and the third layer, that represents the lower crust, has a thickness of 13 km and a velocity of 3.88 km/s. Different values of N, between 1 and 20 km, were tested to evaluate the influence of this parameter in the geometric spreading behavior.
The parameters b and Q can be estimated for each frequency by linearizing equation (4), taking logarithms at both sides. Thus, for a given frequency f we can rewrite equation (4) as a set of linear equations of the form: Where d i = log A(f,r i ) -logN is the normalized amplitude at distance r i , m i = −πf(r i -N)loge/υ, υ = 3.4 km/s is the average S-wave velocity and c i = −log(r i ). 1/Q and b are estimated by solving equation (5) with a leastsquares inversion scheme.

Attenuation functions
The spectral amplitudes modeled with equation (1) describe the S-wave attenuation in northeastern Sonora, in the region close to the rupture zone of the 1887 (M W = 7.5) earthquake. We determined 23 spectral attenuation functions that represent the decay of the S-wave energy for discrete frequencies defined between 0.4 and 63.1 Hz. Figure 6 shows a sample of these attenuation functions and the observed spectral amplitudes for both horizontal (black circles) and vertical (open circles) components. We observe that spectral amplitudes at low frequencies (up to 5.0 Hz) tend to attenuate less with hypocentral distance than amplitudes at higher frequencies. Figure 7 compares the observed spectral amplitudes (circles) from a M L = 3.5 earthquake (event 9 in Table 1) with the nonparametric attenuation functions (dashed lines) scaled with the corresponding values of S i (f) (equation (1)) resulting from solving equation (2). We also plotted in Figure 7 A(f, r) for S i (f) = 1(solid lines) as a reference to show that for a given frequency the shape of the attenuation function is the same regardless of the event size. In other words, the rate of decay of the spectral amplitudes with distance is assumed to be the same for all the events.
To evaluate the effect of site conditions of the recoding stations, we calculated spectral attenuation functions separately for both vertical and horizontal ground-motion components using the whole distance range (Figure 6). Under the assumption that the vertical component of acceleration is less susceptible to site-amplification effects, a comparison of the attenuation functions obtained using the vertical component provides a reference to evaluate the effect of local site response on the attenuation functions determined with the horizontal components. We found no substantial differences between the observed spectral data or between the calculated attenuation functions ( Figure 6).

Estimates of Q
The attenuation functions can be modeled with three parameters: N, b and Q (equation 4) . N is a reference distance chosen as the closest source-station distance, b define the rate of decay due to geometrical spreading, and Q accounts for the attenuation due to anelasticity, friction losses, multipathing and scattering effects. Figure 8 shows estimates of b and Q, calculated for N = 1, obtained using the attenuation functions for the hypocentral distance range between 10 and 140 km, and the horizontal acceleration spectra. The values of Q (Figure 8), obtained making N = 1, increase with frequency from 27 at 0.4 Hz to 3547 at 63.1 Hz, while the values of b show a weak dependence with frequency. We also calculated Q using values of N of 10 and 20 finding negative values of Q for most frequencies. Thus, we consider that N = 1 is an adequate reference distance, since provides positive values of Q and a good fit to the observed spectral amplitude decay.
Although the idealized far-field body wave geometrical spreading in a homogeneous whole space is given by G (r) = r − 1 , some studies (Ou and Herrmann 1990;Burger et al., 1987;Chapman andGodbee, 2012 Morozov, 2010) indicate that the apparent geometrical spreading in velocity structures other than a homogeneous whole space will generally involve more complex amplitude decay than 1/r. The free parameter b in equation (4) can take different values at different frequencies and accounts for a geometrical spreading that may depart from the value expected for a homogeneous medium. If we adopt a constant geometric spreading model with b = 0.21, that corresponds to the average value of b in the frequency band of 0.5 -63.1 Hz, and recalculate Q (Figure 9 and Table 3 Although the model they used to describe the sources of attenuation involve a term that accounts for the nearsurface attenuation, their Q-frequency relation (equation (7)) follows the general trend of Q obtained in this study (circles in Figure 9). However, between 0.63 Hz and 4.0 Hz equation (7) underestimates the values of Q here determined ( Figure 9). The data sets used in previous studies (Castro et al., 2008(Castro et al., , 2009 and that used in this paper have important differences. First, we used a greater number of events located in the epicentral area of the 3 May 1887 M W 7.5 earthquake; second, we used hypocenters relocated by Castro et al. (2010) and the precise locations guarantee a better control on the source-station distances; third, the volume sampled by the source-station paths in this paper follows the strike of the faults closer than previous studies. Jeon and Herrmann (2004) estimated a Q-frequency relation similar to equation (6) with a geometric spreading coefficient b defined by sections depending on the hypocentral distance (b = 1.2 for 0 < r ≤ 50 km; b = 0.55 for 50 ≤ r ≤ 60 km; b = 0.2 for 60 ≤ r ≤ 90 km; b = 0.1 for 90 ≤ r ≤ 140 km and b = 0.5 for 140 < r ≤ 500 km) in the Basin Figure 9 Estimations of Q for S waves in the region close to the fault zone of the 1887 earthquake. Black continuous line is the mean-square fit of the Q values (circles) using b = 0.21 obtained in this work. Grey dotted line corresponds to the regional estimation proposed by Castro et al. (2008). Discontinuous black line corresponds to the values of Q obtained by Jeon and Herrmann (2004) for the Basin and Range province in Utah (USA). and Range Province in Utah (USA). This function overestimates slightly the values of Q that we estimated between 3 and 16 Hz (Figure 9), but in general their relation: follows a similar trend of our Q estimates. The parameter b depends on the hypocentral distance, and for distances greater than 50 km has values below the theoretical value of one.

Attenuation in the north-south direction
To evaluate the attenuation of S waves in the north-south direction, along the strike of the faults that ruptured during the 1887, we calculated attenuation functions using data from the same site to avoid the possible influence of site amplification effects. We selected station OAX (Figure 1) for having the largest amount of records and for being located inside the fault zone. We formed a dataset with horizontal spectral records of OAX from 33 events with magnitudes between 1.0 and 3.5 that occurred in the 20-90 km distance range from the station. We determined attenuation functions for the same frequencies and procedure (equation 1) as before, and then we used them to estimate b and Q with equation (4). Figure 10 displays b and Q obtained with the attenuation functions determined for station OAX, using values of N = 20 (the closest source-station distance in the data set). We found that Q increases with frequency following the function Q ¼ 17:6 AE 1:4f 1:15AE0:12 ð9Þ while b is weakly dependent of frequency, having a value close to 1.4 for the whole frequency band analyzed (1.3-63.1 Hz). Some studies of body wave attenuation shows that the coefficient b associated with the geometrical spreading is more complex than b = 1 and they attribute this to complex velocity structures (Ibáñez et al., 1993;Olafsson et al., 1998;Castro et al., 1999;Akıncı et al., 2006;Padhy, 2009). Ou and Herrmann (1990) found that G(r) may depend upon source depth in layered structures, while Burger et al., (1987) remarked the importance of post-critical reflections from mid-lower crustal velocity discontinuities and the Figure 10 Estimates of b (left frame) and Q (right frame) for N = 20 obtained with records from station OAX for hypocentral distance in the range of 20-90 km. Table 3 Values of Q (Figure 9) estimated at different frequencies using N = 1 and b = 0.21 in equation (4) Figure 11 Comparison of the nonparametric attenuation functions obtained with station OAX (continuous lines) with expected attenuation calculated using equation (4) with geometric spreading G(r) = r − b (dotted lines) and theoretical geometric spreading G(r) = r − 1 (dashed lines). The circles are the observed horizontal spectral amplitudes used to obtain the attenuation functions.
Moho in controlling amplitudes in certain distance ranges. Chapman and Godbee (2012) found values of the geometrical spreading in eastern North-America up to r −4 for the vertical components and up to r -1.5 for the horizontal components, which exceeds notably the theoretical standards for body waves. In our case, we are finding similar values of the parameter b (for N = 20), for the OAX data set, to those reported by Chapman and Godbee (2012) for the horizontal components ( Figure 10). Malagnini et al., (2002) introduced a slightly frequency-dependent geometrical spreading to model ground motion in northeastern Italy, a region where the attenuation parameters vary considerably. Akıncı et al. (2006) proposed geometric spreading functions that are frequency dependent on a stepped way for frequencies below and above 1.0 Hz and for a different range of distances in the Marmara (Turkey) region. The results obtained in this paper (Figures 8 and 10) show that b is weakly dependent of frequency and that there are certainly other factors, beyond simple theoretical models, to fully account for all the processes affecting the attenuation of S-waves in the Sonora region.
To verify the consistency of our results, we compare in Figure 11 the nonparametric attenuation functions (continuous lines), calculated with equation (2) using the observed spectral amplitudes (circles) from station OAX, with the expected amplitude decay (grey dashed lines) calculated using equation (4) with the estimated values of Q and b, using N = 20 ( Figure 10) for that site. Figure 11 shows also, for comparison, the geometrical spreading function G(r) = r − 1 (dashed lines). The curves calculated with equation (4), with the geometrical spreading model G (r) = r − b , follow the nonparametric functions very closely at all frequencies, while the function G(r) = r − 1 (dashed lines) underestimate the observed amplitudes at short distances (r < 40 km) and overestimates the observations at larger distances (r > 40 km).

Summary and conclusions
We used precise hypocenter locations from local events to find nonparametric spectral attenuation functions for the region that ruptured with the M w = 7.5 1887 earthquake of Sonora. The attenuation functions determined decay more slowly with hypocentral distance at low frequencies (f < 4 Hz) than results reported previously in the area using more distant recordings. Consequently the Q-frequency relation (eq. (7)) obtained by Castro et al. (2008) tends to underestimate Q between 0.63 and 4.0 Hz. The values of Q estimated for the S-waves (Figure 9) show a clear dependence of Q with the frequency, and agree with the models proposed by Castro et al. (2008) for the region, and with the values of Q calculated by Jeon and Herrmann (2004) for the Basin and Range Province in the state of Utah (USA).
The values of Q estimated with station OAX (eq. (9) and Figure 10) for S-wave paths traveling along the strike of the fault system located near the rupture of the 1887 event, in the north-south direction, are considerably lower than the average Q estimated using sourcestation paths from multiple stations and directions (eq. (6) and Figure 9). These results indicate that near the fault zone S waves attenuate considerably more than at regional scale, particularly at low frequencies. For instance, at 0.5 Hz Q in the north-south direction, along the strike of the faults (eq. (9)), is 10.5 times smaller than the average Q (eq. (6)). This may be the result of strong scattering near the faults due to the fractured upper crust and higher intrinsic attenuation due to stress concentration near the faults.
The geometric spreading models found have weak frequency dependence (Figures 8 and 10) and can be approximated as G(r) = 1/r 0.21 for the average Q. This spreading function predicts slower amplitude decay with hypocentral distance than the r − 1 theoretical model of body waves.