On Godunov-type finite volume methods for seismic wave propagation

Abstract

The computational complexity of simulating seismic waves demands continual exploration of more efficient numerical methods. While Finite Volume methods are widely acclaimed for tackling general nonlinear hyperbolic (wave) problems, their application in realistic seismic wave simulation remains uncommon, with rare investigations in the literature. Furthermore, seismic wavefields are influenced by sharp subsurface interfaces frequently encountered in realistic models, which could, in principle, be adequately solved with Finite Volume methods. In this study, we delved into two Finite Volume (FV) methods to assess their efficacy and competitiveness in seismic wave simulations, compared to traditional Finite Difference schemes. We investigated Gudunov-type FV methods: an upwind method called wave propagation algorithm (WPA), and a Central-Upwind type method (CUp). Our numerical analysis uncovered that these finite volume methods could provide less dispersion (albeit increased dissipation) compared to finite differences for seismic problems characterized by velocity profiles with abrupt transitions in the velocity. However, when applied to more realistic seismic models, finite volume methods yielded unfavorable outcomes compared to finite difference methods, the latter offering lower computational costs and higher accuracy. This highlights that despite the potential advantages of finite volume methods, such as their conservative nature and aptitude for accurately capturing shock waves in specific contexts, our results indicate that they are only advantageous for seismic simulations when unrealistic abrupt transitions are present in the velocity models.

1 Introduction

Seismic wave propagation, a cornerstone of computational seismology, involves numerically solving elastic wave propagation problems. This process, often referred to as the direct problem in seismic imaging, lies at the core of various applications within computational seismology. These applications span a broad spectrum, including geophysical exploration, regional wave propagation analysis, global or planetary seismology studies, investigations into ground motion and dynamic rupture phenomena, seismic topography (utilizing techniques like full waveform inversion), among others (Igel 2017).

Different mathematical models have been developed to simulate wave propagation in the complex geological structure of the earth, using various numerical techniques, such as finite differences, finite elements, finite volumes, spectral or boundary element methods, discontinuous Galerkin method, and hybrid methods, which have been constantly improved to obtain an accurate and reliable solution of the wave equation for various media (Igel 2017). Still, due to its low computational cost and simplicity, the finite difference method remains the most widely used method. However, although the finite difference method is currently found in large applications in the analysis of seismic wave propagation, in complex earth structures, including those with large velocity contrasts, strong heterogeneity, relief, or topographic attenuation, the methodology can fail (Kalyani et al. 2014). This is mainly due to the assumption often used in formulating finite difference methods that both the wave propagation medium and the solution function of the problem are smooth functions. For problems with discontinuities, the accuracy of finite difference methods deteriorates and may fail to represent propagation through regions of abrupt transitions.

The finite volume method offers a great advantage by generally not requiring smoothness of the propagation medium or the solution, and is a natural choice for modeling in heterogeneous media. In this type of method, each mesh cell can be given different material properties, the propagation of information between cells is done by analyzing the physics imposed on the cell boundaries, and is suitable and popular for hyperbolic conservation laws and equilibrium law systems; See, for example, Godlewski and Raviart (1996), Kröner (1997), Hesthaven (2017), Kurganov (2016), Kurganov (2018),Toro (2013), LeVeque (2002b).

In this numerical study, we investigated two Godunov-type finite volume methods, based on the REA algorithm (Reconstruct - Evolve - Average), to explore their potential in modeling seismic waves across heterogeneous media.

The first method, known as the Wave Propagation Algorithm (WPA), is an Upwind type, originally developed by LeVeque (1997). This method leverages the resulting waves from each Riemann solution in conjunction with limiter flux functions to achieve sharp resolution of discontinuities without inducing numerical oscillations. It exhibits second-order accuracy where the solution remains smooth and has demonstrated efficacy in numerically solving various wave propagation problems, for example, advection equation, linear acoustics equations (in homogeneous and heterogeneous media), Burger’s equation, Euler equations for gas dynamics, shallow water equations, and even the elasticity equations, as well as examples of applications in tsunami simulation.

The second method is a Central-Upwind (CUp) type method pioneered by Kurganov and Lin (2007), building upon the one-dimensional fully and semi-discrete central-upwind schemes introduced in earlier work (Kurganov et al. 2001). This method aims to reduce numerical dissipation, providing a more accurate solution of the evolved quantities on the original grid. The scheme’s simplicity and generality make it an attractive alternative for solving wave propagation problems.

Its applications can be seen, for example: in Kurganov and Lin (2007) where the scheme is applied to solve Euler equations; in Kurganov (2018) a central-upwind method for shallow water equations is developed and, most recently, in Kurganov et al. (2021) an improvement on the numerical dissipation of Central-Upwind schemes with implementation on a wide variety of complex numerical examples.

In the existing literature, very few different finite volume methods have been applied to seismic wave propagation or full waveform inversion problems. Some notable articles address the problem for unstructured grids. For instance, in Dumbser et al. (2007), they present the development of arbitrary high-order finite volume methods on unstructured meshes for the seismic wave propagation problem in 2D and 3D settings. Glinsky-Olivier et al. (2006) proposed a comprehensive reanalysis of the finite-volume approach based on unstructured triangular meshes. More recently, Wang et al. (2022) presented a cell-centered finite volume scheme for the diffusive–viscous wave equation on general polygonal meshes. Additionally, Brossier et al. (2008) developed a finite volume method in the frequency domain for 2D problems. A certain common knowledge that Finite Volume methods are expected to be more expensive than Finite Differences for regular cartesian grids is sometimes heard, but rarely documented.

Alongside, while numerical methods for solving elasticity equations, often in their first-order form as coupled hyperbolic systems, have been documented in the literature using techniques such as WPA (LeVeque 2002a; c) or modified Central-Upwind methods (Kurganov and Pollack 2011), there is a notable absence of literature that evaluates their performance when applied to seismic problems. This research gap motivated our numerical study, which provides valuable insights into the effectiveness and limitations of these numerical methods. We analyze their performance when applied to seismic problems and compare them with the finite difference methods that predominate in the literature.

The paper is structured as follows: In Section 2, we introduce the equations governing the problem under study. Section 3 provides a detailed exposition of the finite volume methods implemented. Finally, in Section 4, we present the numerical results, including their application to synthetic and realistic seismic problems.

2 Modeling equations

In this section, we present the equations governing our study problem (for further details see, e.g. Toro (2009) or LeVeque (2002b)).

2.1 2-D Elastic wave equations

