A 19-Point Average-derivative Optimal Scheme for 3D Frequency-domain Scalar Wave Equation

Forward modeling is an important part of FWI. In line with FWI, forward modeling can be divided into two categories: timedomain modeling and frequency-domain modeling. Compared to time-domain modeling [8], frequency-domain modeling has its advantages: convenient manipulations of a single frequency, multishot computation based on a direct solver, and easy implementation of attenuation [9]. In addition, no wavefield-storage issue occurs when constructing the gradient of FWI for frequency-domain modeling. This is not the case when constructing the gradient of FWI for time-domain modeling. This is because the forward source wavefield and backward receiver wavefield are computed in opposite time direction [10,11].

Forward modeling is an important part of FWI. In line with FWI, forward modeling can be divided into two categories: timedomain modeling and frequency-domain modeling. Compared to time-domain modeling [8], frequency-domain modeling has its advantages: convenient manipulations of a single frequency, multishot computation based on a direct solver, and easy implementation of attenuation [9]. In addition, no wavefield-storage issue occurs when constructing the gradient of FWI for frequency-domain modeling. This is not the case when constructing the gradient of FWI for time-domain modeling. This is because the forward source wavefield and backward receiver wavefield are computed in opposite time direction [10,11].
The main disadvantage of frequency-domain modeling is that it can only be performed in an implicit way by solving a set of linear equations. In comparison with the time-domain modeling, this disadvantage is particularly obvious when it comes to 3D computation. Therefore, reducing the number of grid points per shortest wavelength is in great demand in particular when direct solution techniques are employed. Using a rotated coordinate system, Jo et al. [9] developed a 9-point optimal scheme for 2D scalar wave equation. This 9-point scheme reduces the number of grid points per wavelength to approximately 4, and leads to significant reductions of computer memory and CPU time. Hustedt et al. [12] and Operto et al. [13] generalized the rotatedcoordinate method to variable density case and 3D case, respectively.
To overcome the disadvantage of the rotated optimal 9-point scheme, Chen [14] developed a new 9-point finite-difference scheme for 2D scalar wave equation based on an average-derivative approach [15,16]. This new scheme imposes no restriction of equal directional sampling intervals, and reduces the number of grid points per shortest wavelength to approximately 4 for both equal and unequal directional sampling intervals. Chen [17] further generalized the average-derivative method and developed a 27-point optimal scheme for 3D scalar wave equation. Compared to the classical 7-point scheme, the number of grid points per shortest wavelength is reduced from approximately 13 to approximately 4 by this 27-point optimal scheme for both equal and unequal directional sampling intervals within the relative phase error of one percent.
The 27-point scheme Chen [17] includes 9 optimization coefficients, and has a relatively high complexity. By reducing the number of optimization coefficients from 9 to 5, a simplified 19-point averagederivative optimal scheme for 3D scalar wave equation can be achieved. In the next section, I will present this new 19-point scheme. This is followed by the optimization of coefficients and a numerical dispersion analysis. Numerical examples are then presented to demonstrate the theoretical analysis.

