Source model of the 11th July 2004 Zhongba earthquake revealed from the joint inversion of InSAR and seismological data

1.   Introduction

The 11th July 2004 Zhongba earthquake occurred about 100 km to the north of Yarlung Zangbo suture belt, within the Lhasa block which is located in south of Tibet. Generally speaking, the collision of several continental fragments accreted to the southern margin of Eurasia led to the thick crust of the Tibetan Plateau. A widely held view is that normal faulting on the plateau, which extracts potential energy from the crust beneath the plateau, started from the time when maximum elevation of the plateau had been reached (Blisniuk et al., 2001; Molnar and Tapponnier, 1978). There are numerous NNE-SSW or near NS oriented valleys extended from Himalaya to Lhasa block. Many normal faults, on one of which this earthquake took place, are active in this region (Blisniuk et al., 2001). GPS observations across the plateau estimate ENE-WSW stretching of 21.6±2.5 mm·a^-1 between 79ºE and 93ºE across the central plateau, with 10–14 mm·a^-1 NNE-SSW shortening occurring between north of the Himalaya and south of the Altun and Kunlun faults (Zhang et al., 2004). Indentification of the real fault plane and slip history of this earthquake can provide useful information on mechanics of lithosphere (geodynamic processes) in Tibetan Plateau, as well as the stress field. This earthquake is useful for understanding seismotectonics of the normal faulting mechanisms, including issues such as centroid depth of the event, rupture velocity, etc. Rupture velocity appears to be an indicator of fault maturity, for example, the high rupture velocity of the 2009 Inglewood event is compatible with mature fault of Newport-Inglewood fault (Luo et al., 2010).

M 6.0–6.5 earthquakes happen much more frequently than larger events (M> 7.0) and can excite strong ground motion that causes economy and life loss. But details of the source process are relatively difficult to determine just from teleseismic waveform as our limited knowledge of 3D Earth structure which could not be ignored when studying the seismograms at higher frequency. As seismometers are more sensitive to the temporal than to the spatial distribution of slip, there are trade-offs among model parameters (namely between slip, rupture velocity and rise time) in source models derived from seismological records. This is especially true if no, or poorly distributed, near-source strong-motion data are available. Fortunately, geodetic data can provide more constraint on the slip distribution of this kind of earthquake. Recently developed geodetic and remote sensing techniques, such as InSAR (interferometric synthetic aperture radar), GPS (Global Position System) and optical images could provide extra information of the static ground deformation (Ji et al., 2001, 2002a, b). But static data for earthquake with magnitude around M 6, such as InSAR, alone cannot determine the area of the fault plane independent of magnitude of slip. Thus, the kinematic slip model might be better constrained from complementary data sets of seismological waveform records, field measurement of fault slip and measurements of static deformation using geodetic and remote sensing techniques.

Accurate fault geometry is important for finite fault inversion as it not only affects radiation pattern of seismic wave-fields but also the static displacement on the free surface (Konca et al., 2010). The mechanism of this event was reported by different agencies as listed in Table 1, where we can see diversity of the fault orientation, especially in strike and depth. Thus, it is necessary to resolve the ambiguity of fault geometry before doing finite fault inversion. The regional broadband seismic data are first used to determine the point source mechanism, which is 171º/42º/-83º/6.16/11 (strike/dip/rake/slip/M_W/depth). Then we fixed the fault plane and allowed the earthquake to move around so as to determine its best location by fitting the InSAR data. Finally, we use this fault geometry to obtain a distributed source model by the joint inversion of InSAR and teleseismic body waves.

Table  1.  Focal mechanism solutions of the 11th July 2004 Zhongba M_W6.2 earthquake

 | Show Table

DownLoad: CSV

2.   InSAR and seismic waveform data

Two descending SAR data sets required by ENVISAT satellite of the European Space Agency (ESA) on 17th March and 8th September 2004 have a small baseline and can fully cover the deformation field caused by the 2004 Zhongba earthquake. We applied ROI PAC InSAR software and DORIS orbit data to produce a deformation interferogram, and removed the topographic signal from the interferogram using a digital elevation model with ∼90 m resolution generated from National Aeronautics and Space Administration (NASA)'s Shuttle Radar Topography Mission (Farr et al., 2000).

At regional distance, this earthquake was well recorded by a temporary seismic array from PASSCAL (Program for Array Seismic Studies of the Continental Lithosphere), which consists of 68 broadband stations (Figure 1). The near north-south trend profile covers the epicentral distance from 130 km to 360 km with the maximum azimuth gap of 220º. Seismic data at these distance ranges are usually dominated by Pnl waves (a P-wave train), direct (Sg or Sn) and Mohoreflected shear waves as well as large-amplitude surface waves. Pnl waves contain depth phases such as sPn, sPg, therefore they are sensitive to the depth of the earthquake, and we enhance the weight for Pnl waves relative to surface waves using the cut-and-paste (CAP) method (Zhao and Helmberger, 1994). Because the absolute amplitude waveform data are used instead of normalized ones, the improved CAP method can better take advantage of relative ratio of different phases and avoid singularities in the source parameter space at those points where source orientation generates nodal synthetics (Zhu and Helmberger, 1996).

