Je Tiens À Remercier Madame

Monsieur Soumaya Latour

Ludovic Bodet D'avoir Accepté

Artiom, Gaétan, Moussa, Clément, Rémi, Javier, Salia, Juraj Mohamed Steven

Ali Zaynab

Hassan Kassem

Bilal, Hussein, Hadi, Khaled, Abbas, Mohamed Zamat

Hasan André

II

Keywords: Near surface geophysics, Wave propagation, numerical modeling, finite volume method, Riemann solver, Optimized PML boundary conditions, dispersion analysis, surface waves, porous and granular media

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 interest 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 unconsolidated 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 classical 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.

Résumé

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

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és 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

Chapter 1 Introduction (English)

The Earth has always aroused the curiosity and interest of humans and scientists in particular. Indeed, the subsurface, which contains many natural resources (water, gas, etc.), can also constitute a natural risk because of its lithological and topographical characteristics. These characteristics conjugated to site/topographic effects can be at the origin of an accentuation of the destructive effects of an earthquake. Besides, another interest to define the properties of the near surface is to estimate the rate of saturation of fluids in complex media and monitor in time and space those fluid contents at different scales. This is of particular importance in gas and fluid recovery (hydrocarbon or gas industry), or, more important, to estimate the water resources in the context of climate change. All these reasons are at the center of the concerns of the community of geophysicists who seek to better understand the near surface and therefore to characterize it.

In this thesis, I use seismic waves as a physical phenomenon allowing to determine the physical parameters of the subsurface. Seismic waves have natural sources, for example earthquakes that take place every day more or less deeply in the lithosphere, as well as natural surface sources (interactions with the hydrosphere or atmosphere) or anthropogenic sources (intended or not). The mechanical waves which propagate inside the globe, from the source to the receivers (accelerometers or velocimeters, located on the surface or in borehole), diffract, reflect, convert, and permit one to characterize the medium in which they have traveled. By analyzing, processing and/or inverting the seismograms, we can therefore reconstruct the internal structure of the medium crossed by the waves and quantify the constitutive parameters of this medium. However, natural media are complex and heterogeneous, and their interpretation is often difficult and incomplete. Often we explain only the arrival times of the first waves (refraction, tomography), measure the amplitude of certain reflected waves (AVO analyzes) or invert the seismic modes associated to the surface waves that are present in the diagrams of dispersion.

To better understand acoustic and seismic measurements in sediments and soils, many studies on unconsolidated granular media have been conducted in situ, and on analogue experiments at the laboratory scale where theoretical models have been developed. In my thesis, I 1 will focus on granular media which are a type of complex medium difficult to characterize.

A granular medium is an assembly of discrete particles in contact with a liquid or a gas [START_REF] Castellanos | The relationship between attractive interparticle forces and bulk behaviour in dry and uncharged fine powders[END_REF] and whose size is very variable: fine powders (∼ 10 -6 m), desert dunes (∼ 10 -4 m), ballast railways (∼ 10 -1 m)... Many industries deal with grains in agriculture, civil engineering and pharmaceutical production manufacturing. They also play an important role in geotechnical processes such as landslides or erosion and, on a greater scale, the seismic plate tectonics which determines the morphology of the Earth. Despite their apparent simplicity, granular media are difficult to classify among the three usual states of matter, namely solid, liquid and gaseous and are sometimes regarded as an additional state of matter [START_REF] Jaeger | Physics of the granular state[END_REF]. Depending on its properties and the external conditions, a granular medium can behave as a solid, a liquid or a gas [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. Moreover, a granular medium has properties that are not found in any of the other three states, among which we can cite the vault effect, dilatancy, segregation and an overall elasticity depending on the pressure applied to the inhomogeneous network of inter-grain contacts [START_REF] Dantu | Contribution à l'étude mécanique et géométrique des milieux pulvérulents[END_REF][START_REF] Makse | Packing of compressible granular materials[END_REF][START_REF] Radjaï | Contact dynamics study of 2d granular media: critical states and relevant internal variables[END_REF]. We can distinguish two large groups: the granular cohesive and non-cohesive stacks. They are cohesive when the attractive force is much greater than the weight of the particle and non-cohesive when the forces between grains are only repulsive and the shape of the material is determined by the conditions at the limits (applied stress or gravity). Understanding the effects of cohesion and developing methods to characterize this type of medium are attracting more and more interests. As fine powders are very complex to characterize because of the dependence of their macroscopic state to many variables such as inter-grain adhesion, anelasticity and surface irregularities, geophysicists choose to study ideal granular materials which have no cohesion and have a simple geometry (smooth spherical balls) [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF][START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][START_REF] Tournat | Acoustics of unconsolidated "model" granular media: An overview of recent results and several open problems[END_REF]].

The granular media are complex media, because of their heterogeneity associated with the weak consolidation of the grains. They have voids and the related pores are filled with different fluids. The presence of fluids modifies the overall mechanical characteristics, via the type of fluid (which has its own mechanical characteristics), its quantity, its spatial distribution (depending on porosity) and its flow capacity (related to permeability). The response of these media to a mechanical excitation, at different scales, can be predicted by seismic models such as V p and V s power-law models depending on the depth of the medium [START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, part ii: Velocities in shallow unconsolidated sand[END_REF][START_REF] Bachrach | Seismic velocities and poisson's ratio of shallow unconsolidated sands[END_REF][START_REF] Gassmann | Elastic waves through a packing of spheres[END_REF]. In this thesis, in order to model this media, we consider on one hand an elastic rheology (simple) and on the other hand we add complexity to the rheology by considering a poro-elastic model inspired from the Biot theory [Biot, 1956a,b]. The need to simulate correctly the propagation of seismic waves in complex media (elastic and poroelastic) are then necessary. However, we only know analytical or semi-analytical solutions of wave propagation problems in extremely simple cases: homogeneous models or models formed by plane layers. Since the 1970s, sophisticated numerical methods have emerged with the spectacular increase in the capacities (computing power and memory size) of computers. Some of these solve the elastodynamics equations in the frequency domain and other calculate the entire wave field in the time domain for complex 2D and 3D structures. In a non-exhaustive way, we can cite the finite differences, the finite elements, the spectral elements, the finite volumes, the discontinuous Galerkin method, ...

Scientific problem and objectives

The main objective of this thesis is to model unconsolidated granular media. To achieve this objective, several sub-objectives have been defined as follows:

• Develop a numerical tool which calculates the entire wave field of a two-(or three-)dimensional elastic model with complex structures. In three-dimensional configurations, a parallelized version of the codes must be used.

• Develop signal processing tools to calculate phase velocity curves and to detect first arrival travel times and wave propagation modes of the experimental and simulated data.

Several challenges should be addressed when implementing a code that solves the elastodynamics equations in the time domain for purely elastic or poroelastic media (with or without attenuation):

• The implementation of absorbing boundary conditions for the discretized domain to simulate an open domain. An exhaustive review of the main techniques used is detailed in [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]] and [ [START_REF] Meza-Fajardo | A nonconvolutional, split-field, perfectly matched layer for wave propagation in isotropic and anisotropic elastic media; stability analysis[END_REF]]. In the context of this thesis, convolutional perfectly matched layer (C-PML) absorbing conditions have been adapted to the finite volume method and implemented [START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF]. This type of absorbing condition constitutes an optimized version of the classical PMLs initially developed for electromagnetic wave propagation applications [START_REF] Bérenger | A Perfectly Matched Layer for the absorption of electromagnetic waves[END_REF].

• Taking into account free surface conditions (air-ground). The precise modeling of wave fields close to the surface is an imperative for subsurface geophysical applications (due to the effect of Rayleigh waves). In order to adapt and validate the free surface conditions, we used an antisymetric condition combined with a high order scheme to calculate variables close to the free surface. The antisymetric condition is also called the method of images using by [START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF] for 'standard staggered grid' (SSG) stencils, whose principle consists in canceling the vertical components of the field of constraint and compute the other components at the level of the boundary [Levander, 1988a;[START_REF] Robertsson | A numerical free-surface condition for elastic/viscoelastic finitedifference modeling in the presence of topography[END_REF]. We achieve more accuracy by implementing a non-centered fourth-order spatial scheme close to the free surface proposed by [START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF]. Besides, the Dirichlet boundary conditions (rigid boundary) consist in setting the whole velocity components to zero at the boundary.

• The implementation of a finite volume method and its comparison with classical staggered grid finite differences and finite elements like the spectral element method. Indeed, we propose here 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 the frequency shift convolutional (C-PML) or the non convolutional formulations (ADE-PML). But here, they are applied to the spatial fluxes derivatives, which is a different formulation than classical PMLs that are generally applied to the spatial derivatives of the primitive variables (particle velocities and stresses). They can be even further improved by introducing sponge layers-like or paraxial absorbing conditions. The finite volume is able to compute the physical variables (particle velocities and stresses) and take into account the physical properties of the medium (density and seismic velocities) at a same grid point, which allows more flexibility and facilities to treat numerically the free surface and solid-solid or fluid-solid interfaces. Furthermore, this FVM version is able to introduce coarsening of the grids more easily when increasing seismic velocities with depth (as it is usually the case in many geophysical applications): in this case mesh cells could be coarsened by simple interpolations at fine-to-coarse grid interfaces located at different depths.

Beyond this thesis, it will also be much easier with this RFV, when compared to staggered finite differences, to introduce distorted meshes, topographies, bathymetry, and distorted interfaces between geological layers. This RFV method can be an alternative to staggered grid finite differences (SSG-FD) and a compromise between SSG-FD and high-order finite element methods like the spectral finite element method (SPECFEM).

Applications

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.

The validation of all the points mentioned in the previous section was done through theoretical and experimental cases. Indeed, more specifically, we consider some synthetic models composed of two materials (checkerboard-like), and others composed of three materials among which water (three-layer model). In these cases, we place the source and the receivers sometimes in the middle of the medium to check if solutions at the interface are well calculated, and sometimes close to the free surface to check if we correctly calculate the surface waves. More particularly 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 as in [START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF].

In addition to the synthetic cases considered, and given the difficulty to characterize a model at the scale of a natural medium, the validation of the numerical modeling code was carried out on a more realistic model calculated in the context of unconsolidated granular media. Indeed, the use of perfectly controlled physical models is better suited. Ultrasonic techniques (measurements by laser interferometry) were used because they allow non-contact measurements and high sampling capacity, making it possible to simulate classical seismic recordings, on a centimeter scale and with excellent reproducibility. This approach has been shown to be successful for the analysis of surface waves [START_REF] Bodet | Near-offsets effects on Rayleigh-wave dispersion measurements: Physical modelling[END_REF][START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF][START_REF] Bodet | Laser-Doppler acoustic probing of granular media with in-depth property gradient and varying pore pressures[END_REF][START_REF] Bodet | Surface-wave inversion limitations from laser-Doppler physical modeling[END_REF]. The experimental data recorded on unconsolidated granular media are compared qualitatively and quantitatively to synthetic seismograms. This work was carried out in collaboration with researchers of the laboratory METIS/Paris Sorbonne (CNRS, Pierre et Marie Curie University) (Laboratoire des Milieux environnementaux, transferts et interactions dans les hydrosystèmes et les sols) who provided us with the experimental data obtained from a laboratory experiment performed on unconsolidated granular media. This experience with more realistic models open the way for other experiments, namely the injection of air and/or fluid and the introduction of a non-horizontal topography.

Plan of the thesis

The Chapter 1 is dedicated to theoretical reminders of wave equations in an elastic (3D and 2D) and poro-elastic medium (2D). The different types of waves, the source implementation, the different boundary conditions, and the time/dispersion analysis tools are also developed.

In Chapter 2, we present a state of the art of seismic modeling by making an inventory of the different numerical methods used in recent years to model wave propagation. We also develop the theoretical aspects of the finite volume method used for the development of the elastic code. We describe a finite difference method (to compare with the finite volume method), particularly the Standard Staggered Grid (SSG) discretization schemes, and also the spectral finite elements method used in SPECFEM (taken as reference solution) [START_REF] Casarotti | Cubit and seismic wave propagation based upon the spectral-element method: An advanced unstructured mesher for complex 3d geological media[END_REF][START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Spectral and spectral-element methods for the 2D and 3D elastodynamics equations in heterogeneous media)[END_REF][START_REF] Komatitsch | Advanced numerical modeling in geophysics (Modélisation numérique avancée pour la géophysique)[END_REF][START_REF] Komatitsch | Comment on "Multidomain Pseudospectral Time-Domain (PSTD) method for acoustic waves in lossy media[END_REF][START_REF] Rietmann | Forward and adjoint simulations of seismic wave propagation on emerging large-scale gpu architectures[END_REF][START_REF] Wang | Numerical modeling of seismic wave propagation and ground motion in underground mines[END_REF]. After developing various points related to absorbing conditions and the implementation of the free surface, we conclude this chapter with a study in which we consider three main models: a checkerboardlike heterogeneous medium (to validate the absorbing and free surface conditions), a fluid-solid model with water at the top of the medium and two solid layers below (to test the stability of the PMLs conditions for acoustic-elastic coupling simulations) and a realistic model (power-law model) determined in the context of unconsolidated granular media.

In Chapter 3, we revisit a study carried out in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] on unconsolidated granular media at the laboratory scale. We consider different models with different rheologies (elastic or poro-elastic), different dimensions (3D or 2D), various boundary conditions (PML or Dirichlet) and different numerical modeling of the source (vibrating stick or point source) in order to reproduce the experimental data at the laboratory scale. The localization of the source was important and crucial, and it allowed us to obtain signals much more similar to the measured signals in terms of amplitude in both time and frequency (dispersion diagrams) domains. The 3D and 2D elastic models are respectively simulated using UNISOLVER and Seismic_CPML [Martin et al., 2008c,d], two finite differences codes which use staggered grids. The poroelastic models simulations are done with a version of Seismic_CPML for porous media. Then, a qualitative and quantitative study was carried out by comparing the simulated data with the experimental data. In this context, the signal processing tools implemented numerically during this thesis are used in order to make a time and dispersion analysis.

Introduction (Français)

La structure de la Terre a toujours suscité la curiosité et l'intérêt des hommes et des scientifiques en particulier. En effet, le sous-sol, qui contient de nombreuses ressources naturelles (eau, gaz, etc.), peut également constituer un risque naturel en raison de ses caractéristiques lithologiques et topographiques. Ces caractéristiques conjuguées aux effets de site/topographie peuvent être à l'origine d'une accentuation des effets destructeurs d'un séisme. Par ailleurs, un autre intérêt à définir les propriétés de la proche surface est d'estimer le taux de saturation des fluides dans des milieux complexes et de suivre dans le temps et l'espace les contenus de ces fluides à différentes échelles. Ceci est particulièrement important pour la récupération de ressources en gaz et de fluides (industrie des hydrocarbures ou du gaz), ou, plus important encore, pour estimer les ressources en eau dans le contexte du changement climatique. Toutes ces raisons sont au centre des préoccupations de la communauté des géophysiciens qui cherchent à mieux comprendre la proche surface et donc à la caractériser.

Dans cette thèse, j'utilise les ondes sismiques comme un phénomène physique permettant de déterminer les paramètres physiques de la subsurface. Les ondes sismiques sont générées par des sources naturelles, comme par exemple les tremblements de terre qui ont lieu chaque jour plus ou moins profondément dans la lithosphère, ainsi que par des sources naturelles en surface (interactions avec l'hydrosphère ou l'atmosphère) ou anthropiques (intentionnelles ou non). Les ondes mécaniques qui se propagent à l'intérieur du globe, de la source aux récepteurs (accéléromètres ou vélocimètres, situés en surface ou dans des puits de forage), se diffractent, se réfléchissent, se transforment, et permettent de caractériser le milieu dans lequel elles ont voyagé. En analysant, traitant et/ou inversant les sismogrammes, on peut donc reconstruire la structure interne du milieu traversé par les ondes et quantifier les paramètres constitutifs de ce milieu. Cependant, les milieux naturels sont complexes et hétérogènes, et leur interprétation est souvent difficile et incomplète. Souvent on n'explique que les temps d'arrivée des premières ondes (réfraction, tomographie), on mesure l'amplitude de certaines ondes réfléchies (analyses AVO) ou encore on inverse les modes sismiques associés aux ondes de surface présentes dans les diagrammes de dispersion.

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 sur des expériences analogiques à l'échelle du laboratoire où des modèles théoriques ont été développés. Dans ma thèse, je me concentrerai sur les milieux granulaires qui sont un type de milieu complexe difficile à caractériser.

Un milieu granulaire est un assemblage de particules discrètes en contact avec un liquide ou un gaz [START_REF] Castellanos | The relationship between attractive interparticle forces and bulk behaviour in dry and uncharged fine powders[END_REF] et dont la taille est très variable : poudres fines (∼ 10 -6 m), dunes du désert (∼ 10 -4 m), ballast ferroviaire (∼ 10 -1 m)... De nombreuses industries traitent les grains dans l'agriculture, le génie civil, la fabrication de produits pharmaceutiques, etc. Ils jouent également un rôle important dans les processus géotechniques tels que les glissements de terrain ou l'érosion et, à plus grande échelle, dans la tectonique des plaques sismiques qui détermine la morphologie (lithosphère en partiulier) de la Terre. Malgré leur apparente simplicité, les milieux granulaires sont difficiles à classer parmi les trois états habituels de la matière, à savoir solide, liquide et gazeux. Ils sont parfois considérés comme un état supplémentaire de la matière [START_REF] Jaeger | Physics of the granular state[END_REF]]. En fonction de ses propriétés et des conditions extérieures, un milieu granulaire peut se comporter comme un solide, un liquide ou un gaz [START_REF] Jaeger | Granular solids, liquids, and gases[END_REF]. De plus, un milieu granulaire possède des propriétés qui ne se retrouvent dans aucun des trois autres états, parmi lesquelles on peut citer l'effet de voûte, la dilatance, la ségrégation et une élasticité globale dépendant de la pression appliquée au réseau inhomogène de contacts intergrains [START_REF] Dantu | Contribution à l'étude mécanique et géométrique des milieux pulvérulents[END_REF][START_REF] Makse | Packing of compressible granular materials[END_REF][START_REF] Radjaï | Contact dynamics study of 2d granular media: critical states and relevant internal variables[END_REF]]. On peut distinguer deux grands groupes : les empilements granulaires cohésifs et non cohésifs. Ils sont cohésifs lorsque la force d'attraction est beaucoup plus grande que le poids de la particule et non cohésifs lorsque les forces entre les grains sont uniquement répulsives et que la forme du matériau est déterminée par les conditions aux limites (contrainte appliquée ou gravité). La compréhension des effets de cohésion et le développement de méthodes pour caractériser ce type de milieu suscitent de plus en plus d'intérêt. Comme les poudres fines sont très complexes à caractériser en raison de la dépendance de leur état macroscopique à de nombreuses variables telles que l'adhésion entre les grains, l'inélasticité et les irrégularités de surface, les géophysiciens choisissent d'étudier des matériaux granulaires idéalisés qui n'ont aucune cohésion et ont une géométrie simple (boules sphériques ou ellipsoidales lisses) [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF][START_REF] Jaeger | Physics of the granular state[END_REF][START_REF] Jaeger | Granular solids, liquids, and gases[END_REF][START_REF] Tournat | Acoustics of unconsolidated "model" granular media: An overview of recent results and several open problems[END_REF].

Les milieux granulaires sont des milieux complexes, en raison de leur hétérogénéité liée à la faible consolidation des grains. Ils présentent des vides et ces pores peuvent être remplis de différents fluides. La présence de fluides modifie les caractéristiques mécaniques globales, par l'intermédiaire du type de fluide (qui a ses propres caractéristiques mécaniques), sa quantité et sa distribution spatiale (dépendant de la porosité) et sa capacité d'écoulement (liée à la perméabilité). La réponse de ces milieux à une excitation mécanique, à différentes échelles, peut être prédite par des modèles sismiques tels qu'un modèle de loi de puissance V p et V s en fonction de la profondeur du milieu [START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, part ii: Velocities in shallow unconsolidated sand[END_REF][START_REF] Bachrach | Seismic velocities and poisson's ratio of shallow unconsolidated sands[END_REF][START_REF] Gassmann | Elastic waves through a packing of spheres[END_REF]. Dans cette thèse, afin de modéliser ce milieu, nous considérons d'une part une rhéologie élastique (simple) et d'autre part nous complexifions la rhéologie en considérant un modèle poro-élastique inspiré de la théorie de Biot [Biot, 1956a,b]. La nécessité de simuler correctement la propagation des ondes sismiques dans des milieux complexes (élastiques et poro-élastiques) s'impose alors. Cependant, nous ne connaissons des solutions analytiques ou semi-analytiques des problèmes de propagation des ondes que dans des cas extrêmement simples : modèles homogènes ou formés de couches planes. Depuis les années 1970, des méthodes numériques sophistiquées sont apparues avec l'augmentation spectaculaire des capacités (puissance de calcul et taille mémoire) des ordinateurs. Certaines d'entre elles permettent de résoudre les équations de l'élastodynamique dans le domaine fréquentiel et d'autres dans le domaine temporel pour des structures complexes en 2D et 3D. De manière non exhaustive, on peut citer les différences finies, les éléments finis, les éléments spectraux, les volumes finis, la méthode de Galerkin discontinue, ...

Problème scientifique et objectifs

L'objectif principal de cette thèse est de modéliser les milieux granulaires non consolidés. Pour atteindre cet objectif, plusieurs sous-objectifs ont été définis comme suit :

• Développer un outil numérique dédié au calcul de l'ensemble du champ d'onde d'un modèle élastique géométrique bidimensionnel (ou tridimensionnel) avec des structures complexes. Dans les configurations tridimensionnelles, une version parallélisée des codes devra être utilisée.

• Développer des outils de traitement du signal pour calculer les courbes de vitesse de phase et pour détecter les temps de parcours des premières arrivées et les modes de propagation des ondes des données expérimentales et simulées.

Plusieurs défis peuvent être rencontrés lors de l'implémentation d'un code qui résout les équations de l'élastodynamique dans le domaine temporel pour un milieu purement élastique ou poroélastique (avec ou sans atténuation) :

• L'implémentation de conditions aux limites absorbantes pour le domaine discrétisé afin de simuler un domaine ouvert. Une revue exhaustive des principales techniques utilisées est détaillée dans [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF] et [START_REF] Meza-Fajardo | A nonconvolutional, split-field, perfectly matched layer for wave propagation in isotropic and anisotropic elastic media; stability analysis[END_REF]. Dans le cadre de cette thèse, les conditions absorbantes de type PML qui utilisent des formulations convolutives avec décalage en fréquence (C-PML) ont été adaptées à la méthode des volumes finis et implémentées [START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF]. Ce type de conditions absorbantes constitue une version optimisée des PML classiques initialement développées pour des applications de propagation d'ondes électromagnétiques [START_REF] Bérenger | A Perfectly Matched Layer for the absorption of electromagnetic waves[END_REF].

