Yder Masson

Chao Lyu

Peter Moczo

Yann Capdeville

Barbara Romanowicz

Jean Virieux

2D seismic wave propagation using the distributional finite-difference method: further developments and potential for global seismology

Keywords: Computational seismology, Numerical modelling, Wave propagation, Numerical solutions, Seismic anisotropy, Structure of the Earth

We present a time-domain distributional finite-difference scheme based on the Lebedev staggered grid for the numerical simulation of wave propagation in acoustic and elastic media. The central aspect of the proposed method is the representation of the stresses and displacements with different sets of B-splines functions organized according to the stag-2 Masson et al. gered grid. The distributional finite-difference approach allows domain-decomposition, heterogeneity of the medium, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grid, and fluid-solid interfaces. Numerical examples show that the proposed scheme is suitable to model wave propagation through the Earth, where sharp interfaces separate large, relatively homogeneous layers. A few domains or elements are sufficient to represent the Earth's internal structure without relying on advanced meshing techniques.

We compare seismograms obtained with the proposed scheme and the spectral element method, and we show that our approach offers superior accuracy, reduced memory usage, and comparable efficiency.

INTRODUCTION

Over time, seismology has become increasingly dependent on high-performance computing resources to compute the seismic wavefield, and efficient algorithms for solving the wave equation have been developed. Different methods have been adopted within the earthquake seismology community and in exploration geophysics.

Many earthquake seismologists has adopted the Spectral Element Method (SEM), beginning with pioneering modeling of the crust [START_REF] Tape | Adjoint tomography of the southern california crust[END_REF]) and uppermost mantle [START_REF] Fichtner | Full seismic waveform modelling and inversion[END_REF]) at local and regional scales, since then extended to many regions of the world (e.g. [START_REF] Magnoni | Adjoint tomography of the italian lithosphere[END_REF]). Notably, the latest global tomographic models that unveil the structure of our planet's deep interior typically required millions of CPU hours for their construction [START_REF] French | Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography[END_REF][START_REF] Bozdag | Global adjoint tomography: first-generation model[END_REF][START_REF] Lei | Global adjoint tomography-model GLAD-M25[END_REF]. Indeed, until the advent of the SEM, seismic waveform modeling of the Earth's mantle at the global scale relied on the computation of synthetics using normal mode perturbation theory [START_REF] Woodhouse | Mapping the upper mantle: Three-dimensional modeling of earth structure by inversion of seismic waveforms[END_REF][START_REF] Li | Global mantle shear velocity model developed using nonlinear asymptotic coupling theory[END_REF][START_REF] Mégnin | The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms[END_REF], the validity of which required strong assumptions on the strength and smoothness of the 3D heterogeneity.

The introduction of the SEM to global seismology [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]) opened up a new era, providing a comparatively efficient method for the accurate numerical computation of the entire seismic wavefield, in the presence of topography and potentially arbitrary lateral variations of 3D structure.

Given the computational challenges associated with full 3D simulations, where the computational Distributional finite-difference modeling 3 cost increases with the fourth power of frequency, several approaches have been introduced to approximate the wavefield under specific conditions. For example, when assuming an axisymmetric or partly axisymmetric model is a good enough approximation to compute the wavefield generated by localized structures. This led to the development of C-SEM, a "coupled mode-SEM" code (Capdeville et al. 2003a) in which SEM computations in the mantle were coupled through a DtN operator to fast 1D normal mode computations in the core, which was used in the development of the first SEM-based global radially anisotropic shear velocity models of the upper mantle [START_REF] Lekić | Inferring upper-mantle structure by full waveform tomography with the spectral element method[END_REF][START_REF] French | Waveform tomography reveals channeled flow at the base of the oceanic asthenosphere[END_REF]) and of the whole mantle [START_REF] French | Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography[END_REF]. Capdeville et al. (2003b) extended the C-SEM method to an annulus of SEM coupled above and below with 1D modes, enabling the modeling of ultra-low velocity zones at the base of the mantle in 3D, using diffracted S waves at periods down to 8 s [START_REF] Cottaar | An unsually large ULVZ at the base of the mantle near hawaii[END_REF] or even shorter [START_REF] To | 3D effects of sharp boundaries at the borders of the african and pacific superplumes: Observation and modeling[END_REF]. Hypothesizing that it is more important to accurately represent the heterogeneity in the vertical plane containing the source and the receiver, [START_REF] Nissen-Meyer | AxiSEM: broadband 3-D seismic wavefields in axisymmetric media[END_REF] proposed an efficient global 2.5D solver, AxiSEM, in which cylindrical symmetry is assumed around each source. The requirement of cylindrical symmetry in AxiSEM was later released in AxiSEM3D [START_REF] Leng | AxiSEM3D: broad-band seismic wavefields in 3-d global earth models with undulating discontinuities[END_REF], which allows for a smooth variation of the structure in the direction perpendicular to the great circle plane. AxiSEM3D has already received much attention in the global seismological community, in particular for the modeling of complex structures in D" [START_REF] Li | Kilometer-scale structure on the core-mantle boundary near hawaii[END_REF]. Meanwhile, SPECFEM3D GLOBE [START_REF] Tromp | Spectral-element and adjoint methods in seismology[END_REF], which makes no approximations on the structure, is open-source and has been significantly optimized over time, has now become the standard for full seismic wavefield computations at regional and global scales and has been used in particular for the development of global radially anisotropic shear velocity models GLAD-M15 [START_REF] Bozdag | Global adjoint tomography: first-generation model[END_REF]) and GLAD-M25 [START_REF] Lei | Global adjoint tomography-model GLAD-M25[END_REF]. Other groups have developed their own SEM codes, notably in the framework of the CSEM (Collaborative Seismic Earth Model, e.g. [START_REF] Fichtner | The collaborative seismic earth model: generation 1[END_REF] or in the framework of regional and crustal seismic imaging [START_REF] Trinh | Efficient 3D time-domain elastic and viscoelastic Full Waveform Inversion using a spectral-element method on flexible Cartesian-based mesh[END_REF][START_REF] Lu | High-resolution surface wave tomography of the european crust and uppermost mantle from ambient seismic noise[END_REF]). These codes are however not open source. Although very higher-order spectral element methods, such as the 12th-order method, are more computationally efficient than the classical 4th-order SEM [START_REF] Lyu | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation[END_REF], the overall computational efficiency of SEM decreases significantly when its order exceeds 20 due to the extremely small available time step, which limits the usefulness of super-high-order (≥ 20th-order) SEM for global-scale simulations.

Still, there is a need for more efficient methods to reach the higher frequencies required for the study of complex structures in the deep Earth, as the computational burden of the SEM is prohibitive for global waveform modeling at periods shorter than 10 s. In this study, we investigate an alternative to SEM for modeling seismic wave propagation at the global scale called the distributional finitedifference method (DFDM) introduced in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF]. The proposed algorithm shall allow for do-main decomposition, medium heterogeneity, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grids, free surfaces, and fluid-solid interfaces.

SEM is the favored approach to model wave propagation at the planet's scale; on the other hand, the finite-difference method (FDM) is predominant in exploration geophysics and earthquake ground motion simulations in local surface sedimentary structures (e.g. [START_REF] Virieux | A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging[END_REF][START_REF] Moczo | The finite-difference modelling of earthquake motions: Waves and ruptures[END_REF][START_REF] Moczo | Numerical wave propagation simulation[END_REF]. This apparent discrepancy may be due to communities of practice, software availability, and, most likely, the nature of the seismic data and the model geometry. Active seismic imaging relies primarily on body wave analysis with a focus on reflections, while surface waves represent an essential part of the data processed in seismology and global tomography. Thus, SEM, especially suited for surface wave modeling, is preferred in seismology. FDM needs to be adapted to model surface waves accurately. Otherwise, it is remarkably simple and efficient: an attractive feature in exploration geophysics. Fluid-solid interfaces also need to be modeled accurately and influence the method to be chosen [START_REF] De Basabe | A comparison of finite-difference and spectral-element methods for elastic wave propagation in media with a fluid-solid interface[END_REF]. There are, of course, numerous exceptions to the general statements above. FDM has been successfully used to model wave propagation at the global and regional scale [START_REF] Igel | SH-wave propagation in the whole mantle using high-order finite differences[END_REF]Igel et al. 1995;[START_REF] Chaljub | Sensitivity of ss precursors to topography on the upper-mantle 660-km discontinuity[END_REF][START_REF] Jahnke | Global SH-wave propagation using a parallel axi-symmetric finite-difference scheme[END_REF][START_REF] Kawai | Complete synthetic seismograms up to 2 Hz for transversely isotropic spherically symmetric media[END_REF][START_REF] Li | Global synthetic seismograms using a 2-d finitedifference method[END_REF]. Summation-by-part finite-difference schemes, named SBP-SAT approaches, which naturally account for the free surface, have reached recently an improved accuracy [START_REF] Sjögreen | A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation[END_REF][START_REF] Petersson | Wave propagation in anisotropic elastic materials and curvilinear coordinates using a summation-by-parts finite difference method[END_REF]. Such an advanced SBP-SAT software is openly available to the seismology community (SW4, Petersson & Sjogreen 2017). SEM has also been employed for exploration and geotechnical applications [START_REF] Zheng | Spectral-element simulations of elastic wave propagation in exploration and geotechnical applications[END_REF]. The literature in FDM is abundant, and we refer the reader to review papers [START_REF] Virieux | A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging[END_REF][START_REF] Moczo | The finite-difference modelling of earthquake motions: Waves and ruptures[END_REF][START_REF] Moczo | Numerical wave propagation simulation[END_REF], and the references therein for a comprehensive review.