To model compression waves (P-Waves) in 2-D one can use the following system

$$\begin{split}\varepsilon_{t}-u_{x}-v_{y}&=0,\\ (\rho(x,y)u)_{t}-\sigma_{x}(K(x,y);\varepsilon)&=0,\\ (\rho(x,y)u)_{t}-\sigma_{y}(K(x,y);\varepsilon)&=0,\end{split}$$ (1)

where, for brevity, the partial derivatives are denoted by subscripts. Here $u$ and $v$ are the velocities in the $x$ and $y$ directions, respectively, $\varepsilon$ is the strain, $\sigma$ is the stress, and $\rho$ is the density.

The system 1 can be written in the form of a conservation law

$$\textbf{q}_{t}+f(\textbf{q},x,y)_{x}+g(\textbf{q},x,y)_{y}=0,$$ (2)

subject to the initial condition, $\textbf{q}(x,0)=\textbf{q}_{0}(x)$. Here

$$\textbf{q}=\left(\begin{array}[]{c}\varepsilon\\ m_{u}\\ m_{v}\end{array}\right),\quad f(\textbf{q})=\left(\begin{array}[]{c}m_{u}/\rho\\ -\sigma\\ 0\end{array}\right),\quad g(\textbf{q})=\left(\begin{array}[]{c}m_{v}/\rho\\ 0\\ -\sigma\end{array}\right),$$ (3)

where $m_{u}=\rho u$ and $m_{v}=\rho v$ denote the corresponding momenta.

For small strains, we can assume a linear stress-strain relationship of the form

$$\sigma(\varepsilon,x,y)=K(x,y)\varepsilon,$$

where $K(x,y)$ is the volumetric modulus of compressibility. In this case, the linear hyperbolic system $\textbf{q}_{t}+(\textbf{A}(x,y)\textbf{q})_{x}+(\textbf{B}(x,y)\textbf{q})_{y}=0$ has coefficient matrices

$$A(x,y)=\left[\begin{array}[]{ccc}0&-1/\rho(x,y)&0\\ -K(x,y)&0&0\\ 0&0&0\end{array}\right],B(x,y)=\left[\begin{array}[]{ccc}0&0&-1/\rho(x,y)\\ 0&0&0\\ -K(x,y)&0&0\end{array}\right].$$ (4)

This system is a simplification such that only P-waves are modeled, not S-waves (LeVeque 2002a).

3 Numerical Methods

In one spatial dimension, the Finite Volume Method (FVM) involves partitioning the domain into discrete intervals known as finite volumes or control volumes and tracking an approximation of the flow integral q over each of these volumes. At each time step, these values are updated by approximating the flow through the boundaries of the intervals.

The mesh contains control volumes of type $\mathcal{C}_{i}=(x_{i-1/2},x_{i+1/2})$, of fixed width $\Delta x=x_{i+1/2}-x_{i-1/2}$, and a time step $\Delta t=t^{n+1}-t^{n}$, as shown in Figure 1. The value $\textbf{Q}^{n}_{i}$ will approximate the average value over the $i$-th interval in time $t_{n}$:

$$\textbf{Q}^{n}_{i}\approx\frac{1}{\Delta x}\int_{x_{i-1/2}}^{x_{i+1/2}}\textbf{q}(x,t_{n})dx\equiv\frac{1}{\Delta x}\int_{\mathcal{C}_{i}}\textbf{q}(x,t_{n})dx,$$ (5)

and

$$\textbf{F}^{n}_{i-1/2}=\mathcal{F}(\textbf{Q}^{n}_{i-1},\textbf{Q}_{i}^{n}),\qquad\textbf{F}^{n}_{i+1/2}=\mathcal{F}(\textbf{Q}^{n}_{i},\textbf{Q}_{i+1}^{n}),$$ (6)

where $\mathcal{F}$ is a numerical flow function on the respective interfaces.

\begin{overpic}[width=227.62204pt,height=85.35826pt,tics=5]{fig1} \end{overpic}

Figure 1: Illustration of the variables in a finite volume method to update the value $\textbf{Q}^{n}_{i}$ by the flux at the interfaces of the control volume.

3.1 High-resolution wave propagation algorithms (WPA)

In this section, we introduce the Wave Propagation Algorithm, an upwind Godunov-type FVM developed by LeVeque (2002b). This algorithm is formulated as a wave propagation scheme for systems of hyperbolic conservation laws and represents a second-order method based on the Lax-Wendroff approach. In fact, it reduces to the classic Lax-Wendroff method in the absence of flux limiters. By incorporating flux limiters (referred to as wave limiters within the context of the Wave Propagation Algorithm), the method achieves high resolution while mitigating non-physical oscillations near discontinuities or steep solution gradients. For a more comprehensive understanding of the Wave Propagation Algorithm, readers are referred to LeVeque (2002b).

This method requires solving a Riemann problem at cell edges to update the cell average 5 in each time step. For an autonomous system, the Riemann problem at $x_{i-1/2}$ is a particular case of a Cauchy problem and consists of the hyperbolic equation with constant coefficients,

$$\textbf{q}_{t}+\textbf{A}\textbf{q}_{x}=0,$$ (7)

together with a special initial condition that is piecewise constant,