• La prise en compte des conditions de surface libre (air-sol). La modélisation précise des champs d'ondes proches de la surface est un impératif pour les applications géophysiques de subsurface (en raison de l'effet des ondes de Rayleigh). Afin d'adapter et de valider les conditions de surface libre, nous avons utilisé une condition d'antisymétrie combinée à un schéma d'ordre élevé pour calculer les variables proches de la surface libre. La condition d'antisymétrie est aussi appelée méthode des images utilisée par [START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF] Au-delà de cette thèse, il sera également beaucoup plus facile avec cette RFV, par rapport aux différences finies en quinconce, d'introduire des maillages déformés, des topographies, de la bathymétrie, et des interfaces déformées entre couches géologiques. Cette méthode RFV peut être une alternative aux différences finies à grille en quinconce (SSG-FD) et un compromis entre les SSG-FD et des méthodes d'éléments finis d'ordre élevé telles que la méthode des éléments finis spectraux (SPECFEM).

Applications

La méthode des volumes finis et les différents types de conditions aux limites sont testés et validés sur différents cas synthétiques hétérogènes. La méthode des volumes finis est comparée à d'autres techniques comme les différences finies et les éléments finis d'ordre élevé. Enfin, nous appliquons notre méthode à une configuration de couplage fluide-solide et à certains 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.

La validation de tous les points mentionnés dans la section précédente a été faite à travers des cas théoriques et expérimentaux. En effet, plus précisément, nous considérons certains modèles synthétiques composés de deux matériaux (modèle en damier), et d'autres composés de trois matériaux dont l'eau (modèle à trois couches). Dans ces cas, nous plaçons la source et les récepteurs parfois au centre du milieu étudié pour vérifier si les solutions à l'interface sont bien calculées, et parfois près de la surface libre pour vérifier si nous calculons correctement les ondes de surface. Plus particulièrement, nous portons notre attention sur l'implémentation de la résolution numérique des ondes de surface comme les ondes de Rayleigh, ce qui n'est pas trivial avec les différences finies en quinconce classiques. Nous avons donc implémenté un schéma spatial d'ordre 4 non centré à la surface libre pour atteindre une plus grande précision comme dans [START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF].

En plus des cas synthétiques considérés, et étant donné la difficulté de caractériser un modèle à l'échelle d'un milieu naturel, la validation du code de modélisation numérique a été effectuée sur un modèle plus réaliste calculé dans le contexte des milieux granulaires non consolidés. En effet, l'utilisation de modèles physiques pour des milieux parfaitement contrôlés est mieux adaptée. Les techniques ultrasonores (mesures par interférométrie laser) ont été utilisées car elles permettent des mesures sans contact et une grande capacité d'échantillonnage, permettant de simuler des enregistrements sismiques classiques, à l'échelle centimétrique et avec une excellente reproductibilité. Cette approche s'est avérée fructueuse pour l'analyse des ondes de surface [START_REF] Bodet | Near-offsets effects on Rayleigh-wave dispersion measurements: Physical modelling[END_REF][START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF][START_REF] Bodet | Laser-Doppler acoustic probing of granular media with in-depth property gradient and varying pore pressures[END_REF][START_REF] Bodet | Surface-wave inversion limitations from laser-Doppler physical modeling[END_REF] [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] sur les milieux granulaires non consolidés à l'échelle du laboratoire. Nous considérons différents modèles avec différentes rhéologies (élastique ou poro-élastique), différentes dimensions (3D ou 2D), différentes conditions aux limites (PML ou Dirichlet) et différentes modélisations numériques de la source (pot vibrant ou source ponctuelle) afin de reproduire les données expérimentales à l'échelle du laboratoire. La localisation de la source était importante et cruciale, elle nous a permis d'obtenir des signaux beaucoup plus similaires aux signaux mesurés en termes d'amplitude dans les domaines temporel et fréquentiel (diagrammes de dispersion). Les modèles élastiques 3D et 2D sont respectivement simulés à l'aide de UNISOLVER et Seismic_CPML [Martin et al., 2008c,d], deux codes de différences finies qui utilisent des grilles en quinconce. Les simulations des modèles poroélastiques sont réalisées avec une version de Seismic_CPML pour les milieux poreux. Ensuite, une étude qualitative et quantitative a été réalisée en comparant les données simulées avec les données expérimentales. Dans ce contexte, les outils de traitement du signal implémentés numériquement au cours de cette thèse sont utilisés afin d'effectuer une analyse temporelle et de dispersion.

Chapter 2

Seismic wave propagation: elastic and poro-elastic equations and data analysis tools

Introduction

In this chapter we describe theoretical aspects related to the wave propagation in elastic or poroelastic media. We will develop the equations of elastodynamics from the principles of continuum mechanics (CM). We will present the different types of waves, sources and boundary conditions. Finally, to model unconsolidated granular media, we use qualitative and quantitative analysis tools that we will present here.

Elastodynamics equations

The equations of elastodynamics are derived from the principles of CM. They ensure conservation of momentum and they link stresses, strains and displacements (or velocities). In this section, we present conservation equations and Hooke's laws, the three dimensional (3D) and two dimensional (2D) elastic wave equations for a linear isotropic medium. In order to model unconsolidated granular media partially saturated with fluid, a poroelastic model in two dimensions (2D) is also considered. We will present more particularly the Biot model without attenuation effects (see [Biot, 1956a], [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: Higher-frequency range[END_REF]). For further details, we invite the reader to consult the following books: [Aki and Richards, 2002], [START_REF] Spencer | Continuum mechanics[END_REF], [START_REF] Pride | Seismic attenuation due to waveinduced flow[END_REF], [START_REF] Semblat | Waves and vibrations in soils[END_REF].

Conservation laws: conservation of mass and momentum

Consider a body with a mass M (t) of volume V at an instant t in a fixed frame of reference in space. The law of mass conservation for this body is written as:

13 dM (t) dt = 0 (2.1)
Let u be the displacement of the body M from position x 0 at time t 0 to a position x at the instant t. The velocity vector of mass displacement is then written in the form ∂u ∂t and the acceleration vector ∂ 2 u ∂t 2 . Applying the fundamental principle of dynamics, we obtain the equation of conservation of momentum:

V ρ ∂ 2 u ∂t 2 dV = V F dV + S T dS (2.2)
This equation relates the momentum contained in the volume V and the sum of two types of forces:

• external forces F applied to the volume V,

• internal forces T resulting from the traction exerted on the surface S of the body M.

In equation (2.2), ρ represents the density in Kg/m 3 , F corresponds to the volumetric mass density of external forces in N/m 3 and T the surface force in N/m 2 .

Stress tensor

The i -th component of the surface forces T (tractions) is:

T i = σ ij n j (2.3)
where σ ij is the stress tensor and n j the j -th component of the vector normal to S. The expression (2.3) implies an implicit summation over the indices j (Einstein convention). Using the divergence theorem, the surface integral of equation (2.2) is transformed into a volume integral:

S T i dS = V ∂ j σ ij dV (2.4)
where ∂ j is the partial derivative according to the j -th direction. Thus from expressions (2.2) and (2.4), we obtain the equation for the conservation of momentum:

ρ ∂ 2 u i ∂t 2 = F i + ∂ j σ ij (2.5)

Strain tensor

In general, the relation between strain tensor (ϵ ij ) and displacements is written as follows:

ϵ ij = 1 2 ∂u i ∂x j + ∂u j ∂x i + ∂ 2 u k ∂x i ∂x j (2.6)
In the elastic case (small deformations), the hypothesis of small perturbations is applied in the last part of equation (2.6). We then get the strain tensor for infinitesimal strains:

ϵ ij = 1 2 ∂u i ∂x j + ∂u j ∂x i = 1 2 (u i,j + u j,i ) (2.7)

Hooke's law: strain-stress relation

The rheological law of the medium is formalized by the strain-stress relation for elastic, visco-elastic, fluid media... Its generalized expression is written in a Taylor expansion form of the stress tensor:

σ ij = σ 0 ij + C ijklmn ϵ kl + D ijklmn ϵ kl ϵ mn + O(ϵ 3 ) (2.8)
where σ 0 ij represents the initial stress state of the system when it is at equilibrium (before perturbation), C ijklmn and D ijklmn are tensors of order 4 and 6 respectively. O(ϵ 3 ) defines the order of the Taylor expansion given in equation (2.8).

In the case of an elastic medium, we assume that the stress state of the system at the equilibrium is zero, and the generalized Hooke's law is used and represented by the tensor:

C klij = C ijlk = C ijkl = C jikl (2.9)
Hence the classical stress-strain relationship in elasticity is:

σ ij = C ijklmn ϵ kl (2.10)
For an isotropic elastic linear medium (properties of the medium identical according to all directions), the elastic tensor is written as:

C ijklmn = λδ ij δ kl + µδ ik δ jl + µδ il δ jk (2.11)
where λ and µ are constants called the Lamé parameters and δ ij is the Kronecker symbol defined as follows: δ ij = 1 if i = j and δ ij = 0 if i ̸ = j otherwise. By replacing expression (2.11) in equation (2.10), the strain tensor is written as:

σ ij = λδ ij ϵ kk + 2µϵ ij (2.12)
λ and µ are related to physical parameters (measurable). Other parameters such as Young's modulus E, the bulk modulus K and the Poisson's ratio v can be expressed using Lamé parameters:

• The Bulk modulus K = λ + 2 3 µ • Young's modulus E = µ(3λ+2µ) λ+µ • the Poisson's ratio ν = λ 2(λ+µ)
• The compressional wave or P -wave velocity V p = ( λ+2µ ρ )

• The shear wave or S-wave velocity V s = ( µ ρ ) Some useful properties of these parameters can be mentioned. The Poisson's ratio is such that:

-1 < ν ≤ 1 2
(2.13)

The limit ν = 1 2 corresponds to a very weakly consolidated medium. For most rocks, the coefficient of Poisson is between 0.20 and 0.35. P and S waves velocities are such that:

V s < V p √ 2 (2.14)
Hence the name of these waves (P = Primary, S = Secondary), the direct P waves always arriving first, before the S waves. We also attribute the meanings P = Pressure and S = Shear to these waves, the P waves being waves of compression and the S waves shear waves.

3D elastic wave equations

In summary, by taking the stress formulation given previously, the elastodynamics equations can be formulated at the second order in displacement as:

ρ ∂ 2 u i ∂t 2 = ∂ j σ ij + s i , ϵ ij = 1 2 (u i,j + u j,i ) , σ ij = λδ ij ϵ kk + 2µϵ ij ,
However, the first-order formulation (velocity-stress) of the 3D elastic wave equations for a linear, isotropic medium submitted to external forces are given by [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggeredgrid finite differences[END_REF]

2D elastic wave equations

In the 2D particular case, 2D elastic wave equations for an isotropic medium submitted to external forces can be written using a velocity-stress formulation such as the following linear and hyperbolic system [START_REF] Dumbser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-II. The three-dimensional isotropic case[END_REF] where λ and µ are Lamé parameters, ρ is the density, and s x and s z are the space dependent source terms in x and z directions. The compressional stress components are given by σ xx and σ zz , and the shear stress is σ xz . The components of particle velocities in direction x and z are denoted by v x and v z , respectively.

In a 3D homogeneous medium, the geometric attenuation is in 1 r while in 2D (as is the case here) it decreases in 1 r . r being the source-observation position distance. This is a property that will be used in Chapter 3 to reduce time consumption in our calculations and our signal analysis by rescaling our seismograms from 3D to 2D.

2D poro-elastic wave equations

Porous materials are made of a solid phase (called the frame) and of a fluid phase, and can be considered as an interconnected network of pores inside the solid [START_REF] Pride | Seismic attenuation due to waveinduced flow[END_REF]. When a fluid flow is able to cause the solid to deform, the material is called poroelastic. Unconsolidated granular media can be seen as a poroelastic material in which air or water can play the role of the fluid and, grains the solid. Poroelastic materials are most of the time modelled using the Biot theory [Biot, 1956a,b]. The differential or " strong" formulation of the poroelastic wave equations can be written as [Carcione, 2007a[START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF]:

ρ∂ 2 t u s + ρ f ∂ 2 t w = ∇. C : ∇u s -αP f I ρ f ∂ 2 t u s + ρ w ∂ 2 t w = -∇P f -K∂ t w P f = -αM ∇.u s -M ∇.w (2.17) (2.18) (2.19)
where

u s = (u s i ) i=1,2 , w = Φ u f -u s and u f = (u f i ) i=1,2
are respectively the solid, relative, and fluid displacement vectors; Φ is the porosity; and C is the stiffness tensor of the isotropic elastic solid matrix, defined as:

σ s ij = (C : ϵ) ij = λ s δ ij ϵ kk + 2µϵ ij ϵ ij = 1 2 ∂u s i ∂x j + ∂u s j ∂x i P f = -αM ∇.u s -M ∇.w (2.20) (2.21) (2.22)
where indices i and j can be here 1 or 2 in 2D and with the Einstein convention of implicit summation over a repeated index. P f is the pressure in the fluid. σ s and ϵ are respectively the stress and strain tensors of the isotropic elastic solid frame. The stress tensor of the fluid-filled solid matrix is σ = σ s -αP f I, and ρ = Φρ f +(1 -Φ) ρ s is the density of the saturated medium, where ρ s and ρ f are the solid and fluid densities, respectively. The apparent density is

ρ w = a ρ f Φ
with a the tortuosity. The shear modulus is µ , and λ s = λ -α 2 M is the Lamé coefficient in the solid matrix, where λ is the Lamé coefficient of the saturated matrix. The α and M variables are functions of the porosity and bulk moduli of the fluid and solid components of the porous medium and are given by the following expressions:

α = 1 - K f r K s M = 1/ [Φ/K f + (Φ -α)/K s ] (2.23) (2.24)
where K f r is the incompressibility modulus of the porous frame, K s is the incompressibility modulus of the solid matrix, and K f is the incompressibility modulus of the fluid. The viscous damping coefficient is:

K = κ/η (2.25)
where κ is the permeability of the solid matrix and η is the fluid viscosity. Equations (2.17) to (2.20) can be written using a first order velocity-stress formulation:

ρ w ρ -ρ 2 f ∂ t v s = ρ w ∇.σ + ρ f ∇P f + ρKv f ρ w ρ -ρ 2 f ∂ t v f = -ρ f ∇.σ -ρ∇P f -ρ f Kv f ∂ t σ = C : ∇v s -α∂ t P f I ∂ t P f = -αM ∇.v s -M ∇.v f (2.26) (2.27) (2.28) (2.29)
where

v s = (v s i ) i=1,2 and v f = ∂ t w = (v f i ) i=1,2
are the solid and filtration velocity vectors, respectively. σ is the effective stress tensor of the porous medium. As in [START_REF] Zeng | A staggered-grid finite-difference method with perfectly matched layers for poroelastic wave equations[END_REF], using the trace of the strain tensor T r(ϵ) = ϵ ii and an auxilliary variable ξ, we rewrite the system as:

ρ w ρ -ρ 2 f ∂ t v s i = ρ w ∂ j σ ij + ρ f ∂ i P f + ρKv f i ρ w ρ -ρ 2 f ∂ t v f i = -ρ f ∂ j σ ij -ρ∂ i P f -ρ f Kv f i ϵ ij = 1 2 ∂ j v s i + ∂ i v s j ∂ t ξ = -∂ i v f i P f = -αM T r(ϵ) + M ξ σ s ij = λ s δ ij T r(ϵ) + 2µϵ ij σ ij = σ s ij -αP f δ ij (2.30) (2.31) (2.32) (2.33) (2.34) (2.35) (2.36)
This system of equations has seven wave eigenvalues related to seven wave velocity modes (instead of five for the elastic case). Those wave velocities are ±V pF AST , ±V pSLOW , ±V s and 0. The fast and slow wave velocities V pF AST and V pSLOW can be expressed as [START_REF] Sidler | A pseudospectral method for the simulation of 3-D ultrasonic and seismic waves in heterogeneous poroelastic borehole environments[END_REF] :

V pF AST = -b 1 + √ ∆ 2a 1
(2.37)

V pSLOW = -b 1 - √ ∆ 2a 1 (2.38)
where

a 1 = ρ 11 ρ 22 -ρ 2 12 b 1 = -Sρ 22 -Rρ 11 + 2gaρ 11 c 1 = SR -ga 2 ρ 11 = ρ + ρ f Φ(a -2)
ρ 12 = Φρ f (1 -a) ρ 22 = aρ f Φ S = λ + 2µ R = M Φ 2 ∆ = b 2 1 -4a 1 c 1 ga = M Φ(α -Φ)
In table 2.1, a summary of the different parameters of the poroelastic model is provided. Biot's characteristic frequency f c defines the transition between two poroelastic regimes (with or without attenuation) and is given by [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: Higher-frequency range[END_REF]; [Carcione, 2007a]; [START_REF] Morency | Spectral-element simulations of wave propagation in porous media[END_REF]:

f c = min( ηΦ 2πaρ f κ ) (2.39)
where η, Φ, a, ρ f and κ are given above in the parameter table.

In our study of Chapter 3, our maximum frequency range f max of the source is such that

f max < f c .
Therefore, in the experimental and numerical modeling of unconsolidated granular media under study, we will choose to stay in the poroelastic regime without attenuation.

Types of waves

We distinguish two main types of waves: body and surface waves. The first propagate throughout the volume of the Earth and pass through its inner layers. And the second propagate only on the surface of the Earth.

• Body waves include P-waves and S-waves (Figure 2.1). P-waves or Primary (Pressure) waves are faster than all other waves (including S-waves or Secondary (Shear)-waves). .

Then they are the first to be recorded by seismometers. P-waves are polarized in the wave propagation direction while S-waves are polarized perpendicularly to the direction of propagation but propagating in the same direction as the P-waves (see Figure 2.1).

• Surface waves include Love waves and Rayleigh waves (Figure 2.2). Surface waves are able to travel only within a few seismic wavelengths of the surface of a solid. Often, the term ground roll is used to describe Rayleigh waves that are trapped in the shallowest layers close to the surface.

Body waves

The wave equation in a homogeneous isotropic elastic medium can be written in a vector form as follows:

ρ ∂ 2 ⃗ U ∂t 2 = (λ + 2µ) ⃗ ∇ div ⃗ U -µ ⃗ rot ⃗ rot ⃗ U (2.40)
The equation (2.40) can be rearranged to give: .

∂ 2 ⃗ U ∂t 2 = λ + 2µ ρ ⃗ ∇ div ⃗ U - µ ρ ⃗ rot ⃗ rot ⃗ U (2.41)
The displacement vector ⃗ U can be written in a decomposed form such as:

⃗ U = ⃗ ∇Φ + ⃗ rotΨ (2.42)
where Φ and Ψ are the dilatational and rotational potentials respectively (Helmholtz potentials). Substituting the expression (2.42) into the equation (2.41) and applying the following formula of the Laplacian operator ∇ 2

∇ 2 ⃗ U = ⃗ ∇ div( ⃗ U ) -⃗ rot ⃗ rot ⃗ U , (2.43) we obtain ∇ λ + 2µ ρ ∇ 2 Φ - ∂ 2 Φ ∂t 2 + ⃗ rot µ ρ ∇ 2 Ψ - ∂ 2 Ψ ∂t 2 = 0.
(2.44)

For P-waves, S-waves are not present ( ⃗ rotΨ = 0) and the displacement vector ⃗ U is then written as ⃗ ∇Φ. In this case, the equation (2.44) becomes the following wave equation for P-waves (2.45): . It can also be expressed as :

λ + 2µ ρ ∇ 2 Φ = ∂ 2 Φ ∂t 2 , (2.45)
∇ 2 Φ = 1 V 2 p ∂ 2 Φ ∂t 2 (2.46)
where the velocity of the P waves is formulated by :

V p = ( λ + 2µ ρ ) (2.47)
By following the same idea, for a pure rotational motion ( ⃗ ∇Φ = 0), the vector of displacement is written as ⃗ U = ⃗ rot ⃗ Ψ, and the wave equation for S-waves becomes:

µ ρ ∇ 2 Ψ = ∂ 2 Ψ ∂t 2 . (2.48)
This equation can be expressed as :

∇ 2 Ψ = 1 V 2 s ∂ 2 Ψ ∂t 2 (2.49)
From equation (2.49), we can easily deduce the velocity of the S waves:

V s = ( µ ρ ) (2.50)
Another parameter describing the elastic properties is the Poisson's ratio ν which can be expressed as a function of P and S velocities as follows:

ν = V 2 p -2V 2 s 2 V 2 p -V 2 s (2.51)
Another parameter is the ratio Vp Vs which is used in numerical modeling for geotechnical or hydrogeophysical applications to identify the presence of fluids [START_REF] Hamada | Reservoir fluids identification using vp/vs ratio?[END_REF]. Values of this ratio are generally smaller than 5 in sediments, above or around water tables (phreatic nappes).

Surface waves a) Rayleigh waves

Discovered by Lord Rayleigh [START_REF] Rayleigh | On waves propagated along the plane surface of an elastic solid[END_REF], Rayleigh waves are surface seismic waves that propagate in the direction of propagation following an ellipsoidal motion in the plane perpendicular to the surface. They are the result of P and SV waves interference and carry most of the seismic energy during an earthquake. The further away from the surface, the more their amplitude decrease (exponentially). Considering a plane (O, x, z) with the z axis oriented downwards and the free surface (zero normal traction at the surface) positioned at z = 0 (O being the origin of the coordinate system), the dilatational and rotational potentials ν Vs Vp V R Vs 0.00 0.707 0.862 0.25 0.507 0.919 0.33 0.500 0.932 0.40 0.408 0.941 0.50 0.000 0.955 are expressed as (section 7.4 of [Aki and Richards, 2002]):

Φ = A exp -rkz exp ik(x-V t) Ψ = B exp -skz exp ik(x-V t) (2.52)
where r and s are two positive constants and k is the wave-number. By injecting these expressions respectively into equations (2.49) and (2.50), we search for a solution propagating at speed V. This solution must be progressive in the horizontal direction (Ox) and stationary in depth (Oz). The search for solutions A and B of the system obtained with the free surface boundary condition thus leads to the Rayleigh wave dispersion equation:

2 - V 2 V 2 s 2 = 4 V 2 V 2 p V 2 V 2 s (2.53)
where V = V R the Rayleigh wave velocity which is lower than the surface waves velocities (0

< V R < V s < V p ).
[ [START_REF] Landau | Hydrodynamic fluctuations[END_REF] proposes another Rayleigh wave equation that shows that Rayleigh waves travel at a speed always lower than the S wave:

V R V s 6 -8 V R V s 4 + 8 V R V s 2 3 -2 V 2 s V 2 p -16 1 - V 2 s V 2 p = 0 (2.54)
As the ratio Vs Vp can be expressed as a function of Poisson's ratio : [START_REF] Viktorov | Rayleigh and Lamb waves: physical theory and applications[END_REF] proposes an approximation of the ratio V R Vs as a function of the Poisson ratio ν:

V s V p = 1 -2ν 2 (1 -ν) (2.55)
V R V s = 1.12ν + 0.87 1 + ν (2.56
) [START_REF] Harker | Velocity and attenuation of ultrasound in suspensions of particles in fluids[END_REF] give in table 2.2 the typical values of the Poisson's ratio and the corresponding Vs Vp and V R Vs ratios. In a layered medium (like the crust and the upper mantle) the velocity of the Rayleigh waves depends on the frequency and wavelength content of the medium. Here are some impor-tant properties of the Rayleigh wave:

• The Rayleigh wave at the surface of a flat homogeneous medium is non-dispersive (independent of the frequency). The wave will only be dispersive if heterogeneities of elastic properties are present in the thickness where the Rayleigh wave propagates. It is the case with a thin layer of slow speed located between the free surface and a homogeneous half-space of faster speed located below.

• The presence of a topography can induce a dispersion of Rayleigh waves even if the medium is homogeneous [START_REF] Pilant | On s-wave directivity patterns[END_REF]. This dispersion depends strongly on the curvature of the surface which can, when it is very strong, convert Rayleigh waves into body waves [START_REF] Rulf | Rayleigh waves on curved surfaces[END_REF].

• The Rayleigh waves penetrate over a certain thickness close to the free surface and this is the reason why they give some information about the V s properties close to the surface. Methods like fundamental modes inversions can provide such information on V s .

• The Rayleigh wave is easily recognizable by its retrograde elliptical polarization in the propagation plane. A three-component sensor placed on the ground surface will therefore be able to record both the horizontal component parallel to the direction of propagation, and the vertical component. And this will allow us to detect Rayleigh waves more easily and to invert them to obtain V s properties.

b) Love waves

Love waves discovered by Augustus Edward Hough Love [START_REF] Love | Some problems of geodynamics[END_REF] are horizontally polarized shear waves (SH waves). These waves only appear in a semi-infinite heterogeneous medium where a low velocity layer covers one or more higher propagation velocity layers. Among seismic waves, Love waves generate the most energy during an earthquake and are therefore the most destructive.

c) Stoneley waves

