Evidence of Non-Linear Elasticity of the Crust from the Mw7.6 Manyi (Tibet) Earthquake

Satellite synthetic aperture radar (SAR) interferometry shows that the magnitude 7.6 Manyi earthquake of 8 November 1997 produced a 170 km-long surface break with up to 7 m of left-lateral slip, reactivating a N76E quaternary fault in western Tibet. The radar interferometric map reveals asymmetric, along-strike displacement profiles between the two sides of the surface rupture, a pattern that cannot be explained with linear elastic theory. This observation suggests that the elastic moduli of the crust in tension and in compression are different due to the presence of cracks in the crust at shallow depth. A model indicates that a ratio of 2 between compressive and tensile elastic moduli can account for the observed asymmetry, a ratio consistent with laboratory and borehole measurements.

Geodetic measurements of the static surface displacement field produced by large earthquakes provide information about their source mechanism. Linear models of dislocation in an elastic half space (1) are widely used to represent earthquake faults and generally provide a satisfactory fit to near and far field displacement data obtained with conventional geodetic techniques (2). However, laboratory experiments and in situ measurements in boreholes have shown that many crustal rocks exhibit a non-linear elastic behavior in compression and tension with a dependence of the Young's modulus on the minimum principal stress (3). The effects of non-linear elastic properties of crustal rocks have never been observed in geodetic signals associated with elastic crustal strain because of the low sampling density of techniques that make use of ground-based instrument arrays (4). Here we have used SAR interferometry (5) to analyze coseismic surface displacement in the near field of the 1997 Manyi, Tibet, earthquake rupture.

The magnitude (Mw ) 7.6 Manyi earthquake occurred on 8 November, 1997 in the western part of Tibet with a hypocenter depth of 15 km. The region is flat with an elevation of 5000 m near the western end of the Kunlun fault (6). The focal mechanism (7) indicates left-lateral, strike-slip motion on a fault plane parallel to a mapped quaternary fault (8) (Fig. 1). The aftershocks define an east-northeast trend approximately aligned with the fault (9). To map the earthquake displacement, we combined SAR images acquired by the European Remote Sensing satellite (ERS-2) on three adjacent tracks before and after the event (10) (Fig. 1). The topographic signal was removed from interferograms (11) using a digital elevation model (12). The resulting interferogram shows the surface displacement produced by the earthquake in the direction parallel to the radar line of sight (Fig. 2A, 13).

The three-track mosaic (Fig. 2A) covers the 170 km -long section of the fault that ruptured during the earthquake. The two half-lobe patterns north and south of the rupture are consistent with a left-lateral motion on the fault. The maximum range displacement is observed immediately east of the center of the rupture and exceeds 1.2 m on each side of the fault (Fig. 2B). Assuming left-lateral slip on a N76E fault, this indicates that the total relative offset reached 7 m. The smoothness of the contours of the surface displacement field suggests a smooth fault geometry and a smooth distribution of slip on the fault at shallow depth. The displacement profiles perpendicular to the fault indicate a relatively shallow fault depth. Assuming a uniform slip distribution with depth, elastic dislocation modeling (1) shows that a fault depth of ~8 km is required to fit the observed displacement profiles (Fig. 2B).

Surface condition changes probably related to strain, ground shaking, and to the presence of water in ponds are responsible for the loss of phase coherence near the fault (Fig. 2A, 3A). To map the surface rupture associated with the earthquake we used the two-component offset field obtained by sub-pixel registrations of the radar amplitude images acquired before and after the earthquake (14). Unlike the interferometric phase, which is prone to coherence loss due to changes in surface conditions, the offset field provides surface displacement estimates in azimuth and range in the zone bordering the surface break. The dislocation produced by the earthquake was readily visible in the offset fields so that the 170 km-long surface break could be mapped with a precision of a few image pixels (Fig. 3B).