$$\textbf{q}(x,0)=\left\{\begin{array}[]{ll}\textbf{Q}_{i-1},&if\;x<x_{i-1/2},\\ \textbf{Q}_{i},&if\;x>x_{i-1/2}.\end{array}\right.$$ (8)

The Riemann solution for a system of $m$ equations typically consists of $m$ wave that we denote by $\mathcal{W}_{i-1/2}^{p}$, for $p=1,2,...,m$, propagating with speeds $s^{p}_{i-1/2}$. The problem can be easily solved in terms of the eigenvectors $\textbf{r}^{p}_{i-1/2}$. The standard approach is to decompose the jump in Q as a linear combination of the eigenvectors to define waves $\mathcal{W}_{i-1/2}^{p}$

$$\textbf{Q}_{i}-\textbf{Q}_{i-1}=\sum_{p=1}^{m}\alpha^{p}_{i-1/2}\textbf{r}^{p}_{i-1/2}\equiv\sum_{p=1}^{m}\mathcal{W}^{p}_{i-1/2}.$$ (9)

The coefficients $\alpha_{i-1/2}^{p}$ are given by

$$\alpha^{p}_{i1/2}=\textbf{R}^{-1}_{i-1/2}(\textbf{Q}_{i}-\textbf{Q}_{i-1}),$$ (10)

where $\textbf{R}_{i-1/2}$ is the matrix or right eigenvectors.

The WPA method takes the form

$$\textbf{Q}_{i}^{n+1}=\textbf{Q}_{i}-\frac{\Delta t}{\Delta x}(\mathcal{A}^{+}\Delta\textbf{Q}_{i-1/2}+\mathcal{A}^{-}\Delta\textbf{Q}_{i+1/2})-\frac{\Delta t}{\Delta x}(\tilde{F}_{i+1/2}-\tilde{F}_{i-1/2}),$$ (11)

where the terms

$$\begin{array}[]{c}\mathcal{A}^{-}\Delta\textbf{Q}_{i-1/2}=\sum\limits_{p=1}^{m}(s^{p})^{-}\mathcal{W}_{i-1/2}^{p},\\ \\ \mathcal{A}^{+}\Delta\textbf{Q}_{i-1/2}=\sum\limits_{p=1}^{m}(s^{p})^{+}\mathcal{W}_{i-1/2}^{p},\end{array}$$ (12)

are called fluctuations and the flux correction term is defined by

$$\tilde{F}_{i-1/2}=\frac{1}{2}\sum_{p=1}^{m}|s^{p}|\left(1-\frac{\Delta t}{\Delta x}|s^{p}|\right)\tilde{\mathcal{W}}^{p}_{i-1/2},$$ (13)

are based on the waves $\mathcal{W}_{i-1/2}^{p}$ and speeds $s_{i-1/2}^{p}$ resulting when solving the Riemann problem for any two states $\textbf{Q}_{i-1}$ and $\textbf{Q}_{i}$ and

$$\tilde{\mathcal{W}}^{p}_{i-1/2}=\psi(\theta^{p}_{i-1/2})\mathcal{W}^{p}_{i-1/2},$$ (14)

is a limited version of the original wave, where

$$\theta^{p}_{i-1/2}=\frac{\mathcal{W}^{p}_{I-1/2}\cdot\mathcal{W}^{p}_{i-1/2}}{\mathcal{W}^{p}_{i-1/2}\cdot\mathcal{W}^{p}_{i-1/2}}\qquad\mathrm{with}\;I=\left\{\begin{array}[]{cc}i-1&if\;s^{p}>0,\\ i+1&if\;s^{p}<0,\end{array}\right.,$$ (15)

and $\psi$ is a wave limiter. In this paper, WPA is implemented using the SuperBee (SB) limiter in the flux limiter version

$$\psi(\theta)=\max(0,\min(1,2\theta),\min(2,\theta)).$$ (16)

As the objective of the study is to verify the behavior of this type of method when applied to the elasticity equations in heterogeneous media, we used the f-wave formulation, proposed in Bale et al. (2003), which is often convenient for equations with spatially varying flux functions. In this case, the fluctuations are defined as

$$\begin{array}[]{ccc}\mathcal{A}^{-}\Delta\textbf{Q}_{i-1/2}&=&\sum_{p:s^{p}_{i-1/2}<0}\mathcal{Z}^{p}_{i-1/2},\\ \mathcal{A}^{+}\Delta\textbf{Q}_{i-1/2}&=&\sum_{p:s^{p}_{i-1/2}>0}\mathcal{Z}^{p}_{i-1/2},\end{array}$$ (17)

where we sum only over the $p$ for which $s^{p}_{i-1/2}$ is negative or positive and replace the correction flow by

$$\tilde{\textbf{F}}_{i-1/2}=\frac{1}{2}\sum_{p=1}^{M_{w}}sgn(s^{p}_{i-1/2})\left(1-\frac{\Delta t}{\Delta x}|s^{p}_{i-1/2}|\right)\tilde{\mathcal{Z}}^{p}_{i-1/2},$$ (18)

where $\tilde{\mathcal{Z}}^{p}$ is a limited version of the wave $\mathcal{Z}^{p}$, that are called f-waves, obtained in the same way that $\tilde{\mathcal{W}}^{p}$ would be obtained from $\mathcal{W}^{p}$.

In the two-dimensional case, the value $\textbf{Q}^{n}_{ij}$ represents an average of cells over the grid cell $(i,j)$ at time $t_{n}$.

$$\textbf{Q}_{i,j}\approx\frac{1}{\Delta x\Delta y}\int_{y_{j-1/2}}^{y_{j+1/2}}\int_{x_{i-1/2}}^{x_{i+1/2}}\textbf{q}(x,y,t_{n})dxdy,$$ (19)

where $\Delta x=x_{i+1/2,j}-x_{i-1/2,j}$, $\Delta y=y_{i,j+1/2}-y_{i,j-1/2}$.

\begin{overpic}[width=170.71652pt,height=142.26378pt]{fig2} \end{overpic}

Figure 2: two-dimensional finite volume grid, where $Q_{i,j}$ represents the average of cells.

In the special case of a rectangular grid of the form $\mathcal{C}_{ij}=[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$, as shown in Figure 2, where $\Delta x=x_{i+1/2}-x_{i-1/2}$ and $\Delta y=y_{j+1/2}-y_{j-1/2}$, we have two different ways to extend the wave propagation algorithm to two dimensions, the first leads us to the Fully Discrete Wave Propagation Algorithm (FWPA) given by

$$\begin{array}[]{ccl}\textbf{Q}_{ij}^{n+1}&=&\textbf{Q}_{ij}-\frac{\Delta t}{\Delta x}(\mathcal{A}^{+}\Delta\textbf{Q}_{i-1/2,j}+\mathcal{A}^{-}\Delta\textbf{Q}_{i+1/2,j})\\ \\ &&-\frac{\Delta t}{\Delta x}(\mathcal{B}^{+}\Delta\textbf{Q}_{i,j-1/2}+\mathcal{B}^{-}\Delta\textbf{Q}_{i,j+1/2})\\ \\ &&-\frac{\Delta t}{\Delta x}(\tilde{\textbf{F}}_{i+1/2,j}-\tilde{\textbf{F}}_{i-1/2,j})-\frac{\Delta t}{\Delta y}(\tilde{\textbf{G}}_{i,j+1/2}-\tilde{\textbf{G}}_{i,j-1/2}),\end{array}$$ (20)

where $\tilde{\textbf{F}}$, and $\tilde{\textbf{G}}$ are calculated via an algorithm that uses information from neighboring cells as described in LeVeque (2002b).

The second way leads us to the Dimensional Splitting Wave Propagation Algorithm (DSWPA), where given the linear system of two-dimensional constant-coefficients

$$\textbf{q}_{t}+\textbf{A}\textbf{q}_{x}+\textbf{B}\textbf{q}_{y}=0,$$ (21)

can be divided into two steps

$$\begin{array}[]{cccc}Step\ 1\ (x-sweeps):&\textbf{q}_{t}+\textbf{A}\textbf{q}_{x}&=&0,\\ Step\ 2\ (y-sweeps):&\textbf{q}_{t}+\textbf{B}\textbf{q}_{y}&=&0.\end{array}$$ (22)

In the $x-sweeps$ we have

$$\textbf{Q}_{i,j}^{*}=\textbf{Q}_{i,j}-\frac{\Delta t}{\Delta x}(\textbf{F}_{i+1/2,j}^{n}-\textbf{F}_{i-2/2,j}^{n}),$$ (23)

where $\textbf{F}^{n}_{i-1/2,j}$ is the WPA numeric flux for the one-dimensional problem between the cells $\mathcal{C}_{i-1,j}$ and $\mathcal{C}_{ij}$ and in the $y-sweeps$ we have

$$\textbf{Q}_{i,j}^{n+1}=\textbf{Q}_{i,j}^{*}-\frac{\Delta t}{\Delta y}(\textbf{G}^{*}_{i,j+1/2}-\textbf{G}^{*}_{i,j-1/2}),$$ (24)

where $\textbf{G}^{n}_{i,j-1/2}$ is the WPA numeric flux for the one-dimensional problem between the cells $\mathcal{C}_{i,j-1}$ and $\mathcal{C}_{ij}$.

3.2 Semi-discrete Central-Upwind Scheme with Kurganov-Lin flux (CUp)

In this section, we describe the semi-discrete Central-Upwind scheme with Kurganov-Lin flux, developed by Kurganov and Lin (2007), which is based on the REA algorithm.

For simplicity, a uniform grid is considered, $x_{i}:=i\Delta x$. It is also assumed that at the current time level, $t=t_{n}$, the average solution cells, $\textbf{Q}^{n}_{i}$, are available. Then, the evolution of the solution computed for the next time level, $t=t_{n+1}$, can be presented by the one-dimensional semi-discrete scheme CUp that has the form of flux difference

$$\frac{d}{dt}\textbf{Q}_{i}(t)=-\frac{\textbf{H}_{i+1/2}(t)-\textbf{H}_{i-1/2}(t)}{\Delta x},$$ (25)

choosing the numerical flux $\textbf{H}_{i-1/2}$ , as

$$\textbf{H}_{i-1/2}(t):=\frac{a^{+}_{i-1/2}\textbf{f}(\textbf{q}^{-}_{i-1/2})-a^{-}_{i-1/2}\textbf{f}(\textbf{q}^{+}_{i-1/2})}{a^{+}_{i-1/2}-a^{-}_{i-1/2}}+a^{+}_{i-1/2}a^{-}_{i-1/2}\left[\frac{\textbf{q}^{+}_{i-1/2}-\textbf{q}^{-}_{i-1/2}}{a^{+}_{i-1/2}-a^{-}_{i-1/2}}-\textbf{q}^{int}_{i-1/2}\right],$$ (26)

we obtain the Central-Upwind Kurganov-Lin scheme. Where the one-sided local velocities, $a^{\pm}_{i-1/2}$, are given by

$$a_{i-1/2}^{+}:=\max_{w\in C(\textbf{q}_{i-1/2}^{-},\textbf{q}_{i-1/2}^{+})}\left\{\lambda_{m}\left(\frac{\partial\textbf{f}}{\partial\textbf{q}}(w)\right),0\right\}\geq 0,$$ (27)

$$a_{i-1/2}^{-}:=\min_{w\in C(\textbf{q}_{i-1/2}^{-},\textbf{q}_{i-1/2}^{+})}\left\{\lambda_{1}\left(\frac{\partial\textbf{f}}{\partial\textbf{q}}(w)\right),0\right\}\leq 0,$$ (28)

where $\lambda_{1}<\cdots\lambda_{m}$ are the eigenvalues of the Jacobian $\frac{\partial\textbf{f}}{\partial\textbf{q}}$, and $C(\textbf{q}_{i-1/2}^{-},\textbf{q}_{i-1/2}^{+})$ is a curve in phase space connecting the corresponding values on the left and right sides of the linear reconstruction at $x=x_{i-1/2}$,

$$\textbf{q}_{i-1/2}^{-}:=\textbf{Q}_{i}-\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i},$$ (29)

$$\textbf{q}_{i-1/2}^{+}:=\textbf{Q}_{i-1}+\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i-1},$$ (30)

respectively, where the numerical derivatives

$$\partial\textbf{q}_{i}^{n}=\left.\frac{\partial\textbf{q}}{\partial x}\right|_{x=x_{i},t=t_{n}}+\mathcal{O}(\Delta x),$$ (31)

are calculated using flux limiters.

The flux limiter used in this work was the SuperBee (SB) defined by

$$\begin{split}\partial\textbf{q}_{i}^{n}=maxmod&\left(minmod\left(2\frac{\textbf{Q}_{i}^{n}-\textbf{Q}_{i-1}^{n}}{\Delta x},\frac{\textbf{Q}_{i+1}^{n}-\textbf{Q}_{i}^{n}}{\Delta x}\right)\right.,\\ &\quad\left.minmod\left(\frac{\textbf{Q}_{i}^{n}-\textbf{Q}_{i-1}^{n}}{\Delta x},2\frac{\textbf{Q}_{i+1}^{n}-\textbf{Q}_{i}^{n}}{\Delta x}\right)\right),\end{split}$$ (32)

where

$$maxmod(x_{1},x_{2},...)=\left\{\begin{array}[]{cc}\max_{i}\{x_{i}\},&if\ x_{i}>0\ \forall i,\\ \min_{i}\{x_{i}\},&if\ x_{i}<0\ \forall i,\\ 0,&other\ case\end{array}\right..$$ (33)

Finally, the correction term in (26) (which is an embedded anti-diffusion term that corresponds to the reduced value in numerical dissipation compared to the original semi-discrete CUp scheme of Kurganov and Tadmor (2000)) is

$$\textbf{q}^{int}_{i-1/2}=minmod\left(\frac{\textbf{q}_{i-1/2}^{+}-\textbf{w}_{i-1/2}^{int}}{a_{i-1/2}^{+}-a_{i-1/2}^{-}},\frac{\textbf{w}_{i-1/2}^{int}-\textbf{q}_{i-1/2}^{-}}{a_{i-1/2}^{+}-a_{i-1/2}^{-}}\right),$$ (34)

with

$$\textbf{w}_{i-1/2}^{int}=\frac{a_{i-1/2}^{+}\textbf{q}_{i-1/2}^{+}-a_{i-1/2}^{-}\textbf{q}_{i-1/2}^{-}-\left\{\textbf{f}(\textbf{q}_{i-1/2}^{+})-\textbf{f}(\textbf{q}_{i-1/2}^{-})\right\}}{a_{i-1/2}^{+}-a_{i-1/2}^{-}}.$$ (35)

For the two-dimensional Central-Upwind Semi-Discrete Scheme with Kurganov-Lin flux, similar to case 1-D, we consider a uniform grid $x_{i}=i\Delta x$, $y_{j}=j\Delta y$, and $\Delta t:=t_{n+1}-t_{n}$. The resulting semi-discrete two-dimensional CUp scheme will then be obtained in the following flow difference form:

$$\frac{d}{dt}\textbf{Q}_{i,j}(t)=-\frac{\textbf{H}^{x}_{i+1/2,j}(t)-\textbf{H}^{x}_{i-1/2,j}(t)}{\Delta x}-\frac{\textbf{H}^{y}_{i,j+1/2}(t)-\textbf{H}^{y}_{i,j-1/2}(t)}{\Delta y},$$ (36)

choosing second-order numeric flux as being

$$\begin{split}\textbf{H}^{x}_{i-1/2,j}(t):=&\frac{a^{+}_{i-1/2,j}\textbf{f}(\textbf{q}^{E}_{i-1,j})-a^{-}_{i-1/2,j}\textbf{f}(\textbf{q}^{W}_{i,j})}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}}\\ &+a^{+}_{i-1/2,j}a^{-}_{i-1/2,j}\left[\frac{\textbf{q}^{W}_{i,j}-\textbf{q}^{E}_{i-1,j}}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}}-\textbf{q}^{x-int}_{i-1/2,j}\right],\end{split}$$ (37)

