Abderrahmin Ben Jazia
email: abderrahim.ben-jazia@centrale-marseille.fr

Cédric Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

C Bellis
email: bellis@lma.cnrs-mrs.fr

Wave propagation in a fractional viscoelastic Andrade medium: diffusive approximation and numerical modeling

Keywords: Viscoelasticity, Andrade model, Fractional derivatives, Transient wave propagation, Finite differences

Wave propagation in a fractional viscoelastic Andrade medium: diffusive approximation and numerical modeling

Introduction

There is a long history of studies discussing or providing experimental evidences of frequency-dependent viscoelastic attenuations, as observed in e.g. metals [1], acoustic media [2, 3] and in the Earth [START_REF] Flanagan | Attenuation of broadband P and S waves in Tonga: Observations of frequency dependent Q[END_REF][START_REF] Lekić | Measurement and implications of frequency dependence of attenuation[END_REF]. Such a behavior is classically modeled using a fractional derivative operator [START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF][START_REF] Carcione | Wave Fields in Real Media: Wave propagation in Anisotropic, Anelastic and Porous Media[END_REF], a mathematical tool generalizing to real parameters the standard derivatives of integer orders [START_REF] Podlubny | Fractional Differential Equations[END_REF]. While fractional calculus is now a mature theory in the field of viscoelasticity [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF], some issues remain commonly encountered. They mostly revolve around the two questions of: (i) Incorporating fractional dissipation into viscoelastic models that both fit experimental data and have a theoretical validity regarding, e.g., causality properties [START_REF] Holm | A causal and fractional all-frequency wave equation for lossy media[END_REF][START_REF] Waters | Causality-imposed (Kramers-Kronig) relationships between attenuation and dispersion[END_REF] or the Kramers-Kronig relations [START_REF] Waters | On the applicability of Kramers-Krönig relations for ultrasonic attenuation obeying a frequency power law[END_REF].

(ii) Implementing numerically these fractional models to perform wave propagation simulations. This latter problem is commonly tackled using standard approaches [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] for modeling constant-law of attenuation over a frequency-band of interest, i.e. with the fractional viscoelastic model being approximated by multiple relaxation mechanisms [START_REF] Näsholm | Linking multiple relaxation, power-law attenuation, and fractional wave equations[END_REF].

Bearing in mind the issue (i) discussed above, it is chosen to anchor the present study to a specific, yet prototypical, physically-based viscoelastic model, namely the Andrade model. Initially introduced in [1] to fit experimental flow laws in metals, it has been further investigated in [START_REF] Andrade | On the validity of t 1/3 law of flow of metals[END_REF]. It is now used as a reference in a number of studies [START_REF] Gribb | Low-frequency shear attenuation in polycrystalline olivine: Grain boundary diffusion and the physical significance of the Andrade model for viscoelastic rheology[END_REF][START_REF] Jackson | Grainsize-sensitive viscoelastic relaxation in olivine: Towards a robust laboratory-based model for seismological application[END_REF][START_REF] Sundberg | A composite viscoelastic model for incorporating grain boundary sliding and transient diffusion creep: Correlating creep and attenuation responses for materials with a fine grain size[END_REF][START_REF] Bellis | Sensitivity of seismic measurements to frequency-dependent attenuation and upper mantle structure: an initial approach[END_REF] for the description of observed frequency-dependent damping behaviors in the field of geo-physics and experimental rock mechanics. Moreover, the Andrade model creep function, as written, can notably be decomposed as the sum of a fractional power-law and a standard Maxwell creep function, therefore corresponding rheologically to a spring-pot element arranged in series with a spring-dashpot Maxwell model. Therefore, while being physically motivated and rooted in experiments, this model gives leeway to cover the spectrum from a conventional rheological mechanism to a more complex fractional model, and this with only a few parameters.

This study focuses on the issue (ii), namely the numerical modeling of wave propagation within an Andrade medium that exhibits fractional attenuation. The objective is to develop an efficient approximation strategy of the fractional term featured in this viscoelastic model in view of the investigation and simulation of its transient dynamical behavior.

A model-based approach is explored in the sense that one aims at a direct approximation of the original constitutive equation. Therefore, the latter is not intended to be superseded by another viscoelastic model that would be designed to fit only a given observable. For example, the usual approach that employs a multi-Zener model typically approximates the quality factor only.

The article aim and contribution are twofold:

(i) Deploy an approximation of the fractional derivative featured in the constitutive equation considered. A direct discretization of this term, that is associated with a non-local time-domain convolution product [START_REF] Podlubny | Fractional Differential Equations[END_REF] requires the storage of the entire variables history, which is out of reach for realistic simulations. The Grünwald-Letnikov approximation of fractional derivatives constitutes a tractable approach, commonly used in viscoelasticity [START_REF] Galucio | Finite element formulation of viscoelastic sandwich beams using fractional derivative operators[END_REF]. Its main drawback concerns the stability analysis to be performed for the numerical scheme so-obtained. Indeed, Von-Neumann stability of multistep schemes requires to bound the characteristic roots of the amplification matrix, which may be a difficult task.

We do not follow this approach here. Alternatively, a so-called diffusive representation is preferred [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF], as it allows to recast the equations considered into a local-in-time system while introducing only a limited number of additional memory variables in its discretized form [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF]. Following later improvements of the method in [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF][START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF][START_REF] Haddar | Efficient solution of a wave equation with fractional-order dissipative terms[END_REF][START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF], an efficient quadrature scheme is investigated in order to obtain a satisfactory fit of the reference model compliance.

(ii) Implement the resulting approximated model into a wave propagation scheme.

While the available literature on the numerical simulation of transient wave propagation within fractionally-damped media is relatively scarce, see e.g. [START_REF] Wismer | Finite element analysis of broadband acoustic pulses through inhomogenous media with power law attenuation[END_REF][START_REF] Caputo | Wave simulation in biologic media based on the Kelvin-Voigt fractional-derivative stress-strain relation[END_REF], the aim is here to demonstrate the efficiency of the proposed approach. For the sake of simplicity, the viscoelastic medium considered is assumed to be unidimensional and homogeneous. After discretization of the dynamical system at hand, a Strang splitting approach [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] is adopted, both to reach an optimal stability condition and to enable the use of an efficient time-marching scheme coupled with an exact integration step. Moreover, deriving a semianalytical solution for the configuration considered, as a baseline, a set of numerical results is presented to assess the quality of the numerical scheme developed. The overall features and performances of the diffusive representation are finally discussed to compare the Andrade model with its diffusive approximated counterpart.

Notably, this study demonstrates that the behavior of fractional viscoelastic models such as the Andrade model can be correctly described using a diffusive approximation.