To produce the along strike slip distribution curve we first constructed two profiles of displacement with respect to the fixed, far-field reference south and north of the fault. For each point along the rupture, displacements observed in the nearest, coherent areas south and north of the fault were extrapolated to the fault (Fig. 3A). The line of sight projection of the slip distribution curve is the difference of the north and south displacement curves (Fig. 4A). In the regions of overlap between adjacent tracks, the observed range displacement is systematically larger in the eastern track than in the western track (Fig. 4A). This discrepancy is due to the difference in the line of sight incidence angle between adjacent orbits from where the ground is imaged (15). If we assume that the displacement of the ground is horizontal, it is possible to correct for the effect of line of sight projection by dividing the observed range change at each point by the sine of the local incidence angle (16). The result is the projection of the ground displacement on the horizontal axis perpendicular to the satellite track (cross-track) (Fig. 4B). The agreement between the cross-track displacement curves observed from adjacent orbits confirms that the displacement vector is horizontal along the fault trace (Fig. 4B).

The observed slip distribution is bimodal. The main event ruptured the eastern 150 km-long section of the fault with a maximum strike-slip offset of 7 m. A sub-event ruptured the western 40 km-long section of the fault with a maximum strike-slip offset of 2.6 m. If we assume a uniform slip distribution over a depth of 8 km and an elastic shear modulus of 3.3 x10¹⁰ N/m², the main event implies a geodetic moment of 1.4 x10²⁰ Nm, consistent with the seismic data (17).

The north and south along-strike displacement profiles show a marked asymmetry with respect to the fault. The profiles have a triangular shape between 50 and 150 km but the maximum displacement of the southern side occurs 100 km from the western end of the fault, 20 km west of where it occurs along the northern side (Fig. 4, A and B). The north and south displacement curves are the symmetric images of one another with respect to the point located at 110 km, 0 cm on the plots of Figures 4 A and B. Such a pattern cannot be accounted for by models based on linear elastic theory (18).

A possible explanation for this asymmetry is the presence of cracks in the crust. The effects of cracks in an elastic medium is to reduce its rigidity (19). This effect is non-linear and depends on the sign of the effective normal stress resolved on the crack plane. Cracks would open in a tensile stress regime to a greater extent than they would close under a compressive regime of the same magnitude. Furthermore, new cracks can be created in tension but not in compression. Thus, the rigidity of the crust in compression should be larger than its rigidity in tension. Such a non-linear elastic property has been observed on many crustal rock samples in laboratory tests and in borehole measurements (3).

To test the effect of non-linearity on the static displacement field produced by a strike-slip dislocation in an elastic medium, we used a three-dimensional finite element code (20) and assumed different values of the Young's modulus in compression and tension as a first approximation of the dependence of the Young's modulus on the minor principal stress. We modeled the earthquake dislocation as a 100 km long and 15 km wide, vertical fault plane. The imposed strike-slip distribution has a trapezoidal shape with a constant maximum value of 7 m along the 20 km-long central section, linearly decreasing to 0 m toward the two ends of the fault. As a double couple mechanism, the planar fault geometry and the symmetry of slip distribution define four quadrants in the elastic half-space, separated by the fault plane and the auxiliary plane. For an east-west, left-lateral, strike-slip fault, the northeast and southwest quadrants are dilational and the northwest and southeast quadrants are compressional. We assigned different values of the Young's modulus to the compressional and dilational quadrants and compared 3 cases (21). The first case is the linear case where compressional (Ec) and dilational (Ed) Young's moduli are equal, as in the analytic solution (1). The two other cases are non-linear and assume compressional Young's modulus values twice and four times larger than the tensional Young's modulus value. The linear case shows a slight asymmetry between the two curves with the maximum range displacement in the south occurring ~10 km east of where it is occurring on the north side (Fig. 4C). This small asymmetry is in the opposite direction from the observed asymmetry and is due to the effect of projection of the 3-component displacement vector on the radar line of sight (22). The non-linear cases (Ec/Ed=2 and 4) lead to a larger asymmetry of the fault-parallel profiles, similar to the asymmetry of the observed profiles (Fig. 4, A and C). North and south of the fault, the location of the maximum displacement is shifted toward the tensile, weaker quadrant.