A Stoneley wave [START_REF] Stoneley | Elastic waves at the surface of separation of two solids[END_REF] or interface wave is a type of boundary wave that propagates along a solid-fluid or solid-solid boundary. Amplitudes of Stoneley waves have their maximum values at the boundary between the two contacting media and decay exponentially towards the depth of each of them.

All the different types of waves (body and surface/interface waves) are summarized in Figure 2.3.

Source type

The complexity of the characterization of the source is closely related to the technique used to generate excitement, but also to the heterogeneity of the medium, which makes its modeling very difficult. A seismic source can be generated from nature (earthquake, ambiant noise, etc.) or men (hammer on a plate, mass drops, explosives, industrial explosion ...). Numerically, a source of excitation can be introduced into the medium either as a collocated force (then located at a point of the grid), or like an explosion (pure point compression or moment tensor). In theory, the time dependence of the source can be any signal with a limited frequency band. Most of the time, as is often the case in geophysics during seismic surveys, a Ricker function (second derivative of Gaussian in time) can be used because it is close to the source signal of a hammer hit pulse. We recall here the Gaussian and the Ricker functions:

• Gaussian:

g(t) = A exp -π 2 f 2 0 (t-t 0 ) 2 (2.57) • Ricker: R(t) = A 1 -2π 2 f 2 0 (t -t 0 ) 2 exp -π 2 f 2 0 (t-t 0 ) 2 (2.58)
with f 0 the dominant frequency of the source spectrum and t 0 = c f 0 a numerical time delay allowing the source to trigger with almost zero amplitude, c being a coefficient close to 1.

For earthquakes or explosive sources, the source can be introduced as a seismic moment tensor M(x, t), tensor of rank three in 3D (two in 2D) having all the properties of a stress tensor, and for which the equivalent force is given by:

s(x, t) = -div M(x, t) (2.59)
where

M =      M xx M xy M xz M yx M yy M yz M zx M zy M zz      in 3D and M =   M xx M xy M yx M yy   in 2D.
For an explosion, all the off-diagonal terms M ij (i ̸ = j) are vanished and P-waves are thus generated around the isotropic source.

In this thesis, we focus on point sources: we test more particularly the velocity implementation (through the terms s x , s y and s z in the equations (2.15)).

Boundary conditions

Free surface

To capture surface waves, i.e Rayleigh waves, a free surface must be implemented through Neumann-type conditions. In the case of the spectral elements, it is naturally implemented thanks to the variational formulation. Nevertheless, in the case of finite differences and finite volumes, it can be implemented by considering anti-symmetric conditions or non-centered space derivatives at the surface. This corresponds to a zero normal traction on the surface. In 2D elastic wave equations, it can be expressed as follows:

σ.⃗ n = 0, (2.60)
where ⃗ n is the external normal vector at the free surface.

This corresponds to σ zz = σ xz = 0 if the free surface is considered as the upper edge of the computational domain.

Absorbing conditions

To prevent reflections at the outer-boundaries of the domain, absorbing conditions are considered. Paraxial absorbing conditions like those of [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF] (first order) or that of [START_REF] Stacey | Improved transparent boundary formulations for the elastic wave equation[END_REF] (second order) which is an improvement of the previous one are generally used. The latter essentially absorb the waves arriving perpendicular to the edge. However, it is not very satisfactory. High-order paraxial methods [START_REF] Bécache | Stability of Perfectly Matched Layers, group velocities and anisotropic waves[END_REF]) or PMLs (Convolutional Perfectly Matched Layer C-PMLs or Auxiliary Differential Equations ADE-PMLs) for example are preferred.

C-PML formulations have been developed for different media to optimize the absorption of the waves impinging the outer PML layers at grazing incidence. For more details we refer to [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]; [START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF]; [START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave equation in poroelastic media[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF].

The non-convolutional ADE-PML [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF] has also been applied and is a more flexible generalization of the CPML formulation. Indeed, it is applicable to higher order timestepping schemes more easily. We adapted the CPMLs and ADE-PMLs to the finite volume technique by applying the stretching functions and the memory variables directly to the spatial derivatives of the numerical fluxes (see equations (3.36) and (3.37)) instead of the derivatives of the particle velocities and stresses. This formulation is more well suited and adapted for the finite volume method when compared to the classical CPML and ADE-PML formulations, because our formulation is able to deal with the stretching of the numerical fluxes derivatives and their related PML memory variables that are both based on left and right Riemann-based state vectors and variables. Besides, a null flux at the outer grid points is assigned instead of classical Dirichlet conditions. In section 2.5.3, we summarize the ADE-PML and C-PML formulations for the 2D elastic wave equations. The Stacey conditions are also mentioned in section 2.5.4. In all the following numerical tests, we only show the simulated solutions with CPML conditions because ADE-PML as well as hybrid Stacey/C-PML (or ADE-PML) are not bringing substantial improvements of the solutions.

ADE-PML and C-PML absorbing boundary conditions for 2D elastic wave equations

In this section we recall the main steps of the design of the C-PML and ADE-PML conditions. The C-PMLs are generally firstly designed in the frequency domain with the elastic wave equations (see eq. 3.1) expressed as:

iωρv x = ∂ x σ xx + ∂ z σ xz , iωρv z = ∂ x σ xz + ∂ z σ zz , iωσ xz = µ (∂ x v x + ∂ z v z ) , iωσ xx = (λ + 2µ) ∂ z v z + λ∂ x v x , iωσ zz = λ∂ z v z + (λ + 2µ) ∂ x v x .
(2.61)

Hereafter we explain how the C-PML and ADE-PML are formulated both in the frequency and time domain.

a) C-PML formulation:

The CPML technique consists in replacing all the spatial derivatives ∂ ∂x by 1 sx ∂ ∂x along x-direction (same thing along z-direction) and finding a better choice of the stretching function s x by introducing variables α x ≥ 0 and K x ≥ 1 such that:

s x = K x + d x α x + iω . (2.62)
By sake of simplicity we take K x = 1 because such a choice is satisfactory for many seismic wave propagation problems.

We denote by sx (t) the inverse Fourier transform of 1 sx . Then, ∂ x can be expressed as:

∂ x = sx (t) * ∂ x .
(2.63)

After calculating the expressions of sx (t) and after some algebraic operations we define

ζ x (t) = -d x H(t)e -(dx+αx)t , (2.64)
where H is the Heaviside distribution, and we thus obtain:

∂ x = ∂ x + ζ x (t) * ∂ x .
(2.65)

Since we have null initial conditions, we can approximate the convolution term at time step n following the recursive convolution method of [START_REF] Luebbers | FDTD for Nth-order dispersive media[END_REF] by:

ψ n ≃ (ψ x ∂ x ) n ≃ N -1 m=0 Z x (m)(∂ x ) n-m , (2.66) CHAPTER 2. SEISMIC WAVE PROPAGATION with Z x (m) = (m+1)∆t m∆t ζ x (t)dτ, = -d x (m+1)∆t m∆t e -(dx+αx)τ dτ, = a x e -(dx+αx)m∆t.
(2.67)

Setting

b x = e -(dx+αx)∆t , (2.68)
and

a x = d x d x + α x (b x -1) , (2.69)
the convolution term ψ x acts as a memory variable on a function f (velocities and stress components) updated at each time step n as:

ψ n x (f ) = b x ψ n-1 x (f ) + a x (∂ x f ) n-1/2 , (2.70)
which means that the unsplit CPML formulation can easily be implemented in an existing finite-difference code by simply replacing the spatial derivatives ∂ x with ∂ x + ψ x and advancing ψ x in time using the same time evolution scheme as for the other (existing) variables.

b) ADE-PML formulation:

As for the CPML, the spatial derivatives along the axis perpendicular to the PML layer, say x, are rewritten in a stretched coordinate x [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]) :

∂ x = 1 s x ∂ x , (2.71)
where

s x = K x + d x α x + iω . (2.72)
Then, following ?, we can express 1 sx as follows :

1 s x = 1 K x - d x K 2 x 1 dx Kx + α x + iω , (2.73) with d x = d 0 x L N , K x = 1 + (K max -1) m , α x = α max 1 - x L p , (2.74)
where L is the thickness of the PML layer, N = 2, d 0 = -(N +1)Cpmaxlog(Rc)

2L

, C pmax the maximum speed of the pressure wave and R c being the target theoretical reflection coefficient, chosen here as 0.1% (see [START_REF] Collino | Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]).

We also take p = 1 and α max = πf 0 , where f 0 is the dominant frequency of the seismic source. K max usually lies between 1 and 20 [START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF]) and we take here K max = 1 and m = 2. Therefore, as an example and without loss of generality, the spatial derivative

∂ x σ xz is transformed in the CPML layer into 1 s x ∂ x σ xz = 1 K x ∂ x σ xz - d x K 2 x 1 dx Kx + α x + iω ∂ x σ xz . (2.75)
Let us denote by Q σxz x the auxiliary memory variable associated with ∂ x σ xz which is expressed as :

Q σxz x = - d x K 2 x 1 dx Kx + α x + iω ∂ x σ xz , (2.76)
which leads to

d x K x + α x + iω Q σxz x = - d x K 2 x ∂ x σ xz .
(2.77)

In the time domain this equation now becomes

∂ t Q σxz x + d x K x + α x Q σxz x = - d x K 2 x ∂ x σ xz .
(2.78) Similar formulations are derived for the x and z derivatives of v x , v z , σ xx and σ zz in 2D and consequently the whole system of equations for elastodynamics with ADE-PML conditions is given by:

∂ t σ xx = (λ + 2µ) 1 K x ∂ x v x + Q vx x + λ 1 K z ∂ z v z + Q vz z , ∂ t σ zz = (λ + 2µ) 1 K z ∂ z v z + Q vz z + λ 1 K x ∂ x v x + Q vx x , ∂ t σ xz = µ 1 K x ∂ x v z + Q vz x + 1 K z ∂ z v x + Q vx z , ρ∂ t v x = 1 K x ∂ x σ xx + Q σxx x + 1 K z ∂ z σ xz + Q σxz z , ρ∂ t v z = 1 K z ∂ x σ xz + Q σxz x + 1 K z ∂ z σ zz + Q σzz z .
(2.79)

The differential equations associated to the memory variables are discretized using here a second order time stepping scheme. For instance, the equation associated to the memory variable Q σxz x is expressed as follows

∂ t Q σxz x + 1 τ x Q σxz x = - d x K 2 x ∂ x σ xz , (2.80)
where τ x = 1 dx Kx +αx . After a second-order time stepping scheme, we have

Q σxz x (t n ) = b x Q σxz x (t n-1 ) + a x ∂ x σ n-1/2 xz , (2.81) where b x = exp -∆t τx and a x = τ x dx K 2 x (1 -b x ).
We can retrieve the CPML formulation (see equation (2.70)) from the previous equation (2.81). Indeed, we can deduce a x and b x of the C-PML formulation by using second-order Taylor expansion of the exponentials introduced in the a x and b x coefficients involved in the time-stepping scheme (2.81) applied to the variables Q. In other words, the a x and b x of ADE-PML are second-order Taylor expansions of a x and b x of C-PML formulation. The difference between ADE-PML and C-PML formulations lies in the fact that the ADE formulation can be extended to a higher-order time scheme.

Stacey conditions

As an example, on the Model B in Chapter 3 (see Figure 3.5b) with a free surface, we have applied Stacey absorbing boundary conditions. It improved a little bit the solutions but we preferred to apply PML conditions because they were accurate enough for our studies. In this purpose, the normal tractions are applied at the different outer boundaries (left: L, right: R, top: T, bottom: B) and reformulated as

σ.⃗ n = T (2.82)
where T is the traction, σ the stress matrix and ⃗ n the outgoing normal given by

σ =   σ xx σ zx σ xz σ zz   , ⃗ n = n x n z T , T = T x T z T and T = V p (ρV ) T + V s (ρV ) N , (2.83)
where (ρV ) N , (ρV ) T are respectively the normal and tangential components of the momentum on the absorbing interface, V p and V s being the P and S-wave velocities. At the different boundaries, we thus obtain:

at the left boundary L:

⃗ n L = 0 -1 T , σ xx = -ρV p v x , σ xz = -ρV s v z .
(2.84) at the right boundary R:

⃗ n R = 0 1 T , σ xx = ρV p v x , σ xz = ρV s v z .
(2.85) at the bottom boundary B:

⃗ n B = -1 0 T , σ zz = -ρV p v z , σ zx = -ρV s v x .
(2.86) at the top boundary T:

⃗ n T = 1 0 T , σ zz = ρV p v z , σ zx = ρV s v x .
(2.87)

First arrivals travel-time computation in stratified media

In seismology, ray theory is a method used to approximate the wave propagation in some complex media. In particular, this method is applied to model high-frequency seismic body waves in a stratified medium whose physical properties vary continuously. Indeed, the approximation obtained by ray theory at all distances and depths in the case of unconsolidated sands (our object of study) are satisfactory [START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, part ii: Velocities in shallow unconsolidated sand[END_REF]], [START_REF] Bachrach | Seismic velocities and poisson's ratio of shallow unconsolidated sands[END_REF]. Theoretical first arrival travel-time can then be predicted [START_REF] Bachrach | Seismic velocities and poisson's ratio of shallow unconsolidated sands[END_REF] to determine a P-wave velocity model varying with depth. This is what we will use in Chapter 3 for experimental cases with power-law variations of V p with depth.

The travel distance X and traveltime T of the first arrivals can be expressed as functions of the ray parameter p and the velocity V = V (Z) [Aki and Richards, 2002]; [START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, part ii: Velocities in shallow unconsolidated sand[END_REF]]:

X = 2 Zmax 0 pV dZ 1 -p 2 V 2 (Z) (2.88) T = 2 Zmax 0 dZ V 1 -p 2 V 2 (Z) (2.89)
where Z max is the maximum depth at which the ray penetrates. The ray parameter p, during the travel time of the seismic wave, can be expressed as:

p = 1 V (Z max ) = sin i 0 V (Z = 0)
(2.90) where i 0 is the incidence angle, V (Z max ) and V (Z = 0) are the velocities of the deepest layer (highest velocity) and the surface layer (lowest velocity) respectively (Figure 2.4). The ray parameter p is calculated using Snell's Law (Figure 2.5) which is deduced from Fermat's principle:

sin

θ inc V pA = sin θ t V pB (2.91)
If the angle of incidence is greater than a critical angle θ c , the wave is totally reflected and is not refracted at all [Aki and Richards, 2002]. If the incident angle is lower than the critical angle, then the wave is totally reflected. If the incident angle is equal to the critical angle, then the wave propagates along the interface.

Multi-channel analysis of surface waves (MASW)

In the near surface context, and among the different seismic methods for determining the structure and physical properties of a medium, the MASW technique is used to study the propagation of surface waves (Rayleigh and Love). In a stratified medium (vertically heterogeneous), the propagation of Rayleigh wave depends on the frequency (dispersive). Therefore, by extracting dispersion information of the medium, one can determine the variation of the shear wave speed (V s ) as a function of the depth. The data acquisition profile of this method consists of a seismic source (weight drop or seismic hammers for example) and receivers placed at the surface along the propagation direction in order to record the propagation speeds of the waves (usually the vertical component).

For instance, [START_REF] Nazarian | Nondestructive testing of pavements using surface waves[END_REF][START_REF] Nazarian | Use of spectral analysis of surface waves method for determination of moduli and thicknesses of pavement systems[END_REF], set up an acquisition device composed of only two receivers and a seismic source generated by a sledge hammer. This technique consists in calculating the Fourier Transforms (FFT) of phase spectra of surface waves and is called Spectral Analysis of Surface Waves (SASW). The corresponding wavenumber related to the configuration of the acquisition device is too low (insufficient) to calculate accurately all the propagation modes of surface waves which may occur.

Then, [START_REF] Mcmechan | Analysis of dispersive waves by wave field transformation[END_REF]; [START_REF] Gabriels | In situ measurements of shear-wave velocity in sediments with higher-mode rayleigh waves[END_REF]; [START_REF] Park | Predictions of shear-wave velocities in southern california using surface geology[END_REF]]; [START_REF] Foti | Multistation methods for geotechnical characterization using surface waves[END_REF] propose to increase the wavenumber resolution by increasing the number of receivers and by reducing the distance between them. This gave birth to MASW: Multichannel Analysis of Surface Waves; a technique that makes the data processing much faster and more robust than the SASW method. The dispersion curve of Rayleigh wave phase velocity can be computed in different ways:

• the f-k transform [START_REF] Gabriels | In situ measurements of shear-wave velocity in sediments with higher-mode rayleigh waves[END_REF]; [START_REF] Foti | Multistation methods for geotechnical characterization using surface waves[END_REF].

• the τ -p transform [START_REF] Mcmechan | Analysis of dispersive waves by wave field transformation[END_REF]].

• the phase-difference method called also the slant-stack transform [START_REF] Mokhtar | Seismic velocity and q model for the shallow structure of the arabian shield from short-period rayleigh waves[END_REF]; [START_REF] Park | Predictions of shear-wave velocities in southern california using surface geology[END_REF]].

• the linear Radon transform [START_REF] Luo | Rayleigh-wave dispersive energy imaging using a high-resolution linear radon transform[END_REF].

In our study, we choose the slant-stack transform to extract the dispersion curve in frequencyvelocity phase domain. This method optimizes the τ -p transform and provides a high spectrum resolution with an optimized number of receivers [Park et al., 1998]; [START_REF] Socco | Surface-wave analysis for building near-surface velocity models-established approaches and new perspectives[END_REF]; [START_REF] Xia | Estimation of near-surface shear-wave velocities and quality factors using multichannel analysis of surface-wave methods[END_REF].

Once the dispersion curves have been calculated, they are inverted to estimate S-wave velocity as a function of the depth of the medium. With MASW methods, we can reconstruct a medium up to 30 m deep [START_REF] Park | Masw horizontal resolution in 2d shear-velocity (vs) mapping[END_REF][START_REF] Park | Optimum field parameters of an masw survey[END_REF]. This depth is determined by the the configuration of the acquisition device chosen. Now, we develop the slant-stack transform. The particle motion is recorded at each time t and at each receiver position x, giving D(t, x) a data set in the time-offset domain, which is also called seismogram. Then, a 1D Fourier Transform in time is applied at each receiver position x, giving a 2D f -x spectrum D(f, x):

D(f, x) = +∞ -∞ D(t, x) exp (i2πf t) dt (2.92)
which can be rewritten as:

D(f, x) = P (f, x)A(f, x) (2.93)
where P (f, x) and A(f, x) = | D(f, x)| are the phase and the amplitude spectrum respectively.

As in Park et al. [1996] (Dynamic Linear Move Out/DLMO approach), we apply an operator Φ(f, x), depending on the offset, to the equation (2.93) :

Φ(f, x) = exp -i2πf x V ph (2.94)
where V ph is the phase velocity calculated from a shot gather obtained at a reference location. This allows to correct at the same time the possibly existing offset effects and the linearly sloping events of waves trapped in the very near surface layers (ground roll). Thus, the phase-difference method [START_REF] Park | Predictions of shear-wave velocities in southern california using surface geology[END_REF]] calculates the dispersion diagram D(f, V ph ) as follows:

D(f, V ph ) = x D(f, x) | D(f, x)| Φ(f, x) dx = x P (f, x) Φ(f, x) dx (2.95)
Finally, the discrete form of equation (2.95) for a discrete offset x i is

D(f, V ph ) = ∆x N x i D(f, x i ) | D(f, x i )| Φ(f, x i ) = ∆x N x i P (f, x i ) Φ(f, x i ) (2.96)
where N x is the number of receivers and ∆x is the uniform receivers' space. The amplitude of the dispersion diagram is normalized by the maximum value at each frequency. The phase velocity dispersion curves corresponding to different modes are obtained by selecting these maxima (see appendix A.5).

Conclusion

