何玉梅---中国科学院地质与地球物理研究所.pdf

Bulletin of the Seismological Society of America, Vol. 93, No. 5, pp. 2253–2263, October 2003 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space by Yu-Mei He, Wei-Min Wang, and Zhen-Xing Yao Abstract Based on the generalized reflection and transmission coefficient matrix method, formulations for surface static displacements in a layered half-space are extended to include tensile and inflation point sources from a point pure shear dislocation source. Equations for calculating internal displacement fields from these sources are also derived. The validity of the formula and precision of the new method are illustrated by comparing the consistency of our results and the analytical solutions given by Okada’s (1985, 1992) code in a homogenous half-space and Wang et al.’s (2003) numerical solutions in a multilayered half-space. We also study the effect of a layered half-space on the surface displacement created by various finite faults. Several typical velocity structures in reality are selected. For strike-slip, reverse dipslip, and tensile finite-fault models, the focal depth is very sensitive to the presence of the layered model. The slip displacement is more sensitive to the layered model in the case of the normal dip-slip sources. More numerical tests show that the sensitive slip is mainly due to the ultralow-velocity topsoil. For inflations, the source depth and volume change also altered due to the layered model. Introduction Static deformation data can be used to investigate geological fault movements and related stress distributions. Maruyama (1964) gave the analytical solution of a static elastic dislocation in an infinite or semi-infinite media. The analytical solutions for surface and internal deformations generated by shear and tensile faults in a half-space were derived by Okada (1985, 1992) and were widely used. However, a homogeneous half-space model may oversimplify the real Earth. For example, the existence of a soft top layer tends to generate surface displacements that are larger near the epicentral area but decrease rapidly with distance (Sato and Matsu’ura, 1973). Numerical tests also revealed that, for a pure shear thrust fault with a 30⬚ dip angle, the response generated using a homogeneous half-space model can underestimate the displacement by up to 10% if the observation point is located in the hanging wall, and it will overestimate the displacement by up to 30% if the observation point is located in the footwall (Ji et al., 2001). Therefore, when dealing with geodetic data, a more realistic layered Earth model is often preferred (Hearn et al., 2002). Using the Thomson–Haskell matrix method, BenMenahem and Singh (1968), Singh (1970), Sato (1971), and Sato and Matsu’ura (1973) presented algorithms for calculating static displacement for a point shear dislocation in a layered half-space model. After improving the method, Wang et al. (2003) gave an efficient algorithm for three different types of point double-couple dislocations. Roth (1990) and Ma and Kusznir (1994, 1995) extended this method to deal with internal deformation and strain. Kennett (1974) partitioned the Thomson–Haskell propagator matrix into submatrices and introduced them as the generalized reflection-transmission coefficient matrices. Based on the generalized reflection-transmission coefficient matrix and the discrete wavenumber summation method (Kennett, 1983; Luco and Apsel, 1983; Yao and Harkrider, 1983), Xie and Yao (1989) developed an approach for dealing with the static problem. They derived expressions of static deformation for shear dislocation sources in a layered half-space. In this article, following Xie and Yao (1989), we expand their approach to include tensile and inflation sources. Formulations for both surface and internal displacement fields are derived. Numerical tests are conducted to examine these equations. To validate our method, we compare our point source solutions to Okada’s (1985) analytic results and Wang et al.’s (2003) numerical solutions. Using the new method, the responses for various finite faults buried in different layered models are calculated and compared to Okada’s rectangular fault solutions in a homogenous halfspace. The effects of focal depth, slip, dip, and rake angles are investigated. Theory In the frequency domain, the homogeneous elastostatic equation (without the source term) can be derived from the 2253 2254 Y.-M. He, W.-M. Wang, and Z.-X. Yao elastodynamic equation by making the circular frequency x → 0 (Zhu and Rivera, 2002). Since the Lamé theorem remains true even for static problems, displacement potentials can be adopted and the vector elastostatic equation can be turned into simpler scalar Laplacian equations. For an isotropic layered half-space, the solutions of displacement potentials are similar to those in the elastic dynamic problem. It is convenient to derive source functions following methods used in wave motion problems (Yao and Harkrider, 1983). Although there is no wave propagation in a static problem, the existence of an interface still distorts static displacement fields on both sides of it. We will borrow ideas of reflection and transmission from the wave problem but rename them as static reflection and transmission. Similarly, without frequency dependency, terms related to exp(ⳮkz) and exp(Ⳮkz) are no longer linked with up- and downgoing waves. They are actually different eigenfunctions representing static displacement fields. However, for the ease of description, we sometimes still call them up- and downgoing disturbances of static displacements. Following these processes, the original generalized reflection and transmission coefficient matrix method for the wave propagation problem can be modified for dealing with the static problem. x y U3 TENSILE z Figure 1. where S is the area of the fault plane, h is the z coordinate of the point source, Jm(kr) is the ordinary Bessel function of order m, and Tensile and Inflation Source Functions Xie and Yao (1989) derived source terms for static shear dislocations. Similarly, we can obtain source functions and related source coefficients for an arbitrarily oriented tensile dislocation. According to Steketee (1958), the displacement field ui(x1, x2, x3,) due to a dislocation Duj(n1, n2, n3) across a surface S in an isotropic medium is given by ui ⳱ 冮冮 S  Duj kdjk ⳵uni ⳵uij ⳵uki Ⳮl Ⳮ n dS, ⳵nn ⳵nk ⳵nj k 冢 冣 2 S u⳱ U A 4p m⳱0 3 m 兺 2 ⬁ 冮 Pe m 0 ⳮk|zⳮh| ⬁ Jm (kr)kdk, w⳱ S U3 Am 4p m⳱0 兺 冮 SV e 2 v⳱ S U A 4p m⳱1 3 mⳭ3 ⬁ 兺 m 0 ⳮk|zⳮh| Jm (kr)kdk, 冮 SH e 0 m ⳮk|zⳮh| Jm (kr)kdk, A0P0 ⳱ A0(1)P0(1) Ⳮ A0(2)P0(2), A0SV0 ⳱ A0(1)SV0(1) Ⳮ A0(2)SV0(2), A0(1) ⳱ 1, 2 A(2) 0 ⳱ ⳮsin d, A1 ⳱ sin(2d) sin h, A2 ⳱ ⳮsin2d cos(2h), A4 ⳱ ⳮsin(2d) cosh, A5 ⳱ sin2d cos(2h) (1) where djk is the Kronecker delta, k and l are Lamé constants, nk is the direction cosine of the normal to the surface element dS, and uij is the ith component of the displacement at (x1, x2, x3,) due to a point force at (n1, n2, n3) acting in the jth direction. For a fault with dip d and strike h and buried in an infinite space (shown in Fig. 1), a tensile dislocation can be expressed as U3m. Here m is the direction of the rupture vector of the fault plane, and U3 is the scalar displacement perpendicular to the fault. Then, the displacement potential functions ␾, w, and v of a static tensile point source can be written as (2) Geometry of the dislocation model. (3) are orientation factors. Pm , SVm , and SHm are source coefficients and can be expressed as P(1) 0 ⳱ ⳮD , 1 Ⳮ D P(2) 0 ⳱ 1 ⳮ 4D ⳮ3 , SV0(2) ⳱ , 2(1 Ⳮ D) 2(1 Ⳮ D) P1 ⳱ ⳮeD , 1 Ⳮ D SV1 ⳱ e , 1Ⳮ D P2 ⳱ 1 , 2(1 Ⳮ D) SV2 ⳱ ⳮ1 , SH2 ⳱ ⳮ1, 2(1 Ⳮ D) SV(1) 0 ⳱ ⳮ1 , 1Ⳮ D (4) SH1 ⳱ e, ⳮ ⳮ ⳮ where Pm, SVm, and SHm represent PⳭ m , Pm , SVm , SVm ⳮ and SHⳭ m , SHm , plus and minus superscripts denote up- and downgoing static disturbances, e is ⳮ1 for the minus superscript and equal to 1 for the plus superscript, and D ⳱ (k Ⳮ l)/(k Ⳮ 3l). Compared with static point shear dislocation source functions (Xie and Yao, 1989), we have the same source coefficients except for P0 and SV0. For a tensile dislocation, 2255 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space o its rupture vector parallels the direction of the fault-plane normal, and its rake k disappears, which means orientation factors for shear and tensile sources are different. For the inflation point source, we have a simple expression that u⳱ V 4p ⬁ 冮 0 P0eⳮk|zⳮh|J0 (kr)kdk, r z1 1 z3 zs (5) j where V is the volume change, and P0 ⳱ _ Dⳮ 1 . 1 Ⳮ D zn z +n (6) Static Displacement Solutions for Shear, Tensile, and Inflation Sources in a Layered Half-Space Consider n ⳮ 1 parallel, homogeneous, and isotropic layers overlaying a semi-infinite medium (Fig. 2). A righthanded cylindrical coordinate system is used, and the z axis is vertically downward. The layers and interfaces are labeled according to their distances away from the free surface. Following the generalized reflection and transmission coefficient matrix method (shortened to generalized R/T coefficient matrix method [Kennett, 1983]), using the static shear dislocation source given by Xie and Yao (1989) and tensile source derived in the last section, the surface and internal displacement fields can be expressed as 3 wz ⳱ 2 S Uj Ajm 4p j⳱1 m⳱0 兺 兺 3 2 S qr ⳱ Uj Ajm 4p j⳱1 m⳱0 兺 兺 3 mh ⳱ Figure 2. The geometry of the layered model. The asterisk denotes the source that is located at depth zs and on the z axis. The source layer is divided into two ⳮ layers. zⳭ n and zn denote points just below and above interface zn. Similar to these expressions, the surface and internal displacement fields due to a point inflation source can be written as wz ⳱ 冮 w J (kr)kdk, m m ⬁ 冮 0 2 S Uj AjmⳭ3 4p j⳱1 m⳱0 兺 兺 m [qm J˙ m(kr) ⳮ mm Jm (kr)]kdk, kr ⬁ 0 m m where U1, U2, and U3 correspond to the strike-slip, dip-slip, and tensile components. The wm, qm, and mm, m ⳱ 1, 2, 3 represent the displacement fields under vector surface harmonic base, and their detailed expressions will be given later. The orientation factors are A10P0 ⳱ A10P(2) 0 , A20P0 ⳱ A20P(2) 0 , (1) (1) (2) A30P0 ⳱ A30 P0 Ⳮ A(2) 30 P0 , (1) (2) (2) A30SV0 ⳱ A(1) 30 SV0 Ⳮ A30 SV0 , A20 ⳱ sind cosd, A21 ⳱ ⳮcos(2d) sinh, A22 ⳱ sind cosd cos(2h), A24 ⳱ ⳮcos(2d) cosh, A25 ⳱ ⳮsind cosd sin(2h), A10 ⳱ 0, A12 ⳱ sind sin(2h), A14 ⳱ ⳮcosd sinh, A15 ⳱ sind cos(2h), (2) A30 ⳱ ⳮsin2d, A31 ⳱ sin(2d) sinh, V 4p ⬁ 冮 w J (kr)kdk, 0 0 0 (9) ⬁ 冮 q J˙ (kr)kdk. 0 0 0 冢wq 冣 ⳱ R (I ⳮ R R) m ev 1S D ⳮ1 SL FS ⳮ1 T1S U (I ⳮ RD RU ) m Ⳮ Pm RSL D Ⳮ SVm Pⳮ Ⳮ mⳮ , SVm 冤 冢 冣 冢 冣冥 1S ⳮ1 1S FS ⳮ1 mm ⳱ 2(I ⳮ RD,L ) TU,L(I ⳮ RSL D,LRU,L) Ⳮ ⳮ (RSL D,LSHm Ⳮ SHm ). A11 ⳱ cosd cosh, A(1) 30 ⳱ 1, V 4p For calculating the surface displacement, wm, qm, and mm can be expressed as m m qr ⳱ (7) 冮 [q kr J (kr) ⳮ m J˙ (kr)]kdk, m zn n z ⬁ 0 z2 2 (8) (10) (11) Then the surface displacements are obtained by putting equations (10) and (11) into equations (7) and (9). Here I is a unit matrix, the subscript L refers to the SH problem. Source coefficients Pm, SVm, and SHm are given by equation (4). For the P-SV problem, (i.e., equation 10), A32 ⳱ ⳮsin2d cos(2h), A34 ⳱ ⳮsin(2d) cosh, A35 ⳱ sin2d cos(2h). R⳱ ⳮ 0 D 冢1/D 0冣 1 1 (12) 2256 Y.-M. He, W.-M. Wang, and Z.-X. Yao is the static reflection coefficient matrix at the free surface, and ⳮ1 冢1/D 1/D 1 冣 Rev ⳱ (1 Ⳮ D1) 1 (13) 冢 冣 qm RS FR ⳮ1 FS ⳮ1 ⳱ Rev(I ⳮ RD RU ) TURS(I ⳮ RSL D RU ) wm 冤R 冢SVP 冣 Ⳮ 冢SVP 冣冥 , (14) Ⳮ m Ⳮ m SL D ⳮ m ⳮ m FR ⳮ1 RS FS ⳮ1 mm ⳱ Rev,L(I ⳮ RRS TU,L (I ⳮ RSL D,L RU,L) D,L RU,L) (15) Ⳮ ⳮ RSL D,L SHm Ⳮ SHm , ( ) where 1 ⳮ1 R 冢11 ⳮ11 冣 Ⳮ 冢ⳮ1 ⳮ1冣 (16) Rev,L ⳱ 1 Ⳮ RFR U,L (17) FR U and are static receiver functions at zR. RFR U is the upgoing generalized reflection coefficient matrix between the free surRS RS face z1 and zⳮ R . RD and T U are the down- and upgoing static generalized R/T coefficient matrixes between zⳭ R and zⳮ . Apparently, z ⳱ 0 makes equations (14) and (15) deS R generate to the surface displacement equations (10) and (11). For receivers below the focal depth (i.e., zR ⬎ zS), we have 冢wq 冣 ⳱ R (I ⳮ R R ) m ev m SR RL ⳮ1 SL ⳮ1 TDSR I ⳮ RFS U D U RD ( 冤 Ⳮ ⳮ (ⳮRFS U,L SHm ⳮ SHm ), (19) where 1 is the receiver function matrix. Their detailed expressions are given in the Appendix. RFS U is the upgoing static generalized reflection coefficient matrix between the free surface 1S 1S z1 and source depth zⳮ s . RD and TU are down- and upgoing static generalized R/T coefficient matrixes between zⳭ 1 and SL zⳮ , and R is the downgoing static generalized reflection s D Ⳮ coefficient matrix between zⳭ s and zL (see Fig. 2). Detailed expressions of these static R/T coefficient matrixes at interfaces are also given in the appendix. For the SH problem, the static receiver function and reflection coefficient at the free surface are 2 and 1, respectively. Similar to the wave motion problem (Kennett, 1983), the surface static solution (equations 10 and 11) can be extended to include internal displacement fields. For receivers above the focal depth, that is, zR ⬍ zS, we have Rev ⳱ RL ⳮ1 SR SL ⳮ1 mm ⳱ Rev,L(I ⳮ RSR TD,L(I ⳮ RFS U,LRD,L) U,LRD,L) ) 冢 冣 ⳮ 冢 冣冥, (18) Ⳮ FS Pm ⳮRU SVⳭ m Pⳮ m SVⳮ m Rev ⳱ 1 ⳮ1 冢11 ⳮ11 冣R Ⳮ 冢ⳮ1 ⳮ1冣 RL D (20) Rev,L ⳱ RRL D,L Ⳮ 1 (21) and are static receiver functions for P-SV and SH problems at zR, SR respectively. RSR D and TU are down- and upgoing static genⳮ eralized R/T coefficient matrixes between zⳭ S and zR . Substituting equations (14) and (15) or (18) and (19) into equations (7) and (9), we can calculate internal displacement fields in a layered model. Numerical Examples After the derivation of the source terms for tensile and inflation sources, we obtain the solution of the surface static displacements in a layered half-space due to shear, tensile, and inflation point sources. As numerical examples, we investigate three aspects. At first, we test the accuracy of the formulations by comparing our solutions with Okada’s (1985) analytical results in a homogenous half-space and Wang et al.’s (2003) numerical solutions in a multilayered half-space. Secondly, because in most cases we need to calculate the response of finite faulting using a point source approximation, we will discuss the size of the proper subfault area that can be substituted by a point source with various focal depths. Finally, because a uniform half-space Earth model is usually used in inverting geodetic data to estimate the coseismic source mechanism, we are generally concerned how a homogenous half-space model can bias the estimation of fault parameters compared to a layered halfspace model. Here we make a system analysis on the effect of various medium models and source geometries on parameters estimated by static deformation data. To test the accuracy of equations derived in the last section, the static responses of a point source in a homogenous half-space are calculated using our method and compared with Okada’s analytical solutions in Cartesian coordinates. In these numerical tests strike-slip, dip-slip, tensile sources with various dip angles as well as inflation sources are used. Their focal depths are fixed at 5 km. We indicate the discrepancy between our results with Okada’s analytical solutions using the relative error value defined as du uthis paper ⳮ uokada ⳱ , uokada ⬆ 0. u uokada (22) The relative error values of results with different point 2257 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space source models are displayed in Figure 3, in which the open circles indicate relative error values and the crosses are the vertical displacements from Okada’s solutions. Generally, uthis paper would reach to Ⳳ10ⳮ7 cm when uokada equals 0. We also calculate the relative error values of vertical and horizontal displacements adopting various models with various focal depths, dip, and rake angles at different observation points (x, y). In these numerical tests, the maximum discrepancy is less than 2% of the displacement u generally. Furthermore, we compare our point source solutions with Wang et al.’s numerical results in a multilayered half-space. The velocity structure of the layered half-space model is given in Table 1. The discrepancy between our results and Wang et al.’s is presented in Figure 4. Strike-slip and dipslip sources with various dip angles are used. Their focal depths are 5 km. From the figure, we can see that larger discrepancies mainly correspond with lower displacements u and the maximum discrepancy is less than 5% of the displacement usually. The two kinds of comparisons validate our method. 1 0.25 0.15 0 0.05 -1 -20 -10 0 10 0.5 0.6 0 0.4 -0.5 0 0.1 0.05 -0.05 -1 -20 -10 0 10 0.8 0.5 0.6 0.6 0 0.4 -0.5 0.2 -1 0 -1.5 0 -2 -0.2 -20 -10 0 10 0.4 0.06 0.5 0.02 0 -0.02 0.5 0.2 -1 -20 -10 0 10 distance(km) left strike-slip 20 0.15 0 0.05 -0.5 0 10 -0.05 20 -1 -20 -10 0 10 20 1 0.35 1 0.5 0.25 0.5 0 0.15 0 -0.5 0.05 -0.5 -1 -0.05 0.2 -10 0 10 0.3 0.2 -1 -20 20 -10 0 10 20 1 0.35 1 0.5 0.25 0.5 0 0.15 0 -0.5 0.05 -0.5 -1 -0.05 0.1 0 -0.2 -0.5 -0.4 -0.1 -2 -10 0 0 -0.5 -0.06 0.5 -1.5 -20 1 0.25 0.4 20 0.6 1 -1 -20 1 -0.2 1 0.1 -0.5 0.4 1 0.35 0 0.8 0.8 20 1.00 ⬁ 1.2 -0.4 -0.5 0 1.37 3.00 0.5 20 0.5 0.15 1.15 3.58 1 -2 10 Th (km) 1.6 0 0 q (g/cm3) Vp, Vs, q, and Th are the P-wave velocity, S-wave velocity, density, and thickness of each layer, respectively. The values of model-layer rigidity are derived from the velocity models. -1.5 -10 Vs (km/sec) 2.00 6.20 -1 -20 1 Vp (km/sec) 0 -0.2 0.2 Table 1 Velocity Structure 0.2 20 0.25 dip angle 60(degree) 1 -0.5 -0.05 dip angle 90(degree) Uz(cm) 1 0.8 0.5 error(%) dipangle 30(degree) 0.35 Since static deformations decay rapidly with distance from the source, ground deformation modeling is made ordinarily near the fault, so a point source approximation is almost never valid (Lay and Wallace, 1995). In our method we divide the finite fault into several subfaults that are substituted by the point sources. Beyond all doubt, the subfault would be a point source when its area decreases infinitely. However, utilization of the very small area of subfaults will require abundant computer time. The proper size of subfaults is a trade-off between precision and efficiency. In the fol- -0.6 -1 -20 -10 0 10 reverse dip-slip 20 -0.1 -0.2 -20 -10 0 10 20 -1 -20 tensile Figure 3. Comparison of the surface displacements from our result and Okada’s created by point source. The homogenous half-space model is used. From left to right are strike-slip, dip-slip, tensile, and inflation source, respectively. The upper row is for a dip of 30⬚, the middle row is for a dip of 60⬚, and the lower row is for a dip of 90⬚, and their focal depths are 5 km. There are 41 receivers located perpendicular to strike. The plotted crosses are vertical components from Okada’s, and the open circles indicate discrepancy between our results and Okada’s. -10 0 inflation 10 20 Y.-M. He, W.-M. Wang, and Z.-X. Yao 0.4 5 0.3 3 0.2 1 0.1 -1 0 -3 3 0.8 1 -1 0 -5 -20 -10 0 10 -5 -20 -10 0 10 20 1.2 5 3 0.2 -3 -0.4 20 0.3 dip angle 60(degree) 5 0.4 -0.1 5 3 0.8 1 1 0.1 0.4 -1 0 -1 0 -3 -0.1 -5 -20 dip angle 90(degree) Uz(cm) 1.2 -10 0 10 5 0.1 3 1 0 -1 -0.05 -3 -0.1 -0.15 -5 -20 -10 0 10 distance(km) left strike-slip 20 -5 -20 0.15 0.05 -3 -0.4 20 -10 0 10 20 5 0.8 3 0.4 error(%) dip angle 30(degree) 2258 1 0 -1 -0.4 -3 -0.8 -5 -20 -10 0 10 20 reverse dip-slip Figure 4. Comparison of the surface displacements from our result and Wang et al.’s in a layered half-space. From left to right are strike-slip and dipslip source, respectively. The plotted crosses are vertical components from Wang et al.’s, and the open circles indicate the discrepancy between our and Wang et al.’s results. The other arrangement of the figure is same as that in Figure 3. lowing numerical test, we investigate the relationship between the area of subfault and focal depth and source mechanism. Our results are based on a point source, and Okada’s analytical solutions are from a planar rectangular fault. By analysis of the discrepancy between our results and Okada’s analytical solutions, we can estimate the lower limit of the subfault area. Various subfaults, focal depths, and source mechanisms are investigated. Our study suggests that the proper area of the subfault would decrease with shallower focal depth. This result is also displayed in Figure 5, in which various symbols mark different source depths. From this figure we can see clearly that the discrepancy increases with the larger subfault area and decreases with the deeper focal depth. All these facts suggest the 1 km ⳯ 1 km subfault is small enough when the focal depth is greater than 3 km. At shallower depth, we need to divide the finite fault into much smaller subfaults such as 0.3 km ⳯ 0.3 km for a focal depth of 1 km. The development of this article is the derivation of the source terms for tensile and inflation point sources, as well as the use of a layered half-space model. We will now discuss how a homogenous half-space model can bias the estimation of fault parameters than a layered half-space model. In the numerical experiments we calculate the static displacement as “data” based on a prescribed source in a layered half-space model, then deduce its source parameters by fitting this data based on a homogenous half-space model. In the forward calculation, we set up a 6 km ⳯ 4 km planar rectangle fault buried in a layered half-space and divided the finite fault into 24 1 km ⳯ 1 km subfaults (see Fig. 6). Each subfault has the same slip dislocation of 1 m. Substituting point source for subfault, we obtain the vertical and horizontal ground displacement data by our method. In the inversion we utilize a hybrid global search algorithm (Liu et al., 1995) and calculate the responses using Okada’s analytical solutions of a planar rectangular fault. In Okada’s code a single 6 km ⳯ 4 km fault with one slip displacement is applied. The strike, length, and width of the fault are fixed, and only the focal depth, slip displacement, rake, and dip angles are inverted. Nevertheless, before our investigation of the effect of the layered model, we need to test the efficiency of our forward and inverse methods first. We construct the data using the homogenous half-space model and our method and invert it with Okada’s solutions. Various source mechanisms are investigated. Table 2 shows that in any cases the inverted focal depth, slip displacement, rake, and dip angels are quite close to the real values. All these facts validate the reliability of our forward and inverse method. Hereafter we will probe the effect of medium model. Several typical fault models and corresponding velocity structure in reality are investigated. The general structures for a strike-slip fault (e.g., San Andreas fault) and reverse dip-slip earthquake (e.g., Chi-Chi earthquake), the typical rift structure (e.g., Corinth Gulf) for normal and tensile faulting earthquakes, and the structure of an active volcano (e.g., Mt. Etna) for inflation sources are selected. The layered earth models of different earthquake zones are given in Tables 3, 4, 5, and 6, respectively. For all kinds of designed fault models, the focal depth, rake, and slip displacement are fixed and the dip angle is changed from 30⬚ to 90⬚. The final inverse results including focal depth, slip displacement, dip, and rake angles are given in Table 7. Since the 1989 Loma Prieta earthquake ruptured a part of the San Andreas fault, we choose the Loma Prieta velocity structure model (Table 3; Wald, et al. 1991) to discuss the strike-slip fault. From Table 7 we observe that the inverted rake and dip angles almost match the “true” values of the layered models. The two parameters sensitive to the existence of the layered half-space model are the focal depth and slip displacement. The largest discrepancy with the real value in focal depth and slip can reach 12% and 7%, respectively. For reverse dip-slip model, we select the layered central Taiwan Earth model (Table 4) given by Ma et al. 2259 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space 8 40 strike slip focal depth of 1km focal depth of 2km focal depth of 3km focal depth of 10km 6 dip slip focal depth of 1km focal depth of 2km focal depth of 3km focal depth of 10km 30 20 tensile 16 focal depth of 1km focal depth of 2km focal depth of 3km focal depth of 10km (%) 12 4 20 8 2 10 0 4 0 0 0.2 0.4 0.6 0.8 width of square finite fault(km) 1 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figure 5. The discrepancy of maximum vertical displacements created by different finite fault area with that by point source. Various focal depths are investigated and illustrated by different symbols. Here the horizontal axis gives the width of square finite fault, and the vertical axis show the discrepancy. x 5km y 1km Table 2 1km A point source z Figure 6. 3D dislocation model in a layered halfspace used to calculate the coseismic displacements with our method. (1996). The character of the inverse parameters is similar to that of the strike-slip source when the dip angle equals 30⬚. The largest discrepancy in depth and slip can reach 16% and 4%, respectively. Unexpectedly, the inverted focal depth and slip are almost unchanged when we investigate the reverse dip-slip source with dip angles greater than 50⬚. As an active extension regime, the refined velocity structure of Corinth Gulf (Table 5; Zahradnik, 1999) is applied to the normal dip-slip source. Different from the inverted results for strike-slip and reverse dip-slip sources, the results presented in Table 7 show the dip angle is also sensitive to the presence of the layered half-space model as well as focal depth and slip. Moreover, the slip is the most sensitive parameter, and the maximum discrepancy in slip can reach 27%. More numerical tests show that the sensitive slip is mainly due to the ultralow-velocity topsoil. The discrep- Real and Inverted Results using Homogenous Half-Space Model for 4 by 6 Subfault Plane. Focal Depth (km) Dip (⬚) Rake (⬚) Slip (m) Left strike-slip Real Inverted Real Inverted Real Inverted 5.00 5.28 5.00 5.05 5.00 5.01 50.00 51.98 70.00 70.05 90.00 90.01 0.00 0.04 0.00 0.02 0.00 0.02 1.00 1.05 1.00 1.01 1.00 1.00 Reverse dip slip Real Inverted Real Inverted Real Inverted 5.00 5.05 5.00 4.99 5.00 5.01 30.00 30.07 50.00 50.12 70.00 70.17 90.00 90.05 90.00 90.52 90.00 90.32 1.00 1.00 1.00 1.01 1.00 1.01 Normal dip slip Real Inverted Real Inverted Real Inverted 5.00 4.99 5.00 4.97 5.00 4.94 30.00 29.88 50.00 50.01 70.00 69.96 ⳮ90.00 ⳮ89.76 ⳮ90.00 ⳮ90.06 ⳮ90.00 ⳮ89.95 1.00 1.00 1.00 1.00 1.00 1.00 Tensile Real Inverted Real Inverted Real Inverted 5.00 4.99 5.00 5.01 5.00 5.04 50.00 50.05 70.00 69.99 90.00 89.88 – – – – – – 1.00 1.00 1.00 1.00 1.00 1.00 Inflation Real Inverted 5.00 4.90 – – – – 1.00 0.99 The forward calculation used our method. Okada’s (1985) solution is applied to invert focal depth, dip, rake, and slip. For inflation source, only source depth and volume change are investigated. 2260 Y.-M. He, W.-M. Wang, and Z.-X. Yao Table 3 Table 6 Loma Prieta Velocity Structure Velocity structure of Mt. Etna, Italy Vp (km/sec) Vs (km/sec) q (g/cm3) Th (km) Vp (km/sec) Vs (km/sec) q (g/cm3) Th (km) 1.73 3.38 4.29 4.80 5.37 5.74 6.15 6.25 6.27 6.67 8.00 1.00 1.95 2.48 2.77 3.10 3.31 3.55 3.61 3.62 3.85 4.62 1.50 1.55 1.85 2.05 2.26 2.45 2.58 2.62 2.63 2.77 3.28 0.1 0.4 0.5 2.0 2.0 2.0 2.0 4.0 5.0 7.0 50.0 3.50 4.00 4.50 5.00 6.00 6.50 8.05 2.02 2.31 2.60 2.89 3.46 3.75 4.65 2.06 2.19 2.33 2.46 2.73 2.87 3.30 2.0 5.0 3.0 3.0 3.0 15.0 50.0 Model from Wald et al. (1991). Table 4 Central Taiwan Crustal Model Vp (km/sec) Vs (km/sec) q (g/cm3) Th (km) 3.50 3.78 5.04 5.71 6.05 6.44 6.83 2.00 2.20 3.03 3.26 3.47 3.72 3.99 2.0 2.3 2.5 2.6 2.6 2.6 3.0 1.0 3.0 5.0 4.0 4.0 8.0 6.0 Model from Ma et al. (1996). Table 5 Corinth Gulf Crustal Model q (g/cm3) Th (km) 0.82 1.54 2.57 3.29 3.70 4.56 1.49 1.83 2.31 2.65 2.84 3.25 0.5 0.5 1.0 3.0 34.0 50.0 Modified Model 1.54 2.57 3.29 3.70 4.56 1.83 2.31 2.65 2.84 3.25 1.0 1.0 3.0 34.0 50.0 Vp (km/sec) Vs (km/sec) 1.42 2.67 4.45 5.70 6.40 7.90 2.67 4.45 5.70 6.40 7.90 Model Model from Zahradnik (1999). ancy in slip would decrease when the P-wave velocity of the topsoil increases (see Table 5 and Table 7: inverted b). The same velocity structure of the Corinth Gulf is also used for the tensile source. Table 7 shows that the parameter most sensitive to the layered half-space model is the focal depth, and its discrepancy ranges from 18% to 33%. The modified velocity struc- Modified from Villasenor et al. (1998). ture of Mt. Etna, Italy (Table 6; Villasenor et al., 1998) is employed for inflation source where rake and dip is not exit and slip displacement is substituted by volume change. The source depth decreases to 4.19 km and the volume changes to 0.91. Another set of inversions is made with only the horizontal displacements since these are often the only components available. The results shown in Table 8 are similar to those in Table 7. Discussions and Conclusion Over the last two decades, several workers have investigated the effect of Earth layering on fault displacements for both strike-slip and dip-slip sources (e.g., Sato, 1971; Sato and Matsu’ura, 1973; Savage, 1987; Ma and Kusznir, 1994; Savage, 1998; Cattin et al., 1999). Most of them simply constructed the model of a superficial layer overlying a half-space to represent the layered half-space model. In this article, based on the generalized reflection and transmission coefficient matrix method, formulations for surface static displacements in a layered half-space are extended to include tensile and inflation point sources from a point pure shear dislocation source. We also select several typical velocity structures in reality and various source mechanisms including tensile and inflation and obtain newly quantified results. Generally speaking, using the layered half-space model, the calculated ground displacements would be greater (less) than that applying homogenous half-space model depending mainly on rigidity increase (decrease) with depth (Ma and Kusznir, 1994). In our article, we choose several velocity structures as approximations of the Earth. Without question the rigidity would increase with depth in these models and the calculated displacements would be amplified. For fitting the higher displacements with a homogenous half-space model, one has to decrease the depth and change the slip and alter the dip angle sometimes. Our results show, for strike-slip sources and reverse dip-slip sources with a dip of 30⬚, that the focal depth is the main sensitive parameter to the stratified media. The neglect of the layering, especially the low-velocity layers at shallower depth, would cause one to underestimate the depth of the dislocation (Cervelli et al., 2002). Two sets of inversions, the first with the vertical and 2261 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space Table 7 Table 8 Real and Inverted Results using Layered and Homogenous HalfSpace Model for 4 by 6 Subfault Plane, Respectively Real and Inverted Results using Layered and Homogenous HalfSpace Model for 4 by 6 Subfault Plane Focal Depth (km) Dip (⬚) Rake (⬚) Slip (m) Focal Depth (km) Dip (⬚) Rake (⬚) Slip (m) Left strike slip* Real Inverted Real Inverted Real Inverted 5.00 4.64 5.00 4.40 5.00 4.44 50.00 50.12 70.00 69.51 90.00 90.23 0.00 ⳮ0.16 0.00 0.00 0.00 ⳮ0.02 1.00 1.04 1.00 1.06 1.00 1.07 Left strike slip Real Inverted Real Inverted Real Inverted 5.00 4.40 5.00 4.50 5.00 4.37 50.00 49.48 70.00 69.82 90.00 90.04 0.00 0.02 0.00 ⳮ0.01 0.00 ⳮ0.06 1.00 1.07 1.00 1.06 1.00 1.04 Reverse dip slip† Real Inverted Real Inverted Real Inverted 5.00 4.22 5.00 4.94 5.00 4.85 30.00 29.01 50.00 50.33 70.00 70.04 90.00 90.14 90.00 90.07 90.00 89.87 1.00 1.04 1.00 1.02 1.00 1.01 Reverse dip slip Real Inverted Real Inverted Real Inverted 5.00 4.27 5.00 4.95 5.00 4.86 30.00 29.42 50.00 50.15 70.00 69.89 90.00 90.13 90.00 90.13 90.00 89.89 1.00 1.11 1.00 1.01 1.00 1.03 Normal dip slip‡ Real Inverted a Inverted b Real Inverted a Inverted b Real Inverted a Inverted b 5.00 4.66 4.68 5.00 4.41 4.47 5.00 4.59 4.57 30.00 31.76 31.65 50.00 48.43 48.98 70.00 68.25 68.22 ⳮ90.00 ⳮ89.82 ⳮ90.03 ⳮ90.00 ⳮ89.92 ⳮ89.68 ⳮ90.00 ⳮ89.90 ⳮ90.09 1.00 1.15 0.99 1.00 1.11 1.05 1.00 0.73 1.00 Normal dip slip Real Inverted a Inverted b Real Inverted a Inverted b Real Inverted a Inverted b 5.00 4.68 4.80 5.00 4.47 4.76 5.00 4.90 5.04 30.00 31.85 31.72 50.00 48.20 49.83 70.00 67.26 65.85 ⳮ90.00 ⳮ90.06 ⳮ89.68 ⳮ90.00 ⳮ89.81 ⳮ90.00 ⳮ90.00 ⳮ89.81 ⳮ89.65 1.00 1.15 1.14 1.00 1.21 1.12 1.00 0.53 1.04 Tensile§ Real Inverted Real Inverted Real Inverted 5.00 4.10 5.00 4.06 5.00 3.34 50.00 47.35 70.00 69.09 90.00 90.02 – – – – – – 1.00 0.98 1.00 1.10 1.00 1.01 Tensile Real Inverted Real Inverted Real Inverted 5.00 3.35 5.00 4.07 5.00 3.33 50.00 43.66 70.00 69.15 90.00 90.09 – – – – – – 1.00 0.90 1.00 1.12 1.00 1.01 Inflation| Real Inverted 5.00 4.19 – – – – 1.00 0.91 Inflation Real Inverted 5.00 4.19 – – – – 1.00 0.91 Vertical and horizontal ground displacements are used simultaneously. * The layered parameters are shown in Table 3. † The layered parameters are shown in Table 4. ‡ The layered parameters are shown in Table 5. § The layered parameters are shown in the top half of Table 5. | The layered parameters are shown in Table 6. horizontal displacements and the second with only the horizontal displacements, have similar results. However, the effects of the layered half-space model are also varied with different source mechanisms. For example, for a reverse dipslip source with steep dip angle (⬎50⬚), the effect of the layered half-space model is small and can be ignored. Cattin et al. (1999) analyzed the effect of a layer at the top of the upper crust on a shallow normal-faulting earthquake using analytical and numerical modeling. They concluded that using a homogenous half-space to interpret the coseismic displacements leads to shallower faults than in reality in all cases and tends to overestimate the slip by 10%– 20% when horizontal displacements alone are considered. For normal dip-slip source, we found the ultralow-velocity topsoil would cause slip much sensitive even though the vertical and horizontal components are used simultaneously. When the P-wave velocity of the topsoil increases to a cer- Horizontal ground displacements alone are considered. The layered parameters are same as table 7. tain value (e.g., 2.67 km/sec), our results agree with Cattin et al.’s (1999) roughly. However, similar with Ji (2002), our results show that the dip angle is also sensitive to the nonuniformity of the medium, which differs from the results of Savage (1998) and Cattin et al. (1999). In conclusion, in most cases the neglect of the layered half-space model, especially the low-velocity layers at shallower depth, would cause one to underestimate the depth of a source. Moreover, using a homogenous half-space to simulate the coseismic displacements also leads to a misestimate of the slip and dip angle in the case of tensile and normal faulting. Only for the reverse dip-slip source with steep dip angle (⬎50⬚) can one neglect the layered half-space model. Acknowledgments We thank Dr. Xiao-Bi Xie and Prof. Tian-yu Zheng for helpful discussions. We are also grateful to Dr. R. Bürgmann and two reviewers for valuable comments. This research is supported by the National Science Foundation of China (49634150) and the National Key Basic Research Program (1999433). 2262 Y.-M. He, W.-M. Wang, and Z.-X. Yao References Ben-Menahem, A., and S. J. Singh (1968). Multi-polar elastic fields in layered half-space, Bull. Seism. Soc. Am. 58, 1519–1572. Cattin, R., P. Briole, H. Lyon-Caen, P. Bernard, and P. Pinettes (1999). Effects of superficial layers on coseismic displacements for a dip-slip fault and geophysical implications, Geophys. J. Int. 137, 149–158. Cervelli, P., P. Segall, K. Johnson, M. Lisowski, and A. Miklius (2002). Sudden aseismic fault slip on the south flank of Kilauea volcano, Nature 415, 1041–1018. Hearn, E. H., R. Bürgmann, and R. E. Reilinger (2002). Dynamics of İzmit earthquake postseismic deformation and loading of the Düzce earthquake hypocenter, Bull. Seism. Soc. Am. 92, 172–193. Ji, C. (2002). Analysis of complex faulting: wavelet transform, multiple datasets, and realistic fault geometry, Ph. D. Thesis, Caltech, Pasadena, California. Ji, C., V. H. Donald, T. A. Song, K. F. Ma, and D. J. Wald (2001). Slip distribution and tectonic implication of the 1999 Chi-Chi, Taiwan, Earthquake, Geophys. Res. Lett. 28, 4379–4382. Kennett, B. L. N. (1974). Reflections, rays, and reverberations, Bull. Seism. Soc. Am. 64, 1685–1696. Kennett, B. L. N. (1983). Seismic Wave Propagation in Stratified Media, Cambridge U Press, New York. Lay, T., and T. C. Wallace (1995). Modern Global Seismology, Second Ed., Academic, New York. Liu, P., S. Hartzell, and W. Stephenson (1995). Non-linear multiparameter inversion using a hybrid global search algorithm: applocation in reflection seismology, Geophys. J. Int. 122, 991–1000. Luco, J. E., and R. Apsel (1983). On the Green’s functions for a layered half-space. I. Bull. Seism. Soc. Am. 73, 909–929. Ma, K. F., J. H. Wang, and D. P. Zhao (1996). Three-dimensional seismic velocity structure of the crust and uppermost mantle beneath Taiwan, J. Phys. Earth 44, 85–105. Ma, X. Q., and N. J. Kusznir (1994). Effects of rigidity layering, gravity, and stress relaxiation on 3-D subfurface fault displacement fields, Geophys. J. Int. 118, 201–220. Ma, X. Q., and N. J. Kusznir (1995). Coseismic and postseismic subsurface displacements and strains for a dip-slip normal fault in a three-layer elastic-gravitational medium, J. Geophys. Res. 100, 12,813–12,828. Maruyama, T. (1964). Static elastic dislocations in an infinite and semiinfinite medium, Bull. Earthquake Res. Inst. Tokyo Univ. 42, 289– 368. Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space, Bull. Seism. Soc. Am. 75, 1135–1154. Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space, Bull. Seism. Soc. Am. 82, 1018–1040. Roth, F. (1990). Subsurface deformation in a layered elastic half-space, Geophys. J. Int. 103, 147–155. Sato, R. (1971). Crustal deformation due to dislocation in a multi-layered medium, J. Phys. Earth 19, 31–46. Sato, R., and M. Matsu’ura (1973). Static deformations due to the fault spreading over several layers in a multilayered medium. I. Displacement, J. Phys. Earth 21, 227–249. Savage, J. C. (1987). Effect of crust layering upon dislocation modeling, J. Geophys. Res. 92, 10,595–10,600. Savage, J. C. (1998). Displacement field for an edge dislocation in a layered half-space, J. Geophys. Res. 103, 2439–2446. Singh, S. J. (1970). Static deformation of a multi-layered half-space by internal sources, J. Geophys. Res. 75, 3257–3263. Steketee, J. A. (1958). On Volterra’s dislocation in a semi-infinite elastic medium, Can. J. Phys. 36, 192–205. Villasenor, A., H. M. Benz, L. Filippi, G. De Luca, R. Scarpa, G. Patanè, and S. Vinciguerra (1998). Three-dimensional P-wave velocity structure of Mt. Etna, Italy, Geophys. Res. Lett. 25, 1975–1978. Wald, D. J., D. V. Helmberger, and T. H. Heaton (1991). Rupture model of the 1989 Loma Prieta earthquake from the inversion of strong motion and broadband teleseismic data, Bull. Seism. Soc. Am. 81, 1540–1572. Wang, R., F. L. Martin, and F. Roth (2003). Computation of deformation induced by earthquakes in a multi-layered elastic crust: FORTRAN programs EDGRN/EDCMP. Comp. Geosci. 29, 195–207. Xie, X. B., and Z. X. Yao (1989). A generalized reflection-transmission coefficient matrix method to calculate static displacement field of a stratifield half-space by dislocation source, Chinese J. Geophys. 32, 191–205. Yao, Z. X., and D. G. Harkrider (1983). A generalized reflection transmission coefficient matric and discrete wavenumber method for synthetic seismograms, Bull. Seism. Soc. Am. 73, 1685–1699. Zahradnik, H. (1999). Focal mechanisms from the Corinth Gulf BB stations, final report for the Inco-Copernicus COME project, Charles University, Prague. Zhu, L., and L. Rivera (2002). A note on the dynamic and static displacements from a point source in multilayered media, Geophys. J. Int. 148, 619–627. Appendix Static Reflection and Transmission Coefficient Matrix For the P-SV problem, the downgoing and upgoing static reflection and transmission coefficient matrixes are defined as RD ⳱ 冢rr Ⳮ pp Ⳮ ps Ⳮ Ⳮ rsp tⳭ tsp pp Ⳮ , TD ⳱ Ⳮ Ⳮ , rss tps tss RU ⳱ rⳮ pp rⳮ ps ⳮ ⳮ rsp tⳮ tsp pp ⳮ , TU ⳱ ⳮ ⳮ , rss tps tss 冢 冣 冣 冢 冢 冣 冣 (A1) where rⳭ pp ⳱ 2D1kd1(l1 ⳮ l2)/(l1D1 Ⳮ l2), 2 2 rⳭ sp ⳱ ⳮ[4D1k d1(l1 Ⳮ l2D2)(l1 ⳮ l2) Ⳮ (l1D1 Ⳮ l2)(l1D1 ⳮ l2D2)]/(l1 Ⳮ l2D2)(l1D1 Ⳮ l2), Ⳮ rps ⳱ (l2 ⳮ l1)/(l1D1 Ⳮ l2), rⳭ ss ⳱ 2D1kd1(l2 ⳮ l1)/(l1D1 Ⳮ l2), (A2) tⳭ pp ⳱ l1(1 Ⳮ D1 Ⳮ D2 Ⳮ D1D2)/(1 Ⳮ D2)(l1D1 Ⳮ l2), tⳭ sp ⳱ 2l1D1kd1(1 Ⳮ D1 Ⳮ D2 Ⳮ D1D2)/(1 Ⳮ D2)(l1D1 Ⳮ l2), tⳭ ps ⳱ 0, tⳭ ss ⳱ l1(1 Ⳮ D1 Ⳮ D2 Ⳮ D1D2)/(1 Ⳮ D2)(l1D1 Ⳮ l2), (A3) rⳮ pp ⳱ 0, rⳮ sp ⳱ (l1D1 ⳮ l2D2)/(l1D1 Ⳮ l2), rⳮ ps ⳱ (l1 ⳮ l2)/(l1 Ⳮ l2D2), rⳮ ss ⳱ 0, (A4) tⳮ pp ⳱ l2(1 Ⳮ D2)/(l1 Ⳮ l2D2), tⳮ sp ⳱ 2D1kd1 l(1 Ⳮ D2)/(l1D1 Ⳮ l2), tⳭ ps ⳱ 0, tⳭ ss ⳱ l2(1 Ⳮ D2)/(l1D1 Ⳮ l2). (A5) 2263 Static Deformation Due to Shear and Tensile Faults in a Layered Half-Space Here subscripts 1 and 2 denote parameters in upper and lower media, d is the thickness of the upper medium, and k is the wavenumber. For the SH problem, the static reflection and transmission coefficients are RD,L ⳱ (l1 ⳮ l2)/(l1 Ⳮ l2), TD,L ⳱ 2l1 /(l1 Ⳮ l2), RU,L ⳱ (l2 ⳮ l1)/(l1 Ⳮ l2), TU,L ⳱ 2l2 /(l1 Ⳮ l2). of phase in a static problem, the existence of a layer with a width of d1 still distorts static displacement fields on adjacent layers. We use the amplitude decay to describe the phenomenon. The relationship between the static generalized reflection and transmission coefficient matrix at the tops of two neighboring layers and the static reflection and transmission coefficients at an interface is Ⳮ ˜Ⳮ Ⳮ ⳮ ⳮ ˜ⳮ ⳮ r̃Ⳮ pp ⳱ exp(ⳮ2kd1)rpp, t pp ⳱ exp(ⳮkd1)tpp, r̃pp ⳱ rpp, t pp ⳱ exp(ⳮkd1)tpp, (A6) Ⳮ ˜Ⳮ Ⳮ ⳮ ⳮ ˜ⳮ ⳮ r̃Ⳮ sp ⳱ exp(ⳮ2kd1)rsp , t sp ⳱ exp(ⳮkd1)tsp , r̃sp ⳱ rsp , t sp ⳱ exp(ⳮkd1)tsp , Ⳮ ˜Ⳮ Ⳮ ⳮ ⳮ ˜ⳮ ⳮ r̃Ⳮ ps ⳱ exp(ⳮ2kd1)rps , t ps ⳱ exp(ⳮkd1)tps , r̃ps ⳱ rps , t ps ⳱ exp(ⳮkd1)tps , Ⳮ ˜Ⳮ Ⳮ ⳮ ⳮ ˜ⳮ ⳮ r̃Ⳮ ss ⳱ exp(ⳮ2kd1)rss , t ss ⳱ exp(ⳮkd1)tss , r̃ss ⳱ rss , t ss ⳱ exp(ⳮkd1)tss . Static Reflection and Transmission Coefficient Matrix with Amplitude Decay For the elastic problem, the generalized reflection and transmission coefficient matrix at the tops of two neighboring layers is defined as reflection and transmission coefficient matrix with phase delay. Although there is no concept (A7) Institute of Geology and Geophysics Chinese Academy of Science P.O. Box 9825 Beijing 100029, P.R. China Manuscript received 19 June 2002.