The north-south asymmetry of the displacement field can also be seen in the interferometric map (Fig. 2 A and B). The core of the contour lines (fringes) in the south appears to be centered 20 km west of where it is centered on the north side of the fault. This asymmetry tends to disappear as the distance to the fault increases suggesting that non-linear elasticity is mostly effective in the shallow crust (~3-5 km).

Crustal fluids at lithostatic pressure are required to maintain Mode-1 deformation of cracks effective at depth. Hydrological studies have shown that earthquake strain can mobilize crustal fluids to depths of 5 km or more (23). However, processes involving crustal fluid flow are known to have relaxation times of the order of a few days to a year, depending on cracks connectivity and the porosity and permeability of the rocks (23, 24). The SAR data acquired on ERS tracks 76, 305, and 33 used in this study cover 8, 24, and 40 days of the post-seismic period, respectively (10). The agreement between sections of the displacement curves observed from adjacent orbits in overlapping parts of the satellite swaths (Fig. 4B) indicates that the observed displacement was not significantly affected by any relaxation process during those periods of observation. Thus, the observed asymmetry of the displacement field near the fault has to be explained by a short term process associated with the dislocation. These findings suggest that, at shallow depth (<3km), efficient interconnection between cracks and the fact that some cracks may not be completely filled with fluids lead to a short-term, non-linear elastic response of the crust to coseismic stress changes (25). However, the short-term response does not preclude longer-term relaxation of crustal fluid pressure gradients contained within less connected cracks, especially at greater depth (>3 km). The geodetic signature of such processes is not observed in the along strike displacement profiles because of the short period of observation after the earthquake and the fact that processes occurring at greater depth have a small influence on the surface displacement field at short distance from the fault.

We conclude that non-linear elasticity related to the presence of cracks in the shallow crust provides a plausible explanation for the observed asymmetric displacement pattern across the Manyi earthquake fault. Simple modeling suggests that a ratio of 2 to 4 between compressional and extensional Young's moduli in the shallow part of the seismogenic crust can account for the radar observations. Although extrapolating values of elastic parameters derived from laboratory measurements to large volumes of crust must be done with caution, a contrast of 2-4 between compressional and tensional rigidity constants is commonly observed in laboratory test samples (3).

REFERENCES AND NOTES

1. Y. Okada, Bull. Seismol. Soc. Am. 75, 1135 (1985).

2. See for example, Hudnut et al., Bull. Seismol. Soc. of Am., 84, 625 (1994), J. Freymuller, N.E. King, and P. Segall, ibid., p. 646.

3. Many crustal rocks exhibit non-linear elastic properties with a Young's modulus dependent on the confining pressure. Tests run on laboratory samples showed that rocks such as sandstone, marble, limestone, quartzite, and granite have a Young's modulus in tension smaller than the Young's modulus in compression in a ratio ranging from 1/2 to 1/10 [J.C. Jaeger and N.G. Cook, Fundamentals of Rocks Mechanics (Chapman and Hall, London, ed. 3, 1979), chap.10 and 12; B.C. Haimson and T.M. Tharp, Soc. Petrol. Eng. J., 14, 145 (1974); F.K. Khan and F.K. Irani, Mech. of Mater., 6, 271 (1988); B. Stimpson and R. Chen, Can. Geotech. J., 30, 338 (1993); E.T. Brown, J.W. Bray, and F.J. Santarelli, Rock Mech. and Rock Eng., 3, 189 (1989)].

4. Current geodetic techniques are based on arrays of a few tens of benchmarks with station spacing ranging from tens of meters for small aperture trilateration arrays [G.A. Sylvester, Geophys. Res. Lett., 20, 1079 (1993)] to ~10 km, or more, for Global Positioning System (GPS) networks [e.g., Z.K. Shen et al., Bull. Seismol. Soc. of Am. 84, 780 (1994); F.K. Wyatt, D.C. Agnew, and M. Gladwin, ibid., p. 768].

