1. Introduction
The general knowledge of seismicity of the Mexican Pacific Coast has been described in terms of its relationship within the regional context with the Middle America Trench, the convergence rate observed and predicted by plate tectonics processes, the use of seismic gap theory for forecast earthquakes on the region (Reyes [1]; Eissler and McNally [2]; Singh [3]; as well as Heaton and Kanamori [4]; among others) and prediction of strong ground motions using empirical relationship and seismic source scaling procedures (Somerville [5] and Somerville [6]).
The work here presented, mainly deals with the Colima-Jalisco region (ruptures areas of 1932, and 1973 earthquakes) as well as the remained gap of the Tecoman (2003) earthquake. The conduction of studies in this region to simulate seismic scenarios acquires special significance because of the hazard of an eventual large magnitude earthquake. For the large earthquakes recording in this region last and present century the rupture areas of the 1932, 1973 earthquakes and the small gap remainder from Tecoman earthquake the seismic convergence period have been exceeded or is in its limit.
Let us look into this with detail. Firstly, for the earthquakes of 3 and 18 June 1932, Nishenko and Singh [7] found that the average displacement on the fault of 1932 earthquakes as 155 cm. Pardo and Suarez [8] found that subduction rate of Rivera plate in this region is estimated to be from 2 to 5 cm/yr, take the small estimation of convergence rate of 2 cm/yr yields a recurrence period of about 77 years assuming that the convergence is entirely taken up by the seismic slip (Nishenko and Singh [7]). However considering a rate larger than 2 cm/yr the recurrence period could be much lower. The 79 years elapsed since 1932 indicate that the recurrence period has expired. Although the rupture area of Manzanillo 1995 Mw 8.0 earthquake invades the rupture area of 1932 earthquakes, but represent less than 15 % of rupture area. Secondly, for the January 30 1973 earthquake (Mw 7.3) Reyes [1] estimate an average slip of 144 cm. For this region Minster [9] estimate a convergence rate of 5.6 cm/yr, this suggest a recurrence time period of about 25 years Reyes [1]. This recurrence period consider the various uncertainties and is in very good agreement with the interval between this earthquake and the preceding 1941 earthquake (32 years). The 38 years trascurred since 1973 indicate that the recurrence period have expired Reyes [1].
Thirdly, is the existence of a seismic gap in the region located between the rupture areas of the 1973 and the Tecoman earthquakes. In figure 1 Quintanar [12] shows that the aftershocks location of the Tecoman earthquake lies north of El Gordo graben and the aftershock area encompass part of the rupture area of the 1932 and 1995 earthquakes. The area between the limits of the rupture areas of the 1995 and the1973 earthquakes is what has been called the Colima seismic gap. The northwest area of this gap ruptured with the Tecoman earthquake in 2003. The other half of the gap, roughly to the southeast, remains quiet Quintanar [12].
The later show the existence of 3 different zones where is necessary to make a simulation of strong ground motions to estimate the response spectra. From these 3 zones the more important (by magnitude and area) is the region broken by earthquakes of 1932; however, the possibility of conduct a simulation in this area is limited by the poor instrumentation in the past and therefore the absence of data the region.
In the case of small gap remainder from Tecoman and the rupture area of 1973 earthquakes, the Mw 5.3 earthquake of August 13, 2006 (element event) offered the opportunity to generate a seismic scenario for the studied area that constitutes a potential seismic source in this region in the near future. We apply: (i) the empirical Green’s function method (EGFM) proposed by Irikura [13] to estimate the peak ground acceleration (PGA) and the acceleration time histories using the Somerville [6] relationship allowing us scaling from moderate to major earthquake. Traditionally the EGFM is well known technique commonly used to simulate an already occurred earthquake. The scientific community thinks that would be of great benefit if this technique can be applied to forecast an expected earthquake in Mexico. Under this premise, an immediate doubt arise because the absence of observed records (event we intent simulate) to compare with synthetics. Conduct a study to simulate strong ground motions under this situation makes difficult to constrain and validate the results, but at the same time is one of the more important challenge simulating a future earthquake. In order to overcome the absence of observed records two methodologies were applied: (i) The first one was Somerville [6] relationship together with the EGFM, and (ii) the second one was ground motions predictive equations (GMPE) relationships to make accurate estimation of strong motions. Aguirre and Irikura [14] used acceleration records of Mexican subduction earthquakes to validate the Somerville [5, 6] relationship for subduction earthquakes and found that the asperities size were well predicted by this relationship. The EGFM is a well established methodology used in the field of earth sciences to estimate the ground motions. What the EGFM requires to estimate the ground motions according to the tectonic conditions of the region, are the fault parameters determination. Works like the one done by Irikura [15] address the study of the factors involved with the inner and outer fault parameters of the source to do more accurate estimations of ground motions. Those outer and inner fault parameters are estimated from the inversion of the waveform in studies of rupture processes using strong ground records. Using the results of many kinematic inversions, Somerville [6] obtained some relationships that synthesize the main characteristics of the earthquakes, stated as follows: (i) the seismic moment, (ii) the rupture area, (iii) the slip average, (iv) the combined area of asperities, and (iv) the area of the largest asperity, among others. In this study we are using these relationships because of their great utility in the estimation of seismic ground motions. The above statement provides an efficient way to work when a limited number of parameters to be considered in numerical simulations is available.
Finally from both estimations; it is said the PGAs, and the curves obtained with the two GMPE here used: (Ordaz [16], and Young’s [17]), an trial and error iterative process of residuals minimization was conducted to identify the result that in the statistical sense better matched with the two GMPE here used.
The main contribution of this method is that it reflects a model that considers the source, the path, and the site effects. Another important contribution of the method is that reliable estimations about the energy distribution can be achieved in the high frequency band (between 0.1- and up to 10-Hz). This frequency range is of engineering interest because of the following reasons: (i) Many structures, including tall buildings and long bridges have their natural frequencies in the above frequency range, and (ii) 8 of the 10 major cities of this state are located in the sedimentary basins of the Colima graben and could amplify the ground motions in the frequency range of 0.1 to 10 Hz. It is therefore important to investigate how ground motions up to 10 Hz are generated from great subduction-zone earthquakes. This kind of investigations play a vital role in the effort to propose an scenario of strong ground motions from future large subduction earthquakes in the area in study and to evaluate the performance of structures subject to ground motions.
2. Tectonic
The Colima state is located in Mexico’s Pacific coast. The tectonic of the region is complex, in which the Rivera, the Cocos, and the North American plates converge. In addition to the above, the existence of a microplate has also been proposed by DeMets and Stein [18], and Bandy [19]. There are significant changes in the parameters of the subduction process along the subduction zone on the Pacific coast of Mexico, which has been divided in four sections by Pardo and Suárez [8]. Although the dip of the interplate contact geometry is constant to a depth of 30 km, lateral changes in the dip of the subducted plate are observed once it is decoupled from the overriding plate. In front of the Jalisco block, the Rivera plate has a dip of 45° and its subduction rate below the North American Plate is estimated to be from 2 to 5 cm/yr. The Cocos plate below Colima shows a similar dip to that of the Rivera but the subduction rate below the North American Plate is estimated to be from 4 to 6 cm/yr. To the south, the dip of the Cocos plate decreases gradually and is almost sub-horizontal at Guerrero (where it subducts with a velocity from 6 to 7 cm/yr) before increasing again farther south to the large values observed in Central America. Pardo and Suárez [8] explained the observed no parallelism between the volcanic belt and the subduction zone by these large lateral variations.
3. Data
To simulate acceleration time histories of the Tecomán earthquake, we used two sources of data. Firstly, records of the 13 August 2006 earthquake. These records were obtained from previous temporal campaign in that area and from the support of other institutions with permanent instrumentation in that zone. Its epicenter, focal mechanism, and seismic moment were obtained from Centroid Moment Tensor Project [20]. The location of instruments that recorded this event (figure 1) is next described. The instruments were from permanent seismic networks: 15 Etna episensor wideband accelerographs from d.c. to 200 Hz at 200 samples per second from the national accelerations network of Instituto de Ingeniería (IINGEN) of Universidad Nacional Autonoma de Mexico (UNAM); two Guralp CMG40T-DM24 flat response wideband velocity type seismographs from 0.5 to 100 Hz at 100 samples per second from the network Red Sismica del Estado de Colima (RESCO). Secondly, data from temporal networks installed in the region as part of this project as follow: (i) four Altus Etna wideband accelerographs from d.c. to 100 Hz at 100 samples per second, four (ii) Geosig strong-motion recorder model 18 with analogue-digital converter, wideband accelerometers from d.c. to 100 Hz recording at 100 samples per second. Because 2 of the 25 records used in this study were velocity records, it was necessary to transform them to acceleration. Also, it was necessary to remove the instrumental response of each of the different instruments.
4. Method
The method used to model the target event requires a small magnitude event (earthquake of 13 August 2006) called element event, with hypocenter in close proximity to the earthquake that we want to simulate. For this particular case, the magnitude, location, focal mechanism and source parameters of the element event was reported by Harvard CMT (Mw 5.3, 18.45°N latitude and -103.63°W longitude, depth 23.5 km, strike 38°, dip 23°, and slip 96°, seismic moment 1.12e24 dyne-cm). For the simulated event (target event) and taking in consideration that the area in study is the region between the limits of the rupture areas of the Tecoman (2003) and the 1973 earthquakes, the hipocentral location was proposed just inside of the area in study and near to element event (18.45°N latitude and -103.75°W longitude). Considering the rupture area of 1973 earthquake and area of remainder gap of Tecoman earthquake we proposed 70 km along strike of fault area. Along the dip, we propose 80 km considering an intermediate value of dip length of neighbors earthquakes. Based on the above considerations the proposed effective rupture area is of 5600 km^{2}. Using equation (1), Somervile [6] we estimate a seismic moment of 1.1091e27 dine-cm.
Where A is the rupture area and Mo is the seismic moment.
Using equation (2) by Kanamori [21], the maximum estimated M_{W} magnitude is 7.3.
Where Mw is the moment magnitude and Mo is the seismic moment.
We use the relationships of Somerville [6] to characterize the source parameters as follows: (i) equation (3) to estimate the combined area of asperities (A2), (ii) equation (4) to estimate area of largest asperity (A1), (iii) equation (5) to estimate the hipocentral distance to center of closest asperity (CA), and (iv) equation (6) to estimate the rise time (Rt) that is related to seismic moment of a small earthquake. For the S-wave propagation velocity we used 3.4 km/s.
The fault plane was defined considering: an azimuth of 38°, a dip of 23° and a slip of 96°. These parameters were taken assuming that the mainshock will have the same focal mechanism as the element earthquake.
To estimate the number of sub-events, we applied the ω^{-2} spectral model, Aki [22], obtaining the number of sub-events necessary and estimate N^{3} by using the relationship between seismic moments of the target event (M_{0}), and the element event (m_{0}) that is used as empirical Green's function. N^{3} is equal to the number of sub-faults in direction of the strike (Nx), the dip (Nw) and the time (Nt).
The above description clearly states that it is necessary to find the parameter N, which will be used to scale the fault area for the event to simulate. Since it is divided into N x N subfaults, N^{3} is obtained using the equation (7), and the relationship between these parameters is stated through equations (7), (8), and (9).
where Ū_{0} and ū_{0} is the flat level of the displacement Fourier spectrum for the target and element events respectively. On the other hand M_{o} and m_{o} are the seismic moments of the target and element events respectively. The relationship for high frequency is given by:
where Ā_{0} and ā_{0} are the flat level of the acceleration Fourier spectrum of the target and element events respectively.
Then the synthetic motion of the target event A(t), is given by the element event a(t) using the equations 10 and 11.
where n’ is an appropriate value to eliminate spurious periodicity, r is the distance from the station to the element event hypocenter, r_{ij} is the distance from the station to the sub-fault (i, j), t_{ij} is the sum of the delay times due to the rupture propagation and the differences of distances between the location of the element event and the location of the target event (i,j) at the observed site.
During the simulation process, we use the Somerville [6] relationships to assume and vary the inner and outer fault parameters in order to simulate acceleration records. Such parameters are: the rupture velocity, the rise time and the point where the rupture starts among others.
The PGAs were estimated for the three orthogonal components. On the other hand, earthquake magnitude, focal depth, hypocentral distance, and site characteristics (rock or soil) were the controlling parameters to estimate curves of ground motion by using the GMPE from Ordaz [16] and Young’s [17] used in this study. To compare our results with respective GMPE the PGAs of two horizontal components was computed according with table 1.
GMPE | Horizontal components | Equation |
Youngs et al. (1997) | Geometric mean ^{*} | |
Ordaz et al. (1989) | Cuadratic mean |
After the above steps were completed, we proceeded to generate and compare the mean value of the residual between the PGAs and each one of the GMPEs by applying the definition of mean residuals as the weighted sum of the residuals of the logarithmic values between observed and estimated. The above step was applied to identify, in the statistical sense trough the estimation of the residuals, how realistic our PGA estimations are.
5. Results and discussion
We applied the EGFM to generate a lot of models using the Somerville [6] relations to characterize the parameters of the source. The total rupture area was 5,571.90 km^{2} with an average slip of 3.2 m and a combined area of asperities of 1,269.11 km^{2}, in which the area of largest asperity was of 812.33 km^{2}. The hipocentral distance to the center of the closest asperity was estimate at 18 km. It should be pointed out that the Somerville [6] relationships do not define an azimuth to specify the location of asperity, this mean the largest asperity can be located within an azimuth of between 0^{o} to 360^{o}. Then in the modeled process the position of largest asperity was varied from 0 to 360^{ o}.
The main goal of this study was in the context of obtaining the most probable scenario of the ground response in the major cities of region upon the occurrence of a M_{W} 7.3 earthquake generated in the area in study. In this work, all the simulations were carried out to obtain a statistical sample that may represent the most probable scenarios of the ground response at the studied sites. For that purpose, the fault parameters and rupture process (i.e., (i) the azimuth of the closest asperity to hypocenter, (ii) the rise time, (iii) the rupture velocity, and (iv) the location of SMGA within the fault plane) were varied in an iterative process in order to generate a statistical sample of the most probable ground response scenario.
According to Somerville [6] an earthquake with a moment magnitude of 7.3 should have around 2.4 asperities. Based on the above, our procedure was divided in two stages. In the first stage, we used 2 SMGA and 3 SMGA in the second stage. These SMGA were positioned inside of different points into the dislocation area adopting the Somerville [6] criteria. In order to evaluate which scenario is the most probable we compare the PGA obtained from each scenario to those predicted by two GMPE, from Ordaz [16] and from Young´s [17].
It should be pointed out that both GMPE were obtained considering earthquake data from subduction tectonic environments as follows: (i) Ordaz [16] used subduction earthquakes from the Mexican pacific coast, and (ii) Young’s [17] subduction earthquakes from around the world. The process of finding the best fit of the simulated PGA with respect to Ordaz [16], and Young’s [17], GMPEs was based on the smallest residual criteria between the simulated PGA and the estimated GMPEs. The above, allowed us to identify 2 different source models that provided PGA values that better matched with the used GMPE as follows: (i) the first one with 2 SMGA, and (ii) the second one with 3 SMGA.
In order to find a best residual in each interaction, we varied the position of SMGA, the rupture velocity, rise time, radiation pattern, and the size of SMGA changes in the strike or dip directions, according to the azimuth of the station or stations that had poor adjustment with GMPE. In the process of modeling we found little sensitivity of the synthetics to the rise time variations. On the other hand, we found high sensitivity of the synthetics to the rupture velocity variation, the size of the SMGA and its location inside of fault plane. The parameters with major weight in the modeled are the number, the size and the location of SMGA. The optimal model is a combination of all these parameters. The best model for each stage was determined by minimizing the residual between synthetic and observed PGA.
Figures 2a and 2b show the comparison between both GMPE versus our results for the lowest residual case of the source models (two SMGA). Figure 2a shows the Young’s [17] GMPE curves for rock and soil.
Our data showed three clusters that according to their hypocentral distances are distributed as follows: (i) The first group was distributed within the distance from 35 to 60 km, in this group 5 of the simulated PGA were located nearby of both curves; (ii) The second group is defined for distances range from 60 to 120 km, on this case the PGA are distributed almost evenly below the GMPE curves; (iii) Finally, the third group is defined for distances range from 120 to 500 km, on this particular situation the PGA values are over-estimated by the GMPE and show a clear tendency to attenuate faster than the pattern showed on the GMPE. Figure 2b shows the Ordaz [16] GMPE, the author uses thrust subduction earthquakes from Mexico (such as event simulated in this study). In general the comparison shows that 90% our results are located above this curve.
For the GMPE of Ordaz [16], we compare only the 19 stations seated on rock. For the GMPE of Young’s [17] we compare the 19 stations seated on rock, and 6 stations seated on soil sites (table 3), each of these groups with the respective curves for sites on rock and for sites on soil. From the modeling process of the target event, when three SMGA were used the lowest residuals we obtained between the PGA and the GMPE were: (i) 0.011 for Young’s [17] GMPE (Rock), (ii) 0.101 Young’s [17] GMPE (Soil), and (iii) 0.252 for Ordaz [16] GMPE for the sites shown in Figure 2c and 2d. The best fit between the estimated values of PGA and GMPE, estimated by means of the lowest residual, was obtained for the source model with two SMGA. The respective estimated residuals were: (i) 0.009 for Young’s [17]GMPE (Rock), (ii) 0.024 for Young’s [17] GMPE (Soil) and (iii) 0.245 for Ordaz [16].
Ramirez-Gaytán [23] simulate PGA and acceleration time histories for Tecoman earthquake located in the adjacent area at the event here simulated. The particularity and importance of this study is that the authors used observed strong motion records as a comparison reference to adjust the synthetics and found that from the comparison between PGA synthetics with respect to Young’s [17] GMPE for rock and soil curves, and the Ordaz [16] GMPE curve (Figures 2e and 2f) behave very similar to the description previously provided for the same ranges of distances obtained in this study.
We made an additional comparison, for Tecoman earthquake Singh [10] used records of stations located at distances larger than 50 km and comparing the PGA versus Ordaz [16]. In figure 2g it can be seen that the behavior is the same: observed PGAs lie above the curve. The comparison made are important because in the first case Ramirez-Gaytán [23] generate PGA based in a previous model Ramirez-Gaytán [24] when using observed records to compare with synthetics. In the second case Singh [10] use real data, in both cases the results are very similar with the obtained in this study. All mentioned studies correspond at the same tectonic environment. Figure 2g show than for Singh [10] major PGA locate at distances above 120 km this explained because he use data from regional stations. In the case of Ramirez-Gaytan [23] major PGA are located at distances from 10 to 100 km this explain why authors only use data from local networks.
In this study, we compare our PGA's with those obtained by the GMPE of Ordaz [16] and because they used data from Mexico. However, all of these studies considered all of the data available at the time of their analyses. This means that the data used were essentially strong motion data from the southern part of the country Tejeda- Jácome and Chávez-García [25]. None of these studies included data from western Mexico, along the northern section of the subduction zone. It is uncertain that ground-motion prediction equations developed for Guerrero in a very different tectonic setting can be applied to Colima Tejeda-Jácome and Chávez-García [25]. For this reason, it is important to compare our results with GMPE of other parts of the world, such as the GMPE of Young’s [17]. They used some similar parameters to those of the event simulated in this study (tectonic environment or hypocentral distance magnitude, etc.). This seems to be justified when we observed that the minor residual is reached when compared with the GMPE of Young’s [17] with minor residual (0.009) than Ordaz [16] with residual of 0.245 for the model with best adjust (model with two SMGA). However, we expect that any of the two GMPEs compared in this paper satisfy the detailed similarities with the event simulated in this study for the Colima region. For this reason, our intention is to use local and world parameters in order to validate our results and to give a degree of confidence when applying this methodology to future earthquakes for determination of acceleration time histories and PGA.
The main contribution of the process detailed above is not obtaining PGA, for this case is sufficient to consult GMPE. The purpose of this paper is to prepare acceleration time histories to be used by structural engineers on the analysis and design of structures. In modern seismic design approaches the quality of a structural solution frequently depends on the detailed knowledge that designers have on the characteristics of the seismic ground motion that the structure will suffer at a site if an earthquake of a given magnitude occur at a given location. Figures 3a-3d show the acceleration time histories for the model with minor residual after applying the process previously described. The largest acceleration was obtained in the near-to-source station SJAL (rock site) with 0.47 g. It is important to point out that stations CUP and SCT (figure 3c) located in Mexico City with hipocentral distances larger than 470 km still produced considerable values of peak accelerations (3.53 and 3.55 gal) similar of those of the GDLC station (with 1.71 and 5.3 gal) located to 273.48 km from the epicenter. This is due to the fact that seismic waves are enormously amplified at lake-bed sites respect to hill-zone sites in Mexico City, although it has been suggested that even hill-zone sites suffer amplification and Singh [26].
On the other hand, with acceleration time histories it is possible to generate a response spectrum, which considers forces related to parameters of maximum response like spectral acceleration. Response spectra are essential for the seismic design, any effort to accurately predict these should be done. The synthetic acceleration response spectra for an equivalent viscous damping of 5 percent were calculated and compared with the elastic acceleration design spectra for structures of group B (standard occupancy) of the Manual of Civil Structures MOC-2008 [27], a model design code in Mexico and seismic provisions for current Mexico’s Federal District Code NTCS-2004 [28], Sites CUP and SCT. In the MOC-2008 [27] code, seismic hazard in Mexico is defined as a continuum function where peak accelerations in rock are associated with return periods that were obtained using an optimization design criterion to define the seismic coefficients for the plateaus of the elastic design spectra for standard occupancy structures Tena-Colunga [29].
It can be also seen that only in 1 of 25 stations with distances comparable to the source dimensions the seismic ordinates are underestimated respect to the elastic design spectra. This station is SJAL located in rock site having seismic ordinate of 1.97 g in structural period of T= 0.3 sec. This is possibly because radiation patterns and source heterogeneity, which are taken into account in this work, causes spatial variations in ground motion around the fault Somerville [30]. These changes are very clear in the EGFM, but may be unnoticed in the GMPE. Therefore, obtaining future acceleration time histories and response spectra is a matter of essential interest for the present seismic design to achieve more efficient and safer structures.
6. Conclusions
We use the 25 records of August 13 2006 earthquake Mw 5.3 as element event to simulate strong ground motions for an eventual earthquake Mw 7.3 in the studied area. To reach this objective we integrate the advantages of three methodologies (EGFM, Somerville [6] relations and GMPE) to estimate the possible PGA, acceleration time histories and response spectra for an eventual earthquake in the studied area. We apply the empirical Green’s function method (EGFM) whose main contribution is to reflect a model that considers the source, the path, and the site effects. In Mexico this method has been used to simulate an event that already occurred. In this study we applied it to predict some probable earthquake which may be expected in the region. To overcome the absence of observed records we made use of Somerville [6] relations to be able to make more accurate predictions of strong motions and two GMPE adequate for region to compare our results. The process of finding the best adjustment generated 2 different models (2 and 3 SMGA). This process of minimizing the residual between synthetics and observed PGA clearly shows that the mean residual for 25 stations is obtained when comparing with GMPE of Young’s [17] and modeled with 2 SMGA. Ramirez-Gaytan [23] simulate PGA for Tecoman earthquake, whit difference that in this case the earthquake had occurred and exist observed records to compare and adjust synthetics, results are similar to those obtain in this study. Singh [10] comparing real PGA for Tecoman earthquake versus GMPE of Ordaz [16], results are similar to those obtain in this study. The purpose of this paper is to rescue the acceleration time histories of simulated event prepared to be used by structural engineers to analyze and design structures. Response spectrum show that for 1 of 25 stations (this station is near the source or with distance comparable with the source dimensions) the seismic ordinates are underestimated with the design spectra of the MOC-2008 [27] due possibly to radiation patterns and source heterogeneity, which is still to be confirmed by future records. For any practical evaluation of the seismic hazard in terms of response spectra is possible to integrate the advantages of three methodologies aboard in this study to estimate the possible PGA, acceleration time histories and response spectra for an eventual future earthquake.