Perfectly matched layers simulation in 2-D VTI media
PPerfectly matched layers simulation in 2-D VTI media
P. Contreras, G. Larrazabal, and C. Florio
Multidisciplinary Center for Scientific Visualization and Computation(CEMVICC) Science Faculty Universidad de Carabobo, Carabobo, Venezuela. (Dated: November 12, 2019)Implementing computational boundary conditions, such as perfectly matched layers
PML doeshave advantages for forwarding modeling of the earth’s crust. The mathematical modeling of manyphysical problems encountered in industrial applications often leads to a system of linear partialdifferential equations
PDEs . It considerably improves the visualization of seismic events relevantto oil and gas exploration. In this work, we present an efficient numerical scheme for hyperbolicpartial differential equations where a computational technique takes care of reflections at the bordersdomain using a linear elastic-wave system of decoupled equations with
PML -type boundaryconditions. The key idea is to introduce a layer that absorbs the reflections from the bordersimproving images visualization. Anisotropy has been reported to occur in the Earth’s three mainlayers. The crust, the mantle, and the core; but this implementation refers to the case of a verticaltransversal anisotropic medium
VTI in the crust-layer. Images and screen-shots of the longitudinal P impulse-response and the SV transverse impulse-response are obtained at different times. Thiscomputational method enables us to achieve images for the P and SV response-impulses, and toobtain high quality synthetic seismograms for the PP and PS reflection events in a 2-D VTI two-layer model.
Keywords:
Seismic anisotropy,
VTI medium, computational modeling, perfectly matched layer.
PACS numbers: 07.05.Tp 94.20.Bb 91.30.Cd 91.30.Ab
Introduction
The wave equation plays an important role in seismic oil and gas exploration since it allows modeling the earthstructure. In the past decades, seismic anisotropy has been gaining attention from academic and industry, in partthanks to advances in anisotropy parameter estimation
A material is said to be anisotropic if the value of one or more of its properties varies with direction, consequentlyanisotropy consideration improves the subsurface imaging. Seismic Anisotropy can be defined as the dependence ofseismic velocity on a direction or upon an angle. Anisotropy is described by a 4 th order elasticity tensor with 21independent components for the lowest-symmetry case . However in practice, observational studies are unable todistinguish all 21 elements and anisotropy considerations are usually simplified. For seismic exploration, the mostcomplicated case occurs in fractured monoclinic media, with 9 elastic constants . In general, two or more sets ofvertical non-corrugated and not perpendicular fractures produce an effective monoclinic medium with a horizontalsymmetry plane. The second most important application occurs for the orthorhombic model. The orthorhombicmodel describes a layered medium fracture in two orthogonal directions , However in the simplest form that seismicexploration uses it, there are two main kinds of transverse anisotropy TI . One is called horizontal transverse isotropy- HTI which is a common model in shear-wave studies of fractured reservoirs that describes a system of parallel verticalpenny-shaped cracks embedded in an isotropic host rock , henceforth this kind of anisotropy is associated with cracksand fractures. The other one is called vertical transversal anisotropy- VTI and it is associated with layering andshales. Sometimes people call it vertical polar anisotropy.In the beginning, forward modeling was done by simulating the scalar P wave field obtained from the acoustic waveequation. However, the earth-crust is elastic and anisotropic and all propagation modes should be considered inorder to observe anisotropy effects. Forward modeling and parameter estimation are almost the most fundamental toall other anisotropy applications in oil exploration . In part, this can be done by numerically solving the hyperbolicelastic wave equation WE . The vector and tensor fields which represent the WE improves the information about the3D earth-crust geology. This translates into a better subsurface imaging and it is of extreme importance for the oiland gas seismic exploration.The most common type of anisotropy that occurs in the earth’s crust is the vertical transverse anisotropy VTI whichis observed in sedimentary rocks . To achieve 3D seismic modeling with VTI anisotropy requires the knowledge offive elastic constants. However in two dimensions only four constants are involved, namely C , C , C and C (in this work we use Voigt notation for the elastic constants) . Henceforth shear wave splitting is not consideredin a 2-D modeling because of the lack of the elastic constant C . Abundant geological evidence of shales shows theimportance of VTI models in seismic reservoir characterization .In this work, we implement an algorithm using a linear finite difference decoupled
PDEs in terms of the velocities a r X i v : . [ phy s i c s . g e o - ph ] J a n and stress tensor components. To implement such a model we use central-difference second order approximationfor space variables and time partial derivatives. Henceforth we use the finite-difference time-domain FDTD method . In addition, we use staggered stencils for computational storage . Staggered cells grant the calculation ofdifferent physical quantities at different mesh-grid points (see Figure 2), reducing the computational cost storage. Ad-ditionally, with this algorithm we can introduce different source-receiver configurations, making it ideal for simulationin oil and gas seismic modeling.To implement a FDTD solution of the WE , a computational frame must first be established. The computationaldomain is the physical region over which the simulation is performed. To achieve the condition of non-reflectiveborders, we use the artificial technique of perfectly coupled layers PML . Consequently, the effect due to thereflection of the waves at the borders of the domain is automatically reduced.
PML implementation in seismic exploration requires a reformulation of the linear WE to eliminate the unwantedreflections. Thus, this work is structured in the following way: Section I introduces the topic of seismic anisotropyand its relevant application in R&D Oil industry. Section II briefly describes the PDEs for a 2-D medium with
VTI shale anisotropy, and the theoretical implementation of the
PML obtaining a decoupled
EWs (elastic equationsystem) including the new border computational conditions. Section
III describes the numerical implementation thatwill allow CPU time reduction. Subsequently, we find the correct decoupled finite-difference scheme in terms of thestaggered-stencil mesh taking into consideration the
PML (see Figure 2). Finally, in section IV the outcome showsthe applicability of this technique. We were able to achieve very neat and sharp subsurface images without reflectionevents at the edges of the 2-D earth-crust two dip layer shale model. As we stated before, the formulation of a linear
WEs in terms of temporal derivatives of the velocity and the stressfields can be used to propagate waves. The
FDTD technique reproduces the fields forward in time-domain. This iscalled forward modeling. This procedure is designed through a 3D staggered finite-difference grid . This formulationis widely used in the literature for seismology purposes . In 2-D VTI anisotropy, it involves the transformation offive coupled first-order
PDEs , namely the two general equations of motion for the vector velocity ( v x , v z ) field whichare: ρ ∂v x ∂t = f x + ∂σ xx ∂x + ∂σ xz ∂z , ρ ∂v z ∂t = f z + ∂σ xz ∂x + ∂σ zz ∂z , where ρ is a constant density, f x and f z are external forces driven by the source. Since we use the FDTD time-domaintechnique, a seismic pulse is used as the source, then the response of the system over a wide range of frequencies canbe obtained with a single simulation.And the three equations for the stress σ ij second order tensor field σ ij = (cid:20) σ xx σ xz σ zx σ zz (cid:21) where only four elastic constants are used as the input variables for the WEs : ∂σ xx ∂t = c ∂v x ∂x + c ∂v z ∂z , ∂σ zz ∂t = c ∂v z ∂z + c ∂v x ∂x , ∂σ xz ∂t = c (cid:104) ∂v x ∂z + ∂v z ∂x (cid:105) , the longitudinal P and the transverse SV modes given by the previous five PDEs are coupled . In this subsection, the
PML technique is briefly explained and applied to the 2-D modeling. In order to make thefinite difference method more stable and convergent, we follow the arguments presented previously for the stabilityand the dispersion conditions in a 2-D
EWs (we should mention that a 3D anisotropy implementation has not beenachieved with this technique, due to failures in the stability and dispersion conditions when considering anisotropy ).Even in works dedicated to seismology, the use of this technique helps to improve 2-D imaging resolution at largerscales. The seismological domains have borders as it happens in seismic exploration. PML in this case, reducesthe reflections at the borders in a significant amount compared with the use of traditional mathematical boundaryconditions .The
PML mathematical implementation consist in the introduction of a modified coordinate system, where theexpansion coefficient is a complex number with an evanescent imaginary part. This generalization is achieved throughthe following substitution ∂∂x →
11 + i p ( x ) ω ∂∂x , and ∂∂z →
11 + i p ( z ) ω ∂∂z , where p ( x ) and p ( z ) are the coefficients of the PML . They are given by the following expressions p ( x ) = p ( x/L ) N , and p ( z ) = p ( z/L ) N , with p = − v p logR c / (2 L ), where L is the thickness of the perfectly coupled layer, N the size of the problem and p has an approximate value of 341.9When the above coordinate system is replaced in the equations for the VTI medium of the previous section, thesystem becomes a new linear decoupled
PDE - EWs . Because of that replacement appears new equations for v x and v x new decoupled velocities: ρ (cid:104) ∂∂t + p ( x ) (cid:105) v x = f x + ∂σ xx ∂x , ρ (cid:104) ∂∂t + p ( z ) (cid:105) v x = ∂σ xz ∂z , For the decoupled v z and v z velocities the equations are: ρ (cid:104) ∂∂t + p ( x ) (cid:105) v z = f z + ∂σ xz ∂x , ρ (cid:104) ∂∂t + p ( z ) (cid:105) v z = ∂σ zz ∂z . These decoupled velocities are related to the coupled ones by the relation v x = v x + v x and v z = v z + v z .The stress tensor field is mathematically treated in the same way, obtaining a decoupled system for the new stressequations: The new decoupled σ xx and σ xx stress equations are: (cid:104) ∂∂t + p ( x ) (cid:105) σ xx = c ∂v x ∂x , (cid:104) ∂∂t + p ( z ) (cid:105) σ xx = c ∂v z ∂z , and the new equations for the σ zz and σ zz decoupled stress tensors components equal to: (cid:104) ∂∂t + p ( x ) (cid:105) σ zz = c ∂v x ∂x , (cid:104) ∂∂t + p ( z ) (cid:105) σ zz = c ∂v z ∂z . The decoupled shear σ xz and σ xz stress tensor components follow the equations: (cid:104) ∂∂t + p ( x ) (cid:105) σ xz = c ∂v z ∂x , (cid:104) ∂∂t + p ( z ) (cid:105) σ xz = c ∂v x ∂z . where now σ xx = σ xx + σ xx , σ zz = σ zz + σ zz , and σ xz = σ xz + σ xz . These 5 linear PDE equations will be used for aforward modeling in the following section.
Computational spacetime implementation
The computational model that we propose in this simulation is divided into a mesh of N x × N z points. Fig. (1) thefinite-difference scheme shows the computational edge domain where the PML are applied according to the two-layerdip model. ∆ x and ∆ z are defined as the distances between the points in such a way that x = n x ∆ x y z = n z ∆ z with n x = 1 ...N x y n z = 1 ...N z . For the step in time ∆ t we have that t = n ∆ t , being n the step time (in Fig. 2) thestencil time is not showed .The five variables v x , v z , σ xx , σ zz , and σ xz are discretize into a two dimensional staggered grid mesh, see Fig. (2) fora better explanation, where the velocity components are stored in both stencils, the normal components of the stresstensors are stored in one of the stencils, while the shear stress components are located at another stencil. However,the major difference of a staggered-grid scheme is that the velocity and stress components are not known at the samegrid point, as it can be seen from Fig. (2), henceforth we use a previous scheme proposal with little modifications where different stencils are used for normal and shear stress field computations. FIG. 1: Finite-difference simplified scheme showing the computational 2 two-layer dip model. The vertical, horizontal andcrossed lines near the borders show where the
PLM are appliedFIG. 2: Staggered finite-difference grid for the velocity and stress updates, both stencils show the grids points where differentfields components are calculated. Velocity components are calculated in both stencils (black circles), the shear stress componentsare calculated in the blue stencil (red starts), and the cinnamon stencil is used to calculate the normal stress components (greensquares).
Following Fig. (2) in the space domain v x is calculated at the points ( i ± / , k ), v z is calculated at the points( i, k ± / σ xx and σ zz are calculated at the points ( i, k ) and finally the shearstress σ xz is calculated at the points ( i ± / , k ± / ρ ( x, y, z ) is considered a constant density ρ (seethe last line of Table I for its correspondent values and units) Physically it means that the VTI model is homogeneousand that is does not take into consideration earth-crust heterogeneities.An explosive source is simulated using the amplitude of a Ricker wavelet with a peak frequency f ∼
30 Hertz.Amplitudes are added to the velocities v x and v z at each n time-step. The stability limit of the discrete system isgiven by the following condition:∆ t < .
606 ∆ xV p , where V p account s for the horizontal longitudinal speed . The dispersion condition for the discrete system is givenby the inequality: ∆ x < V s f , where V s is the transverse wave speed. However some authors showed that the stability of the classical PML modeldepends upon the physical properties of the anisotropic medium and that it can be intrinsically unstable. Howeverthe elastic constants in table where previously proven to meet Becache stability criteria. This algorithm could beextended to 3-D situations by further studying the stability and dispersion conditions.Henceforward we proceed to obtain the discretization of the velocity fields as follows for the decoupled finite-difference system: the v x and the v x components are computed at time-steps points t n +1 and t n (the stencil-time isnot shown here ) with the following finite-difference equations: v x | n +1 i +1 / ,k = v x | ni +1 / ,k e − p ( x )∆ t + 1 ρ ∆ t ∆ x e − . p ( x ) ∆ t (cid:104) σ xx (cid:12)(cid:12)(cid:12) n +1 / i +1 ,k − σ xx (cid:12)(cid:12)(cid:12) n +1 / i,k (cid:105) , (1) v x | n +1 i +1 / ,k = v x | ni +1 / ,k e − p ( z )∆ t + 1 ρ ∆ t ∆ z e − . p ( z ) ∆ t (cid:104) σ xz (cid:12)(cid:12)(cid:12) n +1 / i +1 / ,k +1 / − σ xz (cid:12)(cid:12)(cid:12) n +1 / i +1 / ,k − / (cid:105) . (2)The v z and v z components are calculated at the times-step t n +1 and t n points with the help of equations: v z | n +1 i,k +1 / = v z | ni,k +1 / e − p ( x )∆ t + 1 ρ ∆ t ∆ x e − . p ( x ) ∆ t (cid:104) σ xz (cid:12)(cid:12)(cid:12) n +1 / i +1 / ,k +1 / − σ xz (cid:12)(cid:12)(cid:12) n +1 / i − / ,k +1 / (cid:105) , (3) v z | n +1 i,k +1 / = v x | ni,k +1 / e − p ( z )∆ t + 1 ρ ∆ t ∆ z e − . p ( z ) ∆ t (cid:104) σ xx (cid:12)(cid:12)(cid:12) n +1 / i,k +1 − σ xx (cid:12)(cid:12)(cid:12) n +1 / i,k (cid:105) . (4)For the stress fields the discretization is as follows: First, the decoupled equations for the normal stress tensorcomponents σ xx and σ xx are computed at the mesh step-time points t n ± / using the equations: σ xx | n +1 / i,k = σ xx | n − / i,k e − p ( x )∆ t + C ∆ t ∆ x e − . p ( x ) ∆ t (cid:104) v x (cid:12)(cid:12)(cid:12) ni +1 / ,k − v x (cid:12)(cid:12)(cid:12) ni − / ,k (cid:105) , (5)and σ xx | n +1 / i,k = σ xx | n − / i,k e − p ( z )∆ t + C ∆ t ∆ z e − . p ( z ) ∆ t (cid:104) v z (cid:12)(cid:12)(cid:12) ni,k +1 / − v z (cid:12)(cid:12)(cid:12) ni,k − / (cid:105) . (6)Second, the decoupled equation for the normal stress components σ zz and σ zz are computed at the time points t n ± / by means of the following expressions: σ zz | n +1 / i,k = σ zz | n − / i,k e − p ( x )∆ t + C ∆ t ∆ x e − . p ( x ) ∆ t (cid:104) v x (cid:12)(cid:12)(cid:12) ni +1 / ,k − v x (cid:12)(cid:12)(cid:12) ni − / ,k (cid:105) , (7)and σ zz | n +1 / i,k = σ zz | n − / i,k e − p ( z )∆ t + C ∆ t ∆ z e − . p ( z ) ∆ t (cid:104) v z (cid:12)(cid:12)(cid:12) ni,k +1 / − v z (cid:12)(cid:12)(cid:12) ni,k − / (cid:105) . (8)Third, the decoupled equations for the shear σ xz and σ xz components are computed at time points t n ± / by meansof: σ xz | n +1 / i +1 / ,k +1 / = σ zz | n − / i +1 / ,k +1 / e − p ( x )∆ t + C ∆ t ∆ x e − . p ( x ) ∆ t (cid:104) v z (cid:12)(cid:12)(cid:12) ni +1 / ,k − v z (cid:12)(cid:12)(cid:12) ni − / ,k (cid:105) , (9)and σ xz | n +1 / i +1 / ,k +1 / = σ zz | n − / i +1 / ,k +1 / e − p ( z )∆ t + C ∆ t ∆ z e − . p ( z ) ∆ t (cid:104) v x (cid:12)(cid:12)(cid:12) ni,k +1 / − v x (cid:12)(cid:12)(cid:12) ni,k +1 / (cid:105) . (10)The elastic constants used for the two-layer model are given according to table . These values correspond to differenttypes of VTI media and they are able to accomplish the stability condition using a
PML computer domain . TABLE I:
VTI elastic constants for model Fig. elastic constants ( × N m − ) upper layer bottom layer C C C C ρ ( Kg m − ) 7.100 3.200 Analysis and Results
The present paper describes a new space-time 2-D methodology with time-step and space-step control to solve 2-D
VTI linear decoupled
EWs with
PML conditions using staggered-grids Fig. 1 and Fig. 2. Henceforth, we developa computational model for an elastic wave propagation simulation for 2-D
VTI anisotropic media, by combining a
FDTD method on staggered grids with a
PML boundary condition. For that, a High-Performance Linux applicationin C (over 1,000 lines of code) using Make, GDB and Valgrinds memcheck, to generate and visualize 2-D response-impulses and synthetic seismograms was developed.We establish that visualization of the impulse-response of the P and SV modes are improved by the PML . As aconsequence, unwanted reflections from the borders in Fig. 3 are totally eliminated according to Virieux .In Fig. 3 simulations were performed with an explosive source (Ricket-wavelength) centered in the middle of themodel to eliminate reflections from the edges. In Fig. 3 time screen-shots for t = 0 . , .
3, and 0 . PML boundary conditions absorb the reflections from the border,implying an exceptional numerical performance of the computational
PML technique. Another conclusion is that the P SV response-impulse presentstriplications. Moreover, the effects of the VTI elastic constants on SV response-impulses affects the direction ofpropagation. We observe from the snapshots in Fig. 3 how the value of the elastic constant C affects the directionof the SV . Triplications in the SV wavefront are modeled. However, one question is not solve yet in this work, theasymmetric behavior of the SV impulse-response in Fig. 3 with respect to the vertical axis.Subsequently, the explosive source is relocated to the upper central part of the 2-D model sketched in Fig. 4 A .On the surface, fifty receivers are uniformly distributed on both sides of the source. We see how the velocity fieldscancel out at the edges of the model in Fig. 4 A . Then we obtain the synthetic seismograms for the two layers model,each of them has different elastic constants values. The synthetic seismograms of the vertical component V z (Fig. 3 B ) and the horizontal component V x (Fig. 3 C ) are shown with the primary reflections events due to the dip layer.Henceforth, it remarkably shows the effectiveness of the PML using the decoupled linear system of equations(1)-(10) for 2-D
VTI media simulation in Oil and Gas R&D.
Acknowledgments
One of the authors, P. Contreras wishes to express his gratitude to Dr. V. Grechka for pointing out an invariantsymmetric question for the SV impulse-response. The authors also acknowledge Drs. E. Sanchez, and Prof. J.Moreno for several discussions regarding the numerical implementation of this work. Perez M., Grechka V., and Michelena R., Geophysics. (1999) 1266. Grechka V., and Tsvankin I., Geophysics. (1998) 1079. Musgrave, M. Crystal Acoustics. Holden-Day. (1970). Helbig, K. Foundations of anisotropic for exploration Geophysics. Pergamon Press. (1994). L. D. Landau and E. M. Lifshitz. Elasticity Theory. Pergamon Press, 1982. Dellinger, J. Anisotropic seismic wave propagation. Ph.D. thesis. Stanford University, (1991). Grechka V., Contreras P. and Tsvankin I., Geophysical Prospecting. (2000) 577. Contreras P., Grechka V. and Tsvankin I., Geophysics. (1999) 1219. Tsvankin I., 66th Ann. Internat. Mtg., Soc. Expl. Geophys. (1996b) 1850. Contreras P., Klie H., and Michelena R., 68th Ann. Internat. Mtg. Soc. Expl. Geophys. (1998) 1491. L. Thomsen, Geophysics. (1986) 1954. D. F. Winterstein, Geophysics. , (1990) 1070. Faria E. y Stoffa P., Geophysics. (1994) 282. Virieux J., Geophysics. (1986) 889. Festa G. and Nielsen S. Bulletin of the Seismological Society of America. (2003) 891. Komatitsch D. and Martin R., Geophysics. (2007) SM155. Kane Y., IEEE Transactions on Antennas and Propagation. (1966) 302. Becache E., Fauqueux S., and Joly P., J. Comput. Phys. (2003) 399.
P SVP SVSVP
FIG. 3: Screen-shots at 0.2, 0.3 and 0.5 sec. are displayed for the vertical velocity wave-field component propagating througha homogeneous medium with elastic constants from the upper layer given in Table I. The triplication of the SV wave aroundthe intermediate angles are observe to be asymmetrical for the SV impulse-response. The absorbing effect of the PML at theedges of the computational domain is notable.
C: VxB: VzA: MODELOFUENTE PP PPPSV PSV
FIG. 4: Two-layers