Location and moment tensor inversion of small earthquakes using 3D Green’s functions in models with rugged topography: application to the Longmenshan fault zone

With dense seismic arrays and advanced imaging methods, regional three-dimensional (3D) Earth models have become more accurate. It is now increasingly feasible and advantageous to use a 3D Earth model to better locate earthquakes and invert their source mechanisms by ﬁtting synthetics to observed waveforms. In this study, we develop an approach to determine both the earthquake location and source mechanism from waveform information. The observed waveforms are ﬁltered in different frequency bands and separated into windows for the individual phases. Instead of picking the arrival times, the traveltime differences are measured by cross-correlation between synthetic waveforms based on the 3D Earth model and observed waveforms. The earthquake location is determined by minimizing the cross-correlation traveltime differences. We then ﬁx the horizontal location of the earthquake and perform a grid search in depth to determine the source mechanism at each point by ﬁtting the synthetic and observed waveforms. This new method is veriﬁed by a synthetic test with noise added to the synthetic waveforms and a realistic station distribution. We apply this method to a series of M W 3.4–5.6 earthquakes in the Longmenshan fault (LMSF) zone, a region with rugged topography between the eastern margin of the Tibetan plateau and the western part of the Sichuan basin. The results show that our solutions result in improved waveform ﬁts compared to the source parameters from the catalogs we used and the location can be better constrained than the amplitude-only approach. Furthermore, the source solutions with realistic topography provide a better ﬁt to the observed waveforms than those without the topography, indicating the need to take the topography into account in regions with rugged topography.

between synthetic and observed body and surface waveforms (e.g., Dziewonski et al. 1981;Nabelek 1984;Dreger and Helmberger 1991;Kikuchi and Kanamori 1991;Zhao and Helmberger 1994;Zhu and Helmberger 1996). To date, most earthquake source studies utilize 1D Earth models to calculate synthetic waveforms or Green's functions. To reduce the effects of the arrival time difference caused by 3D structural heterogeneities on source mechanism inversion, Zhao and Helmberger (1994) developed a ''cut and paste'' (CAP) method that breaks seismograms into segments of body and surface waves and allows time shift when fitting data. In places with highly heterogeneous Earth structures, this CAP approach may not fully account for the complexity of waveform phenomena due to the 3D structures. One of the major sources of bias in moment tensor inversion is phase skipping between the observed and model-predicted waveforms (Liu et al. 2004). When traveltime delays due to 3D velocity heterogeneities are comparable to the wave period, phase skipping becomes a challenging problem. Hjorleifsdottir and Ekstrom (2010) estimated the effects of the 3D Earth model on the global centroid moment tensor (GCMT) solutions using synthetic data and found only small errors on average in the scalar moment and tensor elements, while the centroid depths and times were biased. Another method to correct the 3D effect using 1D velocity models is the calibration event approach developed by Tan and Helmberger (2007) and Chu et al. (2009), which utilizes a ground true event to derive various corrections for locating events. This approach can improve the accuracy of moment tensor inversion especially for very small earthquakes (M \ 3.5).
With dense arrays and advanced imaging methods, regional 3D Earth models have become more accurate. It is now possible to better locate earthquakes and invert the source mechanisms with a 3D Earth model by fitting observed waveforms. Accurately determining earthquake locations and inverting source mechanisms with a 3D Earth model should also be an important component of full-3D waveform tomography, as the bias of the earthquake locations may lead to serious bias in velocity structure in local seismic tomography (Thurber 1992). Liu et al. (2004) developed a waveform inversion technique using synthetic waveforms with the spectral-element method and Fréchet derivatives to determine the source parameters of small to moderate earthquakes in southern California. Since the derivatives are determined numerically by differentiating synthetics with respect to the source parameters (six moment tensor components, longitude, latitude, depth, and the reference location), massive forward simulations are needed for each event inversion. If the initial location is far from the true location, the derivatives of the location parameters may not be accurate enough to represent the non-linear variation of waveforms as a function of the location. Zhao et al. (2006) developed a strain Green's tensor (SGT) source inversion technique based on 3D reference models by applying source-receiver reciprocity to reduce the number of calculations. In their study, the earthquake location is assumed to be known and only the source mechanism is inverted. Shen et al. (2015) adapted the CAP approach into the SGT method with a 3D velocity model by fitting multifrequency and time-shifted waveforms. The method includes performing grid search in the vicinity of the reference location to find the event location, and solving a linear inverse problem to obtain the source mechanism that provides the best waveform fit. Because of the time shift, the event location is actually determined only by waveform amplitude, with little or no traveltime constraints.
In this study, we develop an approach to determine both the earthquake location and source mechanism from waveform information. We filter the observed waveforms in different frequency bands and select the time windows of individual phases. Instead of picking arrival times, we measure the traveltime difference by cross-correlation between the synthetic waveforms based on the 3D Earth model and observed waveforms. The earthquake location is determined by minimizing the cross-correlation traveltime difference. We then fix the horizontal location of the earthquake and perform a grid search along the vertical direction and invert the source mechanism at each point by fitting the synthetic and observed waveforms in the CAP sense. With no time shift between the synthetic and observed waveforms in the location step, we assume phase skipping is not an issue. Thus, this approach is appropriate in places with well-determined 3D velocity models for the periods of the waves used in the source inversion. In the following, we describe the methodology, verification, and application in the Longmenshan fault (LMSF) zone.

