Gregory A. Newman Earth Sciences Division Lawrence Berkeley National Laboratory

1. A New Approach to Joint Imaging of Electromagnetic and Seismic WavefieldsInternational Symposium on Geophysical Imaging with Localized WavesSanya, Hainin Island ChinaJuly 24-28, 2011 Gregory A. Newman Earth Sciences Division Lawrence Berkeley National Laboratory

2. PRESENTATION OVERVIEW • Controlled Source EM and Magnetotellric Data Acquisition • Motivation for Joint Imaging • Formulation of the Joint Inverse Problem • Large Scale Modeling Considerations • Need for High Performance Computing • Joint MT/CSEM Imaging Results • Gulf of Mexico (synthetic example) • Joint Imaging - EM & Seismic Data • Issues & Proposed Methodology • Laplace-Fourier Wave Field Concepts • Similarities and Differences to the EM Wave Field • Imaging Across Multiple Scale Lengths • 3D Elastic Wave Fields • Preconditioners and solvers for 3D Seismic Wave Field Simulations • Conclusions

3. Marine CSEM & MT Surveying • CSEM • Deep-towed Electric Dipole transmitter • ~ 100 Amps • Water Depth 1 to 7 km • Alternating current 0.01 to 3 Hz • ‘Flies’ 50 m above the sea floor • Profiles 10’s of km in length • Excites vertical & horizontal currents • Depth of interrogation ~ 3 to 4 km • Sensitive to thin resistive beds • MT • Natural Source Fields • Less than 0.1 Hz • Measured with CSEM detectors • Sensitive to horizontal currents • Depth of interrogation 10’s km • Resolution is frequency dependent • Sensitive to larger scale geology

4. j = a S {(djobs - djp)/ej}2 + b S {(Zjobs - Zjp)/pj}2 + lhmhWT Wmh+ lvmvWT Wmv s.t. mv mh dobs and dp are N observed and predicted CSEM data Zobs and Zp are M observed and predicted MT impedance data e & p = CSEM and MT data weights mh = horizontal conductivity parameters mv = vertical conductivity parameters W = Ñ2 operator; constructs a smooth model lh&lv = horizontal & vertical tradeoff parameters a &b = scaling factors for CSEM and MT data types JOINT 3D INVERSE MODELING Minimize: N M j=1 j=1

5. LARGE-SCALE 3D MODELINGCONSIDERATIONS • Require Large-Scale Modeling and Imaging Solutions • 10’s of million’s field unknowns (fwd problem) • Solved with finite difference approximations & iterative solvers • Imaging grids 400 nodes on a side • Exploit gradient optimization schemes, adjoint state methods • Parallel Implementation • Two levels of parallelization • Model Space (simulation and inversion mesh) • Data Space (each transmitter/MT frequency - receiver set fwd calculation independent) • Installed & tested on multiple distributed computing systems; 10 – 30,000 Processors • Above procedure satisfactory except for very largest problems • To treat such problems requires a higher level of efficiency • Optimal Grids • Separate inversion grid from the simulation/modeling grid • Effect: A huge increase in computational efficiency ~ can be orders of magnitude

6. Optimal Grids 100 km 10 km 10 km 100 km Ωm imaging grid Ωs simulation grid sail lines

7. GRID SEPARATION EFFICIENCIES • Advantages • Taylor an optimal simulation grid Ωs for each transmitter-receiver set • Inversion grid Ωm covers basin-scale imaging volumes at fine resolution • Simulations grids much smaller, a subset of the imaging grid • Faster solution times follow from smaller simulation grids • What’s Required • A mapping of conductivity from Ωm to Ωs & Ωs to Ωm • Conductivity on Ωs edged based • Conductivity on Ωm cell based • An appropriate mixing law for the conductivity mappings