Figure  1.  Regional (left) and teleseismic (right) broadband stations (triangles) used in this study for the July 11th 2004 Zhongba earthquake in western Tibet of China, the stars indicate the epicenter and the solid square is the Zhongba County. Observation of stations within the gray rectangle are plotted in Figure 2. The concentric circles in the right panel are shown every 30° of azimuthal distance from the epicenter.

DownLoad: Full-Size Img PowerPoint

Attempts to fit Pnl and surface waves in absolute time are difficult since the synthetic seismograms are computed with an approximated 1D crustal model. Therefore we allow different time shifts to align Pnl waves and surface waves respectively; we also filter the data and the synthetics so that they only reflect the long wavelength smoothed structure of the crust. Pnl waveforms are filtered to include higher frequencies (0.02–0.2 Hz) than shear waves and surface waves (0.01–0.1 Hz) to better identify the depth phases that constitute Pnl. We compute the synthetics using a frequency-wavenumber double integration (Zhu and Rivera, 2002) in a layered elastic crustal model (see Table 2).

Table  2.  Layered crustal velocity model of Tibetan Plateau (from Zhu et al., 2006)

 | Show Table

DownLoad: CSV

Teleseismic P and SH waves are usually used for finite fault inversion of large earthquakes because they travel through relatively simple lower mantle and are less affected by shallow structure and triplications. We collected 23 teleseismic P-waves with good azimuthal coverage (Figure 1) for finite fault inversion. As this was a normal event, radiation pattern of SH waves are relatively weak at teleseismic distance, thus we do not have good SNR (signal to noise ratio) records for SH waves. Instrument responses are removed from original data and converted to displacement. A band-passed filter (0.002–1.0 Hz) is applied for these records after manually picked onset of P and SH arrivals. Higher corner frequency ( > 1.0 Hz) in bandpass-filtering teleseismic P waveforms is not helpful because teleseismic P wave amplitude at frequency above 1 Hz can be different from theoretical prediction by a factor of 10 (Ni et al., 2010; Ruff and Helmberger, 1982).

3.   Inversion and results

3.1   Point source mechanism

Firstly, we used the CAP technique on the regional broadband waveform records to obtain a point source mechanism. This method applies a direct grid search through all possible solutions to find the global minimum of misfits between the observations and the synthetics, allowing time shifts between portion of seismograms and the synthetics. The synthetic displacement for a double couple-source could be written as

here, i=1, 2, 3 corresponds to three fundamental faults, i.e., vertical strike-slip, vertical dip-slip, and 45º dipslip. G_i's are the Green's functions, A_i's are the radiation coefficients, and φ is the station azimuth. M₀ is scalar moment; θ, δ and λ are strike, dip and rake of the source that we want to determine from data. They are estimated by fitting the data in L₂ norm manner, with time shifts allowed between synthetic and data to get the maximum cross-correlation coefficients (Zhu and Helmberger, 1996), as shown in Figure 2.

Figure  2.  Regional displacement of data (black) and synthetic (red) calculated by the best point source solution (see Figure 3 for details). Waveforms at each station are segmented into Pnl (vertical and radial) waves and surface waves (vertical, radial and tangential). Station name is displayed at the left side of each wave train, with epicentral distance above (in km) and azimuth below (in degree). The first number below every waveform pair is the time shift necessary to align them up and the second number is cross-correlation coefficient in percentage. Pnl waves are filtered between 5 s and 50 s and surface waves are filtered between 10 s and 100 s. Here we show only a portion of fits to data, see the supplement for all waveform fits.

DownLoad: Full-Size Img PowerPoint

We obtained a strike of 171º/342º, dip of 42º/48º, rake of -83º/-97º, moment magnitude of 6.2 and a depth of about 11 km (Figure 3), which is almost a pure dip-slip (normal) event. Here we displayed waveform fits for only a portion of representative stations, fits for all the stations could be found in Figure A1. The synthetics generated by the preferred point source mechanism fit the data pretty well for almost all stations, both for Pnl waves and surface waves. By taking advantages of relative amplitude ratio between different phases, such as Rayleigh wave/Love waves, Pnl/surface waves, CAP can usually produce stable and robust mechanism even for this case in which the station coverage has a maximum azimuth gap of 220º. As stations are densely distributed along a near northsouth line and have azimuth coverage of about 140º, we can observe some detail changes of radiation pattern for Love waves (Figure 4). We observed a nodal plane along azimuth of about 75º, where the data show exactly the same feature (station H1240 and H1230 in Figure 2). The azimuth of this Love wave-nodal plane is most sensitive to the strike of the fault. The Love wave radiation pattern with strike 10º different from preferred mechanism is shown in Figure 4, from which we can see clearly the nodal plane is no longer along the azimuth of station H1240 as shown in the data.