Method
To utilize a 3D model in both location and moment inversion, one approach is to do a simultaneous inversion of the location and moment tensor. However, waveforms are strongly nonlinear to the source location. To reduce the difficulty of the nonlinear problem in a linearized iterative inversion, we first use the time delay measured by crosscorrelation between the observation and synthetics to locate the event, then utilize the amplitude information between the observation and synthetics to invert the moment tensor. In the second step, we also allow the event to move in depth to find the best-fit moment tensor and event depth.

Constructing synthetic waveform by receiver strain
Green's tensor We need a way to efficiently calculate the synthetic waveform from a trial source location with a trial moment tensor to receivers on a 3D velocity model to locate events and invert source moment tensors. One straightforward approach is to perform one numerical simulation for each source location and moment tensor, but it requires extensive computation to produce synthetic waveforms if there are many potential locations and moment variations. A more efficient way is to utilize the reciprocity property of the wavefield between the source and receiver to calculate the Green's function from a receiver location to all potential source locations (Eisner and Clayton 2001), reducing the number of simulations to be proportional to the number of receivers. This approach has been used in waveform tomography (e.g., Zhao et al. 2005;Chen et al. 2007) and source parameters inversion (e.g., Zhao et al. 2006;Lee et al. 2011) and we briefly summarize the formulation here. From the representation theorem (Aki and Richards 2002, Eq. 3.23), the displacement of a point moment source can be written as where Green's function G nj (r, t, r s ) relates a unit impulsive force at location r s with direction ê j to the displacement response at location r in direction ê n . M is the moment tensor and the superscript S denotes the spatial gradient operator on the source coordinates. Since M is symmetrical and Green's function is reciprocal, Eq. (1) can be reformatted using the Green's function from the receiver to the source (Zhao et al. 2006) u n ðr; t; r s Þ ¼ where q s i G jn (r s , t, r) is named as the receiver strain Green's tensor (SGT) in Zhao et al. (2006) and can be obtained from the tensor outputs of three numerical simulations (three single forces, one simulation per force) by applying a pseudo delta source time function (STF) at the receiver location. In the following, the synthetic waveforms are assumed to be calculated by Eq. (2) from a precalculated SGT database, which contains all the 3D and topographic effects of the velocity model.

