Impact of Source Modelling and Poroelastic Models on Numerical Modelling of Unconsolidated Granular Media: Application at the Laboratory Scale

Abstract

The near surface is characterized by using different numerical techniques, among them seismic techniques that are non-destructive. More particularly, for a better understanding of acoustic and seismic measurements in unconsolidated granular media that can constitute the near surface, many studies have been conducted in situ and also at the laboratory scale where theoretical models have been developed. In this article, we want to model such granular media that are difficult to characterize. At the laboratory scale, dry granular media can be modelled with a homogenized power-law elastic model that depends on depth. In this context, we validate numerically a similar power-law elastic model for such media by applying it to a homogenized elastic medium or to the solid frame of a poroelastic medium that consists of solid and air components. By comparing the response of both rheologies, we want to highlight what poroelastic media can bring to better reproduce the experimental data in the time and frequency domains. To achieve this objective, we revisit studies carried out on unconsolidated granular media at the laboratory scale and we compare different models with different rheologies (elastic or poroelastic), dimensions (2D or 3D), boundary conditions (perfectly matched layer/PML, or Dirichlet) and locations of the source (modelled as a vibratory stick or a point force) in order to reproduce the experimental data. We show here that a poroelastic model describes better the amplitudes of the seismograms. Furthermore, we study the sensitivity of the seismic data to the source location, which is crucial to improve the amplitude of the signals and the detection of the different seismic modes.

Appendices

Appendix A. 3D Elastic Wave Equations

The elastodynamics equations can be formulated at the second order in displacement as:

$$\begin{aligned}&\rho \frac{\partial ^2 u_i}{\partial t^2} = \partial _j \sigma _{ij} + s_i ,\\&\epsilon _{ij} = \frac{1}{2} \left( u_{i,j} + u_{j,i} \right) , \\&\sigma _{ij} = \lambda \delta _{ij} \epsilon _{kk} + 2 \mu \epsilon _{ij}. \end{aligned}$$

However, the first-order formulation (velocity-stress) of the 3D elastic wave equations for a linear, isotropic medium submitted to external forces is given by Graves (1996):

$$\begin{aligned}&\rho \frac{\partial v_x}{\partial t} = \frac{\partial \sigma _{xx}}{\partial x} + \frac{\partial \sigma _{xy}}{\partial y} + \frac{\partial \sigma _{xz}}{\partial z} + s_x \nonumber ,\\&\rho \frac{\partial v_y}{\partial t} = \frac{\partial \sigma _{xy}}{\partial x} + \frac{\partial \sigma _{yy}}{\partial y} + \frac{\partial \sigma _{yz}}{\partial z} + s_y \nonumber ,\\&\rho \frac{\partial v_z}{\partial t} = \frac{\partial \sigma _{xz}}{\partial x} + \frac{\partial \sigma _{yz}}{\partial y} + \frac{\partial \sigma _{zz}}{\partial z} + s_z \nonumber ,\\&\frac{\partial \sigma _{xx} }{\partial t} = \left( \lambda + 2 \mu \right) \frac{\partial v_x}{\partial x} + \lambda (\frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z}) \nonumber ,\\&\frac{\partial \sigma _{yy} }{\partial t} = \left( \lambda + 2 \mu \right) \frac{\partial v_y}{\partial y} + \lambda (\frac{\partial v_x}{\partial x} + \frac{\partial v_z}{\partial z}) \nonumber ,\\&\frac{\partial \sigma _{zz} }{\partial t} = \left( \lambda + 2 \mu \right) \frac{\partial v_z}{\partial z} + \lambda ( \frac{\partial v_x}{\partial x} + \frac{\partial v_z}{\partial z} ) \nonumber ,\\&\frac{\partial \sigma _{xy}}{\partial t} = \mu \left( \frac{\partial v_y}{\partial x} + \frac{\partial v_x}{\partial y}\right) \nonumber ,\\&\frac{\partial \sigma _{xz}}{\partial t} = \mu \left( \frac{\partial v_z}{\partial x} + \frac{\partial v_x}{\partial z}\right) \nonumber ,\\&\frac{\partial \sigma _{yz}}{\partial t} = \mu \left( \frac{\partial v_z}{\partial y} + \frac{\partial v_y}{\partial z}\right) . \end{aligned}$$

(2)

In these equations, $v_x,$ $v_y,$ $v_z$ are the velocity components; $\sigma _{xx},$ $\sigma _{yy},$ $\sigma _{zz},$ $\sigma _{xy},$ $\sigma _{xz},$ $\sigma _{yz}$ are the stress components; $s_x,$ $s_y,$ $s_z$ are the body-force components; $\rho$ is the density; $\lambda$ and $\mu$ are Lamé parameters.

Appendix B. 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 (Dumbser and Käser 2006):