Figure  3.  Depth sensitivity of the cut-and-paste (CAP) inversion for the July 11th 2004 Zhongba earthquake. The errors are scaled by the minimum error at the best depth, and the point source mechanism we obtain is 171°/42°/−83°/6.16/11 km (strike/dip/rake/moment/magnitude/depth).

DownLoad: Full-Size Img PowerPoint

Figure  A1.  Regional displacement of data (black line) and synthetic (thick gray line) calculated by the best point source solution, the dashed lines denote the components not used in the inversion. See Figure 2 for detail of description.

DownLoad: Full-Size Img PowerPoint

3.2   Joint inversion of InSAR and teleseismic P-waves

Determining the actual fault plane for this M6- class earthquake is not easy because the earthquake did not rupture on the free surface. Fortunately, we have a descending InSAR image for this event, which shows elongated bulls-eye pattern with maximum lineof-sight displacement of 34 cm corresponding to hanging wall subsidence, corresponding to the fault plane with strike of 171º and dip of 42º (Figure 4). Furthermore, to reduce the uncertainty of earthquake location, we searched for the location of the fault plane by moving it east-westward. The fault plane length along strike is set long enough to make sure all possible slip will not reach the edge of the fault plane. At each longitude grid point, a static inversion is done to find the best slip model. The scaled misfits are plotted versus longitude of epicenter (Figure 5) when the depth of hypocenter is fixed at 11 km, and a minimum error is found at longitude of 83.73º.

Figure  4.  Normalized radiation pattern of direct SH waves for strike of 171° (left), 161° (middle) and 181° (right), while dip (42°) and rake (−83°) angles are fixed. White triangles are regional broadband stations, the small black box indicates the stations used in Figure 2. Epicenter of the earthquake is shown as white star.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Static inversion to find the best fault plane location by moving fault plane (strike 171°) east-westward. Misfits are plotted versus longitude of epicenters, scaled by the minimum error. We choose a fault plane which is long enough to handle the error of event latitude.

DownLoad: Full-Size Img PowerPoint

To derive finite source kinematic models, we use the approach developed by Ji et al.(2002a, b) which allows the joint inversion of seismic waveforms and coseismic static displacements. As teleseismic and InSAR data provide complementary constraints on the spatiotemporal evolution of the rupture, we used both data sets to recover the slip history.

The determination of a source model for a given fault geometry, is an underdetermined problem due to numerous trade-off among model parameters including rise time and rupture velocity. Following the observation that most seismic ruptures seem to result from the propagation of slip pulse (Heaton, 1990), the inversion procedure can be regularized by assuming that the rupture consists of the propagation of a rupture front with a finite pulse width (Ji et al., 2002a, b; Wald et al., 1991). However, even in this case, trade-offs remain between the distributions of rupture velocity, rise time, and slip if only seismological data are inverted, especially when seismic station distribution is sparse and the rupture dimension is small. It is thus generally found that, in absence of geodetic constraints, a wide range of kinematic models, can fit the seismological observations equally well (Konca et al., 2007). The trade-offs are significantly reduced if geodetic observations of coseismic deformation are available and inverted jointly with the seismological observations. Even with these assumptions the determination of a source with finite dimensions remains generally underdetermined if the fault discretization is too fine. One way to regularize the inversion is setting some constraints on the roughness of the slip distribution (Ji et al., 2002a), which is the approach adopted here.

For each subfault, we solve for the slip amplitude and direction, rise time and rupture velocity. The rise time indicates the duration of time taken for the fault patch to slip to final static amount and is prescribed as a modified cosine function (Hartzell et al., 1996). The rupture velocity specifies the speed of the local rupture front. For each parameter, we specify extremal bounds and a discretization interval. The InSAR data is uniformly down sampled to about 1 000 points, and we solve for quadratic ramps in the InSAR data to correct for orbital errors not removed through baseline reestimation and interseismic deformation.

We define the best fit model as having the lowest objective function, given as E_wf +W_IE_I +W_SS +W_MM where E_wf is the waveform misfit (calculated for each wavelet channel with different weighting), E_I is the geodetic misfit, S is a normalized, second derivative of slip between adjacent patches (smoothing), M is a normalized seismic moment, and W_I, W_S and W_M are the relative weighting applied to the geodetic misfit, smoothing, and moment, respectively. The least squares misfits are calculated for the teleseismic and geodetic data. Here we test different values of W_I, and we found that equal weighting between the waveform and geodetic misfits did not significantly degrade the fits to the teleseismic or geodetic data between the individual and joint inversions given the normalizations employed by this method. The static Green's functions at free surface are calculated by following Xie and Yao (1989), using the same 1D crustal model as used in regional waveform modeling and teleseismic body-wave calculation.