5. The advantage of SAR interferometry over ground based geodetic techniques is to provide spatially continuous maps of surface displacement over broad areas [H. Gabriel, R. Goldstein, and H. Zebker, J. Geophys. Res., 94, 9183 (1989); Massonnet et al., Nature, 364, 138 (1993); H. Zebker et al., J. Geophys. Res., 99, 19617 (1994)]. The interferometric maps we generated here cover three swaths of 100 km in width and have a pixel size of ~90 m on the ground, after data averaging and geocoding.

6. Landsat images reveal quaternary faulting in the area of Tibet extending west of the western end of the Kunlun fault and south of the Altyn Tagh fault including east-northeast strike-slip faults and north-south normal faults [(8); R. Armijo et al., J. Geophys. Res., (1989)].

7. The Harvard centroid moment tensor solution is: strike=262, dip=79, rake=4, depth=15km, and moment=2.8 x 10²⁷ Dcm; the National Earthquake Information Center (NEIC) focal mechanism is: strike=71, dip=90, rake=4, depth=29 km, and moment=1.5 x 10²⁰ Nm. Both mechanisms indicate mostly purely left-lateral, strike-slip motion on an east-northeast striking fault.

8. Tapponnier and P. Molnar, J. Geophys. Res., 82, 2905 (1977).

9. The seismicity catalog of the National Earthquake Information Center (NEIC) includes 80 aftershocks in 1997 and 1998, 16 of which have magnitude grater than 4.5. The largest event listed is of M 5.3 and occurred on 9 November, 1997, ~150 km south of the surface rupture. We identified in the radar data the surface displacement of a shallow earthquake that could not be associated with any of the aftershock locations listed in the NEIC catalog (26) (Fig. 1).

10. ERS-2 SAR data were acquired on 16 March and 16 November, 1997 for track 76, 19 August and 2 December, 1997 for track 305, and 22 May and 18 December, 1997 for track 33 (Fig. 1). Each pair of raw SAR images was processed to an interferogram using the JPL radar processing software ROI_PAC.

11. In each SAR image pixel, the phase value represents a measure modulo-2p of the distance between the radar antenna and the ground. An interferogram is obtained by averaging the product of one complex SAR image and the complex conjugate of another image, after sub-pixel coregistration. The phase of each pixel in an interferogram is the difference of the phase of the corresponding pixels in the two SAR images. The interferometric phase depicts variations of the antenna-ground path length difference between the two images.

12. The topographic signal can be removed from an interferogram by simulating the topographic phase using a digital elevation model and the geometric parameters defining the configuration of the orbits at times of acquisitions of the data [Massonnet et al., Nature, 364, 138 (1993); Murakami et al., J. Geophys. Res., 101, 8605 (1996)]. For this study, we used the 90 m resolution Digital Terrain Elevation Data (DTED) [National Imagery and Mapping Agency].

13. For ERS, the satellite line of sight is nearly perpendicular to the orbit. On descending passes, the observation is made from the southeast with an incidence angle of 18.5 in near range (east edge of image) and 27 in far range (west edge of image) [European Space Agency, ERS-1 System (ESA Publication Division, ESTEC, Noordwijk, Netherlands, 1992)].

14. Interferometric processing of SAR images requires sub-pixel registration of the two images. This is achieved by searching for local peaks of amplitude correlation between the two images. This procedure provides us with azimuth (parallel to satellite track) and range (perpendicular to satellite track) offset fields that reflect the geometric shift and distortions between the two images. Typical noise in offset fields in the present data is approximately 1/20th of the pixel size. The size of an ERS SAR single look pixel on the ground is 4 m in azimuth and 20 m in range. This technique is similar to the approach developed by R.E. Crippen [Episodes, 15, 56 (1992)] using visible satellite images (SPOT).

