Applied Seismic Simulation

1. 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

2. Introduction 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)

3. Contents: Theory, Applications, Areas of Interest for: • 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)

4. Applied Seismic Simulation Codes (by L. Engell-Sørensen) nc=non commercial=academic use

5. Method, Finite Difference 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).

6. Representation Result 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.

7. General TI Medium 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)

8. FD Scheme: Higher Order Diff and Shift Operators (Staggered Grid) Newton’s law: Hooke’s law:

9. Parallel Method 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)

10. Data Exchange • 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

11. Applied Finite Difference: Basalt Model in 3D Grid 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.

12. P Velocity

13. Shot Geometry and Source Signal

14. Performance • 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.

15. Common Offset Gathers of Traction (z) for 80 shots of Displacement Source (z) 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

16. 1025 m

17. 1275m

18. 1775 m

19. Finite Difference, Areas of Interest • 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

20. Method, Surface Waves 1 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

21. Method, Surface Waves 2 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):

22. Normal Mode Summation 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)

23. Outward Radial P-SV Component

24. Upward Vertical P-SV Component

25. Counter-Clockwise SH Component

26. Rayleigh Wave Eigenfunctions

27. Rayleigh Wave Parameters

28. Love Wave Eigenfunctions

29. Love Wave Parameters

30. Applied Surface Waves 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

31. Surface Waves, Areas of Interest Ocean-bottom registrations above an oil-field: • Near-surface P- and S-velocities (anisotropy) • Near-surface attenuation • General seismic-source signature

32. Method, Ray-Tracing 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.

33. Radiation Pattern of Point Sources 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.

34. Green’s Function Green’s function is defined such that uj = Gjm Fm (Einstein summation convention). Substitution yields the wave equation of motion for Gjm

35. Ray-approximation 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).

36. Geometrical Spreading 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/q1r/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)

37. Dipolar Point Source With Moment 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.

38. Ray Approximation Displacement 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:

39. Ray Approximation Displacement for Seismic Explosive Source For seismic explosive sources we have Mmk=M·δmk and obtain the displacement component

40. Applied Ray-Tracing 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).

41. Recursive Ray Tracer 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.

42. Grid Locations West-East

43. Ray-Traced Attributes Source (T1), Vertical Cross-Section along triangles, Light Blue: Structural Model

44. Initial P-Polarity in X-Direction Ea ------- Eb Eb ------- Ec * Ec ------- Ed Ekofisk Case

45. P-Polarity in X-Direction Ea ------- Eb Eb ------- Ec * Ec ------- Ed Ekofisk Case

46. Ray-Tracing, Areas of Interest • 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)

47. Applied 1D Methods • 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)

48. 1D Methods, Areas of Interest • 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

49. Other Methods for Full Waveform Modelling • Reflectivity method (Kind) • Spectral methods (Bouchon, 1982) • Pseudo Spectral Methods (Carcione) • Finite element • Spectral element methods (Caltech)

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