Introduction

Seismic hazard in Mexico is mainly related to the subduction zone in the Pacific Ocean coast. Significant effort has been devoted to understand plate motion in this complex geological region and the characterization of its dynamic behavior (^{Pérez-Campos, 2008}). However, estimation of ground motion variations due to local site conditions is of great engineering interest, with an important contribution to studies of seismic hazard, where the main purpose is to reduce human and economic losses that ground motions may lead to civil structures (^{Kramer, 1996}).

In order to estimate local site-response contribution to the recorded seismic ground motion, it is necessary to remove the source and path effects. Several source and path effects inversion methods have been developed (^{Field and Jacob, 1995}), each involving assumptions in order to constrain source and/or path effect estimations. Most of these methods include a site response term, but their main goal is the enhancement of source and path effects (^{Field and Jacob, 1995}). In this work, the interest is focused on an accurate estimation of the site response at southern Mexico by using high quality seismic records from a temporal and permanent very broad band array, but also in source and path effects to a minor degree.

In this paper a generalized-inversion scheme to estimate local site response effects was used (^{Field and Jacob, 1995}), where shearwave amplitude spectra is represented by a parameterized source and path effects model and a frequency dependent site response term for each site. Estimating the local site response using this non-reference site dependent technique involves solving a strongly nonlinear multiparameter problem. Solutions to this parameterized inversion scheme were introduced by using two global optimization methods, also called heuristic methods: genetic algorithms (GAs) and simulated annealing (SA). This methods allow to explore a wider range of non-linear potential solutions, as compared to other conventional linear methods (Rodriguez Zúñiga et al., 1997; ^{Ortiz-Alemán et al., 2004}).

An inversion problem solution by using heuristic techniques, such as GAs or SA, requires setting parameters of initial models, including upper and lower bounds to their feasible values. In the case of the sourceeffect term, parameter values and their bounds were taken from the scalar seismic moment as reported at Global Centroid-MomentTensor (CMT) and SSN catalogues. For the path effect term, the theoretical arrival time of the S-wave and the attenuation function for earthquakes traveling from the coast towards the continental region (^{García et al., 2009}), were considered. An initial local site response term was set by using the non-reference site dependent technique proposed by ^{Lermo and Chávez-García (1993)}, which is a variation of the spectral ratio method. Initially proposed by ^{Nakamura (1989)} to analyze Rayleigh waves from microtremor records, this technique consists in taking spectral ratios between horizontal and vertical components of the ground motion records. ^{Field and Jacob (1995)} compared various site response estimation techniques, including this one, finding the basic assumption of local geology as being relatively transparent to the motion observed on the vertical component supported by the revealed frequency dependence of site response. Although they pointed out that investigations based only on the vertical component may underestimate the local site response at sedimentary sites.

Finally, the source, path and site parameters corresponding to the S wave amplitude spectra of 219 velocity seismic records (162 from MASE and 57 from SSN) were estimated by means of the previously described generalized inversion approach. The simultaneously inverted parameters from the records of 55 stations showed a significant agreement between the observed and calculated spectra, which supports the reliability of such estimates. A comparision between local site response terms estimated through both techniques, H/V ratios and the global inversion approach, as introduced in this work, is discussed at the conclusions section.

Data acquisition and processing

Two different seismic networks provided the velocity records used in this work. The first one was the temporary network deployed as part of the MASE project. In the first stage of this project, one hundred broadband seismometers were located in a profile across the southcenter of Mexico (Figure 1), from December 2004 to July 2007. The second seismic network was the permanent Mexican SSN, formed by an array of 23 seismic stations covering the whole country.

Figure 1 Geographic location of the 55 seismic stations (14 SSN - 41 MASE), and epicentral location of the six seismic events considered in this study. 

Original data sets consisted in 68 seismic events, with at least one recording in any of the 123 stations from the total array (100-MASE and 23-SSN). We carefully selected seismic records, discarding those with undesirable (less than 2) signal-to-noise ratios.