A 19-Point Scheme for 3D Wave Equation
The 3D frequency-domain scalar wave equation can be written as where P is the pressure wavefield, ω is the angular frequency, and v(x, y, z) is the velocity.
An average-derivative optimal 19-point scheme for equation (1) can be obtained as follows ( Figure 1): where From equation (4) where ph V is the phase velocity, , θ is the propagation angle, ϕ is the azimuth angle, and The coefficients α, β, γ, c, and d are obtained by minimizing the phase error: The optimization coefficients for different 1  Tables 2 and 3, respectively. Now I perform numerical dispersion analysis. First, I consider the case where r1=1 and r2=1, which corresponds to the equal directional intervals ∆x=∆y=∆z. Figure 2 shows normalized phase velocity curves of the classical 7-point scheme (3) and the average-derivative optimal 19-point scheme (2) for fixed azimuth angle ϕ and different propagation angles θ. Figure 3 shows normalized phase velocity curves of the classical 7-point scheme (3) and the average-derivative optimal 19-point scheme (2) for fixed propagation angle θ and different azimuth angles ϕ. From these figures, one can conclude that within the phase velocity error of 1%, the classical 7-point scheme (3) requires approximately 13 grid points per shortest wavelength, while the average-derivative optimal 19-point scheme (2) requires approximately 4 points.
In order to obtain an overall estimation of the phase velocity errors varying with ϕ and θ, the following relative phase velocity error is considered:   In equation (2), the approximation of the derivative in one direction involves an average of wavefield values from remaining two directions. In this way, a family of approximations to the derivative is obtained, which depends on a free parameter (α, β, or γ). The free parameters c and d from the average of the mass acceleration term play the same role. Therefore, based on optimization techniques, the optimization approximation can be chosen to reduce dispersion errors. The average-derivative optimal 19-point scheme (2) also includes the classical 7-point scheme as a special case, because when α=0, β=0, γ=0, c=1, and d=0, the scheme (2) reduces to the classical 7-point scheme: . , Ex cos k y cos k z (9) Figure 4 shows the relative phase velocity errors for the classical 7-point scheme and the average-derivative optimal 19-point scheme for different 1. With increasing 1, the relative phase velocity errors for the classical 7-point scheme increase. On the other hand, the relative phase velocity errors for the classical average-derivative optimal 19-point scheme are within 1% as 1 G increases. Figures 5-7 show the results for the case where r1=1 and r2=2. In this case, one can draw the same conclusion: within the phase velocity error of 1% and for equal and unequal directional sampling intervals, the classical 7-point scheme (3) requires approximately 13 grid points per shortest wavelength, while the average-derivative optimal 19-point scheme (2) requires approximately 4 grid points. For other cases on r 1 and r 2 , similar analysis can be made.

Numerical Examples
In this section, two numerical examples are presented to verify the theoretical analysis on the classical 7-point scheme (3) and the average-derivative optimal 19-point scheme (2). First, I consider a homogeneous velocity model with a velocity of 4000 m/s (Figure 8a). Horizontal and vertical distances are x=2 km, y=2 km, and z=1 km, respectively. A Ricker wavelet is placed at the center of the model as a source. The receivers are placed at a depth of 250 m. PML absorbing boundary conditions are applied along the edges of the model. In this numerical example, a monochro-matic wavefield of 20 Hz is computed. According to the criterion of 4 grid points per smallest wavelength, horizontal sampling interval is determined by dx=4000/20/4 m=50 m. Set dy=dx and dz=dx/2. Accordingly, horizontal and vertical samplings are nx=ny=nz=41. In this case, the optimization coefficients for the average-derivative optimal 19-point scheme (2) are α=0.098871, β=0.098869, γ=0.095057, c=0.458533, and d=0.090244. In this example, analytical solution is available to make comparisons with numerical solutions. For the analytical solution, the following formula is used where  is the Fourier transformation of the Ricker wavelet, and       Second, I consider a heterogeneous velocity model. Figure 9a shows a salt dome velocity model. The velocity of the salt dome is 5000 m/s, and the velocity of the overburden is 4000 m/s. The distances, sampling numbers, the Ricker wavelet, and the frequency used in this example are the same as those used in the homogeneous velocity model. For the top boundary, a free-surface condition is used, and for other edges of the model, PML absorbing boundary conditions are applied. The source is located at a depth of 200 m, and the receivers are placed immediately under the free surface. In this example, analytical solution is not available. Therefore, I consider the 27-point optimal scheme in Chen [17] for a comparison. Figure 9b shows the real part of the 20 Hz monochromatic wavefield computed with the 27-point optimal scheme, the classical 7-point scheme (3) and the average-derivative optimal 19-point scheme (2). The simulation result with the averagederivative optimal 19-point scheme (2) is in good agreement with the result of the 27-point scheme. On the other hand, the result with the classical 7-point scheme (3) exhibits large discrepancies with the result of the 27-point scheme due to numerical dispersion.

Conclusions
To simplify the 27-point optimal scheme, a 19-point averagederivative optimal scheme for 3D frequency-domain scalar wave equation is developed. The optimization coefficients are obtained by minimizing the phase velocity errors, and they vary with the ratios