$$\begin{split}\textbf{H}^{y}_{i,j-1/2}(t):=&\frac{b^{+}_{i,j-1/2}\textbf{f}(\textbf{q}^{N}_{i,j-1})-b^{-}_{i,j-1/2}\textbf{f}(\textbf{q}^{S}_{i,j})}{b^{+}_{,j-1/2}-b^{-}_{i,j-1/2}}\\ &+b^{+}_{i,j-1/2}b^{-}_{i,j-1/2}\left[\frac{\textbf{q}^{S}_{i,j}-\textbf{q}^{N}_{i,j-1}}{b^{+}_{i,j-1/2}-b^{-}_{i,j-1/2}}-\textbf{q}^{y-int}_{i,j-1/2}\right],\end{split}$$ (38)

we get the CUp scheme with Kurganov-Lin flux, unilateral local propagation speeds can be estimated, for example, by

$$\begin{array}[]{ccc}a_{i-1/2,j}^{+}&:=&\max\left\{\lambda_{1}\left(\textbf{A}(\textbf{q}^{W}_{i,j})\right),\lambda_{m}\left(\textbf{A}(\textbf{q}^{E}_{i-1,j})\right),0\right\},\\ b_{i,j-1/2}^{+}&:=&\max\left\{\lambda_{1}\left(\textbf{B}(\textbf{q}^{S}_{i,j})\right),\lambda_{m}\left(\textbf{B}(\textbf{q}^{N}_{i-1,j})\right),0\right\},\\ a_{i-1/2,j}^{-}&:=&\min\left\{\lambda_{1}\left(\textbf{A}(\textbf{q}^{W}_{i,j})\right),\lambda_{m}\left(\textbf{A}(\textbf{q}^{E}_{i-1,j})\right),0\right\},\\ b_{i-1/2,j}^{-}&:=&\min\left\{\lambda_{1}\left(\textbf{B}(\textbf{q}^{S}_{i,j})\right),\lambda_{m}\left(\textbf{B}(\textbf{q}^{N}_{i-1,j})\right),0\right\},\end{array}$$ (39)

