Modeling seismic wave propagation in a coal-bearing porous medium by a staggered-grid finite differe

Mining Science and Technology( China)21(2011)727-731Contents lists available at SciVerse Science DirectMining Science and Technology( China)ELSEVIERjournalhomepagewww.elsevier.com/locate/mstcModeling seismic wave propagation in a coal-bearing porous mediumby a staggered-grid finite difference methodZou Guangui", Peng Suping, Yin Caiyun, Deng Xiaojuan, Chen Fengying, Xu YanyongState Key Laboratory of Coal Resources and Mine Safety, China University of Mining g Technology, Beijing 100083, ChinaARTICLE INFOA BSTRACTReceived 1o januarv 2011A staggered-grid finite difference method is used to model seismic wave records in a coal bearing.medium. The variables analyzed include the order of the difference calculations, the use of amatch layer to provide absorbing boundary conditions, the source location, the stability condnditions.Accepted 12 March 2011Available online 6 November 2011and dispersion in the medium the results show that the location of the first derivative of the dynamicvariable with respect to space is coincident with the location of the first derivative of the kinematic variable with respect to time. Outgoing waves are effectively absorbed and reflection at the boundary is veryweak when more than 20 perfect match layer cells are used. Biot theory considers the liquid phase to beorous mediaus so the ratio of liquid to solidStaggered-grid finite differencerosity Numerical dispersion and generation of false frequencies is reduced by increasing theof the difference calculations and by reducing the grid size and time step. Temporal second orracy, a tenth order spatial accuracy, and a wavelength over more than ten grid points gave acceptablenumerical results. Larger grid step sizes in the lateral direction and smaller grid sizes in the vertical direion allow control of dispersion when the medium is a low speed body this provides a useful way to sim-ulate seismic waves in a porous coal bearing mediumo 2011 Published by Elsevier B V on behalf of China University of Mining Technology1 IntroductionP-wave equation in a porous medium and then applied waveletmulti-scale simulation to these P-wave equations[8 Pei unfoldeda porous medium comprises a solid matrix and a fluid within the high order staggered-grid finite difference when simulatingthe pores. Examples are coal with gas in the coal pores, limestone elastic wave propagation in a three dimensional anisotropic porouswith water in the karst pores, or sandstone with oil in the pores. medium [9]The medium is traditionally considered as a single phase but a por-The actual medium for the seismic propagation is a semious medium with two phases is more consistent with the actual infinite space but a numerical simulation must be constructed assituation. Understanding the characteristics and properties of wave a finite geological model. Boundary reflections then occur in thepropagation in complex media, like porous ones, require the use of simulated wave propagation. many scholars have worked to elim-numerical simulation, a basic research approach. The staggered- inate these boundary reflections. The perfect match layer(PML)grid finite difference method is considered to have higher precision was used to absorb boundary reflections by Collin et al who ap-and to require less computation when modeling such media(1-2). plied finite difference modeling( FDM)and finite element modelingUsing Biot theory Guo et al. considered the P-wave response in (FEM)to several examples [101. Zeng gave the PMl absorbingporous media and proposed a finite difference algorithm as a solu boundary conditions for a porous medium [111on to this problem 3. Dai et al. applied a finite difference meth-Coal-bearing strata are at a shallow depth and show large propod to obtain the speed versus stress equation for a porous medium erty differences with respect to deep rock: The propagation speed[4]. Carcione et al. performed a numerical simulation of the P-wave and material density, for example. Seismic wave propagationin a Biot poro-elastic medium [5]. Ozdenvar et al. used the stag. coal-rock is also distinct from that seen in an oil fieldgered-grid finite difference method to simulate a seismic wave inSo the question is how wave propagation may be numericallyfinite difference method for the P-wave equations in porous isotro- a careful attention Soal-bearing porous media: A problem needingpic media[7]. Liu used a matrix to express the finite difference of aWe apply herein a staggered-grid finite difference approach tothe solution of the stress-velocity elastic wave equations in a porous medium. The variables are located based on the principles ofCorresponding author. Tel. +86 10 62331305.E-mailaddress:cumtzgg@foxmail.com(GZou).the staggered-grid finite difference. The PML absorbing boundary中国煤化工1674-5264/5- see front matter o 2011 Published by Elsevier B V on behalf of China University of Mining Techdoi:10.1016mst201103.002YHCNMHG728G Zou et aL/ Mining Science and Technology(China) 21(2011)727-731conditions are used to eliminate any boundary reflections. Simu-lated seismic records are analyzed under different parametersincluding boundary conditions, source settings, stability condi-tions, and dispersion. Finally, a simulated example is shown as◆-an example.2. Staggered-grid, finite difference method for modeling theelastic wave equation in porous mediumThe geometric equation of strain-displacement in the porousmedium is cast as a generalized Hooke ' s law. then this result iscombined with a generalized Darcys law to give the elastic waveFig. 1. Staggered-grid diagram: left,equation in a porous medium when the time derivative of the dis- time.it, space gridplacement is given by the speed the results shown in Eqs. (1)-(8)areobtained(Pfa-PuPa) Dr=Pu ax Pa( x+82then vz and i are set at location 4 by (i+0.5.j+0.5. n). Then ooax+6z +(pu+pa (b 4) om, and s are set at location 2 by(i+0.5. j, n+0. 5)and tx is set(Pi2-PuP2l8t-Pu2 BZ Pz( ax+82)+(Pu2+Pu)b 6- (u) at location 3 by (i,j+0.5. n+0.5). The locations of the othervariables may be deduced by analogy.Second, the first derivative is converted to a difference using the(2) scheme)+a(-C:)(3pn(,)++)h(4) where the subscript of variable U is the spatial location, the superC+C亚+Q(+亚(5) script is the time location, Lx the x derivative of the variable withan accuracy of order 2m, and am the difference weight coefficient(6) This equation shows that the first derivative of the staggered-gridhas an accuracy of the order 2m For coefficients under a different(7) order of accuracy one may refer to the existing literature[12].Fromthe truncation error of the calculations it can be shown that the er(8) ror decreases gradually as the difference order becomes higher. Forexample, the error has reached 0.000218 with order 10 accuracywhere v or vz is the solid-phase velocity and the subscript repre-Following these steps the difference scheme of an isotropic por-sents either the x or z direction. The fluid-phase velocity is given ous elastic wave may easily be derived Taking Eqs. (1)and(5)forby ix or ir in the x or z directions. The solid-phase stresses are example. with a second order time accuracy and a tenth orderoxx, z and Txx. The fluid-phase stress component is given by s. space accuracy we findThe dynamic, mass density of the medium is given by P12. pu. v+l-w(Piz-PuPa(2124"yand p22 and bu is the dissipation factor of the medium. The solidphase elasticity coefficients are Gu. C13. C31. C33. and Css. whichare usually represented in terms of A and N for an isotropic, porous+(n2+2)(a12时-bhn2(1medium as, e.g.C1=A+2N,C13=C3=A,C3=C55=Ndt*{Cn1∞4+C13lFinally, Q is the solid-fluid coupling elastic parameter and R ishe fluid-phase elastic parameter.+Q1(tb+l°+)(11)The staggered-grid finite difference of these equations is foundin two steps. First, the kinematic variables such as speed and the where the subscript of each variable is the spatial location, thedynamic variables such as stress are located on a grid. The projec superscript is the time location, L and ly the x and the z deriva-tion of each variable is found by an algorithm Taking Eq (1). for tives of the variable to a tenth order accuracy. The differenceexample, the first derivative of the stress variables s and o on scheme of the other equations may be deduced by analogythe right hand side are located with respect to x. the stress variableox is located with respect to z coincident with the speed variable 3. Key issues in the analyses of the numerical simulationv on the left hand side located with respect to t. the rest can bededuced by analogy3. 1. Boundary conditionsFollowing this rule the arrangement of the variables is obtainedas shown in Fig. 1. The left part is the time-space staggered-grid Eliminating the boundary reflection is one important part ofwhere variables are represented by the same color when they seismic numerical modeling. The PML absorbing boundary condi-occur at the same time. The locations in directions x, y, and t are tion is one maturerepresented in terms of the variables i, j and n the rightmost part boundary condition中国煤化工 blished by first.is the projection of the different spatial variables at a certain time. based on the coneIf u and i, are set at location 1, and then represented by(i.j. n). the stress versus veCNMHG decomposingand introducingG. Zou et aL. /Mining Sciencelogy( China)21(2011)727-73用4100ource posoonm4=04≠0200,80)d≠0,d0N3.509 GPa, /a35711 GPa,Interior0=1.0513 GPa, R-0 8105 GPa,Pu=1688.5 kg/m, Pu"194.0 kg/,,00200300400bur 7. 605 MParsmFig 4 An X component snapshot in the solid- and fluid-phases(leftFig. 2. Perfect match layer diagram and corresponding physical parametersIn porous media the seismic source is not only loaded by a solidbut also by a fluid: Deciding the ratio of solid to fluid is a problemPorous medium theory suggests the liquid phase is homogeneouslydistributed so the ratio of fluid loading on the seismic sourceshould be decided by the porosity. If B is the porosity and ws W号0.5are the weighting factors for solid and fluid motion, respectively.AA八A~we may writeWs=1-B,W=p(12)The seismic source can act on both the vertical stress compd10 PMLnent and the fluid stress component. then it seems as a true"stress-1.030 PMIsource."The seismic source can also act on the velocity componentof the solid particles and a fluid element whereupon it appears as a"velocity source". In this paper the Ricker wavelet is regarded as a0.102"stress source". The dominant frequency of the seismic source is aparameter related to the model grid and the dispersion. For fixedT(s)medium parameters an increase in the dominant frequency ofFig 3. x-component waveforms for differentthe source requires that the spatial and temporal steps must be re-duced in sizea damping factor in the X and Z directions. Second, these equations 3.3. Stability analysisare then transformed into the frequency-domain thereby changingthe Cartesian coordinates into the complex plane. Finally, a reverseHigher order differences are used in the finite-difference meth-Fourier transformation gives the PMl absorbing boundary condiboundary condi- od. These calculated results show that the truncation error is verytion required for the porous media in the time domain. A practical small. However, when the equations involve successive differencesexample using this is given here.cumulative effects from errors become obvious. Taking the timeAs shown in Fig. 2, a single-layer isotropic porous medium mod- step, for example, the n- and(n +0.5)step results are used to cal-el may be designed to test the effect of the PMl absorbing bound- culate the(n+ 1yth result. The errors in the n-and the(n+0.5)ary conditions. the blank space around the source is the internalsteps affect the(n+1)step result The error adds up step by stepregion, the edges are the PMl, and the physical parameters are de- An analysis of error propagation is required: this is referred to asnoted in the right, lower cormer of Fig. 2. The boundary is dividedstability analysis.perpendicular to the x-axis, and at each of the four corners. the h the stability analysis in a porous medium is similar to that forinto three parts: In a direction parallel to the x-axis, in a directionelastic medium. The stability condition can be expressed asattenuation factors dx and dy of each part are denoted in the upper [11, 131right hand corner of Fig. 2. The model grid here is an 81*81 onewith a grid size of 5 m. The p-wave stress source, with a frequencyy(v}2+≤△h(13)of 20Hz, is located at position(200, 80)while the receiving point is where Ah is the grid step in the x- and the z-directions, v theat(0: 400, 280). The time sampling is done every 0.0005 ms and thevelocity of the fast P-wave in the saturated porous medium, andtime accuracy is of order 2 while the space accuracy is of order 10. vs the velocity of the S-wave. The stability condition of the regu-The number of perfect match layers is set to 10, 20, or 30. Theresulting X-component of the seismic wave in the solid-phase is lar interior region is not critical but the stability condition at theshown in Fig. 3. It can be seen that the boundary reflection is elim- PML absorbing boundary region needs a more rigorous standardinated to good effect in the case of 20 or 30 layers. Fig 4 shows the When the parameter b is negligible, in practical numerical mod-wave field snapshots of the seismic waves at 400 ms in the case of eling, the PML boundaries must meet the following conditions20 layers. From this it can be seen that the strong boundary reflec- [13,141tion is effectively eliminated△V2v≤min(△x,△z)(14)3.2. Source implementation3. 4. Wave dispersionSource implementation is another key factor when modeling中国煤化工wave propagation In Table 1 some typical wavelet functions usedDifferent frequenYHropagate at dif-ferent velocities AsCN MH Spread out: this730Zou et aL/ Mining Science and Technology( China )21(2011)727-731Table 1Typical wavelet functions.Wavelet type/avelet function and the main parametersf(r)=(1-2m26r2)elr, fo is the center frequency of seismic wavesZero-phase band-pass wavelet f(r)=2 sin(2rAft)cos(2rfor)/rt, fo is the center frequency of seismic waves, Af is half the difference between high and low cutoffCosine waveletf(r)=cos(2rfore-ar, a is the attenuation coefficient. fo is the center frequency ofLimited bandwidth pulse point f(o=e-ait- tor, a is the attenuation coefficient, to is the center frequency of seismic wavesphenomenon is called seismic wave dispersion. During the numer-1000ical modeling seismic event splitting, waveform distortion, orThe plysical parameter in the first layerwaveform ghosting may all be seen. All of these are caused by 400№7.125GPa.A=5.2GPaseismic wave dispersion [15]0-25575 GPa R=3.62 GPa.There are two reasons for this. One is the properties of the med- 1000L-P=188kgm3,p2=0.131kgm3,um. For example, a wave propagating in the reservoir sand showsPz=0.349 kg/m, bu=0.52 uPas/m2wave dispersion at ultrasonic frequencies. This phenomenon wasThe physical parameter in the first bayeproved by numerical simulation. the other is that unreasonable=6535 GPa. +11. 37 GPa.model parameters can cause seismic wave dispersion during20.895 GPa, R=0.321 GPa,numerical modeling. For example factors that affect seismic wavep1=207kg3,p0.085kgh3dispersion include the time step, the grid size, the propagation dis- 2000Pz=0. 181 kg/m, b1 3.2 HPas/tance, and the speed [16]. The latter cause of dispersion should beavoided in the finite-difference method. there are two methods to Fig. 6. Perfect match layer diagram and physical parameters of the example model.weaken this seismic wave dispersion. One is raising the differenceaccuracy in the spatial and also the temporal regime the other isdecreasing both the grid and the time stepsTaking Model 1 for example, temporal second-order accuracyand spacial tenth-order accuracy are used and the dominant fre-quency of the seismic source is 20 Hz. The spatial and temporalsteps are either Ax=Az=10mand△t=0001s,or△x=△z=5mand At=0.0005 s. The resulting waveforms of the X componentin the solid are shown in Fig. 5. A grid step of 10 m gives an unreasonable waveform after 0. 28 s. However, the waveform is verygood when the grid step is 5 m. Generally speaking a second-order5001000150020000500100015002000accuracy in the temporal domain, a tenth-order accuracy in thespatial domain, and a wavelength over more than ten grids pointswill give acceptable numerical resultsFig. 7. A Z component snapshot in the solid or fluid phases (left or rightA large grid step in the lateral components and a small grid stepespectively) at 400 ms, PML boundary conditions.in the vertical components can control dispersion when the medium is a low speed body [17]. The velocity in the coal bed is usually1800 m/s, but the velocities at the seam top and floor are larger Forexample, a typical top and floor velocity would be more than3500 m/s. This large velocity variation makes the coal bed thelow speed body and dispersion easily occurs in this situation, EIncreasing the grid density in the lateral direction improves results i 60in this regard; however the stability may be affected by the grid04008001200160020000400800120016002000x(m)0.5Fig& A single-shot record of the X (left)and Z (right) components in the solidand the time step size. For this reason changing these parametersshould be comprehensively investigatedΔx=△z=10m,△r=0.001s-△x=A=5m,△=0.005s4. Results and discussion0.2030.4The geological mIsotropIc porous m中国煤化工 interface depthis 1000 m. The comFig. 5. Waveform of the X component in the solid-phase10 m grid step andHCNMHGOnodes with a-dary. A P-waveG Zou et al/ Mining Science and Technology(china)21(2011)727-731source with a dominant frequency of 20 Hz is located at spacial po-control dispersion when the medium is a low speed bodysition(400, 1000). Receivers are located at(0: 2000, 600)and theThis is a useful way to simulate wave propagation in asample interval in time is 1 ms Second order temporal accuracycoal-bearing porous medium.and tenth order spatial accuracy are used. The P-wave velocity assigned to layer 2 is 2000 m/s, which is typical for porous coal.The model parameters are summarized in Fig. 6.AcknowledgmentsFig. 7 shows snapshots of the solid particle, X component at atime of 400 ms the fast P-wave is marked Pi and the slow p-waveThanks for the opinions of the judges and the editors. This paperis marked P2. The s-wave is just designated as such. this figure is supported by the National Basic Research Program of China( Nos.shows the P1 wave the P1-S wave the p2-P1 wave and the p2-s 2009CB219603 and 2006CB202209) the National Natural Sciencewave from the interface. Transmitted waves are also apparent from Foundation of Special Equipment(no 50727401 )and the Nationalthe second layer, namely: the Pl wave, the Pl-s wave the Science Technology Pillar Program in the Eleventh Five-Year PlanP2-Plwave, and the P2-s wave. Fig. 8 shows the shot gather of Period (no. 2007BAK28B03)the solid particle X- and Z-components. Here, the direct P1 wave,the direct P2 wave, and the reflected P1, P2, P1-S, and p2-s waves Referencesare seen[1] Ma ZT. Introduction to Geophysics Computation. Shanghai: Tongji University5 Conclusions[2] Wang XC, Xia CL Liu xw. Two-dimensional P-wave field continuation ofundulating surface. Min Sci Technol 2009: 38(5): 112-6A finite-difference simulation method has been studied for mod- 13] Guo J P-wave field in Biot medium simulated by finite difference method, oileophys Prospect 1992: 12(2): 515-20ling seismic waves propagating in a coal bearing porous medium[4]Dai N, Vafidis A, Kanasewich ER. Wave propagation in heterogeneous, porousThe simulation is based on finite-difference theory and PMl absorb-velocity-stress, finite- dce method, Geophysicsing boundary conditions. The key factors considered are grid space1995:602):327-40step, variable location, order of the calculations, and the PML absorb [51 Carcione JM. Quiroga-Goode G. Some aspects of the physics arodeling of Biot P-waves. J Compute Acoust 1996: 3: 262-8ing boundary conditions. We draw the following conclusions:[6] Ozdenvar T, McMechan G. Algorithms for staggered-grid computations forporoelastic, elastic. acoustic and scalar wave equation. Geophys Prosp1)The variables are divided into dynamic and kinematic vari-199:45403-20[7] Zhang HX, He BS, Ning SN. High-order finite diffeables. Their positions must be mutually staggered. The pro-4:4:307-13jection of variables has a position function: The location of 18] Liu KA. zhang xM, Liu jQ The wavelet multithe first derivative of the dynamic variable with respect to 19) Pei ZL. a staggered-grid high-order finitemethod for modelinrespect to time. the difference accuracy increases for higherNat Sci)]uation in 3-D dual-phase anisotropic media. J China Univ Pet16-20.order differencing.[10] Collin F. Tsogka C. Application of the perfectly matched absorbing layermodelcan effectively attenuate outgoing seismic waves. the [111 Zeng YQ He jQ, Liu QH. The application of the perfectly matched layer innumerical results show that outgoing waves are effectivelyal modeling of wave propagation in poroelastic media, Geophysics200166(4):1258-66.absorbed and that the reflection is very small when more (12 Mou YG. Pei ZL Seismic numerical modelingthan 20 cells of Pml are usededia. Beijing: Oil and Industry Press: 2007.(3)Biot theory considers the liquid-phase as a homogeneous 131 Zhao HB. Wang xM, Wang Di chen H Applications of the boundardistribution so the solid to liquid ratio at the seismic sourcedepends upon the porosity.reservoirs using Biot theory. Geophysics 1991: 56(3):328-39 poroelastic[14] Zhu X, McMechan GA. Numerical simulation of seism(4)Wave dispersion is controlled by the grid and the time stepize smaller sized steps can weaken wave dispersion. Gen[15]Deng Jx. Yu J. The research of the characteristics of ultrasonic velocityerally speaking, second-order accuracy in the temporal and2003:39:835-4tenth-order accuracy in the spatial domains, and a source [161 Sun cY, Gong T]. Zhang Y Zhang w. Analysis on dispersion and alias inwavelength covering more than ten grid points will providefinite-difference solution of wave equatIon. Oil Geophys Prospect2009:44:43-8acceptable numerical results. a large grid step in the lateral [17 Li S). Sun CY. Gao JH. Analysis of dispersion suppression in wave equationdirection and a small grid step in the vertical direction willnumerical simulation. Geophys Prospect Pet 2008: 47(5): 444-8H中国煤化工CNMHG