Optimized Finite Difference Methods for Seismic Acoustic Wave Modeling

The finite difference (FD) methods are widely used for approximating the partial derivatives in the acoustic/elastic wave equation. Grid dispersion is one of the key numerical problems and will directly influence the accuracy of the result because of the discretization of the partial derivatives in the wave equation. Therefore, it is of great importance to suppress the grid dispersion by optimizing the FD coefficient. Various optimized methods are introduced in this chapter to determine the FD coefficient. Usually, the identical staggered grid finite difference operator is used for all of the first-order spatial derivatives in the first-order wave equation. In this chapter, we introduce a new staggered grid FD scheme which can improve the efficiency while still preserving high accuracy for the first-order acoustic/elastic wave equation modeling. It uses different staggered grid FD operators for different spatial derivatives in the first-order wave equation. The staggered grid FD coefficients of the new FD scheme can be obtained with a linear method. At last, numerical experiments were done to demonstrate the effectiveness of the introduced method

Modelling Seismic Wave Propagation for Geophysical Imaging

International audienceThe Earth is an heterogeneous complex media from the mineral composition scale (10−6m) to the global scale ( 106m). The reconstruction of its structure is a quite challenging problem because sampling methodologies are mainly indirect as potential methods (Günther et al., 2006; Rücker et al., 2006), diffusive methods (Cognon, 1971; Druskin & Knizhnerman, 1988; Goldman & Stover, 1983; Hohmann, 1988; Kuo & Cho, 1980; Oristaglio & Hohmann, 1984) or propagation methods (Alterman & Karal, 1968; Bolt & Smith, 1976; Dablain, 1986; Kelly et al., 1976; Levander, 1988; Marfurt, 1984; Virieux, 1986). Seismic waves belong to the last category. We shall concentrate in this chapter on the forward problem which will be at the heart of any inverse problem for imaging the Earth. The forward problem is dedicated to the estimation of seismic wavefields when one knows the medium properties while the inverse problem is devoted to the estimation of medium properties from recorded seismic wavefields

Parsimonious finite-volume frequency-domain method for 2D P-SV-wave modeling

International audienceA new numerical technique for solving the 2D elastodynamic equations based on a finite volume approach is proposed. The associated discretization is through triangles. Only fluxes of required quantities are shared between cells, relaxing meshing conditions compared to finite element methods. The free surface is described along the edges of the triangles which may have different slopes. By applying a parsimonious strategy, stress components are eliminated from the discrete equations and only velocities are left as unknowns in triangles, minimizing the core memory requirement of the simulation. Efficient PML absorbing conditions have been designed for damping waves around the grid. Since the technique is devoted to full waveform inversion, we implemented the method in the frequency domain using a direct solver, an efficient strategy for multiple-source simulations. Standard dispersion analysis in infinite homogeneous media shows that numerical dispersion is similar to those of O(¢x2) staggeredgrid finite-difference formulations when considering structured triangular meshes. The method is validated against analytical solutions of several canonical problems and with numerical solutions computed with a well-established finite-difference time-domain method in heterogeneous media. In presence of a free surface, the finite-volume method requires ten triangles per wavelength for a flat topography and fifteen triangles per wavelength for more complex shapes, well below criteria required by the staircase approximation of finite-difference methods. Comparison between the frequency-domain finite-volume and the O(¢x2) rotated finite-difference methods also shows that the former is faster and less-memory demanding for a given accuracy level. We developed an efficient method for 2-D P-SV-wave modeling on structured triangular meshes as a tool for frequency-domain full-waveform inversion. Further work is required to assess the method on unstructured meshes

Optimized viscoelastic wave propagation for weakly dissipative media

The representation of viscoelastic media in the time domain becomes more challenging with greater bandwidth of the propagating waves and number of travelled wavelengths. With the continuously increasing computational power, more extreme parameter regimes become accessible, which requires the reassessment and improvement of the standard ‘memory variable' methods to implement attenuation in time-domain seismic wave-propagation methods. In this paper, we propose a method to minimize the error in the wavefield for a fixed complexity of the anelastic medium. This method consists of defining an appropriate misfit criterion based on a first-order analysis of how errors in the discretized medium propagate into errors in the wavefield and a simulated annealing optimization scheme to find the globally optimal parametrization. Furthermore, we derive an analytical time-stepping scheme for the memory variables that encode the strain history of the medium. Then we develop the coarse grained memory variable approach for the spectral element method (SEM) and benchmark it using the 2.5-D code AxiSEM for global body waves up to 1 Hz. Showing very good agreement with a reference solution, it also leads to a speedup of a factor of 5 in the anelastic part of the code (factor 2 in total) in this 2.5-D approach. A factor of ≈15 (3 in total) can be expected for the 3-D case compared to conventional implementation