The resulting model is shown to be well characterized mathematically while being easily tractable numerically in view of performing simulations in the time domain. This article is organized as follows. The Andrade model is presented and discussed in Section 2. Considering the featured fractional derivative, a corresponding diffusive approximated (DA) version of the former is subsequently formulated and referred to as the Andrade-DA model. The evolution problem is investigated in Section 3, with the derivation and analysis of the first-order hyperbolic system associated with the Andrade-DA model. Section 4.1 is concerned with the definition and computation of an efficient quadrature scheme for the diffusive approximation, while the implementation of the fully discretized system is described in Section 4.2. Corresponding numerical results are presented and discussed in Section 5.

Fractional viscoelastic model

Preliminaries

The causal constitutive law describing the behavior of a 1D linear viscoelastic medium can be expressed in terms of the time-domain convolution

ε(t) = t 0 χ(t -τ ) ∂σ ∂τ (τ ) dτ, (1)
with creep function χ, stress field σ and strain field ε = ∂u/∂x associated with unidimensional displacement u, time t and space coordinate x.

Next, for parameters satisfying 0 < β < 1, the so-called Caputo-type fractional derivative [START_REF] Carcione | Wave Fields in Real Media: Wave propagation in Anisotropic, Anelastic and Porous Media[END_REF][START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF][START_REF] Podlubny | Fractional Differential Equations[END_REF] of a causal function g(t) is defined as

d β g dt β (t) = 1 Γ(1 -β) t 0 (t -τ ) -β dg dτ (τ ) dτ, ( 2
)
where Γ is the Gamma function. Defining the direct and inverse Fourier transforms in time of a function g(t) as

ĝ(ω) = +∞ -∞ g(t)e -iωt dt, g(t) = 1 2π +∞ -∞ ĝ(ω)e iωt dω,
where ω is the angular frequency and i = √ -1, then the frequency-domain counterpart of equation (2) reads

d β g dt β (ω) = (iω) β ĝ(ω), (3)
so that definition (2) is a straightforward generalization of the derivative of integer order.

Andrade model

The Andrade model [1] is characterized by the creep function given by

χ(t) = J u + t η + A t α H(t), 0 < α < 1, (4)
with Heaviside step function H(t), unrelaxed compliance J u , viscosity η and two positive physical parameters A and α. Usual fits with experimental data correspond to 1 3 ≤ α ≤ 1 2 [START_REF] Andrade | On the validity of t 1/3 law of flow of metals[END_REF][START_REF] Gribb | Low-frequency shear attenuation in polycrystalline olivine: Grain boundary diffusion and the physical significance of the Andrade model for viscoelastic rheology[END_REF]. The composite law (4) can be additively decomposed into a standard Maxwell rheological mechanism with creep function t → J u + t/η and a relaxation time τ Mx = η J u , together with a power law dependence in time t → A t α which constitutes its main feature.

Examples behaviors of the creep function (4) are illustrated in Figure 1a.

The Fourier transforms ε and σ of the strain and stress are linked by ε = N σ with the complex compliance N being defined as N(ω) = iω χ(ω). The latter can be deduced from the Fourier transform χ of the creep function (4) as

N(ω) = J u + (iη ω) -1 + A Γ(1 + α) (i ω) -α . (5)
Straightforward manipulations on (3), ( 4) and ( 5) lead to the following constitutive equation in differential form for the Andrade model

∂ε ∂t = J u ∂σ ∂t + 1 η σ + A Γ(1 + α) ∂ 1-α ∂t 1-α σ. (6)

Dispersion relations

The complex wave number k(ω) satisfies

k(ω) = √ ρ ω [N(ω)] 1/2 := ω c(ω) -iζ(ω), (7)
where the phase velocity c and the attenuation ζ are given by

c(ω) = 2 ρ(|N| + Re[N]) , ζ(ω) = ω ρ(|N| -Re[N]) 2 . ( 8
)
Owing to equations ( 5) and ( 8), the following limits hold:

lim ω→0 c(ω) = 0, lim ω→+∞ c(ω) = 1 √ ρJ u := c ∞ , lim ω→0 ζ(ω) = 0, lim ω→+∞ ζ(ω) = +∞. (9)
Moreover, when A > 0, the creep function ( 4) is an increasing and concave function.

As a consequence, owing to the theoretical developments in [START_REF] Hanyga | On wave propagation in viscoelastic media with concave creep compliance[END_REF] and [START_REF] Hanyga | Wave propagation in linear viscoelastic media with completely monotonic relaxation moduli[END_REF], the attenuation The other physical parameters are: ρ = 1200 kg/m 3 , c ∞ = 2800 m/s and η = 10 9 Pa.s.

(a) Creep function χ(t) 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3 f (Hz) Q A=0 A=1E-10 A=2E-10 1/3 2/3 (b) Quality factor Q(f ) 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3
The horizontal solid line in panel (c) denotes the high-frequency limit c ∞ .

ζ(ω) for the Andrade model turns out to be sublinear in the high-frequency range, i.e.

ζ(ω) = ω→+∞ o(ω). (10)
This key property confirms the relevance of the choice of the Andrade model as a prototypical example of fractional viscoelastic media.

The quality factor Q is defined as the ratio

Q(ω) = - Re[k 2 ] Im[k 2 ] = - Re[N] Im[N] . (11)
According to [START_REF] Lekić | Measurement and implications of frequency dependence of attenuation[END_REF] and in the low and high-frequency regimes, the frequency-dependent behavior follows

Q(ω) ∼ ω→0 Q 0 ω 1-α with Q 0 = η A Γ(1 + α) cos απ 2 , Q(ω) ∼ ω→+∞ Q ∞ ω α with Q ∞ = J u A Γ(1 + α) sin απ 2 -1 . ( 12
)
Sample behaviors of the Andrade model for α = 1/3 and a varying parameter A are sketched in Figure 1. Notably, the case A = 0 corresponds to the standard Maxwell model.

The corresponding attenuation curve shows that, within the frequency range considered, the associated high-frequency regime ζ(ω) ∼ ω→+∞ 1 2η ρ Ju is rapidly attained. Alternatively, when A = 0, one observes in Fig. 1b the slopes 2/3 and 1/3 of the quality factor in log-log scale at low and high frequencies respectively, as expected from [START_REF] Waters | On the applicability of Kramers-Krönig relations for ultrasonic attenuation obeying a frequency power law[END_REF]. The attenuation ζ is represented as a function of the frequency f and displayed in linear scale in Fig. 1d to emphasize the sublinear high-frequency behavior [START_REF] Holm | A causal and fractional all-frequency wave equation for lossy media[END_REF].