In the following chapters, we develop the numerical methods that solve the wave equations in an elastic and poroelastic medium. A finite volume code based on a Riemann solver is validated on one hand by considering different complex synthetic cases, and on the other hand by comparing the solutions with those of SPECFEM (spectral elements) and seismic CPML (finite differences) that are two reference codes for wave propagation simulations. Once the numerical tools have been validated, we use them to reproduce an experiment carried out in the laboratory on unconsolidated granular media and partially saturated with fluid. The signal processing tools developed at the end of this chapter are used to detect the first arrival traveltimes and the principal modes of propagation of the waves, and to determine the best model characterizing the media.

Chapter 3

Modeling of elastic wave propagation using a Riemann-based Finite Volume method: comparisons with Finite Difference and Spectral Element methods

Introduction

Interpretation of seismic data requires a numerical resolution of two-or three-dimensional seismic wave propagation equations for different rheologies : elastic, poroelastic, viscoelastic media etc. The choice of an appropriate numerical method with regards to its accuracy, flexibility, performance and ease of implementation depends essentially on the physical problem to be solved. More particularly, we are interested here in geophysical problems, for which it is not easy to solve accurately strong gradients of the seismic properties with depth as well as both the surface waves on topographies and the absorbing boundaries.

There are many numerical methods to solve the wave equations, of which mainly including, the finite difference method/FDM, the pseudospectral method/PSM, the finite element method/FEM, the spectral element method/SEM, the discontinuous Galerkin method/DG, and the finite volume method/FVM. Each of these numerical methods has its own advantages and disadvantages. In the following, we first give a brief review on these methods.

• The finite difference method/FDM is one of the most used numerical methods for solving wave propagation equations. The mesh is composed of points and the solutions are calculated using a partial derivative approximation (Taylor expansion). Depending on the arrangement of the points (mesh), and the order of accuracy (in time and space), different finite difference schemes can be derived, including one of the most famous: the staggeredgrid scheme [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggeredgrid finite differences[END_REF][START_REF] Sei | A family of numerical schemes for the computation of elastic waves[END_REF][START_REF] Virieux | SH wave propagation in heterogeneous media: velocity-stress finite-difference method[END_REF][START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF]. High order schemes are easy to construct and to implement. The FDM also demands a lower computing cost and fits parallel computing. All of these reasons make this method very useful. However, the inherent restriction of using regular meshes limits its application to complex topography: for instance surface waves such as Rayleigh waves are difficult to compute with good resolution. Furthermore, the accurate calculation of wavefield propagation solutions that involve interfaces, such as free surface topography or fluid -solid interface, is a challenge for finite difference [Gibson et • The pseudospectral method/PSM commonly uses regular meshes and can be viewed as the limit of the FDM with infinite-order accuracy in space [START_REF] Fornberg | High-order finite differences and the pseudospectral method on staggered grids[END_REF][START_REF] Fornberg | A practical guide to pseudospectral methods[END_REF] [START_REF] Seriani | [END_REF]] is also commonly used to solve the wave propagation equations. Since it uses quadrangle/hexahedra or triangle/tetrahedra elements, the FEM is flexible and is thus well suited to model complex topography, heterogeneous media and interfaces between various materials. It allows also a natural treatment of free surface conditions. However, its solutions are essentially continuous which is not always the case in practice, particularly when the medium's properties change sharply. To circumvent this problem, a discontinuous Galerkin method may be used. However, the latter is difficult to implement, particularly for a 3D problem. Solving 3D problems with the FEM method is very expensive in terms of computational ressources/CPU time and computer memory because a global mass matrix must be inverted at each time step.

• To overcome this mass matrix issue, the spectral element method/SEM originally introduced by Patera [1984a] for computational fluid mechanics has been developped for wave propagation modeling [START_REF] Bernardi | Mixed spectral element approximation of the Navier-Stokes equations in the stream-function and vorticity formulation[END_REF] The FVMs which can be seen as a generalization of the FDM, are typically first order in space. However, in the last decade, high order FVM schemes have been developped on unstructured grids using a high-order polynomial reconstruction technique [START_REF] Qian | Coupled THINC and level set method: A conservative interface capturing scheme with high-order surface representations[END_REF]. Furthermore, this method is able to preserve the cell averages for various orders of interpolation. Besides, the Runge-Kutta (RK) schemes are currently used for time discretization. In order to simulate wave propagation on unstructured meshes efficiently, the FVM is particularly interesting due to its good adaptivity to complex geometries, unlike the FD methods.

We thus opted for a finite volume method with Riemann solvers and flux limiters which is able to meet all the needs mentioned above: this offers a compromise between unstructured FVM methods and FD methods. Furthermore, we adapted the formulation of perfectly matched layer/PML absorbing conditions to the FVM, which allows us to avoid the reflection of the waves at the outer boundaries. These PMLs are inspired from the frequency shift convolutional (CPML) or non-convolutional (ADEPML) PMLs of [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF]]. Here we introduce a fluxes-based formulation instead of the primitive variable formulation of the CPMLs or ADEPMLs to take into account the Riemann state vector-based fluxes at the cell interfaces.

In a first step, we solve the two-dimensional wave equation in an elastic (linear) medium with a flat topography using the FVM. Furthermore, we implement and validate a fluxes-based formulation of the frequency shift convolutional (CPML) or non-convolutional (ADEPML) perfectly matched layer absorbing conditions at the outer boundaries [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF], and also a free surface boundary condition that physically represents the ground. The validation of the developed numerical tool consists in comparing the solutions obtained in different heterogeneous checkerboard-like test cases (with or without free surface) to those obtained with codes such as SPECFEM (spectral finite elements) or Seismic_CPML (finite differences) [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF][START_REF] Martin | Seismic wave propagation in nonlinear viscoelastic media using the auxiliary differential equation method[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF], SPECFEM being considered as a reference code in the geophysical community. In a second step, we implement a fluid-solid coupling and compare FVM to FD and SEM solutions. In a third step, we suggest a more complex velocity structure, e.g. induced by a strong gradient of properties with depth as involved in unconsolidated granular media [START_REF] Bergamo | P-and S-wave velocity models of shallow dry sand formations from surface wave multimodal inversion[END_REF]

Governing equations and their finite volume discretization

The 2D elastic wave equations for an isotropic medium submitted to external forces can be written using a velocity-stress formulation such as the following linear and hyperbolic system [START_REF] Dumbser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-II. The three-dimensional isotropic case[END_REF]):

ρ ∂v x ∂t - ∂σ xx ∂x - ∂σ xz ∂z = s x , ρ ∂v z ∂t - ∂σ xz ∂x - ∂σ zz ∂z = s z , ∂σ xx ∂t -(λ + 2µ) ∂v x ∂x -λ ∂v z ∂z = 0, ∂σ zz ∂t -λ ∂v x ∂x -(λ + 2µ) ∂v z ∂z = 0, ∂σ xz ∂t -µ ∂v z ∂x + ∂v x ∂z = 0, (3.1)
where λ and µ are Lamé parameters, ρ is the density, and s x and s z are the space dependent source terms in x and z directions. The compressional stress components are given by σ xx and 3.2. GOVERNING EQUATIONS AND THEIR FINITE VOLUME DISCRETIZATION 43 σ zz , and the shear stress is σ xz . The components of particle velocities in direction x and z are denoted by v x and v z , respectively.

The compressional wave or P -wave velocity and the shear wave or S-wave velocity are given by V p = ( λ+2µ ρ ) and V s = ( µ ρ ), respectively. In this section, we propose a finite volume method that is designed for hyperbolic equations and which uses Riemann solvers and wave (or flux) limiters whose main role is to reduce numerical oscillations that occur at the discontinuities (see LeVeque [2002a]). We first present the method in the case of a 1D hyperbolic system, and then we extend it to the 2D case by applying dimension by dimension or split methods.

1D hyperbolic system in conservation form

A 1D hyperbolic system in conservation form can be written as follows:

∂ ∂t Q(x, t) + ∂ ∂x F (Q(x, t)) = 0, (3.2)
where Q is a vector of m components (of the conserved variables) and F the physical flux vector. If F is linear then we can rewrite the above system in a compact form:

∂ ∂t Q(x, t) + A ∂ ∂x Q(x, t) = 0, (3.3)
where A = ∂ ∂Q F is the Jacobian matrix. If F is non-linear then we can linearize the system locally by using a Roe matrix (see for example [START_REF] Toro | Godunov methods: Theory and applications[END_REF]).

In particular, in the case of 1D wave equations, we have m = 2, Q = v x σ T and

A =   0 -1 ρ -(λ + 2µ) 0   . (3.4)

a) Riemann problem

The Riemann problem associated to the 1D linear hyperbolic system of m equations (3.3) is the following initial value problem:

       ∂ ∂t Q(x, t) + A ∂ ∂x Q(x, t) = 0, Q(x, 0) = Q 0 (x), (3.5)
such that at each interface defined along x we denote by Q L and Q R the left and right states vectors of Q. To find a Riemann solver, we study the eigenstructure of the problem: we decompose the Jacobian matrix A by diagonalization as follows: where R x and (R x ) -1 are matrices that contain the right and left eigenvectors, respectively and Z x is a diagonal matrix containing the eigenvalues ζ p , p = 1, ..., m, of matrix A.

A = R x Z x (R x ) -1 , (3.6)
The exact Riemann solver is obtained by decomposing the jump ∆Q = Q R -Q L at the cell interfaces into eigenvectors of A and is expressed as:

α = (R x ) -1 ∆Q. ( 3.7)
We denote by r p the right p-th eigenvector contained in R x and by l p the left p-th eigenvector contained in (R x ) -1 . The components of α can be written as:

α p = l p ∆Q. ( 3

b) Time and space discretizations

In one space dimension, the mesh of the finite volume method is composed of intervals (surfaces in 2D and volumes in 3D). At each time step, we calculate the flux F at the endpoints of the intervals to approach the integral of Q on each of these cells.

We assume that the grid is uniform. We consider an interval [a, b] subdivided into N subintervals and let us denote the i-th cell by C i = x i-1/2 , x i+1/2 (see Figure 3.1). The unknown wavefield computed at the cell located at point x i is denoted by Q i .

The value Q n i approximates the average value over the i-th cell interval at time t n :

Q n i ≈ 1 ∆x x i+1/2 x i-1/2 Q(x, t n )dx ≡ 1 ∆x C i Q(x, t n )dx, (3.11)
where ∆x = x i+1/2 -x i-1/2 is the length of the cell.

By using the properties of the conservation laws (mass and momentum conservation), one can build numerical methods which calculate quite precisely the real solutions (shock waves for example). In our case, the integral of Q we want to approximate, corresponds to the value of Q at the midpoint of the interval to O(∆x 2 ). Taking the conservation law, we can say that the discrete sum N i=1 Q n i ∆x which approximates the integral of Q over the entire interval [a, b], will change only due to fluxes at the boundaries x = a and x = b.

The discretized integral form of the conservative equations can thus be expressed as:

d dt C i Q(x, t n )dx = F (Q(x i-1/2 , t) -F (Q(x i+1/2 , t)).
(3.12)

We integrate (3.12) in time from t n to t n+1 :

C i Q(x, t n+1 )dx - C i Q(x, t n )dx = t n+1 t n F (Q(x i-1/2 , t)dt - t n+1 t n F (Q(x i+1/2 , t))dt. (3.13)
Rearranging (3.13) and dividing by ∆x:

1 ∆x C i Q(x, t n+1 )dx = 1 ∆x C i Q(x, t n )dx+ 1 ∆x t n+1 t n F (Q(x i-1/2 , t)dt - t n+1 t n F (Q(x i+1/2 , t))dt .
(3.14) It is difficult to evaluate the right-hand side of (3.14) because of the variation of Q(x i±1/2 , t) with time. The following form of the numerical schemes is thus suggested:

Q n+1 i = Q n i - ∆t ∆x F n i+1/2 -F n i-1/2 , ( 3.15)
where F n i-1/2 is an approximation of the average flux at the interface located at x = x i-1/2 :

F n i-1/2 ≈ 1 ∆t t n+1 t n F (Q(x i-1/2 , t))dt. (3.16)
At time t n+1 , the cell averages Q n+1 i , i = 1, ..., N are approximated after a time step of length ∆t = t n+1 -t n using Q n i , i = 1, ..., N , the cell averages at time t n . Then, a fully discrete method can be obtained by approaching the numerical flux using the values Q n . In the case of a hyperbolic problem where wave propagates with finite speed, a good approximation of the flux F n i-1/2 can be calculated using only the cell averages Q n i-1 and Q n i on either side (left and right) of the interface [LeVeque, 2002a]. Here the left state

Q n L i-1/2 is Q n i-1
and the right one

Q n R i-1/2 is Q n i .
Thus let us introduce an expression of the numerical flux as follows:

F n i-1/2 = F(Q n L i-1/2 , Q n R i-1/2 ) = F(Q n i-1 , Q n i ), (3.17)
where F is some numerical flux function. The scheme then becomes

Q n+1 i = Q n i - ∆t ∆x F(Q n i , Q n i+1 ) -F(Q n i-1 , Q n i ) . (3.18)
The specific scheme depends on how we choose the expression of F. We compute flux F at each interface by solving a Riemann problem for (Q L , Q R ) states.

In general the flux F at the interface x i-1/2 depends on the local Jacobian matrix

A i-1/2 associated to 1 2 (Q n i-1 + Q n i ), and the jump ∆Q i-1/2 given in equation (3.10) with Q n L i-1/2 = Q n i-1 and Q n R i-1/2 = Q n i .
We can suppose that the solution of the Riemann problem consists of m waves W p traveling at speeds ζ p , each of which may be positive or negative (LeVeque [2002a]). Then the cell average is updated by:

Q n+1 i = Q n i - ∆t ∆x [A -∆Q] i+1/2 + [A + ∆Q] i-1/2 , ( 3.19)
where

[A + ∆Q] i-1/2 = m p=1 ζ p i-1/2 + W p i-1/2 , (3.20) [A -∆Q] i+1/2 = m p=1 ζ p i+1/2 - W p i+1/2 . ( 3

.21)

A + and A -are respectively the positive and negative parts of matrix A evaluated at x i-1/2 and x i+1/2 , and (ζ p i-1/2 ) + = max(ζ p i-1/2 , 0) and (ζ p i+1/2 ) -= min(ζ p i+1/2 , 0) are the respective associated eigenvalues. More exactly, A + and A -are defined as:

A + = R x (Z x ) + (R x ) -1 , (3.22) A -= R x (Z x ) -(R x ) -1 , ( 3.23)
with (Z x ) + and (Z x ) -the matrices containing the eigenvalues (ζ p ) + , (ζ p ) -for p = 1, ..., m respectively.

Equations (3.20) and (3.21) correspond to the right-going fluctuation from the left edge and the left-going fluctuation from the right edge of the cell C i , respectively. From this formulation, a first order expression of the flux F (Godunov method) at the different interfaces x i-1/2 and x i+1/2 can be extracted (see LeVeque [2002a]):

F i-1/2 = A + i-1/2 Q i-1 + A - i+1/2 Q i , ( 3.24)
F i+1/2 = A + i-1/2 Q i + A - i+1/2 Q i+1 .
(3.25)

GOVERNING EQUATIONS AND THEIR FINITE VOLUME DISCRETIZATION 47

In order to introduce a higher resolution method, we add a numerical artificial viscosity term F to the fluxes F given in equations (3.24) and (3.25) which depends also on the Jacobian A and the jump ∆Q = m p=1 W p :

F i-1/2 = A + i-1/2 Q i-1 + A - i+1/2 Q i + Fi-1/2 , ( 3.26)
F i+1/2 = A + i-1/2 Q i + A - i+1/2 Q i+1 + Fi+1/2 (3.27) where Fi-1/2 = 1 2 m p=1 |ζ p i-1/2 | 1 - ∆t ∆x |ζ p i-1/2 | W p i-1/2 . (3.28)
The final scheme is thus given by (see LeVeque [2002a] and [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]):

Q n+1 i = Q n i - ∆t ∆x [A + ∆Q] i-1/2 + [A -∆Q] i+1/2 - ∆t ∆x Fi+1/2 -Fi-1/2 . (3.29)

c) Fluxes and wave limiters

To reduce numerical oscillations, LeVeque [2002a] proposes to multiply the p-waves W p by a limiter function ϕ as follows:

W p i-1/2,j = αp i-1/2,j r p i-1/2,j , (3.30) αp i-1/2,j = ϕ(θ p i-1/2,j )α p i-1/2,j , (3.31)
where θ p i-1/2,j evaluates the regularity of the solution at the interface. The θ function is given by

θ p i-1/2,j = α p I-1/2,j α p i-1/2,j , (3.32) with I =    i -1 if λ p > 0, i + 1 if λ p ≤ 0.
Here are some examples of limiter functions (LeVeque [2002a], van [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. a second-order sequel to Godunov's method[END_REF], [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]) that we will consider in the next sections:

• SMART: ϕ(θ) = max (0, min (0.25 + 0.75θ, 4, 2θ)) (see [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]);

• KOREN: ϕ(θ) = max 0, min 2θ, min 1+2θ 3 , 2 (see [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]));

• MC: ϕ(θ) = max 0, min 1+θ 2 , 2, 2θ (see LeVeque [2002a]);

• Van Leer:

ϕ(θ) = θ+|θ| 1+θ (see LeVeque [2002a]).

2D formulation

In the 2D linear case, we can rewrite equations (3.1) in a compact form :

∂ ∂t Q(x, z, t) + A ∂ ∂x Q(x, z, t) + B ∂ ∂z Q(x, z, t) = s(x, z), ( 3.33)
where

Q = v x v z σ xx σ zz σ xz T , s = s x s z 0 0 0 T , A =            0 0 -1/ρ 0 0 0 0 0 0 -1/ρ -(λ + 2µ) 0 0 0 0 -λ 0 0 0 0 0 -µ 0 0 0            , ( 3.34)
and

B =            0 0 0 0 -1/ρ 0 0 0 -1/ρ 0 0 -λ 0 0 0 0 -(λ + 2µ) 0 0 0 -µ 0 0 0 0            . ( 3.35)
We apply a split or a dimension by dimension method to extend the scheme to the 2D case which is similar to a space and time predictor-corrector scheme.

We denote by F and G the numerical fluxes respectively in x and z directions using wave limiters. We first compute an intermediate solution along the x-direction (as similarly done in the previous 1D case):

Q * i,j = Q n i,j - ∆t ∆x [A + ∆Q] i-1/2,j + [A -∆Q] i+1/2,j - ∆t ∆x Fi+1/2,j -Fi-1/2,j (along x-direction).
(3.36) Then, we inject Q * along the z-direction:

Q n+1 i,j = Q * i,j - ∆t ∆z [B + ∆Q * ] i,j-1/2 + [B -∆Q * ] i,j+1/2 - ∆t ∆z Gi,j+1/2 -Gi,j-1/2 (along z-direction) (3.37) where [A + ∆Q] i-1/2,j = m p=1 ζ px i-1/2,j + α p i-1/2,j r px i-1/2,j , (3.38) [A -∆Q] i+1/2,j = m p=1 ζ px i+1/2,j - α p i+1/2,j r px i+1/2,j , (3.39)
3.2. GOVERNING EQUATIONS AND THEIR FINITE VOLUME DISCRETIZATION 49

[B + ∆Q] i,j-1/2 = m p=1 ζ pz i,j-1/2 + β p i,j-1/2 r pz i,j-1/2 , (3.40) [B -∆Q] i,j+1/2 = m p=1 ζ pz i,j+1/2 - β p i,j+1/2 r pz i,j+1/2 , (3.41) Fi-1/2,j = 1 2 m p=1 |ζ px i-1/2,j | 1 - ∆t ∆x |ζ px i-1/2,j | W px i-1/2,j , (3.42) Gi,j-1/2 = 1 2 m p=1 |ζ pz i,j-1/2 | 1 - ∆t ∆z |ζ pz i,j-1/2 | W pz i,j-1/2 . (3.43)
Here, α and β are the Riemann solutions along the x and z directions respectively, and are given by

α = (R x ) -1 (Q R -Q L ) , ( 3.44)
and

β = (R z ) -1 (Q T -Q B ) , (3.45)
where Q R , Q L , Q T , Q B are respectively the right, left, top and bottom state vectors.

For the 2D elastic wave equations, the expressions of matrices A, B, A

+ , A -, B + , B -, R x , (R x ) -1 , R z , (R z ) -1 , Z + and Z -are presented in appendix A.1.
As in the 1D case, A + and A -are the positive and negative parts of matrix A evaluated at x i-1/2 and x i+1/2 along x-direction (same thing for the eigenvalues ζ px+ and ζ px-of matrices Z + and Z -). Similarly is done for matrices B + , B -, Z + and Z -(eigenvalues ζ pz + and ζ pz -) along z-direction.

Reference solutions

To estimate the accuracy of the proposed FVM method, we choose as reference solutions the spectral element (SPECFEM) solutions that are fourth-order in space and second-order in time (Newmark scheme). And we compare also the FVM to FD/finite difference solutions to see if it gives better results. We recall hereafter the numerical method used in SPECFEM and the second-order [START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF] and fourth-order [START_REF] Graves | Simulating seismic wave propagation in 3D elastic media using staggeredgrid finite differences[END_REF][START_REF] Levander | Fourth-order finite-difference p-sv seismograms[END_REF][START_REF] Moczo | 3-D fourth-order staggered-grid finite-difference schemes: stability and grid dispersion[END_REF]] staggered grid finite difference discretization schemes in space used for different applications.

a) SPECFEM

The SPECFEM code [START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Spectral and spectral-element methods for the 2D and 3D elastodynamics equations in heterogeneous media)[END_REF] uses the spectral element method which is an integral method based on solving the wave equation in its variational formulation (also called weak formulation). The spectral element method is a variant of the finite element method that combines the mesh flexibility of finite elements and the high order spatial precision of spectral methods (use of high order polynomial interpolation) [START_REF] Faccioli | 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method[END_REF][START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF]Patera, 1984b]. The finite and spectral elements are represented in Figure 3.2. This method is based on a variational (or weak) formulation of the wave propagation equations. The solutions are calculated over a Gauss-Lobatto-Legendre (GLL) mesh points, which are roots of the Legendre polynomial. Numerical quadrature rules are designed for this particular mesh. This makes the mass matrix to be diagonal. This is very important because no matrix system has to be solved at each time step, which makes the spectral element method very flexible and easy to parallelize. Furthermore, interfaces and free surface conditions are naturally taken into account in the variational formulation. Hereafter we recall the finite element method for wave propagation equations.

For a bounded domain Ω, the formulation of the conservation of momentum is written:

