domain: sdu.stu.gp info math phys

text:
Yder Masson
email: yder.masson@gmail.com

Distrib utional finite-differ ence modelling of seismic w av es

Keywords: Numerical modelling, Computational seismology, Wave propagation

This study introduces a distributional finite-difference method (DFDM) for modelling the propagation of elastic waves in heterogeneous media in the time domain. DFDM decomposes the modelling domain into multiple elements that can have arbitrary sizes. When large elements are used, the proposed method closely resembles the finite-difference method because the wavefield is updated using operations involving band diagonal matrices only. Thus DFDM is computationally efficient. When smaller elements are used, DFDM looks closer to the finite-element or the spectral element methods and permits to mesh complicated structures. A complete multidomain algorithm for modelling elastic wave propagation in arbitrarily heterogeneous media is presented. The algorithm's stability is discussed, and the usual Courant condition governs the stability of the proposed scheme. Numerical examples show that the proposed algorithm accurately accounts for free surfaces, solid-fluid interfaces and accommodates non-conformal meshes in their basic for m. Seismog rams obtained using the proposed method are compared to those computed using analytical solutions and the spectral element method. To achieve comparable accuracy, DFDM requires fewer points per wavelength than the spectral element method, for example.

(SEM; [START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF][START_REF] Nissen-Meyer | Axisem: broadband 3-D seismic wavefields in axisymmetric media[END_REF], boundary element methods (BEM; [START_REF] Bouchon | Calculation of synthetic seismograms in a laterall y v arying medium b y the boundary element-discrete w avenumber method[END_REF][START_REF] Papageorgiou | A discrete wavenumber boundary element method for study of the 3-D response 2-D scatterers[END_REF][START_REF] Ba | Soil-structure interaction in transversely isotropic layered media subjected to incident plane SH waves[END_REF], finite-volume methods (FVM; [START_REF] Dormy | Numerical simulation of elastic wave propagation using a finite volume method[END_REF][START_REF] Dumbser | Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D[END_REF] and Discontinuous Galerkin methods (DGM; de la Puente et al. 2008 ;[START_REF] Hermann | Quantitative accuracy analysis of the discontinuous Galerkin method for seismic wave propagation[END_REF][START_REF] Étienne | An hp-adaptive discontinuous galerkin finite-element method for 3-D elastic wave modelling[END_REF][START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF][START_REF] Bonnasse-Gahot | Hybridizable discontinuous Galerkin method for the 2-D frequency-domain elastic wave equations[END_REF] ). In the elements methods, the discrete algorithm is constructed from the wave equation's weak formulation, making it easier to account for various boundary conditions accurately, model discontinuities and free surfaces with complex shapes, and perform domain decomposition. One should note that SEM has enabled important advances in seismology [START_REF] Tape | Adjoint tomography of the souther n Califor nia Cr ust[END_REF][START_REF] French | Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography[END_REF][START_REF] Bozda G | Global adjoint tomography: firstgeneration model[END_REF] ) and earthquake physics [START_REF] Ampuero | Upper limit on damage zone thickness controlled b y seismo genic depth, in Fault Zone Dynamic Processes: Evolution of Fault Properties During Seismic Rupture[END_REF][START_REF] Weng | The dynamics of elongated earthquake ruptures[END_REF], 2020 ). In recent years, DGM and Arbitrary high-order DERi v ati v e DGM (ADER-DG) hav e become increasingly popular to model dynamic earthquake rupture (Wilcox et al. 2010 ;[START_REF] Breuer | Sustained P etascale P erformance of Seismic Simulations with SeisSol on SuperMUC[END_REF]. This is because DGM do not impose any field continuity across their boundaries and allows include non-linear interface conditions. It is also possible to implement local time stepping [START_REF] Uphoff | Extreme scale multi-physics simulations of the Tsunamigenic 2004 Sumatra megathrust earthquake[END_REF] and to use adaptive mesh refinement [START_REF] Kunkel | p4est: scalable algorithms for parallel adaptive mesh refinement on forests of octrees[END_REF].

In this work, a distributional finite-difference method (DFDM) is introduced. The principal moti v ation behind this approach is to obtain a method that leads to codes with a simple structure similar to FDM while taking advantage of the accurate treatment of the boundary conditions of fered b y SEM or DGM. In DFDM, the discrete operators e v aluating the partial deri v ati ves are obtained using the weighted residuals method [START_REF] Salih | Weighted residual methods[END_REF][START_REF] Reddy | Introduction to the Finite Element Method[END_REF] which is in many ways similar to the weak-formulation approach; they can be seen as discrete equi v alents of the distributional deri v ati ve

f ( d) ( x) = +∞ -∞ f ( d) ( l ) δ( l -x) d l . (1a) = ( -1) d +∞ -∞ f ( l ) δ ( d ) ( l -x)dl , ( 1b
)
where B j ( x ) ( j = 1, N B ) denote the N B basis functions with support inside and f B j ( j = 1, N B ) are the expansion coefficients stored in the v ertical v ector f B . The d th deri v ati ve of f ( x ) is gi ven b y

f ( d) ( x ) = N B j= 1 f B j B ( d) j ( x ) , (3)
where B

( d) j ( x ) is the deri v ati ve of order d of the basis function B j ( x ). Let the desired approximation ˜ f ( d) ( x ) of the exact deri v ati ve f ( d ) ( x ) be represented in an arbitrary basis β as the linear combination ˜ f ( d)

( x ) = N β j= 1 f ( d) β j β j ( x ) , ( 4
)
where β j ( x ) denote the N β basis functions with support inside and f ( d) β j are the expansion coefficients stored in the vertical vector f ( d) β . To obtain ˜ f ( d) ( x ), one needs to determine the expansion coefficients in the vector f ( d) β . When the approximation is equal to the true deri v ati ve, w e ha ve ˜ f ( d)

( x ) = f ( d) ( x ) ∀ x ∈ .
(5)

Assuming that the expansion coefficients f B j and f ( d) β j in eqs ( 3) and ( 4) are known, it is not straightforward to check that eq. ( 5) is satisfied because the two functions are represented in different bases. One way to test the equality in eq. ( 5) is to measure f ( d ) ( x ) and ˜ f ( d) ( x ) by multiplying both sides of eq. ( 5) by the same arbitrary test function φ i and integrating over the domain . This gives

x + x - φ i ˜ f ( d) dx = x + x - φ i f ( d) dx (6a) = ( -1) d x + x - φ ( d ) i f dx + d-1 p= 0 ( -1) d-p-1 φ ( d -p-1) i f ( p) x + x - , (6b)
where the right-hand side of eq. ( 6b) is obtained by repeated integration by part, assuming that φ and f are continuously differentiable. When the test function φ i is thought as a Dirac delta function δ( xx 0 ) with the fundamental property δ( xx 0 ) f ( x )dx = f ( x 0 ), eq. ( 6) simply verifies that ˜ f ( d) ( x ) and f ( d ) ( x ) have the same value at point x 0 . If this is true for any point x 0 in : [ x -, x + ], then eq. ( 5) is satisfied. Thus, a strategy for computing the approximate deri v ati ve ˜ f ( d) ( x ) is to search for the optimal expansion coefficients f ( d) β j that satisfy the constraint in eq. ( 6) for a set of N φ test functions φ i ( x ) ( i = 1, N φ ). Because ˜ f ( d) ( x ) is represented using a finite number of unknown coefficients f ( d) β j , one only needs to use a finite number of test functions to determine these coefficients. This general approach is called the method of weighted residuals. By substituting f ( d ) ( x ) and ˜ f ( d) ( x ) in eq. ( 6) with their expansions in eqs ( 3 ) and ( 4 ), respecti vel y, we obtain

N β j= 1 f ( d) β j x + x - φ i β j dx = N B j= 1 f B j x + x - φ i B ( d) j dx (7a) = ( -1) d N B j= 1 f B j x + x - φ ( d ) i B j dx + d-1 p= 0 ( -1) d-p-1 N B j= 1 f B j φ ( d -p-1) i B ( p) j x + x - (7b)
or, in matrix form

Mf ( d) β j = K ( d) f B = K * ( d) + Q -( d) + Q + ( d) f B , (8)
where 25). (e) The dual bases obtained from the Hermite bases in (b) using eq. ( 25). (f) The dual bases obtained from the B-spline bases in panel (c) using eq. ( 25). (g) The dual bases obtained from the polynomial bases in panel (a) using eq. ( 20). (h) The dual bases obtained from the Hermite bases in panel (b) using eq. ( 20). (i) The dual bases obtained from the B-spline bases in panel (c) using eq. ( 20). and f, g

M = M φ,β = ( M ) i, j = φ i , β j x + x - (9a)
K ( d) = M ( d) φ,B = K ( d) i, j = φ i , B ( d) j x + x - (9b) K * ( d) = K * ( d) φ,B = K * ( d) i, j = ( -1) d φ ( d ) i , B j x + x - (9c) Q ± ( d) = Q ± ( d) φ,B = Q ± ( d) i, j = ± d-1 p= 0 ( -1) d-p-1 ± ( d-p-1) q ± ( p) T = ± d-1 p= 0 ( -1) d-p-1 φ ( d -p-1) i ( x ± ) B ( p) j ( x ± ) (9d)
± ( p) = ± ( p) φ = ± ( p) i = φ ( p) i ( x ± ) (9e) q ± ( d) = q ± ( d) B = q ± ( d) i = B ( d) i ( x ± ) (9f) K ( d) = K * ( d) + Q -( d) + Q + ( d) , ( 9g
) (a) (b) (c) (f) (d) (e) (i) (g) (h)
x + x - = x + x - f ( x ) g ( x )dx . ( 10
)
In eqs ( 8) and ( 9), the matrix M is often referred to as the mass matrix, the matrices K ( d) and K * ( d) involving the deri v ati ves of the basis functions are usually called stiffness matrices and the matrices Q - ( d) and Q + ( d) account for the left and the right boundary values, respectively. The linear system in eq. ( 8) can be solved to obtain the expansion coefficient vector

f ( d) β = D ( d) φ,β,B f B , (11)
where the operator

D ( d) φ,β,B = M -1 K ( d) = M -1 K * ( d) + Q -( d) + Q + ( d) . ( 12
)
is the DFD operator of order d . Provided that the inverse matrix M -1 (or some generalized inverse of the matrix M ) exists, the DFD operator in eq. ( 12) can be used to compute approximations of the d th deri v ati ve of f ( x ). The DFD operator D ( d) φ,β,B acts on the e xpansion v ector f B representing f ( x ) in basis B and returns the expansion vector f ( d) β representing ˜ f ( d) ( x ) in basis β. In this study, the DFD operator D ( d) φ,β,B will be constructed using the piecewise polynomial bases presented in the Appendix and pictured in Figs 1 (a)-(c) and (A1 ). All the basis functions considered have compact support (i.e. equal to zero outside a finite interval). It is crucial to obtain efficient algorithms with a finite-difference-like structure. Indeed, when those bases are used, the matrices K ( d) and K * ( d) are band diagonal (i.e. similar to finite-difference matrices but with different coefficients). The matrix M is band diagonal and symmetric. It can be factorized using Cholesky factorization. Thus the multiplications by the matrices K ( d) , K * ( d) , and M can be implemented very ef ficientl y. They have a comparable computational cost. The elements of the vectors ± ( p) and q ± ( d) are zero except the first one or the last one. The elements of the matrices Q ± ( d) are zero except the first one (upper-left) or the last one (bottom-right). It simplifies the treatment of the boundary conditions. Prior ending this section, one should make two informative remarks:

(i) F irst, w hen the basis functions B j ( x ) and β j ( x ) are interpolant, that is, when the expansion coefficients are equal to the values of the functions at some interpolation or collocation points, that is we have

f B j = f ( x i )( x i = 1 , N B ) and f ( d) β j = ˜ f ( d) ( x j ) ( x j = 1 , N β )
, then the operator in eq. ( 12) is a true finite-difference operator.

(ii) Secondly, when the test functions φ i are taken as the dual basis functions

δ i ( x ) = N β j= 1 ( M -1 β ) i j β j ( x ) (e.
g. in F ig. 1 ), w here the ( M -1 β ) i j are the components of the inverse of the matrix

M β = ( M β ) i, j = β i , β j x + x -
, the expansion coefficients f ( d) β j in eq. ( 4) can be obtained directly using f ( d) d) dx . Further, when the basis functions β i are also interpolant, we have

β j = x + x - δ j f (
f ( d) β j = f ( d) ( x j ) = x + x - δ j f ( d) dx , ( 13
)
which is a discrete analogue to Dirac's formula in eq. ( 1b) where the Dirac delta function δ( x ) has been replaced with its approximation δ i ( x ) defined above.

These two observations outline the relationship between the DFD operators and classic finite-difference operators. In the next section, the DFD operator D ( d) φ,β,B is modified to account for prescribed boundary conditions.

DFD operators with prescribed boundary values

In this section, the DFD operator in eq. ( 12) is reformulated in a way that is more convenient to impose prescribed boundary conditions and to perform domain decomposition. At the limits x -and x + of domain , we may express the boundary values of f ( d ) ( x ) using the weighted averages

f ( d) ( x -) = ( 1 -γ -) f ( d) ( x -) + γ -f -( d) (14a) f ( d) ( x + ) = ( 1 -γ + ) f ( d) ( x + ) + γ + f + ( d) , ( 14b
)
where γ -and γ + are the left and right weights, respecti vel y, f ( d) ( x -) and f ( d) ( x + ) are internal values taken inside domain that can be computed using eq. ( 2), and, f - d) and f ( d) are external values taken from the domains -: [ . . . , x -] and + : [ x + , . . . ] neighbouring domain to the left and to the right, respecti vel y. Assuming that the function f ( x ) has representations similar to that in eq. ( 2) in domains -and + , the external values may be computed using

(
+
f -( d) = f ( d) ( x + ) -= N B j= 1 f B j B j ( x + ) - = q + ( d) T f B - (15a) f + ( d) = f ( d) ( x -) + = N B j= 1 f B j B j ( x -) + = q -( d) T f B + , (15b)
where | . . . | specifies that the quantities within the brackets are associated with domain , in particular, x -= x + -and x + = x -+ .

Starting back from eq. ( 6) and following the steps thereafter using the substitutions in eq. ( 14) shows that the expansion coefficients in eq. ( 11) can be computed using the modified expression

f ( d) β = D ( d) φ,β,B,γ-,γ+ f B + d-1 p= 0 b -( p) f -( p) + b + ( p) f + ( p) , ( 16
)
where the original DFD operator D ( d) φ,β,B in eq. ( 12) transforms to

D ( d) φ,β,B,γ-,γ+ = M -1 K * ( d) + ( 1 -γ -) Q -( d) + ( 1 -γ + ) Q + ( d) = M -1 K ( d) -γ -Q -( d) -γ + Q + ( d) (17)
and the boundary operators are

b ± ( p) = b ± ( p) φ,β,γ± = ±γ ± ( -1) d-p-1 M -1 ± ( d-p-1) . ( 18
)
When setting γ -= γ + = 0 in eqs ( 17) and ( 18), eq. ( 16) reduces to eq. ( 11) and the DFD operator in eq. ( 17) is the same as the one in eq. ( 12). In this case, the boundary values f ( d) ( x -) and f ( d) ( x + ) are implicitly specified by the expansion coefficients in vector f B . When taking γ -= 1 and/or γ + = 1, the DFD operator in eq. ( 17) returns the expansion coefficients vector f ( d) β assuming that f ( d) ( x -) = 0 and/or f ( d) ( x + ) = 0, respecti vel y. The left and right boundary values can be imposed independently by adding the products b ± ( p) f ± ( p) to the vector f ( d) β . The external values f -( d) and f + ( d) may be set arbitrarily or computed using eq. ( 14) knowing the expansion coefficients vectors f B - and f B + in the neighbouring domains -and + , respecti vel y.

Adjoint DFD operators

In this section, we construct adjoint pairs of DFD operators acting on functions represented in the first basis ˆ B 1 and returning their deri v ati ve represented in the second basis ˆ B 2 , and inversely. Such pairs of DFD operators allow to repeatedl y dif ferentiate the field variables (e.g. velocity, displacement or stresses, etc.) represented in different bases, going back and forth from one variable to another in a consistent manner. Further, they maintain the self-adjoint nature of the wave equation when it is discretized, which is essential to achieve numerical reciprocity and obtain numerical schemes that are conditionally stable (see Section 2.1). This is because adjoint operators lead to a block-antidiagonal first-order linear system that preserves the intrinsic staggered property [START_REF] Virieux | P-sv wave propagation in heterogeneous media: velocitystress finite-difference method[END_REF] ) using dedicated basis function.

Adjoint DFD operators (i.e. where the second operator is the transpose of the first one or its opposite) can be obtained directly from the general expressions in eqs ( 17) and ( 18) by using specific basis functions for φ, β and B . Let f ( x ) and g ( x ) be two functions represented in the bases B 1 i ( i = 1, N 1 ) and B 2 i ( i = 1, N 2 ) through their expansion coefficient vectors f 1 and g 2 , respecti vel y. By constructing and factorizing (e.g. using Cholesky factorization) the mass matrices associated with the bases B 1 and B 2 as

M 11 = M B 1 ,B 1 = L 1 L T 1 (19a) M 22 = M B 2 ,B 2 = L 2 L T 2 , ( 19b
)
we can construct the two sets of orthogonal basis functions (illustrated in Fig. 1 ) 12 acts on the vector ˆ g 2 representing g ( x ) in basis ˆ B 2 and returns the vector ˆ g ( d) 1 representing the d th deri v ati ve of g ( x ) in basis ˆ B 1 . We have:

ˆ B 1 i = N 1 j= 1 ( L -1 1 ) i j B 1 j ∀ i = 1 , N 1 (20a) ˆ B 2 i = N 2 j= 1 ( L -1 2 ) i j B 2 j ∀ i = 1 , N 2 . ( 20b
ˆ f ( d) 2 = ˆ D ( d) 21 ˆ f 1 + d-1 p= 0 ˆ b -( p) 2 f -( p) + ˆ b + ( p) 2 f + ( p) (21a) ˆ g ( d) 1 = ˆ D ( d) 12 ˆ g 2 + d-1 p= 0 ˆ b -( p) 1 g -( p) + ˆ b + ( p) 1 g + ( p) . ( 21b
)
If needed, the expansion vectors in the original bases B 1 and B 2 may be obtained using

f 1 = ( L -1 1 ) T ˆ f 1 f ( d) 2 = ( L -1 2 ) T ˆ f ( d) 2 g 2 = ( L -1 2 ) T ˆ g 2 g ( d) 1 = ( L -1 1 ) T ˆ g ( d) 1 . ( 22
)
To enforce the adjoint property ˆ D

ˆ D

( d)

21 = D ( d) ˆ B 2 , ˆ B 2 , ˆ B 1 ,α-,α+ = L -1 2 D * ( d) ( L -1 1 ) T (23a) ˆ D ( d) 12 = D ( d) ˆ B 1 , ˆ B 1 , ˆ B 2 , ( 1 -α-) , ( 1 -α+ ) = ( -1) d L -1 1 D * ( d) T ( L -1 2 ) T = ( -1) d ( ˆ D ( d) 21 ) T (23b) ˆ b ± ( p) 2 = b ± ( p) ˆ B 2 , ˆ B 2 ,α± = L -1 2 b * ± ( p) 2 (23c) ˆ b ± ( p) 1 = b ± ( p) ˆ B 1 , ˆ B 1 , ( 1 -α±) = L -1 1 b
Figure 2. Illustration of the geometrical mapping between the physical space coordinates ( x , y ) and the computational space coordinates ( x , y ).

respecti vel y. Note that choosing the complementary boundary weights 24a) and (24b). We shall see that when modelling wave propagation in a single domain, the condition ˆ D

γ -= α -, γ + = α + for ˆ D ( d) 21 and γ -= (1 -α -), γ + = (1 -α + ) for ˆ D ( d) 12 enforces the property D * ( d) B 2 , B * 2 ,B 1 = ( -1) d D * ( d) B 1 , B * 1 ,B 2 in eqs (
( d) 12 = ( -1) d ( ˆ D ( d)

21

) T is sufficient to obtain stable numerical schemes (see Section 2.1), Neumann or Dirichlet boundary conditions can be accounted for by setting α ± = 0 or α ± = 1, while the boundary values f ± ( p) and g ± ( p) may be zeroed out or used to impose an external source. When modelling wave propagation in multiple domains

i ( i = 1, N d ), the additional constrain | α + | i = (1 -| α -| i+ 1 ) which translates the continuity condition | f ( d) ( x + ) | i = | f ( d) ( x -) | i+ 1 ,
is needed to ensure that the global DFD operator (that e v aluates the deri v ati ves in all domains) stays adjoint). The operators in eq. ( 23) are the basic bricks of the proposed algorithm, they can be combined to compute partial deri v ati ve in higher dimensions as discussed in Section 2.4.

DFD operators in higher dimensions

In this section the DFD operators obtained in the previous sections are extended to higher dimensions. A general formula for the DFD operators in two dimension is presented, the operators in higher dimensions can be derived in a similar manner. The objective is to compute the deri v ati ves of a 2-D function f ( x , y ) in the rectangular domain (e.g. as pictured in Fig. 2 )

: = ( x , y ) : x -≤ x ≤ x + , y -≤ y ≤ y + . ( 26
)
Given the 1-D sets of basis functions

X B i ( x ) ( i = 1 , N X B , X β i ( x ) ( i = 1 , N X β ), X φ i ( x ) ( i = 1 , N X φ ) defined on the 1-D domain x : [ x -, x + ] and, Y B i ( y ) ( i = 1 , N Y B ), Y β i ( y ) ( i = 1 , N Y β ), Y φ i ( y ) ( i = 1 , N Y φ ) defined on the 1-D domain y : [ y -, y + ],
we can construct the 2-D sets of basis functions

B i, j ( x , y ) = X B i ( x ) Y B j ( y ) 1 ≤ i ≤ N X B 1 ≤ j ≤ N Y B ( 27a
)
β i, j ( x , y ) = X β i ( x ) Y β j ( y ) 1 ≤ i ≤ N X β 1 ≤ j ≤ N Y β (27b) φ i, j ( x , y ) = X φ i ( x ) Y φ j ( y ) 1 ≤ i ≤ N X φ 1 ≤ j ≤ N Y φ . ( 27c
)
Following the approach in Section (2.1), we represent the function

f ( x , y ) = N X B i= 1 N Y B j= 1 F B i j B i j ( x , y ) = N X B i= 1 X B i ( x ) N Y B j= 1 Y B j ( y ) F B i j (28)
using N X B × N Y B basis functions B ij ( x , y ) with support inside the rectangular domain in eq. ( 26) and expansion coefficients F B i j stored in matrix F B . The exact partial deri v ati ves of f ( x , y ) are

f ( d x ,d y ) ( x , y ) = N X B i= 1 N Y B j= 1 F B i j B ( d x ,d y ) i j ( x , y ) = N X B i= 1 X ( d x ) B i ( x ) N Y B j= 1 Y ( d y ) B j ( y ) F B i j , ( 29
)
where

B ( d x ,d y ) i j ( x , y ) = ∂ d x + d y ∂ x d x ∂ y d y B i j ( x , y ) = ∂ d x ∂ x d x X B i ( x ) ∂ d y ∂ y d y Y B j ( y ) = X ( d x ) B i ( x ) Y ( d y )
B j ( y ) denotes the partial deri v ati ve of the basis function B ij ( x , y ), and, d x and d y are the orders of the partial deri v ati ves in direction x and y , respecti vel y. We w ant to obtain the approximation

˜ f ( d x ,d y ) ( x , y ) = N X β i= 1 N Y β j= 1 ˜ f ( d x ,d y ) β i j β i j ( x , y ) = N X β i= 1 X β i ( x ) N Y β j= 1 Y β j ( y ) ˜ f ( d x ,d y ) β i j (30) of f ( d x ,d y ) ( x , y ) represented using N X β × N Y β basis functions β ij ( x , y
) with support inside the rectangular domain in eq. ( 26) and unknown expansion coefficients ˜ f

( d x ,d y ) β i j stored in matrix ˜ f ( d x ,d y ) β . Starting with the target equality ˜ f ( d x ,d y ) ( x , y ) = f ( d x ,d y ) ( x , y ) , ( 31
)
multiplying both sides of eq. ( 31) by the test function φ kl = φ (0 , 0) kl in eq. ( 27c) and integrating over the rectangular domain in eq. ( 26)

x + x - y + y - φ kl ˜ f ( d x ,d y ) d x d y = x + x - y + y - φ kl f ( d x ,d y ) d x d y , ( 32
)
and replacing f ( d x ,d y ) and ˜ f ( d x ,d y ) with their expansions in eqs ( 29) and ( 30), respecti vel y, we obtain

N X β i= 1 x + x - X φ k X β i d x N Y β j= 1 y + y - Y φ l Y β j d y ˜ f ( d x ,d y ) β i j = N X B i= 1 x + x - X φ k X ( d x ) B i d x N Y B j= 1 y + y - Y φ l Y ( d y ) B j d y F B i j , (33)
or, in matrix form

M x ˜ f ( d x ,d y ) β ( M y ) T = K ( d x ) x F B ( K ( d y ) y ) T , ( 34
)
where the matrices

M x = M X φ ,X β , M y = M Y φ ,Y β , K ( d x ) x = K ( d x )
X φ ,X B and K

( d y ) y = K ( d y )
Y φ ,Y B are similar to those in the 1-D case and can be computed using eqs ( 9a ) and ( 9b ). The expansion coefficients ˜ F β i j can be computed by solving the linear system in eq. ( 34), and the 2-D DFD operator

D ( d x ,d y )
φ,β,B can be defined using

˜ f ( d x ,d y ) β = D ( d x ,d y ) φ,β,B F B . ( 35
)
By defining the action of the tensor product between two matrices A and B on a matrix X as

( A ⊗ B ) X = AXB T = B ( AX ) T T , ( 36
)
where ⊗ denotes the tensor product, and, using the identity

( A ⊗ B ) = ( A ⊗ I ) ( I ⊗ B ) = ( I ⊗ B ) ( A ⊗ I ) , ( 37
)
where I is the identity matrix, we can express the 2-D DFD operator in eq. ( 35) as

D ( d x ,d y ) φ,β,B = M x -1 K ( d x ) x ⊗ M y -1 K ( d y ) y (38a) = D ( d x ) X φ ,X β X B ⊗ D ( d y ) Y φ ,Y β ,Y B (38b) = D ( d x ) X φ ,X β X B ⊗ I I ⊗ D ( d y ) Y φ ,Y β ,Y B (38c) = I ⊗ D ( d y ) Y φ ,Y β ,Y B D ( d x ) X φ ,X β X B ⊗ I (38d)
which is simply the tensor product between the 1-D DFD operators D

( d x )
X φ ,X β X B and D

( d y ) Y φ ,Y β ,Y B .
When acting on a matrix X the operator ( A ⊗ I ) in eq. ( 37) applies the matrix A to all columns of matrix X while the operator ( I ⊗ B ) applies the matrix B to all rows of matrix X . Thus the 2-D DFD operators simpl y appl y the 1-D DFD operators D

( d y ) Y φ ,Y β ,Y B and D ( d x )
X φ ,X β X B to all the rows and to all the columns of the matrix F B , respecti vel y. When isolating the boundary values as in Section 2.2, in place of eq. ( 35), we obtain the following expression for computing the partial deri v ati ves of f ( x , y )

F ( d x ,d y ) = d x -1 p x = 0 d y -1 p y = 0 ( b -( p x ) x ⊗ b -( p y ) y )f --( p x , p y ) + d x -1 p x = 0 ( b -( p x ) x ⊗ D -o ( d y ) y ) f -o ( p x , 0) + d x -1 p x = 0 d y -1 p y = 0 ( b -( p x ) x ⊗ b + ( p y ) y )f -+ ( p x , p y ) (39a) + d y -1 p y = 0 ( D o-( d x ) x ⊗ b -( p y ) y ) f o-(0 , p y ) + ( D o o ( d x ) x ⊗ D o o ( d y ) y ) F oo + d y -1 p y = 0 ( D o+ ( d x ) x ⊗ b + ( p y ) y ) f o+ ( p y , 0) (39b) + d x -1 p x = 0 d y -1 p y = 0 ( b + ( p x ) x ⊗ b -( p y ) y )f +-( p x , p y ) + d x -1 p x = 0 ( b + ( p x ) x ⊗ D + o ( d y ) y ) f + o (0 , p x ) + d x -1 p x = 0 d y -1 p y = 0 ( b + ( p x ) x ⊗ b + ( p y ) y )f ++ ( p x , p y ) , ( 39c
)
where F oo = F B and the superscripts ( --, -o , -+ , o -, oo , o + , + -, + o , + + ) identify the domains from which the values are computed, as illustrated in Fig. 2 . The 1-D operators involved in the computations are

b ± ( p x ) x = b ± ( p x ) X φ ,X β ,γ x ± (40a) b ± ( p y ) y = b ± ( p y ) Y φ ,Y β ,γ y ± (40b) D ab ( d x ) x = D ± ( p x ) X φ ,X β ,X ab B ,γ x -,γ x + (40c) D ab ( d y ) y = D ± ( p y ) Y φ ,Y β ,Y ab B ,γ y -,γ y + . ( 40d
)
The external values of f ( d x ,d y ) ( x , y ) at the bottom-left ( --), bottom-right ( +-), top-left ( -+ ) and top-right ( ++ ) corners of domain are defined as

f ( d x ,d y ) ( x -,y -) = f --( d x ,d y ) f ( d x ,d y ) ( x -,y + ) = f -+ ( d x ,d y ) (41a) f ( d x ,d y ) ( x + ,y -) = f +-( d x ,d y ) f ( d x ,d y ) ( x + ,y + ) = f ++ ( d x ,d y ) . ( 41b
)
The values of

f ( d x ,d y ) ( x ,
f ( d x , 0) ( x -,y ) = N -o Y j= 1 f -o ( d x , 0) j Y -o j ( y ) f ( d x , 0) ( x + ,y ) = N + o Y j= 1 f + o ( d x , 0) j Y + o j ( y ) (42a) f (0 ,d y ) ( x ,y -) = N o- X j= 1 f o-(0 ,d y ) j X o- j ( x ) f (0 ,d x ) ( x ,y + ) = N o+ X j= 1 f o+ (0 ,d x ) j X o+ j ( x ) . ( 42b
)
Note that the bases representing the external values do not need to be the same as those representing f ( x , y ) inside domain which is useful for non-conformal multidomain modelling. In higher dimension the DFD operators can be obtained directly by taking the tensor product of the 1-D operators, for example, in three dimension, we have

D ( d x ,d y ,d z ) φ,β,B = D ( d x ) X φ ,X β X B ⊗ D ( d y ) Y φ ,Y β ,Y B ⊗ D ( d z ) Z φ ,Z β Z B . ( 43
)
where the 1-D operators D ( d)

A,B,C are given by eq. ( 12).

Accounting for material heterogeneity using weighted DFD operators

This section introduces an approximation to account for heterogeneity in the medium when modelling wave propagation without modifying the DFD operators introduced in the previous sections. To solve the wave equation numerically in heterogeneous media, we need to compute products of the form

p( x ) = h ( x ) f ( x ) (44)
between some spatially varying physical property h ( x ) (e.g. density, elastic moduli, etc.) and some field variable f ( x ) (e.g. velocities, displacements, stresses, etc.). Let p ( x ), h ( x ) and f ( x ) be represented as the linear combinations

f ( x ) = N ˆ B j= 1 ˆ f ˆ B j ˆ B j ( x ) (45a) p( x ) = N ˆ B j= 1 ˆ p ˆ B j ˆ B j ( x ) (45b) h ( x ) = N H j= 1 h j H j ( x ) , ( 45c
)
where p ( x ) and f ( x ) are represented in the same orthogonal basis ˆ B (i.e. as defined in eq. 20) with expansion coefficients ˆ p ˆ B j and ˆ f ˆ B j stored in vectors ˆ p ˆ B and ˆ f ˆ B , respecti vel y. h ( x ) is represented in an arbitrary basis H with expansion coefficients h j stored in vector h . Our objective is to determine the expansion coefficient vector ˆ p ˆ B that gives an accurate approximation of the product in eq. ( 44). To do so, we proceed as we did for deriving the DFD operators. By multiplying both sides of eq. ( 44) by the same test function

φ i ( x ) = ˆ B i ( x ) ( i = 1 , N ˆ B
) and inte grating ov er the domain of interest we obtain the linear system

M ˆ p ˆ B = I ˆ p ˆ B = ˆ p ˆ B = M h f ˆ B , (46)
where

M = I = M i, j = ˆ B i , ˆ B j x + x - = δ i j (47a) M h = ( M h ) i, j = h ˆ B i , ˆ B j x + x - = ˆ B i , h ˆ B j x + x - . ( 47b
)
Note that an obvious advantage when working with the orthogonal basis in eq. ( 20) is that the mass matrix in eq. ( 47a) is the identity matrix.

In this case, the expansion coefficient vector p ˆ B is obtained b y multipl ying the discrete field variable f ˆ B by the weighted mass matrix M h . This matrix may be pre-computed and stored prior to running a simulation. Though this is not a problem in 1-D, storing the matrix M h requires lots of memory, and it dramatically increases the computational burden in 2-D and 3-D. To overcome this issue, one can substitute the mass matrix M h by the diagonal matrix M h obtained by summing all line coefficients of M h onto the diagonal. We have:

M h = diag ( h ) , ( 48
)
where the vector

h = W h (49)
can be obtained directly from the vector h representing the physical properties using the matrix

W = W ˆ B ,H = W ˆ B ,H i j = x + x - ˆ B i ( x ) H j ( x ) N ˆ B k= 1 ˆ B k ( x )dx (50)
that is independent of the physical properties. It follows that the vector ˆ p ˆ B j can be approximated using

ˆ p ˆ B ≈ h • ˆ f ˆ B , (51)
where • denotes the element-wise or Hadamard product. The approximation M h ≈ M h is a common practice in finite-element modelling called mass lumping. Note, ho wever , that performing the mass lumping on the matrix M h constructed from the orthogonal basis ˆ B guarantees that one obtains an exact result when the weight function h ( x ) is constant. This is not true when the mass lumping is applied to the regular mass matrices M 11 M 22 in eq. ( 19), for example. That aspect is key in the proposed approach and gives excellent approximations when the physical properties are varying smoothly and the basis functions in ˆ B i are spatially localized. As for the DFD operator, we can extend this approach to higher dimension using the tensor product. For example, the 2-D counterpart of the vector in eq. ( 49) is

H = W x ⊗ W y H , ( 52
)
where the matrix H contains the expansion coefficients of h ( x , y ) represented in basis H ( x , y ) = X H ( x ) Y H ( y ), and W x and W y can be computed using eq. ( 50) with the substitutions ˆ

B i = ˆ X i , ˆ H j = ˆ X H j and ˆ B i = ˆ Y i , ˆ H j = ˆ Y H j ,
assuming that ˆ X i and ˆ Y i are used to represent the field variables.

D F D M O D E L L I N G I N T H E T I M E D O M A I N

1-D w av e pr opagation modelling and stability analysis

This section presents a 1-D DFD scheme for modelling wave propagation in the time domain. It is shown that when the adjoint DFD operators introduced in Section 2.3 are used to discretize the wave equation, we obtain a numerical scheme that is stable under the classic Courant (CFL) condition.

In one dimension the wave equation can be written using the displacement-stress formulation

ρ ü = ∂ ∂ x σ (53a) σ = k ∂ ∂ x u. ( 53b
)
where u is the displacement, ü is the acceleration, σ is the 1-D stress, ρ is the density and k is an elastic constant. Eq. ( 53) can be discretized on a single domain by representing the displacement u and the stress σ in the orthogonal bases ˆ B 1 and ˆ B 2 , respecti vel y, and, the density ρ and the elastic constant k in an arbitrary basis H . In this case, the discrete wave equation is obtained directly by substituting the spatial partial deri v ati ves with the operators in eqs ( 23a ) and ( 23b ) computed using d = 1. Remember that those operators encompass the boundary conditions which are specified by the coefficients α ± . The time derivative can be approximated using the explicit second order finite-difference scheme and the material properties can be replaced by the matrices [ ρ 1 ] and [ k 2 ] computed using eq. ( 49). We obtain the following numerical scheme for solving the wave equation in the time domain

[ ρ 1 ] ˆ u n + 1 -2 ˆ u n + ˆ u n -1 t = -ˆ D T 21 ˆ σ n (54a) ˆ σ n = [ k 2 ] ˆ D 21 u n . (54b)
At each time step, the stress is computed first using eq. ( 54b). Then, knowing the displacement at the time steps n and n -1, eq. ( 54a) is used to obtain the displacement at time step n + 1. When the operators

ˆ D 21 = L -1 2 D * (1) ( L -1 1 ) T and ˆ D T 21 = L -1 1 ( D * (1) ) T ( L -1
2 ) T are computed using the bases presented in the Appendix, the matrix D * (1) is band diagonal and the matrices L 2 and L 2 are lower triangular band diagonal matrices. The computational cost associated with the operator D * (1) is the same as that of a standard finite-difference matrix (i.e. when the bandwidths of the matrices are the same). The operators L -1 1 , L -1 2 , ( L -1 1 ) T and ( L -1 2 ) T can be ef ficientl y implemented using forward and backward substitution algorithms. Their computational cost is about half that of the operator D * (1) . Thus the total computational cost of the algorithm in eq. ( 54) is about twice that of a comparable finite-difference algorithm (when the bandwidth of the FD matrix is the same as that of D * (1) ).

To discuss the stability of the numerical scheme, it is convenient to write it using the compact form

M u n + 1 -2 u n + u n -1 t = -Ku n , ( 55
)
where

M = [ ρ 1 ] (56a) K = ˆ D T 21 [ k 2 ] ˆ D 21 , ( 56b
)
that is the general form of the discrete wave equation [START_REF] Fichtner | Full Seismic Waveform Modelling and Inversion[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]. It is known that the time domain numerical scheme in eq. ( 55) is conditionally stable if M is a diagonal positive definite matrix and K is a symmetric positive definite matrix [START_REF] Samarskii | The Theory of Difference Schemes[END_REF][START_REF] Mattsson | Summation by parts operators for finite difference approximations of second derivatives[END_REF][START_REF] Do Vgilo Vich | High-accuracy finite-difference schemes for solving elastodynamic problems in curvilinear coordinates within multiblock approach[END_REF]. In this case, the numerical scheme is stable when it satisfies an analogous to the classic Courant-Friedrichs-Lewy (CFL) convergence condition

t ≤ CFL t ≈ C min x v , ( 57
)
where CFL t denotes the largest time step for which the scheme is stable, C is a small constant, x is the grid spacing and v is the wave velocity. Because the DFD operators are adjoint, the matrix K has the general form A T BA where B = [ k 2 ] is a diagonal matrix by construction. Therefore, K is positive semidefinite if B is positive. The matrix M = [ ρ 1 ] is diagonal by construction. Consequently, the numerical scheme in eq. ( 54) is conditionally stable if [ ρ 1 ] and [ k 2 ] are positive matrices. This is usually the case for realistic smoothly varying distributions, further, the positiveness of these matrices can be easily verified prior running a simulation.

So far, we have seen that the adjoint DFD operators in Section 2.3 directly lead to a stable numerical scheme when a single domain is used. When modelling wave propagation using multiple domains, the boundary terms (i.e. the sum terms in eq. 21) need to be accounted for when computing the partial deri v ati ves of the wavefield. In this case, the boundary coefficients α ± in eq. ( 24) need to be set consistently to ensure that the total DFD operators (that e v aluate the deri v ati ve in the whole computational domain) stay adjoint. For example, if the computational domain is divided in N intervals i = [ x i , x i+ 1 ] ( i = 0, N ), a continuity condition of the following type needs to be enforced

α -[ i] = 1 -α + [ i-1] = α i for all i = 1 , N . ( 58
)
Here, the subscript [ i ] is used to indicate the domain i in which the parameters α ± are defined. From now on, a similar notation is used for the other variables. The condition in eq. ( 58) translates the fact that the values of wavefield are assumed to be uniquely defined at the boundaries x i between the different domains, that is

u ( x i ) = u ( x + ) [ i-1] = u ( x -) [ i] = α i u + [ i-1] + (1 -α i ) u - [ i] (59a) σ ( x i ) = σ ( x + ) [ i-1] = σ ( x -) [ i] = (1 -α i ) σ + [ i-1] + α i σ - [ i] , (59b)
where

u ± [ i] = ( q ± 1 ) T [ i] ( L -1 1 ) T [ i] ˆ u [ i] (59c) σ ± [ i] = ( q ± 2 ) T [ i] ( L -1 2 ) T [ i] ˆ σ [ i] (59d)
are the boundary values computed in the domain i from the expansion coefficients vector ˆ u [ i] and ˆ σ [ i] . It is important to note that the continuity conditions u

( x + ) [ i-1] = u ( x -) [ i] and σ ( x + ) [ i-1] = σ ( x -) [ i]
are not satisfied by the numerical solution that is discontinuous between two consecutive elements. In general, we have u

+ [ i-1] = u - [ i] and σ + [ i-1] = σ - [ i]
. Ho wever , eq. ( 59) should be used when evaluating the wavefield at the interfaces between the elements. When e v aluating the partial deri v ati ves in eq. ( 53) using eq. ( 21) with the substitutions

ˆ f 1 → ˆ u , ˆ g 2 → ˆ
σ , f ± → u ± and g ± → σ ± , the numerical scheme in eq. ( 54) becomes

[ ρ 1 ] [ i] ˆ u n + 1 [ i] -2 ˆ u n [ i] + ˆ u n -1 [ i] t = ˆ D T 21 [ i] ˆ σ n [ i] -( B + [ i-1] ) T ˆ σ n [ i-1] + ( B - [ i+ 1] ) T ˆ σ n [ i+ 1] (60a) ˆ σ [ i] = [ k 2 ] [ i] ˆ D 21 [ i] ˆ u n [ i] -B - [ i] ˆ u n [ i-1] + B + [ i] ˆ u n [ i+ 1] , (60b)
where

B - [ i] = α i L -1 2 [ i] q - 2 [ i] ( q + 1 ) T [ i-1] ( L -1 1 ) T [ i-1] (60c) B + [ i] = (1 -α i+ 1 ) L -1 2 [ i] q + 2 [ i] ( q - 1 ) T [ i+ 1] ( L -1 1 ) T [ i+ 1] . ( 60d
)
In each domain i , the wavefield can be updated by first evaluating eq. ( 60b) knowing the displacement in the neighbouring domains, and then e v aluating eq. ( 60a) knowing the stress in the neighbouring domains. Note that the continuity condition in eq. ( 58) ensures that the boundary operators B ±

[ i] in eq. ( 60a) are the transpose of those in eq. ( 60b). Practically, to compute the stress we take the following steps in all domains:

(i) Solve L T 1 [ i] a [ i] = ˆ u n [ i] . (ii) Compute b [ i] = D * (1) [ i] a [ i] .
(iii) Update the first and last elements of b using b

1 [ i] ⇐ b 1 [ i] + α i b N [ i-1] b N [ i] ⇐ b N [ i] + (1 -α i+ 1 ) b 1 [ i+ 1] . (iv) Solve L 2 [ i] c [ i] = b [ i] . (v) Compute ˆ σ n [ i] = [ k 2 ] [ i] c [ i] .
Similarly, to update the displacement we take the following steps in all domains:

(i) Solve L T 2 [ i] a [ i] = ˆ σ n [ i] . (ii) Compute b [ i] = D * (1) [ i] T a [ i] .
(iii) Update the first and last elements of b using b

1 [ i] ⇐ b 1 [ i] + (1 -α i ) b N [ i-1] b N [ i] ⇐ b N [ i] + α i+ 1 b 1 [ i+ 1] . (iv) Solve L 1 [ i] c = b [ i] . (v) Compute ˆ u n + 1 [ i] = 2 ˆ u n [ i] -ˆ u n -1 [ i] + t [ ρ 1 ] -1 [ i] c [ i] .
By assembling the linear system in eq. ( 60), we obtain the update equation for the complete computational domain. We have:

M g ˆ u n + 1 g -2 ˆ u n g + ˆ u n -1 g t = -D T g K g D g ˆ u g , ( 61a
)
where the aggregated displacement vector is

ˆ u g = ˆ u [0] . . . . . . . . . ˆ u [ N ] T , ( 61b
)
the global DFD operator is and the matrices that contain the distributions of physical properties in the different domains are

We see that the global DFD operators D g and D T g are still adjoint when the continuity condition in eq. ( 58) is satisfied. Therefore, the multidomain scheme in eq. ( 61a) is also stable under the classic Courant condition in eq. ( 57). Practically, the linear system in eq. ( 61a) does not need to be assembled, and eq. ( 60) should be used to update the wavefield in the different domains. From a computational perspective, in eq. ( 60), it is better to factorize the multiplications by the matrices L 1 / 2 and L 1 / 2 T as detailed in the next section. The computational cost of the boundary operators is negligible as they simply communicate the first and last expansion coefficients between neighbouring elements (when the operators are constructed using the bases described in the Appendix). Thus the computational cost of the multidomain algorithm in eq. ( 60) is marginally higher than that of the single-domain algorithm in eq. ( 54) for a given problem dimension. The 1-D algorithm is extended to 2-D in the next section.

2-D w av e pr opagation modelling

The governing equations for modelling wave propagation in a 2-D linear isotropic elastic media can be written using the displacement-stress formulation

ρ ü x = ∂ x σ xx + ∂ y σ xy + f x (62a) ρ ü y = ∂ x σ xy + ∂ y σ yy + f y (62b) σ xx = ( λ + 2 μ) ∂ x u x + λ∂ y u y (62c) σ yy = λ∂ x u x + ( λ + 2 μ) ∂ y u y (62d) σ xy = μ( ∂ x u y + ∂ y u x ) , ( 62e
)
where u i is the displacement in direction i , ü i is the acceleration in direction i , σ ij are the stresses, ρ is the density, λ and μ are the Lam é parameters and ∂ i denotes the partial deri v ati ve with respect to the spatial direction i . Using the smooth coordinate transformation

x = x ( x , y ) y = y ( x , y ) (63)
illustrated in Fig. 2 , we can write eq. ( 62) in the conserv ati ve form [START_REF] Thompson | Numerical Grid Generation[END_REF]

; Appel ö & Petersson 2009 ; Do vgilo vich & Sofrono v 2015 ) | J | ρ 0 0 ρ ü x ü y = ⎡ ⎢ ⎣ ∂ ∂ x ∂ ∂ y 0 0 0 0 ∂ ∂ x ∂ ∂ y ⎤ ⎥ ⎦ | J | ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂ x ∂ x 0 ∂ x ∂ y ∂ y ∂ x 0 ∂ y ∂ y 0 ∂ x ∂ y ∂ x ∂ x 0 ∂ y ∂ y ∂ y ∂ x ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎣ σ xx σ yy σ xy ⎤ ⎥ ⎦ + | J | f x f y (64a) ⎡ ⎢ ⎣ σ xx σ yy σ xy ⎤ ⎥ ⎦ = ⎡ ⎢ ⎣ λ + 2 μ λ 0 λ λ+ 2 μ 0 0 0 μ ⎤ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂ x ∂ x ∂ y ∂ x 0 0 0 0 ∂ x ∂ y ∂ y ∂ y ∂ x ∂ y ∂ y ∂ y ∂ x ∂ x ∂ y ∂ x ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂ ∂ x 0 ∂ ∂ y 0 0 ∂ ∂ x 0 ∂ ∂ y ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ u x u y , ( 64b
)
where | J | denotes the determinant of the Jacobian matrix. We have:

J ( x , y ) = ⎡ ⎢ ⎣ ∂ x ∂ x ∂ x ∂ y ∂ y ∂ x ∂ y ∂ y ⎤ ⎥ ⎦ J -1 ( x , y ) = ⎡ ⎢ ⎢ ⎣ ∂ x ∂ x ∂ x ∂ y ∂ y ∂ x ∂ y ∂ y ⎤ ⎥ ⎥ ⎦ | J | = ∂ x ∂ x ∂ y ∂ y - ∂ y ∂ x ∂ x ∂ y (65a)
In order to discretize eq. ( 64), within each domain (e.g. domain in Fig. 2 ) the fields variables, the physical properties and the geometric variables are first represented as linear combinations similar to that in eq. ( 28) using the 2-D basis functions

B 1 i, j ( x , y ) = X 1 i ( x ) Y 1 j ( y ) 1 ≤ i ≤ N X 1 1 ≤ j ≤ N Y 1 (66a) B 2 i, j ( x , y ) = X 2 i ( x ) Y 2 j ( y ) 1 ≤ i ≤ N X 2 1 ≤ j ≤ N Y 2 (66b) H i, j ( x , y ) = X H i ( x ) Y H j ( y ) 1 ≤ i ≤ N X K 1 ≤ j ≤ N Y K , (66c)
where X 1 , X 2 , X H and, Y 1 , Y 2 , Y H are 1-D bases defined on domain [ x -, x + ] and [ y -, y + ], respecti vel y. The displacements u x , u y are represented in basis B 1 with their expansions coefficients stored in the matrices U x and U y , respecti vel y, and, the stresses σ xx , σ yy and σ xy are represented in basis B 2 with their expansions coefficients in the matrices σ xx , σ yy and σ xy , respecti vel y. We assume that all other variables are represented in some arbitrary basis H that may be different for each variable. By computing and factorizing the 1-D mass matrices [i.e. using eq. ( 9a) followed by a Cholesky factorization]

M x 1 = M X 1 ,X 1 = L x 1 L T x 1 M x 2 = M X 2 ,X 2 = L x 2 L T x 2 M x 1 = M Y 1 ,Y 1 = L y 1 L T y 1 M y 2 = M Y 2 ,Y 2 = L y 2 L T y 2 (67a)
we can construct the orthogonal bases

ˆ B 1 i, j ( x , y ) = ˆ X 1 i ( x ) ˆ Y 1 j ( y ) ˆ X 1 i ( x ) = N X 1 j= 1 ( L -1 x 1 ) i j X 1 j ( x ) ˆ Y 1 i ( y ) = N Y 1 j= 1 ( L -1 y 1 ) i j Y 1 j ( y ) (68a) ˆ B 2 i, j ( x , y ) = ˆ X 2 i ( x ) ˆ Y 2 j ( y ) ˆ X 2 i ( x ) = N X 2 j= 1 ( L -1 x 2 ) i j X 2 j ( x ) ˆ Y 2 i ( y ) = N Y 2 j= 1 ( L -1 y 2 ) i j Y 2 j ( y ) (68b)
and the dual bases

B * 1 i, j ( x , y ) = X * 1 i ( x ) Y * 1 j ( y ) X * 1 i ( x ) = N X 1 j= 1 ( M -1 x 1 ) i j X 1 j ( x ) Y * 1 i ( y ) = N Y 1 j= 1 ( M -1 y 1 ) i j Y 1 j ( y ) (69a) B * 2 i, j ( x , y ) = X * 2 i ( x ) Y * 2 j ( y ) X * 2 i ( x ) = N X 2 j= 1 ( M -1 x 2 ) i j X 2 j ( x ) Y * 2 i ( y ) = N Y 2 j= 1 ( M -1 y 2 ) i j Y 2 j ( y ) . ( 69b
)
Using these bases, the displacements u x , u y may be equally represented in the bases By approximating the time deri v ati ves ü i using a second order centred finite-difference, and, by substituting the spatial partial derivatives in eq. ( 64) with the discrete operator in eq. ( 39), where the 1-D DFD operators are computed using eq. ( 23), we obtain the discrete linear system

ˆ σ xx = λ + 2 μ 2 • ∂ x ∂ x 2 • ˆ U n (1 , 0) x + ∂ y ∂ x 2 • ˆ U n (0 , 1) x + λ 2 • ∂ x ∂ y 2 • ˆ U n (1 , 0) y + ∂ y ∂ y 2 • ˆ U n (0 , 1) y (70a) ˆ σ yy = λ 2 • ∂ x ∂ x 2 • ˆ U n (1 , 0) x + ∂ y ∂ x 2 • ˆ U n (0 , 1) x + λ + 2 μ 2 • ∂ x ∂ y 2 • ˆ U n (1 , 0) y + ∂ y ∂ y 2 • ˆ U n (0 , 1) y (70b) ˆ σ xy = μ 2 • ∂ x ∂ y 2 • ˆ U n (1 , 0) x + ∂ y ∂ y 2 • ˆ U n (0 , 1) x + μ 2 • ∂ x ∂ x 2 • ˆ U n (1 , 0) y + ∂ y ∂ x 2 • ˆ U n (0 , 1) y (70c) ˆ S xx = | J | 2 • ∂ x ∂ x 2 • ˆ σ xx + ∂ x ∂ y 2 • ˆ σ xy (71a) ˆ S xy = | J | 2 • ∂ y ∂ x 2 • ˆ σ xx + ∂ y ∂ y 2 • ˆ σ xy (71b) ˆ S yx = | J | 2 • ∂ x ∂ y 2 • ˆ σ yy + ∂ x ∂ x 2 • ˆ σ xy (71c) ˆ S yy = | J | 2 • ∂ y ∂ y 2 • ˆ σ yy + ∂ y ∂ x 2 • ˆ σ xy (71d) ˆ U n + 1 x = 2 ˆ U n x -ˆ U n -1 x + 2 t ρ -1 1 • | J | -1 1 • ˆ S (1 , 0) xx + ˆ S (0 , 1) xy + F x (72a) ˆ U n + 1 y = 2 ˆ U n y -ˆ U n -1 y + 2 t ρ -1 1 • | J | -1 1 • ˆ S (1 , 0) yx + ˆ S (0 , 1) yy + F y , ( 72b
)
where 1 and 2 denote the averaged properties of the quantities within the brackets computed using eq. ( 52), with ˆ B = ˆ B 1 and ˆ B = ˆ B 2 , respecti vel y.

The solution to the discrete linear system is obtained recursi vel y. At each time step, knowing the displacements ˆ U n -1 and ˆ U n at times ( n -1) t and n t , the displacement ˆ U n + 1 at time ( n + 1) t is computed by evaluating eqs ( 70), ( 71) and ( 72) successi vel y. In that sequence, most of the computational burden lies in the computation of the partial deri v ati ves ˆ U n (0 , 1)

x

, ˆ U n (1 , 0) x , ˆ U n (0 , 1) y , ˆ U n (1 , 0) y , ˆ S (1 , 0) xx , ˆ S (0 , 1) xy , ˆ S (1 , 0)
yx and ˆ S (0 , 1) yy . These partial deri v ati ves are obtained by applying two types of operators to the field variables:

(i) Operators of the form ( A x ⊗ I ) , ( A x ⊗ I ) T = A x T ⊗ I , I ⊗ A y and I ⊗ A y T = I ⊗ A y T which apply either a matrix A x (or its transpose) to all columns of the matrix variables or a matrix A y (or its transpose) to all rows of the matrix variables. When the matrices A x and A y are band diagonal matrices (as in our examples), these operators can be thought as classic finite-difference operators. They perform similar operations and have the same computational cost when the bandwidth of the matrices is the same as the length of the finite-difference stencil. These operators are constructed from the matrices (vectors are taken as single column matrices)

q ± x 1 = q ± (1) X 1 q ± x 2 = q ± (1) X 2 q ± y 1 = q ± (1) Y 1 q ± y 2 = q ± (1) Y 2 (73a) b * ± x 1 = ±(1 -α x ± ) q ± x 1 b * ± x 2 = ±α x ± q ± x 2 b * ± y 1 = ±(1 -α y ± ) q ± y 1 b * ± y 2 = ±α y ± q ± y 2 (73b) D * (0) x = K (0) X 2 ,X 1 D * o-(0) x 1 = K o-(0) X 1 ,X o- 2 D * o-(0) x 2 = K o-(0) X 2 ,X o- 1 (73c) D * (0) y = K (0) Y 2 ,Y 1 D * o+ (0) x 1 = K o+ (0) X 1 ,X o+ 2 D * o+ (0) x 2 = K o+ (0) X 2 ,X o+ 1 (73d) D * (1) x = K (1) X 2 ,X 1 -α x -Q -(1) X 2 ,X 1 -α x + Q + (1) X 2 ,X 1 D * -o (0) y 1 = K -o (0) Y 1 ,Y -o 2 D * -o (0) y 2 = K -o (0) Y 2 ,Y -o 1 (73e) D * (1) y = K (1) Y 2 ,Y 1 -α y -Q -(1) Y 2 ,Y 1 -α y + Q + (1) Y 2 ,Y 1 D * + o (0) y 1 = K + o (0) Y 1 ,Y + o 2 D * + o (0) y 2 = K + o (0) Y 2 ,Y + o 1 (73f)
that are computed using eq. ( 9) and stored prior running the simulation.

(ii) Operators of the form

L -1 x ⊗ I , L -1 x ⊗ I T = ( L -1 x ) T ⊗ I , I ⊗ L -1 y and I ⊗ L -1 y T = I ⊗ ( L -1 y )
T which apply either an inverse matrix L -1

x (or its transpose) to all columns of the matrix variables or an inverse matrix L -1 y (or its transpose) to all rows of the matrix variables. These operators can be thought as recursive pre-and post-filters that are applied prior and after differentiating the field variables using the operators belonging to the first class abov e. The y can be implemented ef ficientl y using forward and backward substitution and their computational cost is half that of the operators of the first class because they involve upper and lower triangular band diagonal matrices. These operators are constructed from the matrices in eq. ( 67) which are computed and stored prior running the simulation.

Practically, the partial derivatives ˆ

U n (1 , 0) x , ˆ U n (0 , 1) x , ˆ U n (1 , 0)
y and ˆ U n (0 , 1) y in eq. ( 70) are computed as follows. The displacements ˆ U n x and ˆ U n y represented in basis ˆ B 1 are first transformed using 1-D backward substitutions to obtain the displacements U n x and U n y in basis B 1 . We have:

U n x = L -1 x 1 ⊗ I T I ⊗ L -1 y 1 T ˆ U n x U n y = L -1 x 1 ⊗ I T I ⊗ L -1 y 1 T ˆ U n y . ( 74a
)
At this point, because the wavefield is represented in the basis B 1 that has compactly supported basis functions, it can be easil y e v aluated at arbitrary receivers locations. Then, the boundary vector variables

u -o x = q - x 1 T ⊗ I U n x u -o y = q - x 1 T ⊗ I U n y (75a) u + o x = q + x 1 T ⊗ I U n x u + o y = q + x 1 T ⊗ I U n y (75b) u o- x = I ⊗ q - y 1 T U n x u o- y = I ⊗ q - y 1 T U n y (75c) u o+ x = I ⊗ q + y 1 T U n x u o+ y = I ⊗ q + y 1 T U n y (75d)
are computed and sent to/exchanged with the neighbouring domains/elements (in the case of multidomain modelling). Note that, when the inner basis functions vanish at points x -and x + (i.e. when X and1 ), the equations above simply pick up the first/last row/column of the matrices U n . As the boundary vector variables become available, the partial derivatives can be computed in the dual basis B * 2 using the finite-difference-like scheme

1 i ( x -) = X 1 i ( x + ) = 0 for i = 2 , N X 1 -1 and Y 1 i ( x -) = Y 1 i ( x + ) = 0 for i = 2 , N Y 1 -1, as in Figs A1
U * n (1 , 0) x = D * (1) x ⊗ I I ⊗ D * (0) y U n x + b * - x 2 ⊗ I I ⊗ D * -o (0) y 2 u -o x + b * + x 2 ⊗ I I ⊗ D * + o (0) y 2 u + o x (76a) U * n (0 , 1) x = D * (0) x ⊗ I I ⊗ D * (1) y U n x + D * o-(0) x 2 ⊗ I I ⊗ b * - y 2 u o- x + D * o+ (0) x 2 ⊗ I I ⊗ b * + y 2 u o+ x (76b) U * n (1 , 0) y = D * (1) x ⊗ I I ⊗ D * (0) y U n y + b * - x 2 ⊗ I I ⊗ D * -o (0) y 2 u -o y + b * + x 2 ⊗ I I ⊗ D * + o (0) y 2 u + o y (76c) U * n (0 , 1) y = D * (0) x ⊗ I I ⊗ D * (1) y U n y + D * o-(0)
x 2

⊗ I I ⊗ b * - y 2 u o- y + D * o+ (0) x 2 ⊗ I I ⊗ b * + y 2 u o+ y . ( 76d
)
Note that one may start to compute the first term to the right-hand side in the above equation prior receiving the boundary vector variables. When using conformal basis functions (i.e. that are the same in all domains/elements), we have D

* (0) x = D * o-(0) x = D * o+ (0) x and D * (0) y = D * -o (0) y = D * + o (0)
y . Thus, if the inner basis functions vanish at the boundaries [i.e. when

X 2 i ( x -) = X 2 i ( x + ) = 0 for i = 2 , N X 2 -1 and Y 2 i ( x -) = Y 2 i ( x + ) = 0 for i = 2 , N Y 2 -1],
the operators involving the vectors b * simply add the vector u to the first/last column/row of the result matrix. Finally, the partial deri v ati ves are transformed using 1-D forw ard substitutions to obtain their representation in basis ˆ B 2 . We have:

ˆ U n (1 , 0) x = L -1 x 2 ⊗ I I ⊗ L -1 y 2 U * n (1 , 0) x ˆ U n (1 , 0) y = L -1 x 2 ⊗ I I ⊗ L -1 y 2 U * n (1 , 0) y (77a) ˆ U n (0 , 1) x = L -1 x 2 ⊗ I I ⊗ L -1 y 2 U * n (0 , 1) x ˆ U n (0 , 1) y = L -1 x 2 ⊗ I I ⊗ L -1 y 2 U * n (0 , 1) y . ( 77b
)
The partial deri v ati ves ˆ S (1 , 0) xx , ˆ S (0 , 1) xy , ˆ S (1 , 0) yx and ˆ S (0 , 1) yy in eq. ( 72) are computed in a similar manner. The variables ˆ S xx , ˆ S xy , ˆ S yx and ˆ S yy represented in basis ˆ B 2 are first transformed using 1-D backward substitutions to obtain their representation in basis B 2 . We have

S xx = L -1 x 2 ⊗ I T I ⊗ L -1 y 2 T ˆ S xx S xy = L -1 x 2 ⊗ I T I ⊗ L -1 y 2 T ˆ S xy (78a)
S yx = L -1 x 2 ⊗ I T I ⊗ L -1 y 2 T ˆ S yx S yy = L -1 x 2 ⊗ I T I ⊗ L -1 y 2 T ˆ S yy . (78b)
Then, the boundary vector variables

s -o xx = q - x 2 T ⊗ I S xx s o- xy = I ⊗ q - y 2 T S xy s -o yx = q - x 2 T ⊗ I S yx s o- yy = I ⊗ q - y 2
T S yy (79a)

s + o xx = q - x 2 T ⊗ I S xx s o+ xy = I ⊗ q - y 2 T S xy s + o yx = q - x 2 T ⊗ I S yx s o+ yy = I ⊗ q - y 2 T S yy . ( 79b
)
are computed and sent to/exchanged with the neighbouring domains/elements. Here, we assume that the domains neighbouring domain in Fig. 2 are oriented similarly. When the neighbouring elements are rotated, the boundary vectors in eq. ( 79) need to be rotated accordingly. In turn, the partial deri v ati ves are computed in the dual basis B * 1 using

S * (1 , 0) xx = -D * (1) x ⊗ I T I ⊗ D * (0) y T S xx + b * - x 1 ⊗ I I ⊗ D * -o (0) y 1 s -o xx + b * + x 1 ⊗ I I ⊗ D * + o (0) y 1 s + o xx (80a) S * (0 , 1) xy = -D * (0) x ⊗ I T I ⊗ D * (1) y T S xy + D * o-(0)
x 1

⊗ I I ⊗ b * - y 1 s o- xy + D * o+ (0)
x 1

⊗ I I ⊗ b * + y 1 s o+ xy (80b) S * (1 , 0) yx = -D * (1) x ⊗ I T I ⊗ D * (0) y T S yx + b * - x 1 ⊗ I I ⊗ D * -o (0) y 1 s -o yx + b * + x 1 ⊗ I I ⊗ D * + o (0) y 1 s + o yx (80c) S * (0 , 1) yy = -D * (0) x ⊗ I T I ⊗ D * (1) y T S yy + D * o-(0)
x 1

⊗ I I ⊗ b * - y 1 s o- yy + D * o+ (0)
x 1

⊗ I I ⊗ b * + y 1 s o+ yy . (80d)
Finally, the partial derivatives are transformed using 1-D forward substitutions to obtain their representation in basis ˆ B 1 . We have:

ˆ S (1 , 0) xx = L -1 x 1 ⊗ I I ⊗ L -1 y 1 S * (1 , 0) xx ˆ S (0 , 1) xy = L -1 x 1 ⊗ I I ⊗ L -1 y 1 S * (0 , 1) xy (81a) ˆ S (1 , 0) yx = L -1 x 1 ⊗ I I ⊗ L -1 y 1 S * (1 , 0) yx ˆ S (0 , 1) yy = L -1 x 1 ⊗ I I ⊗ L -1 y 1 S * (0 , 1) yy . (81b)
Overall, from a computational perspective, the proposed algorithm is very similar to classic finite-difference schemes; because the main part of the computations in eqs ( 76) and ( 80) are finite-difference like operations (with dif ferent coef ficients). The pre-processing steps in eqs ( 74) and ( 78) and the post-processing steps in eqs ( 77) and ( 81) modestly add up to the computational burden and give the DFD operator an enhanced accuracy. The explicit treatment of the boundary values in eqs ( 75) and ( 79) permits to specify boundary conditions easily and to account for sharp discontinuities as well as solid-fluid interfaces, as discussed in Section 4.

N U M E R I C A L E X A M P L E S

In this section, numerical examples demonstrating the flexibility and accuracy of the DFD algorithm are presented. Domain decomposition, surface wave modelling and the solid-fluid interface are addressed in separate examples. The following comments apply to all the examples discussed after that:

All the DFD simulations use the algorithm described in Section 3.2. The displacements u x and u y are represented in the same basis B 1 defined in eq. ( 66a). The 1-D bases X 1 and Y 1 are al wa ys of the same kind. They can be either one of the left-skewed bases pictured in black in Figs A1 and 1 (a)-(c) (with a different N ). The stresses σ xx , σ yy and σ xy are represented in the same basis B 2 defined in eq. ( 66b). The 1-D bases X 2 and Y 2 are al wa ys of the same kind as X 1 and Y 1 . They can be either one of the right-skewed bases pictured in red in Figs A1 and 1 (a)-(c). The usage of the left-and right-skewed basis functions is analogous to the MacCormack-type scheme in the finite-difference method (see e.g. [START_REF] Hixon | Evaluation of a high-accuracy maccormack-type scheme using benchmark problems[END_REF][START_REF] Hixon | Compact implicit maccormack-type schemes with high accuracy[END_REF][START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF]. It ef ficientl y pre vents grid decoupling phenomena. All the physical properties λ, μ, ρ and the quantities associated with the curvilinear mapping are represented in the same basis H defined in eq. ( 66c ). That basis is a 2-D second-order B-spline basis. The parameters α x/y ± in eq. ( 73) have been set as follows:

(i) α x/y ± = 0 at all interfaces surrounding the computational domain to model a free surface. (ii) α x/y ± = 1 / 2 at all interfaces separating two domains that are both solid/elastic ( μ = 0) or fluid/acoustic ( μ = 0). (iii) α x/y ± = 0 inside the solid/elastic domain and α x/y ± = 1 inside the fluid/acoustic domain, at all solid/fluid interfaces.

Other values could be used for the solid-solid and fluid-fluid interfaces. It was observed that taking α x/y ± = 0 or 1 at the solid-fluid interfaces ef ficientl y pre vent spurious modes from appearing. As a shortcut, domains' dimensions are given as the number of gridpoints (instead of the number of basis functions or expansion coefficients in the different spatial directions). It makes the comparison with other methods easier. By extension, the average number of expansion coefficients per wavelength is given as the number of gridpoints per wavelength. Because the principal focus of the paper is introducing the DFDM algorithm, we adopt a conserv ati ve approach and present differenceseismograms when testing the numerical solution against analytical solution and other modelling methods. An in-depth accuracy analysis of DFDM should rely on more advanced techniques (e.g. [START_REF] Kristekov Á | Misfit criteria for quantitative comparison of seismograms[END_REF][START_REF] Kristekov Á | Time-frequency misfit and goodness-of-fit criteria for quantitative comparison of time signals[END_REF].

Y. Masson

We first present numerical examples showing wave propagation in one dimension. The accuracy of DFDM is tested against SEM, FDM and an analytical solution. We then illustrate the flexibility of DFDM for domain decomposition by modelling 2-D wave propagation using elements with different shapes and sizes. The last examples investigate the relative accuracy of DFDM and SEM. The results are compared with analytical solutions. We first discuss surface wa ve modelling, follow ed by the case of a solid-fluid interface with and without topography. The size of the domains can be adapted, and small elements can be used to mesh complex structures in the velocity model.

Wav e pr opagation modelling in one dimension

In Fig. 5 In all simulations, the length of the velocity model is 6 km. The wave speed is v = 3000 m s -1 and the density is ρ = 2500 kg m -3 . The time step t = 3.85 × 10 -2 s is the same for the five simulations and is small enough to minimize the error associated with time discretization. A Ricker wavelet with peak frequency f 0 = 10 Hz is generated by imposing the displacement's initial values. The grid dimension is the same in all simulations; there are 157 gridpoints or basis functions in all simulations. There are 3 points per wavelength at the estimated maximum frequency f max = 25 (Hz). Dirichlet and Neumann boundary conditions are imposed to the left and to the right of the computational domain, respecti vel y. (For the FD simulation, the boundary conditions are imposed using the image method, that is the displacement is assumed to be antisymmetric at the left boundary and symmetric at the right boundary.) The numerical is estimated by taking the difference between the computed wavefield and an analytical solution.

One first emphasizes that the discontinuous nature of the wavefield is not visible in the DFDM simulations because the jump discontinuities between neighbouring are very small.

The DFD algorithm based on the B-spline basis in Fig. 1 (c) has the smallest numerical error followed by the DFD algorithm based on the Hermite basis in Fig. 1 (b), the DFD algorithm based on the Polynomial basis in Fig. 1 (a), the SEM algorithm and the FD algorithm. This is likely due to the higher convergence rate of the B-spline basis. The error is about similar in the DFD simulations based on the Hermite and the Polynomial bases in Fig. 5 (b) and in Fig. 5 (c). Interestingly, the DFD simulation based on the Polynomial basis in Fig. 5 (c) exhibits some slight stationary noise at the interfaces between the DFD domains located around 2 and 4 km. This noise is not present in the DFD simulations based on the Hermite and the B-spline bases. This is attributed to the Polynomial basis's lower continuity and suggests that highly continuous bases such as the B-spline basis are desirable. All the DFD simulations have a smaller numerical error than the SEM and the FD simulations. This suggests that the DFD operators require fewer points per wavelength to achieve an accuracy comparable to the SEM or the FD approaches.

Wav e pr opagation modelling using different element/domain sizes

To highlight the flexibility of the DFD algorithm, three similar simulations with different domain decompositions are presented. The same elastic model and seismic source are used in all simulations while the elements's sizes vary. The results are illustrated in Fig. 6 . The velocity model is a 1 m 2 domain with P -wave velocity v p = 3400 m s -1 , S -wave velocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The source time function is a Ricker wavelet similar to that in Fig. 3 with peak frequency f 0 = 60 kHz and an estimated maximum frequency f max = 150 kHz. The time step is t = 1.18 × 10 -7 s. Fig. 6 (a) presents the first simulation. It uses a single domain, and the basis function used is the fourth-order polynomial basis in Fig. 1 (a). The grid dimension is 250 × 250 = 62 500. There are 5.6 points per P -wavelength and 3.3 points per S -wavelength. Note that this is well below the recommended value of five points per wavelength for a spectral element simulation (with polynomial order p = 4).

Fig. 6 (b) presents the second simulation. The velocity model is divided into 25 domains with random shapes. Each domain uses a different basis chosen randomly (from those in Figs A1 and1 a-c). The grid dimensions of each domain are adapted to have 5.6 points per P -wavelength and 3.3 points per S -wavelength. The grid dimension is 7629. Most interfaces are non-conformal. The number of basis functions, their type, and polynomial order can be different in neighbouring elements. More precisely, there are three conformal interfaces where the bases exactly match and 37 non-conformal interfaces (i.e. where the basis functions are different on each side of the interface). There are 5, 3, 3, 5, 4 and 5 domains using the basis in Figs A1 (a When comparing the three simulations in Fig. 6 , we observe no differences in the modelled wavefield. The DFDM accurately handles the non-conformal interfaces, mixed polynomial orders and multiple basis types. Despite the relatively small number of points per wavelength (for the S wa ve), w e observe no visible dispersion; this emphasizes the method's accuracy. The apparent advantage of the DFDM is that the domains' sizes can be adapted depending on the complexity of the velocity model. In large regions where the elastic properties are smoothly varying, one can proceed as in Fig. 6 (a), and use one or a few large domains. When needed, one can use much smaller elements as in Fig. 6 (c). In this case, we can mesh complex structures and, the DFDM resembles more to the finite/spectral element method. It is also possible to transition from large to small elements using standard mesh refinement techniques.

Surface w av e modelling on a fine grid

In this series of e xamples, wav e propagation in a 2-D homogeneous isotropic half-space model with a free surface is addressed. The results are presented in Figs 7 (a), (b) and 8 . It is important to note that the DFD algorithm allows to explicitly include the stress conditions at the interface and put the shear velocity to zero in the fluid domain. This is possible because the proposed algorithm is based on the first-order PDEs.

The computational domain has dimensions 6400 m by 3200 m, the P -wav e v elocity is v p = 3400 m s -1 , the S -wave velocity is v s = 1963 m s -1 and the density is ρ = 2500 kg m -3 . An e xplosiv e source is placed at the centre of the free surface, and the wavefield is recorded at fiv e receiv ers equally spaced at the surface. The source time function is a Ricker wavelet similar to that in Fig. 3 with peak frequency f 0 = 10 Hz and an estimated maximum frequency f max = 25 Hz. The distance from the source to the receivers ranges from 320 to 1920 m.

A DFD simulation has been run for each one of the fourth-order bases in Figs 1 (a)-(c), to compare their accuracy. In these simulations, the model is divided into 15 domains with dimensions 1280 m by 1067 m. Each domain contains 116 × 97 gridpoints. In total, the wavefield is represented using 15 × 116 × 97 = 168 780 gridpoints. All the interfaces are conformal in these simulations. the fiv e receiv ers represented by b lues circles in F ig. 7 . The seismograms obtained using the DFDM (dashed blue lines) match perfectly the analytical solution (solid black line). The seismograms obtained using the SEM (dashed red lines) exhibit some spurious ringing in the surface waves because the number of points per wavelength is insufficient. Panels (c) and (d) show the difference (multiplied by 2) between the computed seismograms and the analytical solution.

The difference in panels (c) and (d) are cut when the analytical solution predicts no more waves. The maximum difference between the analytical solution and the computed seismograms is about 8 and 25 per cent of the signals' amplitude for the DFD and SEM seismograms, respectively.

A spectral element simulation has been run using the Specfem code [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF]. This simulation divides the model into 144 × 72 elements with polynomial order 4, producing a grid with 166 753 points.

The number of points per wavelength is similar in the DFD and the SEM simulations. We have 12 points per wavelength for the P wave and 7 points per wavelength for the S wave (these estimates are based on f max = 25 Hz). For a fair comparison, the time step is fixed to t = 8.5 × 10 -4 s in all simulations.

The three DFD simulations produced nearly indistinguishable results and only the wavefield obtained using the fourth polynomial basis is presented in Fig. 7 (a). The wavefield obtained from the SEM simulation is pictured in Fig. 7 (b). We observe no visual difference between the two snapshots. The P wave, S wa ve, Ra yleigh-wa ve and Head-wave emerge as expected, as indicated in Fig. 7 (b). One notes no significant numerical dispersion. The seismograms recorded at the fiv e receiv ers for the DFD and the SEM simulation are shown in Figs 8 (a) and(b). These seismograms are compared to an analytical solution e v aluated using the Gar6more software [START_REF] Ezziani | Analytical solution for waves propagation in heterogeneous acoustic/porous media part I: the 2D case[END_REF]. This software implements the Cagniard-de Hoop method. The computed seismograms perfectly match the analytical solution because there is a sufficient number of points per wavelength. Not sur prisingly, the er ror pictured in Figs 8 (c) and (d) is increasing with distance because of the accumulation of the numerical dispersion. The Ra yleigh wa ve that has the smallest number of points per wavelength is the most affected. Conversel y, the P w ave is sampled by the largest number of gridpoints shows almost no error. The maximum error observed in these examples is about 5 per cent of the signals' amplitude for the SEM simulation and 2.5 per cent for the DFD simulations. Thus, the DFDM appears to be more accurate than SEM when using the same polynomial order and the same number of points per wavelength.

Surface w av e modelling on a coarse grid

Examples similar to those in the previous section but using a coarser grid are presented. The grid spacing has been divided by two, and the time step was adapted accordingly; we have t = 1.7 × 10 -3 s. The parameters that are different from the fine-grid case are listed in Fig. 4 . In this situation, we expect to observe a lot more numerical dispersion. The results are presented in Figs 7 (c), (d) and9 . In the DFD simulations, the 15 domains have dimensions 58 × 49 gridpoints and there is 15 × 58 × 49 = 42 630 gridpoints in total. In the SEM simulation, the model is divided into 72 × 36 elements. This produces a grid with 41 905 points. There is the same number of points per wavelength in the SEM and the DFD simulations. We have 6.1 points per P -wavelength and 3.5 points per S -wavelength.

The seismograms obtained for the three DFD simulations are almost identical, and only the results obtained with the B-spline basis are presented. The error in the computed seismograms is more prominent than in the fine grid case because there are fewer points per wavelength. The maximum error is about 25 per cent of the signals' amplitude for the SEM simulation and 8 per cent for the DFD simulations. Most importantly, the increase in the error is more considerable for the SEM simulation. A significant amount of numerical dispersion affects the surfaces wave (as indicated by the red arrows in Fig. 7 d). This numerical dispersion appears as high-frequency oscillations in the SEM seismograms (the dashed red lines in Figs 9 a and b). On the contrary, the DFD seismograms (the dashed blues lines in Figs 9 a and b) match remarkably well the analytical solution. These examples confirm that the DFDM is more accurate when using the same number of points per wavelength and polynomial order for the basis functions.

Wav e pr opag ation near a fluid-solid interf ace using a fine grid

A situation that is notoriously delicate to model is when waves propagate through media containing fluid-solid interfaces. In this case, spurious modes may arise, and the modelling algorithms often need to be adapted (see, e.g. [START_REF] Komatitsch | Wave propagation near a fluidsolid interface: a spectral-element approach[END_REF]. The case of two horizontal layers separated by a flat interface [START_REF] Komatitsch | Wave propagation near a fluidsolid interface: a spectral-element approach[END_REF] is presented to demonstrate that the DFD is accurate in this delicate situation. The results are shown in Figs 10 (a), (b) and 11 .

The velocity model has dimensions 6400 m by 3200 m. The bottom layer is a 1600 m high elastic region with P -wave velocity v p = 3400 m s -1 , S -wave velocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The top layer is a 1600 m high acoustic region with P -wave velocity v p = 1500 m s -1 and density ρ = 1020 kg m -3 . The waves are generated by an e xplosiv e source placed in the acoustic region. The source time function is a Ricker wavelet with peak frequency f 0 = 10 Hz and an estimated maximum frequency f max = 25 Hz. The wavefield is recorded using five receivers in the acoustic region and five receivers in the elastic region. The receivers in the acoustic region measure the direct, reflected and refracted/head waves. The receivers in the elastic region record the converted and transmitted waves. Two simulations have been run.

The first simulation is a DFD simulation using the fourth-order B-spline basis in Fig. 1 (a). The velocity model is divided into 12 domains with dimensions 193 × 109 = 21 037 gridpoints. The total grid has 252 444 points. The second simulation is using the SEM software Specfem [START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF], the velocity model is divided into 144 × 108 = 15 552 elements with polynomial order 4 and has 249 841 gridpoints in total.

In both simulations, the time step is t = 8.5 × 10 -4 and the number of points per wavelength is the same. There are 12.2 points per P -wavelength and 7 points per S -wavelength in the elastic layer. There are 5.4 points per P -wavelength in the acoustic layer. In the DFD simulation, the solid-fluid interface is simply handled by setting α y + = 0 in eq. ( 73) inside the three elastic domains neighbouring the interface and α y -= 1 inside the three acoustic domains. In the SEM simulation, the solid-fluid interface is handled using the method in [START_REF] Komatitsch | Wave propagation near a fluidsolid interface: a spectral-element approach[END_REF] where problem is formulated in terms of displacement in elastic regions and a velocity potential in acoustic regions.

In both simulations, we observe no numerical artefact, and no spurious modes emerge. The DFD and the SEM snapshots in Figs 10 (a) and (b) are visually identical. The corresponding seismograms in Fig. 11 perfectly agree with the analytical solution [START_REF] Ezziani | Analytical solution for waves propagation in heterogeneous acoustic/porous media part I: the 2D case[END_REF] ). The error does not exceed 5 per cent of the signals' amplitude. The error is almost similar for DFD and SEM and is more significant for the receivers placed in the acoustic region because it has fewer points per wavelength. These examples demonstrate that the DFDM accurately models wave propagation in the vicinity of a solid-fluid interface. No modification to the basic algorithm is needed in that situation.

Wav e pr opag ation near a fluid-solid interf ace using a coarse grid

The following examples are similar to those in the previous section but use a coarser grid to investigate the robustness of the DFDM. All the parameters are identical except the grid dimensions and the time step, as summarized in Fig. 4 . The results are presented in Figs 10 (c),(d) and 12 .

In the DFD simulation, the dimensions of the 15 domains are 97 × 55 = 5335 gridpoints. The complete grid has 64 020 gridpoints. In the SEM simulation, the velocity model is divided into 72 × 54 = 3888 elements. The grid contains 62 713 points. In both simulations, the time step is t = 1.7 × 10 -3 s, and the number of points per wavelength is the same. There are 6.1 points per P -wavelength and 3.5 points per S -wavelength in the elastic layer, and 2.7 points per P -wavelength in the acoustic layer.

When comparing the snapshots in Figs 10 (c) and (d) to those in Figs 10 (a) and (b), we see that the wavefield obtained using the SEM simulation exhibits a significant level of high-frequency noise in the acoustic region where the number of points per wavelength is the smallest. For the DFDM, the snapshot in Fig. 10 (c) looks almost similar to that in Fig. 10 (a). It emphasizes the superior accuracy of the DFDM when a few points per wavelength are being used. This effect is visible in the seismograms pictured in Fig. 12 . The SEM seismograms (plotted using dashed red lines) contain high-frequency oscillations which are not visible in the DFD seismograms (represented using dashed blue lines). These latter match the analytical solution relati vel y well. In Figs 12 (c) and (d), we see that the measured error is about 20 per cent of the signals' amplitude for the DFDM and goes as high as 60 per cent for SEM. Even though the DFD seismograms are also affected by a small amount of numerical dispersion, all the wave packets observed correspond to the arri v al of physical waves propagating through the model.

Wav e pr opag ation near a fluid-solid interf ace with surf ace topography

As a tribute to Dimitri Komatitsch, the last example considers wave propagation near a fluid-solid interface with surface topography and is adapted from the famous example in [START_REF] Komatitsch | Wave propagation near a fluidsolid interface: a spectral-element approach[END_REF]. This example is significant to demonstrate the ability of the DFDM to accurately model wave propagation when the domains are deformed, for example to account for curved discontinuities. It also confirms that the DFDM can accurately account for heterogeneity within the domains. Because the same approximation is used to account for the curvilinear mapping and the heterogeneity in the material properties. The results are presented in Figs 13 and 14 . The velocity model is almost similar to the one introduced in the previous section. It has dimensions 6400 m by 3200 m and consists of two horizontal layers separated by a sinusoidal interface. The sinusoidal interface is centred vertically. Its height varies between 1420 and 1780 m, and its periodicity is 1067 m. The properties of the two layers are the same as in the previous example. The source parameters are also identical except its position. The wavefield is recorded using eight receivers in the acoustic region.

A DFD simulation and a SEM simulation have been run. The DFD simulation uses the fourth-order B-spline basis in Fig. 1 (a), and the velocity model is divided into 12 domains. The domains have different dimensions in the elastic and acoustic regions. Thus, the mesh is non-conformal at the interface. In the acoustic region, the domains have dimensions 193 × 109 = 21 037 gridpoints, and in the elastic region, the domains have dimensions 147 × 83 = 12 201 gridpoints. The total grid has 199 428 points. There are 9.3 points per P -wavelength and 5.4 points per S -wavelength in the elastic layer and 5.4 points per P -wavelength in the acoustic layer. For the SEM simulation, the velocity model is divided into 144 × 108 = 15 552 elements with polynomial order 4. The grid has 249 841 points. There are 12.2 points per P -wavelength and 7 points per S -wavelength in the elastic layer and 5.4 points per P -wavelength in the acoustic layer. In both simulations the time step is t = 8.5 × 10 -4 .

The snapshots in Fig. 13 show that the presence of the sinusoidal interface produces a comple x wav efield where multiple reflected and transmitted waves emerge and pass at the receiver's locations at different times. There is no visual difference between the wavefields computed using DFDM and SEM. The seismograms recorded at the receivers are plotted in Fig. 14 and match perfectly. Since no analytical solution was readily available for an interface with topography, the signals are compared by subtracting the DFD seismograms from the SEM seismo grams. The maximum dif ference does not exceed 6 per cent of the seismograms' amplitude. In this example, both DFDM and SEM gi ve suf ficientl y accurate results. As discussed in the previous section, other tests using the same configuration but a coarser grid have been performed. These tests led to the same conclusion.

D I S C U S S I O N

Prior ending our discussion, let's outline the principal characteristics of DFDM:

(i) The computational domain is decomposed into multiple subdomains or elements similarly to what is done in FEM, SEM or DGM. In this study, rectangular/hexahedral elements are considered. Ho wever , the proposed approach could be extended to account for triangular/tetrahedral elements.

(ii) The subdomains can be arbitrarily large, and w e ma y use multiple elements depending on the velocity model. This is different from most element methods, where the element size is usually small and fixed. PSM and SBP-FDM also allow for arbitrary element size. Thus in DFDM, we can do the domain decomposition by hand when a few large elements are used, as in FDM. Using an external meshing tool is also possib le w henev er comple x structures and interfaces need to be represented accurately, as in FEM, SEM or DGM.

(iii) Within a subdomain, the structure of the proposed algorithm is relati vel y simple and similar to FDM because the DFD operators are constructed from 1-D basis functions having compact support using a tensor product. This makes DFDM efficient, because all the computations rely on 1-D operations involving band diagonal matrices.

(iv) In DFDM, the CPU timescales with the problem's dimension as in FDM and is largely independent of the element size. This is different from spectral methods such as PSM or SEM, for example.

(v) Inside the subdomains, DFDM can account for smooth and discontinuous material heterogeneity. In this work, smooth heterogeneity is incorporated by lumping the mass matrix. The latter is diagonal in DFDM when the material is homogeneous. As opposed to traditional mass lumping in FEM, this does not significantly deteriorate the accuracy of DFDM. It is also possible to adapt some FDM schemes with a subcell resolution to account for internal discontinuities. This approach is treated in a separate publication.

(vi) The free surface can be accounted for naturally in DFDM because the discrete differential operators are derived using a weighted residuals method that is analogous to taking the weak form of the wave equation. Thus we obtain an accuracy similar to FEM, SEM, DGM and SBP-FDM. Surface topography can be accommodated using boundary conforming curvilinear meshes.

(vii) In the case of multidomain modelling, the global linear system does not need to be assembled similarly to FEM, SEM and DGM. The boundary values are exchanged between neighbouring elements, allowing for efficient parallel computations.

(viii) As opposed to FEM and SEM and similarly to DGM, the numerical wavefield is discontinuous. Both the displacement and the stress can be explicitly defined at the surface of the elements. In this respect, DFDM is closer to DGM, which could make it possible to include non-linear interface conditions, arising, for example in dynamic earthquake rupture problems. This is not considered in this paper.

(ix) Nonconformal interfaces, where the basis functions are different in two neighbouring subdomains, are accounted for naturally in DFDM. Thus, the discontinuous grid is a natural feature of DFDM. This includes solid-fluid interfaces that require specific treatment in SEM, for example.

(x) The proposed DFDM algorithm relies on the first-order PDEs and uses the same equations in the fluid and the solid regions. This is similar to FDM on staggered grids. It differs from SEM which is based on the second-order PDEs and uses acoustic and elastic equations in fluid and solid domains, respecti vel y. Acoustic w ave propagation is simply achieved by putting the shear velocity to zero in the fluid domains.

(xi) In DFDM, the field variables are represented using or thonor mal bases constr ucted from basis functions having compact support. The orthogonal basis functions are not compactly supported and span the entire width of the elements. This gives the discrete differential operators in DFDM a spectral aspect that is reminiscent of PSM and ensures good accuracy.

(xii) An important aspect of DFDM is that the stress and the displacement components are represented in different bases. This gives some flexibility to annihilate non-physical, parasitic modes. In this paper, left-and right-skewed basis functions are used to mimic the MacCormack-type scheme used in FDM (see e.g. [START_REF] Hixon | Evaluation of a high-accuracy maccormack-type scheme using benchmark problems[END_REF][START_REF] Hixon | Compact implicit maccormack-type schemes with high accuracy[END_REF][START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF]. One may also use shifted basis functions analogous to the staggered grids in FDM (see e.g. [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF][START_REF] Virieux | P-sv wave propagation in heterogeneous media: velocitystress finite-difference method[END_REF]; this is addressed in a separate publication.

(xiii) In this paper, a displacement-stress formulation is used w hich av oids recording the stress versus time (see e.g. [START_REF] Luo | Parsimonious staggered grid finitedifferencing of the wave equation[END_REF]. Ho wever , the velocity-stress formulation that is necessary for simulating rupture propagation, for example, can also be used and will be detailed in a separate publication.

(xiv) The DFD operators are adjoint, and the wave equation's self-adjoint nature is preserved at the discrete level. This critical property ensures the numerical reciprocity and the stability (subject to the CFL condition) of the proposed algorithm.

The points above put DFDM in perspective with other modelling methods. It should help e v aluating the adequacy of DFDM for addressing specific geophysical problems.

C O N C L U S I O N

In summary, an explicit time-stepping DFD scheme for modelling wave propagation in heterogeneous elastic media has been introduced. The originality of this work is that the DFD differential operators are derived independently from the governing equation describing wave propagation, using a variational approach. This permits discretizing the wave equation directly from its strong form by substituting the partial deri v ati ves with the DFD differential operators. This is similar to what is done in the finite-difference method.

Adjoint pairs of DFD operators have been constructed and are needed to ensure the numerical stability of the DFD algorithm. Another important aspect of this work is the introduction of an averaging scheme to account for heterogeneity in the elastic properties of the medium. This without increasing the computational complexity of the algorithm. The computational structure of the proposed algorithm is similar to that of the finite-difference method.

Because the modelling domain can be decomposed into domains or elements having arbitrary sizes, the DFDM is very flexible and can be used like a finite-difference method by using a few large elements, or, like a finite/spectral element method where complex structures in the velocity model are meshed using small elements. The proposed formalism also allows representing the wavefield using different basis functions, grid spacing (i.e. non-conformal meshes), and polynomial order in the various domains or elements, which adds to its flexibility.

Finall y, the DFDM accuratel y models the most common situations, including free surf aces and solid-fluid interf aces, and can thus be used to address a wide variety of geophysical problems.

A C K N O W L E D G M E N T S

I would like to thank the grant I-SITE E2S-ANR-THERMapp and the France-Berkeley Fund (Project #10-2019) for their support. I wish to thank the editor Sidao Ni as well as Peng Yong and Peter Moczo for their insightful and constructive comments that helped me improve the quality of this paper. I thank Julien Diaz for providing his code and, Jean Virieux and Yann Capdeville for enthusiastic discussions about DFDM.

D ATA AVA I L A B I L I T Y

The data underlying this paper will be shared on reasonable request to the corresponding author (yder masson@berkeley.edu).

R E F E R E N C E S

Figure 1 .

 1 Figure 1. (a) The fourth order polynomial bases obtained following the approach in Section A1 (see Fig. A1 a for details). (b) The fourth order Hermite bases obtained following the approach in Section A2 (see Fig. A1 b for details). (c) The fourth order B-spline bases obtained following the approach in Section A3(c) (see Fig. A1 c for details). (d) The dual bases obtained from the polynomial bases in panel (a) using eq. (25). (e) The dual bases obtained from the Hermite bases in (b) using eq. (25). (f) The dual bases obtained from the B-spline bases in panel (c) using eq. (25). (g) The dual bases obtained from the polynomial bases in panel (a) using eq. (20). (h) The dual bases obtained from the Hermite bases in panel (b) using eq. (20). (i) The dual bases obtained from the B-spline bases in panel (c) using eq. (20).

)

  By representing the functions f ( x ) and g ( x ) in the orthogonal bases ˆ B 1 and ˆ B 2 , respectively, we can construct a pair ˆ D f ( x ) in basis ˆ B 1 and returns the vector ˆ f ( d) 2 representing the d th deri v ati ve of f ( x ) in basis ˆ B 2 , and, ˆ D ( d)

  eqs ( 17 ) and ( 18 ) using the bases substitutions B = ˆ B 1 , φ = β = ˆ B 2 with weights γ -= α -, γ + = α + , and, the operators ˆ D eqs ( 17 ) and ( 18 ) using the substitutions B = ˆ B 2 , φ = β = ˆ B 1 with complementary weights γ -= (1 -α -), γ + = (1 -α + ).

  y ) along the left ( -o ), right ( + o ), bottom ( o -) and top ( o + ) faces of domain represented in the 1-D bases Y -o , Y + o , X o -and X o + , respecti vel y, are represented as

Figure 3 .

 3 Figure3. The source-time function (Ricker wavelet) used in the numerical simulations (left-hand panel) and its Fourier transform (right-hand panel). f 0 denotes the peak frequency at which the amplitude of the Fourier transform is maximum and f max ≈ 2.5 f 0 is the estimated maximum frequency above which the amplitude of the Fourier transform becomes negligible. In this particular example f 0 = 10 Hz and f max = 25 Hz.

Figure 4 .

 4 Figure 4. Summary of the parameters that take a different value in the simulations in Figs 7 and 10 .

Figure 5 .

 5 Figure 5. 1-D snapshots (i.e. showing the displacement as a function of distance at a given time) in five different simulations (3 DFD simulations with the fourth-order bases in Figs 1 (a)-(c), 1 FD simulation, and 1 SEM simulation). The dashed red lines show the computed wavefield in the top panels, and the solid black lines show the analytical solution. The difference between the numerical and analytical solutions (multiplied by 2) is represented in the bottom panels. Panel (a) shows the result obtained for a DFD simulation using 3 subdomains having the same dimension and the B-spline basis in Fig.1 (c). Panel (b) shows the result obtained for a DFD simulation using 3 subdomains having the same dimension and the Hermite basis in Fig.1 (b). Panel (c) shows the result obtained for a DFD simulation using 3 subdomains having the same dimension and the Polynomial basis in Fig.1 (a). Panel (d) shows the result obtained for an SEM simulation using 39 elements having the same dimension and polynomial order p = 4. Panel (e) shows the result obtained for an FD simulation with 4th-order precision in space and 157 gridpoints.

Figure 6 .

 6 Figure 6. Snapshots showing the horizontal displacement in three different simulations at time t = 1.4 × 10-4 . The outermost wave is a P wave, while the innermost wave is an S wave. The velocity model is a square meter domain with P -wave velocity v p = 3400 m s -1 , S -wave velocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The source is a horizontal point force, and the red diamond marks its location. The source time function is a Ricker wavelet with peak frequency f 0 = 60 kHz. The time step is the same for the three simulations, and the grid dimension (i.e. number of degrees of freedom) is almost identical. We have about 5.6 points per P -wavelength and 3.3 points per S -wavelength. In panel (a), the simulation is performed using a single domain and resembles a finite-difference simulation. The simulation uses the fourth-order polynomial basis. Panel (b) is similar to panel (a) but the velocity model is divided in 25 domains (delimited by the solid black lines). The wavefield is represented using a different basis in each domain. Panel (c) is similar to panel (b) but the velocity model is divided in 625 domains. Note that the three snapshots are visually identical. No numerical dispersion is observed despite the small number of points per -wavelength and the variety of bases used to represent the wavefield in the different domains. These simulations illustrate the ability of the DFDM to adapt/refine the meshing depending on the complexity of the velocity model. In the large and relati vel y homo geneous region, one can use one or a few big domains. The size of the domains can be adapted, and small elements can be used to mesh complex structures in the velocity model.

  , 1-D simulations in a bounded homogeneous medium are presented to highlight the amount of numerical dispersion obtained when modelling wave propagation by implementing the DFD algorithm in Section 3.1 with the three bases in Figs 1 (a)-(c). A small number of points per wavelength is used to emphasize the dispersion of the different numerical schemes. The results correspond to three simulations based on the DFD algorithm with the bases in Figs 1 (a)-(c) plus one FD and one SEM references simulations. The FD simulation is fourth-order accurate in space and second-order accurate in time (Le v ander 1988 ).

Figure 7 .

 7 Figure 7. The snapshots show the magnitude of the displacement vector for four different simulations at time t = 0.85 s. The magnitude of the displacement is raised to the power 1/3 to enhance waves having a small amplitude. The velocity model has dimensions 6400 m by 3200 m with P -wave velocity v p = 3400 m s -1 , S -wave velocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The waves are generated by an e xplosiv e source placed at the position marked by a red diamond. The source time function is a Ricker wavelet with peak frequency f 0 = 10 Hz. The receivers are represented using blue circles. The seismograms recorded in the simulations (a, b) and (c, d) are plotted in Figs 8 and 9 , respecti vel y. The dif ferent seismic waves are indicated in panel (b). Panel (a) shows the wavefield obtained for a DFD simulation with 15 domains. The grid has 168 780 gridpoints. The wavefield is represented with the fourth-order polynomial basis in Fig. 1 (a). There are 12.2 points per P -wavelength and and 7 points per S -wavelength. The time step is t = 8.5 × 10 -4 . Panel (b) is similar to (a) but for a SEM simulation with 10 368 elements and 166 753 gridpoints. The number of points per wavelength and the time step are the same in panels (a) and (b). Panel (c) is similar to panel (a) but the grid is coarser. The grid has 42 630 points. The wavefield is represented using the fourth-order B-spline basis in Fig. 1 (c); There are 6.1 points per P-wavelength and and 3.5 points per S-wavelength. The time step is t = 1.7 × 10 -3 . Panel (d) is similar to panel (c) but for a SEM simulation with 2592 elements and 41 905 gridpoints. The number of points per wavelength and the time step are the same in panels (d) and (c).In panels (a) and (b), there are enough points per w avelength. The P w ave and the S wave are modelled accurately. There is no visible numerical dispersion, and the wavefields obtained using DFDM and SEM are identical. In panels (c) and (d), the number of points per wavelength and the time step are divided by two with respect to panels (a) and (b). In panel (d), we see that the SEM simulation exhibits significant numerical dispersion for the surface wave, as indicated by the red arrows. In panel (c), we observe no numerical dispersion with the DFDM, and the wavefield looks identical to panel (a).

  ), 1 (a), A1 (b), 1 (b), A1 (c) and 1 (c), respecti vel y. Fig.6(c) presents the third simulation which was designed similarly to the second simulation. The only difference is that the model is divided into 625 elements. The grid dimension is 7472. There are 124 conformal interfaces and 1076 non-conformal interfaces. There are 108, 110, 121, 99, 94 and 93 domains using the bases in FigsA1 (a), 1 (a), A1 (b), 1 (b), A1 (c) and 1 (c), respecti vel y.

Figure 8 .

 8 Figure 8. Comparison between the seismograms obtained using DFDM, SEM and the analytical solution when modelling surface wave propagation on a fine grid, as detailed in Figs 7 (a) and (b). Panels (a) and (b) show the horizontal and the vertical displacements, respecti vel y. The recordings correspond to the five recei vers represented b y b lues circles in F ig. 7 . Panels (c) and (d) show the difference (multiplied by 5) between the computed seismograms and the analytical solution.The differences in panels (c) and (d) are cut when the analytical solution predicts no more waves. As expected, the numerical error increases with the distance to the source because of the accumulation of the numerical dispersion. The maximum difference between the analytical solution and the computed seismograms is about 2.5 and 5 per cent of the signals' amplitude for the DFD and SEM seismograms, respectively.

Figure 9 .

 9 Figure 9. Comparison between the seismograms obtained using the DFDM, SEM and the analytical solution when modelling surface wave propagation on a coarse grid, as detailed in Figs 7 (b) and (c). Panels (a) and (b) show the horizontal and the vertical displacements, respecti vel y. The recordings correspond tothe fiv e receiv ers represented by b lues circles in F ig. 7 . The seismograms obtained using the DFDM (dashed blue lines) match perfectly the analytical solution (solid black line). The seismograms obtained using the SEM (dashed red lines) exhibit some spurious ringing in the surface waves because the number of points per wavelength is insufficient. Panels (c) and (d) show the difference (multiplied by 2) between the computed seismograms and the analytical solution. The difference in panels (c) and (d) are cut when the analytical solution predicts no more waves. The maximum difference between the analytical solution and the computed seismograms is about 8 and 25 per cent of the signals' amplitude for the DFD and SEM seismograms, respectively.

Figure 10 .

 10 Figure 10. The snapshots show the magnitude of the displacement vector for four different simulations at time t = 1, 7 s. The magnitude of the displacement is raised to the power 1/3 to enhance waves with a small amplitude. The velocity model has dimensions 6400 m by 3200 m and is divided into two horizontal la yers. The low er la yer is a 1600 m high elastic region with P -wav e v elocity v p = 3400 m s -1 , S -wav e v elocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The upper layer is a 1600 m high acoustic region with P -wave velocity v p = 1500 m s -1 and density ρ = 1020 kg m -3 . The waves are generated by an e xplosiv e source placed at the position marked by a red diamond. The source time function is a Ricker wavelet with peak frequency f 0 = 10 Hz. The receivers are represented using blue circles and the seismograms recorded in the simulations (a, b) and (c, d) are plotted in Figs 11 and 12 , respecti vel y. The principal seismic waves are indicated in panel (b). Panel (a) shows the wavefield obtained for a DFD simulation with 12 domains. The grid has 252444 points. The wavefield is represented using the fourth-order B-spline basis in Fig. 1 (a). There are 12.2 points per P -wavelength and 7 points per S -wavelength in the elastic layer. There are 5.4 points per P -wavelength in the acoustic layer. The time step is t = 8.5 × 10 -4 s. Panel (b) is similar to panel (a) but for a SEM simulation with 15 552 elements and 249 841 gridpoints. The number of points per wavelength and the time step are the same in panels (a) and (b). Panel (c) is similar to panel (a) but is using a coarser grid. The grid size is 64 020 points. the wavefield is represented using the fourth-order B-spline basis in Fig.1 (c). There are 6.1 points per P -wavelength and 3.5 points per S -wavelength in the elastic layer. There are 2.7 points per P -wavelength in the acoustic layer. The time step is t = 1.7 × 10 -3 s. Panel (d) is similar to panel (c) but for a SEM simulation with 3888 elements and 62 713 gridpoints. The number of points per wavelength and the time step are the same in panels (d) and (c). There are enough points per wavelength in panels (a) (b) to model wave propagation accurately. No numerical dispersion is observed, and the wavefields obtained using DFDM and SEM are identical. In panels (c) and (d), the number of points per wavelength and the time step are di vided b y two with respect to panels (a) and (b). In panel (d), we see that the SEM simulation exhibits lots of numerical noise in the acoustic region because there is not enough points per P -wavelength. In panel (c), we observe no numerical dispersion with the DFDM, and the wavefield looks identical to panel (a).

Figure 11 .

 11 Figure 11. Comparison between the seismograms obtained using DFDM, SEM and the analytical solution when modelling wave propagation near a solid-fluid interface on a fine grid, as detailed in Figs 10 (a) and (b). Panels (a) and (b) show the horizontal and the vertical displacements, respecti vel y. The recordings correspond to the ten receivers represented by blues circles in Fig.10. The first fiv e receiv ers are placed in the acoustic region, and the last five receivers are placed in the elastic re gion. P anels (c) and (d) show the difference (multiplied by 10) between the computed seismograms and the analytical solution. The differences in panels (c) and (d) are cut when the analytical solution predicts no more waves. The arrivals observed later in the computed seismograms come from reflections on the free surfaces at the sides of the domain. The maximum difference between the analytical solution and the computed seismograms is about 5 per cent of the signals' amplitude in both simulations.

Figure 12 .

 12 Figure 12. Comparison between the seismograms obtained using DFDM, SEM, and the analytical solution when modelling wave propagation near a solid-fluid interface on a coarse grid, as detailed in Figs 10 (b) and (c). Panels (a) and (b) show the horizontal and the vertical displacements, respecti vel y. The recordings correspond to the 10 receivers represented by blues circles in Fig.10. The first fiv e receiv ers are placed in the acoustic region, and the last fiv e receiv ers are placed in the elastic region. The seismograms obtained using the DFDM (dashed blue lines) match the analytical solution (solid black line). The seismograms obtained using SEM (dashed red lines) exhibit some spurious high-frequency noise generated in the acoustic region where the number of points per wavelength is insufficient. Panels (c) and (d) show the difference (multiplied by 5) between the computed seismograms and the anal ytical solution. The dif ference in panels (c) and (d) are cut when the analytical solution predicts no more waves. The maximum difference between the analytical solution and the computed seismograms is about 20 per cent of the signals' amplitude for the DFD seismograms and 60 per cent for SEM seismograms.

Figure 13 .Figure 14 .

 1314 Figure13. The snapshots show the magnitude of the displacement vector for four different simulations at time t = 1.2 s. The magnitude of the displacement is raised to the power 1/3 to enhance waves having a small amplitude. The velocity model is that of[START_REF] Komatitsch | Wave propagation near a fluidsolid interface: a spectral-element approach[END_REF]. It has dimensions 6400 m by 3200 m and is divided into two layers separated by a sinusoidal interface. The lower layer is an elastic region with P -wave velocity v p = 3400 m s -1 , S -wave velocity v s = 1963 m s -1 and density ρ = 2500 kg m -3 . The upper layer is an acoustic region with P -wave velocity v p = 1500 m s -1 and density ρ = 1020 kg m -3 . The waves are generated by an e xplosiv e source placed at the position marked by a red diamond. The source time function is a Ricker wavelet with peak frequency f 0 = 10 Hz. The receivers are represented using blue circles The seismograms recorded in the simulations are plotted in Fig. 14 . The principal seismic waves observed are indicated in panel (b). Panel (a) shows the wavefield obtained for a DFD simulation with 12 domains. The six domains in the elastic region contain 147 × 83 = 12 201 gridpoints. The six domains in the acoustic region contain 193 × 109 = 21 037 gridpoints. Thus, the grid at the sinusoidal interface is non-conformal. The total grid has 199 428 points. The wavefield is represented using the fourth order polynomial basis in Fig. 1 (a). There is 9.3 points per P -wavelength and 5.4 points per S -wavelength in the elastic region. There are 5.4 points per P -wavelength in the acoustic region. The time step is t = 5 × 10 -4 s. Panel (b) is similar to panel (a) but for a SEM simulation with 15 552 elements and 249 841 gridpoints. There are 12.2 points per P -wavelength and 7.1 points per S -wavelength in the elastic region. There are 5.4 points per P -wavelength in the acoustic region; The grid is conformal at the interface. The time step is the same in panels (a) and (b). The wavefield in panels (a) and (b) look identical.(a) (b)

  , respecti vel y, and, the stresses σ xx , σ yy and σ xy may be equally represented in the bases ˆ

	and ˆ U y or U * coefficient matrices ˆ σ xx , ˆ σ yy and ˆ σ xy or σ *	ˆ B 1 or B *	1 using the expansion coefficient matrices ˆ U x B 2 or B * 2 using the expansion

x and U * y xx , σ * yy and σ * xy , respecti vel y.

  Albin, N. & Klarmann, J., 2016. An algorithmic exploration of the existence of high-order summation by parts operators with diagonal norm, J. Sci. Comput., 69 (2), 633-650. Alford, R., K ell y, K. & Boore, D.M., 1974. Accuracy of finite-difference modeling of the acoustic wave equation, Geophysics, 39 (6), 834-842. Alterman, Z. & Karal, F. Jr, 1968. Propagation of elastic waves in layered media by finite difference methods, Bull. seism. Soc. Am., 58 (1), 367-398.

= M 11 ˆ g( d) 1 of the d th deri v ati ve of f ( x ) and g ( x ) in the dual bases (illustrated in Fig.1)

A P P E N D I X : C O N S T RU C T I O N O F T H E B A S I S F U N C T I O N S

In this section, I detail the construction of the three different kinds of bases used to represent the wavefield and construct the DFD operators. The three bases kinds are piece wise pol ynomials bases where the basis functions have compact support. Using basis function with compact suppor t is cr ucial to building up efficient DFD algorithms because we obtain band diagonal matrices from eq. ( 10). Also, when representing functions with the three bases considered, the end values to the left and the right only depend on the first and the last basis functions, respecti vel y. In other words, the basis functions are zero at the limits of the definition interval except for the first one and the last one. This makes it easier to impose the boundary condition when modelling wave propagation. The first basis kind is a simple polynomial basis; it is based on polynomial interpolation, its basis functions are continuous with discontinuous first derivative ( C 0 continuity). The second basis kind is a cubic Hermite basis; it is constructed using Hermite interpolation, its basis functions, and their first deri v ati ve are continuous ( C 1 continuity). The third basis kind is a B-spline basis; its basis functions and their first few derivatives are continuous ( C p -1 continuity). Both the polynomial and the Hermite bases are interpolating. The expansion coefficients representing a function in those bases are the values of the function at some interpolation points. The B-spline basis is not interpolating. All the bases considered have the partition of unity property.

A1 Polynomial basis

The polynomial basis is constructed from the piecewise interpolating function

with polynomial order p defined on the interval [ a , b ]. The function g ( x ) takes a set of N points

together with their attached data y = [ y 1 , . . . , y N ] and a knot vector

as an input and returns the interpolated value g ( x ) = y ( x ) of the data y at point x . Within each knot interval [ k i , k i + 1 ] the interpolating function takes the value of the interpolating polynomial P i ( x ) passing by the subset of data points with coordinates [( x i , y i ), . . . , ( x i + p , y i + p )]. We have:

The polynomial P i ( x ) may be e v aluated, for example, using Lagrange interpolation. Provided that the interpolating function g ( x ) can be computed, a family of basis functions B i ( i = 1, N ) can be constructed by interpolating the discrete vector basis

In this study, all the examples are computed using e venl y spaced points, that is we have

for i = 1, N . Further, to ensure continuity of the basis functions B i , the knot vector k is taken as a subset of the point vector x . We have:

where the shift parameter l allows to control the degree of symmetry or the skewness of the basis functions B i . The second order bases obtained with parameters p = 2, l = 2, N = 17 and p = 2, l = 3, N = 17 are pictured in Fig. A1 (a) in black and red, respecti vel y, and, the fourth order bases obtained with parameters p = 4, l = 3, N = 17 and p = 4, l = 4, N = 17 are pictured in Fig. 1 (a) in black and red, respecti vel y.

A2 Cubic Hermite basis

The cubic Hermite basis is constructed from the Hermite piecewise interpolating function

where

The three primary kinds of bases used to construct the DFD operators. (a) The black and red lines represent the basis functions of two distinct polynomial bases (1 and 2, respecti vel y) constructed following the approach in Section A1 with polynomial order p = 2 and N = 17; they are continuous with discontinuous first deri v ati ve. The basis functions in black are left-skewed and have been constructed using the shift parameter l = 2. The basis functions in red are right-skewed and have been constructed using the shift parameter l = 3. In this particular example, the basis functions in red overlap the black ones when flipping their left-and right-hand sides. When modelling wave propagation, the displacement is represented using the left-skewed basis while the stress is represented using the right-skewed basis, for example, to prevent spurious modes from appearing. (b) Same as in (a) but using the Hermite bases described in Section A2, the basis functions in black and red have been constructed using the shift parameters l = 0 and l = 2, respecti vel y, they are continuous with continuous first deri v ati ve. (d) Same as in (a) but using the B-spline bases described in Section A3, the basis functions in black and red have been constructed using the knot shift parameters α k = 1/8 and α k = -1/8, respecti vel y, they are continuous as well as their first deri v ati ves.

In this study, the deri v ati ves in vector y are computed from the vectors y and x using

where P i ( x) is the the deri v ati ve of the polynomial P i ( x ) passing by the subset of data points with coordinates [( x i , y i ), . . . , ( x i + p , y i + p )]. As for the polynomial basis, the shift parameter l allows to control the skewness of the basis functions H i . Note that the order p of the polynomials P i ( x ) can be greater than that of the Hermite polynomials h ij ( p = 3). When increasing the order of P i ( x ) the basis functions away from the interval's boundaries will converge toward the ideal interpolant sinc( x ), this allows to reduce the computational error even though the polynomial order of the basis function is constant ( p = 3). As for the polynomial basis, a set of basis functions H i ( i = 1, N ) are constructed by interpolating the discrete vector basis in eq. ( A3 ). We have:

where e i is computed from x i and e i using eq. ( A13).

A3 B-spline basis

Given the knot vector with p + 1 repeated end knots

where a ≤ k i < k i + 1 ≤ b , the B-spline basis functions B

The B-spline basis functions B p i ( x) have C p -1 continuity when the inner knots are all different (i.e. there is no repeated knots). In this study, we construct the B-spline basis using the shifted knot vector k α with components

where the knot shift parameter α k produces basis function with a right shift when α k > 0 and a left shift when α k < 0. The components of the zero shift knot vector k 0 are

The