Research Paper Journal of Engineering and Technological Sciences Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave Energy Syamsurizal Rizal1. Awali Priyono2,*. Andri Dian Nugraha2. Mochamad Apri3. Mochamad Agus Moelyadi4 & David P. Sahara2 \\\\\\\\ Graduate Program of the Geophysical Engineering Department. Faculty of Mining and Petroleum Engineering. Institut Teknologi Bandung. Jalan Ganesa No. Bandung 40132. Indonesia. Global Geophysics Research Group. Faculty of Mining and Petroleum Engineering. Institut Teknologi Bandung. Jalan Ganesha No. Bandung 40132. Indonesia Industrial and Financial Mathematics Research Group. Institut Teknologi Bandung. Jalan Ganesa No. Bandung 40132. Indonesia. Department of Aerospace and Aeronautical Engineering. Institut Teknologi Bandung. Jalan Ganesa No. Bandung 40132. Indonesia Corresponding author: awali_p@yahoo. Abstract The numerical dispersion of 2D acoustic wave modeling has become an interesting subject in wave modeling in producing better subsurface images. Numerical dispersion is often caused by error accumulation with increased grid size in wave Wave modeling with high-order finite differences was carried out to reduce the numerical error. This study focused on variations in the numerical order to suppress the dispersion due to numerical errors. The wave equation used in modeling was discretized to higher orders for the spatial term, while the time term was discretized up to the second order, with every layer unabsorbed. The results showed that high-order FD was effective in reducing numerical dispersion. Thus, subsurface layers could be distinguished and observed clearly. However, from the modeling results, the wave energy decreased with increasing distance, so the layer interfaces were unclear. To increase the wave energy, we propose a new source in modeling. Furthermore, to reduce the computational time we propose a proportional grid after numerical dispersion has disappeared. This method can effectively increase the energy of reflected and transmitted waves at a certain depth. The results also showed that the computational time of high-order FD is relatively low, so this method can be used in solving dispersion Keywords: acoustic wave. forward modeling. high-order finite difference. numerical dispersion. proportional grid method. Taylor series. Introduction For the last two decades, the finite difference (FD) method has been an important and relevant topic of study in geophysics . The application of the finite difference method in seismic modeling is done to generate wave propagation in the imaging subsurface . and to determine the location of microseismic sources . using a forward modeling approach. However, wave propagation in seismic modeling often results in numerical dispersion generated by numerical errors because of the grid size and subsurface complexity . Due to numerical dispersion in forward modeling . , the subsurface image will become unclear. The low-order finite difference method can only be applied to smaller grid sizes. In addition, in layers with smaller acoustic impedance, the seismic wave energy is low. To overcome this problem, order variations of the FD method were carried out, ranging from low order to high order and grid size variations based on velocity. Reference . solved the numerical dispersion problem by using a high-order FD and employing the rhombus and Lax-Wendorof schemes with uniform grid sizes. By using these two schemes, it was possible to reduce numerical dispersion to the smallest possible size. Explicit method application to improve spatial and temporal accuracy was carried out by . using a staggered grid method and a uniform grid size. In their study, the Copyright A2023 Published by IRCS - ITB ISSN: 2337-5779 Eng. Technol. Sci. Vol. No. 3, 2023, 402-418 DOI: 10. 5614/j. Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. Taylor expansion approach and a combination of Taylor expansion and least-squares optimization were applied to increase numerical stability. Based on the research, numerical dispersion can be suppressed as much as possible, but it still requires complicated mathematical formulation. Reference . developed a high-order FD with the Lax-Wendroff discretization method. This resulted in high modeling accuracy, where numerical dispersion was reduced optimally but required complicated mathematical formulation. The FD method with an implicit scheme and a uniform grid size, developed by . ,19-. , uses discretization with a high-order cross By using this method, better accuracy was obtained than the lower-order FD. However, the method has not yet been explored in terms of grid size variation based on velocity. The FD method developed in the previous studies used cross-discretization . ,19,. and a staggered-grid scheme . ,24-. Although this method can increase modeling accuracy to the order of 2M, it still shows numerical dispersion for high frequencies. For this reason, an FD method was developed with a hybrid discretization model . , combining the cross and rhombus models, to increase modeling accuracy and reduce numerical dispersion. The results showed that the combined discretization scheme was more accurate and more stable than the cross-discretization scheme. However, this application does require increased computation time because it involves calculation not only at . , . grid points but also at . grid points. Besides that, the wave energy for deeper targets is still low. To reduce numerical dispersion in modeling acoustic waves . , the same discretization was obtained by Reference . for a 3D case. The key idea of the scheme proposed here is to apply a rhombus stencil in 3D space with a constant grid size for the whole domain space and to apply the least square algorithm to generate optimized FD coefficients. Generally. FD modeling is performed using a uniform grid with a small cell size to get clean seismic sections. retain the quality of the seismic cross-section and a relatively short computational time, modeling with a discontinuous grid was carried out by . Reference . generalized the same concept in modeling the 2D P-SV model with staggered-grid FD. Reference . applied a discontinuous grid method to a rough topography. Reference . used various grid FD methods that can solve any integer number for a grid size ratio. Reference . applied a nonuniform grid size to simulate seismic wave propagation using a 3D elastic wave. Reference . applied a discontinuous grid to model 3D finite differences. Reference . executed a discontinuous grid in the staggered-grid FD method using the Lanczos downsampling filter to suspend the high-frequency noise transformed from finer grids to coarser grids. Reference . used a discontinuous grid to solve computational time problems due to near-surface wave propagation using 2D acoustic waves in the frequency domain. In this study, we applied a high-order FD with the cross-discretization scheme. Furthermore, we propose a norm of moment tensor as new a source in wave modeling to increase the wave energy and the proportional grid method to reduce the computational time in propagating the acoustic waves in homogeneous, heterogeneous, and complex mediums. This method was chosen because of its efficiency and effectiveness in increasing wave During the propagating process, we assumed that the waves were unabsorbed in every layer because the waves propagate through several dense layers with very small porosity, so the wave energy is constant in every layer. We also simulated the effect of FD order on the numerical error, so a clean seismogram with a sharp reflector is produced. Applying this method, a clean seismic cross-section with the layer interfaces was obtained. Furthermore, the calculation of computational time was also carried out by evaluating the computational time of low-order to high-order schemes and calculating the cumulative time difference between high-order and loworder. Furthermore, the computational time was evaluated by using code developed in Matlab R2021b. Method The 2D acoustic wave equation is as follows . ,39-. ( ) where P = P. ih, z jh, t nA) is the pressure wavefield, v is the wave velocity in the medium, t is time, x is the coordinate in the horizontal direction, z is the coordinate in the vertical direction/depth, and S. is a scalar function of time and source strength, i is the grid point index in the x-direction, j is the grid point index in the zdirection, and n is the grid point index for the time step, h is the grid size. A is the time step. Eq. is applied to Syamsurizal Rizal et al. a simple source, namely explosions. To increase the wave energy, we propose a new source in 2D acoustic wave equation as follows: | ( )| M( ) = ( ) is the moment tensor representing the source mechanism. | ( )| is a norm of M. Moment tensor is further restructured to form three basic source types: the isotropic (ISO) source, double-couple (DC) sources, and the compensated linear vector dipole (CLVD) source . Coefficient Determination Determining the FD coefficient has been thoroughly described in . Generally, two spatial terms of Eq. above can be described in a high order, while a temporal term can be described in the second order. The terms are then substituted into Eq. Every term is changed to a trigonometric function and further expanded by using the Taylor series. Thus, we obtain the coefficient as follows: 2Oc Oc . Eq. is a special matrix, commonly called the Vandermonde matrix. Proportional grid for layered velocity model with contact impedance. Proportional grid for layered velocity model with uncontact impedance. Figure 1 Illustration of proportional grid method. Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. Proportional grid for complex velocity model. Figure 1 Continued. Illustration of proportional grid method. A Discrete Model of Wave Differential Equation Waves are generated by forward modeling with wave sources at a certain depth in the subsurface. To simulate forward modeling of the wavefield at each grid point at a certain time, wavefield calculations in discrete form need to be made. The 2D acoustic wave differential equation can be written in discrete form as follows: =2 1Oe where co and cm are FD coefficients, v is the wave velocity in the medium, t is the time grid, and AEx is the spatial grid in the x-direction, while the spatial grid in the z-direction is AEz, which is equal to AEx. Previous studies on numerical dispersion usually used a uniform grid size in wave modeling. In this paper, we propose a proportional grid method . to reduce the computational time, where grid size change is set based on the velocity interval, so the grid size will increase with increasing depth. This is done to avoid numerical instability at high velocity. Figure 1 illustrates the proportional grid methods. For velocity v O a, the grid size is OIx1. for a < v O b, the grid size is OIx2. while for v > b, the grid size is OIx3, where the values of a and b are arbitrary velocity values selected from the velocity that exist in the modeling domain. Results To determine the subsurface conditions, forward modeling was carried out, solving the 2D acoustic equation with an isotropic source, so Eq. was changed to: Then, the norm of Eq. was substituted to Eq. ( )| Syamsurizal Rizal et al. is the moment tensor representing the source mechanism. is the moment tensor of an isotropic source. | ( )| is a norm of M. , which is a scalar. In acoustic wave the source is a scalar. Where P = P. ih, z jh, t nA) is the pressure wavefield, v is the wave velocity in the medium, t is time, x is the coordinate in the horizontal direction, z is the coordinate in the vertical direction/depth. In this section, we simulate wave propagation on homogeneous, heterogeneous, and complex mediums, the source of which was located at a certain depth. Homogeneous Velocity Model To generate wave propagation in an isotropic homogeneous medium, we solved the second-order wave differential equation using the FD method, going from the second order to higher orders. In this study, the modeling accuracy only went up to the fourteenth order because modeling with this order has given good results where the seismic records no longer showed any numerical dispersion. If modeling is done at a lower order, numerical dispersion is still visible, while if it is done at a higher order, a longer computation time is required. However, on the other hand, this is adjusted to the numerical dispersion conditions. When modeling at the fourteenth order still appears a numerical dispersion, then it is necessary to do modeling at a higher order. Figure 2. is an illustration of a homogeneous layer where the wave velocity in the medium is 2,700 m/s. Based on the Figure 2. , the wave source is placed at a depth of 1,300 m with a frequency of 25 Hz. To record the signals, receivers . are placed on the surface with a distance of 15. 625 m between one geophone and another, assuming a 2D earth model with a grid size of 15,625 m. The simulation results of wave propagation from the source to the receiver using the model in Figure 2. are shown in Figure 3. ourth orde. and Figure 3. ourteenth orde. Based on the recording in Figure 3. , it is clear that very strong numerical dispersion appeared during the wave propagation simulation. To reduce this, we performed simulations in the sixth (M = . and the tenth order (M = . however, the results of the seismic recordings still showed numerical dispersion. Numerical dispersion was reduced optimally at the fourteenth order (M = . The phenomenon vanished after performing a high-order simulation, as shown by the black arrow in Figure 3. This shows that it was greatly reduced, giving increased modeling accuracy. The phenomenon will cause a low resolution of the seismic section, which ultimately results in uncertainty in the interpretation of seismic data, especially in determining the locations of microseismic sources . Figures 3. show plots of the numerical dispersion in wave modeling in a homogeneous medium. Figure 3. is the numerical dispersion for all receivers and Figure 3. is the numerical dispersion for one receiver. Homogeneous Velocity Model (Vp = 2700 m/. Heterogeneous Velocity Model with Contrast Impedance Consisting of 6-Layers. Complex Velocity Model (Marmousi Velocity Mode. Heterogeneous Velocity Model with Uncontrast Impedance Consisting of 6-Layers Figure 2 Velocity model. Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. Fourth order forward propagation. Shot gather of fourth order. 14th order forward propagation. Shot gather of 14th order. Numerical Error for All Receiver. Numerical Error of Single Seismogram. Figure 3 Forward modeling for homogeneous velocity model with source location at 1,300 m depth. AEx = AEz = AEh = 15,625 m, source frequency 25 Hz. Black arrows show the numerical dispersion. Heterogeneous Velocity Model In this subsection, we apply the numerical method to simulate microseismic propagation in a heterogeneous The subsurface model is shown in Figure 2. This model assumes that every layer of the earth is flat, whose topographical effects can be ignored. The velocity in every layer increases with an increase in depth: Syamsurizal Rizal et al. 2,200 m/s, 2,700 m/s, 3,500 m/s, 4,100 m/s, 4,700 m/s, and 5,300 m/s, respectively. The geophone distance was made equal to the size of the computing grid, which was 15,625 m. The source was placed in a subsurface depth of 1,300 m with a frequency of 25 Hz. The simulation results using the model in Figure 2. are shown in Figure Based on Figure 4. , it can be seen that there was a large numerical dispersion, as shown by the black arrow, which is confirmed by the recording display in Figure 4. , wherein the numerical error obtained in the experiment was 0. 0219 to 3. To reduce the numerical error, simulations were then performed using the sixth, eighth, and tenth order. However, the results showed that at M = 2, 3, and 5, the numerical error was still The numerical error highly decreased in the experiment with the fourteenth order FD (M = . , so the numerical dispersion has disappeared in Figure 4. The numerical error in this simulation was about 2 x 10-9 to 6 x 10-3. Figure 4. shows the shot gather in a heterogeneous medium for the FD of the fourteenth order. Figure 4. ourteenth orde. shows better results when compared to Figure 4b . ourth orde. When compared with the wavefield in Figures 4. , the wavefield in Figures 4. and 4d increased by about 93. Figures 4e and 4. show plots of the numerical dispersion at wave modeling in a homogeneous medium. Figure 4. show the numerical dispersion for all receivers and Figure 4. shows the numerical dispersion for some To observe the effect of increasing modeling accuracy on numerical dispersion in a complex heterogeneous medium, a simulation of wave propagation was carried out on the Marmousi velocity model . , as shown in Figure 2. The results of the forward modeling on the complex heterogeneous medium are shown in Figures 5a and 5. The simulation at the fourth order shows a very large numerical dispersion, as can be seen from the shot gather in Figure 5. For this reason, an increase in modeling accuracy was carried out with M = 7 . ourteenth orde. Based on this simulation, a clean shot gather (Figure 5. was obtained, but the numerical dispersion had not completely disappeared. Figures 5. show plots of the numerical dispersion at wave modeling in a homogeneous medium. Figure 5. shows the numerical dispersion for all receivers and Figure 5. shows the numerical dispersion for some receivers. Discussion Impact of Modeling Parameter As explained in the introduction, numerical dispersion is affected by the grid size and subsurface complexity. this section, we explain the effect of increasing the grid size on numerical dispersion. The magnitude of the numerical dispersion is proportional to the grid size. Modeling with a very large grid size will produce numerical dispersion because the stability factor is very small compared to the stability requirements, namely O 1 Oo2 (A is the stability facto. Therefore, eliminating numerical dispersion still requires a large M value. Meanwhile, modeling with a very small grid size will cause numerical instability because the stability factor is greater than the standardized stability criterion with > 1 . The next simulation was performed with reducing If the Oo2 number of layers but the grid size remains 15. 625 m, the numerical dispersion phenomenon still appears as in the 6 layers model so that eliminating numerical dispersion still requires a fairly large M value. However, if the grid size gets smaller, the phenomenon of numerical dispersion is still visible, and to eliminate numerical dispersion does not need high modeling accuracy. Some different results acquired with Changes in the FD accuracy in complex mediums as in low-order FD, the subsurface layer cannot be distinguished due to overlapping wavefields affected by numerical dispersion but the high-order FD reduces wavefield overlapping caused by numerical dispersion so that the subsurface layers can be seen clearly. Numerical dispersion makes wave energy dispersed into several dispersion waves so that the amplitude of the wave appears to widen. After the finite difference order is increased, the accuracy increases because the wave energy is focused. Based on Figures 4 and 5, it can be seen that seismic wave energy depends on the distance from the wave source . due to layer heterogeneity and complexity . When the distance increases, the energy of seismic wave energy decreases. In order to increase the energy at a great distance, we propose a new source, as mentioned in Eq. The small grid size is very effective for reducing the numerical dispersion phenomenon because the grid spacing is constrained by the lowest velocity in the model . However, this requirement will lead to an Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. increase in computational time. To reduce the computational cost, a nonuniform grid size has been widely applied using fine and coarse grids to discretize low and high-velocity areas, respectively . ,50-. Fourth order forward propagation. Shot gather of fourth order. 14th order forward propagation. Shot gather of 14th order. Numerical Error for All Receiver. Numerical Error of Single Seismogram. Figure 4 Forward modeling in heterogeneous velocity model . , source location at 1,300 m depth. AEx = AEz = AEh = 15,625 m. The source frequency is 25 Hz. Black arrows show the numerical dispersion. In addition to reducing the computational time, the use of a discontinuous grid can reduce multiples that arise due to the subsurface complexity. Some factors that complicate the near-surface conditions are strong heterogeneity, topographic relief, and strong attenuation. This problem can be overcome by using a finer grid Syamsurizal Rizal et al. size on the near-surface and a coarser grid for deeper areas . For higher velocity, wave modeling with a very small grid size will cause numerical dispersion. Conversely, for a low velocity, modeling with the largest grid size will cause numerical instability. Therefore, a balance between velocity and grid size is needed that satisfies the stability criteria. To reduce the computational time, we propose a proportional grid method. In this method, the grid size increases with increasing velocity. 4th order forward propagation. Shot gather of 4th order. 14th order forward propagation. Shot gather of 14th order with 25 Hz frequency. Numerical Error for All Receiver. Numerical Error of Single Seismogram. Figure 5 Forward modeling in Marmousi velocity model source location at 350 m depth. AEx = AEz = AEh = 15,625 m. The source Frequency is 25 Hz. The black ellipse shows the numerical dispersion. Figure 6 shows the seismic recordings before and after the new source and the proportional grid were applied. The velocity model used is the heterogeneous velocity model shown in Figure 2b, namely v = 2,200 m/s. 2,700 Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. m/s. 3,500 m/s. 4,100 m/s. 4,700 m/s. 5,300 m/s from top to bottom. In Figure 6a, a seismic recording can be seen with some unclear layer boundaries due to amplitude attenuation, even though the layers have the largest acoustic impedance. During propagation, the waves get an energy loss so that the wave amplitude decreases because of repeated reflections by layer boundaries to a certain depth. Furthermore, for the heterogeneous velocity model . ix layer. , the proportional grid method (Figure 1. ) was applied, with the following grid < 3000 AE Ie OI = OIEa = 5,625 O 4500 AE Ie OI = OIEa = 10,625 > 4500 AE Ie OI = OIEa = 15,625 Modeling results using this method show that the energy of transmitted and reflected waves increased even though the depth increased because the waves had an increase in amplitude, as shown in Figure 6. Therefore, increasing the amplitude will increase the wave energy. From Figure 6. it can be seen that there was a slight increase in numerical dispersion, but this can be solved by increasing the order of finite differences (M valu. Figure 6. was modeled at M = 8, and the result was that the numerical dispersion was reduced slightly. Figure 6. , the numerical dispersion has vanished when modeling was done at M = 10. Before Applied New Source and Proportional Grid Method, 14th Order. After Applied New Source and Proportional Grid Method 14th Order. After Applied New Source and Proportional Grid Method 16th Order. After Applied New Source and Proportional Grid Method 20th Order. Figure 6 Seismic recording of forward modeling on heterogeneous models . ix layer. with a source frequency of 25 Hz. P wave velocity from top to bottom is 2,200 m/s. 2,700 m/s. 3,500 m/s/. 4,100 m/s. 4,700m/s. 5,300 m/s. Source Location at 1,300 m depth from the surface. Black arrows show the numerical dispersion, while red arrows show the amplitude gain. Furthermore, subsurface modeling was carried out on layers with low acoustic impedance. The velocity model used a heterogeneous model consisting of six layers (Figure 2. The parameters of the P wave velocity from Syamsurizal Rizal et al. the top to the bottom layer are shown in Figure 2d, namely 2,200 m/s, 2,400 m/s, 2,600 m/s, 2,800 m/s, 3,000 m/s, and 3,200 m/s, respectively. The source location was at a depth of 650 m. Figure 7. shows the seismic record of forward modeling with a conventional source and a uniform grid size (OIh = 15,625 . From Figure 7. it can be seen that the energy of the reflected wave was smaller than the energy of the reflected wave, as shown in Figure 6. due to the smaller acoustic impedance. To increase the energy of reflected wave, modeling was carried out by applying a new source, with the following grid conditions (Figure 1. < 2400 AE Ie OI = OIEa = 5,625 O 2800 AE Ie OI = OIEa = 10,625 > 2800 AE Ie OI = OIEa = 15,625 The boundaries between layers are clearly visible in Figure 7. To reduce the numerical dispersion that appears in Figure 7. , the modeling accuracy was improved by increasing the FD order up to M = 6 and M = 8. Figure 7. is the seismic recording at the twelfth order using a new source and a proportional grid. From Figure 7. , it can be seen that the numerical dispersion has disappeared. Modeling using this method was continued at the sixteenth order, and a clearer seismic section was obtained, as shown in Figure 7. The new source plays a role in reducing wave energy dissipation from the source to all directions. Variation in grid size play a role in reducing numerical dispersion and computational time. Before Applied New Source and Proportional Grid Method, 8th Order. After Applied New Source and Proportional Grid Method, 8th Order. After Applied New Source and Proportional Grid Method, 12th Order. After Applied New Source and Proportional Grid Method, 16th Order. Figure 7 Seismic recording of forward modeling on the heterogeneous model . , with a source frequency of 20 Hz. P wave velocity from top to bottom was 2,200 m/s. 2,400 m/s. 2,600 m/s/. 2,800 m/s. 3,000m/s. 3,200 m/s. Source location at 650 m depth from the surface. Black arrows show the numerical dispersion, while red arrows show the amplitude gain. Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. Figure 8 is a comparison of the seismic records before and after applying a new source and the proportional grid size in a complex velocity model (Marmousi mode. Figure 8. is the seismic record before applying the new source and the proportional grid method, while Figure 8. is the simulation after applying both methods. using a new source and a variety of grid sizes based on velocity, the boundaries between layers become more visible, but it creates a slight numerical dispersion effect, as shown in Figure 8. Thus, to eliminate the numerical dispersion, we carried out modeling at order 12 (M = . , as shown in Figure 8. To get the best results, the order of finite difference needed to be increased again, namely, to order 16 (M = . The modeling results are shown in Figure 8. The grid conditions (Figure 1. ) for the Marmousi model were as follows: < 2200 AE Ie OI = OIEa = 2,625 O 2600 AE Ie OI = OIEa = 5,625 O 3200 AE Ie OI = OIEa = 7,625 O 4000 AE Ie OI = OIEa = 10,625 O 5000 AE Ie OI = OIEa = 12,625 > 5000 AE Ie OI = OIEa = 15,625 . Before Applied New Source and the Proportional Grid Method, 8th Order. After Applied New Source and the Proportional Grid Method 8th Order. After Applied New Source and the Proportional Grid Method 16th Order. After Applied New Source and the Proportional Grid Method 12th Order. Figure 8 Seismic recording of forward modeling on the Marmousi model with a source frequency of 20 Hz, source location at 650 m depth from the surface. The black ellipse shows the numerical dispersion, while the red ellipse shows the amplitude gain. Syamsurizal Rizal et al. Calculation of Computational Time The effectiveness of FD computation was reviewed by calculating the computational time for the complex velocity model. In this experiment, we simulated low order (M = . to high order (M = . Based on the experimental results, it was found that the computational time for the fourth order was 142. 14757 seconds . 36913 minute. , for the sixth order it was 186. 31678 seconds . 10528 minute. , and so on. Table 1 shows the computation time using our method in a complex medium. When compared to what was done by Jing et al. , the percentage of CPU time from the code developed was smaller. Based on Table 2, the computational time was 1 s smaller than the computational time in Table 1 because the algorithm did not run using Matlab. Table 3 shows the computational time of some high-order FD. From Table 3 it can be seen that the computation time was relatively high. Table 1 Table 2 Computational time of low-order to high-order finite difference for the complex velocity model. Order Grid Step . Time Step . 15,625 15,625 15,625 15,625 15,625 Computational Time Second Minute 142,14757 2,36913 186,31678 3,10528 248,30288 4,13838 317,26189 5,28770 433,17542 7,21959 51,4561 46,9904 41,8448 39,4061 37,1729 The computational time for generating no visible numerical dispersion in a 2D complex model . Method Lax-Wendroff correction Ae LWC 8 Etgen Optimized time-space-domain finite difference Method Ae OpsTS Table 3 Grid Step Meter Time Step . CPU Time Second Computational time of acoustic wave modeling to reduce visible numerical dispersion . Method The Conventional SFD Scheme M = 8, it = 1. 0 ms Temporal Fourth-Order SFD Scheme M = 8, it = 1. 0 ms . Temporal Sixth-Order SFD Scheme M = 8, it = 1. 0 ms . High-Order Temporal and Spatial TE Based SFD Scheme M = 8. N = 4, it = 1. 0 ms . Optimal SFD Scheme M = 8, it = 1. 0 ms . High-Order Temporal and Spatial TE LS Based SFD Scheme M = 8. N = 2, it = 1. 0 ms . Computational time . Conclusions We simulated seismic wave propagation in isotropic homogeneous and isotropic heterogeneous mediums using the high-order FD method. Simulations were carried out on six-layer in a heterogeneous medium. Based on the results of this study, the high-order FD showed better shot gather than the low-order FD. This was due to the higher-order FD which provided higher accuracy so that the numerical error was reduced. The application of a new source was done to increase the energy of reflected and transmitted waves at low acoustic impedance velocity models. This method can increase the energy of reflection and transmission waves so that the boundaries between layers will be clearly visible, especially layers with large depths. In a velocity model with many layers, the energy of the reflected waves will decrease with increasing depth so that the wave amplitude is reduced. To overcome this problem, a method is needed to maintain the wave energy so that a sharp reflector is obtained. The proportional grid method is very effective in reducing the computational time. The computational time between the low-order method and the high-order method was relatively low. This was reviewed both cumulatively and non-cumulatively. Therefore, the high-order difference can be used in solving dispersion problems. Reducing Numerical Dispersion with High-Order Finite Difference to Increase Seismic Wave DOI: 10. 5614/j. Acknowledgments The authors gratefully acknowledge the Indonesian Directorate General of Higher Education (DIKTI) for the research funding. We also thank to Jamhir Safani and La Ode Ngkoimani for their fruitful discussions. References