∇.σ + f = ρa in Ω σ.n = T d on ∂Ω, (absorbing conditions, normal traction imposed, ...) [σ] Σ .n = t = 0 (free surface condition) (3.46) (3.47) (3.48)
where σ is the stress tensor, a the acceleration vector, f the volumic force field, T d the surface force vector, n the normal vector to the surface Σ and t the normal traction vector. The notation [.] represents the jump at a discontinuity interface.

After integration by parts of the stress terms, the conservation equation of the momentum for any virtual displacement (or velocity) test vector û [START_REF] Dhemaied | Modélisation directe de la propagation d'ondes dans les milieux visco-élastiques : approche 2D par différences finies dans le domaine temporel[END_REF][START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Spectral and spectral-element methods for the 2D and 3D elastodynamics equations in heterogeneous media)[END_REF] becomes:

-

Ω σ : ϵ(û)dΩ + Ω f ûdΩ + ∂Ω T d ûd∂Ω + Σ σ.nûdΣ = Ω ρaûdΩ (3.49)
where ϵ is the strain tensor.

A stress tensor expression for an isotropic linear elastic medium is given by the generalized Hooke's law:

σ = C : ϵ (3.50)
where ":" represents the contracted product of two tensors. We apply the natural free surface condition and we replace equation (3.50) in equation (3.51) for any displacement û ∈ T where T is the set of admissible displacements with the initial conditions on u to obtain the variational formulation of the problem:

-

Ω ϵ(û) : C : ϵ(û)dΩ + Ω f ûdΩ + ∂Ω T d ûdΩ - Ω ρaûdΩ = 0 (3.51) u(x, 0) = u 0 (x) u(x, 0) = v 0 (x) (3.52) (3.53)
where u, u ∈ H 1 (Ω), σ ∈ L 2 (Ω) and x ∈ Ω.

Solving the problem is reduced to minimize the whole vector at the left hand side of equation (3.51) denoted V(û) in the space T . To do this, we introduce the approximation by the method of Galerkin on a finite dimensional space denoted T h :

T h = ûh = n i=1 a i N i + u d , a i ∈ R n ⊆ T (3.54)
where ûh can be obtained from any displacement vector u d filling the surface conditions and n functions N i called shape functions. The n scalars a i are obtained by minimizing:

min ûh ∈T h V(û h ) (3.55)
This last expression corresponds to the spatial discretization of the variational formulation. For the finite element method, the domain Ω can be decomposed into a set of triangles or quadrilaterals (in the 2D case) called finite elements.

The spectral element method used in SPECFEM keeps the same variational formulation (equation (3.51)). The displacement u and the velocity u are continuous everywhere in the domain (in the spectral elements and at the interfaces between elements) while the stress tensor σ is only continuous in the spectral elements (and discontinuous at the interfaces). As in finite elements, each integral will be decomposed into a sum of integrals per subdomain. The elements used for the mesh are quadrangles. Between each reference element, transfer functions are established. The two-dimensional reference elements are squares (Figure 3.3). Furthermore shape functions N k are associated with the control nodes. In 2D, the shape functions are bilinear for four nodes (the vertices of the element) and bicubic for nine nodes: for more details, see [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF]]. For the time discretization, a second order explicit Newmark finite difference scheme [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF] is used.

b) Finite differences

For a given variable f , the second-order discretization of the space derivative along x is

∂f ∂x | i+ 1 2 ,j = f i+1,j -f i,j ∆x , ( 3.56)
and the fourth-order discretization of the space derivative along x is

∂f ∂x | i+ 1 2 ,j = -f i+2,j + 27f i+1,j -27f i,j + f i-1,j 24∆x . (3.57)
Similar operators are applied along the z-direction. These operators are used to compute space derivatives of v x , v z , σ xx , σ zz and σ xz .

The fourth-order discretization scheme in space increases accuracy and allows us to use grids of reasonable size for sources with relatively higher frequency content when compared to second-order space discretization. For time integration we use a second-order leapfrog scheme [START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF]. In figure 3.4 we show the elementary grid cell of a 2D staggered grid finite difference discretization for the elastic wave equations.

Boundary conditions a) Free surface

To capture surface waves, i.e Rayleigh waves, a free surface must be implemented through Neumann-type conditions. In the case of the spectral elements, it is naturally implemented thanks to the variational formulation. Nevertheless, in the case of finite differences and finite volumes, it can be implemented by considering anti-symmetric conditions or non-centered space derivatives at the surface. This corresponds to a zero normal traction on the surface. In 2D elastic wave equations, it can be expressed as follows:

σ.⃗ n = 0, (3.58)
where ⃗ n is the external normal vector at the free surface.

This corresponds to σ zz = σ xz = 0 if the free surface is considered as the upper edge of the computational domain. Numerically, the values of σ zz and σ xz above and below the surface are opposite, and are equal to zero exactly on the free surface.

In order to implement the free surface for all x at the N Z level (points), we need to add ghost cells. For a FD staggered grid, the conditions at the free surface are:

σ xz (:, N Z) = 0 σ xz (:, N Z + 1) = -σ xz (:, N Z -1) σ xz (:, N Z + 2) = -σ xz (:, N Z -2) σ zz (:, N Z + 1) = -σ zz (:, N Z) σ zz (:, N Z + 2) = -σ zz (:, N Z -1). (3.59)
For the calculation of fluxes at any interface of the FVM, we need two cells on either side of the interface. Therefore, if the mesh counts N Z + 2 cells in the z direction, the free surface is placed at the N Z level and the free surface conditions are:

σ xz (:, N Z) = 0 σ xz (:, N Z + 1) = -σ xz (:, N Z -1) σ xz (:, N Z + 2) = -σ xz (:, N Z -2) σ zz (:, N Z) = 0 σ zz (:, N Z + 1) = -σ zz (:, N Z -1) σ zz (:, N Z + 2) = -σ zz (:, N Z -2) (3.60)
To compute accurately the stress variable σ xx on the free surface, we need to compute the spatial derivative of v z in the z-direction. To this end, we use the fourth-order non-centered operator proposed by [START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF]. This is done in the second sweep of the finite volume algorithm where the solutions are computed in the z-direction. Thus we get :

∂v z ∂z (:, N Z) ≈ a 4 v z (:, N Z + 1) + a 3 v z (:, N Z) + a 2 v z (:, N Z -1) + a 1 v z (:, N Z -2) + a 0 v z (:, N Z -3) ∆z ,
where a 0 = 0, 04168, a 1 = -0, 3334, a 2 = 1, 233, a 3 = -0, 6326 and a 4 = -0, 30874.

b) Absorbing conditions

To prevent reflections at the outer-boundaries of the domain, absorbing conditions are considered. Paraxial absorbing conditions like those of [START_REF] Clayton | Absorbing boundary conditions for acoustic and elastic wave equations[END_REF] (first order) or that of [START_REF] Stacey | Improved transparent boundary formulations for the elastic wave equation[END_REF] (second order) which is an improvement of the previous one are generally used. The latter essentially absorb the waves arriving perpendicular to the edge. However, it is not very satisfactory. High-order paraxial methods [START_REF] Bécache | Stability of Perfectly Matched Layers, group velocities and anisotropic waves[END_REF]) or PMLs (Convolutional Perfectly Matched Layer C-PMLs or Auxiliary Differential Equations ADE-PMLs) for example are preferred. C-PML formulations have been developed for different media to optimize the absorption of the waves impinging the outer PML layers at grazing incidence. For more details we refer to [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF][START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF][START_REF] Martin | An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave equation in poroelastic media[END_REF][START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF].

The non-convolutional ADE-PML [START_REF] Martin | A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML)[END_REF] has also been applied and is a more flexible generalization of the CPML formulation. Indeed, it is applicable to higher order timestepping schemes more easily. We adapted the CPMLs and ADE-PMLs to the finite volume technique by applying the stretching functions and the memory variables directly to the spatial derivatives of the numerical fluxes (see equations (3.36) and (3.37)) instead of the derivatives of the particle velocities and stresses. This formulation is more well suited and adapted for the finite volume method when compared to the classical CPML and ADE-PML formulations, because our formulation is able to deal with the stretching of the numerical fluxes derivatives and their related PML memory variables that are both based on left and right Riemann-based state vectors and variables. Besides, a null flux at the outer grid points is assigned instead of classical Dirichlet conditions. In a previous section 2.5.3, we summarized the ADE-PML and C-PML formulations for the 2D elastic wave equations. The Stacey conditions are also mentioned in section 2.5.4. In all the following numerical tests, we only show the simulated solutions with CPML conditions because ADE-PML as well as hybrid Stacey/C-PML (or ADE-PML) are not bringing substantial improvements of the solutions.

Stability condition

At the second order in time the CFL stability condition for FD is verified when

CF L = V p ∆t 1 ∆x 2 + 1 ∆z 2 ≤ 0.6,
where ∆t and ∆x, ∆z are respectively the time step and the space steps that must be carefully chosen [START_REF] Virieux | P-SV wave propagation in heterogeneous media: velocity-stress finitedifference method[END_REF], [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]). For spectral element method, the CFL must be less than 0.68 [START_REF] Komatitsch | Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène (Spectral and spectral-element methods for the 2D and 3D elastodynamics equations in heterogeneous media)[END_REF]. And for the split FVM, the CF L condition can be slightly greater than 1 [LeVeque, 2002a].

Numerical tests

The finite volume method and the boundary conditions i.e, free surface and absorbing conditions, are validated via five test cases:

-Model A: In this case, we consider a checkerboard-like heterogeneous medium with Convolutional Perfectly Matched Layer (CPMLs) applied at all boundaries, with a source and 12 receivers placed in the middle. Two cases are considered: the source is excited in the x-direction or in the z-direction. Through this test, the CPMLs conditions are also validated.

-Model B: In this case, we consider the same case as case A but this time we place the source and the receivers on the free surface. Absorbing conditions are defined at the other boundaries. This test aims at validating both the free surface conditions and the PML conditions in presence of a free surface and the related surface waves.

-Model C: In this case, we consider a fluid-solid model with water at the top of the medium and two solid layers below. This test aims at dealing with a physical case of interest for the hydro-geophysical community. Another interesting point is to test the stability of the PMLs conditions for acoustic-elastic coupling simulations.

-Realistic models: We test two cases in the context of unconsolidated granular media: a multi layer discontinuous model and a model varying continuously with depth. With these two models, we aim at modeling unconsolidated granular media at the laboratory scale. The second model was determined by ray tracing and phase diagram inversion [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]). It allows us to compare the solutions of our code with synthetic data and compare our results (time analysis) with those obtained in 2014.

In the sequel, FVM solutions are compared with those of a finite difference code (Seismic_-CPML) and a spectral element code (SPECFEM). This validation is done in three steps:

-qualitative analysis: comparison of seismograms with Seismic_CPML and SPECFEM; -quantitative analysis: calculation of the error according to a reference solution (solutions of the spectral elements);

-energy curves of the system: it allows us to mainly assess the efficiency and stability of the absorbing conditions.

In total we have tested 4 flux limiters: MC, VANLEER, KOREN and SMART. By comparing the results of different limiters with spectral element solutions, KOREN gave us the best results. In the rest of our study we chose this limiter.

Setup a) Model A

We consider a checkerboard model composed of two materials (see Figure 3.5a) which have the following properties:

• material 1 : ρ = 2800 kg m -3 , V p = 3300 m s -1 , V s = 1905.31 m s -1 ; • material 2 : ρ = 2800 kg m -3 , V p = 2600 m s -1 , V s = 1501.15 m s -1 .
The studied physical model has a width of 105 m and a height of 67.5 m. The size of each block making up the checkerboard is 21 m in length and 20.875 m in width. The source is the first derivative in time of a Gaussian with a dominant frequency f 0 = 200 Hz and a time delay

t 0 = 1.2 f 0 s.
It is given by:

s x = -2A sin(θ)a(t -t 0 )exp(-a(t -t 0 ) 2 ), s z = -2A cos(θ)a(t -t 0 )exp(-a(t -t 0 ) 2 ),
where a = π 2 f 2 0 , A = 10 7 is the amplitude and θ is the force angle relative to the vertical direction.

The source position is (x s , z s ) = (52.37, 31.125) m at the center of the model and we consider two cases:

• θ = 0 o which corresponds to an excitation of the source in the z-direction;

• θ = 90 o which corresponds to an excitation of the source in the x-direction.

The total number of receivers is 12 and the receivers line ranges from (10, 31.125) to (95, 31.125) m (see Figure 3.5a).

CPML conditions [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF] are applied at all boundaries of the model to mimick an infinite medium. Each layer is made up of 15 points (control cells, points, ..).

The time step used is ∆t = 0.03 ms and the total number of time steps is 2000 corresponding to a total of 60 ms. We used the same time step for the different numerical methods.

The Finite Volume method is 2-nd order in space while finite differences and spectral elements are 4-th order. For this, we will consider twice as many cells for finite volumes as for the other methods.

For the Finite Volume space discretization, we consider an uniform mesh such that ∆X = ∆Z = 0.125 m and composed of 841 × 501 finite volumes. The total number of control cells of the mesh is 421341. In this case, the CF L condition stability can be larger than 1 and is here equal to 1.1 due to the time split algorithm.

For the Finite Difference discretization, an uniform mesh is considered such that ∆X = ∆Y = 0.25 m and is composed of 421 × 251 points. The total number of points of the mesh is 131524. In this case, the CF L number has to be smaller than 1 and is equal to 0.53.

For the Spectral Element discretization, we consider a mesh composed of 63×105 elements which corresponds to 250 × 420 points (each element is composed of 4 points in each direction). The total number of points of the mesh is 105000. In this case, the CF L number has to be smaller than 1 and is also equal to 0.53.

b) Model B

We consider exactly the same physical configuration (i.e the checkerboard model), the same time and space discretizations as in the previous test case A for the different methods. We only change the position of the source and the receivers. They are placed at the free surface in order to test the free surface condition implemented there (see figure 3.5b). The source is excited at the surface only in the vertical direction. This time, instead of having absorbing conditions at the surface, we consider the free surface condition as a combination of an antisymmetric condition and a condition derived from the non-centered scheme given in subsection a).

c) Model C

We consider a model composed of water and two other materials (see figure 3.5c) which have the following properties:

• water : ρ = 1000 kg m -3 , V p = 1500 m s -1 , V s = 0 m s -1 • material 1 : ρ = 2500 kg m -3 , V p = 2950 m s -1 , V s = 1550 m s -1 • material 2 : ρ = 2800 kg m -3 , V p = 3500 m s -1 , V s = 2020.72 m s -1
The studied physical model has a width of 60 m and a height of 30 m. The size of each region is 10 m in depth and 60 m in length. The source is the first derivative in time of a Gaussian as defined in the two previous cases with a dominant frequency f 0 = 100 Hz, a time delay t 0 = 1.2 f 0 s and θ = 0 o , which corresponds to an excitation of the source in the z-direction. The source position is (x s , z s ) = (30, 15) m at the center of the model. The total number of receivers is 12 and the receivers line ranges from (5, 29) to (55, [START_REF] Luo | Rayleigh-wave dispersive energy imaging using a high-resolution linear radon transform[END_REF]) m (see figure 3.5c).

CPML conditions, where each PML layer is composed of 15 points, are applied at all boundaries of the model except at the top of the model where we consider a free surface condition.

For the Finite Volume space discretization, we consider an uniform mesh such that ∆X = ∆Z = 0.1 m and is composed of 600 × 300 finite volumes. The total number of control cells of the mesh is 180000. The time step used is ∆t = 10 -5 s and the total number of time steps is 20000 corresponding to 200 ms. In this case, the CF L number is equal to 0.49.

For the Finite Difference discretization, we consider an uniform mesh such that ∆X = ∆Y = 0.2 m and composed of 300 × 150 points. The total number of points of the mesh is 45000. The time step used is ∆t = 2.10 -5 s and the total number of time steps is 10000 corresponding to 200 ms. In this case, the CF L number is equal to 0.49.

For the Spectral Element discretization, we consider a mesh composed of 150×75 elements which corresponds to 600×300 points (each element is composed of 4 points). The total number of points of the mesh is 180000. The time step used is ∆t = 10 -5 s and the total number of time steps is 20000 corresponding to 200 ms. In this case, the CF L number is also equal to 0.49.

Results and discussion a) Model A:

In Figure 3.6, the snapshots of the v z component are shown and we can see that waves are well absorbed. By comparing the seismograms of v x and v z (see Figure 3.7) for the different methods in the case of Model A with a source excited in the z-direction, we notice that finite volumes and the spectral elements coincide. They have the same waveform and very similar amplitudes. On the other hand, the horizontal component of the velocity v x of the finite differences does not have the same waveform. This is due to the fact that a staggered grid is used in this case. Indeed, this component of the speed is not calculated at the same point as that of the two other methods due to half spacing shift in the staggered grids. Hence the interest of the finite volume method compared to the finite differences method which uses a staggered grid. When the solution is excited in the horizontal direction, the problem just rotates 90 degrees. We obtain similar results except that this time it is the vertical component of the speed v x which is not calculated at the same point as that in the two other methods. This explains the difference in the shape of the waveform in the case of finite differences with respect to spectral elements and finite volumes (see Figures 1, 2 and 3 in appendix).

The relative errors of the finite volume solutions as well as those of the finite differences presented in Figure 3.8 are less than 3 % on v z and 0.2 % on v x when compared to SPECFEM solutions. Those of finite volumes (FVM) are much smaller compared to those of finite differences (FD), especially when receivers are close to the source.

The total energy presented in Figure 3.9 decreases towards the end of the simulation, which means that the implemented CPML absorbing conditions are stable. We notice that in the case of FVM, we lose more energy than with finite differences which shows that the CPMLs are even more efficient for FVM than for staggered FD.

b) Model B

In Figure 3.10, the snapshots of the v z component are shown and we can see that waves are well absorbed. By comparing the seismograms of v x and v z (see figure 3.11) for the different methods in the case of Model B with a source excited in the z-direction, we notice that they have the same waveform and a very similar amplitude.

Besides, the relative errors of the FVM are varying between 0.14% and 0.3% for v x and between 0.1% and 0.4% for v z . For the FD they are varying between 0.14% and 0.4 % for v x and between 0.25% and 0.85% for v z (see Figure 3.12). As we can see errors of FVM are smaller than those of the FD method.

Finally, we notice that the energy displayed in Figure 3.13 decreases by several orders of magnitudes towards the end of the simulation, which means that the CPML absorbing conditions are stable too. We notice that in the case of the FVM, we lose also more energy than with FD method.

c) Model C

We can see through the snapshots of v x component represented in Figure 3.14 that waves are well absorbed. By comparing the seismograms of v x and v z (see figure 3.15) for the different methods in the case of Model C with a source excited in the z-direction, we notice that they are practically the same. Indeed, this is what the errors in Figure 3.16 confirm: they are smaller than 0.05% for v x and 0.09% for v z compared to spectral element solutions. Thus, very good qualitative and quantitative results are obtained in this case which involves water in the upper layer of the medium.

Realistic case study in unconsolidated granular media context

Now, at the laboratory scale, we consider a power law model deduced experimentally by ray tracing and phase diagram inversion and validated numerically with purely elastic wave equations for unconsolidated granular media in the case of no fluidization (i.e over-pressure equal to zero) [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]). The model is given by (see Figure 3.17 ):

V p = γ p (ρgh) αp (3.61) V s = γ s (ρgh) αs (3.62)
where ρ, g and h are respectively the density, gravitational acceleration on Earth surface equal to 9.81 m.s -2 and the depth of the media. γ and α are parameters estimated experimentally and given as follows:

α p = 0.3, γ p = 21, α s = 0.33 and γ s = 8.2. We test this model with our finite volume code. We consider first a 10-layer model extracted from the power-law given above. Then, in a second case, we integrate the whole continuous model into our finite volume code and check if we are able to retrieve first arrivals (P-waves) similar to the analytical and finite difference arrival times.

10-layer model

We consider a 10-layer model of 1 m in length and 0.2 m in depth which follows the power law given in equations (4.2a) and (4.2b). Each layer is 0.02 m thick and laterally homogeneous. For all the materials we will consider a density ρ = 1610 kg m -3 . Thus for the layers at the top and at the bottom of the computational domain, the materials are given by:

• seismic velocities close to the free surface :

-V p = 21 × (1610 × 9.81 × 0.02) 0.3 m s -1 = 118.052 m s -1 , -V s = 8.2 × (1610 × 9.81 × 0.02) 0.33 m s -1 = 54.784 m s -1 .
• seismic velocities at the bottom:

-V p = 21 × (1610 × 9.81 × 0.2) 0.3 m s -1 = 235.546 m s -1 , -V s = 8.2 × (1610 × 9.81 × 0.2) 0.33 m s -1 = 117.127 m s -1 .
The source wavelet is the same as in the previous cases A, B and C (but with a dominant frequency equal to 1200 Hz) and is located inside the medium at (x s , z s ) = (0.25, 0.175) m. The free surface is implemented at the top of the computational domain and absorbing conditions composed of 15 points (C-PML layer) are imposed on both vertical sides and at the bottom of the model (see Figure 3.18). 100 receivers are placed at the free surface and spaced 10 mm every 10 grid points.

We consider the same time and space discretizations for finite volume and finite difference methods: 1 × 1 mm 2 cells for a mesh of 1000 × 200 points and a time step ∆t = 1µs which corresponds to a CFL = 0.49 < 1. We run our simulation for 20000 time steps which correspond to 20 ms (physical time).

For the spectral element method, we consider 250 × 50 elements and a time step ∆t = 1µs which correspond to a spacing of ∆x = ∆z = 1 mm and a CF L = 0.49 < 1. In this case, to reach 20 ms (physical time), we run also the simulation for 20000 time steps.

In Figure 3.19, we can see that the waves that arrive on the boundaries where PML are implemented are well absorbed. In Figure 3.20, we notice that P-waves and surface waves (Rayleigh waves) are very well computed with the finite volume method compared to the reference solution of SPECFEM. Indeed, this is what the Figure 3.21 confirms: errors are of the order of 0.1% for both FVM and FD methods but finite volume errors are smaller than those of the finite differences. In Figure 3.22, the total energy curve for the finite volume method shows that the PMLs are stable and that they absorb waves better than the finite difference method.

Power-law model

