TY - JOUR
AU - Rejiba,, Fayçal
AB - SUMMARY In previous studies, the auxiliary differential equation (ADE) method has been applied to 2-D seismic-wave propagation modelling in viscoelastic media. This method is based on the separation of the wave propagation equations derived from the constitutive law defining the stress–strain relation. We make here a 3-D extension of a finite-difference (FD) scheme to solve a system of separated equations consisting in the stress–strain rheological relation, the strain–velocity and the velocity–stress equations. The current 3-D FD scheme consists in the discretization of the second order formulation of a non-linear viscoelastic wave equation with a time actualization of the velocity and displacement fields. Compared to the usual memory variable formalism, the ADE method allows flexible implementation of complex expressions of the desired rheological model such as attenuation/viscoelastic models or even non-linear behaviours, with physical parameters that can be provided from dispersion analysis. The method can also be associated with optimized perfectly matched layers-based boundary conditions that can be seen as additional attenuation (viscoelastic) terms. We present the results obtained for a non-linear viscoelastic model made of a Zener viscoelastic body associated with a non-linear quadratic strain term. Such non-linearity is relevant to define unconsolidated granular model behaviour. Thanks to a simple model, but without loss of generality, we demonstrate the accuracy of the proposed numerical approach. Elasticity and anelasticity, Nonlinear differential equations, Numerical modelling, Computational seismology, Seismic attenuation, Wave propagation 1 INTRODUCTION In seismology or near surface geophysics applications, imaging and characterizing the Earth materials at different scales using seismic wave propagation can face attenuation and non-linear behaviour because of the physical nature of the media under study and the frequency content of the seismic sources. It can involve a wide variety of physical, geometrical, kinematical and structural types of non-linearities as well as combination of these non-linearities types (Ostrovsky & Johnson 2001; Delsanto 2006). Attenuation and non-linearity parametrizations in seismic modelling are not an easy task and are subject to many studies in the last two decades. Non-linear dynamics applied to geomaterials can be essentially due to the presence of soft features where damage is primarily observed as in granular media where the physical source of the non-linearities associated to grain-to-grain interactions has been well identified (Tournat & Gusev 2010; Legland et al.2012). At a larger scale, those non-linearities are arising during the wave propagation and are commonly studied in the context of site effect assessment, and related resonance phenomena or strong ground motion (de)amplification responses. In presence of non-linear site effects, viscous damping or non-linear ground behaviours have also been studied (Bonilla et al.2005, 2011; Régnier et al.2013, 2014, 2016b; Sandikkaya et al.2013; Akkar et al.2014). Actually, attenuation and non-linear characteristics of a medium are frequency dependent and are thus not easy to model. They are a consequence of complex mixture of superimposed distinct linear behaviour of each phase composing the media. In practice, viscoelastic parameters are commonly estimated from seismic data and dispersion analysis, using empirical measurement of the quality factor Q, whether considered as quasi constant in the simpliest and common case over a specific frequency range or as frequency-dependent in case of dispersive and non-linear behaviour. But in many geophysical applications, this approximation is too restrictive because the medium under study can have strong non-linear and attenuation behaviour depending on the seismic/acoustic signal amplitude and the frequency range. In the case of unconsolidated media (granular media, sedimentary basins, soften/weak fault systems), the parametrization of the attenuation and non-linearities must be done for different frequency ranges. Several experiments done in this context have shown resonance effects and slow dynamics phenomena, with the appearance of harmonics in the frequency domain, that are evidences of attenuation and non-linearities dependent on frequency and strain amplitude. As an example of application to earthquake dynamics, soften and weakened materials located at the seismic source, more particularly in the fault gouge area, can be at the origin of strong non-linearities and exhibit non-equilibrium and non-linear dynamics. Earthquake triggering can then be explained by fault failure caused by already weakened/jammed fault systems that are at critical and weak state and excited by seismic-waves impinging the fault system with a sufficiently high strain amplitude (Johnson & Jia 2005). The seismic waves produced by a remote earthquake are exciting the fault and make the fault core modulus to decrease strongly, and the fault core weakens further. However, the strains should reach a sufficiently high threshold value to cause triggering phenomena. In such situation, the stress–strain rheological relations characterizing these phenomena are generally non-linear and can be expressed as Taylor expansions of complex rheology expressions or high polynomial (at least cubic) functions of the strain close to the critical state. Such non-linear rheologies can explain how some earthquakes can be triggered remotely by other huge earthquakes (Gomberg et al.2003; Gomberg & Johnson 2005) like the 7.3 Mw 1992 Landers earthquake (Hill et al.1993; Gomberg et al.2001) or 7.1 Mw Hector Mine (Gomberg et al.2001) and 7.9 Mw Denali earthquakes (Gomberg et al.2004). In the context of the 2011 Tohoku earthquake, non-linear viscoelastic wave propagation modelling has been able to reproduce displacement amplifications in the near surface. This has been enabled by adding non-linearities in both the elastic and attenuation parts of the stress tensor attenuation as well as hysteretic effects (d’Avila et al.2013) in the ground response due to strong seismic motion. As another example, several harmonics appearing in the observed accelerations spectra related to the Northridge earthquake (Delépine et al.2009) are evidences of non-linearity effects generated by frequency modulation phenomena. At intermediate scales as for geotechnical applications, handling the interaction between the structure and the ground at the same time is crucial, particularly because of the potentially significant degradation of the shear modulus (Kramer 1996; Delépine et al.2009). It is another reason why using non-linear stress–strain relations is interesting in seismic wave propagation modelling. At local scales, non-linear extensions of standard elastodynamic equations have been suggested to better describe wave propagation in complex porous media. For instance, the poroelastic wave equations given in Biot (1956a,b) have been extended to non-linear poroelastic wave equations (Biot 1973) in the particular context of Hertz–Mindlin rheology. Such approach has been applied to wave processes with strong acoustic non-linearity in porous rubber- or sandstone-like media mainly in one dimension (Ostrovsky 1991) or to porous media containing spherical or cylindrical pores by Donskoy et al. (1997) and to unconsolidated granular media by Dazel & Tournat (2010). This theoretical framework has been developed in the acoustic community in an effort to understand non-linear measured responses when the source amplitude is increased for instance in compact granular materials under gravity. In 1-D, the interest of the geophysical community to non-linearities is obvious as 1-D non-linear response validations and verifications using different numerical solutions have recently focused the attention on canonical/benchmark cases (Régnier et al.2016a, 2018). Non-linear models in pure solid or poroelastic rheologies have been studied mainly by developing 1-D analytical solutions in the frequency domain in the context of dispersive granular media (Hokstad 2004; Tournat et al.2003, 2004; Tournat & Gusev 2010; Legland et al.2012) or in the time domain for non-linear constitutive laws (Berjamin et al.2017) as well as 1-D finite differences (FD) in the time domain (Favrie et al.2014) or recent 2-D lattice approaches (Wallen & Boechler 2017). In the studies of Tournat’s group, non-linear effects have been exhibited via second harmonics generated by, for instance, non-linear Hertzian contact rheologies or shear to longitudinal mode conversions in granular packed media. 1-D analytical models are essentially formulated using power law-based non-linear elastic and acoustic rheologies for granular and unconsolidated dispersive media. Power law degrees equal or lower than 2 are generally used and quadratic Hertz–Mindlin stress–strain relations are introduced with one non-linear parameter homogeneous to the shear or bulk modulus of granular media at the laboratory scale. In Favrie et al. (2014) a 1-D non-linear version of the Zener model is proposed. The generalized Zener viscoelastic model degenerates correctly towards a pure non-linear elastic model when attenuation effects vanish. Besides, the non-linear elastic stress–strain relations under study are Lennard–Jones potential based functions or cubic functions like the Landau’s model which are widely used in non-destructive testing and slow dynamics modelling. High order polynomial stress–strain functions are commonly introduced in the frequency domain to homogenize solid–solid or solid–fluid interactions from pore to microscale. Ostrovsky & Johnson (2001) made a review of different non-linear elastic behaviours of earth materials and their impact in earth and material science (strong ground motion, non-destructive testing). They describe different sources of non-linearities in waves travelling through rocks due for instance to geometry of cracks and solid matrices, grain shapes, fluid content/saturation effects, the nature of soften/weak materials (presence of cracks, grain–grain or grain–fluid contacts, etc.), as well as hysteresis, relaxation or attenuation effects in these materials (slow dynamic effects). In many experiments, second and third harmonic spectrum amplitudes are evidence of non-linear behaviour and are varying according to the strain amplitude of the fundamental mode and with quadratic and cubic frequencies. In this context the authors give a theoretical high order polynomial expression of the stress–strain relation in which non-linear coefficients are appearing and can be modelled as combinations of the elastic moduli. The non-linear rheologies basically are derivations of the second potential invariant and can be expressed by power laws as quadratic or cubic stress–strain relations. The power laws (with power of 3/2, 2 or 3) of the stress–strain relation and the coefficients involved in the different terms of these laws are highly dependent on the shape of the pores or microcracks involved in the solid material (spherical, cylindrical, ellipsoidal, etc.; Nazarov & Sutin 1997; Nazarov 2001; Ostrovsky & Johnson 2001; Yu-Lin et al.2009). In this study, and without loss of generality, we will carry out the numerical verification using a quadratic term in the strain–stress relationship. Concerning the numerical modelling of seismic wave propagation in realistic media, general overviews (e.g. Carcione et al.2002; Moczo et al.2011, 2014) introduce most of the numerical methods dedicated to mechanical dynamics that include different direct methods such as pseudospectral, continuous (SEM)/discontinous (DG) Galerkin finite-element methods, and FD techniques in frequency or in time domains. 2-D or 3-D FD time domain (FDTD) approaches are commonly applied to discretize and solve the wave equations for elastic, viscoelastic, anisotropic and numerous dispersive media (Levander 1988; Robertsson et al.1994; Graves 1996, 1999; Bohlen 2002; Saenger & Bohlen 2004). These techniques have been generally applied using cartesian and staggered mesh discretizations. But more recently, more sophisticated FD techniques using high-order time-stepping and deformed meshes with non-staggered and collocated mesh point discretization (Zhang & Shen 2010; Zhang et al.2012; Sun et al.2018) have been introduced as well as non-deformed meshes using 2-D immersed methods (Lombard et al.2008; Chiavassa & Lombard 2013; Blanc et al.2014) with specific operators applied at physical interfaces and at the free surface (with mirror (Zhang & Chen 2006; Zhang et al.2012) or non-centred schemes (Lombard et al.2008)). Also, recent developments in including Zener body attenuation models have been made in the SEM (Blanc et al.2016) for a constant quality factor. Usually the linear viscoelastic behaviour is modelled by a memory variable formalism in order to ease the numerical implementation of the time convolution operator associated with the stress–strain relation (Moczo et al.1997; Xu & McMechan 1998; Hestholm 1999; Olsen et al.2000; Day 1998; Day & Bradley 2001; Wang et al.2001; Bohlen 2002; Saenger & Bohlen 2004). The excessive memory storage requirements associated to this material-independent memory variable approach needed to be numerically optimized and reduced by using coarse sampling approach for instance (Kristek & Moczo 2003). The aim of this study is to avoid the use of memory variables by extending an auxiliary differential equation (ADE) approach in three dimensions for a simple (quadratic term) but representative non-linear behaviour. To our knowledge, for the first time in 3-D seismic wave propagation modelling, the ADE formalism is used for implementing a non-linear and viscoelastic constitutive law. Such a formulation has been used with success for a linear viscoelastic case by Dhemaied et al. (2011) in an explicit FDTD-PML scheme using a second-order spatial and time stepping accuracy and a classic staggered grid for a single Zener body viscoelastic model. The authors were inspired by an ADE electromagnetic application in ground penetrating radar by Rejiba et al. (2003) and were motivated by extending the ADE approach to the seismic wave equations for a desired viscoelastic model. Attenuation is commonly described mathematically by different viscoelastic models defined as a combination of relaxation mechanisms (series and parallel settings of springs and dashpots) as in Moczo & Kristek (2005) that allow different frequency dependent formulations of the attenuation. Maxwell, Kelvin–Voigt or Zener body models are generally used, and the reader can be referred to Semblat & Pecker (2009) for more details about those models. The Zener body model is used here and provides a simple stress–strain formulation in the frequency domain and is defined as a Kelvin–Voigt cell and a spring in series. The use of ADE formalism in 3-D-FDTD, for non-linear and dispersive material, is widely used for electromagnetic waves propagation and particularly for photonics and active plasmonics applications where Lorentz–Drude, Raman and Kerr non-linear terms are incorporated (Taflove & Hagness 2005; Greene & Taflove 2006; Taflove et al.2013). This approach is very uncommon (almost inexistant) in seismic wave propagation despite its great potential, except in rare modelling studies in 1-D dispersive and attenuated fluid/acoustic biophysical media (Jiménez et al.2016) or non-destructive testing (Lombard & Piraux 2011). The only mention of ADE formalism concerns the implementation of ADE-perfectly matched layer (PML) absorbing boundary conditions (Martin et al.2010; Zhang & Shen 2010; Moczo et al.2014). The main specificity, and probably the best asset of this method, is that it allows a clear separation to be established between the set of propagation equations: the strain–velocity equation, the velocity–stress equation, and the rheological model defined by the stress–strain relation. Then by conviniently using the differentiation theorem for the Fourier transform, an inverse Fourier transformation from the frequency to the time domain can be performed for each term of the stress–strain relation. As for all FDTD techniques the ADE approach rigorously enforces the vector field boundary conditions at interfaces of different media associated to a small fraction of the impinging pulse width : consequently it is an almost completely general approach that permits the modelling of a broad variety of dispersive and even non-linear materials. As a consequence, this method has particularly the advantage of allowing the flexible implementation of any dispersion law based on non-constant Q quality factor law and on the complex expression of the desired viscoelastic modulus. More generally, by using the ADE approach, the stress–strain relations are not limited to expressions only based on linear dispersion laws or quasi-constant quality factors and can also incorporate non-linear rheological laws to take into account hysteretic, friction and damaging effects for instance. In this study, we will explain how the introduction of a strain–stress relation, involving non-linearity and attenuation, can be efficiently solved numerically in 3-D. We will present the numerical implementation and the numerical and physical verification of a simple Zener model associated with a quadratic non-linear term by introducing an ADE-FDTD approach. Our 3-D scheme includes a second order time-stepping and a fourth order staggered-grid spatial discretization to limit the numerical dispersion and convolutional (Komatitsch & Martin 2007; Martin & Komatitsch 2009; Martin et al.2008, 2010; Bodet et al.2014) or non-convolutional (Martin et al.2010) perfectly matched layer boundary conditions designed for elastic, viscoelastic or poroelastic seismic-wave modelling. We incorporate simultaneously these conditions along with viscoelastic properties inside the stress–strain formulations. It is worth noting that, in this study, non-convolutional PML using ADE formalism is implemented (Martin et al.2010) for more efficiency with less memory storage and subsequently faster computations. In the following, we first give the governing equations of the non-linear constitutive law used. For sake of simplicity and without loss of generality for linear or non-linear stress–strain relations, we consider a homogeneous medium and an associated rheology defined by a single Zener body model for the viscoelastic part of the stress–strain relation and an additional non-linear (quadratic) stress–strain relation term in one direction of space. The 3-D medium can thus be considered here as anisotropic in the direction of source polarization. We then show a set of 3-D numerical experiments to check for the accuracy and validity of the present numerical ADE-FDTD scheme, particularly concerning viscoelastic and non-linearity effects. The numerical experiments are designed for three different rheologies of the medium (elastic, viscoelastic and non-linear viscoelastic) to be compared if excited with the same source signal. We concentrate our numerical (i.e. spectral) analysis and verification process on the isotropic non-linear viscoelastic rheology. The study of non-linear aspects is performed thanks to classical spectral and dispersion analysis tools. Different amplitude factors and frequency contents of the source excitation are tested to clearly exhibit the non-linear response of the medium. The influence of the source type on the non-linear effects is also discussed. 2 GOVERNING EQUATIONS OF THE LINEAR AND NONLINEAR VISCOELASTIC EQUATIONS In this study, we aim at directly introducing linear or nonlinear viscoelastic rheologies, but with a true constitutive law relating the strain deformations to the stresses by writing, in the frequency domain: \begin{eqnarray*} \sigma &=& \sum _{i=1,N} \sigma _i^L \nonumber + \sigma _i^{NL} \nonumber \\ &=& \sum _{i=1,N} a_i \varepsilon _i + \sum _{i=1,N} \sum _{j=1,M} b_i c_{ij} \varepsilon _i^j , \end{eqnarray*} (1) where |$\sigma _i^L$| and |$\sigma _i^{NL}$| are, respectively, the linear and non-linear parts of the global effective tensor σ. ai are complex coefficients depending on the physical and geometrical properties of the linear viscoelastic part of the effective tensor. Subscript i refers to one of the partial tensors σi and j is the j-th polynomial degree of each term of the non-linear part of each i-th partial tensor σi . bi and cij are the polynomial coefficients of the non-linear part of the effective tensor. For simplicity, we will consider N = 1 in all cases, M = 1 in the linear case and M = 2 in the non-linear case, the general effective tensor calculation being possibly treated in the same way. The strong form of the linear viscoelastic wave propagation equation in the time domain is given by: \begin{eqnarray*} \rho \, \frac{\partial ^2}{\partial t^2} {\bf U}&=& \nabla \cdot \sigma ,\nonumber \\ \end{eqnarray*} (2) where U = (ui)i = 1, D (D being the space dimension) is the solid displacement vector and σ is the viscoelastic stress tensor. To introduce a viscoelastic rheology, the stress tensor is defined by a convolution product of the modulus tensor and the strain (Carcione 2007) as follows: \begin{eqnarray*} \sigma _{ij}=({\bf C}\!*\!\ {\bf \varepsilon} )_{ij} \end{eqnarray*} (3) \begin{eqnarray*} \varepsilon _{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right), \end{eqnarray*} (4) where C is the modulus tensor of the viscoelastic solid matrix. In the frequency domain, by applying a Fourier Transform, we can write: \begin{eqnarray*} \hat{\sigma }_{ij}&=&{\bf F}{({\bf C}\!*\!\ \bf \varepsilon )}_{ij}=(\hat{{\bf C}} \hat{\bf \varepsilon })_{ij}=\lambda ^* \delta _{ij} \hat{\varepsilon _{kk}} +2\mu ^* \hat{\varepsilon _{ij}} \nonumber \\ \hat{\varepsilon }_{ij}&=&\frac{1}{2}\left(\frac{\partial \hat{u}_i}{\partial x_j}+\frac{\partial \hat{u}_j}{\partial x_i}\right), \end{eqnarray*} (5) where indices i and j can be 1, 2 or 3 here in 3-D, and with the Einstein convention of implicit summation over a repeated index. ε is the strain tensor of the viscoelastic solid matrix, and ρ is the density of the solid material. |$\mu ^* = \mu (\frac{1+i \omega \tau _0 \gamma _{\mu }}{1+i \omega \tau _0})$| is a complex shear modulus and |$\lambda ^* = \lambda (\frac{1+i \omega \tau _0 \gamma _{\lambda }}{1+i\omega \tau _0})$| is a complex Lamé coefficient parameter of the solid matrix. We can notice that μ* and λ* are non-linear functions of the frequency ω. After some algebraic manipulations, we can also compare the numerical and the analytical compressible and shear moduli of the viscoelastic medium, the expressions of the analytical moduli being given by MP = λ* + 2μ* and MS = μ* and the expressions of the numerical moduli by \begin{eqnarray*} M_P&=& \frac{\hat{\sigma _{zz}}(\hat{\varepsilon _{zz}}+\hat{\varepsilon _{xx}})+\hat{\sigma _{yy}}(\hat{\varepsilon _{yy}}+\hat{\varepsilon _{xx}})}{\hat{\varepsilon _{zz}}(\hat{\varepsilon _{zz}}+\hat{\varepsilon _{xx}})+\hat{\varepsilon _{yy}}(\hat{\varepsilon _{yy}}+\hat{\varepsilon _{xx}})} \end{eqnarray*} (6) \begin{eqnarray*} M_S&=& \frac{\hat{\sigma _{yz}}}{\hat{\varepsilon _{yz}}} . \end{eqnarray*} (7) These expressions will allow us in the next section to validate the numerical solutions at least for the linear viscoelastic case. For simplicity, but without loss of generality, we add non-linear behaviour to the stress–strain law and consider that the non-linearity is only present in one direction of space (the longitudinal z direction for instance). We make the assumption that the fluid part is missing in the non-linear (quadratic) poroelastic formulation of Donskoy et al. (1997) or Dazel & Tournat (2010) and that the elastic part is replaced by the Zener linear viscoelastic strain–stress relation of Dhemaied et al. (2011). The stress–strain |$\hat{\sigma }(\omega )=M(\omega )\hat{\varepsilon }(\omega )$| relation in the frequency domain is given by: \begin{eqnarray*} \hat{\sigma }=\frac{(\lambda _{P,S}+\tilde{\lambda }_{P,S} i\omega )\hat{\varepsilon }+ H \hat{\varepsilon ^m}}{1+\tau _0 i\omega }, \nonumber \\ \end{eqnarray*} (8) where H is a non-linear tensor, m is an exponent that generally depends on the nature of the non-linear rheology (due to elliptical/spherical cracks, inclusions or cavities distributions, ...), τ0 is the stress relaxation time, λP, S and |$\tilde{\lambda }_{P,S}$| are, respectively, the unrelaxed (real) and relaxed (imaginary) parts of the Lamé parameters. Once these rheological laws are given in the frequency domain, they are expressed in the time domain after applying the inverse Fourier transform to each term of the stress–strain relation, and consequently the stress–strain equation yields: \begin{eqnarray*} \sigma + \tau _0 \partial _t \sigma = \lambda _{P,S} \varepsilon + \tilde{\lambda }_{P,S} \partial _t \varepsilon + H \varepsilon ^m ,\nonumber \\ \end{eqnarray*} (9) In a less compact form, the system of equations is developed as: \begin{eqnarray*} \sigma _{ij}+\tau _0 \frac{\partial \sigma _{ij}}{\partial t}&=&\lambda \delta _{ij} \varepsilon _{kk} +2\mu \varepsilon _{ij} +\tilde{\lambda } \delta _{ij} \partial _t \varepsilon _{kk} +2\tilde{\mu } \partial _t \varepsilon _{ij} + H_{ij} \delta _{ij} \varepsilon _{ij}^m \nonumber \\ \varepsilon _{ij}&=&\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \nonumber \\ \partial _t \varepsilon _{ij}&=&\frac{1}{2}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right), \end{eqnarray*} (10) where vi (i = 1, 3 in 3-D) are the solid velocity vector components. As stated earlier, we assume that Hxx = Hyy = Hxy = Hxz = Hyz = 0 and Hzz = D where the parameter D can take zero or non-zero value. In the case of D = 0 we retrieve the linear isotropic viscoelastic case and for D ≠ 0 we retrieve a non-linear anisotropic viscoelastic stress tensor in the frequency domain with a non-linear anisotropy in one direction of space. A power m = 2 is chosen now in this study for the additional non-linear term without loss of generality. The relaxed Lamé coefficients |$\tilde{\lambda }$| and |$\tilde{\mu }$| are expressed as follows: \begin{eqnarray*} \tilde{\lambda } &=&\tau _0 \lambda \gamma _{\lambda } \nonumber \\ \tilde{\mu } &=&\tau _0 \mu \gamma _{\mu } \end{eqnarray*} (11) with \begin{eqnarray*} \gamma _{\lambda }&=& \frac{1}{\lambda }\left(\frac{\tau _p}{\tau _0}(\lambda + 2\mu ) -2 \mu \gamma _{\mu }\right) \nonumber \\ \gamma _{\mu } &=& \frac{\tau _s}{\tau _0} , \end{eqnarray*} (12) where τp and τs are the P and S-wave relaxation times. Eq. (10) is the required ADE defining the stress tensor σ(t) in the time domain. It has been easily implemented in an FDTD code using a semi-implicit scheme centred at time step n + 1 wherein yet-to-be computed velocity and displacement fields at time n + 1/2 are used to create an update formula for a field known at time-step n; almost exactly as detailed in Dhemaied et al. (2011) for the 2-D case. This ADE-FDTD integration scheme is detailed in the next section. 3 NUMERICAL EXPERIMENT 3.1 Setup A second order displacement formulation is chosen and a fourth order staggered-grid FD spatial discretization in space (Graves 1996) as well as a second order time-stepping algorithm are used to solve the whole set of equations. The ADE system of equations is discretized in time at the second order using a semi-implicit scheme as follows: \begin{eqnarray*} \rho \, v_i^{n+1/2} &=& \rho \, v_i^{n-1/2} + \Delta t (\partial _j \sigma _{ij}^{n} +F_i ) \nonumber \\ u_i^{n+1/2} &=& u_i^{n-1/2} +\Delta t v_i^{n+1/2} + \Delta t^2 \frac{v_i^{n+1/2}-v_i^{n-1/2}}{\Delta t} \nonumber \\ \varepsilon _{ij}^{n+1/2}&=&\frac{1}{2}\left(\frac{\partial u_i^{n+1/2}}{\partial x_j}+\frac{\partial u_j^{n+1/2}}{\partial x_i}\right) \nonumber \\ \partial _t \varepsilon _{ij}^{n+1/2}&=&\frac{1}{2}\left(\frac{\partial v_i^{n+1/2}}{\partial x_j}+\frac{\partial v_j^{n+1/2}}{\partial x_i}\right) \nonumber \\ \sigma _{ij}^{n+1}(\tau _0 /\Delta t + 0.5)&=& \sigma _{ij}^n (\tau _0 /\Delta t - 0.5)+ \lambda \delta _{ij} \varepsilon _{kk}^{n+1/2} +2\mu \varepsilon _{ij}^{n+1/2} \nonumber \\ &=& +\tilde{\lambda } \delta _{ij} \partial _t \varepsilon _{kk}^{n+1/2} +2\tilde{\mu } \partial _t \varepsilon _{ij}^{n+1/2} \nonumber \\ &=& + H_{ij} \delta _{ij} \varepsilon _{ij}^2\mid ^{n+1/2} , \end{eqnarray*} (13) where ui and vi (i = 1, 3 in 3-D) are the, respectively the solid displacement and velocity vector components. We implement the following FD operator of each gradient of a variable f in a given direction x, similar operators being applied in the other two directions y and z: \begin{eqnarray*} \frac{\partial f}{\partial x}\mid _{i+1/2,j,k}\equiv \frac{-f_{i+1,j,k}+27 f_{i,j,k}-27 f_{i-1,j,k}+f_{i-2,j,k}}{24} . \end{eqnarray*} (14) Optimized convolution (Komatitsch & Martin 2007) or non-convolution (Martin et al.2010) frequency shift perfectly matched layers absorbing boundary conditions complete the implementation to efficiently truncate the computational domain and mimic this way an infinite medium. The advantage of this non-convolution PML is that we have more flexibility to implement different time-stepping schemes when compared to CPML formulation. For the sake of clarity, we use here, respectively, second and fourth order discretization operators in time and space to integrate the equations as well as CPML boundary conditions. Besides, without loss of generality, higher order operators can also be used to solve the elastodynamic equations and add non-convolution PML conditions as in Martin et al. (2010), where we introduced eighth order Holberg space operators and fourth-order Runge–Kutta (RK4) time-stepping integration. The discretization could also be replaced by more sophisticated space discretizations in space and time on deformed meshes and collocated grids as those mentionned in the introduction. But here, the fourth-order staggered grid numerical discretization and second-order time-stepping we choose are enough for the purpose of this study because we just want to exhibit attenuation and non-linear effects in 3-D with reasonable accuracy. In addition, it should be emphasized that, to regularize possible sharp shocks related to the non-linearity introduced, we use non-centred fourth-order spatial discretization of the quadratic term ε2 involved in the strain–stress law and applied in the source polarization direction z only (here only the component Hzz = D is non-zero in the H non-linear tensor, and εzz is denoted by ε to alleviate notations). To avoid possible high order spurious modes, we denote by |$\overline{\varepsilon }_{i,j,k+1/2}=0.5(\varepsilon ^L_{i,j,k}+\varepsilon ^R_{i,j,k})$| the average strain computed at a cell face discontinuity where two strain deformation states L (left) and R (right) are defined. We assume that |$\varepsilon ^2_{i,j,k+1/2}=\overline{\varepsilon }_{i,j,k+1/2} \varepsilon ^{L}_{i,j,k+1/2}$| if |$\overline{\varepsilon }_{i,j,k+1/2} \gt 0$|⁠, and |$\varepsilon ^2_{i,j,k+1/2}=\overline{\varepsilon }_{i,j,k+1/2} \varepsilon ^R_{i,j,k+1/2}$| if |$\overline{\varepsilon }_{i,j,k+1/2} \lt 0$|⁠. Here we choose |$\varepsilon ^L_{i,j,k+1/2/}=\varepsilon _{i,j,k}$| and |$\varepsilon ^R_{i,j,k+1/2}=\varepsilon _{i,j,k+1}$|⁠. A simplified illustration of our numerical setup is presented on Fig. 1(a). In all cases, we perform our simulations using a dominant frequency of 1500 Hz for a Ricker wavelet point-force source, implemented in the longitudinal (z) direction (see Figs 1a and b). We have sampled the spatial and time steps in our simulations in order to have access to a wide range of wavelengths compared to the source signal. The time step is thus chosen as Δt =0.4 μs and the spatial step is Δx = Δy = Δz = 0.0014 m. The grid mesh is discretized by 1800 points along the longitudinal z-axis and, 200 points over the traversal x direction and 200 points over the depth y. A source is located at |$z=0.2324\, \mathrm{m}$| from the left boundary of the thin slice and receivers are located at all discretization points along the longitudinal length and in the middle of the medium (the length of the models along the longitudinal axis z is 2.52 m). We checked the accuracy of the numerical method for pure viscoelastic and non-linear viscoelastic cases for two different sources (a Ricker and a Gaussian-modulated sine pulse source) that will be used in the next sections of this study. For sake of clarity, the reader is referred to the appendix section for more details. In that section P- and S-wave velocities are very well reproduced for both types of sources in the viscoelastic case. Furthermore, in the non-linear viscoelastic case, the σzz component of the stress tensor (computed using the non-centred discretization of the non-liner term at the cell faces described before) is well reproduced in the frequency domain. Figure 1. View largeDownload slide Simplified illustration of our numerical setup (a). A dominant frequency of 1500 Hz for a Ricker wavelet point-force source has been implemented in the longitudinal (z) direction with several excitation amplitudes in the non-linear viscoelastic case (black diamonds correspond to picked maxima) (b). Figure 1. View largeDownload slide Simplified illustration of our numerical setup (a). A dominant frequency of 1500 Hz for a Ricker wavelet point-force source has been implemented in the longitudinal (z) direction with several excitation amplitudes in the non-linear viscoelastic case (black diamonds correspond to picked maxima) (b). 3.2 Numerical verification of the nonlinear viscoelastic regime 3.2.1 Classical dispersion analysis using a wideband pulse Now, we focus our attention to the physical aspects of our study. To illustrate the effect of the current non-linear regime (Fig. 2a), we first present the seismograms obtained for a reference elastic and homogeneous medium characterized by a pressure-wave velocity (Vp) of 245.63 |$\mathrm{m\, s}^{-1}$|⁠, a shear-wave velocity (Vs) of 122.65 |$\mathrm{m\, s}^{-1}$| and a bulk density |$\rho =1610\, \mathrm{kg\, m}^{-3}$|⁠, these values corresponding to experimental observations in unconsolidated granular media recently studied in Bodet et al. (2014). A second simulation considers the viscoelastic case, based on similar homogeneous properties but with stress and strain relaxation times, respectively, set to τ0 =0.28 ms and τp = τs =0.3 ms. A third simulation finally presents the non-linear case. The D parameter described earlier is set to |$5 \times 10^{13} \times \rho V_p^2$| to highlight non-linear regime effects. Figure 2. View largeDownload slide Seismograms (a) and spectrograms (b) of vertical component particle velocity (raw amplitude) recorded along the simulated models for: elastic, linear viscoelastic and non-linear viscoelastic rheologies medium (source excitation amplitude × 1). Figure 2. View largeDownload slide Seismograms (a) and spectrograms (b) of vertical component particle velocity (raw amplitude) recorded along the simulated models for: elastic, linear viscoelastic and non-linear viscoelastic rheologies medium (source excitation amplitude × 1). The seismograms of particle velocity recorded in the z direction are given as a function of offset (source-recording point distance along the longitudinal z-axis) and time (presented on Fig. 2a). They consist in 1634 traces of 50 000 amplitude samples. Except for visible variations in amplitude of particle velocity between the elastic case and the two viscoelastic cases, distance–time data systematically depict two distinct events: a P- and a S-wave train (P and S on Fig. 2a) of similar apparent velocities for each rheology. However, the effect of viscoelasticity is clearly apparent on corresponding spectrograms (Fig. 2b) which present important differences in terms of amplitude and shape. The non-linear case typically presents frequency-up conversions with at least two amplifications at higher frequencies that are linked to the source peak frequency (⁠|$1.5\, \mathrm{kHz}$|⁠) but could not be exactly and accurately associated to a specific harmonic: the fundamental mode is located around a frequency of 1900 Hz, and the third harmonic is roughly 5700 Hz, but the second harmonic that should be located around a frequency of 3800 Hz does not appear clearly. In addition, the wavefield recorded for each case has been transformed into the frequency-phase velocity domain (by a slant-stack in common shot gathers after correction for geometrical spreading) to quantify the phase velocities of the P- and S-wave events cited above. The resulting dispersion curves, presented on Fig. 3 clearly show that the simulations perfectly describe implemented phase velocities: constant velocity in the elastic case and, by definition, dispersive relationship between velocity and frequency in the viscoelastic case. Figure 3. View largeDownload slide Dispersion curves: maxima picked from frequency-phase velocity transforms of seismograms for both P and S waves compared to corresponding theoretical dispersion curves. Figure 3. View largeDownload slide Dispersion curves: maxima picked from frequency-phase velocity transforms of seismograms for both P and S waves compared to corresponding theoretical dispersion curves. 3.2.2 P-wave extraction using advanced preprocessing tools and pulse modulation We will consider the Ricker source described earlier as a reference in term of amplitude. More runs involving exactly the same parameters as the latter non-linear viscoelastic model are presented, but with different amplitudes of the source (from 0.25 × to 2.25 × the excitation amplitude of the first case, as shown on Fig. 1b). We will focus on the longitudinal z direction, and consequently on the P-wave train only. The seismograms, presented on Fig. 4 obtained for each source excitation amplitude, are thus muted as soon as P- and S-wave trains can be separated. Each corresponding spectrogram on Fig. 5 clearly depicts the harmonics observed earlier resulting from implemented non-linearities. Figure 4. View largeDownload slide Seismograms of particle velocity (raw amplitude) recorded along the simulated models for different source excitation amplitudes. S-wave muted. Figure 4. View largeDownload slide Seismograms of particle velocity (raw amplitude) recorded along the simulated models for different source excitation amplitudes. S-wave muted. Figure 5. View largeDownload slide Amplitude spectrum of particle velocity (raw amplitude) recorded along the simulated models for different source excitation amplitudes. S-wave muted. Figure 5. View largeDownload slide Amplitude spectrum of particle velocity (raw amplitude) recorded along the simulated models for different source excitation amplitudes. S-wave muted. In order to clearly identify the non-linear effects, we extracted traces from each seismogram and spectrogram along the line at four different offsets, corresponding to two to five times the dominant wavelength of the source (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠). These traces are presented on top of Fig. 6. They correspond to the signals computed at each different source excitation amplitude (from × 0.25 to × 2.25 the amplitude of the reference run presented in Section 3.2.1; see colour-scale on Fig. 1b). The corresponding spectra are presented on top of Fig. 7, on which the main frequency peak A1 and two other peaks A′ and A″ clearly appear (maxima are picked and shown as black plus signs). These maxima are not in ratio 2 in frequency, but seem to correspond to the third and fifth harmonics, respectively, which poses the question of their correspondence with harmonic wave components. These traces are then normalized to the source maximum amplitude in both time and frequency domains (bottom lines of Figs 6 and 7). When the normalization in time remains equivalent (every normalized traces can be considered superimposed on Fig. 7), the ratio presented in frequency clearly shows opposite behaviour as a function of source amplitude depending on the considered peak. The black arrows, given on top of Fig. 7 as an example, depict a non-monotonous behaviour of these harmonics as a function of source amplitude. Figure 6. View largeDownload slide Traces recorded at four different offsets given as a function of the source-signal dominant wavelength (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠) presented for each different source excitation amplitude (from × 0.25 to × 2.25; see colour-scale on Fig. 1b) with corresponding normalization to the amplitude at the source (A0). A is the amplitude of the trace at corresponding offset. Figure 6. View largeDownload slide Traces recorded at four different offsets given as a function of the source-signal dominant wavelength (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠) presented for each different source excitation amplitude (from × 0.25 to × 2.25; see colour-scale on Fig. 1b) with corresponding normalization to the amplitude at the source (A0). A is the amplitude of the trace at corresponding offset. Figure 7. View largeDownload slide Amplitude spectra and spectral ratio of traces recorded at four different offsets given as a function of the source-signal dominant wavelength (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠) presented for each different source excitation amplitude (from × 0.25 to × 2.25; see colour-scale on Fig. 1b). A is the amplitude of the trace at corresponding offset, A0 the amplitude at source and A1, A′ and A″ are the amplitudes of the three first maxima detected in the spectrum at each offset. The black arrows are just given on top as an example showing their non-monotonous behaviour as a function of source amplitude. Figure 7. View largeDownload slide Amplitude spectra and spectral ratio of traces recorded at four different offsets given as a function of the source-signal dominant wavelength (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠) presented for each different source excitation amplitude (from × 0.25 to × 2.25; see colour-scale on Fig. 1b). A is the amplitude of the trace at corresponding offset, A0 the amplitude at source and A1, A′ and A″ are the amplitudes of the three first maxima detected in the spectrum at each offset. The black arrows are just given on top as an example showing their non-monotonous behaviour as a function of source amplitude. We picked the amplitude of each peak and presented, on Fig. 8, their ratios between the first peak and the other peaks. These spectral ratios were computed at the 4 different offsets corresponding to two to five times the dominant wavelength of the source (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠) described above. The colours on Fig. 8 depend on these offsets and ratio are given of course for every source excitation amplitudes (from × 0.25 to × 2.25 the amplitude of the reference run presented in Section 3.2.1). At far-field offsets, the ratio shows strong variations with source amplitude and tends to remain stable at near-field offsets. As stated earlier, picked maxima are not in ratio 2 in frequency which suggests overlap (interferences) between adjacent harmonics due to the wide frequency band of the Ricker source and dispersion effects (a 1900 Hz fundamental frequency and a 5700 Hz third harmonic frequency). This observation illustrates as well the fact that non-linear effects accumulate over the distance to the source, as these ratios increase at fixed amplitude A1 and increasing distance. Figure 8. View largeDownload slide Spectral ratio computed at different offsets for every Ricker wavelet source excitation amplitudes (the four different offsets, corresponding to two to five times the dominant wavelength of the source (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠), are given with colours depending offset from blue (2λP) to green (5λP)). Figure 8. View largeDownload slide Spectral ratio computed at different offsets for every Ricker wavelet source excitation amplitudes (the four different offsets, corresponding to two to five times the dominant wavelength of the source (⁠|$\lambda _P \approx 0.18\, \mathrm{m}$|⁠), are given with colours depending offset from blue (2λP) to green (5λP)). Consequently, to exhibit and individualize the different peaks more precisely, we introduce a narrow frequency band pulse source implemented as a Gaussian-modulated sine pulse with a dominant frequency of 1500 Hz. Fig. 9 compares normalized traces obtained from both the Ricker pulse (with amplitude of the source 2 × the excitation amplitude of the first case (3.2.1)) and the Gaussian-modulated sine pulse (Fig. 10). From these seismograms (Figs 10a and b), we computed amplitude (Figs 10c and d) and f−k spectra (Figs 10e and f). The spectra clearly show that, within the same propagation ranges and domains, the source modulation allows the identification of all consecutive harmonics (frequency-up conversions in ratio 2 in frequency), corresponding to the expected non-linear effect. In addition, for example in the case of strong non-linearities arising from a high amplitude of the source, the superimposition of all spectra corresponding to all successive offsets along the analysed direction (Fig. 11) highlights the overlapping process associated with wideband pulse. In the present case, all first six successive modes are clearly defined in terms of peak frequency location and amplitude. In this case of a Gaussian-modulated sine pulse, we can observe that the amplitude maxima A1 have similar non-linear effects (Figs 12a and b) as the A1 maxima in the Ricker wavelet source case (Figs 11a and b), and are much better individualized on a wider frequency range. However, the spectral ratio amplitudes A3/A13 have almost the same exponentially increasing patterns as A″/A13 in the Ricker source case at almost similar offsets. In Fig. 11(b), the second and third peaks A′ and A″ are well separated while in Fig. 12(b) the maxima A2 and A3 are almost superimposed. The ratios of second and third harmonic amplitudes A2 and A3 over first peaks A1, A12 and A13 present a more pronounced non-linear behaviour (Figs 12b–d) at offsets closer to the source location when compared to the Ricker wavelet source case. More particularly, the ratios of second and third harmonic amplitudes A2 and A3 over A12 and A13 show an increasing trend with offset (Figs 12c and). Small oscillations of those spectral ratios for the gaussian-modulated sine pulse appear at small offsets probably due to near field effects, and seem to be mainly associated with amplitude oscillations of A1. Finally, in both non-linear cases, we can see that the spectral ratios are increasing (Figs 11e and f, 12e and f) and are reaching similar high values for values of A1 close to zero at very far offsets to the source while they are converging asymptotically to a constant value for high values of A1 close to the source. Figure 9. View largeDownload slide Test comparing normalized traces obtained from the Ricker pulse and a gaussian-modulated sine pulse: (a) at the source location and (b) at near offset (2λP). Figure 9. View largeDownload slide Test comparing normalized traces obtained from the Ricker pulse and a gaussian-modulated sine pulse: (a) at the source location and (b) at near offset (2λP). Figure 10. View largeDownload slide Seismograms (a, b), spectrograms (c, d) and f−k spectra (e, f) obtained for the Ricker pulse source (a, c, e) and the Gaussian-modulated sine pulse source (b, d, f) are compared to show the impact of two different types of pulse source on the seismograms and the spectra. The source amplitude factor is equal to 2 in both pulse source cases. Figure 10. View largeDownload slide Seismograms (a, b), spectrograms (c, d) and f−k spectra (e, f) obtained for the Ricker pulse source (a, c, e) and the Gaussian-modulated sine pulse source (b, d, f) are compared to show the impact of two different types of pulse source on the seismograms and the spectra. The source amplitude factor is equal to 2 in both pulse source cases. Figure 11. View largeDownload slide Offset presentation (a) of spectrograms of Fig. 10 for the Ricker pulse source. In semi-logarithmic scale: (b) Representation of the amplitudes of the three first maxima A1, A′ and A″ detected in the spectrum at all offsets in the simulation. (c) and (d) represent spectral ratio between the three first maxima computed at all offsets in (c) and (d), and versus the first peak A1 in (e) and (f). Figure 11. View largeDownload slide Offset presentation (a) of spectrograms of Fig. 10 for the Ricker pulse source. In semi-logarithmic scale: (b) Representation of the amplitudes of the three first maxima A1, A′ and A″ detected in the spectrum at all offsets in the simulation. (c) and (d) represent spectral ratio between the three first maxima computed at all offsets in (c) and (d), and versus the first peak A1 in (e) and (f). Figure 12. View largeDownload slide Offset presentation (a) of spectrograms of Fig. 10 for the Gaussian-modulated sine pulse source. In semi-logarithmic scale : (b) Representation of the amplitudes of the three first maxima A1, A2 and A3 detected in the spectrum at all offsets in the simulation. (c) and (d) represent spectral ratio between the three first maxima computed at all offsets in (c) and (d), and versus the first peak A1 at all offsets in (e) and (f). Figure 12. View largeDownload slide Offset presentation (a) of spectrograms of Fig. 10 for the Gaussian-modulated sine pulse source. In semi-logarithmic scale : (b) Representation of the amplitudes of the three first maxima A1, A2 and A3 detected in the spectrum at all offsets in the simulation. (c) and (d) represent spectral ratio between the three first maxima computed at all offsets in (c) and (d), and versus the first peak A1 at all offsets in (e) and (f). 3.3 Discussion on non-linear regimes using non-dimensional analysis In order to characterize better the regime (linear, weakly or strongly non-linear) under which the numerical experiments studied in the previous section are performed, we consider below the simplest non-linear model and its associated non-dimensional numbers. The simulated problem corresponds to an oscillating point (dipolar) source along the z-direction in a homogeneous elastic solid medium, showing linear viscoelastic isotropy. The non-linearity, as described earlier, is anisotropic (only acting on the σzz and εzz relation) and quadratic. The non-linear propagation is inspected along the z-direction. To our knowledge, there is no analytical solution in this configuration for the non-linear process of harmonic generation. We therefore recall here the usual parameters for plane wave propagation of longitudinal waves in a homogeneous medium with quadratic non-linearity. First, we define linear and non-linear time or length scales denoted by τa, τnl, La, Lnl. La and τa correspond to the characteristic length and timescales over which the initial amplitude of a wave is divided by e due to dissipation effects, and Lnl and τnl are the characteristic non-linear length and timescales over which non-linear effects develop. If the receivers are located close to the source then non-linearities cannot accumulate enough. If they are at a distance to the source lower than Lnl then the non-linear effects can be observed but are still in a regime of accumulation. If it is equal to Lnl then the non-linearities are fully accumulated and developped and Lnl denotes the shock front distance. If we non-dimensionalize the non-linear equations, they can be expressed as: \begin{eqnarray*} \omega ^2 \rho U_o {\bf U}^{\prime } &=& 1/L_a div( \frac{(\lambda _{P,S}+\tilde{\lambda }_{P,S} i\omega ) U_o/L_a \hat{\varepsilon ^{\prime }}}{1+\tau _0 i\omega }+ D U_o^2/L_a^2 \hat{\varepsilon ^{\prime 2}} ) \nonumber \\ {\bf U}^{\prime }&=& \frac{(\lambda _{P,S}+\tilde{\lambda }_{P,S} i \omega ) U_o / L_a^2}{\omega ^2 \rho U_o (1+\tau _0 i \omega )} div(\hat{\varepsilon ^{\prime }}) + \frac{D U_o^2}{\rho \omega ^2 U_o L_a^3} div(\hat{\varepsilon ^{\prime 2}}) \nonumber \\ \end{eqnarray*} (15) where U = U′ × Uo, |$x_i = x_i^{\prime } \times L_a$|⁠, t = t′/ω, ε = ε′ × Uo/La with ε′ = 0.5(∇′U′ + ∇′U′T). If we denote the modulus of the viscoelastic seismic velocity by |$Vp_a = (\frac{\mid \mid M_P\mid \mid }{\rho })^{1/2}$|⁠, where MP is the complex modulus associated to the P wave, then the (visco-)elastic length scale can be considered homogeneous to the seismic wavelength La = Vpa*T = Vpa/ω. Then the previous equations can be simplified as: \begin{eqnarray*} {\bf U}^{\prime } &=& \frac{\rho Vp_a^2 U_o \omega ^2/Vp_a^2}{\omega ^2 \rho U_o} div(\hat{\varepsilon ^{\prime }}) + \frac{D U_o^2 L_a}{\rho \omega ^2 U_o L_a^4} div(\hat{\varepsilon ^{\prime 2}}) \nonumber \\ &=& div(\hat{\varepsilon ^{\prime }}) + \frac{D U_o^2 L_a \omega ^4}{\rho \omega ^2 U_o Vp_a^4} div(\hat{\varepsilon ^{\prime 2}}) \nonumber \\ &=& div(\hat{\varepsilon ^{\prime }}) + \frac{D U_o L_a \omega ^2}{\rho Vp_a^4} div(\hat{\varepsilon ^{\prime 2}}) \nonumber \\ &=& div(\hat{\varepsilon ^{\prime }}) + \frac{L_a}{L_{nl}} div(\hat{\varepsilon ^{\prime 2}}) \nonumber \\ &=& div(\hat{\varepsilon ^{\prime }}) + \Gamma div(\hat{\varepsilon ^{\prime 2}}), \end{eqnarray*} (16) where we can consider the non-dimensional Goldberg number Γ that measures the ratio between the linear attenuation effects and the non-linear elastic effects. Γ can be expressed as \begin{eqnarray*} \Gamma =\frac{\tau _a}{\tau _{nl}}=\frac{L_a}{L_{nl}} , \end{eqnarray*} (17) where the different lengths and timescales can be formulated as: \begin{eqnarray*} L_{nl} &=&\frac{\rho Vp_a^4}{\omega ^2 D U_o} \nonumber \\ \tau _{nl} &=&\frac{L_{nl}}{Vp_a} \nonumber \\ &=&\frac{\rho Vp_a^3}{\omega ^2 D U_o} \nonumber \\ \tau _a &=& \Gamma \tau _{nl} \nonumber \\ L_a &=& Vp_a /\omega \nonumber \\ L_a &=& L_{nl} \frac{\tau _a}{\tau _{nl}} . \end{eqnarray*} (18) Since we have |$D= D^{\prime } \rho V_p^2$| with D′ = 5.1013 we can compute Γ: \begin{eqnarray*} \Gamma &=& \frac{D^{\prime } U_o L_a \omega ^2}{Vp_a^2}\nonumber \\ \Gamma &=&D^{\prime } \frac{V_o}{Vp_a} . \end{eqnarray*} (19) In our simulations of the previous Section 3.2.2 and for a source amplitude factor of 2.25, the particle velocity and the displacement at the source are Vo = 2.1 × 10−11 m.s−1 and Uo = 1.4 × 10−14 m. The values of these time and length scales are τa ≃ 0.3 ms (maximum viscoelastic relaxation time), τnl ≃ 0.07714 ms, La ≃ 0.18 m (maximum viscoelastic wavelength), Lnl ≃ 0.04628 m. The non-dimensional Goldberg number is approximately equal to a value of ≃4. Furthermore, the non-linear effects can be measured also in terms of the rate of perturbation of the linear stress tensor by the non-linear terms in the P-wave direction: |$\frac{\sigma _{zz}^{nl}- \sigma _{zz}^{l}}{\sigma _{zz}^{l}}=\frac{E^* \varepsilon _{zz}(1+ D/E^* \varepsilon _{zz})}{E^* \varepsilon _{zz}}=D/E^* \varepsilon _{zz}$|⁠, where |$E^* = \rho V_p^2$|⁠. Here D/E* = D′ = 5 × 1013. Therefore if we take the initial strain εzz0 = 7.77 × 10−14 measured close to the source we have a ratio value around 3.88. When compared to plane waves, in the 3-D simulated problem it is expected that La is reduced due to geometrical attenuation, and Lnl is increased due to decreasing wave amplitude along distance. As a result, the Goldberg number of our problem (the non-linear over linear ratio) is necessarily slighlty smaller than 4, and exhibits clearly the non-linear propagation regime. The non-linear effects on signal amplitudes depicted in Figs 11 and 12 correspond to an amplitude source equal to 2.25 and to a very similar Goldberg non-linear number. If the wave problem would be for plane waves, the configuration of viscoelastic parameters and the non-linear parameter would give a Goldberg number of ≃4. This indicates that non-linear effects would appear in the direction of propagation of the viscoelastic plane wave. However, in our specific problem, the three-dimensional character of the oscillating point source introduces additional effects. First, a higher attenuation with distance from the source is expected due to geometric losses. This leads to a localization of the main non-linear effects closer to the source than in the plane wave case. This also allows for the presence of near field effects such as a non-monotonous evolution of the fundamental wave amplitude along the z-axis, probably at the origin of the oscillations observed in Fig. 12. In our configuration, the fact that the non-linearity operates only on εzz, and is essentially anisotropic, necessarily leads to a preferential accumulation direction for the non-linear effects of harmonic generation. The harmonic generation process is thus possible only in the z-direction, which may lead to the generation of a 'beam' of harmonics with less geometrical attenuation than the fundamental wave. This unusually strong non-linear anisotropy, which, to our knowledge, has not been observed in real systems earlier, is only made possible here thanks to the developed simulation method. 4 CONCLUSIONS The simulation of seismic wave propagation in a 3-D viscoelastic and non-linear media is presently performed using a 3-D ADE-FDTD scheme. Comparing to classical recursive convolution schemes, the flexibility and its perfect suit for multi-threaded execution of the ADE method must be weighted with the need of implementing as many functions as constitutive laws to be tested. However it offers the possibility to even define 3-D models made of completely different non-linear behaviours. The choice of the non-linear quadratic term in the constitutive law is adapted to the granular media behaviour under consideration. Such constitutive law represents fairly well the non-linear behaviour due to grain-to-grain interaction which plays a significant role in earthquake strong ground motion modelling. In the present case, the non-linear constant value D, associated with a Gaussian-modulated sine pulse peak frequency of 1.5 kHz, generates harmonics that are multiple of 1.5 kHz. The non-linearities detected in the displacement spectrograms highlight several harmonics at frequencies that are multiples of the dominant frequency of the excitation source only for narrow band pulse excitations. In case of a wide band pulse, (i.e. when the peak frequency is similar to the bandwidth) the different harmonics are not so well described and individualized, and are overlapping each other. Such aspects are rarely discussed, more likely because such wideband sources are rarely encountered in reality whatever the scale and type of the case study (i.e. seismology, near-surface seismic, non-destructive testing). The results show that we are able to accurately assess non-linear and dispersive behaviour in 3-D with an adequate strategy of source excitation, that is going to a Gaussian-modulated sine wave of narrower frequency bandwidth, which constitutes an efficient tool to tackle for example seismic source inversion problems. Moreover, in three-dimensions the orders of magnitude of the amplitudes simulated are lower than their equivalent in 1-D, which make very difficult their interpretation when observed in rare occasions as in laboratory sandbox experiments (Bodet et al.2014). It is worth noting, that the numerical verification process of the current 3-D scheme has been performed considering a non-linear viscoelastic behaviour along a single direction, while others were kept simply associated with a viscoelastic model. It does not show any sign of numerical instability, and thus illustrates also the capabilities of this numerical approach to appraise a whole range of issues associated to the role of anisotropy associated with both dispersion and non-linear behaviour (Saenger & Bohlen 2004). Finally, we have defined non-linear non-dimensional numbers like the Goldberg number and the ratio of non-linear perturbations to linear terms that are representative of the non-linear regimes. The validity of our 3-D ADE-FDTD approach has been illustrated on an anisotropic analog rather than on a heterogeneous model as the main objective of defining a dispersive-non-linear constitutive law is to model the behaviour associated with a naturally heterogeneous medium at different scales. In the future an interesting goal would consist in using such homogenized models (for dispersion and non-linear behaviour) for several types of problems: at the mesoscale like site effect issues (for example by adjusting the constitutive law of the sedimentary filling) or at the laboratory scale for phenomenological issues to analyse the attenuation parameters coming from the data analysis. At the laboratory scale we could then try to see if non-linear effects due to grain-to-grain contacts can explain acoustic modes appearing on wider frequency ranges at different distances from the source and how they are competing with attenuation effects. ACKNOWLEDGEMENTS This work was granted access to the high-performance computing resources of the French supercomputing centre CCRT under allocation #2016-046351 and #2017-046351 awarded by GENCI (Grand Équipement National de Calcul Intensif) and of the regional computing centre CALMIP of Toulouse, France, under allocation #P1135. REFERENCES Akkar S. , Sandikkaya M.A. , Ay B.Ö. , 2014 . Compatible ground motion prediction equations for damping scaling factors and vertical-to-horizontal spectral amplitude ratios for the broader Europe region , Bull. Earthq. Eng. , 12 ( 1 ), 517 – 547 . Google Scholar Crossref Search ADS Berjamin H. , Lombard B. , Chiavassa G. , Favrie N. , 2017 . Analytical solution to the 1D non-linear elastodynamics with general constitutive laws , Wave Motion , 74 , 35 – 55 . Google Scholar Crossref Search ADS Biot M.A. , 1956a . Theory of propagation of elastic waves in a fluid-saturated porous solid. I: low-frequency range , J. acoust. Soc. Am. , 28 , 168 – 178 . Google Scholar Crossref Search ADS Biot M.A. , 1956b . Theory of propagation of elastic waves in a fluid-saturated porous solid. II: higher-frequency range , J. acoust. Soc. Am. , 28 , 179 – 191 . Google Scholar Crossref Search ADS Biot M.A. , 1973 . Non-linear and semilinear rheology of porous solids , J. geophys. Res. , 23 , 4924 – 4937 . Google Scholar Crossref Search ADS Blanc E. , Chiavassa G. , Lombard B. , 2014 . Wave simulation in 2D heterogeneous transversely isotropic porous media with fractional attenuation: a Cartesian grid approach , J. Comput. Phys. , 275 , 118 – 142 . Google Scholar Crossref Search ADS Blanc E. , Komatitsch D. , Chaljub E. , Lombard B. , Xie Z. , 2016 . Highly accurate stability-preserving optimization of the Zener viscoelastic model, with application to wave propagation in the presence of strong attenuation , Geophys. J. Int. , 205 ( 1 ), 427 – 439 . Google Scholar Crossref Search ADS Bodet L. , Dhemaied A. , Martin R. , Mourgues R. , Rejiba F. , Tournat V. , 2014 . Small-scale physical modeling of seismic-wave propagation using unconsolidated granular media , Geophysics , 79 ( 6 ), T323 – T339 . Google Scholar Crossref Search ADS Bohlen T. , 2002 . Parallel 3-d viscoelastic finite difference seismic modelling , Comput. Geosci. , 28 ( 8 ), 887 – 899 . Google Scholar Crossref Search ADS Bonilla L.F. , Archuleta R.J. , Lavallée D. , 2005 . Hysteretic and dilatant behavior of cohesionless soils and their effects on non-linear site response: field data observations and modeling , Bull. seism. Soc. Am. , 95 ( 6 ), 2373 . Google Scholar Crossref Search ADS Bonilla L.F. , Tsuda K. , Pulido N. , Régnier J. , Laurendeau A. , 2011 . Non-linear site response evidence of Knet and KiKnet records from the 2011 off the pacific coast of tohoku earthquake , Earth, Planets Space , 63 ( 7 ), 50 . Google Scholar Crossref Search ADS Carcione J.M. , 2007 . Wave Fields in Real Media: Theory and Numerical Simulation of Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media , 2nd edn , Elsevier Science . Carcione J.M. , Herman G.C. , Ten Kroode A. , 2002 . Seismic modeling , Geophysics , 67 ( 4 ), 1304 – 1325 . Google Scholar Crossref Search ADS Chiavassa G. , Lombard B. , 2013 . Wave propagation across acoustic / Biot’s media: a finite-difference method , Commun. Comput. Phys. , 13 ( 4 ), 985 – 1012 . Google Scholar Crossref Search ADS d’Avila M. M.S. , Semblat J.-F. , Lenti L. , 2013 . Strong ground motion in the 2011 Tohoku Earthquake: a one-directional three-component modeling , Bull. seism. Soc. Am. , 103 ( 2B ), 1394 – 1410 . Google Scholar Crossref Search ADS Day S.M. , 1998 . Efficient simulation of constant Q using coarse-grained memory variables , Bull. seism. Soc. Am. , 88 , 1051 – 1062 . Day S.M. , Bradley C. , 2001 . Memory-efficient simulation of anelastic wave propagation , Bull. seism. Soc. Am. , 91 , 520 – 531 . Google Scholar Crossref Search ADS Dazel O. , Tournat V. , 2010 . Non-linear biot waves in porous media with application to unconsolidated granular media , J. acoust. Soc. Am. , 127 ( 2 ), 692 – 702 . Google Scholar Crossref Search ADS PubMed Delépine N. , Lenti L. , Bonnet G. , Semblat J.-F. , 2009 . Non-linear viscoelastic wave propagation: an extension of nearly constant attenuation models , J. Eng. Mech. , 135 ( 11 ), 1305 – 1314 . Google Scholar Crossref Search ADS Delsanto P.P. , 2006 . Universality of Nonclassical Non-linearity , Springer . Dhemaied A. , Rejiba F. , Camerlynck C. , Bodet L. , Guérin R. , 2011 . Seismic-wave propagation modeling in viscoelastic media using the auxiliary differential equation method , Bull. seism. Soc. Am. , 101 ( 6 ), 413 – 420 . Google Scholar Crossref Search ADS Donskoy D.M. , Khashanah K. , McKee T.G. Jr. , 1997 . Non-linear acoustic waves in porous media in the context of biot’s theory , J. acoust. Soc. Am. , 102 ( 5 ), 2521 – 2528 . Google Scholar Crossref Search ADS PubMed Favrie N. , Lombard B. , Payan C. , 2014 . Fast and slow dynamics in a non-linear elastic bar excited by longitudinal vibrations , Wave motion . Gomberg J. , Johnson P. , 2005 . Dynamic triggering of earthquakes , Nature , 437 , 830 . Google Scholar Crossref Search ADS PubMed Gomberg J. , Reasenberg P. , Bodin P. , Harris R. , 2001 . Earthquake triggering by transient seismic waves following the landers and hector mine earthquakes , Nature , 411 , 462 – 466 . Google Scholar Crossref Search ADS PubMed Gomberg J. , Bodin P. , Reasenberg P.A. , 2003 . Observing earthquakes triggered in the near field by dynamic deformations , Bull. seism. Soc. Am. , 93 ( 1 ), 118 . Google Scholar Crossref Search ADS Gomberg J. , Bodin P. , Larson K. , Dragert H. , 2004 . Earthquake nucleation by transient deformations caused by the m = 7.9 denali, alaska, earthquake , Nature , 427 , 621 – 624 . Google Scholar Crossref Search ADS PubMed Graves R.W. , 1996 . Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences , Bull. seism. Soc. Am. , 86 ( 4 ), 1091 – 1106 . Graves R.W. , 1999 . Three-dimensional computer simulations of realistic earthquake ground motions in regions of deep sedimentary basin , in The Effects of Surface Geology on Seismic Motion , Vol. 1 , pp. 103 – 120 ., eds Irikura K. , Kudo K. , Okada H. , Sasatani T. , Balkema . Greene J.H. , Taflove A. , 2006 . General vector auxiliary differential equation finite-difference time-domain method for non-linear optics , Opt. Express , 14 ( 18 ), 8305 – 8310 . Google Scholar Crossref Search ADS PubMed Hestholm S. , 1999 . Three-dimensional finite difference viscoelastic wave modelling including surface topography , Geophys. J. Int. , 139 ( 3 ), 852 – 878 . Google Scholar Crossref Search ADS Hill D.P. et al. , 1993 . Seismicity remotely triggered by the magnitude 7.3 Landers, California earthquake , Science , 260 , 1617 – 1623 . Google Scholar Crossref Search ADS PubMed Hokstad K. , 2004 . Non-linear and dispersive acoustic wave propagation , Geophysics , 69 ( 3 ), 840 – 848 . Google Scholar Crossref Search ADS Jiménez N. , Camarena F. , Redondo J. , Sánchez-Morcillo V. , Hou Y. , Konofagou E.E. , 2016 . Time-domain simulation of ultrasound propagation in a tissue-like medium based on the resolution of the non-linear acoustic constitutive relations , Acta Acust. United Acust. , 102 ( 5 ), 876 – 892 . Google Scholar Crossref Search ADS Johnson P.A. , Jia X. , 2005 . Non-linear dynamics, granular media and dynamic earthquake triggering , Nature , 437 ( 6 ), 871 – 874 . Google Scholar Crossref Search ADS PubMed Komatitsch D. , Martin R. , 2007 . An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation , Geophysics , 72 ( 5 ), SM155 – SM167 . Google Scholar Crossref Search ADS Kramer S.L. , 1996 . Geotechnical Earthquake Engineering , Prentice Hall . Kristek J. , Moczo P. , 2003 . Seismic-wave propagation in viscoelastic media with material discontinuities: a 3D fourth-order staggered-grid finite-difference modeling , Bull. seism. Soc. Am. , 93 ( 5 ), 2273 . Google Scholar Crossref Search ADS Legland J.B. , Tournat V. , Dazel O. , Novak A. , 2012 . Linear and non-linear biot waves in a noncohesive granular medium slab: transfer function, self-action, second harmonic generation , J. acoust. Soc. Am. , 131 ( 6 ), 4292 – 4303 . Google Scholar Crossref Search ADS PubMed Levander A.R. , 1988 . Fourth-order finite-difference P-SV seismograms , Geophysics , 53 ( 11 ), 1425 – 1436 . Google Scholar Crossref Search ADS Lombard B. , Piraux J. , 2011 . Numerical modeling of transient two-dimensional viscoelastic waves , J. Comput. Phys. , 230 , 6099 – 6114 . Google Scholar Crossref Search ADS Lombard B. , Piraux J. , Gélis C. , Virieux J. , 2008 . Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves , Geophys. J. Int. , 172 , 252 – 261 . Google Scholar Crossref Search ADS Martin R. , Komatitsch D. , 2009 . An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation , Geophys. J. Int. , 179 ( 1 ), 333 – 344 . Google Scholar Crossref Search ADS Martin R. , Komatitsch D. , Ezziani A. , 2008 . An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave equation in poroelastic media , Geophysics , 73 ( 4 ), T51 – T61 . Google Scholar Crossref Search ADS Martin R. , Komatitsch D. , Gedney S.D. , Bruthiaux E. , 2010 . A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML) , Comput. Model. Eng. Sci. , 56 ( 1 ), 17 – 42 . Moczo P. , Kristek J. , 2005 . On the rheological models used for time-domain methods of seismic wave propagation , Geophys. Res. Lett. , 32 , L01306 . Google Scholar Crossref Search ADS Moczo P. , Bystrický E. , Kristek J. , Carcione J.M. , Bouchon M. , 1997 . Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures , Bull. seism. Soc. Am. , 87 , 1305 – 1323 . Moczo P. , Kristek J. , Galis M. , Chaljub E. , Etienne V. , 2011 . 3-D finite-difference, finite-element, discontinuous-galerkin and spectral-element schemes analysed for their accuracy with respect to p-wave to s-wave speed ratio , Geophys. J. Int. , 187 ( 3 ), 1645 – 1667 . Google Scholar Crossref Search ADS Moczo P. , Kristek J. , Galis M. , 2014 . The Finite-Difference Modelling of Earthquake Motions: Waves and Ruptures , Cambridge University Press . Nazarov V. , 2001 . Acoustic non-linearity of cracks partially filled with liquid: cubic approximation , J. acoust. Soc. Am. , 109 ( 6 ), 2642 – 2648 . Google Scholar Crossref Search ADS PubMed Nazarov V. , Sutin A.E. , 1997 . Non-linear elastic constants of solids with cracks , J. acoust. Soc. Am. , 102 ( 6 ), 2249 – 3354 . Google Scholar Crossref Search ADS Olsen K.B. , Nighor R. , Konno T. , 2000 . 3D viscoelastic wave propagation in the upper Borrego Valley, California, constrained by borehole and surface data , Bull. seism. Soc. Am. , 90 ( 1 ), 134 – 150 . Google Scholar Crossref Search ADS Ostrovsky L. , 1991 . Wave processes in media with strong acoustic non-linearity , J. acoust. Soc. Am. , 90 ( 6 ), 3332 – 3337 . Google Scholar Crossref Search ADS Ostrovsky L. , Johnson P.A. , 2001 . Dynamic non-linear elasticity in geomaterials , Rivista del nuovo cimento , 24 ( 7 ). Régnier J. , Cadet H. , Bonilla L.F. , Bertrand E. , Semblat J. , 2013 . Assessing non-linear behavior of soils in seismic site response: statistical analysis on KiKnet strong motion data , Bull. seism. Soc. Am. , 103 ( 3 ), 1750 . Google Scholar Crossref Search ADS Régnier J. , Bonilla L.F. , Bertrand E. , Semblat J. , 2014 . Influence of the VS profiles beyond 30 m depth on linear site effects: assessment from the KiKnet data , Bull. seism. Soc. Am. , 104 ( 5 ), 2337 . Google Scholar Crossref Search ADS Régnier J. et al. , 2016a . International benchmark on numerical simulations for 1D, non-linear site response (PRENOLIN): verification phase based on canonical cases , Bull. seism. Soc. Am. , 106 ( 5 ), 2112 . Google Scholar Crossref Search ADS Régnier J. , Cadet H. , Bard P. , 2016b . Empirical quantification of the impact of non-linear soil behavior on site response , Bull. seism. Soc. Am. , 106 ( 4 ), 1710 . Google Scholar Crossref Search ADS Régnier J. et al. , 2018 . PRENOLIN: international benchmark on 1D non-linear site response analysis-validation phase exercise , Bull. seism. Soc. Am. , 108 ( 2 ), 876 . Rejiba F. , Camerlynck C. , Mechler P. , 2003 . FDTD-SUPML-ADE simulation for Ground-Penetrating Radar modelling , Radio Sci. , 38 ( 1 ), 1005 . Google Scholar Crossref Search ADS Robertsson J.O. , Blanch J.O. , Symes W.W. , 1994 . Viscoelastic finite-difference modeling , Geophysics , 59 ( 9 ), 1444 – 1456 . Google Scholar Crossref Search ADS Saenger E.H. , Bohlen T. , 2004 . Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid , Geophysics , 69 ( 2 ), 583 – 591 . Google Scholar Crossref Search ADS Sandikkaya M.A. , Akkar S. , Bard P. , 2013 . A non-linear site amplification model for the next Pan European ground motion prediction equations , Bull. seism. Soc. Am. , 103 ( 1 ), 19 . Google Scholar Crossref Search ADS Semblat J.F. , Pecker A. , 2009 . Waves and vibrations in soils: earthquakes, Traffic, Shocks, construction works , in Waves and Vibrations in Soils: Earthquakes, Traffic, Shocks, Construction Works , 500 pp, IUSS Press . Sun Y. , Zhang W. , Chen X. , 2018 . 3d seismic wavefield modeling in generally anisotropic media with a topographic free surface by the curvilinear grid finite difference method , Bull. seism. Soc. Am. , 108 ( 3A ), 1287 . Google Scholar Crossref Search ADS Taflove A. , Hagness S.C. , 2005 . Computational Electrodynamics: The Finite-Difference Time-Domain Method , Artech house . Taflove A. , Oskooi A. , Johnson S.G. , 2013 . Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology , Artech House . Tournat V. , Gusev V. , 2010 . Acoustics of unconsolidated ”model” granular media: an overview of recent results and several open problems , Acta Acust. United Acust. , 96 , 208 – 224 . Google Scholar Crossref Search ADS Tournat V. , Castagnède B. , Gusev V. , Béquin P. , 2003 . Self-demodulation acoustic signatures for non-linear propagation in glass beads , C.R. Mecanique , 331 , 119 – 125 . Google Scholar Crossref Search ADS Tournat V. , Gusev V. , Zaitsev V. , Castagnède B. , 2004 . Acoustic second-harmonic generation with shear to longitudinal mode conversion in granular media , Europhys. Lett. , 66 ( 6 ), 798 – 804 . Google Scholar Crossref Search ADS Wallen S. , Boechler N. , 2017 . Shear to longitudinal mode conversion via second harmonic generation in a two-dimensional microscale granular crystal , Wave Motion , 68 , 22 – 30 . Google Scholar Crossref Search ADS Wang Y. , Xu J. , Schuster G. , 2001 . Viscoelastic wave simulation in basins by a variable-grid finite-difference method , Bull. seism. Soc. Am. , 91 ( 6 ), 1741 – 1749 . Google Scholar Crossref Search ADS Xu T. , McMechan G. , 1998 . Efficient 3-D viscoelastic modeling with application to near-surface land seismic data , Geophysics , 63 ( 2 ), 601 – 612 . Google Scholar Crossref Search ADS Yu-Lin F. , Xiao-Zhou L. , Jie-Hui L. , Li M. , 2009 . The non-linear propagation of acoustic waves in a viscoelastic medium containing cylindrical micropores , Chin. Phys. B , 18 ( 9 ), 3909 – 3917 . Google Scholar Crossref Search ADS Zhang W. , Chen X. , 2006 . Traction image method for irregular free surface boundaries in finite difference seismic wave simulation , Geophys. J. Int. , 167 ( 1 ), 337 – 353 . Google Scholar Crossref Search ADS Zhang W. , Shen Y. , 2010 . Unsplit complex frequency-shifted pml implementation using auxiliary differential equations for seismic wave modeling , Geophysics , 75 ( 4 ), T141 – T154 . Google Scholar Crossref Search ADS Zhang W. , Zhang Z. , Chen X. , 2012 . 3-d elastic wave numerical modeling in the presence of surface topography by a collocated-grid finite difference method on curvilinear grid , Geophys. J. Int. , 190 , 358 – 378 . Google Scholar Crossref Search ADS APPENDIX: VERIFICATION AND ACCURACY TESTS OF THE NUMERICAL METHOD In the linear viscoelastic case, it is important to check the accuracy of the numerical method used. For this purpose, we compare, in the frequency domain and at as 0.5 m distance from the source, the analytical and the numerically computed compressional and shear phase velocities defined as the norms of MP/ρ and MS/ρ, where MP and MS are the moduli given in eq. (7) for a 1500 Hz dominant frequency Ricker or Gaussian-modulated sine pulse wavelet source. In Figs A1b–e), we can observe that P and S phase velocities are very well represented numerically, particularly in a frequency range between around 200 Hz and Πf0 ≃ 4712 Hz. Maximum errors of less than 0.3  and 0.2 per cent are obtained on P and S phase velocities for the Ricker and Gauss-modulated sine source cases, respectively, which is very good. In the non-linear case, we consider the strain ε numerically computed at each time step and we use it to compute the stress directly in the frequency domain after Fourier transform of the desired stress–strain defined by eq. (9). We compare this 'theoretical’ solution to the FD numerical stress solution in the frequency domain. We can see in Fig. A1f that the numerical solution of the σzz component is very accurate. Figure A1. View largeDownload slide Similar numerical and analytical velocity component vz computed in the time domain for the Ricker source (a). Numerical and theoretical P- and S-wave phase velocities for the Ricker (b and c) and the Gaussian-modulated sine (d and e) pulse source cases are compared for the linear viscoelastic case at a 0.5 m distance from the source along the z axis. In (f), for the non-linear viscoelastic case and the Gaussian-modulated sine pulse source of 2.25 maximum amplitude, the numerical and analytical values of the σzz stress component are compared in the frequency domain with very good agreement. Figure A1. View largeDownload slide Similar numerical and analytical velocity component vz computed in the time domain for the Ricker source (a). Numerical and theoretical P- and S-wave phase velocities for the Ricker (b and c) and the Gaussian-modulated sine (d and e) pulse source cases are compared for the linear viscoelastic case at a 0.5 m distance from the source along the z axis. In (f), for the non-linear viscoelastic case and the Gaussian-modulated sine pulse source of 2.25 maximum amplitude, the numerical and analytical values of the σzz stress component are compared in the frequency domain with very good agreement. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Seismic wave propagation in nonlinear viscoelastic media using the auxiliary differential equation method
JF - Geophysical Journal International
DO - 10.1093/gji/ggy441
DA - 2019-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/seismic-wave-propagation-in-nonlinear-viscoelastic-media-using-the-70jwxpYOSm
SP - 453
VL - 216
IS - 1
DP - DeepDyve
ER -