Numerical Modeling for Earthquake Source Imaging: Implications for Array Design in Determining the Rupture Process

The aim of this study is to develop a numerical method to evaluate the optimum seismic station distribution for imaging the source rupture process of an earthquake. Based on the beam-forming technique, the source rupture distribution of an earthquake can be reconstructed through theoretical travel time correction and waveform stacking. Numerical tests show that this method successfully reconstructs the main displacement distribution from an assumed source rupture plane. In accordance with assumed fault models, seismic waveforms are numerically generated as input data for further source imaging. From synthetic seismograms, we reconstructed the rupture distribution of these assume earthquake sources, and analyzed error systematically. Results of this study indicate that receiver distribution types really affect the successful reconstruction of the slip distribution. Furthermore, parameters such as dip angles and frequency content also play important roles in reconstructing earthquake sources. The proposed method is simple, inverse efficiently, and no initial condition required. Further applications of this method are suggested to image source rupture from near field strong motion observations and to design seismic array to effectively observe seismic rupture properties.

earthquake.Recent examples are the 1995 Kobe, Japan, earthquake and the 1999 Chi-Chi, Taiwan, earthquake, both induced several thousand human casualties and significant property losses near their source areas.Taiwan is located in the circum-Pacific earthquake zone with a high population density.Earthquake source studies are always important issues in Taiwan to reduce seismic hazards.The study of the earthquake source rupture process can provide us with not only the magnitude and direction of regional tectonic stress, but also the dynamic rupture characteristics of earthquakes.Furthermore, many source physical properties can be determined from the analysis of an earthquake rupture process.In order to observe detailed rupture processes, seismometers are usually deployed as near earthquake faults as possible.However, how to best use limited resources to efficiently image sources, to timely report observations, and to lessen the damage from earthquakes are challenges to seismologists.
Some studies have considered this issue and reported the modeling or observation for earthquake rupture processes using seismic arrays.Iida et al. (1988) used theoretical waveform inversion analysis to give suggestions about the optimum strong-motion array geometry for source inversion.The dynamic source processes of the 1968 Tokach-Oki and the 1969 Kurile Islands earthquakes were investigated by Schwartz and Ruff (1985) based on source time functions deconvolved from teleseismic P waves.Fukuyama (1991) studied the rupture processes of the 1987 east Chiba earthquake using aftershocks as empirical Green's functions.Fukuyama and Mikumo (1993) used the waveform inversion method to determine the distribution of rupture starting times and slip distribution.In their approach, the slip functions were calculated by the initial crack model, or the use of previous crack inversion analysis.The rupture process of the 1989 Macquarie Ridge earthquake was studied by Tada et al. (1993) using a broadband array in Japan.Because the epicentral distances of those stations were near 100 degrees, the plane wave and point source were assumed and the beam-forming technique (Capon 1973) was employed to analyse the source rupture direction at teleseismic range.Das and Suhadolc 1996 used the standard large range inverse method to invert rupture distributions, and performed numerical experiments using artificial data to determine rupture and moment radiated processes.The slip distribution and rupture history of the 1992 Landers, California, earthquake were imaged by Cohee and Beroza (1994) using near-source displacement seismograms.Dreger et al. (1995) used an empirical Green's function inverse approach to estimate kinematic source parameters of the 1993 Klamath Falls earthquake.Hartog and Schwartz (1996) combined body wave and surface wave data, and used an empirical Green's function method to analyze source time functions.Ide and Takeo (1996) reconstructed the dynamic rupture model of the 1993 Kushiro-Oki earthquake using strong motion data.
The 1999 Chi-Chi earthquake was the largest inland earthquake in the past 100 years in Taiwan.This disastrous earthquake was well recorded by more than 400 free field digital accelerometers operated by the Central Weather Bureau (CWB) which has monitored earthquake activities in Taiwan since 1990 (Shin. 1993).The data set represents the most complete strong-motion records to study the disastrous earthquake's source rupture processes and ground motions.Many important discoveries about this centurial earthquake have been studied and reported.The location and source type have been examined by using seismic data (Shin et al. 2000;Wang et al. 2000), waveform inversions (Kikuchi et al. 2000;Lee and Ma, 2000;Ma et al. 2000), and forward simulations (Dalguer et al. 2001;Huang 2000;Huang et al. 2000;Oglesby and Day 2001).Although, it is realized that the seismic data from the 1999 Chi-Chi, Taiwan, earthquake provided a unique opportunity to study the earthquake source properties based on those strong motion records, however, till now, only few reports discussed the source properties based on seismic array recordings (Huang 2001).In this study, we developed a numerical code based on the beam-forming technique and designed a numerical experiment to demonstrate a method of finding the optimum receivers distribution of a seismic array in order to determine source rupture processes.Employed this method, existing or future seismic arrays can be modified or designed to monitor seismically active faults that have high probabilities of inducing large earthquakes.In addition, previous information on crustal models, fault geometries, and local seismicity can be easily integrated into the design of seismic arrays.

