- 237 Views
- Uploaded on
- Presentation posted in: General

Applied Seismic Simulation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Applied Seismic Simulation

by Lisbeth Engell-Sørensen Centre of Integrated Petroleum Research (CIPR)URL: http://www.ii.uib.no/~lisbethe-mail: Lisbeth.Engell-Sorensen@cipr.uib.no

From Reservoir to Geophysics:

Seminar on Integrated Modelling

CIPR, May 9, 2003

The main purpose of this talk is to resume the seismic methods and to help all CIPR employees and students to choose and understand the available codes for their purpose.

Some definitions used in geophysics:

- Waveform modelling (inversion/trial and error/forward)
- Waveform simulation (forward)
- Asymptotic waveform modelling (phases)
- Full waveform modelling (seismograms)

- Finite-difference simulation of equation of motion (10 slides)
- Full-waveform modelling with surface waves (3 slides)
- Asymptotic wave propagation with ray-tracing (10 slides)
- 1D simulation using reflection coefficients (3 slides)
- Other methods for surface waves and body waves: reflectivity method, spectral and pseudo spectral methods, finite element and spectral element methods (1 slide)

nc=non commercial=academic use

Wave propagation in a general anisotropic, inhomogeneous elastic solid initiated by a body force is given by the equation of motion (Ben-Menahem et al, 1991):

Einstein’s summation convention for repeated indices is used. fj is body force density, ui are displacement components, ij is stress tensor, Cijkl is the tensor of elastic moduli, kare spatial partial derivatives, t is temporal derivative, Ti = ijnjis the traction on a surface with normal n, Cijkl = Cjikl = Cijlk = Cklij (21 independent Cijkl coefficients).

The displacement, u, due to body forces, f, throughout V and to boundary conditions, u and T, on S is (Aki and Richards, 1980)

where G is Green’s tensor. The first term is the contribution from the body forces in V and the surface integral gives the contributions from the boundary conditions on the boundary S of V.

The general anisotropic code uses 21 elastic coefficients.

Eight independent coefficients are needed as input parameters for the general TI medium: density, P-, S-velocity, three Thomson parameters: , , , and two angles for the TI symmetry axis.

(for the following FD modelling described here see: Hokstad et al., 2002; Holberg, 1987; Mittet et al., 1994; Petersen, 1999; Engell-Sørensen and Koster, 2002, Engell-Sørensen, 2002)

Newton’s law:

Hooke’s law:

The parallel code uses multiple processors in order to manipulate subsets of large amounts of data simultaneously

1. Model large problems in limited time

2. Larger memory problems on more processors

Message Passing paradigm MPI (’nested’ parallelism)

- To compute a higher-order difference or shift operation several points are needed on each side of every grid point: Partition with overlapping sub-domains
- We define a common boundary with eight points for eight’s-order FD operator (half-length = eight)
- Since every grid point needs contributions from eight neighbouring grid points in forward and backward directions, the local data in the boundary that is shared by neighbouring sub-domains must be exchanged and added to the neighbouring sub-domain data
- 3D case: at most 26 neighbouring sub-domains with which information has to be exchanged

The main purpose of the modelling is to study waveform propagation and generate realistic data for testing of processing and migration tools applied in basaltic regions.

The work is based on a three-dimensional finite difference (FD) code, TIGER, made by SINTEF. The code computes wave propagation in a 3D anisotropic elastic media for micro seismic and traditional seismic sources.

The geological model is a basalt model, which covers from 24 km to 37 km of a real shot line in horizontal direction and from the water surface to 3.5 km depth.

The vertical parameter distribution is obtained from observations in two wells. At the depth of between 1100 m to 1500 m, a basalt horizon covers the whole sub surface layers.

The 2½D model has been constructed using the compound modelling software from Norsk Hydro.

The model is interpolated to a 3D grid. Each shot includes a subset of the global 3D grid in order to minimize computations.

- The computations were done on the IBM, p690 Regatta Turbo system at Parallab, University of Bergen (1.3 GHz Power4 processors), which consists of three 32 processor nodes with 64 Gbyte memory each (a total of 192 Gbyte memory). The system has a peak performance of 500 Gflops/s.
- The models applied here use about 12 Gbyte of memory. Each shot-model has 551 x 121 x 701 grid points. Temporal spacing in FD modelling is .25 ms and total recording time is 3 s.
- The wall-clock times have been measured for all 80 runs of the 3D wave propagation. Total time is smallest for 8 processors.

The seismic sections show clearly the wave propagation within and near the basalt layer.

Diffractions are observed near the fault, possibly due to the fault geometry

1025 m

1275m

- To study waveform propagation and generate realistic data for testing of processing and migration tools applied in e.g. basaltic regions
- The parallel code enables us to model large-scale realistic geological models in reasonable time
- FD does not identify phases