where, $\lambda_{1}<\lambda_{2}<\cdots<\lambda_{m}$ are the $m$ eigenvalues of the corresponding Jacobians, $\textbf{A}:=\frac{\partial\textbf{f}}{\partial\textbf{q}}$ and $\textbf{B}:=\frac{\partial\textbf{g}}{\partial\textbf{q}}$, and the point values of linear piecewise reconstruction are given by

$$\begin{array}[]{cc}\textbf{q}^{E}_{i,j}:=\textbf{Q}_{i,j}^{n}+\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}&\qquad\textbf{q}^{N}_{i,j}:=\textbf{Q}_{i,j}^{n}+\frac{\Delta y}{2}(\textbf{q}_{y})^{n}_{i,j},\\ \textbf{q}^{W}_{i,j}:=\textbf{Q}_{i,j}^{n}-\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}&\qquad\textbf{q}^{S}_{i,j}:=\textbf{Q}_{i,j}^{n}-\frac{\Delta y}{2}(\textbf{q}_{y})^{n}_{i,j},\\ \end{array}$$ (40)

where again the numerical derivatives are calculated using flux limiters.

Similarly to the one-dimensional case, we call $\textbf{q}^{x-int}_{i-1/2,j}$ and $\textbf{q}^{y-int}_{i,j-1/2}$ the dissipation correction terms that are given by

$$\begin{split}\textbf{q}^{x-int}_{i-1/2,j}=minmod&\left(\right.\frac{\textbf{q}^{NW}_{i,j}-\textbf{w}^{int}_{i-1/2,j}}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}},\frac{\textbf{w}^{int}_{i-1/2,j}-\textbf{q}^{NE}_{i-1,j}}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}}\\ &\qquad,\left.\frac{\textbf{q}^{SW}_{i,j}-\textbf{w}^{int}_{i-1/2,j}}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}},\frac{w^{int}_{i-1/2,j}-\textbf{q}^{SE}_{i-1,j}}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}}\right),\\ \end{split}$$ (41)