We consider a model of 1 m in length and 0.215 m in depth. The simulated model is laterally homogeneous, and its velocities follow the power-law trend with depth for the case of zero over-pressure (see Figure 3.17). It is discretised by 0.5 × 0.5 mm 2 elements for a mesh of 2000×430 points for both the finite volume and finite difference methods. We consider the same source as those of test cases A, B and C but with a dominant frequency equal to 1500 Hz. The source is located at the free surface at (x s , z s ) = (0.25, 0.214) m. As in the previous 10-layer Model test case, the free surface is implemented at the top of the computational domain and absorbing conditions composed of 15 points (C-PML layer) are imposed on both vertical sides and at the bottom of the model (see Figure 3.18). The receivers are spaced each 5 mm every 10 grid points at the free surface. Stability and dispersion conditions are verified. Indeed, the time step ∆t = 1µs corresponds to a CF L = 0.6 < 1. We run our simulation for 20000 time steps which correspond to 20 ms (physical time).

In Figure 3.23, we can see clearly the physically dispersive nature of the wave patterns in this kind of unconsolidated compact medium, and the surface waves are particularly well computed as well as the body waves. SPECFEM solutions are not shown here because it is difficult and non trivial to implement the boundary conditions at the surface where seismic velocities are tending towards zero. In the FVM, the seismic velocities at the surface are taken equal to those defined just one grid point below the free surface. Obviously, by construction, these points below the surface are not located at the same place for FVM and SPECFEM due to the non-equidistant spatial distribution of the Gauss-Lobatto-Legendre points of SPECFEM.

Conclusion

In this chapter, we present a RFV-FSP finite volume method adapted to near surface problems. Flux-based CPML-type absorbing conditions are implemented for this method: they are stable and work very well in both free surface and no free surface cases. We succeeded in calculating the surface waves and more particularly the Rayleigh waves by applying an antisymmetry condition or a 4-th order non-centered operator for the stress σ xx to calculate the solutions more accurately close to the free surface. These boundary conditions show better qualitative and quantitative results compared to the finite difference method. Fluid-solid simulations provided also very good results. On the other hand, we succeeded with this RFV-FSP tool to simulate steep gradients of the seismic properties in depth, in the context of more complex models such as heterogeneous and dispersive unconsolidated granular media. Very good results are obtained with RFV-FSP in the "discontinuously" varying model with depth, and comparisons with SPECFEM are very good. Besides, it is also able to reproduce the physical dispersion of waves that are travelling through a sharp seismic profile varying continuously with depth. This will pave the way to future numerical developments for RFV-FSP and spectral element methods to simulate more complex and heterogeneous media in geophysical applications for non trivial free surface conditions and where strong gradients with depth are present.

With this finite volume method it will be possible to deform the mesh which will allow us to deal with realistic problems and to consider non-flat topographies (which is not quite easy with staggered grid finite differences). Besides, we are interested in high dominant frequency sources which implies a fairly fine mesh. However, due to the resolution of a Riemann problem at each interface, the computational time is greater with the finite volumes than with the finite differences. Thus, to deal with even larger problems, we shall consider accelerating the finite volume code by parallelizing it and/or by introducing a mesh unrefinement strategy. We shall also improve its precision by applying high-order Runge-Kutta time stepping methods.

Chapter 4

Modeling of unconsolidated granular media 4.1 Introduction

The near surface is defined as the first hundreds of meters from the ground. It is a particularly complex environment, which may contain natural resources such as water, gas, hydrocarbons ... The near surface can also be considered as the interface between the upper crust and the atmosphere and is the place of fluid exchanges and human and life activities. The weather and human actions can modify the structure and composition of this zone over time. It then becomes necessary to develop tools that will allow us to follow its evolution. Various geophysical methods and techniques have been developed and implemented such as electrical, electromagnetic, gravimetric and seismic methods. The choice of one of these methods depends on the type of physical properties to be studied, the object and/or area to be characterized and finally the environment in which it is located.

A natural medium can be composed of unconsolidated materials such as clay, silt, sand, gravel or materials derived from the erosion of rocks. It is thus important to characterize media made of these materials to better understand acoustic and seismic measurements in soils and their related geological environments. Indeed, many studies of elastic wave velocities in such materials have been conducted in situ, and on analogue experiments at the laboratory scale where theoretical models have been developed. Among the seismic techniques used to determine the physical properties of such a medium, we can cite the tomography techniques which is based on minimising the first arrivals [START_REF] Improta | Seismic imaging of complex structures by non-linear traveltime inversion of dense wide-angle data : application to a thrust belt[END_REF][START_REF] Le Meur | Tomographie tridimensionnelle à partir des temps des premières arrivés des ondes P et S. Application á la rǵion de Patras (Grèce)[END_REF][START_REF] Ravaut | Tomographie sismique haute résolution de la croûte terrestre : inversion combiné des temps de trajet et des formes d'ondes de données sismiques réflexion/réfraction grand angle multitraces[END_REF][START_REF] Zelt | Seismic traveltime inversion for 2d crustal velocity structure[END_REF] or the MASW (Multiple Analysis of Surface Waves) technique [START_REF] Bitri | Caractérisation du proche sous-sol le long de tracés linéaires par profilage sasw[END_REF]Foti, 2000a;[START_REF] Ganji | Detection of underground obsacles by sasw method -numerical aspects[END_REF]Nazarian and Stokoe, 1984a;[START_REF] Park | Predictions of shear-wave velocities in southern california using surface geology[END_REF][START_REF] Xia | Estimation of near-surface shear-wave velocities and quality factors using multichannel analysis of surface-wave methods[END_REF][START_REF] Xia | Estimation of near-surface shear-wave velocity by inversion of rayleigh waves[END_REF] which studies the surface wave dispersion. Toolboxes like SWIP (Surface-wave inversion and profiling) method, presented in [START_REF] Pasquet | Swip: An integrated workflow for surface-wave dispersion inversion and profiling[END_REF], can be used to estimate the P-pressure and S-shear wave velocities (V p and V s , respectively) of such media. Some common resulting models proposed at the near surface scale show an important increase in P and S wave velocity as a function of depth. Indeed, based on Hertz-Mindlin contact theory (in the context of intergrain forces modeling), the velocity structure of such medium can be modelled as [START_REF] Gassmann | Elastic waves through a packing of spheres[END_REF]:

V p,s = γ p,s (ρgz) αp,s (4.1)
where g is the gravity acceleration, γ p,s is a depth-independent coefficient mainly depending on elastic properties of grains, porosity and coordination number of the packed structure, and α p,s is the power-law exponent. The parameter α p,s depends on the dispersion and the form of grains [START_REF] Makse | Why effective medium theory fails in granular materials[END_REF][START_REF] Schön | Physical properties of rocks: Fundamentals and principles of petrophysics[END_REF][START_REF] Tournat | Acoustics of unconsolidated "model" granular media: An overview of recent results and several open problems[END_REF][START_REF] Zimmer | Seismic velocities of unconsolidated sands: Part 1-pressure trends from 0.1 to 20 mpa[END_REF].

In this chapter, we want to model numerically the seismic wave propagation in unconsolidated granular media. In particular, we revisit a study done in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] on these media at the laboratory scale. We consider the case where the media does not contain a fluid. Indeed, before characterizing the (partial or full) saturation of these media with fluid (water, gas), it is necessary to correctly model the dry media. The physical model deduced is a powerlaw model with α p,s close to 1 3 , values that have been found to be greater than the value of 1 6

proposed by [START_REF] Bachrach | High-resolution shallow-seismic experiments in sand, part ii: Velocities in shallow unconsolidated sand[END_REF]] for a shallow sand medium in situ. We want here to show that our power-law reproduces much better the seismograms and phase dispersion diagrams for the granular medium under study at the laboratory scale. On the medium/long term, all the methodology used here to define the power laws can be applied at the laboratory scale, and thus also to data collected at near-surface field (in-situ) scales.

The main objective of this study is to numerically validate the models considered by verifying whether they allow us to better interpret the data recorded in the laboratory. In order to do this, simulations are run to fully replicate the configuration of the laboratory experiments. The physical problem consists in solving the seismic wave equation for different rheologies (pure elastic, porous) in order to determine which of them best explains the laboratory data.

First, 3D simulations are performed in the case of a pure elastic medium. The calculated solutions are then compared with the experimental data recorded in the laboratory. A time and dispersion analysis are done to compare respectively the first arrivals and the dispersion curves of the experimental data with the simulated ones.

Second, since the physical model considered is horizontally stratified, and the source and receivers belong to the same plane, we can reduce the 3D model to a 2D one. The 2D seismograms can reproduce the 3D one by introducing a source-receiver distance rescaling.

Third, we make the rheology more complex (porous medium for instance) in 2D to compare it to the homogenized elastic medium and to the experimental results.

The comparison between simulated and experimental data is done in two domains:

• Time domain: we detect the first travel-time arrivals of the seismic waves.

• Frequency domain: we calculate the dispersion curves using a slant-stack transform.

Finally, the results obtained are discussed and conclusions and recommendations are drawn. The choice of the numerical tools has been made taking into account four main criteria: precision, performance, simplicity of the method, and the nature of the physical problem to be solved. As the experiment is made with a flat topography, a finite difference method is a good choice. Indeed, in this case, it responds favorably to all the criteria mentioned above.

A 3D fourth-order finite difference code called UNISOLVER is used to solve the wave equations in a stratified elastic medium taking into account the different models of high gradient velocities. This code is parallelized and therefore the computation time is optimized which makes it very efficient for the realization of our simulations. Absorbing and free surface conditions are implemented and validated by comparing with reference codes such as SPECFEM (spectral finite elements). The accuracy of its schemes is fourth-order in space and second-order in time. 2D versions of this code for both elastic and porous media are also used to highlight the contribution of the porous medium to the elastic medium.

We also use the 2D finite volume code validated in Chapter 2 to compare the finite volume solutions with the finite difference ones.

We identified two main issues to model this kind of complex granular medium:

• The rheological law (elastic or poroelastic) describing the mechanical behavior of the medium.

• The modelling of the source as a stick as in the experiment [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] or as a virtual point source located at different depths.

The choice of the rheological law and the source model influences the waveforms and the amplitudes as well as the different surface and P modes in the dispersion diagrams. This will be discussed in the following sections.

Section 2 deals with experimental setup and physical model obtained by inversion of first arrivals and dispersion diagrams. In section 3, the numerical setup for the different rheological models is given. And in sections 4 and 5, a qualitative (seismograms, spectograms and dispersion images) and quantitative (first arrivals, dispersion curves) analysis is done on the different numerical models considered. Finally in section 6, we consider different locations of the numerically modelled sources.

Experimental setup and physical model

Experimental setup

We are interested in an experiment carried out in the laboratory on unconsolidated and granular/porous media [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]. The studied physical medium has a length of 1 m, a width of 0.8 m and a height of 0.22 m. The box is filled with 180-300 µm diameter glass beads (GB1) and the bottom of the box consists of a metallic sieve glued on a perforated plate (Figure 4.1).

The source consists of a metal stick of 0.015 m connected to a waveform generator (Figure 4.2) with a dominant frequency f 0 = 1500 Hz. The latter is injected in the vertical yOz-plane with a tilted angle of 20 o from the vertical normal to the free surface. The source position is (x s , y s , z s ) = (0.4, 0.25, 0.215) m at the free surface.

The bulk density of GB1 is given by ρ GB1 = 1610 Kg.m -3 and its corresponding porosity is Φ GB1 = 0.356.

For a given source location, the normal component of the particle velocity (V z ) is recorded in time at the surface of the medium as a "seismogram" using an oscilloscope. Up to 100 traces were recorded (using an oscilloscope) in linear single-channel walkway mode along the Oy direction.

The experimental data at our disposal correspond to 25 receivers placed linearly and equidistantly at the surface of the medium. From the source over a length of 0.5 m in the y-direction, the spacing between the receivers is 0.02 m. We filter the experimental data with a Butterworth filter in the 200-2700 Hz frequency band as shown in Figure 4.3. In Figure 4.4, we show the seismograms of V z component (filtered and not filtered) at receivers 10, 15 and 20 (located respectively at 20, 30 and 40 cm from the source). In the following sections we will compare solutions only with filtered real data.

Physical model

We are interested in a V p /V s model obtained from the experiment in the HOM22 medium in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]. The latter is deduced by ray tracing-based travel time inversion for the V p model and phase velocity diagram inversion for the V s model and is given by a power law as follows (Figure 4.1):

V p = γ p (ρgh) αp (4.2a) V s = γ s (ρgh) αs (4.2b)
where ρ, g and h are respectively the density, the gravitational acceleration at Earth surface (g = 9.81 m.s -2 ) and the depth of the medium. The parameters γ and α are estimated experimentally and given as follows: α p = 0.3, γ p = 21, α s = 0.33 and γ s = 8.2.

In [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF], this model was partially validated numerically using a 3D elastic finite difference code. But there was still a problem in reproducing the amplitudes of higher modes excepting the fundamental mode. First, we reproduce the results of [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] with UNISOLVER, a 3D finite difference code. Secondly, we reduce the 3D model to a 2D model and finally we integrate the V p /V s model into a 2D poro-elastic model according to the Biot elastic model [Biot, 1956a]. One of the main goals is to better reproduce the amplitudes of higher modes.

Numerical setup

We consider three models:

• 3D elastic: one uses UNISOLVER (a finite difference/FD code of fourth-order in space and second-order in time [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]) to reproduce the results of [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF].

• 2D elastic: we reduce the 3D model to the plane where the metal stick is located. In this case, one uses seismic_CPML_elastic (a FD code fourth-order in space and second-order in time [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF]) and we compare the results to the experimental data.

• 2D poro-elastic: poroelastic materials are most of the time modelled using the Biot theory [Biot, 1956a] and [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: Higher-frequency range[END_REF]. The latter takes into account parameters such as porosity, permeability, viscosity and tortuosity (measured during the experiment). The compressional P wave velocities (fast and slow) and shear S wave velocity in the porous medium depend also on the P and S elastic wave velocities of the solid frame. A version of Seismic_CPML for Biot model is used [START_REF] Martin | An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[END_REF].

For each model, we consider three cases:

• Dirichlet: in this case, all the boundaries except the free surface are Dirichlet boundary conditions. Dirichlet conditions are defined to simulate a rigid material on the five other sides of the box.

• PMLs: in this case, we consider absorbing boundary conditions (C-PMLs) on all the boundaries except the free surface and the bottom boundary where a Dirichlet condition is conserved. The PMLs allow to absorb the outgoing waves to mimic an open medium.

• Full PMLs: in this case, on all the boundaries except the free surface, absorbing boundary conditions (C-PMLs) are considered.

The free surface condition is implemented at the top of the computational domain using the zero normal traction assumption for the different cases.

For the different models and cases, the force source signal (Figure 4.2) is injected all along the metal stick. The latter is simulated numerically by injecting the time wavelet source signal over a series of points located along the stick.

3D/2D elastic models

The physical model has a width of 0.035 m, a length of 1 m and a height of 0.220 m. The density is given by ρ = 1610 Kg/m 3 . In 3D, we consider an uniform mesh such that ∆X = ∆Y = ∆Z = 5. 10 -4 m and composed of 70 × 2000 × 430 points. The total number of points of the mesh is 60, 2 millions points. The time step used is ∆t = 10 -6 s and the total number of time steps is 50000 corresponding to 50 ms (time of the experiment). In this case, the CFL number is equal to 0.6. PML layers are composed of 15 points. The simulated model is laterally homogeneous, and its velocities follow the power-law trend with depth (Figure 4.1). The velocity profiles, continuous in theory, are discretized according to the spacing discretization ∆Z in depth.

The receivers are spaced each 5 10 -3 m (10 grid points) at the free surface to record the seismograms of the normal component of the particle velocity (Figure 4.5).

We reduce the 3D model to a 2D model by considering only the plane where the metal stick (source) is located. The x-direction is omitted.

In order to improve the amplitude and the waveform of the recorded signals respect to the experimental data, we simulate different configurations using the 3D parallel FD code UNISOLVER : : with Dirichlet, PMLs (just at the four vertical outer boundary walls) or Full PMLs (everywhere) conditions except at the free surface. For the 2D case, we consider the same different configurations mentioned before and we also use the 2D code seismic_CPML_elastic.

The 4-th order 3D UNISOLVER parallel code has been scaled over different numbers of processors (from 100 up to 400) using MPI (Message-Passing-Interface) libraries. The computational domain has been cut along the longitudinal y-axis axis (from 20 points down to 5 grid points per processor). A buffer overlapping zone of two grid points between subdomains (one subdomain per processor) is used to communicate the particle velocities and stresses and material properties between processors via 'MPI -SEND' and 'MPI -RECV' libraries. As we can see in Figure 4.7, the strong scaling obtained by measuring the CPU time versus the number of processors is very satisfactory, even if a classical synchronous/blocking communication strategy has been introduced. This strategy is a classical one [START_REF] Komatitsch | An unsplit convolutional Perfectly Matched Layer improved at grazing incidence for the seismic wave equation[END_REF] even if other asynchronous strategies could have been used to reduce and hide communication times by overlapping communications by computations in the inner subdomains as depicted in [START_REF] Chau | MPI implementation of parallel subdomain methods for linear and nonlinear convection-diffusion problems[END_REF][START_REF] El Baz | Asynchronous iterations with flexible communication: contracting operators[END_REF], 2001;Martin et al., 2008a;[START_REF] Miellou | A new class of asynchronous iterative methods with order intervals[END_REF]]. The code has been run on different supercomputing machines such as Olympe of CALMIP computing centre of Toulouse (France), or IRENE of the TGCC platform of the French Nuclear Agency (Saclay/Paris/France).

2D poroelastic model

In the poro-elastic case, one considers the same discretization in space and time as in the elastic cases. The length and width of the box are also preserved.

In our simulation, the fluid considered is the air and the solid components are marbles. The density of the solid ρ s = 2500 Kg/m 3 and the density of the fluid ρ f = 1.20 Kg/m 3 give an apparent density ρ w = 6, 46 Kg/m 3 and the density of the satured medium ρ = 1610 Kg/m 3 (using relations in the Biot model). The porosity, the permeability, the viscosity and tortuosity are given by Φ GB = 0.356, K = 9.510 -11 m 2 , ν = 1.907 • 10 -5 and a = 1.91 respectively. The bulk modulus of air is 0.101 MPa. And the bulk moduli of the solid and saturated (porous) medium are computed depending on the variation of P-wave and S-wave velocities which affect the effective Lamé parameters λ and µ of the porous medium and λ s of the solid frame (see Figure 4.8 for the different velocity models computed with this model). V s and V p characterize the pure elastic model, while V s , V p fast and V p slow characterize the porous medium. We use seismic_CPML_poroelastic, a FD code of fourth-order in space and second-order in time to solve the 2D Biot model in the different boundary cases considered.

When we calculate f max = πf 0 and f c (see equation (2.39)), we obtain : f max = 4710 Hz < f c = 4959.592 Hz which means that we can consider a Biot model without attenuation.

Qualitative analysis

Seismograms and spectrograms

Figure 4.9 shows the experimental and numerical seismograms for the different models (3D elastic, 2D elastic and 2D poro-elastic) in the Dirichlet case. We can clearly see the P-wave followed by a train of guided waves, of the P-SV type. The P-wave, itself guided by the gradient of mechanical properties of the medium, corresponds to fast modes considered as essentially longitudinal. We can also see the reflections of the waves on the boundaries due to the Dirichlet conditions considered. However, in the PMLs and Full PMLs cases (Figure 4.10 and 4.11), these reflections disappear. In the spectrograms, in the case of Dirichlet conditions, the reflections contaminate the frequency signal (Figure 4.12). Imposing absorbing conditions (C-PML) on the boundaries absorbs the reflections and makes the spectrogram clearer (Figure 4. 13 and 4.14). By comparing the PMLs and Full PMLs cases, one notices that the Dirichlet conditions imposed at the bottom of the medium in the PMLs case creates interferences represented by horizontal lines which are much more visible on the white seismograms (Figure 4.10g and 4.10h). We will thus continue our analysis by considering the Full PMLs case only.

Dispersion images

As in [START_REF] Park | Predictions of shear-wave velocities in southern california using surface geology[END_REF]], a slant stack transform (an oblique summation of normalized signal amplitudes in the frequency domain) is applied to the experimental and simulated seismograms of Figure 4.11 and gives us the dispersion images shown in Figure 4.15. We identify three principal modes on the image of experimental dispersion over the source-relative frequency band [0.1 -3000Hz]. The maxima (in black/red in Figures 4.15) correspond to four modes: two low velocity P-SV modes (propagation modes) at Low Frequency (LF) in the 0.25-1 kHz frequency band (mode 0 and 1), and two P-modes (one principal: mode 2) at higher frequencies and higher velocities in the 1.25-2.5 kHz frequency band.

The dispersion images of the numerical elastic models (3D/2D) in Figure 4.15 show clearly only two modes (mode 0 and 1). The mode 2 appears more clearly in the poro-elastic model.

Quantitative analysis

First arrivals

The first-arrival (P-mode) time of the numerical data was picked at each trace and compared to both the theoretical (calculated for the unconsolidated granular medium, see Chapter 1 for more details) and experimental travel times. These travel times are shown on the time versus offset curve of Figure 4.16. Computed arrival times are matching experimental and theoretical results within a 5% error range.

Dispersion curves

The dispersion curves picked for the P-SV modes identified on the dispersion images (Figure 4.15) clearly fit experimental and theoretical dispersion (Figure 4.17). The P-SV theoretical dispersion curves are calculated using CPS/Computing Programs for Seismology (Herman et al 2013).

The P-mode picked on the 3D elastic and 2D elastic dispersion images appeared noisy due to its weak amplitude. But the P-mode picked in the 2D poro-elastic case is better represented compared to experimental data and theoretical dispersion. The theoretical dispersion is computed with the CPS code that uses the Thomson-Haskell matrix propagator technique and includes the complex-valued roots of the dispersion equation. Figure 4.17 show the peaks of the energy maxima corresponding to the identified modes on the dispersion images (Figure 4.15). This confirms that the modes identified on the numerical dispersion images correspond to the fundamental mode (mode 0) and to the second propagation mode (mode 1). The first P -mode (mode 2) is not visible enough on the dispersion image of the 3D/2D elastic model but it is very clear in the poro-elastic one.

We have been able to reproduce qualitatively the modes 0 and 1 but not the modes 2 related to P -modes. This is mainly due to the way the source is introduced numerically. Indeed, the way the source is modelled is crucial to detect all the main surface and P-SV modes (0, 1) as well as volumic wave modes (P modes 2 and higher), and to better reproduce numerically the waveforms of the experimental data. The next section is devoted to this aspect and tries to bring some answers about the source modelling.

Conclusion