Rayleigh waves are P-SV-type waves, that exist in a half space (Aki and Richards, 1980; Stein, 1991), and may be written as:

Two conditions must be fulfilled for Rayleigh waves:

The resulting linear system of equations have non-trivial solutions for four values of the apparent (=phase) velocity cx, where only one can be used

Love waves are SH-type waves, that exist in a layer over a half space (Aki and Richards, 1980; Stein, 1991), and may be written as:

Four conditions must be fulfilled for Love waves:

This gives several (multi mode) solutions for the phase velocity cx, and each mode is a function of kx=/cx (dispersion):

Displacement due to P-SV and SH waves is

where Mij and Gij are the moment tensor and the Green’s tensor, respectively, and an index after comma means spatial derivative at the source. For symmetric moment tensor and multi-layered, anelastic, isotropic medium we have (see Haskell, 1953; Ben-Menahem and Harkrider, 1964; Harkrider, 1970; Panza, 1985; Panza and Suhadolc, 1987; Florsch et al., 1991; Engell-Sørensen, 1993; Engell-Sørensen and Panza, 2000)

Love Wave Eigenfunctions

With mode summation of dispersed modes for P-SV and SH-type waves all P-SV type waves with phase velocity less than the P-wave velocity in the half space are modelled.

Strategy for computing ground displacement :

- The dispersion relations for P-SV- and SH-type waves are found from the system of linear equations obtained by the boundary conditions at every interface and solved for the phase velocities (the eigenvalues) for all modes
- The associated eigenfunctions as function of depth, u(z), w(z), R(z), (z), v(z),L(z), are determined by solving the systems of linear equations obtained by the boundary conditions for the two types of waves
- The frequency domain ground displacements ur, uz, and u are found by an integral expression over all wave numbers, and can be obtained by summing the residue contributions, which are all the modes

Ocean-bottom registrations above an oil-field:

- Near-surface P- and S-velocities (anisotropy)
- Near-surface attenuation
- General seismic-source signature

Rays and wave surface in anisotropic, inhomogeneous media: moving eigenvector trihedral on ray (g1,g2,g3) at M (quasi-longitudinal,quasi-shear,quasi-shear). Fixed reference Cartesian system at O. Wave-surface trihedral at P (e1,e2,e3=N). Slowness vector s, wave number k, phase velocity vector V, and τ are parallel to N. Group velocity vector W is tangent to the ray at P.

The vector equation of motion for the spectral displacement, u, in general anisotropic and inhomogeneous elastic solid due to a point force F at r=r0 is (Ben-Menahem et al., 1991)

where is density, Cijkl(r) the tensor of elastic modules at r, angular frequency, uj components of displacement vector.

Green’s function is defined such that

uj = Gjm Fm

(Einstein summation convention). Substitution yields the wave equation of motion for Gjm

Assume the ”ray approximation” of the high-frequency asymptotic-series solution

where the sum is over all three types of waves propagating in the medium with unit polarization vectors,g, and the travel time between source and receiver is given by (). The polarization vectors are determined by an approximate equation obtained, when the ray-approximation Gjmdefined here is inserted in the vector equation of motion (the Christoffel equation).

The complex amplitude is given by (Červený et al., 1977, Aki and Richards, 1980, Hanyga et al.,1998)

V()=(sisi)–1/2 = phase velocity

R()(r0,r) = complete reciprocal reflection/transmission coefficient

(r0,r) = complete phase shift due to caustics along the ray

[J(r0)/J(r)]1/2 = geometrical spreading

J(r) =det ((x1,x2,x3)/(,q1,q2))=|(r/q1r/q2)·(r/)ray| = Jacobian of the transformation from Cartesian coordinates of a point on the ray to its ray coordinates

(r/)ray=group velocity

(,q1,q2) = ray coordinates, = travel time along the ray

(q1,q2) = curvilinear coordinates on the wave surface (=constant)

Let us assume we have a dipolar point source with moment. The medium response is (Aki and Richards, 1980, Ben-Menahem and Singh, 1981)

where the moment tensor M is a function of frequency.

By taking the spatial derivative of the ray approximation at the source the displacement (4) becomes

where s is slowness. By inserting Green’s tensor we finally get:

For seismic explosive sources we have Mmk=M·δmk and obtain the displacement component

Many ray-tracer codes:

- Recursive ray tracer (wave front tracing) (Hanyga et al., 1998, Moser and Pajchel, 1997;Vinje et al., 1993)
- Two-point ray tracing based on the eikonal equation that describes the kinematics of the wave field (e.g. Virieux, 1990, 1996; Hanyga and Pajchel, 1995; Engell-Sørensen, 1991)
- System of two coupled equations (the ray-equations). Can be solved by a Runge-Kutta scheme (Clarke, 1996).