$$\begin{split}\textbf{q}^{y-int}_{i,j-1/2}=minmod&\left(\right.\frac{\textbf{q}^{SW}_{i,j}-\textbf{w}^{int}_{i,j-1/2}}{a^{+}_{i,j-1/2}-a^{-}_{i,j-1/2}},\frac{\textbf{w}^{int}_{i,j-1/2}-\textbf{q}^{NW}_{i,j-1}}{a^{+}_{i-1/2,j}-a^{-}_{i,j-1/2}}\\ &\qquad\left.,\frac{\textbf{q}^{SE}_{i,j}-\textbf{w}^{int}_{i,j-1/2}}{a^{+}_{i,j-1/2}-a^{-}_{i,j-1/2}},\frac{\textbf{w}^{int}_{i,j-1/2}-\textbf{q}^{NE}_{i,j-1}}{a^{+}_{i,j-1/2}-a^{-}_{i,j-1/2}}\right),\end{split}$$ (42)

where

$$\begin{split}\textbf{w}^{int}_{i-1/2,j}&=\frac{a^{+}_{i-1/2,j}\textbf{q}^{W}_{i,j}-a^{-}_{i-1/2,j}\textbf{q}^{E}_{i-1,j}-[\textbf{f}(\textbf{q}^{W}_{i,j})-\textbf{f}(\textbf{q}^{E}_{i-1,j})]}{a^{+}_{i-1/2,j}-a^{-}_{i-1/2,j}},\\ \textbf{w}^{int}_{i,j-1/2}&=\frac{b^{+}_{i,j-1/2}\textbf{q}^{S}_{i,j}-b^{-}_{i,j-1/2}\textbf{q}^{N}_{i,j-1}-[\textbf{f}(\textbf{q}^{S}_{i,j})-\textbf{f}(\textbf{q}^{N}_{i,j-1})]}{b^{+}_{i,j-1/2}-b^{-}_{i,j-1/2}},\end{split}$$ (43)

and $\textbf{q}^{NE}_{i,j}$, $\textbf{q}^{NW}_{i,j}$, $\textbf{q}^{SE}_{i,j}$, $\textbf{q}^{SW}_{i,j}$ are the values of the corresponding corner points of the linear reconstruction by parts in cell $(i,j)$ given by

$$\begin{array}[]{cc}\textbf{q}^{NE}_{i,j}:=\textbf{Q}_{i,j}^{n}+\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}+\frac{\Delta y}{2}(\textbf{q}_{y})^{n}_{i,j},&\qquad\textbf{q}^{SE}_{i,j}:=\textbf{Q}_{i,j}^{n}+\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}-\frac{\Delta y}{2}(\textbf{q}_{y})^{n}_{i,j},\\ \textbf{q}^{NW}_{i,j}:=\textbf{Q}_{i,j}^{n}-\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}+\frac{\Delta y}{2}(\textbf{q}_{x})^{n}_{i,j},&\qquad\textbf{q}^{SW}_{i,j}:=\textbf{Q}_{i,j}^{n}-\frac{\Delta x}{2}(\textbf{q}_{x})^{n}_{i,j}-\frac{\Delta y}{2}(\textbf{q}_{y})^{n}_{i,j}.\\ \end{array}$$ (44)

4 Numerical Results

In this section, we introduce the numerical seismic test cases utilized to assess the performance and competitiveness of the Godunov-type finite volume methods presented earlier, particularly in comparison with the staggered finite difference method. Throughout these examples, we employ reflective or solid wall boundary conditions at the top and bottom and extend the domain on both sides to prevent reflections at these boundaries from interfering with the receivers. This approach allows us to avoid the use of absorbing boundary conditions, which are not the primary focus of our study.

4.1 Test case 1 - Heterogeneous parallel velocity model (Synthetic test case).

In the first test case, our domain spans from $x_{0}=0$ to $x_{1}=1000m$ and $y_{0}=0$ to $y_{1}=1000m$, featuring a velocity profile as depicted in Figure 3. This profile ranges from velocities of 1.5 km/s to 5.5 km/s. We utilize a Ricker wavelet given by:

$$g(t)=(1-2\pi^{2}f_{M}^{2}(t-t_{0})^{2})\exp(\pi^{2}f_{M}^{2}(t-t_{0})^{2}),$$

where $f_{M}=0.015$ kHz represents the peak frequency and $t_{0}=1/f_{M}$. The seismic source is positioned at the domain’s center, at a depth of $20m$, while receivers are placed along the $x$-axis at a depth of $20m$.

In Figure 3, we illustrate the velocity profile along with the positions of the seismic source (highlighted in red) and receivers (depicted in green, spaced every 50m). Additionally, in Figure 4, we present a reference solution of the field $\sigma(x,y)$ at $t=0.3s$ in the top row and the reference seismogram in the bottom row at $t=1s$ using a finite difference method of order 20 implemented with the Devito package (Louboutin et al. 2019, Luporini et al. 2020). This solution was computed on a grid with 6401 points in the $x$ direction and 6401 points in the $y$ direction.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig3.png} \end{overpic}

Figure 3: - Test case 1. Heterogeneous parallel velocity model.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig4_1} \end{overpic}

\begin{overpic}[width=142.26378pt,height=128.0374pt]{fig4_2} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig4_3} \end{overpic}

\begin{overpic}[width=142.26378pt,height=128.0374pt]{fig4_4} \end{overpic}

Figure 4: - Test case 1. Reference solution for the field $\sigma(x,y)$ (Top row) and reference seismogram (Bottom row) computed using a finite difference method of order 20 (DF20) implemented with the Devito package. The numerical grid consisted of 6401 points in both $x$ and $y$ direction, with a spatial resolution of $\Delta x=\Delta y=0.15625m$. The Courant-Friedrichs-Lewy (CFL) number was set to $0.25$

Figure 5 displays the absolute differences between the reference solution and the solutions obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods. These differences are shown for grid spacing of $\Delta x=\Delta y=5m$, at time $t=0.3s$. Figure 6 shows the comparison with the reference solution, calculated using a finite difference method of order 20, the cut was made at $x=330m$. It’s noticeable that the finite difference methods DF2 and DF8 exhibit slight delays compared to the reference solution and the implemented finite volume methods, indicating a slightly higher dispersion. Conversely, the finite volume methods show greater diffusion, as evidenced by Figure 6, where the max-norm errors are higher for finite volume methods, while the 1-norm errors are lower.

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig5_1} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig5_2} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig5_3} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig5_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig5_5} \end{overpic}

Figure 5: - Test case 1. Absolute value of the difference between the reference solution and the solution obtained by implementing the DF2, DF8, DSWPA, FWPA, and CUP methods, using $\Delta x=\Delta y=5m$ at time $t=0.3s$.