Diffusive approximation: Andrade-DA model

When implementing [START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF], the difficulty revolves around the computation of the convolution product in (2) associated with the fractional derivative of order 1 -α, which is numerically memory-consuming. The alternative approach adopted in this study is based on a diffusive representation, and its approximation, of fractional derivatives. Following [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF], then for 0 < α < 1 equation (2) can be recast as

∂ 1-α ∂t 1-α σ = +∞ 0 φ(x, t, θ) dθ, (13)
where the function φ is defined owing to a change of variables as

φ(x, t, θ) = 2 sin(πα) π θ 1-2α t 0 ∂σ ∂τ (x, τ ) e -(t-τ ) θ 2 dτ. ( 14
)
As φ is expressed in terms of an integral operator with decaying exponential kernel, it is referred to as a diffusive variable. From equation ( 14), it can be shown to satisfy the following first-order differential equation for θ > 0:

     ∂φ ∂t = -θ 2 φ + 2 sin(πα) π θ 1-2α ∂σ ∂t , φ(x, 0, θ) = 0. (15)
The diffusive representation [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF][START_REF] Näsholm | Linking multiple relaxation, power-law attenuation, and fractional wave equations[END_REF] amounts to supersede the non-local term in ( 6) by an integral over θ of the function φ(x, t, θ) which obeys the local first-order ordinary differential equation [START_REF] Andrade | On the validity of t 1/3 law of flow of metals[END_REF]. The integral featured in ( 13) is in turn well-suited to be approximated using a quadrature scheme, so that

∂ 1-α ∂t 1-α σ ≃ L ℓ=1 µ ℓ φ(x, t, θ ℓ ) ≡ L ℓ=1 µ ℓ φ ℓ (x, t), (16)
given a number L of quadrature nodes θ ℓ with associated weights µ ℓ . These parameters with unit of s -1/2 and s 1/2 respectively, and whose computations will be returned to in Section 4.1, will be seen to be decided from the fit of the Andrade model complex compliance.

The frequency-domain versions of equations ( 6), ( 15) and ( 16) lead to the approximated complex compliance Ñ, such that ε = Ñ σ and characterizing the model hereafter referred to as the Andrade-DA model, as

Ñ (ω) = J u + (iη ω) -1 + AΓ(1 + α) 2 sin(πα) π L ℓ=1 µ ℓ θ 1-2α ℓ θ 2 ℓ + iω . (17)
A comparison between (5) and its diffusive approximated counterpart [START_REF] Jackson | Grainsize-sensitive viscoelastic relaxation in olivine: Towards a robust laboratory-based model for seismological application[END_REF] shows that the corresponding complex compliances N and Ñ differ only in their third terms. Therefore, based on equation ( 7), the associated dispersion relations read

k 2 = ω c ∞ 2 1 + AΓ(1 + α) J u κ mod (ω) - iρ ω η (18)
with the function κ mod being defined for the two models considered by

κ mod (ω) =          κ(ω) = (iω) -α Andrade, κ(ω) = 2 sin(πα) π L ℓ=1 µ ℓ θ 1-2α ℓ θ 2 ℓ + iω Andrade-DA. ( 19
)
Finally, the diffusive approximated counterparts of the phase velocity and the attenuation function in ( 8) can be immediately deduced using [START_REF] Sundberg | A composite viscoelastic model for incorporating grain boundary sliding and transient diffusion creep: Correlating creep and attenuation responses for materials with a fine grain size[END_REF][START_REF] Bellis | Sensitivity of seismic measurements to frequency-dependent attenuation and upper mantle structure: an initial approach[END_REF]. In particular, the low-frequency and high-frequency limits of the phase velocity c are equal to those in [START_REF] Mainardi | Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models[END_REF]. Moreover, using tables of standard Fourier transforms, the corresponding time-domain creep function χ, defined by Ñ = iω χ, is obtained as

χ(t) = J u + t η + AΓ(1 + α) 2 sin(πα) π L ℓ=1 µ ℓ θ -1-2α ℓ 1 -e -θ 2 ℓ t H(t). (20)

Evolution equations

With the complex compliance ( 17) of the Andrade-DA model at hand, which constitutes the approximated version of the diffusive representation of the Andrade model ( 5), the present section is concerned with the description and analysis of its dynamical behavior.

First-order system

Let define the parameters

γ ℓ,α = 2 sin(πα) π J u θ 1-2α ℓ , Υ ℓ,α = AΓ(1 + α) γ ℓ,α for ℓ = 1, . . . , L. (21)
Combining the conservation of momentum in terms of velocity field v = ∂u/∂t and equations ( 6), ( 15) and ( 16) yields

                     ∂v ∂t - 1 ρ ∂σ ∂x = F v , (22a)
∂σ ∂t - 1 J u ∂v ∂x = - 1 J u η σ - AΓ(1 + α) J u L j=1 µ j φ j + F σ , (22b)
∂φ ℓ ∂t -γ ℓ,α ∂v ∂x = -θ 2 ℓ φ ℓ - γ ℓ,α η σ -Υ ℓ,α L j=1 µ j φ j + J u γ ℓ,α F σ , (22c)
for ℓ = 1, . . . , L and where F v and F σ are introduced to model external sources. Equations ( 22) are completed by initial conditions

v(x, 0) = 0, σ(x, 0) = 0, φ ℓ (x, 0) = 0 for ℓ = 1, • • • , L.
Gathering unknown and sources terms, let the vectors U and F be defined as

U = v, σ, φ 1 , • • • , φ L T , F = F v , F σ , J u γ 1,α F σ , • • • , J u γ L,α F σ T . ( 23
)
Then the system ( 22) can be written in the matrix-form

∂U ∂t + A ∂U ∂x = SU + F , ( 24
)
where A is given by

A =            0 -ρ -1 0 • • • 0 -J u -1 0 0 • • • 0 -γ 1,α 0 0 • • • 0 . . . . . . . . . . . . . . . -γ L,α 0 0 • • • 0            , (25)
and S reads

S =             0 0 0 • • • 0 0 -(J u η) -1 -AΓ(1 + α)J u -1 µ 1 • • • -AΓ(1 + α)J u -1 µ L 0 -γ 1,α η -1 -θ 2 1 -Υ 1,α µ 1 • • • -Υ 1,α µ L . . . . . . . . . . . . . . . 0 -γ L,α η -1 -Υ L,α µ 1 • • • -θ 2 L -Υ L,α µ L             . ( 26
)
Note that this differential system remains valid in the case of a non-homogeneous viscoelastic medium.