The recursive ray tracer (Hanyga et al., 1998, Moser and Pajchel, 1997) can be used.

The geological model used is an n-layered 3D model with constant or linearly changing P- and S-velocities and density in the layers.

The ray tracer finds the attributes in a 3D grid for rays from a source.

Source (T1), Vertical Cross-Section along triangles, Light Blue: Structural Model

Ea

-------

Eb

Eb

-------

Ec

*

Ec

-------

Ed

Ekofisk Case

Ea

-------

Eb

Eb

-------

Ec

*

Ec

-------

Ed

Ekofisk Case

- Pre-stack and post-stack depth migration in complicated geological structures
- Tomography (i.e. inversion for geological structure)
- Ray tracing identify phases (i.e. remembers the history of individual rays)

- Apply impedance logs from a well
- Compute elastic impedance (EI) synthetic seismograms (e.g. 30 degrees) (Connolly, 1999)
- Stack synthetic EI seismograms (e.g. up till 30 degrees)
- Compare with migrated stacked CDP gathers (e.g. with offset up till 30 degrees)

- Lithology and pore fluid identification
- Calibration of migration results:
CDP-Gather Procedure for migration of data (GEOVECTEUR, SU, PROMAX, CREWES, e.t.c.):

- Decon (remove multiples)
- NMO-correction, stacking velocities, data
- Stack , one trace/CDP
- Migration with stacking velocities and stacked data -> migrated data

- Reflectivity method (Kind)
- Spectral methods (Bouchon, 1982)
- Pseudo Spectral Methods (Carcione)
- Finite element
- Spectral element methods (Caltech)

All methods can be applied but for different purposes in connection with depletion of an oil reservoir

The author would like to thank Norsk Hydro, Statoil, GEUS, and SINTEF for very helpful discussions and Parallab for being helpful with using the new IBM, p690 Regatta system

Phillips Petroleum Company, Norway and Norsk Hydro ASA are thanked for financing the work presented here

Aki K., and P. G. Richards (1980). Quantitative Seismology, Theory and Methods, vol. 1-2., W. H. Freeman, San Francisco, 931 pp.

Ben-Menahem A., R. L. Gibson Jr and A. G. Sena (1991). Green's tensor and radiation patterns of point sources in general anisotropic inhomogeneous elastic media, Geophys. J. Int. 107, 297-308.

Boore, D. M. (1972). Finite difference methods for seismic wave propagation in heterogeneous materials. In B. A. Bolt (editor), Seismology: Seismic Waves and Earth Oscillations (Methods in Computational Physics, Vol. 11). New York: Academic Press.

CREWES MATLAB: Margrave, G. F. (2001). Numerical Methods of Exploration Seismology with algorithms in Matlab, Department of Geology and Geophysics, The University of Calgary, http://www.crewes.org/Samples/EduSoftware

Engell-Sørensen, L. (2003). Optimized 3D Finite Difference Modelling of Basaltic Region, Extended Abstract, EAGE 65th Conference & Exhibition, Stavanger, 2003.

Engell-Sørensen, L. and J. Koster (2002). Optimization and Demonstration of 3D Finite Difference Tools for Modelling Seismic Acquisition: 64th EAGE Conference & Exhibition, Abstract.

Hokstad, K., L. Engell-Sørensen, and F. Maaø (2002). 3-D elastic finite-difference modelling in tilted transversely isotropic media: 72nd Ann. Internat. Mtg: Soc. of Expl. Geophys.

Holberg, O. (1987). Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena, Geophys. Prosp. 35, 629-655.

Mittet R., and L. Amundsen and B. Arntsen (1994). Iterative inversion/migration with complete boundary conditions for the residual misfit field, J. of Seismic Exploration, 1994, 3, p. 141-156. Petersen, S. A. (1999). Compound modelling – A geological approach to the construction of shared earth models: 61st EAGE Conference & Technical Exhibition, Abstract.

SU: Stockwell, J. W. and J. K. Cohen (1998). The New SU User’s Manual, Version 2.5, Center for WavePhenomena, Colorado School of Mines, http://www.cwp.mines.edu/cwpcodes/

Aki K., and P. G. Richards (1980). Quantitative Seismology, Theory and Methods, vol. 1, W. H. Freeman, San Francisco, 557 pp.

Ben-Menahem, A., Harkrider, D. G. (1964). Radiation patterns of seismic surface waves from buried dipolar point sources in a flat stratified Earth. J. Geophys. Res. 69, 2605-2620

Engell-Sørensen, L. and G. F. Panza (2000). Inversion for the Moment Tensor and Source Location Using Love- and Rayleigh-Type Waves, in P. C. Hansen, B. H. Jacobsen & K. Mosegaard (Eds.), Methods and Applications of Inversion, Lecture Notes in Earth Sciences 92, Springer, Berlin, 2000.

