TY - JOUR
AU1 - Koene, Erik F, M
AU2 - Robertsson, Johan O, A
AU3 - Andersson,, Fredrik
AB - SUMMARY We present a method to position point sources at arbitrary locations on finite-difference (FD) grids. We show that implementing point sources on single nodes can cause considerable errors when modelling with the FD method. In contrast, we propose to create a spatially distributed source (over multiple nodes) that nonetheless creates the desired point-source response. The spatial point source is formulated in the wavenumber domain and is based on the FD coefficients used for the wave propagation. Using this ‘FD-consistent source’ on 1-D and 2-D examples, we show that we can obtain superior fits to analytical solutions compared to single-node or sinc-function source implementations, and we show that sources can be offset to arbitrary locations from ‘on’ the grid to ‘off’ the grid, while resulting in solutions that are identical to within machine precision Numerical modelling, Computational seismology, Wave propagation 1 INTRODUCTION The finite-difference (FD) method is one of the most efficient and most widely used seismic modelling methods (Robertsson & Blanch 2020). For example, it plays an essential role in wave-equation-based imaging methods such as reverse-time migration (Baysal et al. 1983) and full-waveform inversion (Tarantola 1984). Within the FD method, derivatives in the wave equation are replaced with FD approximations of limited accuracy. The differences between the continuous wave equation and its discretized counterparts, however, give rise to frequency-dependent phase and amplitude errors. We identify five error sources as a result of discretizing the wave equation. First, only band-limited solutions (in space and time) can be accurately recovered from the discrete grid, constraining the type of solutions that can be computed. Secondly, symmetric FD approximations of derivatives can be made free of amplitude errors, but do introduce phase errors. Such errors can be minimized with, for example, optimized FD operators for the spatial derivatives (Holberg 1987) and filtering procedures for the temporal derivatives (Stork 2013; Koene et al. 2018, 2019). Thirdly, fine heterogeneous models must be carefully sampled onto coarse grids to generate the correct amplitude and phase response (Muir et al. 1992; Hobro 2010). Fourth and fifth on the list are the accurate modelling of point sources and point receivers in space. It is common to represent these by means of sinc functions, such that if the desired point coincides with a node on the grid, energy is added to or recorded from that node only (the oscillatory nature of the sinc function results in zero contributions elsewhere). Conversely, if the point does not coincide with the grid, there is a non-zero injection or recording at all points of the grid. The sinc function may be approximated with a small number of nodes only using, for example, a least-squares approximation (Mittet & Arntsen 1999, 2000) or windowing (Hicks 2002). In this paper, we argue that the sinc function is suboptimal for approximating point sources in FD modelling. We illustrate this first with a simple reproducible example. Consider the 1-D scalar wave equation for pressure p(x, t) at position x and time t, $$\begin{eqnarray*} \frac{1}{c^2}\frac{\partial ^2 p(x,t)}{\partial t^2} - \frac{\partial ^2 p(x,t)}{\partial x^2} = S(t)\delta (x), \end{eqnarray*}$$(1) where c represents a (constant) phase velocity, S(t) is a source-time function and δ(x) is a Dirac delta distribution for a point source in space. We approximate eq. (1) with a second-order accurate FD scheme on a grid of temporal steps Δt and spatial steps Δx, $$\begin{eqnarray*} q_n^{m+1}\!=\!2q_n^m\!-\!q_n^{m-1}\!+\!C^2( q_{n-1}^m\! -\! 2q_n^m\! +\! q_{n+1}^m\! + \!\Delta x^2 S^{m} \delta _n^*), \end{eqnarray*}$$(2) where |$q_n^m=q(n\Delta x,m\Delta t)$| approximates the pressure field, Sm = S(mΔt) samples the source wavelet, |$\delta _n^*=\delta ^*(n\Delta x)$| approximates the spatial source and C = cΔt/Δx is the Courant number. It is well known (e.g. Igel 2016) that for C = 1 and without a source-term present, solutions of eqs (1) and (2) satisfy the same dispersion relation. More precisely, the wavefield p(x, t) is then accurately obtained exactly on the discrete nodes |$q_n^m$|⁠. The FD scheme with C = 1 and without a source is thus highly accurate. Our interest here, however, goes out to the case in which a source is introduced. We will test the C = 1 scheme in combination with two different source excitation methods. First, we simulate a point source in eq. (2) by exciting the source wavelet at a single node coincident with the grid. For example, with a source position xs = 250Δx, that gives $$\begin{eqnarray*} \delta _n^{*,\text{classical}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{\Delta x} & \quad \text{if } n = 250, \\ 0 & \quad \text{otherwise}. \end{array}\right. \end{eqnarray*}$$(3) As time-dependent source functions in this paper, we will use $$\begin{eqnarray*} S(t)\equiv S_{f_0,t_0}(t) = \left(1-2\pi ^2 f_0^2 (t-t_0)^2 \right) e^{-\pi ^2 f_0^2 (t-t_0)^2}, \end{eqnarray*}$$(4) which are Ricker wavelets, implicitly parametrized by the peak frequency f0 and a peak amplitude at time t0. We simulate with Δx = 10 m, c = 2000 m s−1, Δt = 5 ms, such that C = 1. For the source wavelet, we choose t0 = 0.08 s and model four different peak frequencies f0 = {20, 30, 40, 50} Hz. Snapshots of the FD result after 100 time steps, that is, |$q_n^{100}$|⁠, are displayed in Fig. 1(a). For increasing frequencies, we observe a persistent ringing originating from the source location, even after the entire source wavelet has been fed into the simulation. Figure 1. Open in new tabDownload slide Results for 1-D scalar wave simulations, using a Ricker wavelet as a source-time function in (a)–(c) and the derivative of a Ricker wavelet in (d)–(f). Panels (a), (b), (d) and (e) show snapshots at t = 0.5 s for sources excited at the dotted black line for varying frequencies. The snapshots corresponding to different central source frequencies are offset along the vertical axis. Panels (c) and (f) show a comparison between the computed time-series and the expected solutions. Figure 1. Open in new tabDownload slide Results for 1-D scalar wave simulations, using a Ricker wavelet as a source-time function in (a)–(c) and the derivative of a Ricker wavelet in (d)–(f). Panels (a), (b), (d) and (e) show snapshots at t = 0.5 s for sources excited at the dotted black line for varying frequencies. The snapshots corresponding to different central source frequencies are offset along the vertical axis. Panels (c) and (f) show a comparison between the computed time-series and the expected solutions. Now consider placing the source in-between two nodes at xs = 250.5Δx, by injecting the average in two neighbouring nodes, $$\begin{eqnarray*} \delta _n^{*,\text{alternative}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{2\Delta x} & \quad \text{if }n=250\text{ or }n=251, \\ 0 & \quad \text{otherwise}. \end{array}\right. \end{eqnarray*}$$(5) We run the same simulation, except for using t0 = 0.08 + Δx/2/c. Snapshots are displayed in Fig. 1(b) and show no ringing. In Fig. 1(c), we compare solutions for source frequency f = 50 Hz, recorded at x = 3000 m, against the analytical solution $$\begin{eqnarray*} p(3000,t) = \frac{c}{2}\int _0^t S\left(t^{\prime }-\frac{\vert {x_s-3000}\vert }{c}\right) {\rm d}t^{\prime }. \end{eqnarray*}$$(6) Because the coarsely sampled source wavelet is integrated discretely by the FD system, we also create a reference solution using discrete integration, $$\begin{eqnarray*} \tilde{p}(3000,m\Delta t) = \frac{c}{2}\sum _{t=0}^m S\left(t-\frac{\vert {x_s-3000}\vert }{c}\right) \Delta t, \end{eqnarray*}$$(7) which is similarly displayed in Fig. 1(c). We observe that the first point source implementation created significant ringing, while the second forms an excellent match to the reference solution (eq. 7). These observations remain valid for other wavelets also, for example when using a derivative of the Ricker wavelet as source wavelet. In Fig. 1(d), we show snapshots where a source is injected on a single node. In Fig. 1(e), conversely, we use the alternative of eq. (5) to inject the wavelet on two nodes. The single-node excitation now does not produce ringing artefacts, but a clear amplitude difference can be observed in the snapshots when comparing the two methods, particularly for the higher frequency wavelets. In Fig. 1(f), we again compare a trace recorded at x = 3000 m and compare it against the analytical and reference solutions of eqs (6) and (7) (but now with the derivative of a Ricker wavelet as the source-time function). It is clear that the one-node excitation results in amplitudes that are nearly twice as large as desired, whereas we obtain an excellent match against the reference solution when we average the source over just two nodes. The implications of these toy examples are significant. If one designs a highly accurate FD scheme but introduces the source inaccurately, significant errors can still manifest themselves in the results. In the remainder of this paper, we generalize the above idea to model point sources and their spatial derivatives (such as dipole sources) at arbitrary locations on FD grids. 2 THEORY We seek solutions to the scalar wave eq. (1), while being restricted to using the FD method. This has two particular implications: (i) quantities are only defined at discrete space and time coordinates, and (ii) derivatives are replaced with FD approximations. The FD equations thus differ from the continuous equations, and do not necessarily produce the same solutions. For example, the FD result can contain numerical dispersion errors, and we can only recover band-limited solutions from the samples on the grid. In other words, the output from the FD solver may be ‘erroneous’ with respect to the desired solutions. We propose a method to obtain more accurate solutions, by changing the input to the FD solver. 2.1 The band-limited system and the sinc source We define the Fourier transform and its inverse as $$\begin{eqnarray*} \widehat{f}(k) = \int _{-\infty }^\infty f(x) e^{-ikx}{\rm d}x, \end{eqnarray*}$$(8) $$\begin{eqnarray*} f(x) = \frac{1}{2\pi } \int _{-\infty }^\infty \widehat{f}(k) e^{ikx} {\rm d}k. \end{eqnarray*}$$(9) In the wavenumber domain, an nth derivative with respect to space becomes $$\begin{eqnarray*} \widehat{f^{(n)}}(k) = (ik)^n\widehat{f}(k), \end{eqnarray*}$$(10) and a δ(x) function becomes $$\begin{eqnarray*} \widehat{\delta }(k) = 1. \end{eqnarray*}$$(11) Assuming that the solution to the wave equation allows a plane-wave decomposition, and that we only have solutions available at regular steps of Δx in space, then we can at best recover continuous solutions with angular wavenumbers up to Nyquist (kn = π/Δx) from those samples. We define a band-limitation operator |$\mathcal {B}$| of form $$\begin{eqnarray*} \mathcal {B}[p](x)=\frac{1}{2\pi }\int _{-k_n}^{k_n} \left[ \int _{-\infty }^\infty p(x^{\prime }) e^{-ikx^{\prime }}{\rm d}x^{\prime } \right] e^{ikx} {\rm d}k, \end{eqnarray*}$$(12) which combines a standard Fourier transform with a band-limited inverse transform. Then, we can define the band-limited pressure field pB(x, t), as obtained from the total pressure field p(x, t), as $$\begin{eqnarray*} p^B(x,t) \equiv \mathcal {B}[p](x,t), \end{eqnarray*}$$(13) for which it also holds that $$\begin{eqnarray*} \mathcal {B}\left[\frac{\partial ^2 p(x,t)}{\partial x^2}\right] = \frac{\partial ^2 \mathcal {B}[p](x,t)}{\partial x^2} = \frac{\partial ^2 p^B(x,t)}{\partial x^2}. \end{eqnarray*}$$(14) Finally, the projection of the Dirac delta to the set of band-limited functions is $$\begin{eqnarray*} \delta ^B(x) \equiv \mathcal {B}[\delta ](x) = \frac{1}{2\pi }\int _{-k_n}^{k_n} e^{ikx}{\rm d}k =\! \frac{1}{\Delta x} \frac{\sin \left( \frac{\pi x}{\Delta x} \right)}{ \frac{\pi x}{\Delta x} }, \end{eqnarray*}$$(15) which is a sinc function multiplied by the inverse of the grid spacing Δx. Thus, application of |$\mathcal {B}$| to the scalar wave eq. (1) gives $$\begin{eqnarray*} \frac{1}{c^2}\frac{\partial ^2 p^B(x,t)}{\partial t^2} - \frac{\partial ^2 p^B(x,t)}{\partial x^2} = S(t)\delta ^B(x). \end{eqnarray*}$$(16) Solutions to this band-limited wave system (16) are pressure fields containing wavenumbers only up to Nyquist. These are the only solutions that can be accurately recovered from a discretely sampled grid. Hence, even in the absence of any numerical modelling errors, the solutions we can recover with the FD method are those that satisfy eq. (16), not eq. (1). We clarify this point with two examples. Consider injecting a source wavelet onto a single node as in eq. (3), and solving the wave equation with the FD or pseudo-spectral method. During the injection of the source wavelet, it may appear that energy is localized in space when inspecting the solution only on the nodes. But if one interpolates this solution in the Fourier domain (e.g. through zero-padding), non-zero energy will appear between all the nodes. Hence, while exciting all the energy on just a single node in the ‘discrete world’, it is as if one introduced an oscillatory and globally acting source in the ‘continuous world’. In this paper, we focus on spatial band limitations, but similar arguments apply to a temporal discretization. In the FD method, a source wavelet is injected only at regular intervals of Δt. The wavelet fed into the simulation thus consists only of the samples S(mΔt). This act of sampling automatically aliases the injected wavelet (any higher frequency component will, on the temporal grid, appear identical to a lower band-limited frequency). The solution that the FD method computes is thus valid for the injected band-limited wavelet, not the original wavelet S(t). Thus, we can only excite and recover band-limited inputs and outputs. We can sample δB(x) from eq. (15) onto a discrete uniform grid as follows. Assume the following notation for discrete data: p[j] contains points defined at locations x[j] = jΔx for integers j. The band-limited point-source in eq. (15) is then sampled onto a uniform grid with js = xs/Δx, as $$\begin{eqnarray*} \delta ^B[j - j_s] = \frac{1}{\Delta x} \frac{\sin \left( \pi (j-j_s) \right)}{ \pi (j-j_s) }. \end{eqnarray*}$$(17) This point-source notion is plotted in Figs 2(a) and (b). For integer values of js, we obtain the notion of the point-source as used in eq. (3), that excites a wavelet onto a single node only: $$\begin{eqnarray*} \delta ^B[j - j_s] = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{\Delta x} & \quad \text{for }j=j_s\text{ and }j_s\text{ an integer}, \\ 0 & \quad \text{otherwise,} \end{array}\right. \end{eqnarray*}$$(18) which corresponds to the case plotted in Fig. 2(a). Figure 2. Open in new tabDownload slide Point sources can be modelled by sampling a continuous function (plotted in dashed grey lines) onto a grid of integer steps of Δx (marked by coloured circles). To generate sources with an offset to the grid, the chosen function is simply translated to the desired source position, and sampled at the grid nodes again. Panels (a) and (b) correspond to a sinc function, while panels (c) and (d) correspond to the FD-consistent source for a second-order accurate scheme. In the Introduction, we considered cases (a) and (d) to introduce point sources. Figure 2. Open in new tabDownload slide Point sources can be modelled by sampling a continuous function (plotted in dashed grey lines) onto a grid of integer steps of Δx (marked by coloured circles). To generate sources with an offset to the grid, the chosen function is simply translated to the desired source position, and sampled at the grid nodes again. Panels (a) and (b) correspond to a sinc function, while panels (c) and (d) correspond to the FD-consistent source for a second-order accurate scheme. In the Introduction, we considered cases (a) and (d) to introduce point sources. 2.2 The FD-consistent point source When using FD modelling, derivatives are not computed accurately up to the Nyquist wavenumber, while a sinc function is typically still assumed for the point source. In this section, we derive an alternative: a source for the FD system that constructs the point-source response for the band-limited eq. (16). Like the sinc function, the newly proposed sources can also be positioned at arbitrary locations by simple translations in space, which is demonstrated in Figs 2(c) and (d). 2.2.1 The FD system Consider a spatial FD operator $$\begin{eqnarray*} \text{D}^n q(x,t) = \sum _j \alpha _j q(x+j,t), \end{eqnarray*}$$(19) with coefficients αj chosen such that the operator Dn approximates the nth derivative of q(x, t) with respect to space. Substitution of the FD operator in the band-limited wave eq. (16) yields the FD system: $$\begin{eqnarray*} \frac{1}{c^2}\frac{\partial ^2 q(x,t)}{\partial t^2} - \text{D}^2 q(x,t) = S(t)\delta ^\text{FD}(x). \end{eqnarray*}$$(20) The source term δFD(x) is the quantity we seek. We use the variable q rather than p, to emphasize that solutions to eq. (20) correspond to the FD system. The spatial Fourier transform of Dn applied to q(x) is $$\begin{eqnarray*} \widehat{\text{D}^nq}(k) & = \sum _j \alpha _j e^{-ikj}\, \widehat{q}(k) = (i\varphi (k))^n\, \widehat{q}(k), \end{eqnarray*}$$(21) where we expect that the FD operator approximates the nth derivative with respect to space, that is, φ(k) ≈ k. We assume that φ(k) is bijective on the interval [ −kn, kn]. This assumption is satisfied for typical FD operators using Taylor coefficients or optimal FD coefficients. 2.2.2 Useful properties of typical FD operators We highlight the two FD operators that are commonly used in practice, to get some explicit examples of φ(k) in typical applications. First, consider the staggered grid FD operator DFDs that approximates first derivatives with respect to space, defined as $$\begin{eqnarray*} \text{D}_\text{FDs} q(x) = \sum _{l=1}^L \alpha _l \frac{q(x\!+\!(l\!-\!\frac{1}{2})\Delta x)\!-\!q(x\!-\!(l\!-\!\frac{1}{2})\Delta x)}{\Delta x}, \end{eqnarray*}$$(22) where the subscript FDs refers to the staggered grid FD system. Transforming eq. (22) to the wavenumber domain gives $$\begin{eqnarray*} \widehat{\text{D}_\text{FDs} q}(k)= i \sum _{l=1}^L 2\alpha _l\frac{\sin \left( k(l-\frac{1}{2})\Delta x \right)}{\Delta x} \widehat{q}(k). \end{eqnarray*}$$(23) Following the formulation in eq. (21) that defines φ(k), we obtain $$\begin{eqnarray*} \varphi _\text{FDs}(k) = \sum _{l=1}^L2\alpha _l\frac{\sin \left( k(l-\frac{1}{2})\Delta x \right)}{\Delta x}, \end{eqnarray*}$$(24) $$\begin{eqnarray*} {\frac{{\rm d}\varphi _\text{FDs}(k)}{{\rm d}k}} = \sum _{l=1}^L 2\alpha _l (l-\tfrac{1}{2})\cos \left(k(l-\tfrac{1}{2})\Delta x\right). \end{eqnarray*}$$(25) Second, consider the conventional (i.e. non-staggered) grid FD operator DFDc to approximate second derivatives with respect to space: $$\begin{eqnarray*} \text{D}_\text{FDc}^2 f(x) = \alpha _0 \frac{f(x)}{\Delta x^2} + \sum _{l=1}^L \alpha _l \frac{f(x+l\Delta x)+f(x-l\Delta x)}{\Delta x^2}. \end{eqnarray*}$$(26) After transforming eq. (26) to the wavenumber domain, we find $$\begin{eqnarray*} \varphi _\text{FDc}(k) = \frac{{\rm sign}(k)}{i\Delta x}\sqrt{\alpha _0 + \sum _{l=1}^L 2\alpha _l \cos (kl\Delta x)}, \end{eqnarray*}$$(27) $$\begin{eqnarray*} \frac{{\rm d}\varphi _\text{FDc}(k)}{{\rm d}k} = \frac{\sum _{l=1}^L \alpha _l l \sin (k l \Delta x)}{\varphi _\text{FDc}(k)\Delta x}, \end{eqnarray*}$$(28) where the subscript FDc refers to the conventional grid FD system. 2.2.3 Derivation of the FD-consistent point source Similar to |$\mathcal {B}$| in eq. (12), we define a transform operator |$\mathcal {G}$| as $$\begin{eqnarray*} \mathcal {G}\left[ q \right](x) &=& \frac{1}{2\pi }\int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} \left[\int _{-\infty }^\infty q(x^{\prime }) e^{-ikx^{\prime }} {\rm d}x^{\prime } \right] e^{i\varphi (k)x}{\rm d}k, \nonumber \\ &=& \frac{1}{2\pi }\int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} \widehat{q}(k) e^{i\varphi (k)x}{\rm d}k, \end{eqnarray*}$$(29) which transforms q into the angular wavenumber domain with a standard Fourier transform, and then takes an inverse transform using φ(k) rather than k. Limits are taken with an inverse function, that is, φ(φ−1(k)) = k. The limits are taken with inverse functions, to simplify an integration by substitution that takes place later in this section. This operator has the following special property $$\begin{eqnarray*} \frac{\partial ^n \mathcal {G}\left[ q \right](x)}{\partial x^n} &=& \frac{\partial ^n}{\partial x^n}\left[ \frac{1}{2\pi } \int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} \widehat{q}(k) e^{i\varphi (k)x} {\rm d}k \right],\nonumber \\ &=& \frac{1}{2\pi }\int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} (i\varphi (k))^n \widehat{q}(k) e^{i\varphi (k)x}{\rm d}k, \nonumber \\ &=& \frac{1}{2\pi }\int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} \widehat{\text{D}^nq}(k) e^{i\varphi (k)x}{\rm d}k, \nonumber \\ &=& \mathcal {G}\left[\text{D}^nq\right](x), \end{eqnarray*}$$(30) where we used eq. (21). Note the similarity between eqs (14) and (30). In particular, an application of |$\mathcal {G}$| to both sides of the FD wave system in eq. (20) yields (where we momentarily omit the space–time dependencies for brevity) $$\begin{eqnarray*} \mathcal {G}\left[\frac{1}{c^2}\frac{\partial ^2 q}{\partial t^2}\!-\! \text{D}^2 q \right] = \mathcal {G}\left[ S \delta ^\text{FD} \right] \!=\! \frac{1}{c^2}\frac{\partial \mathcal {G}[q]}{\partial t^2} \!-\! \frac{\partial ^2 \mathcal {G}[q]}{\partial x^2} = S\mathcal {G}\left[ \delta ^\text{FD} \right]. \end{eqnarray*}$$(31) We can now make a remarkable observation: eq. (31) has the form of the scalar wave equation. That is, if we make the substitutions |$\mathcal {G}[q](x,t)=p^B(x,t)$| and |$\mathcal {G}[\delta ^\text{FD}](x)=\delta ^B(x)$| in eq. (31), we would obtain eq. (16). We focus on the implication of the second relation, that is, can we realize the operation |$\mathcal {G}[\delta ^\text{FD}](x)=\delta ^B(x)$|? We consider the following ansatz in the wavenumber domain, $$\begin{eqnarray*} \widehat{\delta }^\text{FD}(k) = \frac{{\rm d}\varphi (k)}{{\rm d}k}, \end{eqnarray*}$$(32) because application of |$\mathcal {G}$| yields (using integration by substitution with u = φ(k)) $$\begin{eqnarray*} \mathcal {G}[ \delta ^\text{FD} ](x) &=& \frac{1}{2\pi }\int _{\varphi ^{-1}(-k_n)}^{\varphi ^{-1}(k_n)} e^{i\varphi (k)x} \frac{{\rm d}\varphi (k)}{{\rm d}k}{\rm d}k,\nonumber \\ &=& \frac{1}{2\pi } \int _{-k_n}^{k_n} e^{iux}{\rm d}u= \delta ^B(x). \end{eqnarray*}$$(33) Hence, if we apply |$\mathcal {G}$| to the FD system (20), we obtain $$\begin{eqnarray*} \frac{1}{c^2}\frac{\partial \mathcal {G}[q](x,t)}{\partial t^2} - \frac{\partial ^2 \mathcal {G}[q](x,t)}{\partial x^2} = S(t)\delta ^B(x). \end{eqnarray*}$$(34) Thus, if we use the source expression of eq. (32) in an FD simulation to compute an FD solution q(x, t), then |$\mathcal {G}[q](x,t)$| is a solution to the band-limited scalar wave equation with a band-limited point source. We refer to the source expression in eq. (32) as the FD-consistent point source. We remark that we do not advocate computing |$\mathcal {G}[q](x,t)$|⁠. From inspection, that operation just phase shifts the FD solution in space and, hence, would fix the spatial dispersion error (analogous to the time-dispersion transforms, e.g. Koene et al. 2018, 2019). Instead, we suggest to use FD operators that have little spatial dispersion to begin with. So, we will not compute |$\mathcal {G}[q](x,t)$| and, hence, our solution may still contain spatial dispersion. 2.3 Implementation using the fast Fourier transform 2.3.1 Implementing sinc multipoles We can rapidly approximate δB as defined in eq. (15) with an inverse discrete Fourier transform (IDFT). Let p[j] = p[j + N] be a periodic signal with periodic Fourier series |$\tilde{\text{p}}[l]=\tilde{\text{p}}[l+N]$|⁠, then $$\begin{eqnarray*} \text{p}[j] = \text{IDFT}\big [ \tilde{\text{p}} \big ][j] = \frac{1}{N} \sum _{l=0}^{N-1} \tilde{\text{p}}[l] e^{\frac{i2\pi }{N} l j}, \end{eqnarray*}$$(35) which is rapidly computed with an inverse fast Fourier transform. The variable l in eq. (35) is related to angular wavenumbers as $$\begin{eqnarray*} k_l = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{2\pi }{N\Delta x}l & \quad \text{for } l\lt \frac{N}{2}, \\ \frac{2\pi }{N\Delta x}\left(l-N\right) & \quad \text{for } l\ge \frac{N}{2}. \end{array}\right. \end{eqnarray*}$$(36) We thus approximate and sample δB from eq. (15) with $$\begin{eqnarray*} \text{d}_\text{BL,0}[j - j_s] = \frac{1}{\Delta x}\text{IDFT}\Big [ e^{-ik_lj_s\Delta x} \Big ][j]. \end{eqnarray*}$$(37) But whereas eq. (17) samples the sinc function as-is (Fig. 3a), the use of the IDFT in eq. (37) returns its periodic summation with period N (Fig. 3b), $$\begin{eqnarray*} \text{d}_\text{BL,0}[j - j_s] = \frac{1}{\Delta x}\sum _{h=-\infty }^\infty \frac{ \sin (\pi (j-j_s+hN))}{\pi (j-j_s+hN)}. \end{eqnarray*}$$(38) The discrepancy between eqs (17) and (38) is a good example of the subtle differences between continuous and discrete operations, which is exactly what is of such importance in this paper. Figure 3. Open in new tabDownload slide Sampling a sinc function as a point source on the half-open interval [0, 20Δx). (a) The sinc function (dashed black line) is sampled onto the grid. (b) The periodic summation of the sinc (red line), made up from an infinitely periodic train of sinc pulses (soft red line), is sampled onto the grid. Figure 3. Open in new tabDownload slide Sampling a sinc function as a point source on the half-open interval [0, 20Δx). (a) The sinc function (dashed black line) is sampled onto the grid. (b) The periodic summation of the sinc (red line), made up from an infinitely periodic train of sinc pulses (soft red line), is sampled onto the grid. We generalize the principle to obtain n-th order multipoles such as monopoles (n = 0) and dipoles (n = 1), by applying the derivative operation of eq. (10) to eq. (37), $$\begin{eqnarray*} \text{d}_{\text{BL},n}[j-j_s] = \frac{1}{\Delta x}\text{IDFT}\Big [ (ik_l)^n e^{-ik_lj_s\Delta x} \Big ][j]. \end{eqnarray*}$$(39) 2.3.2 Implementing FD-consistent multipoles In Appendices A and B, we show two cases where it is possible to compute an analytical expression for the FD-consistent point source. In Appendix C, we show that it is also possible to approximate the point source through a least-squares optimization. But we can also obtain its periodic summation with period N (just like eq. 38) through an IDFT as follows. The FD-consistent point source at the origin, defined in eq. (32) and as computed with an IDFT, is denoted by dFD, 0(x) and is computed as $$\begin{eqnarray*} \text{d}_{\text{FD},0}[j] = \frac{1}{\Delta x}\text{IDFT}\left[ \frac{{\rm d}\varphi (k_l)}{{\rm d}k} \right][j]. \end{eqnarray*}$$(40) To shift the point source to other positions, we simply translate the solution over the grid in the wavenumber domain as before, $$\begin{eqnarray*} \text{d}_{\text{FD},0}[j-j_s] = \frac{1}{\Delta x}\text{IDFT}\left[ e^{-ik_l j_s\Delta x} \frac{{\rm d}\varphi (k_l)}{{\rm d}k} \right][j]. \end{eqnarray*}$$(41) We obtain the nth multipole source by applying the derivative operation of eq. (21) to the point source in eq. (41), $$\begin{eqnarray*} \text{d}_{\text{FD},n}[j-j_s] = \frac{1}{\Delta x}\text{IDFT}\left[ (i\varphi (k_l))^n e^{-ik_l j_s\Delta x} \frac{{\rm d}\varphi (k_l)}{{\rm d}k} \right][j]. \end{eqnarray*}$$(42) As an example, a MATLAB (2017) script that solves eq. (42) for the staggered grid FD method can be written as follows: phi = 2*sin(k*([1:L]-1/2)*dx)/dx * (a); J = 2*cos(k*([1:L]-1/2)*dx) * (a.*([1:L]’-1/2)); dFD = 1/dx*real(ifft((1i*phi).^ n.*exp(-1i*k*xs).*J)); In this script, k is the angular wavenumber array of eq. (36), L the number of FD coefficients, dx the step size, a the array of FD coefficients, n the derivative order, xs the point-source location and dFD corresponds to the FD-consistent point source. 2.4 Extension to higher dimensions The derivation so far has focused on the 1-D case. The extension to higher dimensions is, however, straightforward. The 2-D extension of the (periodic) sinc source of eq. (39) is given as the product of two separable point sources (e.g. Hicks 2002). Along each coordinate axis, j and h for example, we use $$\begin{eqnarray*} \text{d}_{\text{BL},nr}[j-j_s,h-h_s] = \text{d}_{\text{BL},n}[j-j_s]\text{d}_{\text{BL},r}[h-h_s]. \end{eqnarray*}$$(43) The 2-D extension of the FD-consistent point source is easily derived. It can be verified that we can use the following transform operator |$\mathcal {G}$| for the 2-D case, $$\begin{eqnarray*} \mathcal {G}[q](x,z) = \frac{\iint \widehat{q}(k_x,k_z) e^{i\left[\varphi (k_x)x+\varphi (k_z)z\right]}{\rm d}k_x {\rm d}k_z}{(2\pi )^2} , \end{eqnarray*}$$(44) where kx and kz are the horizontal and vertical angular wavenumbers, respectively. The integration limits in eq. (44) should be between ±φ−1( · ) for the horizontal and vertical Nyquist wavenumbers. For horizontal and vertical FD operators |$D_x^2$| and |$D_z^2$|⁠, the operator, for example, gives $$\begin{eqnarray*} \mathcal {G}[D_x^2 q+D_z^2q](x,z,t) = \frac{\partial ^2 \mathcal {G}[q](x,z,t)}{\partial x^2} + \frac{\partial ^2 \mathcal {G}[q](x,z,t)}{\partial z^2}. \end{eqnarray*}$$(45) We can carry out the same steps as before in comparing the FD system with the continuous wave equation. As ansatz we then choose the following wavenumber representation of the source, $$\begin{eqnarray*} \widehat{\delta }^\text{FD}(k_x,k_z) = \frac{{\rm d}\varphi (k_x)}{{\rm d}k_x}\frac{{\rm d}\varphi (k_z)}{{\rm d}k_z}. \end{eqnarray*}$$(46) This gives |$\mathcal {G}[d](x,z) = \delta ^B(x,z)$|⁠, analogously to eq. (33). A 2-D inverse Fourier transform of eq. (46) then gives the 2-D point source in space. Generalizing the concept to multipoles and using an inverse fast Fourier transform, following eq. (42), we then find $$\begin{eqnarray*} \text{d}_{\text{FD},nr}[j-j_s,h-h_s] = \text{d}_{\text{FD},n}[j-j_s]\text{d}_{\text{FD},r}[h-h_s]. \end{eqnarray*}$$(47) Hence, the 2-D FD-consistent multipole point source, eq. (47), is given as the product of two separable point sources, just as the 2-D sinc multipole point source. 2.5 The FD-consistent point source as a wavenumber filter 2.5.1 The example from the introduction In eq. (2) in the Introduction, we approximated the scalar wave eq. (1) using the FD method as $$\begin{eqnarray*} \frac{p(x,t-\Delta t)\!-\!2p(x,t)\!+\!p(x,t+\Delta t)}{\Delta t^2} \! =\! c^2 \text{D}_\text{FDs}\text{D}_\text{FDs} p(x,t), \end{eqnarray*}$$(48) where we use the staggered first derivative operator DFDs of eq. (22) using just a single FD coefficient of α1 = 1. That is, DFDsp(x, t) = [p(x + Δx/2, t) − p(x − Δx/2, t)]/Δx, which is then followed by another application of DFDs. The FD-consistent point source at xs = 250.5Δx, using eq. (25) for the expression dφFDs(k)/dk = cos (kΔx/2), is then given by $$\begin{eqnarray*} \text{d}_{\text{FD},0}[j-250.5] &=& \frac{1}{\Delta x}\text{IDFT} \left[ e^{-250.5ik_l\Delta x} \cos (\tfrac{k_l\Delta x}{2}) \right]\![j],\nonumber \\ &=& \frac{1}{\Delta x}\text{IDFT} \left[ \frac{e^{-250ik_l\Delta x} + e^{-251ik_l\Delta x} }{2} \right]\![j],\nonumber \\ &=& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{1}{2\Delta x} & \quad \text{if }j=250\text{ or } j=251,\\ 0 & \quad \text{otherwise},\end{array}\right. \end{eqnarray*}$$(49) which corresponds to the source used in the Introduction (eq. 5). The only difference between the sinc and FD-consistent point source, thus, is the inclusion of the factor cos (kΔx/2) in the inverse Fourier transform operation. In Fig. 4(a), we plot the spectrum of a sinc point source as compared to the FD-consistent point source. This gives a heuristic understanding of the FD-consistent point source: it is an amplitude-related filter in the wavenumber domain. In the Introduction, we saw that the sinc point source (i.e. the single-node excitation) resulted in too large amplitudes compared to the reference solution, particularly for high-frequency waves. The FD-consistent point source, with its non-flat spectrum, created better results. We now see that the proposed source modifies the source energy such that the effectively created wavefield has the appropriate amplitudes for all wavenumbers. In Appendix A, we show that this generalization holds for any staggered grid FD model of half-order L: an FD-consistent monopole half-way between the grid can be represented with just 2L non-zero weights. In Appendix B, we show that a similar generalization holds for the conventional grid, where a dipole positioned exactly on the grid can be represented with just 2L non-zero weights. Figure 4. Open in new tabDownload slide A plot of the spectra due to using a sinc point source versus that of the FD-consistent point source. The horizontal interval goes from 0 up to the Nyquist wavenumber. In (a), we consider a second-order accurate FD-consistent point source, as in the Introduction. In (b), we consider a high-order optimal scheme using 15 coefficients from Liu (2014), see the text. Figure 4. Open in new tabDownload slide A plot of the spectra due to using a sinc point source versus that of the FD-consistent point source. The horizontal interval goes from 0 up to the Nyquist wavenumber. In (a), we consider a second-order accurate FD-consistent point source, as in the Introduction. In (b), we consider a high-order optimal scheme using 15 coefficients from Liu (2014), see the text. 2.5.2 High-order optimal FD Consider another example, using high-order FD. Assume that we solve the same staggered grid scheme as in eq. (48) but for operator half-length L = 15, using optimized FD coefficients proposed by Liu (2014) with a maximum phase velocity error of 1 per cent. We provide the list of coefficients with the online Supporting Information. We can compute the expression dφFDs(k)/dk from eq. (25), which represents the spectrum of the FD-consistent point source. The result is shown in Fig. 4(b). We see that this point-source implementation has a spectrum that oscillates around the spectrum of the sinc point source. As a result, this FD-consistent point source will boost energy at certain wavenumbers, and suppress energy at others. 3 EXAMPLES We now compare the accuracy in FD simulations, using the periodic sinc point source of eq. (39) and the FD-consistent point source of eq. (42). In all examples, zero initial conditions and periodic boundaries are assumed. We provide MATLAB scripts to recreate the examples in the online Supporting Information. 3.1 1-D examples 3.1.1 Staggered grid, monopole source Consider solving the staggered grid scheme in eq. (48) for half-length L = 15, using optimized FD coefficients proposed by Liu (2014) with a maximum phase velocity error of 1 per cent. The spectrum of the corresponding FD-consistent point source is given in Fig. 4(b). We use the time-dispersion transforms from Koene et al. (2018) to eliminate the effect of time dispersion. We assess the results using the following two accuracy measures taken from Mittet (2017): $$\begin{eqnarray*} \text{spectral ratio}(\omega ) = \frac{\vert \text{analytical solution}(\omega )\vert }{\vert \text{FD solution}(\omega ) \vert }, \end{eqnarray*}$$(50) $$\begin{eqnarray*} \text{phase error}(\omega ) = \frac{\arg (\text{FD solution}(\omega ))-\arg (\text{analytical}(\omega ))}{\omega }, \end{eqnarray*}$$(51) for angular frequency ω. As noted in the theory section, we expect that the simulation contains spatial dispersion (thus has a phase error) due to the inaccuracy of the FD operator. However, we expect a spectral ratio of 1 when comparing the computed solution to the analytical solution in eq. (6). We choose Δx = 10 m, Δt = 3 ms and c = 2000 m s−1. We use an f0 = 28 Hz Ricker wavelet. A monopole source is implemented using the sinc source as well as the FD-consistent source formulations. We position the source on a node at xs = 250Δx and record the pressure field on the grid at xr = 280Δx. Note that, for the sinc source, that corresponds to a single-node excitation as in Fig. 2(a). In a second experiment, we shift the source to xs = 250.5Δx. For the sinc and FD-consistent source, that corresponds to exciting the source using potentially all grid nodes. We consider the results in Fig. 5. We show separate traces for the sinc and the FD-consistent source, compared against the analytical solution. On close inspection, it is visible that both sets of traces contain an advanced and a trailing artefact compared to the main lobe. These artefacts are related to spatial dispersion of optimal FD operators; this observation was also made in, for example, Koene & Robertsson (2020). The two obtained results appear very similar otherwise, and appear to fit the analytical solution well. However, when using the sinc source, the spectral ratio shows oscillating artefacts around unity (where the number ‘one’ indicates perfect accuracy). Conversely, the results computed with the FD-consistent source give an accurate spectral ratio over the entire frequency range considered, irrespective of the source position (‘on’ or ‘off’ the grid). The phase error is identical regardless of the source implementation, and is caused by the numerical dispersion of the spatial FD operator in the FD updates, as expected. We see that the source position (whether ‘on’ a node or not) does not influence the spectral ratio or phase error much. Recalling the spectrum of this FD-consistent point source in Fig. 4(b), it is now apparent that the FD-consistent point source has acted as a filter such that the source excited the correct amplitude at all wavenumbers. Figure 5. Open in new tabDownload slide Comparison of a sinc monopole and an FD-consistent monopole in high-order 1-D FD simulations on staggered grids. The left-hand column considers a point source coincident with the grid (a single-node excitation), while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. Figure 5. Open in new tabDownload slide Comparison of a sinc monopole and an FD-consistent monopole in high-order 1-D FD simulations on staggered grids. The left-hand column considers a point source coincident with the grid (a single-node excitation), while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. 3.1.2 Conventional grid, dipole source To illustrate that the methodology applies equally well to non-staggered grids and multipole sources, we now model the following FD system on a conventional (non-staggered) grid: $$\begin{eqnarray*} \frac{p(x,t-\Delta t)-2p(x,t)+p(x,t+\Delta t)}{\Delta t^2} = c^2 \text{D}^2_\text{FDc} p(x,t), \end{eqnarray*}$$(52) with |$\text{D}^2_\text{FDc}$| the FD operator as defined in eq. (26), using optimized FD coefficients proposed by Liu (2013) for half-length L = 15 and a maximum phase velocity error of 0.3 per cent. The list of coefficients is provided with the online Supporting Information. We use the time-dispersion transforms from Koene et al. (2018). We use the same model parameters as the previous example, but model dipole point sources by using n = 1 in eqs (39) and (42). The analytical solution is $$\begin{eqnarray*} p(3000,t) = \frac{{\rm sign}( x_s - 3000 )}{2} S_{f,t_0}\left(t-\frac{\vert {x_s-3000}\vert }{c}\right). \end{eqnarray*}$$(53) The result is shown in Fig. 6. Again, the sinc source yields a solution that has a spectral ratio with artefacts. Conversely, the FD-consistent source results in the same amplitudes as the analytical solution, irrespective of the source location. Figure 6. Open in new tabDownload slide Comparison of a sinc dipole and an FD-consistent dipole in high-order 1-D FD simulations on non-staggered grids. The left-hand column considers a point source coincident with the grid, while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. Figure 6. Open in new tabDownload slide Comparison of a sinc dipole and an FD-consistent dipole in high-order 1-D FD simulations on non-staggered grids. The left-hand column considers a point source coincident with the grid, while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. 3.2 2-D staggered grid examples Having demonstrated that the method produces correct time-series in 1-D, we move on to homogeneous 2-D elastic examples. We model the elastic isotropic stress-velocity system with the staggered grid leap-frog formulation of Virieux (1986). We discretize the model with Δx = Δz = 10 m, Δt = 0.5 ms, vp = 2000 m s−1 and vs = 1200 m s−1, we use optimal high-order coefficients from Liu (2014), with L = 10 for a maximum phase velocity error of 0.1 per cent using the staggered-grid FD operator given in eq. (22). We use periodic boundaries and time-dispersion transforms as before. We create an explosive (P-wave only) source by adding a curl-free force density vector onto the particle velocity components (Hicks 2002). In the 2-D case, that corresponds the following particle velocity updates: $$\begin{eqnarray*} \rho {\begin{pmatrix}\frac{\partial v_x(x,z,t)}{\partial t} \\ \frac{\partial v_z(x,z,t)}{\partial t} \end{pmatrix}} = {\begin{pmatrix}\frac{\partial \sigma _{xx}}{\partial x} + \frac{\partial \sigma _{xz}}{\partial z} \\ \frac{\partial \sigma _{zz}}{\partial z} + \frac{\partial \sigma _{xz}}{\partial x} \end{pmatrix}} + {\begin{pmatrix}\delta (z)\frac{\partial \delta (x)}{\partial x} \\ \delta (x)\frac{\partial \delta (z)}{\partial z} \end{pmatrix}}S(t). \end{eqnarray*}$$(54) For the source-wavelet S(t), we use a 33 Hz integrated Ricker wavelet. We position the source at (x, z) = (2255, 2255) m. On the Virieux grid, that source coincides with the location where the shear stress σxz is positioned. 3.2.1 Checking snapshots using ‘global’ sources As a first-order accuracy check, we will see if we can generate P-waves only, because an erroneous source injection could generate undesired S-wave modes. Hence, we create the 2-D source injection array following eqs (43) and (47). Both methods will generate a source using potentially all nodes on the grid. Fig. 7 shows two snapshots using the sinc and FD-consistent point source each. We note that at the early time, t = 0.0375, the sinc source excitation has caused energy to appear at an entire strip of z locations. Conversely, the FD-consistent point source is more local in space. A similar effect could be seen in Fig. 2, which shows that the FD-consistent point source can be inherently more localized in space compared to the ringy sinc point-source. In Fig. 7(b), we observe that the sinc point source has generated ringing that trails the P wave, and a close look even reveals a faint presence of S waves. Conversely, in Fig. 7(d) for the FD-consistent source, ringing is much reduced and linked to just spatial dispersion. Figure 7. Open in new tabDownload slide Snapshots of vz from modelling an explosive source in 2-D elastic simulations. Panels (a) and (b) use the sinc source, cf. eq. (43), while panels (c) and (d) use the FD-consistent point source, cf. eq. (47). Clipping levels are identical for snapshots of the same time level. Figure 7. Open in new tabDownload slide Snapshots of vz from modelling an explosive source in 2-D elastic simulations. Panels (a) and (b) use the sinc source, cf. eq. (43), while panels (c) and (d) use the FD-consistent point source, cf. eq. (47). Clipping levels are identical for snapshots of the same time level. 3.2.2 Checking snapshots using ‘localized’ sources In Fig. 7(a), we see that the source excites the entire vz field in the z-direction, which is consistent with the idea of a sinc source (in eq. 16 derived to be the appropriate source for band-limited solutions). However, we may prefer a point source that is more localized in space. The typical strategy to achieve this is to approximate the point source with only a few nodes. We can follow the method of Hicks (2002) to model an explosive source as in eq. (54), using a 30 × 30 node grid surrounding the desired point-source location (i.e. 15 points above and 15 points below the point-source, and similarly on the left- and right-hand sides). Given a P-wave velocity of 2000 m s−1 and a maximum source frequency of about 2.5 × 33 = 82.5 Hz, we obtain that waves of at least k = 2πfmax/vp ≈ 0.26 rad m−1 will be present. With a Nyquist wavenumber of kn = π/Δx ≈ 0.31 rad m−1, we thus require that our source attains a bandwidth of at least 82 per cent of the Nyquist wavenumber. Therefore, we run four band-limited experiments: two using Hicks’ method for the sinc point source and two using our least-squares approximation for the FD-consistent point source as given in Appendix C. We approximate the point sources up to 80 per cent and 90 per cent of the Nyquist wavenumbers. For the method of Hicks (2002), that means using b = 10.12 for the dipole source and b = 10.05 for the monopole source for the first case, and b = 5.04 for the dipole and b = 4.98 for the monopole source in the second case. Fig. 8 shows snapshots for the simulations with band-limited point-source implementations. We see that the point sources created with the method of Hicks (2002) are more ‘ringy’ compared to those created with the FD-consistent point source. Moreover, Hicks’ method creates slowly propagating S-wave modes that should not be present for the point sources that cover 90 per cent of the Nyquist range (Fig. 8b). Conversely, the FD-consistent point source creates only P waves, as desired. The differences between Figs 7(d) and 8(b) and (d) are absolutely minimal. Hence, it appears that the FD-consistent point source can be well-approximated using the least-squares method. Figure 8. Open in new tabDownload slide Snapshots of vz from modelling an explosive source in 2-D elastic simulations using a spatially limited source array, covering 80 per cent (top) and 90 per cent (bottom) of the Nyquist bandwidth. Panels (a) and (b) are computed with coefficients due to the method of Hicks (2002). Panels (c) and (d) are computed with the coefficients of an L2 optimization of the FD-consistent point source per eq. (C2). The clipping in the figures is identical to that for the t = 0.75 s cases in Fig. 7. Figure 8. Open in new tabDownload slide Snapshots of vz from modelling an explosive source in 2-D elastic simulations using a spatially limited source array, covering 80 per cent (top) and 90 per cent (bottom) of the Nyquist bandwidth. Panels (a) and (b) are computed with coefficients due to the method of Hicks (2002). Panels (c) and (d) are computed with the coefficients of an L2 optimization of the FD-consistent point source per eq. (C2). The clipping in the figures is identical to that for the t = 0.75 s cases in Fig. 7. 3.2.3 Comparison against an analytical solution We place an artificial vz receiver at (x, z) = (2255, 2655), that is, exactly 400 m below the source injection point. The receiver is implemented by sinc interpolation, as the vz node is not coincident with the desired receiver position. We use the time-dispersion transforms of Koene et al. (2018) to remove the effect of time dispersion. We then compare the result with an ‘analytical’ reference using the program EX2DELEL (Berg et al. 1994), that uses the Cagniard-de Hoop method for a simple elastic model. The two modelling cases considered are the same as those shown in Figs 8(b) and (d), that is, the ‘localized’ point sources optimized to 90 per cent of the Nyquist wavenumber. Fig. 9 shows the recorded traces and the comparison against the analytical solution. What we can see is that although the traces appear similar, the spectral ratio for Hicks’ (2002) solution is considerably worse than the spectral ratio for the FD-consistent point source which is nearly flat. Conversely, both solutions have nearly identical phase errors. The results quantify what we could already see, which is that the FD-consistent point source used in Fig. 8(d) emulates the desired wavefield better than the source used to generate Fig. 8(b). Figure 9. Open in new tabDownload slide A comparison of FD simulations against an analytical solution from EX2DELEL (Berg et al. 1994). Hicks’ simulation corresponds to that seen in Fig. 8(b), whereas the L2 FD-consistent simulation corresponds to that seen in Fig. 8(d). Figure 9. Open in new tabDownload slide A comparison of FD simulations against an analytical solution from EX2DELEL (Berg et al. 1994). Hicks’ simulation corresponds to that seen in Fig. 8(b), whereas the L2 FD-consistent simulation corresponds to that seen in Fig. 8(d). 3.2.4 Flexibility of source formalism Previously, we created an explosive source by adding the gradient of a 2-D monopole onto the particle velocity updates, per eq. (54). We can alternatively create the explosive source by adding a monopole point source onto both normal stress components, corresponding to the following normal stress updates: $$\begin{eqnarray*} {\begin{bmatrix}\frac{\partial \sigma _{xx}(x,z,t)}{\partial t} \\ \frac{\partial \sigma _{zz}(x,z,t)}{\partial t} \end{bmatrix}} = {\begin{bmatrix}\lambda +2\mu & \quad \lambda \\ \lambda & \quad \lambda + 2\mu \end{bmatrix}}{\begin{bmatrix}\frac{\partial v_x(x,z,t)}{\partial x} \\ \frac{\partial v_z(x,z,t)}{\partial z} \end{bmatrix}} +\delta (x,z) S^{\prime }(t). \end{eqnarray*}$$(55) In this case, we require not the wavelet S(t) but its time derivative S′(t) to achieve an identical source as with eq. (54). An interesting question, then, is whether the results using the explosive sources of eqs (54) and (55) yield the same results. We run that test using the ‘full’ sinc and FD-consistent sources as were used in Fig. 7, not their localized approximations. Otherwise, all settings are kept identical: Δx = 10 m, Δt = 0.5 ms, vp = 2000 m s−1 and vs = 1200 m s−1, optimal coefficients from Liu (2014) with L = 10 for a maximum phase velocity error of 0.1 per cent, a 33 Hz integrated Ricker wavelet for S(t), a source at (x, z) = (2255, 2255) m and a vz receiver at (x, z) = (2255, 2655) m. The derivative of the wavelet in eq. (55) is required at Δt/2 offset due to the leap-frog time-integration scheme. Given the source wavelet S(t) that has been filtered using the forward time-dispersion transform (e.g. Koene et al. 2018), we simply obtain the derivative wavelet at the appropriate times using discrete differentiation, $$\begin{eqnarray*} S^{\prime }\left(t+\frac{\Delta t}{2}\right) = \frac{S(t+\Delta t)-S(t)}{\Delta t}. \end{eqnarray*}$$(56) The results are presented in Fig. 10. Although the two different sinc sources appear to yield identical traces in Fig. 10(a), the zoomed section in Fig. 10(b) shows that there are differences between the two. Also visible in Figs 10(a) and (b) is the presence of an artefact at 0.15 s. The origin of the artefact lies in the fact that the sinc point source is strongly present at all grid nodes simultaneously – thus, the source excitation leads to a small but measurable signal immediately at the receiver location. Furthermore, in Figs 10(c) and (d), we see that the difference between the two traces is on the order of 0.1 per cent of the amplitude of the original energy. Figure 10. Open in new tabDownload slide Recordings of vertical particle velocity vz in units of μm s−1, due to 2-D elastic FD simulations using two formally equivalent definitions of explosive sources. Panels (a)–(d) use the sinc source formalism. Panels (e)–(h) use the FD-consistent source formalism. Figure 10. Open in new tabDownload slide Recordings of vertical particle velocity vz in units of μm s−1, due to 2-D elastic FD simulations using two formally equivalent definitions of explosive sources. Panels (a)–(d) use the sinc source formalism. Panels (e)–(h) use the FD-consistent source formalism. Conversely, as demonstrated in Fig. 10(g), the traces due to the FD-consistent point source are identical to within machine precision (<10−13 differences). Note furthermore that when we compare the zoomed plots (Figs 10b and f), we see that ringing is reduced for the FD-consistent point source, which is in line with all previously obtained results. The remaining ringing is likely due to spatial dispersion. This ringing is visible because we zoom-in on the results to an unusually high degree. So, although the formalism of eqs (54) and (55) are identical in theory, they do not yield identical results when implemented with the sinc point-source strategy. Conversely, with the FD-consistent point source, both methods yield results that are identical to within machine precision. Hence, the proposed point source is consistent with the used FD operators. 3.2.5 Flexibility of source positioning Finally, we carry out a test where we excite an explosive source at (x, z) = (2250, 2250) using a monopole excitation onto the σxx and σzz components, following eq. (55). As this location is coincident with the location of σxx and σzz on the Virieux grid, this means that the sinc point source excites the explosive source with just a single-node excitation. We then measure the σxx component at (x, z) = (2550, 2550), that is, at an offset of 300 m in both the x- and z-directions, which is again exactly coincident with the grid. We then carry out a second model run where the source is positioned at (x, z) = (2255, 2255) m, and a receiver is positioned at (x, z) = (2555, 2555) m which records the stress using sinc interpolation. That is, the source–receiver offset remains constant, but both source and receiver are offset by 5 m (i.e. half a grid step). We again use the ‘full’ sinc and FD-consistent source as in Fig. 7, not their band-limited approximations. The resulting recorded traces are shown in Fig. 11. Although the sinc source (with source and receiver ‘on’ as well as ‘off’ the grid) appears to result in very similar traces, there are differences between the two that become more obvious once we zoom in, as in Fig. 11(b), or take the difference of the two traces as done in Figs 11(c) and (d). Note the artefact around t = 0.15, whose origin lies again in the fact that the sinc point source can be strongly present at all grid nodes simultaneously. This ringing is not recorded when we model the sinc source ‘on’ the grid. Note, however, the difference in ringing amplitude between Figs 11(b) and (f). The FD-consistent point source generated less ringing, which is in line with all previously obtained results. The FD-consistent point source identically shows a similar artefact around t = 0.15 s, visible in Fig. 11(f), albeit smaller in size. It is striking, moreover, that the FD-consistent point source again leads to results identical to within machine precision. Figure 11. Open in new tabDownload slide Recordings of horizontal stress σxx in units of μPa, due to 2-D elastic FD simulations using two identical source–receiver offsets. Panels (a)–(d) use the sinc source formalism. Panels (e)–(h) use the FD-consistent source formalism. Figure 11. Open in new tabDownload slide Recordings of horizontal stress σxx in units of μPa, due to 2-D elastic FD simulations using two identical source–receiver offsets. Panels (a)–(d) use the sinc source formalism. Panels (e)–(h) use the FD-consistent source formalism. 4 DISCUSSION AND CONCLUSION The typical assumption in FD wave modelling is that a sinc function sampled onto the FD grid accurately represents a point source. However, in this paper, we show that the sinc point source produces suboptimal results and can cause errors in FD computations. We derived an alternative implementation of point sources by relating the source term of the band-limited wave equation to the FD wave system. The proposed ‘FD-consistent point source’ is a function of the FD coefficients used in the wavefield propagation, and may be obtained through a filtering procedure in the wavenumber domain. The point-source spectrum is altered with a factor dφ(k)/dk, where φ(k) represents the effective wavenumber modification due to the spatial FD approximations. Through several examples we have shown that our proposed method outperforms the sinc function in obtaining a good fit between FD and analytical solutions. We have implemented two formally equivalent point-source descriptions for elastic waves and obtained solutions identical to within machine precision using our method. Additionally, we have shown that our methodology makes it possible to position a source at an arbitrary location ‘on’ or ‘off’ the grid, while also resulting in solutions that are identical to within machine precision. We have deliberately restricted the examples to homogeneous models. As discussed in the Introduction, heterogeneous models can cause inaccuracies in the FD computation that are unrelated to point-source implementations; thus, it would be more difficult to untangle the various sources of errors. However, this will be a subject for future research. Another question open for future research is what the implications are for source–receiver reciprocity in FD modelling, if receivers are accurately modelled with the sinc formalism, but sources are accurately modelled by the proposed modification. We hypothesize that the reciprocity relation is the following: if a wavefield is recorded with the sinc function, it must be injected with the FD-consistent point source, and vice versa. This conjectured reciprocity relation could, for example, explain the high accuracy attained in the method of multiple point sources (MPS) for wavefield reconstruction in FD modelling (see e.g. Vasmel & Robertsson 2016). The method of MPS uses both sinc and FD-consistent methods to position receivers and sources on FD grids; and their form changes when the FD coefficients are changed in exactly the way that we derive in this paper. The theory of this paper can thus for example be used to justify the form of the receivers and sources that are required to accurately implement the method of MPS. SUPPORTING INFORMATION Online_supplement.zip Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS This work was supported by SNF grant 2-77220-15. We thank Dirk-Jan van Manen for interesting discussions. We thank Thomas Bohlen, an anonymous reviewer, and editors Ivan Vasconcelos and Joerg Renner for helpful comments that have significantly improved this manuscript. REFERENCES Baysal E. , Kosloff D.D., Sherwood J.W., 1983 . Reverse time migration , Geophysics , 48 ( 11 ), 1514 – 1524 . 10.1190/1.1441434 Google Scholar Crossref Search ADS WorldCat Crossref Berg P. , If F., Nielsen P., Skovgaard O., Helbig K., 1994 . Analytical reference solutions , Model. Earth Oil Explor. , 77 , 421 – 427 . OpenURL Placeholder Text WorldCat Hicks G.J. , 2002 . Arbitrary source and receiver positioning in finite-difference schemes using kaiser windowed sinc functions , Geophysics , 67 ( 1 ), 156 – 165 . 10.1190/1.1451454 Google Scholar Crossref Search ADS WorldCat Crossref Hobro J. , 2010 . Rapid and accurate finite-difference model generation from discontinuous anisotropic velocity models , in SEG Technical Program Expanded Abstracts 2010 , pp. 2961 – 2965 ., Society of Exploration Geophysicists . 10.1190/1.3513461 Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Crossref Holberg O. , 1987 . Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena , Geophys. Prospect. , 35 ( 6 ), 629 – 655 . 10.1111/j.1365-2478.1987.tb00841.x Google Scholar Crossref Search ADS WorldCat Crossref Igel H. , 2016 . Computational Seismology: A Practical Introduction , Oxford University Press . 10.1093/acprof:oso/9780198717409.001.0001 Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Crossref Koene E. , Wittsten J., Robertsson J., Andersson F., 2019 . Eliminating time dispersion from visco-elastic simulations with memory variables , in 81st EAGE Conference and Exhibition 2019 , Vol. 2019 , pp. 1 – 5 ., European Association of Geoscientists & Engineers . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Koene E.F. , Robertsson J.O., 2020 . Optimal finite-difference operators for arbitrarily sampled data , Geophysics , 85 ( 3 ), 1 – 54 . 10.1190/geo2019-0081.1 Google Scholar Crossref Search ADS WorldCat Crossref Koene E.F.M. , Robertsson J.O.A., Broggini F., Andersson F., 2018 . Eliminating time dispersion from seismic wave modeling , Geophys. J. Int. , 213 ( 1 ), 169 – 180 . 10.1093/gji/ggx563 Google Scholar Crossref Search ADS WorldCat Crossref Liu Y. , 2013 . Globally optimal finite-difference schemes based on least squares , Geophysics , 78 ( 4 ), T113 – T132 . 10.1190/geo2012-0480.1 Google Scholar Crossref Search ADS WorldCat Crossref Liu Y. , 2014 . Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling , Geophys. J. Int. , 197 ( 2 ), 1033 – 1047 . 10.1093/gji/ggu032 Google Scholar Crossref Search ADS WorldCat Crossref MATLAB , 2017 . version 9.2.0 (R2017a) , The MathWorks Inc , Natick, MA . Mittet R. , 2017 . On the internal interfaces in finite-difference schemes , Geophysics , 82 ( 4 ), T159 – T182 . 10.1190/geo2016-0477.1 Google Scholar Crossref Search ADS WorldCat Crossref Mittet R. , Arntsen B., 1999 . General source and receiver positions in coarse-grid finite-difference schemes , in SEG Technical Program Expanded Abstracts 1999 , pp. 1855 – 1858 ., Society of Exploration Geophysicists . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Mittet R. , Arntsen B., 2000 . General source and receiver positions in coarse-grid finite-difference schemes , J. Seism. Explor. , 9 , 73 – 92 . OpenURL Placeholder Text WorldCat Muir F. , Dellinger J., Etgen J., Nichols D., 1992 . Modeling elastic fields across irregular boundaries , Geophysics , 57 ( 9 ), 1189 – 1193 . 10.1190/1.1443332 Google Scholar Crossref Search ADS WorldCat Crossref Robertsson J.O.A. , Blanch J.O.A., 2020 . Numerical methods, finite-difference , in Encyclopedia of Solid Earth Geophysics , 3rd edn. Gupta H.K., ed. Springer Science Business Media . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Stork C. , 2013 . Eliminating nearly all dispersion error from FD modeling and RTM with minimal cost increase , in 75th Conference and Exhibition , SEG, Extended Abstracts . Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Tarantola A. , 1984 . Inversion of seismic reflection data in the acoustic approximation , Geophysics , 49 ( 8 ), 1259 – 1266 . 10.1190/1.1441754 Google Scholar Crossref Search ADS WorldCat Crossref Vasmel M. , Robertsson J.O.A., 2016 . Exact wavefield reconstruction on finite-difference grids with minimal memory requirements , Geophysics , 81 ( 6 ), T303 – T309 . 10.1190/geo2016-0060.1 Google Scholar Crossref Search ADS WorldCat Crossref Virieux J. , 1986 . P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method , Geophysics , 51 ( 4 ), 889 – 901 . 10.1190/1.1442147 Google Scholar Crossref Search ADS WorldCat Crossref APPENDIX A: ANALYTICAL FD-CONSISTENT POINT SOURCE FOR STAGGERED GRIDS We mostly resort to using IDFTs as a rapid and flexible method to compute the FD-consistent point sources. In some cases, however, it is possible to compute the band-limited inverse Fourier transform analytically. The inverse transform of eq. (32), with dφFDs/dk from (25), can for example be expressed as $$\begin{eqnarray*} \delta ^\text{FDs}(x) &=& \frac{1}{2\pi }\int _{-\frac{\pi }{\Delta x}}^{\frac{\pi }{\Delta x}} \frac{{\rm d}\varphi _\text{FDs}(k)}{{\rm d}k} e^{ikx} {\rm d}k,\nonumber \\ &=& \sum _{l=1}^L \frac{\alpha _l\left( l - \frac{1}{2}\right)}{\Delta x} \Bigg \lbrace {\rm sinc}_\pi \left[\frac{x}{\Delta x}+\left( l - \frac{1}{2}\right)\right] + {\rm sinc}_\pi \left[\frac{x}{\Delta x}-\left( l - \frac{1}{2}\right)\right] \Bigg \rbrace , \end{eqnarray*}$$(A1) where sincπ(x) = sin (πx)/(πx). Hence, a particular simplification occurs for a source half-way between two nodes. Using a form similar to (18), we find that eq. (A1) for a monopole half-way between grid nodes is sampled onto integer nodes j as $$\begin{eqnarray*} \delta ^\text{FDs}[j-j_s] = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{\alpha _l\left(l-\tfrac{1}{2}\right)}{\Delta x} & \text{with}\ \vert {j-j_s}\vert =\left(l-\frac{1}{2}\right), j_s+\frac{1}{2}\text{ integer}, \\ 0 & \text{otherwise.} \end{array}\right. \end{eqnarray*}$$(A2) We can give two examples: In the second-order accurate FD scheme used in the Introduction with L = 1, we used α1 = 1, and put a source at 250.5Δx. Following eq. (A2), that corresponds to the point-source weights δFDs[250 − 250.5] = δFDs[251 − 250.5] = α1/(2Δx) = 1/(2Δx) and zero elsewhere. This was equally shown in eq. (49). If we use a fourth-order accurate FD scheme of L = 2 with coefficients α1 = 9/8 and α2 = −1/24, the FD-consistent source has weights δFDs[249 − 250.5] = δFDs[252 − 250.5] = 3α2/(2Δx) = −1/(16Δx) and δFDs[250 − 250.5] = δFDs[251 − 250.5] = α1/(2Δx) = 9/(16Δx). The realization is thus that for staggered FD grids, the FD-consistent monopole positioned half-way between the grid contains only 2L non-zero coefficients that are easily obtained from the used FD coefficients. We remark that the FD-consistent dipole for staggered grids is not obtained by taking the analytical derivative of eq. (A1). As clear from eq. (42), we must apply the FD derivative operator to compute the FD-consistent dipole from the monopole expression. Thus, the FD-consistent dipole for staggered grids is obtained as $$\begin{eqnarray*} \delta ^\text{FDs}_\text{dipole}(x)=\text{D}_\text{FDs}\delta ^\text{FDs}(x), \end{eqnarray*}$$(A3) with DFDs from eq. (22) and δFDs(x) from eq. (A1). APPENDIX B: ANALYTICAL FD-CONSISTENT DIPOLE FUNCTION FOR CONVENTIONAL GRIDS In Appendix A, we compute a band-limited inverse Fourier transform to obtain an expression for the continuous FD-consistent monopole on staggered grids. We may want to carry out a similar procedure for the FD-consistent point source for the conventional grid. However, the inverse Fourier transform for the monopole does not allow a simple closed-form solution. The FD-consistent expression for a dipole, however, can be obtained with less effort. Hence, we take a band-limited inverse Fourier transform to write the FD-consistent dipole for the conventional grid as $$\begin{eqnarray*} \delta ^\text{FDc}_\text{dipole}(x) &=& \frac{1}{2\pi }\int _{-\frac{\pi }{\Delta x}}^{\frac{\pi }{\Delta x}} \frac{{\rm d}\varphi _{\text{FDc}}(k)}{{\rm d}k} i\varphi _\text{FDc}(k) e^{ikx} {\rm d}k \nonumber \\ &=& \sum _{l=1}^L \frac{\alpha _l l}{2\Delta x^2} \Bigg \lbrace {\rm sinc}_\pi \left[\frac{x}{\Delta x}+l\right] - {\rm sinc}_\pi \left[\frac{x}{\Delta x}-l\right] \Bigg \rbrace , \end{eqnarray*}$$(B1) with again sincπ(x) = sin (πx)/(πx). A particular simplification thus occurs for a dipole centred on the nodes. Using a form similar to eq. (18), we find that eq. (B1) for a dipole on the grid nodes may be sampled onto integer nodes j as $$\begin{eqnarray*} \delta ^\text{FDc}_\text{dipole}[j-j_s]=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}-\frac{\alpha _l l}{2\Delta x^2} & \quad \text{with }\ j-j_s=l, j_s\ \text{integer},\\ 0 & \quad \text{otherwise}.\end{array}\right. \end{eqnarray*}$$(B2) We can give two examples: Consider that we compute second derivatives of the wave equation using |$\text{D}_\text{FDc}^2$| with a second-order (L = 1) accurate scheme, then α0 = −2 and α1 = 1. If we want to position the FD-consistent dipole on x = 250Δx, then the coefficients are |$\delta ^\text{FDc}_\text{dipole}[249-250]=1/(2\Delta x^2)$|⁠, |$\delta ^\text{FDc}_\text{dipole}[250-250]=0$| and |$\delta ^\text{FDc}_\text{dipole}[251-250]=-1/(2\Delta x^2)$|⁠, and zero elsewhere. Consider that we compute second derivatives of the wave equation using |$\text{D}_\text{FDc}^2$| with a fourth-order (L = 2) accurate scheme, then α0 = −5/2, α1 = 4/3 and α2 = −1/12. If we want to position the FD-consistent dipole on x = 250Δx, then the coefficients of the source are |$\delta ^\text{FDc}_\text{dipole}[248-250]=-1/(12\Delta x^2)$|⁠, |$\delta ^\text{FDc}_\text{dipole}[249-250]=2/(3\Delta x^2)$|⁠, |$\delta ^\text{FDc}_\text{dipole}[250-250]=0$|⁠, |$\delta ^\text{FDc}_\text{dipole}[251-250]=-2/(3\Delta x^2)$| and |$\delta ^\text{FDc}_\text{dipole}[252-250]=1/(12\Delta x^2)$|⁠. The realization is thus that for conventional FD grids, the FD-consistent dipole positioned exactly on the grid requires just 2L non-zero coefficients that can be easily obtained from the used FD coefficients. APPENDIX C: A LEAST-SQUARES APPROXIMATION Similar to the sinc function, the FD-consistent source proposed in eq. (42) can be non-zero at every node in the simulation. This is, in principle, correct as we seek solutions to the band-limited scalar wave eq. (16) with such a sinc source. As a result, energy due to the point-source is (correctly) introduced at the source but is also immediately present at the receiver. These artefacts are for example visible in Fig. 11. Such results may not be desired. Furthermore, when computing the sinc source or FD-consistent source with an IDFT, the computed weights for the source wrap around the boundaries. This second property is unappealing when we want to model on non-periodic domains. We can overcome both these limitations by obtaining a least-squares approximation of the source array with a small ‘footprint’ and no wrap-around. That is, like in Mittet & Arntsen (1999, 2000) and Hicks (2002), we limit the spatial extent of the region that generates the point-source. As an example, consider using 2M + 1 points around an arbitrarily positioned point-source. Denoting x0 as a point on the grid closest to the desired source position, we create the following error function for the nth order sinc multipole point source, $$\begin{eqnarray*} E^\text{BL}(\beta ,K_c) &=& \int _{-K_c}^{K_c} \Bigg \vert \sum _{m=-M}^M \Bigg (\beta _m e^{-ik\overbrace{(x_0-x_s+m\Delta x)}^\text{offset between nodes and true source}}\Bigg ) \nonumber \\ && - \underbrace{\frac{1}{\Delta x}(ik)^n }_\text{sinc point-source}\Bigg \vert ^2 {\rm d}k, \end{eqnarray*}$$(C1) and similarly for the FD-consistent point source, $$\begin{eqnarray*} E^\text{FD}(\beta ,K_c) &=& \int _{-K_c}^{K_c} \Bigg \vert \sum _{m=-M}^M \Bigg (\beta _m e^{-ik\overbrace{(x_0-x_s+l\Delta x)}^\text{offset between nodes and true source}}\Bigg ) \nonumber \\ && - \underbrace{\frac{1}{\Delta x}(i\varphi (k))^n \frac{{\rm d}\varphi (k)}{{\rm d}k}}_\text{FD-consistent point-source}\Bigg \vert ^2 {\rm d}k. \end{eqnarray*}$$(C2) The critical wavenumber Kc must be chosen as equal to or less than the Nyquist wavenumber, that is, |$K_c \le \frac{\pi }{\Delta x}$|⁠. For a small spatial range given by L, it is often not possible to fit the filter target up to the Nyquist wavenumber with good accuracy (e.g. with errors of less than 1 per cent for every wavenumber considered). Hence, the value of M should not be chosen too small. We then minimize the system of eq. (C1) or eq. (C2) in a least-squares sense to obtain the best coefficients βm, which correspond directly to the desired weights in the space domain. We have included a MATLAB implementation of this principle in the online Supporting Information to this paper. As an example, we revisit the example of Section 3.1.1, where we modelled a monopole with a high-order staggered grid FD operator. We keep everything the same, except that we will now approximate the sources with 2M + 1 = 31 points [using M = 15 in eqs (C1) and (C2), equaling the half-order L of the FD operator]. We choose to approximate the point-source up to Kc = 0.8π/Δx, that is, up to 80 per cent of the Nyquist wavenumber. As an example output, we show the space-domain versions of two point-source approximations coincident with the grid in Fig. C1. The sinc source, as expected from eq. (18), is non-zero only for x = xs where it takes on the value 1/Δx. The FD-consistent source is significantly more oscillatory over the entire interval. Figure C1. Open in new tabDownload slide A space domain display of a least-squares approximation (abbreviated to lsqr in the legend) of the two point source notions discussed in this paper using 31 points, for a point source coincident with the grid at xs = 2500 m. The point-source approximation only takes on a value at every circular mark, the line is connected between points only to make it easier to visually distinguish the two data sets. Figure C1. Open in new tabDownload slide A space domain display of a least-squares approximation (abbreviated to lsqr in the legend) of the two point source notions discussed in this paper using 31 points, for a point source coincident with the grid at xs = 2500 m. The point-source approximation only takes on a value at every circular mark, the line is connected between points only to make it easier to visually distinguish the two data sets. Fig. C2 shows the results obtained with the least-squares approximation of both the point-source methods. The results are clearly in agreement with the results shown in Fig. 5. In essence, the time-series obtained with the band-limited FD-consistent source obtains the correct spectral ratio, whereas that time-series due to the least-squares approximation of the sinc function does not. The phase error is identical in both cases. Figure C2. Open in new tabDownload slide Comparison of least-squares approximations of a sinc monopole and an FD-consistent monopole in high-order 1-D FD simulations on staggered grids. The left-hand column considers a point-source coincident with the grid, while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. Figure C2. Open in new tabDownload slide Comparison of least-squares approximations of a sinc monopole and an FD-consistent monopole in high-order 1-D FD simulations on staggered grids. The left-hand column considers a point-source coincident with the grid, while the right-hand column considers a point-source half-way between two nodes. The middle and bottom plots compare the obtained accuracy against the analytical solution. The least-squares method provides a simple way to localize point sources with only a small number of nodes, and requires no periodicity of the domain. Hence, the least-squares method, as proposed here, is a simple way to apply the theory of this paper to bounded domains with only a limited number of nodes. Note that the method of images found in Hicks (2002) can furthermore be used if the point-source array partially crosses a free surface boundary. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A consistent implementation of point sources on finite-difference grids
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa383
DA - 2020-11-01
UR - https://www.deepdyve.com/lp/oxford-university-press/a-consistent-implementation-of-point-sources-on-finite-difference-c9WObkrlHl
SP - 1144
EP - 1161
VL - 223
IS - 2
DP - DeepDyve
ER -