## Section: Scientific Foundations

### Numerical methods

In the 1990s, numerical modeling techniques that made it possible to study large 3D problems started to emerge. However, many research groups continued to carry out two-dimensional calculations because of the computational cost of these methods in 3D, and of the difficulty of implementing them (one needs to run them on large multiprocessor machines, based on specific programming techniques such as multitasking). In the last thirty years, several methods have been used for the numerical calculation of synthetic seismograms in complex geological models, first in two dimensions, and more recently in three dimensions. The finite-difference technique [44] is the most popular, and was applied successfully to local or regional models [38] , [45] . Another largely used technique is the pseudospectral method, which uses global bases of Chebyshev or Legendre polynomials [25] . However, in many cases of practical interest, these traditional methods suffer from limitations such as numerical dispersion or numerical anisotropy. We thus have in recent years developed a new technique, called the spectral element method, which had been introduced initially in fluid dynamics [41] , and that we applied for the first time to the propagation of waves in 3D structures (see for example [35] , [34] . This work showed the superiority of the spectral element method over more traditional numerical techniques in terms of precision, weak numerical dispersion, and geometrical flexibility making it possible to adapt it to large and complex 3D models [32] , [33] . Such precision is a crucial advantage for the resolution of forward or inverse problems in seismology.

The spectral-element technique developed by Dimitri Komatitsch and his coworkers is based on a variational formulation of the wave equation, and combines the flexibility of a finite element method with the precision of a global pseudospectral method. The finite element grid is adapted to all major discontinuities of the geological model. In order to maintain a relatively constant grid resolution in the entire model in terms of the number of grid points per wavelength, and to reduce the computing time, the size of the elements is decreased with depth in a geometrically-conforming fashion, which allows us to preserve an exactly diagonal mass matrix in the method. The effects of attenuation and anisotropy are taken into account in the technique.

We applied this technique to a large number of real geophysical cases of practical interest, for example the study of strong ground motion and of the associated seismic risk in the densely-populated Los Angeles basin region. This area consists of a basin of great dimension (more than 100 km x 100 km) which is one of the deepest sedimentary basins in the world (the sedimentary layer has a maximum thickness of 8.5 km right underneath downtown Los Angeles), and therefore one of the most dangerous because of the resulting amplification of the seismic waves. In the case of a small recent earthquake in Hollywood (of magnitude MW = 4.2 on September 9, 2001), well recorded by more than 140 stations of the TriNet seismic network of Southern California, we managed for the first time to fit the three components of the vector displacement, while most of the previous studies concentrated on the vertical component only [39] , [46] , [40] , while still obtaining a good fit to the recorded data down to relatively short periods (2 seconds). This study clearly showed how useful sophisticated 3D numerical modeling techniques can be in such a context [31] .

Topography also plays a significant role for the characteristics of the "ground roll", i.e., surface waves recorded by the oil industry in field acquisition experiments, and its effects are essentially three-dimensional. We will use the 3D spectral-element code (SPECFEM3D) developed by D. Komatitsch and his coworkers to generate synthetic data for various 2D and 3D configurations. Some recent articles show that the presence of heterogeneities in the subsurface not only contributes to attenuate the recorded signal but also to delay coherent events. It is also difficult, under these conditions, to distinguish between a distribution of heterogeneities and a mere stratification of layers. We will also carry out 3D simulations with complex topographies. Effects of amplification due to conversions of body waves into Rayleigh surface waves in regions of sharp topography have been observed in the field. Amplitude variations of the signal of an order of magnitude can appear along the recording antenna [43] , [37] . Because of diffraction phenomena, complex arrivals are also observed. We will try to reproduce these effects with our 3D numerical simulations of seismic wave propagation.

A few years ago, Dimitri Komatitsch started to work on this topic of the numerical modeling of the effect of topography on seismic wave propagation with Jacques Muller (TOTAL) and Patrice Ricarte (IFP) within the framework of the "Foothills" project of TOTAL. At the time, calculations were carried out exclusively in 2D because of the computer resources available. A very significant problem for the oil industry is indeed to study models located in foothill basins, in which the topography and the presence of the weathered surface layer (Wz) plays a crucial role on the quality of the seismic data recorded. It is now necessary to generalize such calculations, and obviously nowadays to carry them out in 3D using the seismic modeling tools that we have developed. Roland Martin is currently in the process of using such techniques to apply it to another real case in South America in which the very thin weathered zone makes the simulation very difficult to perform at high-frequency if classical numerical simulation techniques such as staggered finite-difference methods are used.

Logo Inria