$$\begin{aligned}&\rho \frac{\partial v_x}{\partial t} = \frac{\partial \sigma _{xx}}{\partial x} + \frac{\partial \sigma _{xz}}{\partial z} + s_x \nonumber ,\\&\rho \frac{\partial v_z}{\partial t} = \frac{\partial \sigma _{xz}}{\partial x} + \frac{\partial \sigma _{zz}}{\partial z} + s_z \nonumber ,\\&\frac{\partial \sigma _{xx} }{\partial t} = \left( \lambda + 2 \mu \right) \frac{\partial v_x}{\partial x} + \lambda \frac{\partial v_z}{\partial z} \nonumber ,\\&\frac{\partial \sigma _{zz} }{\partial t} = \left( \lambda + 2 \mu \right) \frac{\partial v_z}{\partial z} + \lambda \frac{\partial v_x}{\partial x} \nonumber ,\\&\frac{\partial \sigma _{xz}}{\partial t} = \mu \left( \frac{\partial v_z}{\partial x} + \frac{\partial v_x}{\partial z}\right) , \end{aligned}$$

(3)

where ${\lambda}$ and $\mu$ are Lamé parameters, $\rho$ 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 $\sigma _{xx}$ and $\sigma _{zz}$, and the shear stress is $\sigma _{xz}$. The components of particle velocities in direction x and z are denoted by $v_x$ and $v_z$, respectively.

Appendix C. 2D Poroelastic 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 (Pride et al. 2004). 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) and Biot (1956b)). The differential or "strong" formulation of the poroelastic wave equations can be written as (Carcione 2007), (Carcione 2014):

$$\begin{aligned} \rho \partial _t^2 u^s + \rho _f \partial _t^2 w&= \nabla . \left( C : \nabla u^s - \alpha P^f I\right) , \end{aligned}$$

(4)

$$\begin{aligned} \rho _f \partial _t^2 u^s + \rho _w \partial _t^2 w&= - \nabla P^f - K \partial _t w, \end{aligned}$$

(5)

$$\begin{aligned} P^f&= - \alpha M \nabla . u^s - M \nabla . w, \end{aligned}$$

(6)

where $u^s = (u^s_i)_{i=1,2}$, $w = \Phi \left( u^f - u^s\right)$ and $u^f = (u^f_i)_{i=1,2}$ are, respectively, the solid, relative, and fluid displacement vectors; $\Phi$ is the porosity; and C is the stiffness tensor of the isotropic elastic solid matrix, defined as:

$$\begin{aligned} \sigma _{ij}^s&= \left( C :\epsilon \right) _{ij} = \lambda _s \delta _{ij}\epsilon _{kk} + 2 \mu \epsilon _{ij}, \end{aligned}$$

(7)

$$\begin{aligned} \epsilon _{ij}&= \frac{1}{2} \left( \frac{\partial u_i^s}{\partial x_j} + \frac{\partial u_j^s}{\partial x_i}\right) , \end{aligned}$$

(8)

$$\begin{aligned} P^f&= - \alpha M \nabla . u^s - M \nabla . w, \end{aligned}$$

(9)

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. $\sigma ^s$ and $\epsilon$ are, respectively, the stress and strain tensors of the isotropic elastic solid frame. The stress tensor of the fluid-filled solid matrix is $\sigma = \sigma ^s - \alpha P^f I$, and $\rho = \Phi \rho _f + \left( 1- \Phi \right) \rho _s$ is the density of the saturated medium, where $\rho _s$ and $\rho _f$ are the solid and fluid densities, respectively. The apparent density is $\rho _w = a \frac{\rho _f}{\Phi }$ where a denotes the tortuosity. The shear modulus is $\mu$ and $\lambda _s = \lambda - \alpha ^2 M$ is the Lamé coefficient in the solid matrix, where $\lambda$ is the Lamé coefficient of the saturated matrix. The $\alpha$ 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:

$$\begin{aligned} \alpha&= 1 - \frac{K_{fr}}{K_s} , \end{aligned}$$

(10)

$$\begin{aligned} M&= 1/\left[ \Phi /K_f + (\Phi -\alpha )/K_s\right] , \end{aligned}$$

(11)

where $K_{fr}$ 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:

$$\begin{aligned} K = \kappa / \eta , \end{aligned}$$

(12)

where $\kappa$ is the permeability of the solid matrix and $\eta$ is the fluid viscosity. Equations (4) to (7) can be written using a first-order velocity-stress formulation:

$$\begin{aligned} \left( \rho _w \rho - \rho _f^2 \right) \partial _t v^s&= \rho _w \nabla . \sigma + \rho _f \nabla P^f + \rho K v^f, \end{aligned}$$

(13)

$$\begin{aligned} \left( \rho _w \rho - \rho _f^2 \right) \partial _t v^f&= -\rho _f \nabla . \sigma - \rho \nabla P^f - \rho _f K v^f, \end{aligned}$$

(14)

$$\begin{aligned} \partial _t \sigma&= C: \nabla v^s - \alpha \partial _t P^f I, \end{aligned}$$

(15)

$$\begin{aligned} \partial _t P^f&= - \alpha M \nabla . v^s - M \nabla . v^f, \end{aligned}$$