We focus on the staggered-grid finite-difference approach, which is relevant in this study. We highlight the pioneering work of [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF] in electromagnetism, [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF] in seismic source modeling, and [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] in seismic wave propagation. In staggered-grid finite-difference approaches, the field variables are discretized on different grids, leading to centered finite-difference operators with improved accuracy for spatial derivatives of fields. Further, the staggered-grid approach effectively annihilates the spurious (non-physical) parasitic modes that are problematic in FDM (see, e.g., [START_REF] New | Stable 3-level leapfrog integration in numerical relativity[END_REF][START_REF] Dovgilovich | High-accuracy finite-difference schemes for solving elastodynamic problems in curvilinear coordinates within multiblock approach[END_REF]. Thus, this staggered-grid FDM is praised for its simplicity and computational efficiency and is widely used for the 3D modeling of seismic wave propagation [START_REF] Graves | Modeling three-dimensional site response effects in the marina district basin, san francisco, california[END_REF][START_REF] Olsen | Magnitude 7.75 earthquake on the san andreas fault: three-dimensional ground motion in los angeles[END_REF][START_REF] Pitarka | Three-dimensional simulation of the near-fault ground motion for the 1995 hyogo-ken nanbu (kobe), japan, earthquake[END_REF][START_REF] Cotton | Three-dimensional simulation of earthquakes in the grenoble's basin[END_REF][START_REF] Matsushima | 3D simulation of aftershocks of the Hyogoken nanbu earthquake of 1995[END_REF][START_REF] Wald | The seismic response of the los angeles basin, california[END_REF][START_REF] Kristek | The 1995 Kobe mainshock simulated by the 3D finite differences[END_REF]. The former Virieux algorithm [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF]) permits limited anisotropy (i.e. orthotropic) and can accommodate orthogonal curvilinear coordinates. Several staggered-grid FDM schemes have been introduced to account for general anisotropy and non-Distributional finite-difference modeling 5 orthogonal curvilinear grids, notably, partly-staggered grid scheme (e.g., [START_REF] Saenger | Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[END_REF]) and Lebedev scheme (e.g., [START_REF] Lisitsa | Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity[END_REF][START_REF] De La Puente | Mimetic seismic wave modeling including topography on deformed staggered grids[END_REF]. Even though they do not strictly fall in the staggered-grid FDM category, one should also mention the collocated approaches relying on the MacCormack scheme where left and right FD operators are employed to prevent the parasitic modes from appearing [START_REF] Zhang | Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[END_REF][START_REF] Zhu | Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method[END_REF][START_REF] Lan | Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface[END_REF][START_REF] Sun | Seismic-wave modeling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method[END_REF]. [START_REF] Bernth | A comparison of the dispersion relations for anisotropic elastodynamic finite-difference grids[END_REF] emphasizes that the Lebedev scheme is more efficient and accurate than the partly staggered grid scheme. The stability and grid dispersion of the staggered-grid FDM is analyzed carefully in [START_REF] Moczo | 3D fourth-order staggered-grid finite-difference schemes: Stability and grid dispersion[END_REF]. [START_REF] Kristek | Efficient methods to simulate planar free surface in the 3d 4th-order staggered-grid finite-difference schemes[END_REF] and [START_REF] Bohlen | Accuracy of heterogeneous staggered-grid finite-difference modeling of rayleigh waves[END_REF] addressed the difficulty of accurately modeling the free surface using staggered-grid FDM. Nonetheless, several approaches have been introduced to address this issue [START_REF] Lombard | Free and smooth boundaries in 2-d finite-difference schemes for transient elastic waves[END_REF]Zhang et al. 2012a,b;[START_REF] De La Puente | Mimetic seismic wave modeling including topography on deformed staggered grids[END_REF][START_REF] Gao | An immersed free-surface boundary treatment for seismic wave simulation[END_REF][START_REF] Sun | Seismic-wave modeling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method[END_REF][START_REF] Shragge | Tensorial elastodynamics for isotropic media[END_REF][START_REF] Zang | An overset-grid finite-difference algorithm for simulating elastic wave propagation in media with complex free-surface topography[END_REF]). An aspect of special importance in the finite-difference modeling in both seismic exploration and ground motion simulations is the implementation of material interfaces. It is not trivial in the case of a general shape and position in the uniform grid. The problem has been recently addressed by [START_REF] Mittet | On the internal interfaces in finite-difference schemes[END_REF], [START_REF] Jiang | TTI equivalent medium parametrization method for the seismic waveform modelling of heterogeneous media with coarse grids[END_REF], [START_REF] Koene | Finite-difference modelling of 2-D wave propagation in the vicinity of dipping interfaces: a comparison of anti-aliasing and equivalent medium approaches[END_REF] and [START_REF] Moczo | Material interface in the finite-difference modeling: a fundamental view[END_REF]. In this work, we build up the Lebedev approach to construct our DFDM algorithm.

Besides SEM, numerous methods that decompose the computational domain in multiple elements exist. In particular, the classical finite-element method (FEM), [START_REF] Lysmer | A finite element method for seismology[END_REF][START_REF] Smith | The application of finite element analysis to body wave propagation problems[END_REF][START_REF] Marfurt | Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[END_REF]. Boundary-element methods (BEM) [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][START_REF] Bouchon | Calculation of synthetic seismograms in a laterally varying medium by the boundary element-discrete wavenumber method[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) [START_REF] Käser | Quantitative accuracy analysis of the discontinuous galerkin method for seismic wave propagation[END_REF][START_REF] De Basabe | The interior penalty discontinuous galerkin method for elastic wave propagation: grid dispersion[END_REF][START_REF] De La Puente | Discontinuous galerkin methods for wave propagation in poroelastic media[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]. Mixed FEM approaches following partially the staggered-grid strategy have been developed [START_REF] Nédélec | Mixed finite elements in R 3[END_REF][START_REF] Bécache | A new family of mixed finite elements for the linear elastodynamic problem[END_REF]. Over the last decade, Isogeometric Analysis (IGA) [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Auricchio | Isogeometric collocation methods[END_REF]) has garnered significant attention within the field of structural engineering. This approach discretize Partial Differential Equations (PDEs) using the same type of functions employed in Computer Aided Design (CAD): Non-Uniform Rational B-Splines (NURBS). These NURBS not only provide a precise representation of the geometric aspects of the computational domain but also serve as the basis for defining the solution space for PDEs. Element-based approaches allow to mesh complex geological features and are well suited to model free surfaces accurately and discontinuities thanks to their weak formalism. However, they lead to several numerical challenges. They are more difficult to implement than FDM, often more expensive in computational time, and their accuracy depends on the quality of the meshing [START_REF] Virieux | A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging[END_REF]. One should also mention pseudo-spectral methods (PSM) that are highly accurate but less flexible for representing media with complex geometries [START_REF] Fornberg | The pseudospectral method: Comparisons with finite differences for the elastic wave equation[END_REF][START_REF] Carcione | A 2d chebyshev differential operator for the elastic wave equation[END_REF]. [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] introduced a distributional finite-difference method (DFDM) that combines some essential features of FDM, FEM, SEM, IGA, and DGM. DFDM decomposes the space into multiple elements with arbitrary sizes. Thus, a volumetric mesh may be employed to describe the model better whenever needed. Within the elements, DFDM has an algorithmic structure similar to FDM and relies on compactly supported operators, which ensures efficiency. Our algorithm relies on B-spline bases similarly to IGA with two important differences. First, DFDM does not require designing optimal quadrature rules (e.g. [START_REF] Zou | Efficient and robust quadratures for isogeometric analysis: Reduced gauss and gauss-greville rules[END_REF]) for achieving efficiency. The mass and stiffness matrices are assembled independently prior to the computations, while the element shapes and heterogeneity are accounted for approximately using weights associated with virtual grid points. This approach is wellsuited for geophysical applications where the exact structure of the medium in which waves propagate is rarely known precisely. Second, in DFDM, field variables are represented using different basis functions (similarly to staggered grid FDM). This has the dual benefit of reducing the number of degrees of freedom in the problem and eliminating parasitic modes that are notoriously problematic in wave propagation modeling [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF][START_REF] Lisitsa | Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity[END_REF]. Finally, in DFDM, the numerical wavefield is discontinuous between neighboring elements, both the displacement and the stress can be defined explicitly at the boundaries, similarly to DGM.

The algorithm proposed by [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] relies on a MacCormack type of strategy to remove the parasitic modes; all the displacements are represented in a first basis (with the basis function skewed to the left), and all the stresses are represented in a second basis (with basis functions skewed to the right). This algorithm also shares the structure of the partly-staggered-grid schemes with all the displacement components located at one grid position and all the stress components located at another grid position (e.g., [START_REF] Saenger | Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[END_REF]. [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] proposed a DFDM scheme analogous to the classic algorithm from [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF] where field variables are represented in different bases organized according to the staggered grid: specific differential operators switching consistently from one basis to the other one when estimating a field and its spatial derivative are designed as we shall also see in this work (this operation is not performed in mixed FEM approaches). They showed that the DFDM algorithm is more accurate and faster than its FDM counterpart. Interestingly, the staggered approach permits the computation of the partial derivatives by applying 1D operators in a single spatial direction.

In contrast, the method of [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] requires using the 1D operators in the two spatial directions.

The switch between bases is still performed in both strategies. The two latter approaches also differ in the time scheme. [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] employs a centered second-order time integration scheme with a displacement formulation and [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] employ a first-order leap-frog strategy. Last, the Distributional finite-difference modeling 7 two studies account for material heterogeneity in different manners. [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF], includes material heterogeneity using a mass lumping technique and [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] discretizes the material properties on a virtual grid associated with the basis functions. For clarity and to outline the similarities with FDM, [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] only considered a single isotropic element. Nonetheless, their algorithm applies to limited anisotropy (i.e., orthotropic material) and adapts to orthogonal curvilinear grid (e.g., [START_REF] Xie | An explicit fourth-order orthogonal curvilinear staggered-grid fdtd Distributional finite-difference modeling 37 method for maxwell's equations[END_REF].

In this study, we expand the [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] algorithm to the Lebedev grid, which consists of two Virieux grids. It is needed in seismology to model wave propagation in arbitrary anisotropic media and to work with non-orthogonal curvilinear coordinates. Further, it allows the meshing of complex geological features and accurately model interfaces, particularly free surfaces and fluid-solid interfaces with topography. A key aspect of this paper is the evaluation of the proposed algorithm for modeling seismic wave propagation through the Earth at the global scale. This manuscript is organized as follows. The first section starts with the second-order elastodynamic system of equations. We apply the non-orthogonal curvilinear coordinate transformation. Then, we discretize the field variables using the basis functions organized according to the Lebedev grid.

Thereafter, we briefly recall how a dual basis can be employed to approximate functions and their derivatives. We proceed with the derivation of the one-dimensional DFDM operators employed for approximating the partial derivatives. The 1D operators are then extended to 2D using a tensor product. The treatment of boundary conditions and interfaces is discussed together with the introduction of the 2D operators. Finally, we detail our 2D time-domain algorithm that models wave propagation.

In the second section, we present numerical examples demonstrating that the proposed numerical algorithm can accurately and efficiently model wave propagation through the Earth at the global scale.

We compare the results obtained with DFDM to those obtained using SEM. DFDM is significantly more accurate than SEM when using a similar setup. We finish with concluding remarks and future perspectives.

METHODOLOGY

Governing equations

We consider wave propagation in a two-dimensional linear anisotropic elastic medium. The proposed algorithm relies on the second-order elastodynamic system of equations.

ρü x = ∂ x σ xx + ∂ y σ xy + f x (1a) ρü y = ∂ x σ xy + ∂ y σ yy + f y (1b) σ xx = C 11 ∂ x u x + C 12 ∂ y u y + C 13 (∂ x u y + ∂ y u x ) (1c)
σ yy = C 12 ∂ x u x + C 22 ∂ y u y + C 23 (∂ x u y + ∂ y u x ) (1d)
σ xy = C 13 ∂ x u x + C 23 ∂ y u y + C 33 (∂ x u y + ∂ y u x ) (1e)
where

u i is the displacement component in direction i, f i is the component of the source in direction i,
σ ij are the component of the stress tensor, ρ is the density, C ij can be inferred from elastic parameters of the linear Hooke law, ∂ t is the partial derivative with respect to time, ∂ i is the partial derivative with respect to spatial direction i.

Curvilinear transformation

To allow for curvilinear meshes that follow the topography of the interfaces, we apply a smooth, locally invertible coordinate transformation

x = x(x ′ , y ′ ) y = y(x ′ , y ′ ) (2)
where x and y are the cartesian coordinates in the physical space and x ′ and, y ′ are the coordinates in the computational space. We illustrate such a transformation in Fig. 1. Using the transformation in Equation 2, we can write Equation 1 in the conservative form [START_REF] Thompson | Numerical grid generation[END_REF][START_REF] Appelö | A stable finite difference method for the elastic wave equation on complex geometries with free surfaces[END_REF][START_REF] Dovgilovich | High-accuracy finite-difference schemes for solving elastodynamic problems in curvilinear coordinates within multiblock approach[END_REF] Distributional finite-difference modeling 9

x ′ y ′ x ′ - x ′ + y ′ - y ′ + Ω Ω x- Ω x+ Ω y- Ω y+ x y x = x( ( (x ′ , y ′ ) ) ) y = y( ( (x ′ , y ′
) ) )

Computational space

Physical space

Figure 1. Illustration of the geometrical mapping between the physical space coordinates (x, y) and the com-

putational space coordinates (x ′ , y ′ ). |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 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (3a) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ xx σ yy σ xy ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C 11 C 12 C 13 C 12 C 22 C 23 C 13 C 23 C 33 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂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 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (3b)
where

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 ′ (3c)
and |J| is the determinant of the Jacobian matrix J. We need to transform the linear system above into a discrete linear system to obtain a numerical solution.

Discrete representation of the field variables

We discretize the field variables in Equation 3 following the approach introduced in Masson ( 2022)

and [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF]. We represent the wavefield using orthogonal basis functions that are staggered according to the so-called Lebedev grid for the following reasons:

• Orthonormal bases maintain the self-adjoint antisymmetric nature of the wave equation after discretization [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF]. This is essential to achieve numerical reciprocity and ensure the stability of the numerical scheme (subject to a CFL condition).

• Staggered basis functions ensure that the spurious(non-physical) parasitic modes associated with the centered finite-difference operators are effectively removed (see, e.g., [START_REF] New | Stable 3-level leapfrog integration in numerical relativity[END_REF].

In two dimensions, the Lebedev grid consists of two Here, the upper indices indicate the basis in which the field variable are represented in the two spatial directions. For example, u 12 is represented using the basis B 1 in direction x and the basis B 2 in direction y. Each set of variables is organized following the staggered grid [START_REF] Virieux | P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[END_REF]. For the first Virieux grid, we have and for the second Virieux grid we have

u 21 x (x ′ , y ′ , t) = N 2 x ∑ i=1 N 1 y ∑ j=1 U 21 xij (t) B21 ij (x ′ , y ′ ) = N 2 x ∑ i=1 N 1 y ∑ j=1 U 21 xij (t) B2 xi (x ′ ) B1 yj (y ′ ) (4a) u 12 y (x ′ , y ′ , t) = N 1 x ∑ i=1 N 2 y ∑ j=1 U 12 yij (t) B12 ij (x ′ , y ′ ) = N 1 x ∑ i=1 N 2 y ∑ j=1 U 12 yij (t) B1 xi (x ′ ) B2 yj (y ′ ) (4b) σ 11 xx (x ′ , y ′ , t) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 xxij (t) B11 ij (x ′ , y ′ ) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 xxij (t) B1 xi (x ′ ) B1 yj (y ′ ) (4c) σ 11 yy (x ′ , y ′ , t) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 yyij (t) B11 ij (x ′ , y ′ ) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 yyij (t) B1 xi (x ′ ) B1 yj (y ′ ) (4d) σ 22 xy (x ′ , y ′ , t) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 xyij (t) B22 ij (x ′ , y ′ ) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 xyij (t) B2 xi (x ′ ) B2 yj (y ′ ) (4e)
u 12 x (x ′ , y ′ , t) = N 1 x ∑ i=1 N 2 y ∑ j=1 U 12 xij (t) B12 ij (x ′ , y ′ ) = N 1 x ∑ i=1 N 2 y ∑ j=1 U 12 xij (t) B1 xi (x ′ ) B2 yj (y ′ ) (5a) u 21 y (x ′ , y ′ , t) = N 2 x ∑ i=1 N 1 y ∑ j=1 U 21 yij (t) B21 ij (x ′ , y ′ ) = N 2 x ∑ i=1 N 1 y ∑ j=1 U 21 yij (t) B2 xi (x ′ ) B1 yj (y ′ ) (5b) σ 22 xx (x ′ , y ′ , t) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 xxij (t) B22 ij (x ′ , y ′ ) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 xxij (t) B2 xi (x ′ ) B2 yj (y ′ ) (5c) σ 22 yy (x ′ , y ′ , t) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 yyij (t) B22 ij (x ′ , y ′ ) = N 2 x ∑ i=1 N 2 y ∑ j=1 S 22 yyij (t) B2 xi (x ′ ) B2 yj (y ′ ) (5d) σ 11 xy (x ′ , y ′ , t) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 xyij (t) B11 ij (x ′ , y ′ ) = N 1 x ∑ i=1 N 1 y ∑ j=1 S 11 xyij (t) B1 xi (x ′ ) B1 yj (y ′ ) (5e)
where the Bkl ij (k = 1, 2) are the two-dimensional orthonormal B-spline basis functions. They are constructed from the one-dimensional basis functions Bk

x ij (k = 1, 2) and Bl y ij (l = 1, 2) using a tensorial product. The orthonormal bases Bkl ij (k = 1, 2) may be constructed from arbitrary basis functions B kl ij (k = 1, 2) that we will refer to as the canonical bases. In this study, the canonical bases are taken as the two-dimensional B-spline bases detailed in Appendix A and pictured in Figure 2.

Approximation using dual bases

Before constructing the DFD operators, we show that the orthonormal bases employed for representing the field variables can be used to approximate functions and their derivatives. The orthonormal bases Bk (k=1,2) defined in Appendix C are self-dual. Thus, they provide us with a direct way to find the approximation f k (θ) of an arbitrary function f (θ), where θ is a spatial coordinate (see, e.g., Woźny 2013). We have

f k (θ) = N k ∑ i=1 f k i Bk i (θ) ≈ f (θ) (6)
where the approximation f k (θ) is represented in the basis Bk and the expansion coefficients f k i are obtained directly using the integrals

f k i = ∫ θ+ θ- Bk i (θ)f (θ)dθ. (7)
In Equation 6, the function f k (θ) is an optimal approximation of f (θ) in the least square sense; it

minimizes the norm ||f -f k || in the interval [θ -, θ + ].
Notably, the approximation error vanishes and we obtain an exact result (i.e., ||f -f k || = 0) when the function f (θ)=f k (θ) can be represented in the basis Bk . We can approximate the derivative f ′ (θ) of the function f (θ) with a similar approach. We Distributional finite-difference modeling 13 have

f ′ k (θ) = N k ∑ i=1 f ′ k i Bk i (θ) ≈ f ′ (θ), (8)
where f ′ k i are the expansion coefficients representing the approximation f ′ k (θ) of f ′ (θ) in the basis Bk . Let us notice that this basis for the derivative might not be the one used for approximating the function itself. By substituting f (θ) with f ′ (θ) in Equation 7and using integration by parts, we obtain the expansion coefficients

f ′ k i = ∫ θ+ θ- Bk θi (θ)f ′ (θ)dθ (9a) = -∫ θ+ θ- ( Bk θi ) ′ (θ)f (θ)dθ -Bk θi (θ -)f (θ -) + Bk θi (θ + )f (θ + ). ( 9b
)
where ( B) ′ are the derivatives of B. Equations 7 and 9 are employed to approximate the derivatives of the field variables and thus construct the DFD operators. Equations 7 and 9 are central in DFDM;

they allow one to compute the value/derivative of a function in a basis different from the one in which they are represented. It is similar to what is done in staggered-grid FDM, where we take wavefield values on one grid and compute their derivative on the corresponding shifted grid. This is the key point regarding the staggered-grid strategy we promote.

Distributional finite-difference operators in 1D

We now construct the one-dimensional DFD operators that are the basic bricks of our algorithm. Following [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF], they should act on the field variables expressed in the first basis and return their derivatives' representations in the second basis. We assimilate the function f (θ) introduced in the previous section to a discrete field variable f l (θ) represented in the orthogonal basis Bl . We assume that the discrete representation f l of f l (θ) is known. We can express f l (θ) and its

derivative f ′ l (θ) as f l (θ) = N l ∑ j=1 f l j Bl j (θ) (10a) f ′ l (θ) = N l ∑ j=1 f l j ( Bl j ) ′ (θ) (10b)
where f l j are the expansion coefficients stored in the vector f l , Bl j are the basis functions, and ( Bl j ) ′ are their derivatives. Keeping in mind that the upper indices k and l indicate the basis in which the variables are represented; our objective is to obtain the approximations

f k (θ) ≈ f l (θ) and f ′ k (θ) ≈ f ′ l (θ)
expressed in Bk of the function and its derivative in Equation 10. By substituting Bk i (θ) and f l (θ) in Equation ( 7) with their expressions in Equations (C.2) and (10a) we obtain the linear system

f k = T kl ⋅ f l . ( 11
)
where

T kl = (L k ) -1 ⋅ M kl ⋅ (L l T ) -1 . ( 12
)
The mass matrix M kl and its Cholesky factorization L n are defined in Appendix B and Appendix C.

The operator T kl can be thought of as a basis transform; it acts on the vector f l representing the function f l (θ) in the basis Bl and returns the vector

f k representing the approximation f k (θ) ≈ f l (θ)
in the basis Bk .

By substituting Bk i (θ) and f ′ l (θ) in eq. ( 9) with their expressions in eqs. (C.2) and ( 10b) we obtain the linear system associated with the derivative

f ′ k = D kl ⋅ f l -b k f l (θ -) + d k f l (θ + ), (13)
where

D kl = -(L k ) -1 ⋅ K lk T ⋅ (L l T ) -1 . ( 14
)
and

b k = (L k ) -1 p k (15a) d k = (L k ) -1 q k . ( 15b
)
The stiffness matrix K lk and the vectors p k and q k are defined in Appendix B. When the boundary values f l (θ -) and f l (θ + ) are equal to the first and last elements of the vector f l , the linear system in Equation 13becomes

f ′ k = D kl ⋅ f l , (16)
and by definition, the first-order distributional finite-difference operator is

D kl = (L k ) -1 ⋅ [-K lk T -p k ⋅ p l T + q k ⋅ q l T ] ⋅ (L l T ) -1 . ( 17
)
The operator D kl acts on the vector f l representing the function f l (θ) in the basis Bl . It returns the

vector f ′ k representing the approximation f ′ k (θ) ≈ f ′ l (θ) in the basis Bk .
Notice that the factorization employed here is different (but equivalent) to that employed in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] where the DFD operators have an adjoint form and partly encompass the boundary values.

In this paper, the self-adjoint form of the global or assembled linear system is recovered by using an appropriate averaging of the boundary values between neighboring elements, as detailed in the next section. Notice that all the matrices involved in calculating the DFD operators in Equations 12 and 14

Distributional finite-difference modeling 15 are band-diagonal matrices because the B-spline canonical basis functions have compact support. It ensures the efficiency of the proposed algorithm.

Distributional finite-difference operators in 2D

We construct the 2D-DFD operators needed to compute the partial derivative of a function using a tensor product, i.e., by applying the 1D-DFD operators in the spatial dimensions x ′ and y ′ as detailed in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF].

Let f kl (x ′ , y ′ ) be one of the field variables in Equations 4 and 5. It is represented in the basis Bkl through its expansion coefficients F kl mn stored in the matrix F kl . Throughout this paper, we assume that the indices (i,j) of the arrays correspond to (row, column) and are associated with directions (x,y).

We have

f kl (x ′ , y ′ ) = N k x ∑ m=1 N l y ∑ n=1 F kl mn Bkl mn (x ′ , y ′ ) = N k x ∑ m=1 N l y ∑ n=1 F kl mn Bk xm (x ′ ) Bl yn (y ′ ) ( 18
)
where the 2D basis functions Bkl mn (x ′ , y ′ ) = Bk xm (x ′ ) Bl yn (y ′ ) are constructed from the 1D basis functions Bk xm (x ′ ) and Bl yn (y ′ ). According to the staggered grid, one should represent the partial derivative ∂ ∂x f kl (x ′ , y ′ ) in the basis Bil (with l ∈ [1, 2] and i = 3 -l) and the partial derivative ∂ ∂y f kl (x ′ , y ′ ) in the basis Bkj (with k ∈ [1, 2] and j = 3 -k). We have

∂ ∂x f kl (x ′ , y ′ ) ≈ N i x ∑ m=1 N l y ∑ n=1 D il Xmn Bil mn (x ′ , y ′ ) = N i x ∑ m=1 N l y ∑ n=1 D il Xmn Bi xm (x ′ ) Bl yn (y ′ ) (19) ∂ ∂y f kl (x ′ , y ′ ) ≈ N k x ∑ m=1 N j y ∑ n=1 D kj Ymn Bkj mn (x ′ , y ′ ) = N k x ∑ m=1 N j y ∑ n=1 D kj Ymn Bk xm (x ′ ) Bj yn (y ′ ) ( 20
)
where the expansion coefficients D il Xmn and D kj Ymn stored in the matrices D il X and D kj Y represent the approximations ∂ ∂x f kl (x ′ , y ′ ) and ∂ ∂y f kl (x ′ , y ′ ) in the bases Bil and Bkj . Notice that the function f kl (x ′ , y ′ ) and its partial derivatives are represented in the same 1D basis perpendicular to the partial derivative's direction. Following the approach in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF], the expansion coefficients representing the partial derivatives can be computed using

D ij X (F kj ) = D ik x ⋅ F kj -b i x ⋅ f j T (x ′ -) + d i x ⋅ f j T (x ′ + ) (21a) D ij Y (F il ) = F il ⋅ D jl T y -f i (y ′ -) ⋅ b j T y + f i (y ′ + ) ⋅ d j T y , (21b)
where In the domain Ω, the internal boundary values in Equation 21shall be set to impose the boundary conditions or computed from the internal expansion coefficients F ij using

F
f j (x ′ -) = [f j (x ′ -)] Ω = [F kj T ⋅ b k x ] Ω (22a) f j (x ′ + ) = [f j (x ′ + )] Ω = [F kj T ⋅ d k x ] Ω (22b) f i (y ′ -) = [f i (y ′ -)] Ω = [b i T x ⋅ F il ] Ω (22c) f i (y ′ + ) = [f i (y ′ + )] Ω = [d i T x ⋅ F il ] Ω . (22d)
Here, the brackets [] Ω indicate that the expansion coefficients are taken inside the current domain Ω.

To model a free surface, for example, one can take boundary values of the traction equal to zero when evaluating the partial derivatives of the stresses and use the internal values when evaluating the partial derivatives of the displacements. Similarly, for a domain with fixed boundaries, one can take boundary values of displacement equal to zero when evaluating the partial derivatives of displacements and use the internal values in Equation 22 when evaluating the partial derivatives of the stresses.

We average the boundary values between neighboring domains to connect multiple domains or elements. This continuity condition is needed to ensure numerical reciprocity and the stability of the numerical scheme [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF]. Notice however, that the continuity of the wavefield is not satisfied numerically and the values to the left and right of an interface may be different.The partial derivatives in Equation 21 are evaluated in the domain Ω using of the current domain Ω. This operation can be thought of as an approximate basis transform. We have

f j (x ′ -) = (α x ′ -) [f j (x ′ -)] Ω + (1 -α x ′ -) T jj y Ω←Ω x ′ - ⋅ [f j (x ′ + )] Ω x ′ - (23a) f j (x ′ + ) = (α x ′ + ) [f j (x ′ + )] Ω + (1 -α x ′ + ) T jj y Ω←Ω x ′ + ⋅ [f j (x ′ -)] Ω x ′ + (23b) f i (y ′ -) = (α y ′ -) [f i (y ′ -)] Ω + (1 -α y ′ -) T ii x Ω←Ω y ′ - ⋅ [f i (y ′ + )] Ω y ′ - (23c) f i (y ′ + ) = (α y ′ + ) [f i (y ′ + )] Ω + (1 -α y ′ + ) T ii x Ω←Ω y ′ + ⋅ [f i (y ′ -)] Ω y ′ + ( 23d
T nn θ Ω←Ω n = [L n θ ] -1 Ω ⋅ M nn θ Ω←Ω n ⋅ [L n θ T ] -1 Ω n (θ = x, y; n = i, j) (24a) M nn θ Ω←Ω n def = (M nn θ Ω←Ω n ) ij = ∫ θ+ θ- [B n θi (θ)] Ω [B n θj (θ)] Ω n dθ (θ = x, y; n = i, j) (24b)
where the quantities within the square brackets [] Ω and [] Ω n are those associated with the current Ω and the neighboring Ω n domains. The matrices L n θ are defined in Equation C.1. The matrix T nn θ in Equation 24a is similar to the zeroth-order DFD operator in Equation 12 but acts on the boundary values taken from the neighboring domains [] Ω n . Accordingly, the matrix M nn θ in Equation 24b is similar to the mass matrix in Equation B.1 but the second basis functions B n θj (θ) are taken from the neighboring domain Ω n . When the wavefield is represented in the same basis along the interface between two elements (i.e., for conformal interfaces), the matrices T nn θ in Equation 23 may be omitted because T nn θ = I. To obtain a stable algorithm that satisfies numerical reciprocity, the boundary values must be uniquely defined between neighboring domains [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF]. We must apply the same weighting scheme in neighboring elements. For example, one enforces the following relations between the five domains pictured in Figure 1:

[α x ′ -] Ω = 1 -[α x ′ + ] Ω x ′ - ( 25
) [α x ′ + ] Ω = 1 -[α x ′ -] Ω x ′ + (26) [α y ′ -] Ω = 1 -[α y ′ + ] Ω y ′ - ( 27
) [α y ′ + ] Ω = 1 -[α y ′ -] Ω y ′ + (28)
Here, [α a ] Ω denotes the weight associated with the face a of the domain Ω a . The continuity conditions in Equation 25 ensure that the global differential operators (that approximate the partial derivatives in all domains taken together) will be adjoint after applying the boundary conditions [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF].

Thus far, we derived the DFD operators approximating the partial derivatives in Equation 21 independently of the wave equation. We can use these operators as such to solve the elastodynamic system in Equations 3.

2D algorithm

Our algorithm is obtained by substituting the partial derivatives in Equations 3 with the DFD operators in Equation 21and by substituting the physical quantities by their discrete counterparts. In this work, we adopt the virtual point approach introduced in [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF]. The material properties and quantities associated with the curvilinear mapping are simply evaluated at the virtual grid points, as detailed in Appendix D. We approximate the second-time derivative using a centered second-order finite-difference scheme. Thus, Equations 3 are solved using a recursive procedure. For clarity, we

give the update expressions for the variables associated with the first Virieux grid in Equations 4.

The update expressions for the second Virieux grid in Equation 5 are similar and are obtained by performing the index substitutions (1 → 2; 2 → 1) in Equations (29-31) below. At each step, one first evaluates the stresses in each domain at time t using:

Σ 11 xx = [(C 11 11 ○ J -1 xx 11 + C 11 13 ○ J -1 xy 11 ) ○ D 11 X (U 21 x )](t) (29a)
+ M 11 xx (t) (29e)
Σ 11 yy = [(C 11 12 ○ J -1 xx 11 + C 11 23 ○ J -1 xy 11 ) ○ D 11 X (U 21 x )](t) (29f) + [(C 11 12 ○ J -1 yx 11 + C 11 23 ○ J -1 yy 11 ) ○ D 11 Y (U 12 x )](t) (29g) + [(C 11 22 ○ J -1 xy 11 + C 11 23 ○ J -1 xx 11 ) ○ D 11 X (U 21 y )](t) (29h) + [(C 11 22 ○ J -1 yy 11 + C 11 23 ○ J -1 yx 11 ) ○ D 11 Y (U 12 y )](t) (29i) + M 11 yy (t) (29j)
Σ 22 xy = [(C 22 13 ○ J -1 xx 22 + C 22 33 ○ J -1 xy 22 ) ○ D 22 X (U 12 x )](t) (29k) + [(C 22 13 ○ J -1 yx 22 + C 22 33 ○ J -1 yy 22 ) ○ D 22 Y (U 21 x )](t) (29l) + [(C 22 23 ○ J -1 xy 22 + C 22 33 ○ J -1 xx 22 ) ○ D 22 X (U 12 y )](t) (29m) + [(C 22 23 ○ J -1 yy 22 + C 22 33 ○ J -1 yx 22 ) ○ D 22 Y (U 21 y )](t) (29n) + M 22 xy (t). (29o)
The stresses are then transformed according to the curvilinear mapping using: In Equations 29-31, most of the computational burden lies in the evaluation of the partial derivatives D kl X () and D kl Y () using Equation 21. We give detailed expressions for the partial derivatives needed to evaluate Equations 29 and 31 in Appendix E.

S 11 xx = |J| 11 ○ [J -1 xx 11 ○ Σ 11 xx + J -
In order to obtain a stable scheme, the coefficients α in Equation 23 must be complementary when evaluating the partial derivatives D kl X () and D kl Y () in Equations 29 and 31. If the coefficients 31. This rather simple strategy simply enforces numerical reciprocity and ensures that the discrete wave equation remains self-adjoint [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF].

[α x ′ -, α x ′ + , α y ′ -, α y ′ + ] are used to compute the partial derivatives D kl X (U i ) and D kl Y (U i ) in Equation 29, then, the complementary coefficients [1 -α x ′ -, 1 -α x ′ + , 1 -α y ′ -, 1 -α y ′ + ] must be used to compute the partial derivatives D kl X (S ij ) and D kl Y (S ij ) in Equation
Notice that more elaborated yet related approaches have been introduced in DGM (e.g., [START_REF] Grote | Interior penalty discontinuous galerkin method for maxwell's equations: Energy norm error estimates[END_REF]).

VERIFICATION

To evaluate the accuracy of the proposed algorithm, we model 2D seismic wave propagation in the isotropic 1D radial ak135 Earth model [START_REF] Kennett | Constraints on seismic velocities in the earth from traveltimes[END_REF] using the DFDM algorithm. While our model is isotropic, the numerical tests demonstrate our algorithm's capability to effectively handle anisotropic media. This is because the coefficient representing the curvilinear mapping and the anisotropic moduli are the same, leading to comparable apparent numerical anisotropy (i.e., when assuming undeformed elements). We do not consider a cylindrical symmetry as in Nissen-Meyer et al.

(2014); nonetheless, such a symmetry can be envisioned and may be relevant for future applications.

We compare the results to those obtained using the SEM method. We present the setups for our DFDM and the SEM benchmark simulations and then discuss our results.

Example Setup

We use a simplified isotropic version of the ak135 Earth model with no attenuation and a single-layer crust. Inside the crustal layer, we use a linear velocity profile. It interpolates the values of the elastic properties of the ak135 model from the surface to the base of the lower crust.

To perform the DFDM simulation, we decompose the model into multiple elements, following a standard cube sphere approach [START_REF] Komatitsch | Introduction to the spectral element method for three-dimensional seismic wave propagation[END_REF][START_REF] Nissen-Meyer | A 2-d spectral-element method for computing sphericalearth seismograms-ii. waves in solid-fluid media[END_REF], as shown in Figure3b. There are 8 layers of elements in the numerical mesh. They match the principal interfaces of the ak135 model, that is, the MOHO at a depth of 35 Km, the transition zone between the upper and the lower mantle at a depth of 660 Km, the core-mantle boundary (CMB) at a depth of 2891.5 Km and the inner core boundary (ICB) at a depth of 5153.5 Km. To achieve a relatively homogeneous number of grid points per minimum wavelength and improve efficiency, we use two sub-layers in the

ρ = Cov[u, u ref ] σ[u]σ[u ref ]
1.000 0.946 0.982 0.995

Table 1. Summary of the parameters employed in the different simulations, the measured errors, and CPU times.

The SEM-Ref simulation is the reference simulation with the highest accuracy (i.e., at least 11 grid points per wavelength). The SEM-1 simulation has a standard setup (with at least 5.5 grid points per wavelength).

The SEM-2 simulation uses a finer mesh than SEM-1 and is more accurate (with at least 6.9 grid points per wavelength). The DFDM simulation uses a coarse mesh (with at least 3 grid points per wavelength); it is more accurate and faster than SEM-1 and SEM-2. cases, the mesh matches the principal discontinuities, i.e., the MOHO, the 670 discontinuity, the CMB, and the ICB. In the SEM simulations, the mesh has two doubling layers, and the number of elements is adjusted as a function of the desired accuracy. The first doubling layer is just under the 670 Km discontinuity, and the second one is close to the middle of the outer core, between the CMB and the ICB. For the DFDM simulation, the number of elements is independent of the desired accuracy. We adjust the accuracy by varying the number of grid points within the elements. Because DFDM allows for non-conformal interfaces, the number of grid points is adjusted independently inside individual elements. The numbers of elements employed for the SEM-Ref,

SEM-1, SEM-2 and DFDM are indicated in Table 1.

lower mantle, the outer core, and the inner core. In the horizontal direction, the mesh consists of 12 slices (i.e. 3 slices connected to each face of the central cube). Thus, the numerical mesh contains 93 elements in total. We set the number of grid points inside the elements so that there are at least 3 grid points per wavelength (using the grid spacing associated with the virtual grid points in Appendix D).

We shall see that this gives an accuracy comparable to a sufficiently accurate SEM simulation. The enhanced accuracy of DFDM is attributed to different factors. First, as opposed to SEM, no approximation is used in the mass matrix, which gives the differential operators a spectral nature. Second, the numerical solution is highly continuous, thanks to the B-spline bases. Because fewer elements are needed, this reduces the number of interfaces at which the numerical wavefield is discontinuous. Last, mixed conditions are employed to share the boundary values between neighboring elements. Because DFDM allows for non-conformal interfaces, the number of grid points in different directions does not need to be the same in neighboring elements. Finally, we assign the model parameters by evaluating the ak135 model at the virtual grid points in Equation D.1.

For the simulations, we consider a seismic event located at the north pole, at an intermediate depth of 50 km, to mitigate the amplitudes of the surface and body waves. We represent the seismic source using the moment tensor point source in Equation D.3 with

M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (32)
. We take the source time function as the Ricker wavelet

s(t) = (1 -2π 2 f 2 0 t 2 ) e -π 2 f 2 0 t 2 (33)
with peak frequency f 0 = 2.5 × 10 -2 Hz and an estimated maximum frequency f max=6.25×10 -2 Hz.

We recorded the wavefield at 20 receivers at the surface of the Earth. The receivers represented by triangles in Figure 4 are equally spaced and buried at a depth of 10 m. Their epicentral distances vary from 9 to 180 degrees.

Comparison with SEM

To assess the accuracy of the DFDM algorithm, we compare the DFDM solution to SEM solutions.

The SEM method is well established and known to produce sufficiently accurate results for the setup considered. It is very close to the analytical when a sufficient number of grid points per wavelength is employed. We ran three SEM simulations, we used the SPECFEM2D current source [START_REF] Komatitsch | SPECFEM2D v7.0.0 [software[END_REF]) published under the GPL3 license. The first SEM simulation we will refer to as SEM-Ref is our reference simulation. It is highly accurate and relies on a fine mesh with at least 11 grid points per wavelength. The second SEM simulation, which we will refer to as SEM-1, employs a relatively standard setup with at least 5.5 grid points per wavelength. The third simulation, which we will refer to as SEM-2, is more accurate than SEM-1 and employs a finer mesh with at least 6.9 grid points per wavelength.

Similarly to the DFDM simulation, the SEM simulations mesh pictured in Figure 3a matches the principal interfaces of the ak135 model. The three SEM meshes employed for the simulations have two doubling layers to obtain a relatively homogeneous number of grid points per minimum wavelength.

The doubling layer is located just under the 670 Km discontinuity. The second one is close to the middle of the outer core, between the CMB and the ICB. We give the total number of elements for SEM-Ref, SEM-1, and SEM-2 in Table 1. They are adjusted to achieve the desired accuracy with a minimum number of grid points per wavelength. Notably, the SEM-1 mesh has a single layer of elements inside the crust, while SEM-Ref and SEM-2 use two layers of elements. Thus, the crust is modeled much more accurately in the SEM-Ref and SEM-2 simulations. In the SEM simulations, we assign the model parameters by evaluating the ak135 model at the grid points. We use the same time step in all simulations. It is small enough to minimize the numerical dispersion associated with the time integration scheme. It is close to the stable time step of reference simulation SEM-Ref, which has the most restrictive stability condition. The stable time step in DFDM depends on the basis functions employed and may be adjusted by varying the spacing between the basis functions toward the sides of the elements. With the setup we employ here, the stability condition is slightly more restrictive for DFDM than for SEM. However, because DFDM requires less grid points per wavelength to achieve a certain accuracy, the stable time step is larger in the DFDM simulations.

Results and discussion

We first compare the results obtained using DFDM to those obtained for the reference solution SEM-Ref, which is assumed to be nearly perfectly accurate. In Figure 4 we present snapshots showing the wavefield obtained for DFDM and SEM-Ref. Even though DFDM uses as few as 3 points per wavelength, we observe no visible difference between the two wavefields, which highlights the high accuracy of DFDM when using coarse meshes [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF][START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF]. Notice the P to S and S to P wave conversions at the solid-fluid interfaces are accurately modeled by DFDM, which uses the same scheme in the entire domain, while SEM uses a specific solid-fluid coupling approach. This is due to the specific choice of the boundary coefficients α in Equation 23, as detailed in Appendix E. Generally, when dealing with interconnected purely solid or fluid elements, the choice of alpha values is relatively arbitrary. However, in the case of solid-fluid coupling interfaces, the value of alpha is fixed. On the one hand, we use the internal value of the displacement in the solid elements (i.e. α = 1 ). On the other hand, we take the internal value of the pressure in the fluid elements (i.e. α = 0). This specific choice enforces the continuity of the pressure and gives the most accurate results.

In Figures 5a and5b we superimpose the seismograms simulated by DFDM and SEM-Ref at the 20 receivers. There is no noticeable difference between the seismograms obtained using DFDM and SEM-Ref, both for the body waves and for the strong dispersive surface waves. This is valid for all epicentral distances, which is remarkable and confirms the high accuracy of DFDM.

To quantify the error, we measured the following quantities for our DFDM simulation:

• The normalized maximum error

ε = Max[ |u -u ref | ] Max[ |u ref | ] (34)
where u and u ref are the vectors containing the displacement components at all times at all receivers for the considered and the reference solutions.

• The standard deviation of the error normalized by that of the reference solution where σ denotes the standard deviation.

σ = σ[u -u ref ] σ[u ref ] (35)
• The Pearson correlation coefficient

ρ = Cov[u -u ref ] σ[u]σ[u ref ] (36)
where Cov denotes the covariance.

We show values of the errors in Table 1. These errors may appear relatively large, but the measured Pearson correlation coefficient ρ = 0.995 indicates an excellent correlation between the DFDM and the SEM-Ref simulations (ρ = 1 corresponds to a perfect correlation). In the DFDM simulation, the maximum difference and its standard deviation are 18% and 10%.

To investigate the relative accuracy of DFDM with respect to SEM when a standard setup is used, we compare the differences between the seismograms obtained using DFDM, SEM-1, SEM-2, and SEM-Ref in Figures 5c and5c. We observe that SEM-1 has the most prominent error, while DFDM has the smallest error. Notice that the error is principally associated with the surface waves that are notoriously difficult to model accurately. The surface waves are primarily sensitive to the velocity in the crust, and we know that DFDM is more accurate than SEM when using a similar number of points per wavelength [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF]. Thus, we anticipate that DFDM is more accurate than SEM-1 because DFDM and SEM-1 employ one layer of elements in the crust with 5 grid points in the vertical direction. What is more remarkable is the superior accuracy of DFDM with respect to SEM-2, because SEM-2 uses two layers of elements inside the crust (which corresponds to 9 grid points in the vertical direction). In contrast, DFDM uses only one layer of elements (with 5 grid points in the vertical direction). Further, DFDM uses only 3.5 points per wavelength in the crust in the horizontal direction, while SEM-2 uses almost 7. These observations show that DFDM is remarkably accurate and adapted to model surface waves. Such waves are notoriously challenging in global and regional seismological modeling.

To better quantify the numerical error, we also measured the error, its standard deviation, and the Pearson correlation coefficient for SEM-1 and SEM-2. We show values of the errors in Table 1. In SEM-1, we use the recommended setup for SEM, i.e., at least 5 points per wavelength [START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF])), and we expect no significant numerical dispersion. We find a maximum error of 33%, and its standard deviation is 19%. Thus the typical error is twice as large in SEM-1 as in DFDM.

The measured Pearson correlation coefficient ρ = 0.946 also indicates a degraded correlation with the reference solution. It is remarkable given that the grid size (i.e., the number of grid points) in SEM-1 is about 3.5 times larger than in DFDM. For SEM-2, where there are at least 6.9 points per wavelength, we find a maximum error of 19%, a standard deviation is 12%, and the Pearson correlation coefficient is ρ = 0.982. These values show that the error in SEM-2 is slightly larger than in DFDM.

Distributional finite-difference modeling 29

We consider, however, that these two simulations have comparable accuracy. Looking at the grid dimensions in DFDM and SEM-2 suggests that SEM requires roughly 5 times more grid points to achieve an accuracy similar to DFDM in 2D. Because DFDM requires more computational effort than SEM to update the wavefield, comparing the CPU time measured for the different simulations is interesting. The values given in Table 1 are only indicative because the specfem2D code used for the SEM simulation is a well-established, and optimized software while our DFDM code is still at an early development stage. Nonetheless, we observe that the CPU time for DFDM is slightly shorter than for SEM-1 and significantly smaller than for SEM-2. Thus DFDM appears substantially faster than SEM to achieve comparable accuracy.

DISCUSSION AND OUTLOOK

We showed that DFDM is significantly more accurate than SEM in numerical simulations of the global wave propagation. Through measurements of the CPU time, we observed that the gain in accuracy achieved with DFDM also translates into better computational efficiency. It is especially striking that we can model the surface waves with fewer points per wavelength than in SEM.

For further applications and to better estimate the net gain in computational efficiency with respect to SEM, our algorithm needs to be expanded from two to three dimensions. This is the subject of a forthcoming paper (Lyu et al., 2023a, in prep.). Lateral variations in the velocity structure and topography will also need to be accounted for in real applications. In this study, we limited our numerical tests to an axisymmetric model without topography for an easier comparison with SEM. This is because material heterogeneity is accounted for using slightly different approximations in SEM and DFDM, and obtaining a trustworthy reference solution in heterogeneous media is delicate. Nonetheless, [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF] demonstrated that DFDM can accurately account for heterogeneity within the elements, including sharp interfaces. Furthermore, [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF] showed that DFDM is in agreement with SEM in the case of strong topographic variations. Accounting for the 3D structure of the Earth's crust is especially challenging in global seismology. Our numerical examples demonstrate that DFDM allows us to employ non-conformal meshes to efficiently handle the strong contrast between the crust and the upper mantle. This is more delicate in SEM, where doubling layers or similar approaches are needed to refine the numerical grid with respect to depth or as a function of velocity fluctuations. Furthermore, DFDM also allows us to use small elements as SEM, and similar strategies may be adopted. Thus, DFDM brings additional flexibility to accurately account for structures with complicated geometries. In practice, a 3D crust may be implemented by varying the upper element thickness to accommodate Moho topography; additional elements may be added to account for sharp lateral variations.

The data underlying this article will be shared on reasonable request to the corresponding author (yder masson@berkeley.edu). The SPECFEM software can be obtained at the following address https: //geodynamics.org/resources/specfem2d and is described in [START_REF] Komatitsch | The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF]. we obtain the basis functions of the B-spline basis β p with polynomial order p recursively using the Cox-de Boor formula

β p i (x ′ ) = x ′ -k i k i+p -k i β p-1 i (x ′ ) + k i+p+1 -x ′ k i+p+1 -k i+1 β p-1 i+1 (x ′ ). (A.3)
In our examples, we employ pairs B 1 /B 2 of B-spline bases where

B 1 def = β p 1 (x ′ ) with (p 1 = p = 4, N 1 = N ) and, B 2 def = β p 2 (x ′ ) with (p 2 = p -1 = 3, N 2 = N -1).
In this case, the basis functions are shifted similarly to the staggered finite-difference grid, as illustrated in Figure A1a.

The two-dimensional canonical B-spline bases are constructed from the one-dimensional bases

B k x (x ′ ) (k = 1, 2
) and B l y (y ′ ) (l = 1, 2) defined in direction x ′ and y ′ using a tensor product. We have

B kl ij (x ′ , y ′ ) = B k x i (x ′ )B l y j (y ′ ). (A.4)
Such two-dimensional canonical basis functions are illustrated in Figure 2a.

APPENDIX B: DEFINITIONS: MASS, STIFFNESS AND BOUNDARY MATRICES

We define the matrices and vectors needed to construct the orthonormal basis functions and the DFD operators. We build those matrices and vectors from an arbitrary pair of canonical bases B 1 /B 2 . In this paper, we employ staggered B-spline bases, as detailed in A and illustrated in Figure A1.

Distributional finite-difference modeling 39

The mass matrix is:

M kl θ def = (M kl θ ) ij = ∫ θ+ θ- B k θi (θ)B l θj (θ)dθ ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ θ def = x ′ or y ′ (k, l) def = (1, 2) or (2, 1) . (B.1)
Here, B k θi (θ) and B l θj (θ) are the 1D canonical basis functions presented in Appendix A. The subscript θ def = x ′ or y ′ denotes the direction of the computational domain x ′ or y ′ along which M θ is computed.

θ def = x ′ or y ′ denotes the coordinates in the computational domain with bounds [θ -, θ + ] def = [x ′ -, x ′ + ]
and [y ′ -, y ′ + ]. The superscripts kl indicate the bases from which M kl θ is computed ( B 1 or B 2 in Appendix A ). When the basis functions have compact support, the mass matrix is band diagonal, ensuring the computations' efficiency. Further, if the two bases are identical (i.e., k = l), then the mass matrix M kk θ is symmetric. It is employed to construct the orthonormal bases introduced in the next section.

The stiffness matrix is The boundary matrix Q θ in eq. (B.2) is:

K kl θ def = (K kl θ ) ij = ∫ θ+ θ- B k θi (θ)B l (1) θj (θ)dθ (B.2a) = -∫ θ+ θ- B k (1) θi (θ)B l θj (θ)dθ + [B k θi (θ)B l θj (θ)] θ+ θ- (B.2b) def = -(K lk θ ) T + Q θ (B.
Q θ = -p k θ ⋅ p l θ T + q k θ ⋅ q l θ T = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ -1 0 +1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ p k/l θ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 ⋮ 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ q k/l θ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 ⋮ 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (B.3)
Notice that Equation B.3 is specific to the B-spline bases considered in this paper because all the basis functions are equal to zero at the extremities of the interval considered, except the first one and the last one. More general expressions are available in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF].

APPENDIX C: ORTHONORMAL BASES: ORTHONORMAL B-SPLINE BASES

To construct the orthonormal bases Bk θ employed to represent the field variables, we start by computing the mass matrices M kk in eq. (B.1) associated with the canonical bases B k (for k = 1, 2). Then, we compute the Cholesky factorizations L k θ of the mass matrices. we have

M kk θ = L k θ ⋅ L k θ T ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ θ def = x ′ or y ′ k def = 1 or 2 . (C.1)
and we define the basis functions Bk x i (x ′ ) of the orthonormal bases Bk x (x ′ ) as the linear combinations x , v

Bk θi (θ) = N θ ∑ j=1 (L k θ ) -1 ij B k θj (θ) for θ = x ′ , y ′ and k = 1, 2 (
ij (x ′ , y ′ ) = Bk x i (x ′ ) Bl y j (y ′ ) (C.4)
where Bk x (x ′ ) (k = 1, 2) and Bl y (y ′ ) (l = 1, 2) are the one-dimensional bases defined in directions x ′ and y ′ .

APPENDIX D: DISCRETIZATION OF THE MATERIAL PARAMETERS USING VIRTUAL GRIDS

To discretize the material parameters, we adopt the approach introduced in [START_REF] Masson | P-SV wave propagation in heterogeneous media: Velocity-stress distributional finite-difference method[END_REF].

We define virtual grids where each basis function is associated with a unique grid point. The grid coordinates are computed from the 1D basis in the equation C.2; we have

xk i = ∫ x ′ + x ′ - Bk xi x ′ dx ′ ∫ x ′ + x ′ - Bk xi dx ′ ŷl j = ∫ y ′ + y ′ - Bl yj y ′ dy ′ ∫ y ′ + y ′ - Bl yj dy ′ . (D.1)
The coordinates (x k i , ŷl j ) define the virtual grids G kl associated with the basis Bkl in each domain. They can be employed to assign the values to the material parameter arrays. In the proposed algorithm, we use the following:

|J| kl def = [|J| kl ] i,j = |J|(x k i , ŷl j ) (D.2a) J -1 ab kl def = [J -1 ab kl ] i,j = J ab (x k i , ŷl j ) = ∂a ′ ∂b (x k i , ŷl j )
where ab = xx, xy, yx, yy (D.2b)

ρ kl def = [ρ kl ] i,j = ρ(x k i , ŷl j ) (D.2c) C kl mn def = [C mn ] i,j = C mn (x k i , ŷl j ) (D.2d)

MOMENT TENSOR POINT SOURCE REPRESENTATION

We represent the seismic source using the moment tensor point source

M(x, t) = [ M ⋅ δ(x -x s )] ⋅ s(t) M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M xx M xy M xy M yy ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (D.3)
where s(t) is the source time function, M is the moment tensor, δ is the delta function and x s is the source position vector. Practically, we implement the seismic source by adding the components of the moment tensor to the stresses in Equation 29. We have

M kl xx (t) = M xx ⋅ δ kl s ⋅ s(t) (D.4) M kl yy (t) = M yy ⋅ δ kl s ⋅ s(t) (D.5) M kl xy (t) = M xy ⋅ δ kl s ⋅ s(t) (D.6)
where we take the delta function as following: 1) the weights are complementary to those employed for the displacement (this is needed to obtain a stable scheme); 2) If the axes of the computational domain aren't oriented similarly in the two domains, the components of the transformed stress need to be rotated.

Equation E.2, the boundary values are computed in all domains Ω using:

[u j θ (x ′ -)] Ω = [U kj T θ ⋅ b k x ]
Ω (θ = x, y; j = 1, 2; k = 3 -j) (E.3a)

[u j θ (x ′ + )] Ω = [U kj T θ ⋅ d k x ]
Ω (θ = x, y; j = 1, 2; k = 3 -j) (E.3b)

[u i θ (y ′ -)] Ω = [b i T x ⋅ U il θ ]
Ω (θ = x, y; i = 1, 2; l = 3 -i) (E.3c)

[u i θ (y ′ + )] Ω = [d i T x ⋅ U il θ ]
Ω (θ = x, y; i = 1, 2; l = 3 -i).

(E.3d)

The workflow for computing and exchanging the boundary values according to Equation E.2 is sketched in Figure A2.

The partial derivatives of the stress components needed to evaluate the right-hand-side of Equa- where the 2D arrays S ij θθ contain the expansion coefficients representing the stresses inside the current domain Ω and the 1D arrays s i θθ () contain the expansion coefficients representing the stresses at the boundaries of the current domain Ω. We compute the boundary values in Equation E.4 using: Finally, when the parametric coordinates don't have the same orientation in the two domains, one needs to rotate the boundary values of the transformed stresses. For example, using

s j θθ (x ′ -) = (1 -α x ′ -) [s j θθ (x ′ -)] Ω + (α x ′ -) T jj y Ω←Ω x ′ - ⋅ [s j θθ (x ′ + )] Ω x ′ - ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ θθ = xx,
⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ [s kl xx ] i,
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⋅ 1 |J| ⋅ J (E.7)
where J is the Jacobian matrix that aligns the parametric coordinates in the two domains.

Following [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF], we set the coefficients α θ ′ ± in Equations E.2 and E.5 depending on the nature of the interfaces. We have:

• α θ ′ ± = 1 at all interfaces surrounding the computational domain to model a free surface. • α θ ′ ± = 1/2 at all interfaces between domains that are both solid/elastic or fluid/acoustic. • α θ ′ ± = 1 inside the solid/elastic domain and α θ ′ ± = 0 inside the fluid/acoustic domain, at all solidfluid interfaces.

Figure 2 .

 2 Figure 2. (a) Representation of the canonical B-spline bases employed to construct the orthonormal bases that represent the field variables in the proposed DFD algorithm. The stresses σ 11 xx , σ 11 yy and σ 11 xy are represented using the basis pictured in gray. The stresses σ 22 xx , σ 22 yy and σ 22 xy are represented using the basis pictured in red. The displacements u 12 x and u 12 y are represented using the basis pictured in blue. The displacements u 21 x and u 21 y are represented using the basis pictured in green. (b) Representation of the analogous finite-difference staggered grid. The complete grid to the left is called the Lebedev grid. It decomposes into the two staggered grids to the right (i.e., as in Virieux (1986)).

  kl contains the internal expansion coefficients representing the function f kl (x ′ , y ′ ) inside the domain Ω. The indices kl indicates the basis B kl in which f kl (x ′ , y ′ ) is represented inside Ω. f j (x ′ ± ) contains the expansion coefficients representing the boundary value f kl (x ′ ± , y ′ ) along the faces separating the domains Ω x± and Ω in Figure 1. The index i indicates the basis B j y in which f kl (x ′ ± , y ′ ) is represented inside Ω. f i (y ′ ± ) contains the expansion coefficients representing the boundary value f kl (x ′ , y ′ ± ) along the faces separating the domains Ω y± and Ω in Figure 1. The index j indicates the basis B i x in which f kl (x ′ , y ′ ± ) is represented inside Ω. D ij θ , b i θ are the 1D-DFD operators introduced in Section 2.5.

  ) Here, the brackets' subscripts indicate the domain from which the boundary values are taken. [] Ω corresponds to the current domain, and [] Ω x ′ -, [] Ω x ′ + , [] Ω y ′ -and [] Ω y ′ + correspond to the neighboring domains to the left, right, bottom and top, respectively, as pictured in Figure 1. The weighted average coefficients α x ′ -, α x ′ + , α y ′ -and α y ′ + are associated with the left, right, bottom and top faces of the domain Ω, respectively. The matrices T nn θ in Equation 23 are associated with the faces. They transform the boundary values taken from the neighboring domains Ω n and return their representation in the basis

Figure 3 .

 3 Figure 3. Numerical mesh employed for the SEM simulations (left) and the DFDM simulation (right). In both

]Figure 4 .Figure 5 Figure 6 .

 456 Figure 4. Snapshots showing the magnitude of the displacement modeled in the SEM-Ref (left) and DFDM (Right) simulations, at times 400 s, 800 s, and 1200 s. The two wavefields look perfectly identical even though the mesh of the DFDM simulation is much coarser than in the SEM-Ref simulation.