We selected broadband seismograms from six earthquakes (four subduction and two inslab events) which were simultaneously recorded by both seismic networks from August, 2005 to April, 2007 in southern Mexico. The magnitudes of these events range from 5.0 to 5.6 Mw and their focal depths from 10 to 80 km. Also, we constrained the analysis to a sub-array of only 55 stations (41 from MASE and 14 from SSN). Figure 1 shows epicenter locations of these intermediate depth earthquakes, numbered from 1 to 6 (stars), and locations of the 55 broadband stations (triangles for SSN and squares for MASE). Table 1 contains the main features of the six events depicted in Figure 1.

Table 1 Main features of the intermediate-depth seismic events considered in this work. 1,4,5, and 6 correspond to subduction events and 2 and 3 to inslab earthquakes. 

Selected records were properly corrected by removing the mean and linear trends, and avoiding any other filtering. For all the records, the theoretical arrival time of the S-wave was picked in the horizontal component. Using these theoretical picks, for each one of the three components, we selected a window including the intense part of the ground motion (20s). These windows were Fourier-transformed by an FFT routine using 2048 samples, and then smoothed by 1/3 octave-band filter.

An initial site response for each station was estimated by computing the spectral ratio technique with a single station (^{Lermo and Chávez-García, 1993}), that involves dividing horizontal to vertical components of the velocity records, as shown in equation (1). In the case when some seismic station recorded more than one seismic event, site response was computed as the average of the spectral ratios from all recorded events,

1

where H and V are the horizontal and vertical components, respectively. And EW, NS, and Z are the east-west, north-south and vertical components, respectively.

Spectral decay analysis

As mentioned above, a window including the intense part of the seismic motion and beginning at the shear-wave arrival was employed for signal spectra estimation. Fourier amplitudes were computed for both vertical and horizontal components from the records. In Figures 2 and 3, spectra were plotted for 8 frequency values at 0.1, 0.5, 1, 2, 4, 6, 8, and 10 Hz, respectively, as a function of epicentral distances to the MASE array stations. Continuous black line is the expected attenuation curve by just considering geometrical spreading with a (1/r) law.

As can be noticed in Figures 2 and 3, spectral amplitudes do not match the expected attenuation with distance (1/r). In such figures, it is possible to clearly observe two regions with different amplitude behaviours. First, in the region from r = 200 to approximately 400 km, we observe a consistent amplification of amplitudes in all propagated frequencies, above the geometrical spreading; this is due to the presence of the Trans Mexican Volcanic Belt (TMVB) (see Figure 1), where the strata consist mainly in volcanic rock. Farther from r = 400 km to the north, a greater attenuation than 1/r, particularly for higher than 2 Hz frequencies, was found. This might be due to the complexity in the response at each site of the MASE array, maybe as a result of the tradeoff between large scale amplification patterns (^{Ordaz and Singh, 1992}), and the particular attenuation features at each site.

In Figure 2, amplitudes corresponding to events 1, 4, 5 and 6 were plotted against epicentral distance. In these events, seismic waves traveled perpendicular to the coast towards the MASE array. Events numbered as 1 and 5 have similar epicenter locations, and they show similar amplitude decay behavior with distance for all frequency windows. Amplitudes from seimic event 6 exhibit a similar behavior relative to the other three events (1, 4 and 5) up to 4[Hz], but from 6 to 10 [Hz] amplitudes are higher.

Figure 2 Spectral amplitude [cm] as a function of epicentral distance. Eight frequency values in the range from (0.1 - 10 [Hz]) were plotted for each one of the MASE array seismic stations. The continuous black line represents the expected trend for geometric attenuation. Subduction events labelled as 1,4,5,and 6 (see Figure 1), with similar seismic propagation paths towards the MASE array were depicted. 

In Figure 3, amplitudes for events 2 and 3 are plotted versus distance. Both events propagate almost perpendicular to the MASE array. Event 3 was located to the north and event 2 to the south of the MASE array. Both events show similar spectral decay with distance following the expected trend 1/r.