Energy decay

Studying the energy associated with the system ( 22) is required to characterize the stability of the Andrade-DA model and to provide constraints on the diffusive approximation calculation. For an infinite 1D domain, the stored kinetic and elastic energies are defined as

E v (t) = 1 2 +∞ -∞ ρv 2 dx and E σ (t) = 1 2 +∞ -∞ J u σ 2 dx, (27)
together with a coupled term associated with the diffusive approximation

E d (t) = 1 2 +∞ -∞ L ℓ=1 µ ℓ Υ ℓ,α θ 2 ℓ J u σ - φ ℓ √ J u γ ℓ,α 2 dx. ( 28
)
Then, in the absence of any source term, one has the following property

Proposition 1. If µ ℓ > 0 for all ℓ = 1, . . . , L, then the function E(t) = E v (t)+E σ (t)+E d (t)
is a positive definite quadratic form and dE dt < 0 for all time t > 0.

Proof. In the absence of any source term, then multiplying the momentum equation (22a) by the velocity field v and integrating spatially by parts yields

+∞ -∞ ρv ∂v ∂t + σ ∂v ∂x dx = 0,
assuming that the elastic fields vanish at infinity. Likewise, from equation (22b) and multiplying by σ, one obtains

1 2 d dt +∞ -∞ ρv 2 + J u σ 2 dx + +∞ -∞ σ 2 η + AΓ(1 + α) L ℓ=1 µ ℓ φ ℓ σ dx = 0. ( 29
)
Now, using twice differential equation ( 15), one has for ℓ = 1, . . . , L

σ ∂φ ℓ ∂t + θ 2 ℓ φ ℓ σ -J u γ ℓ,α σ ∂σ ∂t = 0 and φ ℓ J u γ ℓ,α ∂φ ℓ ∂t + θ 2 ℓ φ 2 ℓ J u γ ℓ,α -φ ℓ ∂σ ∂t = 0,
which after subtraction and manipulation entails

φ ℓ σ = φ 2 ℓ J u γ ℓ,α + γ ℓ,α 2θ 2 ℓ d dt J u σ - φ ℓ √ J u γ ℓ,α 2 . ( 30
)
Finally, substituting [START_REF] Hanyga | On wave propagation in viscoelastic media with concave creep compliance[END_REF] in [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] leads to the relation

1 2 d dt +∞ -∞ ρv 2 + J u σ 2 + L ℓ=1 µ ℓ Υ ℓ,α θ 2 ℓ J u σ - φ ℓ √ J u γ ℓ,α 2 dx = - +∞ -∞ σ 2 η + L ℓ=1 µ ℓ Υ ℓ,α φ ℓ √ J u γ ℓ,α
which concludes the proof, owing to the definition of the total energy function E from [START_REF] Wismer | Finite element analysis of broadband acoustic pulses through inhomogenous media with power law attenuation[END_REF] and [START_REF] Caputo | Wave simulation in biologic media based on the Kelvin-Voigt fractional-derivative stress-strain relation[END_REF].

In summary, positivity of the quadrature nodes and weights in ( 16) is crucial to ensure the well-posedness of the system [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF]. This issue will be further discussed in Section 4.1.

Properties of matrices

Some properties of the matrices A (25) and S (26) are discussed to characterize the first-order system (24) of partial differential equations.

Proposition 2. The eigenvalues of the matrix A are

sp(A) = {0, ±c ∞ } , with 0 being of multiplicity L.
As A is diagonalizable with real eigenvalues, then equation ( 24) is a hyperbolic system of partial differential equations, with solutions of finite-velocity. It is emphasized that the eigenvalue c ∞ = 1/ √ ρJ u does not depend on the set of quadrature coefficients {(µ ℓ , θ ℓ )}, so that the phase velocity upper bounds for the Andrade and Andrade-DA models are equal.

Proposition 3. Assuming θ ℓ > 0 and µ ℓ > 0 for ℓ = 1, . . . , L then sp(S) ∋ 0 with multiplicity 1. Moreover the L + 1 non-zero eigenvalues λ ℓ of S are real and, ordering the

nodes as 0 < θ 1 < • • • < θ L , satisfy λ L+1 < -θ 2 L < • • • < -θ 2 ℓ < λ ℓ < -θ 2 ℓ-1 < • • • < λ 1 < 0.
Proof. Let P S (λ) denote the characteristic polynomial of the matrix S, i.e. P S (λ) = det(S -λI L+2 ) with I L+2 the (L + 2)-identity matrix. The line i and the column j of the determinant are denoted by L i and C j , respectively. The following algebraic manipulations are performed successively:

(i) L j ← L j -γ α θ 1-2α j L 1 with j = 2, . . . , L + 1 (ii) C 1 ← C 1 L ℓ=1 (-θ 2 ℓ -λ) (iii) C 1 ← C 1 -γ α θ 1-2α ℓ λ C ℓ L i=1 i =ℓ (-θ 2 i -λ) for ℓ = 2, . . . , L + 1.
From [START_REF] Deü | Simulation of fractionally damped mechanical systems by means of a Newmark-diffusive scheme[END_REF] and definition [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF] of parameters γ ℓ,α and Υ ℓ,α , one deduces

P S (λ) = λ (J u η) -1 + λ L ℓ=1 (-θ 2 ℓ -λ) + λ L ℓ=1 µ ℓ Υ ℓ,α L j=1 j =ℓ (-θ 2 j -λ) := λ Q S (λ).
From the above equation, one has P S (0) = 0 while Q S (0) = 0, therefore 0 is an eigenvalue of the matrix S with multiplicity 1. In the limit λ → 0, then asymptotically

P S (λ) ∼ λ→0 (-1) L (J u η) -1 λ L ℓ=1 θ 2 ℓ , so that sgn(P S (0 -)) = (-1) L+1 . (31)
Moreover, using [START_REF] Matignon | Stability properties for generalized fractional differential systems[END_REF] and the assumptions considered, then at the quadrature nodes one has for all k ∈ {1, . . . , L}

P S (-θ 2 k ) = - 2 sin(πα)A Γ(1 + α) πJ u µ k θ 5-2α k L j=1 j =k (θ 2 k -θ 2 j ) ⇒ sgn(P S (-θ 2 k )) = (-1) L-k+1 .
Finally, the following limit holds

P S (λ) ∼ λ→-∞ (-1) L λ L+2 ⇒ sgn(P S (-∞)) = 1. (32)
We introduce the following intervals