Figure A1 .

 A1 Figure A1. Illustration of the basis functions B 1 and B 2 and, their orthogonal counterparts B1 and B2 for p = 5.

  basis functions B θ .

  denotes the elements of the inverse of the matrix L k θ . The basis function in eq. (C.2) are orthonormal in the sense that Mkk θ def = ( Mkk θ ) ij = ∫ θ+ θ-Bk θi (θ) Bk θj (θ)dθ = δ ij def = I, (C.3) where δ ij is the Kronecker Delta and I is the identity matrix. This result is easily verified by inserting Equation C.2 in Equation C.3 and using Equation B.1. Notice that we chose the Cholesky factorization in Equation C.1 for its simplicity and efficiency. However, any factorization of the form M k θ = A k θ ⋅ A k θ T may be employed to construct the proposed algorithm. Equation C.3 shows that the orthonormal bases Bk θ (θ) are their own dual bases, by definition. Finally, we construct the two-dimensional orthonormal B-spline bases employed to represent the discrete field variables (v 21

  kl s ) ij = ∬ ϕ kl ij (x)δ(x -x s ) dx dy = ϕ kl ij (x s ). (D.7)In the examples, the source time function employed is the Ricker wavelet in Equation 33.

Figure

  Figure A2. (a) Workflow for computing and exchanging the boundary values between two neighboring domains Ω 1 and Ω 2 . Here, we consider the case where the right face of domain Ω 1 coincides with the left face of domain Ω 2 . In both domains, we first evaluate the displacement U ij x at the boundary using Equation E.3. We then send the boundary values [u j x (x ′ )] Ω 1 to the neighboring domain. During the exchange, the transformation matrix T ij y is applied to obtain the boundary values in the basis of the receiving domain. Finally, the internal and external boundary values are weighted (according to the weights on the arrows) and summed to obtain the averaged boundary value in Equation E.2. (b) is similar to (a) but for the stresses. The steps are the same except the

  5d) where the internal boundary values [s i θθ ()] Ω are taken from the current domain Ω and the external Distributional finite-difference modeling 45 boundary values [s i θθ ()] Ω θ ′ ± are taken from the neighboring domains. The matrices T ij θ perform basis transformations to represent the external boundary values in the same basis as the displacement inside the current domain Ω. The weight α θ ′± controls whether the boundary values s i θθ () are taken inside (α θ ′ ± = 0), outside (α θ ′ ± = 1)or in between (0 < α θ ′ ± < 1) the current and the neighboring domains. In Equation E.5, the boundary values are computed in all domains Ω using: xx, yy, xy; j = 1, 2; k = 3 -j) xx, yy, xy; j = 1, 2; k = 3 -j) xx, yy, xy; i = 1, 2; l = 3 -i) (E.6c)[s i θθ (y ′ + )] Ω = [d i T x ⋅ S il θθ ]Ω (θ = xx, yy, xy; i = 1, 2; l = 3 -i) (E.6d)

  senting the components of the moment tensor point source M ij in the basis Bkl ,

							1 xy	11	○ Σ 11 xy ](t)	(30a)
	S 22 xy = |J| 22	○ [J -1 yx	22	○ Σ 22 xx + J -1 yy	22	○ Σ 22 xy ](t)	(30b)
	S 22 yx = |J| 22	○ [J -1 xy	22	○ Σ 22 yy + J -1 xx	22	○ Σ 22 xy ](t)	(30c)
	S 11 yy = |J| 11	○ [J -1 yy	11	○ Σ 11 yy + J -1 yx	11	○ Σ 11 xy ](t).	(30d)
	Finally, the displacements are updated to the next time step, at time t + ∆ t , using:
	1 t ∆ 2	[U 21 x (t + ∆ t ) -2U 21 x (t) + U 21 x (t -∆ t )] = [ (D 21 X (S 11 xx ) + D 21 Y (S 22 xy )) ⊘ |J| 21	⊘ ρ 21	](t) (31a)
	1 t ∆ 2	[U 12 y (t + ∆ t ) -2U 12 y (t) + U 12 y (t -∆ t )] = [ (D 12 X (S 22 yx ) + D 12 Y (S 11 yy )) ⊘ |J| 12	⊘ ρ 12	](t). (31b)
	D kl X (F) are 2D arrays with dimensions (N k x × N l y ) containing the values or expansion coefficients repre-
	senting the partial derivative ∂ ∂x D kl Y (F) are 2D arrays with dimensions (N k x × N l y ) containing the values or expansion coefficients repre-
	senting the partial derivative ∂ ∂y C kl ij are 2D arrays with dimensions (N k x × N l y ) containing the values of the parameters C ij (x k i , ŷl j )
	evaluated at the virtual grid points (x k i , ŷl j ) associated with the basis Bkl ,
	J -1 ij	kl are 2D arrays with dimensions (N k x × N l y ) containing the values of the inverse Jacobian matrix's
	component (J -1 ) ij (x k i , ŷl
	|J| kl are 2D arrays with dimensions (N k x ×N l y ) containing the values of the determinant of the Jacobian
	matrix |J|(x k i , ŷl j ) evaluated at the virtual grid points (x k i , ŷl j ) associated with the basis Bkl ,
	∆ t denotes the time step,

Notice that the displacements at the two previous time step needs to be stored in memory while the stresses are only used as a intermediate step to update the displacements. The arrays above are as follow:

Σ kl ij are 2D arrays with dimensions (N k x × N l y )

containing the values or expansion coefficients representing the stresses σ ij (t) in the basis Bkl ,

U kl i are 2D arrays with dimensions (N k x × N l y )

containing the values or expansion coefficients representing the displacements u i (t) in the basis Bkl , ′ (f kl (x ′ , y ′ )) in the basis Bkl , ′ (f kl (x ′ , y ′ )) in the basis Bkl , M kl ij (t) are 2D arrays with dimensions (N k x × N l y ) containing the values or expansion coefficients reprej ) evaluated at the virtual grid points (x k i , ŷl j ) associated with the basis Bkl , ○ denotes the element-wise or the Hadamard product, ⊘ denotes the element-wise or the Hadamard division.

  j [s kl xy ] i,j [s kl xy ] i,j [s kl yy ] i,j

	⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦	←	⎡ ⎢ ⎣ ⎢ ⎢ ⎢ ⎢	[s kl

xx ] i,j [s kl xy ] i,j [s kl xy ] i,j [s kl yy ] i,j