8. Joint CSEM - MT Imaging Mahogany Prospect, Gulf of Mexico • Study: 3D Imaging of oil bearing horizons with complex salt structures present • Simulated Example: 100 m thick reservoir, 1 km depth, salt below reservoir • Model: 0.01 S/m salt, 2 S/m seabed, 0.05 S/m reservoir, 3 S/m seawater • MT Data: 7,436 data points, 143 stations & 13 frequencies 0.0005 to 0.125 Hz • CSEM Data: 12,396 data points, 126 stations & 2 frequencies 0.25 and 0.75 Hz • Starting Model: Background Model without reservoir or salt • Processing Times: 5 to 9 hours, 7,785 tasks, NERSC Franklin Cray XT4 System Survey Layout y=5 km cross section -5 0 x(km) 5 10 -10

9. JOINT CSEM-MT IMAGING:The Benefits

10. Joint CSEM - MT Imaging Mahogany Prospect Gulf of Mexico

11. Joint Imaging of EM and Seismic Data Issues • Rock Physics Model • links attributes to underlying hydrological parameters • too simplistic • difficult or impossible to define robust/realistic model • Differing Resolution in the Data • EM data 10x lower resolution compared to seismic • RTM & FWI of Seismic Data • requires very good starting velocity model • velocity can be difficult or impossible to define • huge modeling cost due to very large data volumes (10,000’s of shots; 100,000’s traces per shot)

12. Joint Imaging of EM and Seismic Data A way forward • Abandon Rock Physics Model • assume conductivity and velocity structurally correlated • employ cross gradients: t =    • t = 0 =>   ;  = 0 and/or  =0 • Equalize Resolution in the Data • treating seismic and EM data on equal terms • Laplace-Fourier transform seismic data – Shin & Cha 2009

13. Acoustic Wave Equation Propagating Wave Fourier/Frequency Domain Time Domain Damped Diffusive Wave Laplace/Fourier Domain At first glance similar physics & similar resolution with EM fields skin depth:

14. Seismic Imaging: Laplace-Fourier Domain 337 shot gathers 151 detectors/shot maximum offsets 15km BP Salt Model Starting Velocity Model s = 10.5 to 0.5 =0.5 Laplace Image Standard FWI Image s = 10.5 to 0.5 f = 6 to 0.5 =0.5 New FWI Image Laplace-Fourier Image Taken from Shin & Cha, 2009

15. Laplace-Fourier Wavefield Modeling tol= There are differences compared to EM fields • wavelength and skin depth are decoupled Meshing Issues to Consider • grid points per wavelength:10 points – 2nd order accuracy • grid points per skin depth: 6 points – 2nd order accuracy Accuracy Issues • wavefield dynamic range extreme ~ 70 orders • iterative Krylov solvers require tiny solution tolerances

16. LAPLACE-FOURIER IMAGING:Mahogany Prospect • Study: 3D Imaging of oil bearing horizons with complex salt structures present • Simulated Example: 100 m thick reservoir, 1 km depth, salt below reservoir • Model: 6 km/s salt, 3 km/s seabed, 2 km/s reservoir, 1.5 km/s seawater • Seismic Data: 24,577 data points, 287 stations at 1 frequency (σ=1, ω=2π) • Starting Model: Background Model without reservoir or salt • Processing Times: 10.3 hours, 6,250 tasks, NERSC Franklin Cray XT4 System * Survey line N 10000 km 287 sources, (σ=1, ω=2π) Source & Receiver & Spacing 1 km & 300 m Max. offsets 17 km 50 m below sea surface Misfit Survey line N 7500 km Survey line N 5000 km Survey line N 2500 km Survey line N 0 km Survey line N -2500 km Survey line N -5000 km West-East -20 km 20 km

17. LAPLACE-FOURIER IMAGE:Mahogany Prospect -0.9 km -0.9 km 0 km North 0 km North 13.5 km 13.5 km 2.5 km North 2.5 km North 5.0 km North 5.0 km North 7.5 km North 7.5 km North 10.0 km North 10.0 km North 12.0 km North 12.0 km North -18.5 km 18.5 km -18.5 km 18.5 km