We use a simulated annealing algorithm (Rothman, 1986) to find the best fitting model parameters for the joint inversions for coseismic slip. This nonlinear, iterative inversion algorithm is designed to avoid local minima by searching broadly through parameter space in initial steps, and then in later iterations to focus on regions that well fit the data (Sen and Stoffa, 1991).

The fault plane is divided into 15×15 subfaults, with subfault size of 2 km×2 km. During inversion, the slip amplitude is allowed to vary from 0 to 3 m, the rake angles varies for -120º to -40º which centered around the rake determined from local seismic waveforms, the average rupture velocity is selected to range from 2.0 km/s to 3.0 km/s at an interval of 0.1 km/s, and the rise time is allowed to range from 0.5 to 6.0 s at a 0.5 s interval. The rupture velocity is set below 3.0 km/s because no super shear rupture is expected for dip-slip earthquakes (Zhang and Chen, 2006). And as mentioned previously, we set equal weight for static and seismic data. As InSAR data is much less sensitive to deeper slip, it is hard to use InSAR data only to recover slip distribution especially when we have only one descending image. Complementally, our seismic data can tightly constrain the mechanism and moment. Thus, we fixed the total moment of the earthquake to the value (2.2×10¹⁸ N·m) we got from seismic waveform modeling during the joint inversion.

Our preferred finite slip model of the joint inversion is shown in Figure 6. The corresponding teleseismic waveform fits and InSAR fits are shown in Figures 7 and 8. The maximum slip amplitude we obtain is about 1.5 m. Most slip are distributed between 5 km and 11 km in depth, and about 10 km along strike. The average rupture velocity is about 2.2 km/s, which means that this is a relatively slow rupture earthquake. The centroid depth estimated from the finite slip model is about 9 km, which is a little shallower than the point source centroid depth. This might be because of the static data we use is more sensitive to shallower slip and prefers the slip closer to the free surface.

Figure  6.  Slip distribution of cumulative slip distribution (showing slip vectors, and amplitude of slip also represented by the color-coded) and isochrons of the seismic rupture determined by the joint inversion. The rupture times are given relative to the origin time, and the red star indicates the epicenter.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Comparison of the observed (black) and synthetic (red) teleseismic P-wave displacement seismograms. Station names are displayed to the left of the traces along with the azimuths (upper) and epicentral distances (lower) in degrees. Peak amplitude in micron of data is indicated above the end of each trace.

DownLoad: Full-Size Img PowerPoint

Figure  8.  Descending InSAR data recorded by the ENVISAT on 17 March and 8 September 2004 (left), line-of-sight (LOS) motion predicted from the preferred model (Figure 5) (middle) as well as the residual between the data and synthetic (right). The black arrow indicates the looking direction of the satellite.

DownLoad: Full-Size Img PowerPoint

4.   Discussion and conclusions

Using regional broadband seismic records, we first obtained the best point source mechanism of the 11th July 2004 Zhongba earthquake, which is 171º(342º)/42º(48º/-83º(-97º)/11 km for strike/dip/rake/depth respectively. With this densely distributed network, we can determine the strike of the mechanism with high accuracy. Our point source orientation parameters are consistent with mechanism from National Earthquake Information Center (NEIC), global CMT and Elliott et al. (2010), and the centroid depth of 11 km is consistent with NEIC (8 km) and Elliott et al. (2010) (10 km). The fault plane with strike of 171º and dip of 42º is identified as the real fault plane with the help of InSAR data. Then the distributed slip model is obtained by joint inversion of teleseismic P-waves and InSAR data.

Overall, the slip distribution pattern is similar to that by Li et al. (2005) and Elliott et al. (2010). Our joint inversion shows much more compact faulting area than the teleseismic inversion only by Li et al. (2005), which is not surprising because of the limitation in teleseismic waveform inversions only for M6 earthquakes. Slip distributed at 5–11 km depth is consistent with weak lower crust, with results by Zhao and Helmberger (1991) and Zheng (1995), who found that most events in this region occured shallower than 25 km. The normal faulting mechanism is probably due to east-west extension of the upper crust, which is either caused by increased potential energy of elevated crust or geodynamic processes that may be unrelated to plateau formation (Blisniuk et al., 2001). Thinning of the upper crust is most plausibly result of potential-energy increases resulting from spatially and temporally heterogeneous changes in thermal structure and density distribution within the crust and upper mantle beneath Tibet.