Accounting for attenuation should pose no particular problems and may be modeled using readily available time-domain algorithms relying on the generalized Maxwell body (GMB-EK) or the generalized Zener body. For regional simulations, the computational domain needs to be truncated using absorbing boundaries. This is discussed in a separate manuscript (Lyu et al., 2023b) where we show that the perfectly matched layer method (PML) can be easily adapted to DFDM and is stable under certain conditions that we outline.

In this study, we employed a second-order displacement time integration scheme which makes it possible to a parsimonious strategy that avoids the storage of stress components within the elements [START_REF] Luo | Parsimonious staggered grid finite-differencing of the wave equation[END_REF]. While the proposed algorithm is not strictly parsimonious, It can be made parsimonious by changing the synchronization, i.e., adjusting the moments when the boundary values are exchanged and the boundary conditions are applied.

CONCLUSION

We introduced a staggered-grid distributional finite-difference method (DFDM) for modeling seismic wave propagation through the 2D Earth. We compared the performance of DFDM and the popular spectral element method (SEM) when modeling wave propagation through the same reference Earth model. We emphasize the following characteristics of the proposed approach i) the simplicity of the algorithm, locally similar to the finite-difference method, globally similar to the discontinuous finiteelement method; ii) the versatility of the implementation allowing for the domain decomposition, heterogeneity of the medium, curvilinear mesh, anisotropy, non-conformal interfaces, or discontinuous grid and fluid-solid interfaces; and iii) the efficiency with superior accuracy, reduced memory usage (G ≈ 3), and potentially higher efficiency than SEM. Thus, we conclude that DFDM is a potentially interesting alternative to SEM for global tomography, which heavily relies on sizeable global scale simulations. The authors would like to express their gratitude to Editor Ebru Bozdag and Reviewer Hom Nath Gharti for their constructive comments and remarks, which helped improve our manuscript.