Several things can be said at the end of this study. First, the V p -V s model has been well integrated in the elastic and poroelastic codes. Indeed the results obtained (from qualitative and quantitative points of view) are improved compared to those of [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] in the case of zero overpressure: numerical first arrivals and dispersion curves are very close to those of the experimental and theoretical ones. And the Finite Volume (RFV) code developed gives us also satisfying results. Secondly, the rheology considered (purely elastic or poroelastic) and the source location play an important role in the modelling of this kind of complex medium. Indeed, poroelastic models give better results and the deeper the source is placed the better the results are. The fundamental mode is less visible than the higher modes (such as the P-modes).

Chapter 5

Conclusions and perspectives (English) Numerical developments

In this work, we were interested in studying numerical models of unconsolidated granular media.

The first part was devoted to the numerical tools that we need for our study. We developed a finite volume method based on a Riemann solver that computes the full two-dimensional wave field in an elastic medium.

To model surface waves (Rayleigh wave), we used the image method to implement the free surface condition. But to achieve more accuracy, we implemented a non-centered fourth-order spatial scheme close to the free surface. This condition was validated by several synthetic test cases which consider a heterogeneous medium (with or without a fluid layer), and where the source and receivers are placed close to the flat topography.

For the modeling of an open domain, the C-PML absorbing conditions were adapted to the finite volume method. They are stable and the "outgoing" fields are well absorbed. The different validations proposed allow to consider much more complex problems (realistic case) to be solved. One of the main numerical developments of this work is the implementation of a Finite Volume method (RFV) as well as the adaptation of the C-PML absorbing boundary conditions and free surface conditions to the RFV technique. The PML conditions are applied to the spatial fluxes derivatives, which is a different formulation than classical PMLs that are generally applied to the spatial derivatives of the primitive variables (particle velocities and stresses). The validated tool can then be used for more realistic cases.

To lead time and dispersion analyses, we have developed tools that allow us to detect the travel time of the first arrivals, to calculate the dispersion curves and to identify the main modes of wave propagation for the different models considered. These tools were mainly used in Chapter 3 during the comparative study between the 3D elastic, 2D elastic and 2D poroelastic models.

111

Geophysical applications

In terms of applications, the idea was to address several spatial scales (from a few centimeters to tens of meters) and different frequency ranges (from 200 Hz to 1500 Hz). We tested the wave propagation on media with two and three materials including water. The calculated solutions compared to reference solutions (seismic_CPML and SPECFEM) were satisfactory from qualitative (snapshots) and quantitative (seismograms) point of views. Indeed, in most cases, finite volumes produce lower errors than finite differences. The surface waves are well calculated and the waves at the boundaries are well absorbed.

We revisited an experimental study made on unconsolidated granular media in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]. For this, we followed three steps:

• Reproduce the results obtained in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] using UNISOLVER, a 3D parallelized elastic code of finite differences, and the different processing tools developed to detect the first travel-time arrivals and to calculate the dispersion curves.

• Reduce the 3D problem to a 2D one which is possible because the model considered varies according to the depth only, and also because the source and receivers in the experiment belong to the same plane. Here we used seismic_CPML and our finite volume code.

• Add complexity to the model by considering a poroelastic model. This is done with a version of seismic_CPML for Biot model.

In this study, we used the absorbing conditions to remove spurious boundary reflections and therefore to clean up the signal. We conducted a time and dispersion analysis with the different models considered. The latter showed the superiority of the poroelastic model over the elastic one.

Then, we looked at the modeling of the source. We modeled it in two different ways: by a vibrating stick (series of points, as in the experiment) or by a point source. We did a sensitivity analysis on the position of the source. At the end of this study, we realized that the results were very sensitive to the position of the source: the more it was placed at depth (without exceeding a certain limit) the more satisfactory the results were. The obtained results show that the problem encountered in the study in [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] came essentially from the modeling of the source. At the end of this study we can say that the poroelastic model is more representative of unconsolidated granular media.

Perspectives

The 2D elastic finite volume code has been validated for heterogeneous media defined by isotropic and linear strain-stress mechanical laws. However, some complex media can be governed by non-linear laws. Indeed, non-linear dynamics are observed in granular media due to grain-to-grain interactions [START_REF] Martin | Seismic wave propagation in nonlinear viscoelastic media using the auxiliary differential equation method[END_REF][START_REF] Tournat | Acoustics of unconsolidated "model" granular media: An overview of recent results and several open problems[END_REF]]. In the future, we could implement non-linear laws in our finite volume code by modifying the Jacobian matrices introduced in the compact form of the equations [START_REF] Leveque | Finite-volume methods for non-linear elasticity in heterogeneous media[END_REF].

On the other hand, in some seismic data, we observe a dispersion in the seismic wave fields that is often associated with attenuation phenomena other than those arising from geometric attenuation [START_REF] Murphy | Effects of partial water saturation on attenuation in massilon sandstone and vycor porous glass[END_REF]. The dispersive phenomena are mainly due to the presence of geometrical waveguides from the interfaces between the different geological formations of the medium. One way to model this dispersion is by considering a viscoelastic rheology [START_REF] Dhemaied | Seismic-wave propagation modeling in viscoelastic media using the Auxiliary Differential Equation method[END_REF]. One can also ask ourselves if the numerical scheme considered can allow us to solve other rheologies (poroelastic, anisotropic or viscoelastic media) as well as more general elastodynamic equations (with nonlinearities ...).

One of the disadvantages of finite volume method compared to finite element methods is the handling of the free surface. Indeed, in finite element methods, free surface conditions are naturally considered in the variational formulation of the problem. This is why, it could be interesting to discuss more the need or not for fine cells around the free surface. Also, it would be interesting to compare the method of images implemented with other methods that can be used for collocated grids.

In this thesis, we have treated only cases with a flat topography. Knowing that the finite volume method is well suited for the deformation of the mesh (due to the fact that all variables waves are considered at the same grid point), a non-flat topography can be considered. And that was originally one of the principal motivations for which we chose this method, which is a good compromise between finite differences and finite elements methods. Consequently, we will be able to deal with even more complex cases which would allow us, for example, to study site effects or fluid-solid interactions (bathymetry) which interest many geophysicists.

In one hand the spatial accuracy of the method can be improved by considering ENO and WENO schemes [START_REF] Zhang | Eno and weno schemes[END_REF] and in the other hand the time precision can be improved by implementing a fourth order Runge Kutta method. In addition to the accuracy of the results, it is necessary to take into account the speed of calculation which constitutes a real challenge in the geophysical community. This is why it is necessary to consider a full parallelization of the code and introduce it in the UNISOLVER code framework.

The extension of 2D modeling to the 3D case is essential for both large and laboratory scales. As pointed out in the Chapter 3, wave interference issues could be removed or attenuated by extending 2D simulations to 3D ones. For example, the 3D effects for the site effect can not be neglected, and simple corrections of the amplitudes can not resolve the problem.

With regard to the numerical study carried out on unconsolidated granular media, once the model has been validated for dry media, we can now embark on a study on media partially saturated with fluid. A more complete study could be considered with a 3D poroelastic model. This will pave the way to study gas or fluid/water detection and also to monitor some areas for CO2 storage, mitigation of flash floods, etc..

Conclusions and perspectives (Français)

Développements numériques

Dans ce travail, nous nous sommes intéressés à l'étude des modéles numériques de milieux granulaires non consolidés.

La première partie a été consacrée aux outils numériques dont nous avons besoin pour notre étude. Nous avons développé une méthode de volumes finis basée sur un solveur de Riemann qui calcule le champ d'onde bidimensionnel complet dans un milieu élastique.

Pour modéliser les ondes de surface (onde de Rayleigh), nous avons utilisé la méthode des images pour mettre en oeuvre la condition de surface libre. Mais pour obtenir plus de précision, nous avons implémenté un schéma spatial décentré d'ordre 4 à proximité de la surface libre. Cette condition a été validée par plusieurs cas tests synthétiques qui considèrent des milieux hétérogènes (avec ou sans couche fluide), et où la source et les récepteurs sont placés près de la topographie plane.

Pour la modélisation d'un domaine ouvert, les conditions absorbantes de type C-PML ont été adaptées à la méthode des volumes finis. Elles sont stables et les champs "sortants" sont bien absorbés. Les différentes validations proposées permettent de considérer des problèmes beaucoup plus complexes (cas réaliste) à résoudre.

Un des principaux développements numériques de ce travail est la mise en oeuvre d'une méthode des volumes finis (RFV) ainsi que l'adaptation des conditions aux limites absorbantes C-PML et des conditions de surface libre à la technique RFV. Les conditions PML sont appliquées aux dérivées spatiales des flux, ce qui est une formulation différente des PML classiques qui sont généralement appliquées aux dérivées spatiales des variables primitives (vitesses et contraintes des particules). L'outil validé peut alors être utilisé pour des cas plus réalistes.

Pour mener des analyses en temps et de dispersion, nous avons développé des outils qui nous permettent de détecter le temps de trajet des premières arrivées, de calculer les courbes de dispersion et d'identifier les principaux modes de propagation des ondes pour les différents modèles considérés. Ceux-ci ont été utilisé principalement dans le Chapitre 3 lors de l'étude comparative entre les modèles 3D élastique, 2D élastique et 2D poroélastique.

Applications géophysiques

En termes d'applications, l'idée était d'aborder plusieurs échelles spatiales (de quelques centimètres à plusieurs dizaines de mètres) et différentes gammes de fréquences (de 200 Hz à 1500 Hz). Nous avons testé la propagation des ondes sur des milieux à deux et trois matériaux dont l'eau. Les solutions calculées comparées aux solutions de référence (seismic_CPML et SPECFEM) étaient satisfaisantes d'un point de vue qualitatif (snapshots) et quantitatif (sismogrammes). En effet, dans la plupart des cas, les volumes finis produisent des erreurs plus faibles que les différences finies. Les ondes de surface sont bien calculées et les ondes aux frontières sont bien absorbées.

Nous avons revisité une étude expérimentale faite sur des milieux granulaires non consolidés dans [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF]. Pour cela, nous avons suivi trois étapes :

• Reproduire les résultats obtenus dans [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] en utilisant UNISOLVER, un code élastique 3D parallélisé de différences finies, et les différents outils de traitement développés pour détecter le temps des premières arrivées des ondes sismiques et pour calculer les courbes de dispersion.

• Réduire le problème 3D à un problème 2D ce qui est possible car le modèle considéré ne varie qu'en fonction de la profondeur, et aussi parce que la source et les récepteurs dans l'expérience appartiennent au même plan. Ici, nous avons utilisé le code seismic_CPML et notre code de volumes finis.

• Ajouter de la complexité au modèle en considérant un modèle poroélastique. Ceci est fait avec une version de seismic_CPML pour le modèle de Biot.

Dans cette étude, nous avons utilisé les conditions absorbantes pour éliminer les réflexions parasites aux limites et donc pour nettoyer le signal. Nous avons effectué une analyse temporelle et de dispersion avec les différents modèles considérés. Cette dernière a montré la supériorité du modèle poroélastique sur le modèle élastique.

Ensuite, nous nous sommes intéressés à la modélisation de la source. Nous l'avons modélisée de deux manières différentes : par un pot vibrant (série de points, comme dans l'expérience) ou par un point source. Nous avons fait une analyse de sensibilité sur la position de la source. A l'issue de cette étude, nous nous sommes rendus compte que les résultats étaient très sensibles à la position de la source : plus elle était placée en profondeur (sans dépasser une certaine limite), plus les résultats étaient satisfaisants. Les résultats obtenus montrent que le problème rencontré dans l'étude de [START_REF] Bodet | Smallscale physical modeling of seismic-wave propagation using unconsolidated granular media[END_REF] venait essentiellement de la modélisation de la source. A l'issue de cette étude, nous pouvons dire que le modèle poroélastique est plus représentatif des milieux granulaires non consolidés.

Perspectives

Le code de volumes finis (RFV) a été validé pour des milieux hétérogènes définis par des lois mécaniques isotropes et linéaires de contrainte-déformation. Cependant, certains milieux complexes peuvent être régis par des lois non linéaires. En effet, des dynamiques non linéaires sont observées dans les milieux granulaires en raison des interactions grain à grain [START_REF] Martin | Seismic wave propagation in nonlinear viscoelastic media using the auxiliary differential equation method[END_REF][START_REF] Tournat | Acoustics of unconsolidated "model" granular media: An overview of recent results and several open problems[END_REF]. Dans le futur, nous pourrions implémenter des lois non-linéaires dans notre code de volumes finis en modifiant les matrices Jacobiennes introduites dans la forme compacte des équations [START_REF] Leveque | Finite-volume methods for non-linear elasticity in heterogeneous media[END_REF]. D'autre part, dans certaines données sismiques, nous observons une dispersion des champs d'ondes sismiques qui est souvent associée à des phénomènes d'atténuation autres que ceux issus de l'atténuation géométrique [START_REF] Murphy | Effects of partial water saturation on attenuation in massilon sandstone and vycor porous glass[END_REF]. Les phénomènes de dispersion sont principalement dûs à la présence de guides d'ondes géométriques provenant des interfaces entre les différentes formations géologiques du milieu. Une des façons de modéliser cette dispersion est de considérer une rhéologie viscoélastique [START_REF] Dhemaied | Seismic-wave propagation modeling in viscoelastic media using the Auxiliary Differential Equation method[END_REF]. On peut également se demander si le schéma numérique considéré peut nous permettre de résoudre d'autres rhéologies (milieux poroélastiques, anisotropes ou viscoélastiques) ainsi que des équations élastodynamiques plus générales (avec des non-linéarités ...). L'un des inconvénients de la méthode des volumes finis par rapport à celle des éléments finis est le traitement de la surface libre. En effet, dans les méthodes par éléments finis, les conditions de surface libre sont naturellement considérées dans la formulation variationnelle du problème. C'est pourquoi, il pourrait être intéressant de discuter davantage de la nécessité ou non de considérer des cellules fines autour de la surface libre. De plus, il serait intéressant de comparer la méthode des images mise en oeuvre avec d'autres méthodes qui pourraient être utilisées pour les grilles colocalisées.

Dans cette thèse, nous avons traité uniquement des cas avec une topographie plane. Sachant que la méthode des volumes finis est bien adaptée à la déformation du maillage (dû au fait que toutes les variables sont considérées au même point de grille), des topographies déformées pourront être envisagées dans certains cas. C'était d'ailleurs à l'origine une des principales motivations pour lesquelles nous avons choisi cette méthode, qui est un bon compromis entre les méthodes des différences finies et des éléments finis. Par conséquent, nous pourrons traiter des cas encore plus complexes qui nous permettront, par exemple, d'étudier les effets de site ou les interactions fluide-solide (bathymétrie) qui intéressent de nombreux géophysiciens.

La précision spatiale de la méthode peut être améliorée en considérant les schémas ENO et WENO [START_REF] Zhang | Eno and weno schemes[END_REF] et la précision temporelle en implémentant une méthode Runge Kutta de quatrième ordre. En plus de la précision des résultats, il est nécessaire de prendre en compte la temps de calcul qui constitue un véritable défi dans la communauté géophysique. C'est pourquoi il est nécessaire d'envisager une parallélisation du code. L'extension de la modélisation 2D au cas 3D est essentielle à la fois pour les grandes échelles et pour les échelles en laboratoire. Comme indiqué dans le chapitre 3, les problèmes d'interférence des ondes pourraient être supprimés ou atténués en étendant les simulations 2D aux simulations 3D. Par exemple, les effets 3D dûs aux effets de site ne peuvent pas être négligés, et de simples corrections des amplitudes ne peuvent pas résoudre les problèmes.

En ce qui concerne l'étude numérique réalisée sur des milieux granulaires non consolidés, une fois le modèle validé pour des milieux secs, nous pouvons maintenant nous lancer dans une étude sur des milieux partiellement saturés en fluide. Une étude plus complète pourrait être envisagée avec un modèle poroélastique 3D. Cela ouvrira la voie à l'étude de la détection de gaz ou d'eau/fluides et aussi au "monitoring" de certaines zones pour le stockage de CO2, l'atténuation des crues soudaines, etc.

Submitted articles

⋆ Kassem Asfour, Roland Martin, Didier El Baz and Ludovic Bodet. Impact of source modeling and poroelastic models on numerical modeling of unconsolidated granular media : application at the laboratory scale. Submitted to Measurements. are respectively the physical state vector, the source vector and Jacobian matrices.

The propagation velocities of the elastic waves are determined by the eigenvalues ζ i of matrices A and B which are given by

ζ 1 = 0, ζ 2 = -v s , ζ 3 = v s , ζ 4 = -v p , ζ 5 = v p
, where v p = ( λ+2µ ρ ), v s = ( µ ρ ) are the compressional wave (P -wave) velocity and the shear wave (S-wave) velocity, respectively.

We rewrite A and B using their diagonalized expressions such as:

A = R x (R x ) -1 , B = R z (R z ) -1 , where R x =           

Figure 2

 2 Figure 2.1: P-waves (top) are polarized in the direction of propagation and S-waves (bottom) are slower than P-waves and are polarized perpendicular to the direction of propagation. Figure taken from https://www.sms-tsunami-warning.com/pages/seismic-waves#.YOw_6VMzbCI.

Figure 2 . 2 :

 22 Figure 2.2: Love waves (top) are horizontally polarized shear waves(SH waves) and Rayleigh waves (bottom) are characterized by its retrograde elliptical polarization in the propagation plane. https: //www.sms-tsunami-warning.com/pages/seismic-waves#.YOw_6VMzbCI.

Figure 2 . 3 :

 23 Figure 2.3: Types of waves: (a) P-waves traveling in a isotropic elastic solid medium; (b) the vertical component of S-waves traveling in a solid medium; (c) SH-waves (horizontal component of S-waves) traveling in a solid medium (d) Rayleigh waves traveling along a section of the Earth's surface; and (e) Love waves traveling along a section of the earth's surface; (a-e adapted from [Pride, 2005]).

Figure 2 . 4 :

 24 Figure 2.4: Seismic ray path in the case of four layers. The velocity (V) increases with depth in each successive layer (in Z direction). Receivers can be placed at the surface to record seismograms. Picture taken from [Panchuk et al., 2018].

  Figure 2.5: Snell's law: reflection (a), refraction (b) and critical angle (c).

  ; Palermo et al., 2018; Vriend et al., 2015].

Figure 3 .

 3 Figure 3.1: 1D finite volume mesh and computation of unknown wavefields Q and fluxes F through time and space.

  .8) Let W p = α p r p , (3.9)be the jump across the p-th wave component in the solution of the Riemann problem (called the p-wave). Then we have :

Figure 3 . 2 :

 32 Figure 3.2: Finite and spectral elements. Picture taken from Semblat [2008].