Effects of seismic anisotropy and attenuation on first-arrival waveforms recorded at the Asse II nuclear repository

For decades, deep geological storage in former salt mines has been a widely recognized strategy for long-term radioactive waste disposal. However, in the case of the Asse II repository in Lower Saxony, groundwater inflow and instabilities in the geological structures rendered the mine unusable as a long-term solution. The nuclear waste needs to be recovered on grounds of safety reasons, hence the need for detailed structural information in order to build a new shaft. In this context, it is essential to use optimized, modern seismic imaging methods, such as, for instance, full-waveform inversion (FWI), to obtain high-resolution, physical parameter models of the Asse salt structure and its surroundings. The goal of this thesis is to draw conclusions on the future application of elastic FWI using first-arrival waveforms at frequencies up to 20 Hz, potentially including anisotropy and attenuation. For this purpose, simple parameter models were created based on previously known geological information and used as references for synthetic forward modeling tests. The objectives were (a) to see if the models are suitable as initial models for FWI, (b) to assess what type of anisotropy needs to be considered, if at all, and (c) to investigate the significance of attenuation. To facilitate the numerical tests, the mathematics of viscoelastic anisotropic wave propagation was studied and a new 2D finite-difference (FD) anisotropic forward solver was implemented. A detailed comparison of wavefield snapshots and seismograms was conducted between isotropic, vertical-transverse isotropic (VTI), and tilted-transverse isotropic (TTI), as well as elastic and viscoelastic modeling. The results demonstrate that, in general, the models are likely to meet the prerequisites for the successful application of first-arrival FWI up to frequencies of about 20 Hz. While attenuation turned out to be only a minor factor, it is, however, essential to incorporate anisotropy. As the Asse salt structure is complex and steeply dipping, TTI modeling is the preferred way to correctly map the subsurface in high resolution and match first-arrival traveltimes. Furthermore, a comparison with field data acquired over the Asse hill shows that many features present in that data can already be explained using the current approach

Modélisation numérique de la propagation d'ondes en milieux complexes : application aux milieux granulaires non consolidés