18. LAPLACE-FOURIER IMAGE:Mahogany Prospect1kmbelow seabed 12.0 km North 12.0 km North -5.0 km North -5.0 km North 18.5 km 18.5 km -18.5 km -18.5 km West – East West – East

19. Joint EM-Seismic ImagingProblem Formulation and are N observed and predicted EM data and are M observed and predicted Laplace-Fourier seismic data • and= EM and seismic data weights = m conductivity parameters • = m acoustic velocity parameters = Ñ2 operator; constructs a smooth model • and = conductivity & velocity tradeoff parameters • and = scaling factors forEM and seismic data types • are cross gradient structural constraints; is a penalty parameter

20. Recipe for Auxiliary Parameters First carry out separate inversions for seismic and EM => choose smoothing parameters (cooling approach) Consider total objective functional : Next balance data funtionals for seismic and EM : => set => rescale accordingly Test out values for => selected balances => trail values tested out over a few inversion iterations

21. Initial Imaging Resultsmarine example CSEM 16 km max. offsets 17 shots 161 detectors/shot conductivity image correlated with velocity; =1011 conductivity image no correlation to velocity; =0 f = 0.25 seismic 12-16 km offsets 85 shots 121-161 detectors/shot s =5,4,3,2,1 f = 0 velocity image correlated with conductivity; =1011 velocity image no correlation to conductivity; =0 Computational Requirements: 4250 cores – Franklin NERSC Processing Time: ~22 hours

22. Elastic Wave Field Simulator First- order system for velocity –stress components Laplace-Fourier Domain - velocity components, - stress components,  - density,  and  - Lame coefficients. Forces are defined via Moment-Tensor components (R. Graves 1996)

23. Boundary and Initial Conditions The ordinary initial conditions for all components are zeros. The boundary conditions are • a) PML absorbing boundary conditions for velocity • b) free surface boundary for normal stress component

24. Solution Realization Iterative Krylov Methods Coupled System - matrices of FD first derivative operators System transformed : solve only for the velocity components D is complex non-symmetric 15 diagonals for 2nd order scheme 4th order scheme 51 diagonals

25. Solution Accuracy

26. Two Half Spaces Model Test

27. SEG/EAEG SALT MODEL TEST real 3D salt body 200x200x130 nodes y=4800 m s=(3+18.85i) sec-1 y=4800 m. imaginary y=4800 m s=(3+18.85i) sec-1 Snap-shots of velocity field y-component Vertical cross section of velocity (p) Point source at r=(800,800,500), m ). Solution Time 1355 sec, 576 Iterations Solution Tolerance 1e-7

28. Laplace-Fourier Transformation Benefits • Wave field simulations • excellent choice for a preconditioner(frequency domain ) On a class of preconditioners for solving the Helmholtzs equation: Erlangga et al., 2004, Applied Numer. Mathematics, 50 409-425. • Imaging • possibilities to image at multiple scales and attributes A consistent Joint EM seismic imaging approach • known to produce robust macro-models of velocity Critical to successful RTM and FWI of seismic reflection data

29. Conclusions Demonstrated Benefits Massively Parallel Joint Geophysical Imaging • Joint CSEM & MT • acoustic seismic (Laplace-Fourier Domain) • HPC essential Future Developments in LF elastic wavefield imaging • massively parallel (MP) LF elastic wavefield simulator • exploit simulator as a preconditioner • frequency domain wavefield modeling • exam MP direct solvers • gradient based 3D LF elastic imaging code Future Plans in Joint Imaging • joint EM and elastic wavefield 3D imaging capability

30. ACKNOWLEDGEMENTS My Colleagues: M. Commer, P. Petrov and E. Um Research Funding US Department of Energy Office of Science Geothermal Technologies Program

31. Computational Details