I L+1 = -∞, -θ 2 L , I ℓ+1 = -θ 2 ℓ+1 , -θ 2 ℓ for ℓ = 1, . . . , L-1 and I 1 = -θ 2 1 , 0 . (33)
Given that λ → P S (λ) is continuous, then equations [START_REF] Hanyga | Wave propagation in linear viscoelastic media with completely monotonic relaxation moduli[END_REF][START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF] show that the polynomial P S changes sign in each of the intervals I ℓ of (33). Consequently, there exist λ ℓ ∈ I ℓ with ℓ = 1, . . . , L + 1 such that P S (λ ℓ ) = 0 and which coincide with the eigenvalues, with multiplicity 1, of the matrix S of size L + 2.

Proposition 3 states that, under suitable conditions on the quadrature coefficients, the matrix S in (26) has eigenvalues with negative or zero real parts. This property is crucial regarding the numerical modeling developed in the forthcoming Section 4.2. As for the energy analysis given in Proposition 1, positivity of quadrature nodes and weights is again the fundamental hypothesis. Lastly, it is possible to use the above proposition to characterize the spectral radius of the matrix S.

Proposition 4. The spectral radius of the matrix S ( 26) is such that

max θ 2 L , (J u η) -1 + L ℓ=1 µ ℓ Υ ℓ,α ≤ ̺(S) ≤ θ 2 L + (J u η) -1 + L ℓ=1 µ ℓ Υ ℓ,α .
Proof. By definition, one has

tr(S) = -(J u η) -1 + L ℓ=1 (θ 2 ℓ + µ ℓ Υ ℓ,α ) ≡ L+1 ℓ=1 λ ℓ . (34)
According to the proof of Property 3, the eigenvalues λ ℓ satisfy

- L ℓ=1 θ 2 ℓ ≤ L ℓ=1 λ ℓ ≤ - L-1 ℓ=1 θ 2 ℓ .
Substitution in [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF] and providing that ̺(S) = |λ L+1 | allows to conclude the proof.

Semi-analytical solutions

Let us consider a homogeneous medium described either by the Andrade model, i.e.

(x, ω) = i F (ω) 2πc 2 ∞ J u +∞ -∞ k k 2 -k 2 0 e ik(x-xs) dk,
with k 0 being defined for the two models considered according to [START_REF] Sundberg | A composite viscoelastic model for incorporating grain boundary sliding and transient diffusion creep: Correlating creep and attenuation responses for materials with a fine grain size[END_REF][START_REF] Bellis | Sensitivity of seismic measurements to frequency-dependent attenuation and upper mantle structure: an initial approach[END_REF] as

k 0 = ω c ∞ 2 1 + AΓ(1 + α) J u κ mod - iρ ω η 1/2 .
Note that choosing κ mod = κ or κ mod = κ yields the solution associated with the Andrade or with the Andrade-DA model respectively. The poles ±k 0 of the integrand are simple and satisfy Im[k 0 ] < 0. Using the residue theorem, one obtains in the time-domain the stress field solution

σ(x, t) = - sgn(x -x s ) 2πc 2 ∞ J u +∞ 0 Re F (ω)e i(ωt-k 0 |x-xs|) dω. ( 35
)
Similarly, the velocity field satisfies

v(x, t) = 1 2π +∞ 0 Re k 0 ω F (ω)e i(ωt-k 0 |x-xs|) dω. (36)
Finally, for the Andrade-DA model, the associated memory variables φ ℓ are expressed as

φ ℓ (x, t) = - sgn(x -x s )γ ℓ,α 2πc 2 ∞ +∞ 0 Re iω θ 2 ℓ + iω F (ω)e i(ωt-k 0 |x-xs|) dω, ℓ = 1, . . . , L. (37)
In the numerical results presented Section 5, the frequency-domain integrals featured in solutions [START_REF] Kappel | An implementation of Shor's r-algorithm[END_REF], [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] and (37) are computed using a standard quadrature rule over the frequency-band considered.

Numerical methods

Quadrature methods

Two different approaches can be employed to determine the set {(µ ℓ , θ ℓ )} of 2L coefficients of the diffusive approximation [START_REF] Gribb | Low-frequency shear attenuation in polycrystalline olivine: Grain boundary diffusion and the physical significance of the Andrade model for viscoelastic rheology[END_REF]. While the most usual one is based on orthogonal polynomials, the second approach is associated with an optimization procedure applied to the model complex compliance. Both lead to positive quadrature coefficients, which ensures the stability of the Andrade-DA model, as shown by propositions 1 and 3.

Gaussian quadrature. Various orthogonal polynomials can be used to evaluate the improper integral [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF] introduced by the diffusive representation of fractional derivatives.

Historically, the first one has been proposed in [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF], where a Gauss-Laguerre quadrature is chosen. Its slow convergence was highlighted and then corrected in [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF] with a Gauss-Jacobi quadrature. This latter method has been lastly modified in [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF], where alternative weight functions are introduced, yielding an improved discretization of the diffusive variable owing to the use of an extended interpolation range. Following this latter modified Gauss-Jacobi approach, while omitting the time and space coordinates for the sake of brevity, the improper integral ( 13) is then recast as

+∞ 0 φ(θ) dθ = +1 -1 1 -θ γ 1 + θ δ φ( θ) d θ ≃ L ℓ=1 μℓ φ( θℓ ), (38)
with the modified diffusive variable φ defined as

φ( θ) = 4 1 -θ γ-1 1 + θ δ+3 φ 1 - θ 1 + θ 2 ,
and where the weights and nodes {(μ ℓ , θℓ )} can be computed by standard routines [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF].

According to the analysis of [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF], Section 4, an optimal choice for the coefficients in [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF] is in the present case: γ = 3 -4α and δ = 4α -1. Following this approach, then by equating the series ( 38) and ( 16) that both approximate the term (13), the quadrature coefficients are chosen to be defined as

µ ℓ = 4 μℓ 1 -θℓ γ-1 1 + θℓ δ+3 , θ ℓ = 1 -θℓ 1 + θℓ 2 . ( 39
)
Optimization quadrature. Alternatively, the quadrature coefficients can be deduced from the model physical observables.

Note that as the quality factor ( 11) is defined as the ratio With reference to the quantities introduced in ( 19), then for a given number K of angular frequencies ω k , one defines the following objective function

Q(ω) = -Re[N]/ Im[N],
J {(µ ℓ , θ ℓ )} ; L, K = K k=1 κ(ω k ) κ(ω k ) -1 2 = K k=1 2 sin(πα) π L ℓ=1 µ ℓ θ 1-2α ℓ (iω k ) α θ 2 ℓ + iω k -1 2 (40)
to be minimized w.r.t parameters (µ ℓ , θ ℓ ) for ℓ = 1, . . . , L.