Similarly, Figure 7 depicts the absolute differences between the reference seismogram and those obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods for the same grid spacing and time $t=1s$. Figure 8 shows the comparison with the reference solution calculated using a finite difference method of order 20 using the receiver located at $x=330m$. Here, the dispersion introduced by DF2 and DF8 is more apparent, highlighting the higher diffusivity of finite volume methods compared to finite difference methods.

Furthermore, Figure 9 (left) illustrates the computational time required by each method to obtain the seismogram for test case 1. As expected, finite volume methods are more computationally expensive than finite difference methods. Specifically, the computational time required by finite volume methods DFWPA and FWPA on a grid with $\Delta x=\Delta y$ is approximately equivalent to that required by the finite difference method of order 8 on a grid with $\Delta x/2=\Delta y/2$. Figures 9 (center and right) show the error obtained by each method in 1-norm and max-norm versus the computational time required. It’s evident that in 1-norm, the time required by the DF8 method to achieve a certain error is comparable to that required by DSWPA and FWPA for test case 1.

In Figures 10 and 11 we introduced the parameter $\lambda$, representing the percentage decrease in discontinuity in velocity at each layer of the velocity profile presented in this test case, in these figures we show the relationship between $\lambda$ (where $\lambda=0$% corresponds to the original velocity profile with increments of 1km/s at each interface) and the error of each method in both the 1-norm and the max-norm for fixed $\Delta x$ and $\Delta y$. We observe that as $\lambda$ increases, indicating a reduction in the discontinuity at each interface, the error obtained with the finite difference method is significantly diminished, This shows that the greater the jump in the discontinuity in the velocity profile, the greater the error introduced by finite difference method.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig6_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig6_2} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig6_3} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig6_4} \end{overpic}

Figure 6: - Test case 1. At the top row, a comparison with the reference solution is presented for the DF2, DF8, DSWPA, FWPA, and CUP methods using $\Delta x=\Delta y=5m$ (left) and $\Delta x=\Delta y=2.5m$ (right). At the bottom, log-log plots depict the error versus the grid size for the DF2, DF8, DSWPA, FWPA, and CUP methods using the max-norm (left) and 1-norm (right).

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig7_1} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig7_2} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig7_3} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig7_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig7_5} \end{overpic}

Figure 7: - Test case 1. Absolute difference between the reference seismogram and the seismograms obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods, using a spatial resolution of $\Delta x=\Delta y=5m$ at time $t=1s$.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig8_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig8_2} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig8_3} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig8_4} \end{overpic}

Figure 8: - Test case 1. At the top row, a comparison with the reference seismogram is presented for the DF2, DF8, DSWPA, FWPA, and CUP methods using $\Delta x=\Delta y=5m$ (left) and $\Delta x=\Delta y=2.5m$ (right). At the bottom, log-log plots depict the error versus the grid size for the same methods using the max-norm (left) and 1-norm (right).

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig9_1} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig9_2} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig9_3} \end{overpic}

Figure 9: - Test case 1. Plot of the computational time (in seconds) vs. grid size (left), Computational time vs. error in 1-norm for each grid size (center), computational time vs. error in max-norm for each grid size (right), for the methods DF2, DF8, DSWPA, FWPA and CUP, final time $t=1s$.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig10_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig10_2} \end{overpic}

Figure 10: - Test case 1 (variation in discontinuity). Plot of the error vs. the percentage of decrease in the discontinuity in velocity in each layer from the original velocity profile ($\lambda$) for the DF2, DF8, DSWPA, FWPA, and CUP method, with a cutoff at $x=330m$ at $t=0.3s$, the Grid spacing is $\Delta x=\Delta y=10m$ (left) max-norm, (right) 1-norm.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig11_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig11_2} \end{overpic}

Figure 11: - Test case 1 (variation in discontinuity). Plot of the error vs. the percentage of decrease in the discontinuity in velocity in each layer from the original velocity profile ($\lambda$) for the DF2, DF8, DSWPA, FWPA, and CUP method, with a cutoff at $x=330m$ at $t=0.3s$, the Grid spacing is $\Delta x=\Delta y=5m$ (left) max-norm, (right) 1-norm.

4.2 Test case 2 - SEG/EAGE salt body velocity profile.

In this test case, the domain spans from $x_{0}=2000$ to $x_{1}=11000m$ and $y_{0}=0$ to $y_{1}=3000m$, with a velocity profile depicted in Figure 12, with velocities ranging from 1.5km/s to 4.5km/s.

The source was positioned at the middle of the domain at a depth of $40m$, while receivers were situated along the $x$ axis at a depth of $40m$.

Figure 13 displays the absolute difference between the reference seismogram for test case 2 and the seismogram obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods with $\Delta x=\Delta y=20$ at $t=2s$. Figure 14 shows the comparison with the reference solution calculated using a finite difference method of order 20 using the receiver located at $x=7000m$ at the top of the figure, at left and center we show the log-log plots of the error versus the grid size and in the right of the figure, we show a plot of the computational time versus grid size.

\begin{overpic}[width=284.52756pt,height=85.35826pt]{fig12.png} \end{overpic}

Figure 12: - Test case 2. SEG/EAGE salt body velocity profile.

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig13_1} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig13_2} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig13_3} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig13_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig13_5} \end{overpic}

Figure 13: - Test case 2 (SEG/EAGE). Absolute difference between the reference seismogram and the seismograms obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods. These computations were conducted with a spatial resolution of $\Delta x=\Delta y=20m$ and a final time $t=2s$.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig14_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig14_2} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig14_3} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig14_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig14_5} \end{overpic}

Figure 14: - Test case 2 (SEG/EAGE). At the top of the figure, a comparison using the receiver located at $x=7000m$ is presented for the DF2, DF8, DSWPA, FWPA, and CUP methods using $\Delta x=\Delta y=20m$ (left) and $\Delta x=\Delta y=10m$ (right). At the bottom, log-log plots depict the error versus the grid size for the same methods using the max-norm (left) and 1-norm (center). A plot of computational time versus grid size is also provided (right).

4.3 Test case 3 - Marmousi velocity profile.

In this test case, we utilized the Marmousi velocity profile, depicted in Figure 15, originally spanning a horizontal extension of 17km and a depth of 3.5km. Simulations were conducted within a subdomain ranging from $x=4000m$ to $x=14000m$, maintaining the same depth as the original profile. The source was positioned at the middle of the domain at a depth of $40m$, while receivers were distributed along the $x$ axis at the same depth. A frequency of $f_{M}=0.005KHz$ was employed.