Locating events by traveltime delays from waveform cross-correlation
The locations of seismic events are usually constrained by the traveltime. There are two approaches to retrieve the traveltime information contained in observed waveforms. The most intuitive and widely used method is onset-time picking, which relies on the high-frequency assumption that the wavelength of the seismic wave is much smaller than the size of heterogeneities in the velocity model. The second approach performs cross-correlation calculation between observed and synthetic waveforms and uses the time shift of the maximum cross-correlation value as the time delay measurement (Shearer 1997), which accounts for the finite-frequency wave phenomena and may take full 3D effects into account if the synthetic waveform is calculated with a 3D model. In this study, to fully utilize the 3D numerical Green's functions and take into account the effects of small-scale heterogeneity and surface topography (Wang et al. 2016), we adapt the waveform cross-correlation approach to use the Green's functions calculated by 3D numerical waveform simulation instead of 3D raytracing to fit the traveltime information between the observation and synthetics. The event location can be obtained by minimizing the cross-correlation traveltime difference between the observed and synthetic waveforms as where w i = w a w d w s is the weighting factor based on the azimuth distribution w a , epicentral distance w d and data signal-to-noise ratio (SNR) w s . t c is the centroid time shift to correct the error of the origin time of the event and the cross-correlation time shift s is defined as the time lag when the cross-correlation function reaches its maximum where N is the total number of the phase waveform windows, t 1 and t 2 are the starting and ending times of the phase window. u s and u o are the synthetics and observed data, respectively. We can derive the sensitivity of s in Eq.
(3) to the location perturbation (e.g., Klein 1994;Gillard et al. 1996;Shearer 1997), then perform a linearized Gauss-Newton inversion to iteratively update the event location. Alternatively, since the initial locations of the events provided by GCMT or the China Earthquake Networks Center Unified Catalog published by the China Earthquake Data Center (CEDC) should not be far from the true locations, we adapt a grid-search approach by comparing the objective function Eq.
(3) at the grid points around the initial locations to find the point with the best waveform fit as the event location. This approach is more straightforward and does not suffer from the local minimal problem. Thus, we apply it to determine the event location in this study.