Le sous-sol, qui contient de nombreuses ressources naturelles (eau, gaz, pétrole, etc.), peut également constituer un risque naturel en raison de ses caractéristiques lithologiques et topographiques. Par ailleurs, dans le contexte du changement climatique, il devient de plus en plus important d'estimer le taux de saturation des fluides dans ces milieux pour prévenir les catastrophes naturelles comme les glissements de terrain ou des inondations. Toutes ces raisons suscitent l'intérêt des géophysiciens qui cherchent à mieux comprendre la proche surface et donc à la caractériser. En Géophysique, différentes techniques sont utilisées pour caractériser le sous-sol parmi lesquelles des techniques sismiques non destructives. Lorsque les ondes sismiques traversent un matériau donné, elles sont diffractées, réfléchies ou converties et contiennent ainsi des informations sur les phases fluide et solide. Pour mieux comprendre les mesures acoustiques et sismiques dans les sédiments et les sols, de nombreuses études sur les milieux granulaires non consolidés ont été menées in situ et aussi à l'échelle du laboratoire où des modèles théoriques ont été développés. Dans cette thèse, nous souhaitons modéliser des milieux granulaires qui sont un type de milieu complexe difficile à caractériser. Pour atteindre cet objectif, nous avons suivi trois étapes. Premièrement, nous avons développé un outil numérique qui calcule l'ensemble du champ d'ondes d'un modèle élastique bidimensionnel avec des structures complexes. Nous proposons une méthode de volumes finis basée sur un solveur de Riemann (RFV-FSP/Riemann Finite Volume-Fluxes Frequency Shift PML method) pour calculer les champs d'ondes sismiques sur des grilles colocalisées ainsi qu'une formulation des conditions absorbantes de type PML spécifiquement conçue pour la méthode des volumes finis. Ces dernières sont optimisées à incidence rasante en utilisant des formulations convolutives avec décalage en fréquence (C-PML) ou non convolutives (ADE-PML). Ici, elles sont appliquées aux dérivées spatiales des flux, ce qui diffère des PML classiques qui sont généralement appliquées aux dérivées spatiales des variables primitives (vitesses et contraintes des particules). La méthode des volumes finis et les différents types de conditions aux limites sont test´es et validés sur différents cas synthétiques hétérogènes. Les volumes finis sont comparés à d'autres techniques comme les différences finies et les éléments finis d'ordre élevé. Nous appliquons aussi notre méthode à une configuration de couplage fluide-solide et à quelques modèles sismiques d'intérêt dans le contexte de milieux granulaires non consolidés présentant de fortes variations de propriétés avec la profondeur. En particulier nous concentrons notre attention sur la résolution numérique des ondes de surface comme les ondes de Rayleigh. Pour obtenir plus de précision, nous avons implémenté un schéma spatial décentré du quatrième ordre proche de la surface libre. Deuxièmement, nous avons mis en place des outils de traitement du signal qui détectent les temps des premières arrivées sismiques, et calculent les courbes de vitesse de phase et les modes de propagation des ondes. Ces derniers outils sont utilisés pour l'analyse de dispersion. Pour finir, nous revisitons une étude réalisée sur des milieux granulaires non consolidés à l'échelle du laboratoire en utilisant les différents outils développés. Nous comparons différents modèles (2D ou 3D) avec différentes rhéologies (élastique ou poro-élastique), différentes conditions aux limites (PML ou Dirichlet) et différentes modélisations numériques de la source (point source ou pot vibrant) afin de reproduire les données expérimentales. L'étude de la sensibilité des données sismiques à l'emplacement de la source était également cru- ciale pour améliorer l'amplitude des signaux et la détection des différents modes sismiques. Cela nous permettra à l'avenir de mieux imager et comprendre ces milieux complexes.The subsurface, which contains many natural resources (water, gas, oil, etc.), can also constitute a natural risk because of its lithological and topographical characteristics. In the context of climate change, it becomes more and more important to estimate the rate of saturation of fluids in these media to prevent natural disasters like landslides or flash floods. All these reasons arouse the in- terest of geophysicists who seek to better understand the near surface and therefore to characterize it. In geophysics, different techniques are used to characterize the subsurface among them seismic techniques which are non-destructive. When seismic waves are crossing a given material, they are diffracted, reflected or converted and thus contain information on fluid and solid phases. To better understand acoustic and seismic measurements in sediments and soils, many studies on unconsol- idated granular media have been conducted in situ, and at the laboratory scale where theoretical models have been developed. In this thesis, we want to model granular media which are a type of complex medium difficult to characterize. To achieve this objective, we followed three steps. First, we developed a numerical tool which calculates the entire wave field of a two dimensional geometric elastic model with complex structures. And we compare its accuracy to other techniques like the classical staggered-fine difference or the high-order spectral element methods. We propose a finite volume method based on a Riemann solver (RFV-FSP/Riemann Finite Volume-Fluxes frequency Shift PML method) to compute seismic wave fields on collocated grids as well as a formulation of perfectly matched layer (PML) absorbing boundary conditions that are more specifically designed to the finite volume method. The PML boundary conditions are optimized at grazing incidence by using frequency shift convolutional (C-PML) or non convolutional formulations (ADE-PML). Here, they are applied to the spatial fluxes derivatives, which is a different formulation than clas- sical PMLs that are generally applied to the spatial derivatives of the primitive variables (particle velocities and stresses). The finite volume method and the different kinds of boundary conditions are tested and validated on different heterogeneous synthetic cases. The finite volume method is compared to other techniques like finite differences and high order finite elements. Finally, we apply our method to a fluid-solid coupling configuration and to some seismic models of interest in the context of unconsolidated granular media presenting sharp property variations with depth. In particular we focus our attention on the implementation of the numerical resolution of surface waves like the Rayleigh waves, which is not trivial with classical staggered finite differences. We thus implemented a non-centered fourth-order spatial scheme at the free surface to achieve more accuracy. Second, we implemented signal processing tools that calculate phase velocity curves and detect first arrival travel times and wave propagation modes of seismic data. These tools are used for dispersion analysis. Third, we revisit a study carried out on unconsolidated granular media at the laboratory scale using the different tools (finite differences or finite volumes). We compare different models with different rheologies (elastic or poro-elastic), different dimensions (3D or 2D), different boundary conditions (PML or Dirichlet) and different numerical modeling of the source (stick or point) in order to reproduce the experimental data. The study of the sensitivity of the seismic data to the source location was also crucial to improve the amplitude of the signals and the detection of the different seismic modes. This will allow us in the future to better image and understand these complex media