Discretization of multipole sources in a finite difference setting for wave propagation problems

Seismic sources are commonly idealized as point-sources due to their small spatial extent relative to seismic wavelengths. The acoustic isotropic point-radiator is inadequate as a model of seismic wave generation for seismic sources that are known to exhibit directivity. Therefore, accurate modeling of seismic wavefields must include source representations generating anisotropic radiation patterns. Such seismic sources can be modeled as multipoles, that is, a time-dependent linear combination of spatial derivatives of the spatial delta function. Since the solutions of linear hyperbolic systems with point-source right hand sides are necessarily singular, standard results on convergence of grid-based numerical methods (finite difference or finite element) do not imply convergence of numerical solutions. We present a method for discretizing multipole sources in a finite difference setting, an extension of the moment matching conditions developed for the Dirac delta function in other applications, along with numerical evidence demonstrating the accuracy of these approximations. Using this analysis, we develop a weak convergence theory for the discretization of a family of symmetric hyperbolic systems of first-order partial differential equations, with singular source terms, solved via staggered-grid finite difference methods: we show that grid-independent space-time averages of the numerical solutions converge to the same averages of the continuum solution, and provide an estimate for the error in terms of moment matching and truncation error conditions. Numerical experiments confirm this result, but also suggest a stronger one: optimal convergence rates appear to be achievedpoint-wise in space away from the source. (C) 2019 Elsevier Inc. All rights reserved.