Engell-Sørensen, L. (1993). North Sea earthquake source parameters, Ph D. Thesis, Institute of Solid Earth Physics, University of Bergen, Norway. March 1993.

Florsch, N., Fäh, D., Suhadolc, P., and Panza, G. F. (1991). Complete Synthetic Seismograms for High-Frequency Multimode SH-waves. Pageph 136, 4, 529-560.

Herrman, R. (2002). Computer Programs in Seismology – Version 3.20, http://www.eas.slu.edu/People/RBHerrmann/CPS32.html

Harkrider, D. G., (1970). Surface waves in multilayered elastic media. Part II. Higher mode spectra and spectral ratios from point sources in plane layered earth models, Bull. Seism. Soc. Am. 60, 6, 1937-1987.

Haskell, N. A. (1953). Dispersion of surface waves on multilayered media, Bull. Seism. Soc. Am. 43,17-34.

Panza, G. F. (1985). Synthetic seismograms: the Rayleigh waves modal summation. J. Geophys. 58, 125-145.

Panza, G. F., Suhadolc, P. (1987). Complete strong motion synthetics. Ed. Bolt, A. B., Academic Press, Orlando, FL., pp. 135-204.

Stein, S. and M. Wysession (2002). Introduction to Seismology, Earthquakes and Earth Structure, Blackwell Science Inc; 1st ed., 498 pp.

Aki K. and P. Richards (1980). Quantitative Seismology, Theory and Methods, vol. 1, W. H. Freeman, San Francisco, 557 pp.

Ben-Menahem A., R. L. Gibson Jr and A. G. Sena (1991). Green's tensor and radiation patterns of point sources in general anisotropic inhomogeneous elastic media, Geophys. J. Int. 107, 297-308.

Ben-Menahem A. and S. J. Singh (1981). Seismic Waves and Sources, Springer-Verlag, New York, 1108 pp.

Clarke R. A. (1996). Two point raytracing in 3D blocky models: KIM 1996 Annual Report, Istitut Français du Pétrole, Pau, France.

Červený, V., I. A. Molotkov, and I. Pšenčik (1977). Ray Methods in Seismology, University Karlova, Praha.

Engell-Sørensen, L. (1991). Inversion of arrival times of microearthquake sources in the North Sea using a 3-D velocity structure and prior information. Part I. Method. Bull. Seism. Soc. Am. 81 1183-1194.

Hanyga A. and Pajchel J.(1995). Point-to-curve ray tracing in complicated geological models, Geophysical Prospecting 43, 859-872.

Hanyga, A., A. B. Druzhinin, A. D. Dzhafarov and L. Frøyland (1998). 3-D Recursive Cell Ray Tracing in Inhomogeneous and Weakly Anisotropic Discontinuous Media. Norsk Hydro Report.

Moser, T. J. and J. Pajchel (1997). Recursive seismic ray modelling application in inversion and VSP, Geophysical Prospecting, 45, 885-908.

NORSAR-2D Ray Modelling (2003). http://www.norsar.no/Seismod/Products/N2D.html

NORSAR-3D Ray Modelling (2003). http://www.norsar.no/Seismod/Products/N3D.html

SEIS88 (2003). Seismic Waves in Complex 3-D Structures (SW3D), Department of Geophysics, Charles University Prague, http://seis.karlov.mff.cuni.cz/software/index.htm

Vinje V., Iversen E. and Gjøystdal H. (1993). Traveltime and amplitude estimation using wavefront construction. Geophysics 58, 1157-1166.

Virieux J. (1990). Propagation in inhomogeneous media: Ray theory workshop, Institut de Géodynamique, Université de Nice.

Virieux J. (1996). Seismic Ray Tracing in Seismic Modelling of Earth Structure, revue OCEANIS, Institut Océanographique, Paris, vol. 23 (2), 153-205.

Aki K., and P. G. Richards (1980). Quantitative Seismology, Theory and Methods, vol. 1., W. H. Freeman, San Francisco, 557 pp.

Conolly, P. (1999). Elastic impedance, The Leading Edge, April 1999.

CREWES MATLAB: Margrave, G. F. (2001). Numerical Methods of Exploration Seismology with algorithms in Matlab, Department of Geology and Geophysics, The University of Calgary, http://www.crewes.org/Samples/EduSoftware

Geovecteur®Plus 7.1. (2003). CGG Data Processing Software, http://www.cgg.com/proserv/software/products/geovecteur/GVT1.html

PROMAX (2002) Seismic Processing Software, Landmark A Halliburton Company, http://www.lgc.com/news/pressreleases/20020819-landmark+releases+seisspace.htm

SU: Stockwell, J. W. and J. K. Cohen (1998). The New SU User’s Manual, Version 2.5, Center for Wave Phenomena, Colorado School of Mines, http://www.cwp.mines.edu/cwpcodes/