Figure 3 Spectral amplitude [cm] as a function of epicentral distance. Eight frequency values in the range from (0.1 - 10 [Hz]) were plotted for each one of the MASE array seismic stations. The continuous black line represents the expected trend for geometric attenuation. Inslab events labelled as 2 and 3 (see Figure 1), with similar seismic propagation paths towards the MASE array were depicted. 

Parameterized source- and path-effects inversion scheme

Many seismological studies consider a siteresponse term primarily to remove this effect, in order to enhance source and/or path effect estimates. In this work the main interest was the site response estimates themselves. For this last reason, a parameterized source and path effects inversion scheme were used.

Following the formulation from ^{Field and Jacob (1995)}, for a network of I sites over which J events were recorded, amplitude spectrum of the j-th event recorded at the i-th site, A_ij(ƒ), can be written as

2

where E_j(ƒ) represents a source term, P_ij(ƒ) a path term and S_i(ƒ) a site effect term.

In this work the spectral method from ^{Brune (1970)} was used, as a quantitative evaluation tool for source parameters. Using Brune’s model, source parameters can be inferred from the low frequency spectral level, Ω₀, and the corner frequency, ƒ_c, measured from a source radiation pattern and Q factor corrected spectrum. In most seismological studies, spectra are corrected by some estimation of Q, and low frequency spectral levels and corner frequencies are obtained by visual inspection of the observed spectra. When using such kind of approach, uncertainties in Q, ƒ_c and Ω₀ estimations can be large, leading to large uncertainties in source parameters. Here a quantitative approach to source, path and site effects inversion based on the use of a simulated annealing algorithm were introduced.

Following ^{Boatwright et al. (1991)}, the source term in the frequency domain is expressed as

3

Equation (3) is ^{Boatwright’s (1978)} approximation for the ^{Brune (1970)} far-field shear-wave spectrum.

It is important to recall the relation between Ω₀ and M₀:

where p is the density, c the wave speed, r the distance from the source, and U_φθ the term on the radiation pattern. Generally, the focal mechanism corresponding to the radiation pattern U_φθ is unknown; nevertheless, the seismic scalar momentum, M₀, can be estimated by replacing U_φθ with a mean radiation term with a value of 0.52 for the P-wave and of 0.63 for the S-wave (^{Shearer, 2009}).

Assuming path effects to be similar for all array sites, except for a geometrical spreading term of the form

4

where R_x was taken as 100 km for this study, the path effect term is modeled as

5

where 1/r_ij is the geometrical spreading factor described in equation (4), T_ij is the travel time of S waves, Q(ƒ) is a frequency-dependent quality factor, taken as Q(f) = 211ƒ^0.46as proposed by ^{García et al. (2009)}, and t_i* is a local attenuation term.

In such a way, 2J source parameters ( Ω₀ and ƒ₀ for each event), I local attenuation parameters (t_i*), and I frequency dependent terms corresponding to the site response, were considered. In this case, observed spectra involved N discrete frequencies, then the model contains 2J + I (N+1) parameters to invert.

Site response estimation using simulated annealing

An initial parameterized model was proposed. Source effect parameters, Ω₀ and ƒ_c, for each event j, were estimated by using the scalar seismic moment, M₀, reported in the Harvard CMT Catalogue and SSN database. Initially, fc_j ∝ M₀ ^−1/3 and Ω₀ ∝ ƒM₀; nevertheless, the following relations where used:

Parameter values are listed in Table 2. Both estimations were considered as lower and upper limits in the inversion scheme, while average values were set as initial parameter values.

Table 2 Scalar seismic moment, M₀, for the six events also listed and described in Table 1. Information was taken from the Seismic Seismologic Network (SSN) and the Global CMT catalogues. Values for f_c and Ω₀ were computed from M₀. 

In order to set initial values of path effect parameters, the relationship between local attenuation and frequency dependent quality factor, as described in equation (6), was considered,

6

Local attenuation t_i*, also characterizes the attenuation of body waves, and might be estimated using the travel time of S waves as well as the quality factor Q (^{Shearer, 2009}).