METHOD
A seismic array, in this study, is a set of seismometers with a common timing system that is distributed over an area of the earth's surface to record seismic signals.Complex seismic wave propagation within the seismic array can be easily identified.Such an array is useful for studying detailed characteristics of wave propagation across the array.Usually, the seismic array is used to enhance the signal to noise ratio of weak seismic waves.However, seismic array data can be used to determine propagation directions and velocities of seismic waves within an array.The beam-forming technique (Capon 1973;Aki and Richards 1980) is a simple but most effective method to process array data and was used in this study.Before introducting the method used to image source rupture processes, we describe the forward computation method used in this study to generate synthetic seismograms.
Generally, a rectangular rupture plane can be considered as a summation of point sources located on individual grid of fault plane (Fig. 1).The observed ground motion U(t) at a station S from a rupture fault plane (shaded area of Fig. 1) can be computed as: , where is the impulse response of a grid point source (each black dot on the fault plane of Fig. 1).and are the weighting factor and time delay of grid point , respectively.In this study, i and j express the grid indexes of horizontal and vertical directions, respectively.The time delay was computed based on the travel time between each and station S, and accounted for the rupture time from the initial point to each .The weighting was considered as slip amount of each .The travel time can be computed by numerical methods (Boore 1972;Huang 1992) or dynamic raytracing (Cerveny 1985) for very complex earth models.However, in this study, we focused our interest on the rupture properties of earthquake faults, and considered the earth model as a simplified half-space model and straightline ray path.The "observed data" employed by numerical experiments of this study were actually generated by synthetic seismograms using above method.
Source imaging for earthquake rupture processes is an inverse problem.The basic assumption used to invert the source rupture process is that the raypath between each point of the rupture plane and any station is different when seismic station is very close to the hypocenter area.It implies that the travel times of seismic waves propagated from the rupture plane to the array stations are different.Furthermore, there is a time delay to radiate seismic waves from different portions of a fault plane when a seismic fault ruptures with finite velocities.After correction for theoretical travel time, all seismograms can be used to solve the energy distribution by waveform stacking.In this study, we used a beam-forming technique to stack array observations to reconstruct the source slip distributions.According to the beam-forming approach, a time delay is inserted in a seismogram U(t).Then, a set of seismograms which have time delays of each jth station are stacked.The stacked seismogram B(t) is called the beam and is given by: where N is the number of stations.According to this approach, the uncorrelated noise may be reduced to 1/N of original signals.To image the source using beam-forming, we used the concept shown in equation ( 2) and define the beam value for any point with spatial coordinate (x, y) on an assumed rupture plane as: where N is the number of stations to be used, i is the time interval to be stacked and j is the station index, is the seismic wave arrival time from a point source with spatial coordinate (x,y) to station j, is the observed seismogram of station j which is shifted by time , is the weighting factor of station j which is usually assumed to vary between 0 and 1 to define the data quality.Herein, the seismic wave arrival time (travel time to be shifted) is calculated from point (x,y) to station j using a suitable earth model and rupture time delay from initial rupture point to point (x,y).Employing equation (3) to image the earthquake source, we must calculate the beam value of all grid points of a source area according to seismograms recorded by each station.Finally, can be considered as the inverted slip distribution of fault plane.