Figure 16 exhibits the absolute difference between the reference seismogram for test 3 and the seismogram obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods with $\Delta x=\Delta y=20$ at $t=2.5s$. Figure 14 shows the comparison with the reference solution calculated with a finite difference method of order 20 using the receiver located at $x=7000m$ at the top of the figure, at left and center we show the log-log plots of the error versus the grid size and in the right of the figure, we show a plot of the computational time versus grid size.

\begin{overpic}[width=284.52756pt,height=85.35826pt]{fig15.png} \end{overpic}

Figure 15: - Test case 3. Marmousi velocity profile.

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig16_1} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig16_2} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig16_3} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig16_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig16_5} \end{overpic}

Figure 16: - Test case 3 (Marmousi). Absolute difference between the reference seismogram and the seismogram obtained through the implementation of the DF2, DF8, DSWPA, FWPA, and CUP methods, using $\Delta x=\Delta y=20m$ at $t=2.5s$.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig17_1} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig17_2} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig17_3} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig17_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig17_5} \end{overpic}

Figure 17: - Test case 3 (Marmousi). At the top of the figure, a comparison using the reference receiver located at $x=7000$m is presented for the DF2, DF8, DSWPA, FWPA, and CUP methods using $\Delta x=\Delta y=20m$ (left) and $\Delta x=\Delta y=10m$ (right). At the bottom, log-log plots depict the error versus the grid size for the same methods using the max-norm (left) and 1-norm (center). Additionally, a plot of computational time versus grid size is provided (right).

4.4 Test case 4 - A typical velocity field of Santos Basin.

For the fourth test, we utilized a typical velocity profile from the Santos Basin (as depicted in Figure 18), an expansive oil field covering 352260 square kilometers situated in the Atlantic Ocean, approximately 300 kilometers southeast of Santos, Brazil. Renowned as one of the largest ocean basins in the country, the Santos Basin harbors numerous oil reserve exploration and exploitation sites.

The velocity profile initially spanned a horizontal distance of 70.3 kilometers with a depth of 9.92 kilometers. However, our simulations focused on a subdomain within $x=[20000m,50000m]$, covering a horizontal extension of 30 kilometers, and $y=[0,8000m]$, with a depth of 8 kilometers. The seismic source was positioned at the midpoint of the subdomain, specifically at $x=35000m$ and a depth of $32m$. Receivers were aligned along the $x$ axis at the same depth. We employed a frequency of $f_{M}=0.005KHz$.

Figure 19 showcases a reference seismogram at $t=5s$, calculated using a finite difference method of order 20 using devito in a grid with 9601 points in $x$ and 4001 points in $y$. Additionally, Figure 19 presents the plot of the receiver located at $x=33350m$ at $t=5s$. Figure 20 illustrates the absolute difference between the reference seismogram for test 4 and the seismogram obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods with $\Delta x=25m$ and $\Delta y=8m$ at $t=5s$. Additionally, it presents a comparison with the reference solution using the receiver located at $x=33350m$.

\begin{overpic}[width=284.52756pt,height=85.35826pt]{fig18.png} \end{overpic}

Figure 18: - Test case 4. A typical velocity field of Santos Basin.

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig19_1} \end{overpic}

\begin{overpic}[width=142.26378pt,height=128.0374pt]{fig19_2} \end{overpic}

Figure 19: - Test case 4. Reference seismogram of the $\sigma(x,y)$ field, computed using a DF20 method on a mesh consisting of 9601 points in $x$ and 4001 in $y$ directions. The spatial resolution was $\Delta x=3.125m$ and $\Delta y=2m$, and the simulation was conducted until a final time of $t=5s$, with CFL number set to 0.25

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig20_1} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig20_2} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig20_3} \end{overpic}

\begin{overpic}[width=85.35826pt,height=85.35826pt]{fig20_4} \end{overpic}

\begin{overpic}[width=99.58464pt,height=85.35826pt]{fig20_5} \end{overpic}

\begin{overpic}[width=156.49014pt,height=128.0374pt]{fig20_6} \end{overpic}

Figure 20: - Test case 4. The first five plots display the absolute difference between the reference seismogram and the seismograms obtained using the DF2, DF8, DSWPA, FWPA, and CUP methods, using a spatial resolution of $\Delta x=25m$ and $\Delta y=16m$ at time $t=5s$. The last plot illustrates the comparison with the reference receiver located at $x=33350m$.

We can conclude from the figures 14, 17 and 20 that when implementing finite volume methods in more realistic scenarios (test cases 2, 3, and 4), the finite difference methods DF2 and DF8 outperform the finite volume methods. Moreover, the finite volume methods require more computational time. These figures highlight that the finite volume methods exhibit greater error compared to the finite difference methods, both in the max-norm and the 1-norm, and necessitate more computational time, aligning with our expectations.

5 Conclusions

In our study, we delve into a comparative analysis of two finite-volume methods, which have demonstrated efficacy in fluid dynamics and hyperbolic conservation laws, when applied to seismic wave propagation. This investigation is particularly significant considering the scarcity of studies examining finite volume methods in the context of seismic problems. We juxtapose these methods with the traditional finite difference techniques prevalent in seismic literature.

Our findings reveal nuanced performance characteristics contingent upon the nature of the velocity profiles. When confronted with synthetic velocity profiles exhibiting pronounced discontinuities at interfaces, finite volume methods excel. Despite their higher computational overhead and propensity for greater dissipation, they yield lower numerical errors in such scenarios. However, finite difference methods outshine finite volume methods when velocity profile discontinuities are less pronounced, offering superior computational efficiency and accuracy in realistic seismic problems.

The two methods analyzed in this study are high-resolution techniques and typically achieve up to second-order accuracy. This choice strikes a balance between computational efficiency and desired resolution, making it particularly suitable for addressing challenges posed by heterogeneous media and intricate structures. However, it’s worth noting that these methods can be extended to higher orders if necessary. For instance, the Wave Propagation Algorithm (WPA) can be extended to higher orders as demonstrated in Ketcheson et al. (2013), while similar advancements for the Central-Upwind method have been explored in Kurganov (2016). These extensions offer increased accuracy and precision, however, they also increase their computational cost.

While finite volume methods demonstrate prowess in handling hyperbolic problems with discontinuous solutions—such as those involving shock waves or rarefaction waves—their performance diminishes in linear seismic problems. This is primarily attributable to their higher computational complexity, implementation intricacies, and relatively lower accuracy stemming from significant dissipation introduced by TVD limiters.