(16)

where $v^s = (v^s_i)_{i=1,2}$ and $v^f = \partial _t w = (v^f_i)_{i=1,2}$ are the solid and filtration velocity vectors, respectively. $\sigma$ is the effective stress tensor of the porous medium. As in Zeng and Liu (2001), using the trace of the strain tensor $Tr(\epsilon )= \epsilon _{ii}$ and an auxiliary variable $\xi$, we rewrite the system as:

$$\begin{aligned} \left( \rho _w \rho - \rho _f^2 \right) \partial _t v^s_i&= \rho _w \partial _j \sigma _{ij} + \rho _f \partial _i P^f + \rho K v^f_i, \end{aligned}$$

(17)

$$\begin{aligned} \left( \rho _w \rho - \rho _f^2 \right) \partial _t v^f_i&= -\rho _f \partial _j \sigma _{ij} - \rho \partial _i P^f - \rho _f K v^f_i, \end{aligned}$$

(18)

$$\begin{aligned} \epsilon _{ij}&= \frac{1}{2} \left( \partial _j v_i^s + \partial _i v_j^s\right) , \end{aligned}$$

(19)

$$\begin{aligned} \partial _t \xi&= -\partial _i v_i^f, \end{aligned}$$

(20)

$$\begin{aligned} P^f&= - \alpha M Tr(\epsilon ) + M \xi , \end{aligned}$$

(21)

$$\begin{aligned} \sigma _{ij}^s&= \lambda ^s \delta _{ij} Tr(\epsilon ) + 2 \mu \epsilon _{ij}, \end{aligned}$$

(22)

$$\begin{aligned} \sigma _{ij}&= \sigma _{ij}^s -\alpha P^f \delta _{ij}. \end{aligned}$$

(23)

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 $\pm V_{pFAST}$, $\pm V_{pSLOW}$, $\pm V_{s}$ and 0. The shear velocity $V_s$ and the fast and slow P-wave velocities ($V_{pFAST}$ and $V_{pSLOW}$) can be expressed as (Sidler et al. 2014):

$$\begin{aligned} V_s&= \sqrt{ \frac{\mu }{a_1}}, \end{aligned}$$

(24)

$$\begin{aligned} {V_p}_{FAST}&= \sqrt{ \frac{-b_1 + \sqrt{\Delta }}{2 a_1}}, \end{aligned}$$

(25)

$$\begin{aligned} {V_p}_{SLOW}&= \sqrt{ \frac{-b_1 - \sqrt{\Delta }}{2 a_1}}, \end{aligned}$$

(26)

where

$a_1 = \rho _{11} \rho _{22}-\rho _{12}^2,$

$b_1 = -S \rho _{22} - R \rho _{11} + 2 ga \rho _{11},$

$c_1 = S R - ga^2,$

$\rho _{11} = \rho + \rho _f \Phi (a - 2),$

$\rho _{12} = \Phi \rho _f (1-a),$

$\rho _{22} = a \rho _f \Phi ,$

$S = \lambda + 2 \mu ,$

$R = M \Phi ^2,$

$\Delta = b_1^2 - 4 a_1 c_1,$

$ga = M \Phi (\alpha - \Phi ).$

In Table 1, a summary of the different parameters of the poroelastic model is provided.

Table 1 Parameters of the poroelastic model

Biot’s characteristic frequency $f_c$ defines the transition between two poroelastic regimes (with or without attenuation) and is given as follows (see Biot 1956b; Carcione 2007 and Morency and Tromp 2008):

$$\begin{aligned} f_c = \min (\frac{\eta \Phi }{ 2 \pi a \rho _f \kappa }), \end{aligned}$$

(27)

see parameters in Table 1.

In our study, the maximum frequency range $f_{max}$ of the source is such that $f_{max}< f_c$. Therefore, in the experimental and numerical modelling of unconsolidated granular media under study, we choose to stay in the poroelastic regime without attenuation.

Appendix D. Performance of UNISOLVER Parallel Code

The fourth-order 3D UNISOLVER parallel code has been scaled over different numbers of processors (up to 400) of Olympe supercomputer at CALMIP computing centre of Toulouse (France) by using MPI (Message-Passing-Interface) library. The computational domain has been cut along the longitudinal y-axis (from 20 points down to 5 grid points per processor along the y-axis). A buffer overlapping corresponding to 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’ communication operations. As we can see in Fig. 22, the strong scaling obtained by measuring the CPU time versus the number of processors is very satisfactory, despite the use of blocking communication operations. This strategy is classical (see (Komatitsch and Martin 2007)).

Fig. 22

Strong scaling of the UNISOLVER code over 200 up to 400 processors on Olympe supercomputer. Ideal and numerical tests scaling curves are shown

We underline that asynchronous iterative schemes of computation which cover communications by computations in the inner subdomains as depicted in Miellou et al. 1998; El Baz et al. 2001; Chau et al. 2007; El Baz et al. 2005; Martin et al. 2008a are also very efficient.