ANALYSIS AND RESULTS
To describe the analytical approach used, both forward and inversion computations from a simple case shown in Fig. 1 are presented.For illustration, a vertical rupture plane with a source area of 20 kilometer square was used.The fault plane was considered as a pure dip-slip fault and its slip distribution was assigned as two major asperities as shown in Fig. 2. The source rupture type was considered as a circular fault with a constant rupture velocity of 3.2 radiating away from the initial rupture point (point A of Fig. 2).The seismograms used to image the source rupture were generated using the procedure described above (equation 1).The seismic array used to image source rupture included 25 stations.These stations were distributed within an area 10 kilometers square above the seismic fault as shown in Fig. 3.A half-space earth model with seismic velocity of 4 was assumed over the source area.In order to obtain seismograms from a finite source rupture, we divided the continuous rupture plane into grids and simulated the rupture source by the summation of impulse response of individual point source.In this study, the rupture plane was divided as 441 uniformly square grids.The coordinate of every grid point was labelled, and slip displacement was assigned from the contour value in Fig. 2. The point source response was simply considered as a Ricker wavelet with time delay determined from its hypocentral distance divided by seismic velocity of model.Synthetic seismograms can be obtained after suitable time shifting and waveform stacking.The method to compute synthetic seismogram can be considered as an analogy to that using Isochrone (Spudich and Frazer 1984).To generate synthetic seismograms for stations, the Ricker wavelet was designed as having major frequency content of 5 Hz.The synthetic array waveforms are shown in Fig. 4. The seismogram of each station is very complex.Based on equation (3), synthetic seismograms in Fig. 4 were used to image finite fault slip distribution using the same rupture geometry and type as forward modeling.According to the theoretical travel time correction and waveform stacking of these point responses, the energy released from the rupture plane was reconstructed as shown in Fig. 5.The major features of the assumed rupture (Fig. 2) were reconstructed.However, since only finite observations were employed in the source imaging, the inverse results show slightly different results from the original assumed displacement distribution.Of course, the degree of error in inverse results were also affected by other factors, such as fault dip angle, rupture velocity, slip type and background noise.In order to classify the error sensitivity of different source parameters, we will design a numerical experiment to verify the sources of errors.
To assess the errors quantitatively, we need to compare the inverse results and slip distribution provided by forward modeling.First, the slip distributions both forward and inverse were normalized individually.Then, the error of each grid point was determined by the squared value of the difference between the forward and inverse slip displacement.Finally, the single value of error of an assumed model was defined as the normalized summation of errors of all  Fig. 7 indicate that when seismic wave frequencies are lower than 2 Hz, the inverse results affected by source rupture plane dips are larger than the array types.However, when frequencies higher than 2Hz, there is less influence from the angle of the source rupture plane dip, and the array configuration plays an important role.
In order to understand the influence of background noise, we selected the case of frequency 5 Hz and source rupture plane dip angle 67.5° from array type 9 for further study (Fig. 8).This special case has the smallest error value among all analyzed cases in this study   (Fig. 7).Two noises, the 30% peak amplitude random noise, and the 5Hz harmonic noise with a 50% peak amplitude were added to the original noise-free synthetic seismograms.The results are shown in Fig. 9. Seismograms shown in Fig. 9 were employed to stack the source image.The effect of background noise to source imaging was compared in Fig. 10.Figures 10a, b represent the inverse slip distributions with and without background noise, respectively.It is found that the error value in Fig. 10a was only 192% of that in Fig. 10b.The background noise has reduced the sharpness of the slip distribution pattern.However, the source images of Fig. 10a still resolved the source slip distribution acceptable.This indicates that when noise ratio was less than 30%, or there is noticeable noise of a different frequency content, such noise has less effect on the source imaging.