^{García et al. (2009)} proposed two attenuation functions Q(ƒ) for the S-wave group, considering trajectories along the coast and towards the continent, and therefore values of Q⁻¹(ƒ), for both types of seismic stations, coastal and inland. Values of Q⁻¹(ƒ)±1∑.D. of the final models for the coastal and inland subsets for a frequency range form 1 to 10 Hz are showed in Figure 4. We considered the maximum and minimum value of Q−1(f) as upper and lower limits, respectively, for the path parameters.

Figure 4 Estimated values of Q⁻¹ for coastal (triangles) and continental seismic paths (squares). The dashed line indicates the attenuation function Q(ƒ) for coastal paths and the continuous line the continental paths, after ^{García et al. (2009)}. 

An initial guess for the site response term was considered from results of horizontal to vertical spectral ratios, described in previous sections. Once the initial model has been set, the inversion process followed the flow chart depicted in Figure 5, which considers an inversion approach based on the SA method.

Figure 5 Schematic drawing illustrating the method of simulated annealing as used in this study (after ^{Ortiz-Alemán et al., 2004}). 

A well documented issue in this class of inverse problems is the non-uniqueness of the solution. In this case, an observed spectra can be matched to a similar approximation order by a pair of small fc and t_i* values, or by large values of the same parameters. Nevertheless, in the case of site response estimation this trade-off is advantageous to some extent. A very precise estimation of the real values of ƒ_c and t_i* is not necessary because both of them have a minor influence on the site response term. Uncertainties associated to these tradeoff can be reduced if several records are simultaneously inverted. In addition, the initial model as proposed in this work constrains source parameters by including all available a priori information, then the trade-off effects between ƒ_c and t*, previously pointed out, may be reduced to some significant extent.

Three types of inversion were performed in order to include all possible combinations between the number of stations, I, and events, J. Data were fitted by means of equation (2), considering the parameterized expressions for source- and path-effects, by minimizing the function in equation (7),

where A_ij(ƒ_n) and Ae_ij(ƒ_n) are the observed and estimated amplitudes for spectrum i, respectively, at the frequency ƒ_n ; N is the number of spectra samples that were simultaneously inverted and J is the total number of records that were used.

Attenuation function estimation

Another inversion was performed considering source and local site-response parameters estimated in the generalized-inversion scheme including all stations and events (55 stations and 6 events) in order to estimate a new attenuation function for trajectories towards the continent.

Only path-effect parameters were assumed as unknown, while source effect parameters and local site response values were directly fixed from the solution previously achieved by using SA. The local attenuation term values, t_i*, and their limits remained similar to the generalized inversion scheme. For the attenuation function of the form Q(ƒ) = aƒ ^b, values of the parameters a and b and their limits were taken from other two attenuation studies, ^{Ordaz and Singh (1992)} and ^{García et al. (2009)}, that are routinely used for estimation of regional attenuation at southern Mexico. Attenuation parameters, a and b, were then inverted, which allowed to accomplish a new attenuation function, Q(ƒ) = 224ƒ ^0.55. This new function can deal with any kind of seismic events (inslab or subduction earthquakes) and their corresponding seismic propagation paths along the coast or towards the inner continental zone. This effect is notably shown in Figure 4, where intraslab and interplate events are depicted; notwithstanding that continental models present a higher attenuation than coastal ones, the differences are not considerable, being able to be adjusted by the linear function herein provided.

Notable differences in attenuation along the coast and towards the continent have been reported by ^{García et al. (2009)}. This study suggests that seismic waves travelling inland suffer considerably less attenuation than those propagating along the coast.