Moment tensor inversion
Once the event location is determined, the components of the moment tensor of the event are linearly related to the observed waveform by the spatial gradient of the Green's function (Aki and Richards 2002) and can be obtained by a linear inversion that minimizes the misfit between the observation and synthetics (e.g., Liu et al. 2004). Another way to find the best-fit moment tensor is to decompose the moment tensor as a linear combination of six elementary moment tensors (Kikuchi and Kanamori 1991;Shen et al. 2015), where the 6 elementary moment tensors M m are and a m is the coefficient. The source mechanism can be written as Eq. (7) using the coefficients a m , M ¼ a 2 À a 5 þ a 6 a 1 a 4 a 1 a 2 À a 6 a 3 a 4 a 3 a 5 þ a 6 2 4 3 5 : Mathematically, using the six elementary moment tensors to fit the waveforms is equivalent to using the six single components of the moment tensor. The advantage of using the six elementary moment tensors over the single moment components is that each of the elementary moment tensors has direct physical meaning in the source mechanism. For example, if we exclude M 6 in the decomposition, the source is a pure-deviatoric moment tensor. If we want to constrain the event to be a pure double-couple source, we can pose an additional constraint of zero det M.
From Eq. (1), the basis displacement component n for an element moment tensor m at station j can be calculated from the station SGT database (Shen et al. 2015 where S(r s , t) is the STF. The displacement u n for a general moment tensor can be represented as a linear combination of the basis displacements, where N e is the number of the elementary moment tensors used.
To find the best fit combination of the elementary moment tensors, we define the least-square misfit function as where m = [a 1 ,…, a n ] is the inversion parameter, namely the coefficients of the elementary moment tensors, N the total number of the phase waveform windows, w i a weighting factor defined as Eq. (3). N r = P N i¼1 w i is used to normalize the misfit, t 1 and t 2 are the starting and ending times of the phase window, d i (t) is the observed seismogram. dt is the remaining time shift between the observation and synthetics that cannot be explained by the event location but is likely caused by errors in the velocity model and event origin time. The presence of this time shift indicates that the velocity model could be refined by a full 3D finite-frequency tomography (e.g., Zhao et al. 2005;Tromp et al. 2005).
For any given event location, the solution of Eq. (10) can be easily obtained by a linear inversion. As the source depth is affected by systematical errors in the velocity model and event origin time and sensitive to the source mechanism, we again adopt a grid-search approach, performing the moment tensor inversion at each depth while fixing the horizontal position. The final depth of the event and the moment tensor solution are determined by the minimal misfit value.

General work flow
The event location and moment tensor are inter-dependent in the above method. In the event location step, we need a preliminary mechanism to calculate synthetic waveforms, before we can use cross-correlation to fit the synthetics to observations to determine traveltime delays; while in the moment tensor inversion step, we need the event location estimation to invert the source mechanism. To address this inter-dependency, we adopt the following workflow. As Fig. 1 shows, we first invert the moment tensors at each grid point in the vicinity of the reference location provided by GCMT, CEDC, or other catalogs. Then, we estimate the scalar moment after doing an initial moment tensor inversion in each grid point without an event-specific STF. After we obtain the scalar moment, we estimate the source half duration based on the scalar moment of the event. The empirical relationship derived from the modeling of broadband moment-rate functions is (e.g., Ekström and Engdahl 1989;Ekström et al. 1992Ekström et al. , 2012: We then invert the moment tensor again with the eventspecific STF.

Synthetic tests and application to the LMSF zone
The LMSF zone is the most active part of the North-South Seismic Belt in China due to the collision between the Tibetan plateau on the west and the Sichuan basin on the east. Elevation rises from about 600 m in the southern Sichuan basin to over 6500 m at elevation peaks in a horizontal distance of less than 50 km. The regional topographic gradients typically exceed 10 %. Both the great Wenchuan earthquake of May 12, 2008 and the Lushan earthquake of April 20, 2013 happened in this region, which led to significant loss of property and lives. Accurate determination of event locations and source mechanisms is important for understanding the seismicity and studying the velocity structure to reveal the tectonics of this region. There have been many studies of the 3D velocity structures of the LMSF region (e.g., Pei et al. 2010;Liu et al. 2014;Yang et al. 2012). In this study, we use the velocity model of Liu et al. (2014) as our 3D reference model (Fig. 2). This model was obtained by a joint inversion of receiver function and surface wave dispersion from a dense seismic array. The original velocity model is given relative to the flat surface at sea level. In this study, we conform the velocity model to the surface topography, thus the thickness of each layer is unchanged, but the geometry of the modified layer has a topographic variation that follows the surface topography. If we directly use the source location and moment tensor solution from GCMT to simulate waveforms (e.g., event 20080226: Lat. 30.03°N, Long. 101.96°E, depth 29.2 km, M W 5.1) with the 3D reference model, the synthetics and observations do not show a high degree of agreement (see Fig. 3 for an example). Later, we will show that the locations and moment tensors re-determined using 3D SGT improve the waveform fit.
In order to construct the 3D SGT database in an area with steep terrains, such as the LMSF, we use a boundary conforming, curvilinear-grid finite-difference (FD) method (Zhang and Chen 2006;Zhang et al. 2012a, b, c;Zhang et al. 2014) to simulate wave propagation from the stations. The topography data are from GTOPO30. To increase the efficiency of the FD algorithm, we use a non-uniform grid to discretize the computational domain. The grid varies continuously with smaller spacing (1 km) in shallower depth near the free surface and larger spacing (1-2 km) at greater depth. The grid spacing is about 1 km in the north and east directions. The computational domain is surrounded by complex-frequency shifted perfect matched layers implemented through auxiliary differential equations (ADE CFS-PML; Zhang and Shen 2010). We select 21 permanent stations within the study area (Fig. 3) to construct the SGT database. The time step of the simulation is 0.389 s, the record time length is 150 s, and the disk storage of the SGT database is about 6 TB.

Numerical tests
We use the true station distribution in the LMSF zone and two hypothetical events at the same location with different source mechanisms to verify our method in locating events and inverting source mechanisms in the presence of strong topography and velocity heterogeneity.
The true locations of the hypothetical events are at (191, 240, -19.39) km with the origin of the coordinates at longitude 100°E (the X direction), latitude 28°N (the Y direction), and the sea level (the positive Z direction is upward). The first event has only a non-zero moment element M xy = 1.0, while the second event has a more general mechanism (Table 1). We use the curvilinear-grid FD method to generate the synthetic waveforms with the 3D velocity model and surface topography. We use nine seismic stations in inversion (Fig. 4) and add 20 % Gaussian noise in the waveforms of the second event. The initial locations of the two events are both set at (196, 243, -22) km, several kilometers away from the true locations. We use the P and surface wave segments of the three-component seismograms in two frequency bands (15-30 s; 7.5-15 s) in the inversion.
Since the LMSF zone has a complex structure and terrain, we set the search box width 8 km larger than the location error (about 5 km). From the inversion result shown in Table 1 and moment beachball comparison in Fig. 4, the source parameter inversion converges to the true solution. The misfit of depths is partly introduced by the discrete mesh. The polarities, peaks, and troughs of the waveforms fit well even in the case with the added noise (Fig. 5). Numerical tests indicate that our method can be

Application to the real seismic events in the LMSF zone
The LMSF zone has abundant and widely distributed seismic activities. The source mechanism solutions show diverse styles of thrust, normal, and strike-slip faulting (e.g., Chen et al. 1981;Pettersen and Doornbos 1987;Wu et al. 2004). We select 30 M W 3.4-5.6 events in the CEDC catalog (China Earthquake Networks Center Unified Catalog) that were recorded by the 21 broadband three-component stations in the study area from 2008 to 2015. We only use data with SNR [ 3 for inversion. The weighting factor in Eq. (10) is calculated by w i = w a w d w s . To account for various qualities of seismograms, we assign w s = 1 for time windows with SNR between 3-5 and w s = 2 for windows with SNR [ 5. The station azimuth weight w a is inversely proportional to the number of stations in 30°azimuthal bins to lessen the likelihood that waveforms from a cluster of closely spaced stations dominate the solution (Shen et al. 2015). The distance weight is inversely proportional to the epicentral distance. Each event source inversion involves at least seven stations.
Other inversion parameters are set as the same as that in the synthetic tests. Since the data quality varies from event to event, the number of stations left out of the inversion for the validation purpose is different for different events. Among the 30 selected events, the four largest ones have the GCMT solutions (Table 2). Our solutions and the GCMT solutions have similar fault planes and auxiliary planes in the beachball plots (Fig. 6). Even for the event with a poor station azimuth coverage like event 20110410 ( Fig. 6; Table 2), our moment tensor inversion is still robust. For the stations not involved in the source inversion, the synthetic waveforms with the source parameters from our FD SGT method fit the observed better than those of the GCMT solutions (Fig. 7). We note that there are significant differences in the CEDC and GCMT locations, especially in depth. For the four events in Table 2, the depth difference of some events is larger than 10 km.
To demonstrate the effects of the event location by minimizing the cross-correlation traveltime difference, we check the results of the event 20100427 and other five events. From Fig. 8, we can see that the cross-correlation traveltime difference from our approach is smaller than that of the amplitude-only approach. The smaller crosscorrelation traveltime difference can be important if we use the source locations with a full-3D waveform tomography method to alternatively refine the velocity structure and source parameters. The depth determined by the amplitude-only method is significantly shallower (Table 3) and our depths are more close to the result provided by CEDC. While the hypocenters of the CEDC solutions are not ground truth, we note that the 1D velocity model used in CEDC is optimized for the region and the true hypocenters and centroid locations for small events should be close.
Most source inversion studies adopt a flat topography to reduce the difficulty of calculating SGT. This assumption is   Also included are the CEDC locations and magnitudes than that in the case of a flat surface at sea level. The misfit as a function of depth is more convergent and less fluctuant in the case with the realistic topography. The synthetic seismograms based on the source solution with the realistic topography provide a better fit to the observed waveforms especially on the vertical component (Fig. 12). For events 20100427 and 20110410, the focal mechanisms and depths in the case with the realistic topography agree well to that of a flat topography. We note that the frequency band used in this study is not very high due to the consideration of the computational cost and the resolution of the reference model. The wavelength of the surface waves in this study is roughly around 20-100 km. At higher frequencies, the effects of the topography may be more significant.
To demonstrate the effects of the model uncertainty on source inversion, we systematically increase the P-velocity model by 5 % and invert the four largest events in Table 2 based on this new modified model. Figure 13 shows the inverted results of FD SGT method using the reference model (red-and-white beachball), modified reference model (blue-and-white beachball), and GCMT (black-andwhite beachball), respectively. The depths of the four events are shallower based on the modified reference model than those of the reference model while horizontal location and fault planes and auxiliary planes in the beachball plots are similar. Although the source mechanisms and event locations by our method do not significantly change with small model perturbation, the model uncertainty indeed affects the solutions, which indicates that we may use a waveform-based tomography to refine both the 3D Earth structure and earthquake source parameters to improve the data fit. In this study, we combine the velocity model of Liu et al. (2014) with the realistic topography into a 3D Earth model to locate events and invert for source mechanisms using waveform information of 30 small to moderate earthquakes in the LMSF zone. We use the curvilinear-grid FD method to accurately simulate seismic wave propagation. To increase computational efficiency, we use the 3D SGT database technique to calculate the waveforms from potential source locations to receivers. The event location is first determined by minimizing the cross-correlation traveltime difference between the synthetic and observed waveforms, then further refined by grid search along the vertical direction with source mechanism inversion at each point. Numerical synthetic tests verify that this approach can recover the true location and mechanism with 20 % added noise in the synthetic data. The re-determined locations of the 30 events are generally consistent with the CEDC hypocenter locations and the moment tensors of the four largest events generally agree with the GCMT Fig. 8 Histogram of cross-correlation time shifts for event 20100427. a The result using the method of Shen et al. (2015) and b the result using our method solutions. We note that the synthetic waveforms from our FD SGT method fit the observed better than those of GCMT.
We analyze the effects of determining event locations by cross-correlation traveltime difference. Compared to Shen et al. (2015)'s amplitude-only approach, our depths are closer to the results provided by CEDC and the crosscorrelation traveltime difference misfits are smaller, indicating the positive effect of minimizing the crosscorrelation traveltime difference. We also investigate the effect of the surface topography on source inversion. Most source inversion studies adopt a flat topography to reduce the difficulty of calculating SGT. This assumption may be appropriate for regions with smooth terrains but is not suitable for the LMSF zone. By comparing the inversion results in models with and without the rugged topography of the study region, we find that the source solutions in the case with the topography provide a better fit to the observed waveforms, especially on the vertical component. The frequency band used in this study is not very high due to the considerations of the computational cost and the resolution of the reference model. At higher frequencies, the effects of topography may be more pronounced and should be more important in forward simulation.
Our method can provide source mechanisms and event locations of small to moderate earthquakes, which can be further used in waveform-based tomography to refine the 3D Earth structure. Both earthquake source parameters and velocity structures affect waveform traveltimes and should be simultaneously or alternately updated in waveformbased tomography. Our method could be an essential component of a full-3D waveform tomography workflow to provide earthquake locations and mechanisms toward minimizing the cross-correlation traveltime and amplitude misfits.