DISCUSSIONS
The source imaging method proposed by this study is a direct and simple method.Comparing results from inverse and forward approaches, it was found that there were a few differences between them, and the main slip distribution could be reconstructed accurately.If the source rupture process can be presented with high resolution by analyzing the data from a limited station array, the estimated slip distribution can be employed to determine the pattern of the source rupture.Furthermore, the same method can be used to design seismic arrays to observe source rupture.One of these cases is applied to an area where large earthquakes may occur.Typically, information about seismicity and fault geometry has previously been collected in high-risk active fault zones.This information can be employed to set the fault models of candidate rupture type and slip distribution.We can integrate existing geological information, fault geometries, and seismicity information to determine the optimum seismic station distribution from simulating the different array pattern in the fault zone.Thus, the method proposed in this study can be used to design a seismic array for a specific fault after a dense numerical survey.Numerical modeling of cases of possible source rupture types and fault geometries, as described in this study, will assist in designing array patterns that allow for obtaining the optimum resolution of potential rupture processes, and avoid the loss of important information from less suitable patterns.Furthermore, this method can be used to examine existing strongmotion arrays to test their capability to resolve earthquake source rupture processes.
The forward modeling approach used in this study was suitable only for case of seismic stations located very close to the source area, thus, the earth model is assumed as a uniform medium and the ray path is considered as a straight line.Of course, there are not the actual conditions for most seismic arrays which are established in regions of complex geology and surface topography.However, due to the complexity of earth's structure at source regions, the point source impulse response used in this study can be replaced by Green's functions for layer media response (Helmberger 1983;Wang and Herrmann 1980;Kennett 1983) or lateral inhomogeneous media (Boore 1972;Smith 1975;Huang 1992) to image the source rupture.Furthermore, it is possible to apply this approach to real cases using small event seismograms as empirical Green's functions (Hartzell 1978).
When comparing the assumed fault models, in some cases, there were a few errors (from 0 to 2 km) in estimating the locations of asperity centers.These errors were reasonable.First, the seismograms were digitized into discrete time series (the sampling rate, in this study, was 100 per second), and the assumed seismic velocity of the crust model was 4 ; i.e., there was 40 meters of uncertainty on the seismograms.We could have increased the sampling rate to reduce those errors, however, this would have required more computation time and computer memory to image source.Second, the wavelengths of seismograms according to different frequencies used in this study were 0.8 km, 1.6 km and 4 km, respectively.Stacking errors should be of the same order as wavelengths.Results of this study have indicated that the seismic images are less affected by added background noise.The possible reasons are: 1.The frequency of random noise used in this study was higher than seismic waves radiated from seismic source.2. The frequency content of 5 Hz harmonic noise was different from the seismic waves with a Ricker wavelet of 5 Hz major frequencies.The seismic waves with Ricker wavelet played the dominator role in influencing the imaged results.Actually, imaged results will be affected by noise; and its important to note that the argument regarding noise testing in this study were made for some special cases only.

CONCLUSIONS
A numerical method based on the beam-forming technique was developed to study the earthquake source rupture processes.Numerical experiments demonstrated that the slip distribution of an earthquake rupture plane can be successfully reconstructed based on the developed method.Results of this study showed that array configuration, seismic frequency content, and the dip angles of fault plane are the main parameters affecting the resolution of source imaging.Applications of this method were suggested to design seismic arrays to observe source rupture and to image source rupture processes with slight modification of Green's functions, which were for this study, simplified as Ricker wavelets.

Fig. 1 .
Fig. 1.Schematic map showing the relationship between the finite fault rupture plane and a near source seismic station (symbol solid triangular and noted as S).The fault plane is divided by equal space grids.i and j express the grid indexes of horizontal and vertical directions, respectively.(or other black dots on the fault plane) indicate the simplified point source for individual grid area.

Fig. 2 .
Fig. 2. The seismic fault plane with area of 20 km x 20 km.The contour values represent the slip displacement.Symbol A represents the initial point of source rupture.

Fig. 3 .
Fig. 3. Station distribution of a seismic array analyzed in this study.The number next to each station (symbol solid triangular) shows the station number.Those stations are distributed on 10 km x 10 km area.

Fig. 4 .
Fig. 4. Synthetic seismograms computation based on a vertical dip-slip fault rupture for stations of seismic array as shown in Fig. 3.Each seismogram was normalized to its maximum amplitude.

Fig. 6 .
Fig. 6.Seismic array maps of other nine types analyzed in this study.

Fig. 7 .
Fig. 7. Summarized error values with respect to frequencies, dip angles and array types from numerical experiments.

Fig. 8 .
Fig. 8. Schematic map showing the spatial relationship between earthquake source area and seismic array.

Fig. 10 .
Fig. 10.(a) The reconstructed source image of fault slip using synthetic seismograms of Fig. 9 which has artificial noise added.(b) The reconstructed source image for the same fault model of Fig. 10a using noise-free synthetic seismograms.