Figure 3 . 3 :

 33 Figure 3.3: Passage from a reference element (square) to a quadrangle type element. (η and ξ represent the local coordinates of the reference element. Figure taken from Komatitsch [1997].

Figure 3 . 4 :

 34 Figure 3.4: Grid cell of a 2D staggered spatial finite-difference method of Madariaga [1976] and Virieux [1986] used classically for the discretization of the elastodynamic equations.

Figure 3 . 5 :

 35 Figure 3.5: In model A and B a checkerboard medium composed of two materials is considered. (a) Model A: the source and receivers are placed in the middle of the medium and CPML-type absorbing conditions are implemented on all boundaries. (b) Model B: the source and receivers are placed at the surface and a free surface condition is implemented on the top of the medium and CPML-type absorbing conditions on the other boundaries. (c) Model C: a medium composed of water on top and two other materials beneath, with a free surface condition implemented on the top of the medium and CPML-type absorbing conditions on the other boundaries. The source is placed in the center of the medium and the receivers are placed close to the free surface.

Figure 3 . 6 :

 36 Figure 3.6: For the Model A, snapshots of the vertical component of velocity (v z ) at (a) 9 ms and (b) 27 ms computed with the finite volume code.

Figure 3

 3 Figure 3.7: (a) Horizontal component of velocity (v x ) and (b) vertical component of velocity (v z ) at the different receivers with a point source force excited in the z-direction for the Model A obtained with Finite Volumes (FV_KOREN), Finite Differences (FD) and Spectral Elements (SE).

Figure 3

 3 Figure 3.8: (a) Relative errors of the horizontal component of velocity (v x ) and (b) relative errors of the vertical component of velocity (v z ) for the Finite Volume (FV_KOREN) and Finite Difference (FD) methods according to Spectral Elements (SE) method with a point source force excited in z-direction for the Model A.

Figure 3 .Figure 3 .

 33 Figure 3.9: For the Model A, energy of the system with a point source force excited in the z-direction obtained with Finite Volumes (FV_KOREN) and Finite Differences (FD) methods.

Figure 3 .

 3 Figure 3.11: (a) Horizontal component of velocity (v x ) and (b) vertical component of velocity (v z ) at the different receivers with a point source force excited in the z-direction for the Model B obtained with Finite Volumes (FV_KOREN), Finite Differences (FD) and Spectral Elements (SE).

Figure 3 .

 3 Figure 3.12: (a) Relative errors of the horizontal component of velocity (v x ) and (b) relative errors of the vertical component of velocity (v z ) for the Finite Volume (FV_KOREN) and Finite Difference (FD) methods according to Spectral Elements (SE) method with a point source force excited in zdirection for the Model B.

Figure 3 . 13 :

 313 Figure 3.13: For the Model B, Energy of the system with a point source force excited in the z-direction obtained with Finite Volumes (FV_KOREN) and Finite Differences (FD) methods.

Figure 3 .Figure 3 .

 33 Figure 3.14: For the Model C, snapshots of the horizontal component of velocity (v x ) at 15 ms, 25 ms and 45 ms computed with the finite volume code (a,c,e) and seismic_CPML FD code (b,d,f) respectively.

Figure 3 .

 3 Figure 3.16: (a) Relative errors of the horizontal component of velocity (v x ) and (b) relative errors of the vertical component of velocity (v z ) for the Finite Volume (FV_KOREN) and Finite Difference (FD) methods according to Spectral Elements (SE) method with a point source force excited in zdirection for the Model C.

Figure 3 .

 3 Figure 3.17: Power law model for the dry unconsolidated granular medium case (zero over-pressure) using parameters estimated experimentally

Figure 3 . 19 :

 319 Figure 3.19: For the 10-layer unconsolidated granular medium model, snapshots of the vertical component of velocity (v z ) at 3 ms, 4 ms and 6 ms computed with the finite volume code (a,b,c).

Figure 3 .

 3 Figure 3.20: (a) Horizontal component of velocity (v x ) and (b) vertical component of velocity (v z ) on 12 receivers from the source and spaced 10 mm apart for the 10-layer Model with a point source force excited in the z-direction and receivers placed at the free surface obtained with Finite Volumes (FV_KOREN), Finite Differences (FD) and Spectral Elements (SE).

Figure 3 .

 3 Figure 3.21: (a) Relative errors of the horizontal component of velocity (v x ) and (b) relative errors of the vertical component of velocity (v z ) for the Finite Volume (FV_KOREN) and Finite Difference (FD) methods according to Spectral Elements (SE) method for the 10-layer Model with a point source force excited in z-direction and receivers placed at the free surface.

Figure 3 .

 3 Figure 3.22: Energy of the system for the 10-layer Model with a point source force excited in the z-direction and receivers placed at the free surface obtained with Finite Volumes (FV_KOREN) and Finite Differences (FD) methods.

Figure 3 . 23 :

 323 Figure 3.23: For the unconsolidated granular medium model, snapshots of the vertical component of velocity (v z ) at 2 ms, 4 ms, 6 ms and 13 ms computed with the finite volume code (a,b,c,d), respectively.

Figure 4

 4 Figure 4.1: The physical model (PM) was prepared with 180-300 10 -6 m diameter GBs (GB1) sieved into a (1000 × 800 × 220) . 10 -3 m box. The bottom of the box consisted of a metallic sieve (dashed lines) glued on a perforated plate. The PM HOM22 was prepared by sieving GB1 directly onto the metallic sieve. The bulk density ρ GB1 = 1610 Kg.m -3 was estimated during the deposition process and on samples. Its value led to a porosity of Φ GB1 = 0.356. The x s and y s give the source location (marked by a red star), and the blue lines show the record lines and acquisition parameters. Power law model for the dry unconsolidated granular medium case (zero over-pressure) using parameters estimated experimentally (right).

Figure 4

 4 Figure 4.2: The force source signal (green line on the left inset) is sent from a waveform generator to a low frequency (LF) shaker exciting a metal stick buried in the granular material. The laser beam is set at the zero-offset position (0.250 m in the length direction) to record the vertical component of velocity V z of the stick (the red dashed line in the right inset).

Figure 4 . 3 :

 43 Figure 4.3: Experimental data not filtered (left) and filtered (right).

Figure 4

 4 Figure 4.4: Experimental trace not filtered and filtered at receivers 10 (a), 15 (b) and 20 (c).

Figure 4

 4 Figure 4.5: The 3D numerical model is discretised by elements of (5 × 5 × 5) 10 -4 m for a mesh of 70 × 2000 × 430 points (0.035 × 1 × 0.215 m 3 ). Absorbing PML layers composed of 15 points each one or Dirichlet conditions are applied on the outer boundaries except the free surface. The free surface is defined at the top of the computational domain, and Dirichlet conditions are defined to simulate the edges and the bottom of the box. The force source (red arrow) is implemented with an angle of 20 o from the normal to the free surface, to meet the experimental configuration (see Figures 4.2 and 4.6). The receivers (dashed line) are spaced each 5 × 10 -3 m (10 grid points) at the free surface.

Figure 4

 4 Figure4.6: The 2D numerical model is discretised by elements of (5 × 5) 10 -4 m for a mesh of 2000 × 430 points (1 × 0.215 m 2 ). Absorbing PML layers composed of 15 points each one or Dirichlet conditions are applied on the outer boundaries except the free surface which is defined at the top of the computational domain. The force source (in red) is implemented with an angle of 20 o from the normal to the free surface and injected in every source point, to meet the experimental configuration and the stick form. The receivers (dashed line in green) are spaced each 5 × 10 -3 m (10 grid points) at the free surface.

Figure 4 . 7 :

 47 Figure 4.7: Strong scaling of the UNISOLVER code over 200 up to 400 processors. Ideal and numerical tests scaling curves are shown.

Figure 4 . 8 :

 48 Figure 4.8: The different velocities computed for the purely and poroelastic models.

Figure 4 . 9 :Figure 4 Figure 4

 4944 Figure 4.9: Comparison between experimental (a,e) and 3D elastic (b,f), 2D elastic (c,g) and 2D poroelastic (d,h) seismograms (Vertical component of particle velocity) in the Dirichlet case.

Figure 4 . 13 :

 413 Figure 4.13: Comparison between experimental (a) and 3D elastic (b), 2D elastic (c) and 2D poroelastic (d) spectrograms in the PMLs case.

Figure 4

 4 Figure 4.14: Comparison between experimental (a) and 3D elastic (b), 2D elastic (c) and 2D poroelastic (d) spectrograms in the Full PMLs case.

Figure 4

 4 Figure 4.15: Comparison between experimental (a,e) and 3D elastic (b,f), 2D elastic (c,g) and 2D poroelastic (c,h) dispersion images in the Full PMLs case.

Figure 4 .

 4 Figure 4.16: First arrival times of the 3D elastic, 2D elastic and 2D poroelastic models compared to the experimental and analytical first arrival times.

Figure 4

 4 Figure 4.17: Dispersion curves of 3D elastic (a), 2D elastic (b) and 2D poro-elastic (c) compared to experimental and analytical dispersion curves.

CHAPTER 4 .

 4 MODELING OF UNCONSOLIDATED GRANULAR MEDIA

Figure 4 .

 4 Figure 4.18: Numerical setup for the 8 cm deep source case for the 2D elastic and 2D poroelastic case.

Figure 4 . 19 :Figure 4

 4194 Figure 4.19: Experimental trace image (a) compared to the 2D elastic Finite Difference (b), 2D elastic Finite Volume (c) and 2D poro-elastic (d) trace images in the 8 cm deep source case.

Figure 4

 4 Figure 4.22: Experimental dispersion image (a) compared to the 2D elastic Finite Difference (b) 2D elastic Finite Volume (c) and 2D poro-elastic (d) dispersion images in the 8 cm deep source case.

Figure 4 . 23 :

 423 Figure 4.23: Dispersion curves of 2D elastic (a) and 2D poro-elastic (b) in the 8 cm deep source case compared to experimental and analytical dispersion curves.

Conference⋆

  Kassem Asfour, Roland Martin, Ludovic Bodet, and Didier El Baz. Numerical tools to model seismic waves in unconsolidated and partially saturated granular media. Pubished during the general assembly of European Geoscience Union (EGU) 2021.

Figure 1 :

 1 Figure 1: (a) Horizontal component of velocity (v x ) and (b) vertical component of velocity (v z ) at the different receivers with a point source force excited in the x-direction for the Model A obtained with Finite Volumes (FV_KOREN), Finite Differences (FD) and Spectral Elements (SE).

Figure 2 :

 2 Figure 2: (a) Relative errors of the horizontal component of velocity (v x ) and (b) relative errors of the vertical component of velocity (v z ) for the Finite Volume (FV_KOREN) and Finite Difference (FD) methods according to Spectral Elements (SE) method with a point source force excited in x-direction for the Model A.

Figure 3 :

 3 Figure 3: For the Model A, energy of the system with a point source force excited in the z-direction obtained with Finite Volumes (FV_KOREN) and Finite Differences (FD) methods.

Figure 4 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 11 :

 4678911 Figure 4: Comparison between experimental (a) and 3D elastic (b), 2D elastic (c) and 2D poroelastic (d) F-K dispersion images in the Dirichlet case.

Figure 12 :Figure 13 :Figure 15 :

 121315 Figure 12: Experimental f-k dispersion image (a) compared to the 2D elastic (b) and 2D poro-elastic f-k dispersion images in the 8 cm case.

Figure 16 :

 16 Figure 16: Comparison of the numerical seismograms (2D elastic and 2D poro-elastic) with the experimental ones at receivers 10, 15, and 20 in the 4 cm case.

Figure 17 :

 17 Figure 17: Comparison of the numerical seismograms (2D elastic and 2D poro-elastic) with the experimental ones at receivers 10, 15, and 20 in the 5 cm case.

Figure 18 :

 18 Figure 18: Comparison of the numerical seismograms (2D elastic and 2D poro-elastic) with the experimental ones at receivers 10, 15, and 20 in the 6 cm case.

Figure 19 :

 19 Figure 19: Comparison of the numerical seismograms (2D elastic and 2D poro-elastic) with the experimental ones at receivers 10, 15, and 20 in the 7 cm case.

EXPERIMENTFigure 21 :

 21 Figure 21: Experimental f-k dispersion images (a,d,g,j) compared to the 2D elastic (b,e,h,k) and 2D poro-elastic f-k dispersion images (c,f,i,l) in the 4 cm (b,c), 5 cm (e,f), 6 cm (h,i), 7 cm (k,l) cases.

EXPERIMENTFigure 23 :

 23 Figure 23: Experimental dispersion images (a,d,g,j) compared to the 2D elastic (b,e,h,k) and 2D poroelastic dispersion images (c,f,i,l) in the 4 cm (b,c), 5 cm (e,f), 6 cm (h,i), 7 cm (k,l) cases.

  champ de contrainte et à calculer les autres composantes au niveau de la frontière ainsi que les différentes composantes du champ de vitesse[Levander, 1988a;[START_REF] Robertsson | A numerical free-surface condition for elastic/viscoelastic finitedifference modeling in the presence of topography[END_REF]. Nous obtenons plus de précision en implémentant un schéma spatial d'ordre 4 non centré près de la surface libre proposé par[START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF]. Par ailleurs, les conditions de Dirichlet (bord rigide) consistent à annuler toutes les composantes du champ de vitesse à la frontière.• La mise en oeuvre d'une méthode de volumes finis et sa comparaison avec les différences finies classiques avec grilles en quinconce ou des éléments finis comme la méthode des éléments spectraux. En effet, nous proposons ici une méthode de volumes finis basée sur un solveur de Riemann (méthode RFV-FSP/Riemann Finite Volume-Fluxes frequency Shift PML) pour calculer les champs d'ondes sismiques sur des grilles colocalisées ainsi qu'une formulation des conditions absorbantes (C-PML) plus spécifiquement conçues pour la

	les
	propriétés physiques du milieu (densité et vitesses sismiques) à un même point de grille,
	ce qui permet plus de flexibilité et de facilités pour traiter numériquement la surface libre
	et les interfaces solide-solide ou fluide-solide. De plus, cette version du FVM est capa-
	ble d'introduire plus facilement le dérafinement des grilles lorsque les vitesses sismiques
	augmentent avec la profondeur (comme c'est généralement le cas dans de nombreuses
	applications géophysiques) : dans ce cas, les cellules de la grille peuvent être agrandies à
	différentes profondeurs par de simples interpolations aux interfaces entre une grille fine
	et une autre plus grossière.

pour les stencils de type 'standard staggered grid' (SSG), dont le principe consiste à annuler les composantes verticales du méthode des volumes finis. Les conditions aux limites de type PML sont optimisées à incidence rasante en utilisant les 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 est une formulation différente des PMLs classiques qui sont généralement appliquées aux dérivées spatiales des variables primitives (vitesses et contraintes des particules). Elles peuvent même être améliorées en introduisant des conditions absorbantes de type éponge ou paraxiales. Les volumes finis sont capables de calculer les variables physiques (vitesses et contraintes des particules) et de prendre en compte

Plan de la thèse

	l'inventaire des différentes méthodes numériques utilisées ces dernières années pour modéliser la
	propagation des ondes. Nous développons également les aspects théoriques de la méthode des
	volumes finis utilisée pour le développement du code élastique. Nous décrivons une méthode de
	différences finies (pour comparer avec la méthode des volumes finis), en particulier les schémas
	de discrétisation Standard Staggered Grid (SSG-FD), ainsi que la méthode des éléments finis
	spectraux utilisée dans SPECFEM (prise comme solution de référence) : [Casarotti et al.,
	2008; Komatitsch, 1997, 2003, 2008; Rietmann et al., 2012; Wang and Cai, 2017]. Après avoir
	développé différents points relatifs aux conditions absorbantes et à la mise en oeuvre de la
	surface libre, nous concluons ce chapitre par une étude dans laquelle nous considérons trois
	modèles principaux : un milieu hétérogène en damier (pour valider les conditions d'absorption
	et de surface libre), un modèle fluide-solide avec de l'eau au sommet du milieu et deux couches
	solides en dessous (pour tester la stabilité des conditions PMLs pour les simulations de couplage
	acoustique-élastique) et un modèle réaliste (modèle en loi de puissance) déterminé dans le
	contexte des milieux granulaires non consolidés.
	Dans le chapitre 3, nous revisitons une étude réalisée dans
	. Les données expérimentales enregistrées sur
	des milieux granulaires non consolidés sont comparées qualitativement et quantitativement à
	des sismogrammes synthétiques. Ce travail a été réalisé en collaboration avec des chercheurs
	du laboratoire METIS/Paris Sorbonne (CNRS, Université Pierre et Marie Curie) (Laboratoire
	des Milieux environnementaux, transferts et interactions dans les hydrosystèmes et les sols) qui
	nous a fourni les données expérimentales issues d'une expérience de laboratoire réalisée sur des
	milieux granulaires non consolidés. Cette expérience avec des modèles plus réalistes ouvre la
	voie à d'autres expériences, à savoir l'injection d'air et/ou de fluides et l'introduction d'une
	topographie non-horizontale.
	Le chapitre 1 est consacré aux rappels théoriques des équations d'ondes dans un milieu
	élastique (3D et 2D) et poro-élastique (2D). Les différents types d'ondes, l'implémentation
	des sources, les différentes conditions aux limites et les outils d'analyse temps/dispersion sont
	également développés.
	Dans le chapitre 2, nous présentons un état de l'art de la modélisation sismique en faisant

  In these equations, (v x , v y , v z ) are the velocity components; (σ xx , σ yy , σ zz , σ xy , σ xz , σ yz ) are the stress components; (s x , s y , s z ) are the body-force components; ρ is the density; λ and µ are Lamé parameters.

	ρ	∂v x ∂t	=	∂σ xx ∂x	+	∂σ xy ∂y	+	∂σ xz ∂z	+ s x ,
	ρ	∂v y ∂t	=	∂σ xy ∂x	+	∂σ yy ∂y	+	∂σ yz ∂z	+ s y ,
	ρ	∂v z ∂t	=	∂σ xz ∂x	+	∂σ yz ∂y	+	∂σ zz ∂z	+ s z ,
	∂σ xx ∂t	= (λ + 2µ)	∂v x ∂x	+ λ(	∂v y ∂y	+	∂v z ∂z	),
	∂σ yy ∂t	= (λ + 2µ)	∂v y ∂y	+ λ(	∂v x ∂x	+	∂v z ∂z	),
	∂σ zz ∂t	= (λ + 2µ)	∂v z ∂z	+ λ(	∂v x ∂x	+	∂v z ∂z	),
	∂σ xy ∂t	= µ	∂v y ∂x	+	∂v x ∂y	,
	∂σ xz ∂t	= µ	∂v z ∂x	+	∂v x ∂z	,
	∂σ yz ∂t	= µ	∂v z ∂y	+	∂v y ∂z	,	(2.15)
										):

Table 2 . 2 :

 22 VsVp and V R Vs speed ratios depending on the Poisson's ratio.

  al., 2014; Kreiss et al., 1967; Lombard et al., 2008; Minkoff, 2002; Sun et al., 2018; Zhang, 1997, 2012; Zhang et al., 2011].

  ]. It is based either on the fast Fourier transform/FFT technique in the frequency domain[START_REF] Kosloff | Elastic wave calculations by the Fourier method[END_REF][START_REF] Kosloff | Forward modeling by a Fourier method[END_REF][START_REF] Zhang | A new high accuracy locally one-dimensional scheme for the wave equation[END_REF] or on the Chebyshev transform in the time domain[START_REF] Carcione | The wave equation in generalized coordinates[END_REF][START_REF] Carcione | Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media[END_REF]; Carcione et al., 2005; Danecek and Seriani, 2008; Klin et al., 2010; Komatitsch et al., 1996; Sidler et al., 2014;[START_REF] Wang | Modelling seismic wave propagation in a two dimensional cylindrical whole earth model using the pseudospectral method[END_REF]. It has been one of the most used numerical methods because of its accuracy, and a minimum number of grid points is needed to reach the Nyquist frequency and to model wave propagation properly with very low numerical dispersion. However, the PSM is a global method and in general it is difficult to apply parallelization and local adaptivity. More particularly, a large amount of forward and inverse FFT reduces its computational efficiency.

• The Finite Element Method/FEM [Bécache et al., 2002; Cohen et al., 2001; Givoli et al., 2006; Givoli and Neta, 2003; Priolo et al., 1994;

  FEM that commonly involves non-diagonal mass matrices. The SEM is now a popular class of continuous Galerkin techniques as well as a class of high-order FEM. It combines the benefits of both the pseudospectral and the spectral methods.• The Discontinuous Galerkin Method/DG originally introduced by Reed and Hill[1973] for solving linear neutron transport equations is widely used in seismology to solve the elastodynamic equations in the time and frequency domains [Bonnasse-Gahot et al., 2018; Chung and Engquist, 2009; de la Puente et al., 2007; Dumbser et al., 2007; He et al., 2015; Käser and Dumbser, 2006; Käser et al., 2007; Terrana et al., 2018]. It is a high order conservative method. Different time-integration methods can be used but one of the most popular is the Arbitrary high-order DERivatives (ADER) approach which is based on upwind fluxes and high order interpolation in space and time [Dumbser and Käser, 2006].High order schemes in space can be constructed easily. The DG method is very flexible, easy to parallelize (the mass matrix is local rather than global), and contrarily to the FEM, the solutions can be discontinuous at the interfaces between the elements. However, the DG method is computationally (CPU time and memory occupancy) more expensive than SEM or more classical FEM due to the computation of high-order surface and volume integrals. A convergence analysis of DG methods show that this technique proved to be accurate and efficient[START_REF] Engsig-Karup | Dg-fem solution for nonlinear wave-structure interaction using boussinesq-type equations[END_REF][START_REF] Lanteri | Convergence of a discontinuous galerkin scheme for the mixed time-domain maxwell's equations in dispersive media[END_REF]. To overcome the computational cost of the DG method, and without inducing any additional difficulty for the parallelization, an Hybridizable Discontinuous Galerkin/HDG method is proposed[START_REF] Terrana | A spectral hybridizable discontinuous Galerkin method for elastic-acoustic wave propagation[END_REF]].• The Finite Volume Method/FVM is an integral method very popular for the resolution of linear and non-linear hyperbolic equations for fluid mechanics (Navier Stokes equations) and seismic wave propagation. Contrary to the FEM methods based on a variational formulation of the problem, the FVM is based on the "strong" formulation. The solutions are calculated using a piecewise constant approximation and by evaluating approximated fluxes at the interfaces between the finite volumes (intervals in 1D, surfaces in 2D and volumes in 3D)[START_REF] Chung | Convergence analysis of a finite volume method for Maxwell's equations in nonhomogeneous media[END_REF]; Chung and Engquist, 2005; Dumbser et al., 2007; Jianfeng and Tielin, 1999, 2002; LeVeque, 2002a].

;

Komatitsch et al., 2001;[START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF], 2001, 2002, 2003;[START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF][START_REF] Seriani | Elastic wave propagation in complex heterogeneous earth structures: numerical modelling by using a poly-grid spectral element method[END_REF]

. In this method, a diagonal mass matrix

[START_REF] Bernardi | Mixed spectral element approximation of the Navier-Stokes equations in the stream-function and vorticity formulation[END_REF][START_REF] Cohen | Mixed spectral finite elements for the linear elasticity system in unbounded domains[END_REF][START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF]

] is introduced by using Gauss-Lobatto-Legendre quadrature rules which makes the SEM easier to parallelize than the

Remerciements

J'ai effectué mes travaux de thèse à Toulouse entre le GET (Géosciences Environnement Toulouse) de l'OMP (Observatoire Midi-Pyrénées) et le LAAS-CNRS (Laboratoire d'Analyse et d'Architecture des Systèmes). Je tiens à remercier Monsieur Étienne Ruellan et Monsieur Sylvain Bonvalot, directeurs successifs du GET, ainsi que Monsieur Liviu Nicu, directeur du LAAS-CNRS, pour m'avoir accueilli au sein de leurs laboratoires et d'avoir mis à ma disposition toutes les ressources nécéssaires au bon déroulement de ma thèse.

Impact of the numerically modelled source

Source configuration

Instead of considering a source defined by a stick with vibratory displacement, we consider a point source at different depths with the same inclination as the stick and the same time signal source injected during the experiment. The different depths considered are 3.76, 4.7, 5.64, 6.58, 7.52 cm. Let us denote these different cases respectively by 4, 5, 6, 7, 8 cm which correspond to the length of the stick as in Figure 4.18. The best results are obtained in the case of a 8 cm stick length in which the source is modelled as a point source. We choose this case because better waveforms and better modes description have been obtained as discussed hereafter. For the other cases, the reader is referred to appendix A.6.

In addition to the 2D finite difference codes (elastic and poroelastic), we use our RFV (finite volume) code with the same setup configuration as the 2D elastic Finite Difference code. The trace and the dispersion image of the Finite Volume method compared to the Finite difference one are very similar as you can see in Figures 4. 19, 4.21 and 4.22. RFV can be considered as another good candidate to solve this kind of near surface problem.

Waveforms comparison

Compared to the seismograms of the experiment, the numerical seismograms (2D elastic and 2D poro-elastic models) in the 8 cm deep source case are more similar (Figure 4. [START_REF] Zhang | Quadrangle-grid velocity-stress finite-difference method for elastic wavepropagation simulation[END_REF]) than with the pure stick case (Figure 4.11). Quantitatively, we can see in Figure 4.20 that the first numerical waveforms of the 2D elastic and 2D poro-elastic models are very similar to the experimental ones. We also noticed that the computed Rayleigh waveforms have overestimated amplitudes when compared to experimental surface waveforms after ≈ 0.012 s. On the contrary, in the 8 cm deep source case, the amplitudes of surface waves are similar to the amplitudes of the experimental data, which is very encouraging.

Dispersion analysis

In Figure 4.22, we can see the very good similarity of the numerical dispersion images with the experimental one. In the previous section, the P-modes (mode 2) were not clearly reproduced in terms of amplitude. On the contrary, here we obtain the same patterns with very pronounced P-modes. These modes are very well represented in both elastic and poro-elastic cases. In the poro-elastic case, results are a little bit more improved. This can be also clearly observed on the dispersion curves of Figure 4.23. All these results show the crucial importance of the data sensitivity to the source modelling and particularly to the source location at depth. The deeper is the source the more the P -modes can be recovered and moved towards lower frequencies. Since these P -modes are more pronounced, the P -SV wave modes are less and less visible even if they are still existing but with less intensity. Of course, if the source location

Appendix

A.1 Compact form of the 2D elastic wave equations using

We rewrite equations (3.1) in the following compact form:

where 2D POROELASTIC