15. At the latitude of 35N on descending orbit passes, the incidence angle difference between lines of sight of adjacent orbits varies from 4.6 to 4.9 between the far range (western edge) and the near range (eastern edge) of the area of swath overlap (13).

16. The line of sight projection of a displacement vector (e,n,u) in the local east, north, up reference frame is r = [e cos(a) - n sin(a)]sin(i) + u cos(i), where a is the angle between the azimuthal direction of the displacement vector and the direction perpendicular to the satellite track, and i is the incidence angle (13). For horizontal displacement fields (u=0), the observed cross-track component of the displacement vector is r/sin(i) and is independent of the incidence angle.

17. The source time function of the earthquake depicts a similar bimodal shape with a main pulse moment of ~1.5 x10²⁰ Nm in the first 30 seconds followed by a smaller pulse of 10 seconds [University of Michigan web site: www.geo.lsa.umich.edu/SeismoObs/STF.html (1997); A. Velasco, personal communication (1998)].

18. Seismological and radar data indicate that only one fault was involved in the earthquake. Shear stress across the fault must be continuous (equilibrium conditions). Observed displacement asymmetry indicates that strain is not continuous across the fault. Thus, elastic properties must be different between the two sides of the fault. The particular form of asymmetry observed in the displacement curves (Fig. 4B) strongly suggests that elastic properties are correlated with volume strain and not related to crustal heterogeneities along the rupture.

19. W.R. Delameter, G. Herrmann, and D.M. Barnett, J. Appl. Mech., 47, 74 (1975); E.R. Ivins and G.A. Lyzenga, Phil. Trans. R. Soc. Lond., A 318, 285 (1986).

20. F. Saucier, Ph.D. thesis, University of Oregon (1991).

21. This approach is akin to a first-order perturbation approach with respect to a linear model since we implicitly assume that the locations of compressional and tensile volumes are unchanged by non-linearity.

22. Given the symmetry of the fault geometry and slip distribution of the model, the along-strike distribution of the fault-parallel, horizontal component of the displacement vector is symmetric with respect to the auxiliary plane and the distribution of the vertical component is anti-symmetric. The line of sight projection of the vector combines all components of the displacement vector (16) and thus depicts the anti-symmetry of its vertical component with respect to the auxiliary plane.

23. R. Muir-Wood and G.C.P. King, J. Geophys. Res., 98, 22,035 (1993).

24. A. Nur and J.R. Booker, Science, 175, 885 (1972); D.L. Anderson and J.H. Whitecomb, J. Geophys. Res., 80, 1497 (1975); K.W. Hudnut, L. Seeber, J. Pacheo, Geophys. res. Lett., 16, 199 (1989); C.H. Scholz, Geology, 2, 551 (1974); J.B. Rundle and W. Thatcher, Bull. Seismol. Soc. Am., 70, 1869 (1980); G. Peltzer, P. Rosen, F. Rogez, and K. Hudnut, Science, 273, 1202 (1996).

25. A drop of seismic velocities at depths shallower than 3-5 km is observed in many places [P. Reasenberg and W.L. Ellsworth J. Geophys. Res., 87, 10,637 (1982); A.W. Walter and W.D. Mooney, Bull. Seismol. Soc. Am., 72, 1,567 (1982); L.D. Dietz and W.L. Ellsworth, Geophys. Res. Lett., 17, 1,417 (1990); C. Dorbath, D. Oppenheimer, F. Amelung, and G. King, J. Geophys. Res., 101, 27,917 (1996)] and is generally attributed to the presence of open cracks, not saturated with fluids in the shallow crust.

