Coupling of Focal Mechanism with finite Difference Wave Propagation in Rotated Staggered GridJunlun LI Department of Earth, Atmospheric and Planetary Sciences, MITIntroductionFinite difference scheme has been widely used to solve elastic wave propagation in heterogeneous and layered media. Traditionally, staggered grid scheme is used for differencing the equations of Newton’s Motions Law and Generalized Hook’s Law describing the elastic wave propagation. In the two dimensional case, the equations can be expressed as follow: )()2()2(xvzvtzvxvtzvxvtzxtvzxtvzxxzzxzzzxxxzzxzzxzxxx∂∂+∂∂=∂∂∂∂++∂∂=∂∂∂∂+∂∂+=∂∂∂∂+∂∂=∂∂∂∂+∂∂=∂∂µτµλλτλµλτττρττρ Where xv, zvis velocity component, ijτis stress component, ρis density, µλ,are Lame constant. For traditional staggered grid method, we would discretize the above equations as in the following scheme:Fig. 1. Traditional Staggered grid schemeBy placing as above, the velocity components and stress components are in staggered grids. However, when we update velocity components xvand zvand shear stressxzσ, we have to use averaged density and shear modulus. For density we simply use arithmetic average and it is fine when wave transmitted from solid into vacuum or some medium with very low density. However, when we average shear modulus µ, we should use harmonic averaging: )1,(1)1,1(1),1(1),(1)21,21(4−+−−+−+=−−jijijijijiµµµµµ Therefore, when wave transmitted from solid into some medium with zero or very small shear modulus, at the boundary point the equation will blow up.In order to solve this problem, Saenger et. al. proposed another approach called xρzvρxvρxvzzxxσσµλ,,xzσµzρzvrotated staggered grid. Their basic idea is to rotate the coordinate system that all the stresses are located at the center of the cell. In doing so, when stresses are updated, no averaging among the neighboring points is needed. Their algorithm can be expressed in the following plot: Fig. 2. Rotated staggered grid schemeThis new rotated staggered grid can be related to the traditional staggered grid with the following coordinate transforms:22~~xzrzrzxrxxzrzxrxz∆+∆=∆∆∆−∆∆=∆∆+∆∆=In traditional coordinate system, we define the following 2nd order operator as:~z~xρzxvvρzxvvρzvρxvρzxvvxzzzxxσσσµλ,,,ρzxvv∆−−∆+∆=∆−−∆+∆=),2,(),2,(1),,(),,2(),,2(1),,(tzzxutzzxuxtzxuDtzxxutzxxuxtzxuDzxThen we can get the 2nd order operator in the new rotated coordinate system as:∆−∆−−∆+∆+∆=∆+∆−−∆−∆+∆=),2,2(),2,2(1),,(),2,2(),2,2(1),,(~~tzzxxutzzxxuxtzxuDtzzxxutzzxxuxtzxuDzxThen we transform them back into the traditional coordinate system:+−∆∆≈∂∂+∆∆≈∂∂),,(),,(2),,(),,(),,(2),,(~~~~tzxuDtzxuDzrtzxuztzxuDtzxuDzrtzxuxzxzxBy doing so, we can obtain a new set of discretized equations for describing the elastic wave propagation.Thanks to Cruz-Atienza etc., we could have an easier insight into this coordinate rotation approach. Besides the above expression of algorithm, we could define )(21)()(21)(2/1,2/12/1,2/12/1,2/12/1,2/12/1,2/12/1,2/12/1,2/12/1,2/1+−−+−−+++−−+−−+++−−=−+−=jijijijiijzjijijijiijxffffhfDffffhfDAnd if we discretize the original differential equations by these differencing operators, we could have the same difference equations as derived from above.Implementations1). Air Wave PropagationThe whole codes are written in Matlab. Perfectly Matched Layers (PML) with 40 grids in thickness are placed on the four boundaries. For simulation, we use grid size dh=2 m. The whole area is 240x300 grids. There are two layers in this model: from grid 1 to 60 is air, which has a density of 1.292 kg/m^3, Vp=340 m/s, and shear modulus u=0; from grid 61 to 300 is rock, which has a density of 3000 kg/m^3, Vp=3464 m/s, Vs =2000 m/s. From the density, Vp and Vs, we can calculate the Lame constants. Although we now only deal with two layers, the numerical scheme actually can handle heterogeneous media as later we will see.For simulation wave propagation, an isotropic explosion source is used, which is added to stresses zzxxσσ, within a 3x3 square area. As pointed out by Cruz-Atienza etc., this scheme could be a better approximation of point source than just using one single point. The source center is located at grid (120, 100). The source time function is a Ricker wavelet, which in time domain can be expressed as below:)2/exp()1(2222,tataAzzxx−−=σσWhere A is amplitude, a=π20⋅f,0f is the central frequency. The wavelet is demonstrated below:0 0.05 0.1 0.15-50510x 104timeamplitude Fig. 3. Wavelet for Simulating Air Wave propagationThe velocity distributions at t=0.120 s are plotted as belowFig. 4. Air Wave Propagation: first row is Vz and second row is VxFrom the figures above we can clearly see that this numerical scheme works well for high contrast media. The wave transmitted into air can be correctly modeled. No Vx component exists in the air as shear wave can not propagate in the air, whose shear modulus u is zero.One thing should be noticed in the implementations of this paper is that we flip the source time function by π. Due to possible different definition in coordinate system, or possible mistakes in staggered grid rotation or Matlab index usage, we find the phases in our results should be reversed by π to match various benchmark tests. This reverse is applied in source time functions in this paper.2). Single Arbitrarily Oriented Fault SourceRealistic seismic fault is a finite length rupture, ranging from hundreds of meters to hundreds of kilometers. Traditionally, people would use single point source to approximate a fault if observation point is far enough, e.g., tens of times of the fault length away from the fault. However, for a large seismic fault like San Andreas, we should use finite length fault for modeling. In order to model finite length, arbitrarily oriented fault, we first model a single point, arbitrarily oriented “fault”.Fig. 5 . San Andreas Fault Source: http://www.physics.unlv.edu/~jeffery/astro/earth/geology/plate/usgs_020_sanandreas.gifAfter Cruz-Atienza etc., we model this problem as