ACKNOWLEDGEMENT

This

APPENDIX A: CANONICAL BASES: STAGGERED B-SPLINE BASES

In order to derive the distributional finite-difference algorithm, one first needs to choose the canonical bases. In this study, we employ B-spline bases (See, e.g. De Boor 1978) that are highly continuous and offer good accuracy. Some alternative bases are given in [START_REF] Masson | Distributional finite-difference modelling of seismic waves[END_REF].

In one dimension, the B-spline basis β p with polynomial order p and basis functions

is constructed as follows. We first define the equally spaced knot vector k with N + p + 1 elements that starts and ends with p + 1 repeated knots

Then, starting with the zero-order B-spline basis β 0 with basis functions

APPENDIX E: COMPUTATION OF THE PARTIAL DERIVATIVES AND COMMUNICATION OF THE BOUNDARY VALUES.

The partial derivatives of the displacement components needed to evaluate the right-hand-side of Equation 29 are:

where the 2D arrays U ij θ contain the expansion coefficients representing the displacements inside the current domain Ω and the 1D arrays u i θ () contain the expansion coefficients representing the displacements at the boundaries of the current domain Ω. We compute the boundary values in Equation E.1 using ± control whether the boundary values u i θ () are taken inside (α θ ′ ± = 1), outside (α θ ′ ± = 0)or in between (0 < α θ ′ ± < 1) the current and the neighboring domains. In