26. At longitude E87.17 and latitude N34.60 the interferogram shows a double lobe, ~10 km-long pattern, with a maximum range change corresponding to ~10 cm of surface uplift and bounded to the north by an east-west, ~8 km-long line of phase discontinuity. We interpret this feature as the surface displacement field associated with a shallow earthquake having occurred between 19 August, 1997 and 2 December, 1997. The observed displacement and phase discontinuity suggest that the event has left-lateral and thrust components and a fault plane steeply dipping to the south (Crampé and Peltzer, manuscript in preparation). We note however, that the geographic location of this feature is 30 km away from the nearest aftershock epicenter determined by the NEIC during the same period. This may be due to either large errors in the localization of aftershocks or the fact that this shallow event occurred during the propagation of the main rupture making its seismic record indistinguishable from that of the main event.

27. We thank P. Bernard, F. Cornet, E. Ivins, and H. Kanamori for fruitful discussions, P. Lundgren for help with the modeling, and two anonymous reviewers for constructive remarks. ERS data were provided by the European Space Agency. The work was done at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

Fig. 1 Active fault map of the Manyi earthquake area over shaded 30 arc seconds topographic map (GTOPO30). Faults in black are from (8). Red line is surface rupture mapped using offset fields (14). Focal mechanisms and epicenter locations determined by Harvard and by the National Earthquake Information Center (NEIC) are depicted (7). Yellow symbols show seismicity between 8 November, 1997 and 31 December, 1998 from NEIC catalog. Circles and stars depict aftershocks with magnitude smaller and greater than 4.6, respectively. Black dot at longitude E87.17 and latitude N34.60 is a shallow aftershock visible in SAR interferogram and nearby red line is associated surface rupture (26). Boxes (solid line) represent area covered by ERS tracks used here: track 76 (left), track 305 (center), and track 33 (right). Box in dashes shows area covered by Fig. 2A. Dashed lines indicate location of profiles shown in Fig. 2B.

Fig. 2 (A) Interferometric map showing the 8 November, 1997 coseismic surface displacement field. One full color cycle (blue-red-yellow-blue) represents 50 cm of ground shift away from satellite along radar line of sight (13). Uncolored areas are zones of low phase coherence that have been masked before phase unwrapping. The map is a mosaic of three interferograms covering ERS tracks 76, 305, and 33 (Fig. 1). Small color discontinuities observed at frame boundaries are due to differences of incidence angles between adjacent tracks in overlapping regions (15). (B) Observed (red and blue) and modeled (gray) line of sight displacement profiles perpendicular to fault. Red profile (left profile on Fig. 1) intersect the fault where maximum displacement is observed south of the fault. Blue profile (right profile on Fig. 1) is where maximum displacement is observed north of the fault. Fault is modeled in an elastic half-space (1) as a vertical plane, 8-km deep, with 7 m of slip uniformly distributed with depth.

Fig. 3 (A) Detail of interferogram shown in Figure 2A centered on surface rupture at longitude 87.3. Black dots indicate surface rupture as mapped using offset fields (14). Range displacement is measured north and south of the fault at white dots and values are extrapolated up to the rupture using displacement gradient estimated along line segments (black lines). (B) Detail of the offset field in azimuth (14) covering same area as in (A). The surface break follows the discontinuity in the offset field. White pixels over dry lake area north of fault are erroneous offset estimates attributable to lack of features in amplitude image.

Fig. 4 (A) Observed range displacement distribution profiles along south (blue) and north (green) sides of fault. Fault slip distribution (relative displacement between north and south sides) is shown in red. (B) Same as (A) for cross-track component of surface displacement. Profiles are obtained by normalizing profiles shown in Figure 4A by the sine of the incidence angle (13). For horizontal surface displacement the normalized quantity is independent of the incidence angle. 0-8 m vertical scale indicates corresponding strike-slip displacement on fault. (C) Modeled, along-strike range displacement profiles assuming linear and non-linear elastic crust. Colors are as in (A). Dashed lines correspond to linear (Ec/Ed=1) case. Dotted lines (Ec/Ed=2) and solid lines (Ec/Ed=4) correspond to non-linear cases (see text). A constant incidence angle of 23 (13) is assumed for 3-component displacement vector projection. Profiles are measured at a distance of 5 km from the fault.