A straightforward linear minimization of ( 40) may lead to some negative parameters [START_REF] Blanc | Numerical modeling of transient poroelastic waves: the Biot-JKD model with fractional derivatives[END_REF][START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF] so that a nonlinear optimization with the positivity constraints µ ℓ ≥ 0 and θ ℓ ≥ 0 is preferred. The additional constraint θ ℓ ≤ θ max is also introduced to avoid the algorithm to diverge. These 3L constraints can be relaxed by setting

µ ℓ = µ ′ ℓ 2 and θ ℓ = θ ′ ℓ 2 and
solving the following problem with only L constraints min

{(θ ′ ℓ ,µ ′ ℓ )} J {(µ ′ ℓ 2 , θ ′ ℓ 2 )} ; L, K with θ ′ ℓ 2 ≤ θ max for ℓ = 1, . . . , L. (41)
As problem ( 41) is nonlinear and non-quadratic w.r.t. abscissae θ ′ ℓ , we implement the algorithm SolvOpt [START_REF] Kappel | An implementation of Shor's r-algorithm[END_REF][START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] based on the iterative Shor's method [START_REF] Shor | Minimization Methods for Non-Differentiable Functions[END_REF]. Initial values µ ′ 0 ℓ and θ ′ 0 ℓ used in the algorithm must be chosen with care; for this purpose we propose to use the coefficients obtained by the modified Jacobi method (39) for ℓ = 1, . . . , L

µ ′ 0 ℓ = 4 μℓ 1 -θℓ γ-1 1 + θℓ δ+3 , θ ′ 0 ℓ = 1 -θℓ 1 + θℓ . ( 42
)
Doing so, the required positivity constraints are satisfied by the initial guesses while it is expected that this choice already yields a satisfactory quadrature scheme as shown in [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF]. Finally, the angular frequencies ω k for k = 1, ..., K in ( 40) are chosen linearly on a logarithmic scale over a given optimization band [ω min , ω max ], i.e.

ω k = ω min ω max ω min k-1 K-1 . ( 43
)
Remark 1. In the proposed optimization method, both set of quadrature coefficients µ ℓ and θ ℓ are computed by minimization of the objective function J . In particular, the nodes θ ℓ are not imposed to be equidistributed according to [START_REF] Näsholm | On a fractional Zener elastic wave equation[END_REF] as it is the case in the commonly used approach [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]. This point will be returned to in Section 5.2.

Numerical scheme

A numerical scheme is proposed to compute the solution of system [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF]. Introducing a uniform grid with mesh size ∆x and time step ∆t, let U n j denote the approximation of the solution U (x j = j∆x, t n = n∆t) with j = 1, . . . , N x and n = 1, . . . , N t . Straightforward discretization of ( 24) typically yields to the numerical stability condition [START_REF] Blanc | Numerical modeling of transient poroelastic waves: the Biot-JKD model with fractional derivatives[END_REF] ∆t

≤ min ∆x c ∞ , 2 ̺(S)
.

As shown by Proposition 4, the usual CFL bound on the time step ∆t ≤ ∆x/c ∞ may be reduced as η decreases or A increases, which turns out to be detrimental to the numerical scheme. Moreover, as ̺(S) depends on the quadrature coefficients of the diffusive variable, the stability condition would in turn not depend only on meaningful physical quantities such as the maximum phase velocity c ∞ .