For comparison purposes, error functions as in equation (7), for each seismic station were computed by using the attenuation function proposed by ^{García et al. (2009)} and the one obtained in this study. Results are depicted in Figure 6, and it can be observed that along the profile across the south-center of Mexico (MASE array) and also in the case of permanent stations of the SSN, errors are lower in general by using the attenuation function from this work, which is also valid for all propagation paths and any kind of seismic event, as mentioned before. In Table 4 the main similarities and differences bertween datasets used by ^{García et al. (2009)} and by this work are listed, for the computation of a frequency dependent attenuation function. This table allows to appreciate the robustness of the present frequency-dependent attenuation function, due to the fact that the database of stations and records used here is of comparable size relative to the one used by ^{García et al. (2009)}.

Table 3 Local atenuation, t*, corner frequency, ƒ_c, and objective function e, from inversions A, B and C, for the seismic stations of the MASE array (HUIT, JIUT, PSIQ and SATA) and SSN seismic network (CUIG, MOIG, PLIG and TPIG). 

Table 4 Comparison between datasets from ^{García et al., (2009)} and this work, for the computation of frequency dependent atenuation function, Q(f). 

Figure 6 Objective function value, e, computed for the attenuation function. In the upper part, stations of MASE array were located considering their distance relative to the coast, being QUEM the nearest and CHIO the farest station. In the lower part, SSN network stations in alphabetical order. 

Results and conclusions

In this paper a parameterizd stochastic global inversion method using of local seismic response was introduced, which also allows to estimate parameters of source and path effects. The technique of horizontal to vertical (H/V) spectral ratio is used to obtain an initial model for the inversion scheme. This type of initial solution approaches qualitatively local seismic response particularly in the case of hard rock sites. However in the case of sedimentary sites, certainty of estimates for site response from H/V spectral ratio is significantly conditioned by the impedance contrast between sedimentary and basement units, as well as by the three dimensional geometry of the sedimentary basin. H/V spectral ratio in the presence of a low impedance contrast presented is not suitable to use as pointed out by ^{Luzón et al. (2001)} and ^{Sánchez-Sesma et al. (2011)}.

Parameterized nonlinear global inversion allowed us to obtain robust estimates of local site response, by simultaneously inverting the whole set of observed seismic records. Parameters were estimated in a frequency range of 0 - 10 [Hz] for each of the seismic stations (41 - MASE and 14 -SSN). Estimates of corner frequency and low frequency spectrum level for the six seismic events were obtained. Local atenuation values and a regional attenuation function were also introduced.

Uncertainty inherent to the non-uniqueness on the solution was restricted by incorporating a priori information about seismic source, path and site parameters. The kind of parameterization used in this work to estimate such seismic effects was quite appropriate as sensitivities from all estimated parameters show significant contributions on the shape of synthetic spectra computed by the forward problem calculation. Observed spectra from stations of MASE array (CUIG, PLIG, HUIT, JIUT, PSIQ and SATA) were fitted to a similar degree when compared to spectra from stations that were located far and even at opposite ends (MOIG, TPIG), as shown in Figure 7.

Figure 7 Observed, A(ƒ), and synthetic spectra, A_c(ƒ), computed by using estimated parameters from inversions A, B and C, for some stations of the MASE array (HUIT, JIUT, PSIQ, and SATA) and SSN network (CUIG, MOIG, PLIG and TPIG). In each case, local atenuation, t*, and objective function value, e, are indicated. 

Three different inversions of site response were computed. In study case A, the seismic response for each one of the 55 seismic stations was separately inverted and all six seismic events were simultaneously considered. In case B, the seismic response was inverted by considering each one of the six seismic events and the whole set of 55 seismic stations simultaneously. In case C, the site response at 55 seismic stations was inverted and seismic records for six events were considered. Results are summarized in Table 3, where estimated source and path effect parameters are listed for a set of MASE and SSN seismic stations.

In Figure 8, site response estimations at some selected seismic stations from MASE and SSN seismic networks were plotted. For comparison purposes H/V average site effects for the same set of stations are also included. Three different inversion efforts (A, B, and C) are include in the same Figure 8.

Figure 8 Estimated site response, S(ƒ), for seismic stations of the MASE array (HUIT, JIUT, PSIQ and SATA) and SSN network (CUIG, MOIG, PLIG and TPIG), as compared to the site effect computed by average H/V ratios. In each case, the objective function value, e, is indicated. 