Splitting. Alternatively, we follow here the splitting approach analyzed in [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. To implement ( 24) numerically, one solves successively the propagative equation

∂U ∂t + A ∂U ∂x = 0 ( 44
)
and the diffusive equation

∂U ∂t = SU + F . ( 45
)
Due to the structure of matrix S, one defines from [START_REF] Diethelm | An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives[END_REF] the subvectors

U = σ, φ 1 , • • • , φ L T , F = F σ , J u γ 1,α F σ , • • • , J u γ L,α F σ T , (46)
and from ( 26) the submatrix

S =          -(J u η) -1 -AΓ(1 + α)J u -1 µ 1 • • • -AΓ(1 + α)J u -1 µ L -γ 1,α η -1 -θ 2 1 -Υ 1,α µ 1 • • • -Υ 1,α µ L . . . . . . . . . . . . -γ L,α η -1 -Υ L,α µ 1 • • • -θ 2 L -Υ L,α µ L          .
Having separated the two source terms, then equation ( 45) is equivalently recast in the form

         ∂v ∂t = F v , (47a)
∂U ∂t = S U + F . ( 47b
)
The discrete operators associated with the discretizations of ( 44) and ( 47) are respectively denoted by H p and H d . The operator H d depends explicitly on time when the forcing terms F v or F σ are non-zero, whereas H p remains independent on t. The socalled Strang splitting approach of [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF] is then used between time steps t n and t n+1 , for n = 0, . . . , N t -1, which requires to solve ( 44) and ( 45) with adequate time increments as, for j = 1, . . . , N x

U (1) j = H d (t n , ∆t/2) U n j , U (2) j = H p (∆t, j) U (1) , U n+1 j = H d (t n+1 , ∆t/2) U (2) j , (48)
with U (1) = U

(1) 1 . . . U

(1) Nx

T . Since the matrices A and S do not commute, an error associated with the splitting scheme is introduced [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. However, provided that H p and H d are at least second-order accurate and stable, then the time-marching scheme (48) constitutes a second-order accurate approximation of the original equation [START_REF] Birk | An improved non-classical method for the solution of fractional differential equations[END_REF].

Diffusive operator. The physical parameters do not vary with time, thus the matrix S does not depend on t. Owing to Property 3, one has 0 / ∈ sp(S) = {λ 1 , . . . , λ L }, and hence det S = 0. Freezing the forcing terms at t k , with k = n or n + 1, yields for a generic vector

U j = [v j , U j ] T H d (t k , ∆t/2) U j = v j + ∆t 2 F v (x j , t k ), e S ∆t 2 U j -I -e S ∆t 2 S -1 F (x j , t k ) T . ( 49
)
If there is no excitation, i.e. F v = F σ = 0, then integration (49) is exact. The matrix exponential entering the definition of the operator H d is computed using the method ♯2 in [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF] based on a (6/6) Padé approximation. Property 3 ensures that the computation of this exponential is stable.

Propagative operator. To integrate [START_REF] Hanyga | Fractional-order relaxation laws in non-linear viscoelasticity[END_REF], we use a fourth-order ADER (Arbitrary DERivative) scheme [START_REF] Schwartzkopff | Fast high-order ader schemes for linear hyperbolic equations[END_REF]. This explicit two-step and single-grid finite-difference scheme writes

U (2) j = U (1) j - ℓ=2 ℓ=-2 4 m=1 ϑ m,ℓ A ∆t ∆x m U (1)
j+ℓ := H p (∆t, j) U (1) , (50)
where the coefficients ϑ m,k are provided in Table 1. It satisfies the optimal stability condi-

tion c ∞ ∆t / ∆x ≤ 1. m = 1 m = 2 m = 3 m = 4 ℓ = -2 1/12 1/24 -1/12 -1/24 ℓ = -1 -2/3 -2/3 1/6 1/6 ℓ = 0 0 5/4 0 1/4 ℓ = 1 2/3 -2/3 -1/6 1/6 ℓ = 2 -1/12 1/24 1/12 -1/24
Table 1: Coefficients ϑ m,ℓ in the ADER-4 scheme (50)

Numerical results

Configuration

The homogeneous domain considered is 400 m-long and it is characterized by the physical parameters provided in Table 2 and which are consistent with experimentally-based values, see [START_REF] Bellis | Sensitivity of seismic measurements to frequency-dependent attenuation and upper mantle structure: an initial approach[END_REF] and the references therein. In this Section, one aims at assessing the overall performances of the proposed approach. In Section 5.2 we analyze the quadrature method in order to evaluate the model error, i.e. the error associated with the approximation of the Andrade model complex compliance and of associated observables. Section 5.3 is concerned with the validation of the numerical scheme for the wave propagation part. To do so, the numerical velocity field solution is compared to the semi-analytical Andrade-DA solution derived in Section 3.4.

ρ (kg/m 3 ) c ∞ (m/s) η (Pa.s) A (Pa -1 .s -α ) α 1200 2800 10 9 2•10 -10 1/3
Finally, we close the loop in Section 5.4 by comparing the semi-analytical Andrade solution to its numerically computed diffusive approximation-based version. Moreover, we provide a comparison between the theoretical phase velocity and its counterpart measured from the propagation simulations made. A similar comparison is made for the attenuation as a function of frequency and distance.

Validation of the quadrature methods

The angular frequency range of interest [ω min , ω max ] is defined by ω min = ω c /100 and ω max = 10 ω c for a given central source angular frequency ω c = 60 π. The choice of ω min is meant to promote the accuracy of the approximated model over long times. We choose K = 2 L while the parameter θ max introduced in ( 41) is set to θ max = √ 10 ω max to ensure a stable computation of the matrix exponential in (49). Observables of the Andrade model ( 5) are then compared to those of the Andrade-DA model [START_REF] Jackson | Grainsize-sensitive viscoelastic relaxation in olivine: Towards a robust laboratory-based model for seismological application[END_REF] on Figure 2 for the two quadrature methods discussed in Section 4.1. Large deviations are observed when the Gaussian quadrature is used, in particular on the attenuation function. On the contrary, an excellent agreement between the Andrade model and its optimized diffusive counterpart is obtained. Only slight differences can be observed at the scale of the figures within the optimization interval.

On Figure 3 are represented the L = 4 and L = 8 quadrature coefficients, i.e. nodes θ ℓ with corresponding weights µ ℓ , for the two methods considered. Note that, according to [START_REF] Pritz | Analysis of four-parameter fractional derivative model of real solid materials[END_REF], the values provided by the Gaussian approach are used as initial guesses in the minimization [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independent, part II[END_REF]. The scaled optimization angular frequencies √ ω k for k = 1, . . . , K are also shown for the purposes of comparison. Remarkably, the computed optimal nodes do not coincide with equidistributed nodes along the optimization frequency-band, a repartition which is prescribed in the commonly employed approach of [START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF]. 2. Comparison with their approximated counterparts for L = 4 memory variables and using either the modified Gauss-Jacobi approach or the proposed optimization method.

1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3 1E-1 1E+ 0 1E+ 1 f (Hz) Q exact Gaussian optimization (b) Quality factor Q(f ) 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3
Vertical dotted lines delimit the optimization frequency-band. The horizontal solid line in panel (c) denotes the high-frequency limit c ∞ .

The corresponding model error defined as | κ(ω) κ(ω) -1| and associated with the minimization problem [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independent, part II[END_REF] is displayed in Figure 4, for L = 4 (Fig. 4a) and L = 8 (Fig. 4b) diffusive variables. For a given quadrature method, the results are clearly improved as L increases. For a given L, the optimization provides more accurate results compared to the Gaussian quadrature over the frequency band of interest which is delimited by vertical dotted lines.

1E-1 1E+ 0 1E+ 1 1E+ 2 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3 nodes weights Gaussian optimization (a) L = 4 1E-2 1E-1 1E+ 0 1E+ 1 1E+ 2 1E+ 3 1E-2 1E+

Validation of the numerical scheme

While F σ = 0 in (22b), the source in (22a) is imposed at point x s as F v (x, t) = F (t) δ(x -x s ) where F (t) is the function with regularity C 6 that is defined by original version, is lower than in the elastic case, which explains the observed delay. Figure 5b shows a seismogram corresponding to the Andrade-DA model in order to highlight attenuation and dispersion of the waveform.

F (t) =        4 m=1 a m sin (b m ω c t) if 0 ≤ t ≤ 1 f c , 0
Considering the computed Andrade-DA model, Figure 6 compares the semi-analytical and the numerical velocity field solutions corresponding to equation [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] where κ mod = κ and to [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF] respectively. Figure 6b presents the relative spatial L 2 -norm error at final time t f between these two solutions for various discretizations, varying the numbers of nodes in the interval N x = 50 to 6400. These convergence measurements show that order 2 is reached, confirming the theoretical results of Section 4.2.

Validation of the overall approach

To assess the performances of the overall approach, we now confront the results from the propagation simulations made to the Andrade model. Firstly, Figure 7 plots the relative spatial L 2 -norm error at final time t f between the velocity field solution [START_REF] Rekik | Optimization of the collocation inversion method for the linear viscoelastic homogenization[END_REF] where κ mod = κ and this obtained from [START_REF] Yuan | A numerical scheme for dynamic systems containing fractional derivatives[END_REF] for various values of the discretization parameter N x shown to be well-suited for a numerical implementation.

The system at hand is investigated and it is demonstrated that its well-posedness requires the positiveness of the weights associated with the quadrature scheme. To compute the quadrature coefficients, two numerical methods are combined: a polynomial Gaussian approach to get an initial guess jointly with a constrained optimization to approximate the Andrade model complex compliance over a frequency-band of interest. It is shown and the efficiency of the proposed approach is highlighted. The main point of this article is that using a diffusive approximation of a fractional derivative term, entering a given viscoelastic constitutive equation, yields a sound mathematical model, that is also easily tractable numerically to perform wave propagation simulations.

To focus on this message, a simple but realistic fractionally-damped viscoelastic model within a unidimensional and homogeneous configuration has been considered. Its dynamical behavior is described by a first-order hyperbolic system which extension to higher spatial dimensions or heterogeneous media is straightforward. The main limitation of this study concerns the numerical scheme employed to solve the wave propagation problem for two reasons: (i) The splitting approach is of order 2 which constitutes an intrinsic limiting factor even if the employed ADER scheme is of order 4. (ii) The numerical scheme has been developed for a (piecewise)-homogeneous body. Yet, efficient numerical methods are currently available and can be directly employed to perform time-domain simulations within a higher-order scheme that is also valid for heterogeneous configurations. Improving the method along these lines constitute the main focus for future work. Alternatively, arbitrary-shaped material discontinuities within piecewise-homogeneous 2D Andrade media can be handled using an immersed interface method [START_REF] Chiavassa | Time domain numerical modeling of wave propagation in 2D heterogeneous porous media[END_REF].

Another line of research concerns extension of the proposed approach to other fractional viscoelastic model, such as the fractional Kelvin-Voigt model [START_REF] Caputo | Linear models of dissipation whose Q is almost frequency independent, part II[END_REF][START_REF] Caputo | Wave simulation in biologic media based on the Kelvin-Voigt fractional-derivative stress-strain relation[END_REF] or the fractional Zener model [START_REF] Pritz | Analysis of four-parameter fractional derivative model of real solid materials[END_REF][START_REF] Näsholm | On a fractional Zener elastic wave equation[END_REF]. More sophisticated models could also be investigated, such as nonlinear fractional viscoelasticity [START_REF] Hanyga | Fractional-order relaxation laws in non-linear viscoelasticity[END_REF] or nonlocal models in space [START_REF] Hanyga | Spatially fractional-order viscoelasticity, non-locality, and a new kind of anisotropy[END_REF].

Figure 1 :

 1 Figure 1: Behaviors of various viscoelastic models derived from (4): Maxwell model (A = 0) and Andrade model (α = 1/3, with A = 10 -10 Pa -1 .s -α and A = 2.10 -10 Pa -1 .s -α ).

  equations (22a) and[START_REF] Torvik | On the appearance of the fractional derivative in the behavior of real materials[END_REF], or by the Andrade-DA model, i.e. equations (22a) and (22b), together with equation[START_REF] Andrade | On the validity of t 1/3 law of flow of metals[END_REF]. Corresponding semi-analytical solutions are sought in order to validate the ensuing numerical simulations of wave propagation. It is assumed F σ = 0 and excitation F v (x, t) = F (t)δ(x-x s ) at source point x s with time evolution F . Applying space-time Fourier transforms and their inverses leads to the stress field solution in the form of σ

  then obtaining a good fit on the latter does not imply a satisfying approximation of the function N itself. In other words, optimizing an objective function based on Q(ω) might yield a poor approximation of the model constitutive equation. Therefore, a direct optimization of the available Andrade model complex compliance N is preferred.

  Creep function χ(t)

Figure 2 :

 2 Figure 2: Exact observables of the Andrade model with physical parameter values provided in Table 2. Comparison with their approximated counterparts for L = 4 memory variables

Figure 3 :

 3 Figure 3: Set of quadrature coefficients for the two approaches considered. The L points are plotted with abscissae and ordinates corresponding respectively to node θ ℓ and weight µ ℓ values for ℓ = 1, . . . , L. Vertical dashed lines are plotted at the abscissae corresponding to the K = 2L scaled optimization angular frequencies values √ ω k .

  f c = ω c /2 π = 30 Hz and parameters b m = 2 m-1 , a 1 = 1, a 2 = -21/32, a 3 = 63/768 and a 4 = -1/512. The associated frequency bandwidth is highlighted in Fig. A.9b. Moreover, the domain is discretized with N x = 400 nodes and the diffusive approximation is computed by constrained optimization with L = 4 memory

Figure 4 :

 4 Figure 4: Computed error | κ(f )κ(f ) -1| quantifying the discrepancy between the Andrade model complex compliance and its diffusive approximation. Comparison between the modified Gauss-Jacobi approach and the proposed optimization method. Vertical dotted lines delimit the optimization frequency-band.

Figure 5 Figure 5 :

 55 Figure 5 displays snapshots of forward propagating waves from the source point x s = 200. The numerical solutions associated with various values of the attenuation parameters in (20) are plotted on Fig. 5a; namely Hooke model (i.e. purely elastic case which may be obtained in the limit η = +∞ and setting A = 0), Maxwell model (A = 0, η = 10 9 ), and Andrade-DA model (A = 2 • 10 -9 , η = 10 9 ). As predicted by the dispersion analysis of sections 2.3 and 2.4, the phase velocity of the Andrade-DA model, as this of its

1EFigure 6 :

 6 Figure 6: Validation of the numerical scheme for the Andrade-DA model. (a) Snapshots of velocity fields at final time t f for the semi-analytical Andrade-DA solution and its numerical counterpart for L = 4. (b) Relative error in spatial L 2 -norm between these two solutions for a varying value of the discretization parameter N x .

Figure 7 :

 7 Figure 7: Discrepancy between the numerical Andrade-DA solution and the semianalytical Andrade solution. The relative error in spatial L 2 -norm at final time t f between the associated velocity fields solutions are plotted for a varying value of the discretization parameter N x and of the quadrature parameter L.

  Attenuation ζ(f )

Figure 8 :

 8 Figure 8: Discrepancy between exact observables of the Andrade model and their counterparts measured from numerical wave propagation simulations using the Andrade-DA model. These quantities are compared for the values L = 4 and L = 8 of the quadrature parameter and represented within the exciting source frequency bandwidth.

  Figure A.9: Source employed in the wave propagation simulations. (a) Time-domain signal and (b) frequency spectrum with vertical dashed line plotted at abscissa f c = 30 Hz.

Table 2 :

 2 Chosen physical parameters in the Andrade model (4).

dx,

Acknowledgements. The work of Abderrahmin Ben Jazia has been funded by the Ecole Centrale de Marseille, France, for which special thanks are addressed to Guillaume Chiavassa. The authors are thankful to Emilie Blanc for fruitful discussions

Appendix A. Exciting source signal

The Fourier transform of the time-domain source signal (51) reads

For the chosen values [2] T. L. Szabo, Time domain wave equations for lossy media obeying a frequency power law, J. Acoust. Soc. Am. 96 (1994) 491-500.