In order to account for a possible trade-off between source and path effects, affecting the site effect estimations, in Figure 9 a comparison is made for the estimated source spectra for seismic event listed as 6 in Table 1. One station (CEME) near the seismic epicentral location was used for the inversion of source spectra. In order to establish the approximation quality, in Figure 9 the observed and estimated global seismic spectra (considering source, path and site effects) is are shown. At the top of Figure 9, the estimated source term for each seismic station is plotted. The effects of attenuation with distance and local amplification can be considered as less significant at station CEME. In this same figure, source effect can be more clearly distinguishable in this station which is closer to the epicenter. A good agreement between the estimated source and the global spectra was found, on one hand, and the observed spectra on the other.

Figure 9 Estimated and observed seismic spectra for CEME and PLIG stations of the MASE array. In the upper part, the inverted source term for event 6, listed in Table 1 is included. Estimated global seismic spectra account for source, path and site effects. Attenuation with distance is less significant at CEME in comparison to PLIG, as expected, and source effect is more distinguishable at CEME station. The synthetic spectra closely resemble the observed spectra. 

Atenuation Function ^{García et al. (2009)} This work Number of records 469 219 Type of records Velocity and acceleration Velocity Number of stations 56 55 Magnitude of events, Mw 5.0≤Mw≤8.0 5.0≤Mw ≤5.6 Depth of events, H, [km] 10≤H≤30 10≤H≤80 Distance from seismic source to stations, r, [km] 20≤r≤400 50≤r≤500

In Figure 6, the objective function value, computed for the attenuation function was plotted. Attenuation function by ^{García et al. (2009)} is considered in case C. Another inversion case study (D) was then computed, by optimizing coefficients a and b of a general frequency dependent attenuation function of the form Q(ƒ)=aƒ ^b. In Figure 10, observed and synthetic spectra, A(ƒ) and A_c(ƒ) were included, computed by using the attenuation function from ^{García et al. (2009)}, labeled as C, and the attenuation function as derived in this work, labeled as D, for stations from the MASE array (HUIT, JIUT, PSIQ and SATA) and SSN network (CUIG, MOIG, PLIG and TPIG).

Figure 10 Observed and synthetic spectra, A(ƒ) and A_c(ƒ), computed by using the attenuation function from ^{García et al. (2009)} , labeled as C, and attenuation function as derived in this work, labeled as D, for stations from the MASE array (HUIT, JIUT, PSIQ and SATA) and SSN network (CUIG, MOIG, PLIG and TPIG). In each case, value of objective function, e,is indicated. 

It is important to recall that pioneering work on the spectral amplification at the hill zone around Mexico City was done by ^{Ordaz and Singh, (1992)}. They made evident a systematic spectral amplification in the frequency range of ƒ = 0.2 − 2.0 Hz. The results herein obtained confirm that geological conditions represented by the TMVB, specially in the case of subduction events, constitute a region of relative spectral amplification in the frequency range of ƒ = 0.5 − 4.0 Hz. For epicentral distances larger than 400 km (in the north side of the MASE array, beyond the TMVB) there exists a strong spectral decay in seismic energy for the frequency range of ƒ = 4.0 − 10.0 Hz. Nonlinear inversion of source, path and site effects from the MASE experiment dataset probed to be an alternative tool to derive meaningful information on relative seismic amplification or attenuation in south-central Mexico.

Data and resources

Two different seismic networks provided the velocity records used in this work. The first one was the temporary network deployed as part of the MASE project. In the first stage of this project, one hundred broadband seismometers were located in a profile across the south-center of Mexico (Figure 1), from December, 2004, to July, 2007. The second seismic network was the permanent Mexican SSN, formed by an array of 23 seismic stations.

Seismic data from the MASE temporal seismic network and from the permanent SSN network were provided by the Insitute of Geophysics of the National University of Mexico (UNAM). Seismic modeling and data inversion were carried out at the Mexican Institute of Petroleum.