Numerical Simulation of Elastic Wave in Frequency Domain Based on Generalized Finite Difference Method With a Multiaxial Convolutional Perfectly Matched Layer Boundary

The meshless generalized finite difference method (GFDM) was used to simulate elastic wave propagation in the frequency domain. Dispersion analysis and numerical experiments were conducted to evaluate its performance. The dispersion analysis revealed that the GFDM requires approximately 5.6 grid points in a shear wavelength to suppress the dispersion error of 1% with a regular node distribution. To handle the artificial boundary reflection problem, a multiaxial convolutional perfectly matched layer (MCPML) was introduced to the meshless numerical simulation for its good performance. Numerical experiments showed that the meshless GFDM could effectively avoid stepped diffraction and reduce the total number of nodes for the calculation, by adopting an adaptive node distribution according to the model parameters. Compared to the 25-points finite difference method (FDM), the proposed method demonstrated high calculation accuracy, less memory consumption, and good boundary absorption effect, making it a reliable approach for elastic wave simulation in the frequency domain.


I. INTRODUCTION
In recent years, full-waveform inversion, a high-precision and high-resolution velocity modeling method, has gained widespread attention.This technique can be implemented in both time and frequency domains, with frequency-domain inversion offering advantages such as efficient calculation and flexible data selection compared to time-domain inversion [1].The forward wavefield calculation methods are also different due to the different forms of observational data in each domain.Frequency-domain forward modeling offers several advantages over the time domain, including: 1) only calculating the impedance matrix for each fre- The associate editor coordinating the review of this manuscript and approving it for publication was Mehmet Alper Uslu.quency once when multiple sources are present, leading to significant improvement in calculation efficiency through parallelization; 2) easy simulation of attenuation effects; and 3) no time cumulative error, making it convenient for longterm simulations [2].Over nearly 30 years of development, the basic theory of wavefield simulation in the frequency domain has continuously improved, capitalizing on these advantages.
Lysmer et al. first employed the finite-element method (FEM) for a frequency-domain numerical simulation of seismic waves [3].Marfurt highlighted that simulating seismic waves in the frequency domain is effective for attenuation effects simulation in anisotropic media, and suitable for parallel multi-source simulations [4].Hustedt et al. studied the efficiency of frequency-domain simulations and emphasizing that the main calculational cost of frequency-domain simulations lies in the decomposition algorithm of the impedance matrix.Specifically, LU decomposition accounts for approximately 88% of the calculation time in a frequency slice wavefield solution [5].Stekl and Pratt utilized an optimized 9-point finite difference method to simulate wave propagation in viscoelastic media, requiring four sampling points within a shear wavelength to maintain a 5% dispersion error [6].To enhance accuracy, Min et al. proposed an optimize the 25-points FDM, which necessitates only about four sampling points within a shear wavelength to ensure a 1% dispersion error.The 25-points FDM is widely adopted in frequency domain elastic wave equation forward modeling and full waveform inversion due to its fast calculation efficiency, high precision and ease of implementation [7].However, the FDM employs a fixed step to create grids, leading to stepped diffraction during irregular interface simulations and consequently impacting the numerical results.
Recently, the meshless method has gained popularity due to its ability to avoid mesh generation and discretize the solution area into a series of independent nodes.This approach constructs approximate functions on these discrete nodes to obtain linear equations [8].Because the node positions can change flexibly with the model parameters, the meshless method can effectively avoid grid diffraction [9].The GFDM is a meshless method with simple principles, flexible node arrangement, and high calculation accuracy.It has been widely used to solve various mathematical and engineering problems [10], [11], [12].Jensen et al. first introduced the Taylor series expansion of adjacent nodes around the center point based on a distance function to solve differential equations [13].Benito et al. later developed an explicit formula for GFDM [14].Ureña et al. applied GFDM to forward modeling of elastic wave equations and proposed an adaptive method to minimize the effect of irregular node distribution on dispersion [15].Ureña et al. discussed the dispersion and stability conditions of elastic wave forward modeling with GFDM under regular and irregular nodes [16].Benito et al. further studied GFDM to solve seismic wave propagation in homogeneous isotropic media with perfectly matched layer (PML) boundaries [17].Benito et al. discussed the effect of node settings on simulation accuracy [18], while Benito et al. analyzed the influence of node clouds discretization types (regular and irregular) when applying GFDM to solve elastic wave propagation problems [19].Salete et al. discussed the stability of PML boundaries when applying GFDM to solve a two-dimensional seismic wave propagation problem [20].However, most of these studies were conducted in the time domain, where the stability of the calculation requires careful consideration.Unreasonable node distribution directly leads to instability in the calculation, and even if the parameters in the meshless finite difference meet the stability conditions, long-time wavefield propagation in the time domain may produce unstable high wave number components due to rapid cumulative error growth, leading to an unstable calculation [21].
Takekawa et al. introduced a meshfree finite-difference method (MF-FDM) in the frequency domain to address elastic wave propagation problem.Achieving a dispersion error within 1% requires 6.3, 7.3 and 11.1 grid points per shear wavelength for Poisson's ratios of 0.10, 0.25 and 0.45, respectively [22].Subsequently, Takekawa et al. implemented free boundary conditions directly through nodes on the free surface in the frequency domain using MF-FDM.By comparing simulation results with analytical solutions and the FEM, it was demonstrated that MF-FDM can accurately and efficiently handle free boundary conditions [23].The primary distinction between Takekawa's MF-FDM and Benito's GFDM is that Benito employs a distance-weight function to compute the difference stencils, while MF-FDM does not.Selecting an appropriate weight function theoretically improves calculation accuracy.
The treatment of artificial boundaries is crucial for elastic wave numerical simulation.The perfectly matched layer (PML) boundary condition proposed by Berenger in electromagnetic wave simulation, remains the most widely used method today [24].Weng and Weedon introduced the complex extension coordinate system and formulated the PML boundary [25].In recent years, the PML boundaries have been extensively applied to acoustic and elastic wave simulations.To improve stability in certain media, Kuzuoglu and Mittra proposed complex frequency shifted PML (CPML) boundary conditions [26].Kristel et al. enhanced absorption by proposing a multiaxial perfectly matched layer (MPML), which absorbs boundary reflection in multiple orthogonal directions [27].Tian et al. combined CPML and MPML to create non-splitting MCPML in the time domain, enhancing CPML stability and MPML absorption efficiency [28].In the frequency domain, the boundary layer thickness directly impacts memory usage and computational efficiency.Dong et al. and Liu et al. studied MCPML's performance in the frequency domain, which achieved a good absorption effect in a thinner absorption layer [29], [30].
This study explores GFDM for elastic wave forward modeling in the frequency domain.The methodology section introduces GFDM theory and performs a dispersion analysis.To enhance boundary reflected wave absorption, MCPML is incorporated into the GFDM for elastic wave modeling.Simulation results demonstrate the superiority of the proposed algorithm.The homogeneous model test validates the algorithm's correctness and efficiency by comparing it with the analytical solutions and FDM.The Corner Edge model test reveals that the simulation efficiency can be improved by adopting different node distributions according to the model parameters and reducing the total number of nodes in the calculation.The undulating interface model test highlights GFDM's ability to effectively address stepped diffraction.The Overthrust model test confirms GFDM's suitability for simulating elastic wave propagation in complex media.

II. METHODOLOGY A. ELASTIC WAVE EQUATION
Two-dimensional elastic wave equation in frequency domain can be expressed as follows.
In equation (1), ω is the angular frequency, ρ is the density of the medium, λ and µ are Lame constants, u x and u z represent displacements in the x and z directions, f x and f z represent source terms in the x and z directions, respectively.

B. GFDM FOR SPATIAL PARTIAL DERIVATIVE APPROXIMATION
In GFDM, the Taylor series expansion and weighted least-squares method are employed to approximate derivatives of unknown variables using function values from neighboring nodes within a star.A regular or irregular cloud of points is generated in the computational domain and along the boundary.For a given node i, which is denoted as a central node, the m nearest nodes surrounding the node i will be found.The concept of a star refers to a group of established nodes relative to the central node, as illustrated by the black circle in Fig. 1.Each node has assigned a star.If u 0 is the function value at the central node x 0 of the star and u i (i = 1, 2, • • • , m) is the function value at one of the remaining nodes, the Taylor series expansion around the central node can be expressed in the following form The coefficients h i = x i − x 0 , l i = z i − z 0 , and (x 0 , z 0 ) are the coordinates of the central node, and (x i , z i ) are the coordinates of the node i in the star.By truncating the Taylor series with the fourth derivative, the residual function B (u) is defined as In equation ( 3), w i = w (h i , l i ) represents the star distance weight function of a star.In all the examples considered in this study, the weighting function is chosen as where d i denotes the distance between nodes (x i , z i ) and (x 0 , z 0 ), d max is the distance between the central node and the farthest node in the star.Let Then, the residual function defined in (3) can be expressed in matrix form as where To minimize the residual function B (u) with respect to the unknown derivatives at the central point yields the following linear equation system. And According to ( 10)-( 12), the partial derivative vector D u can be expressed as In equation ( 13), a = m i=1 a i .Expanding the partial derivative vector D u , we can get By substituting equation ( 14) into equation ( 1), a discrete formula for solving the elastic wave equation can be obtained.
In equation (15), u x,0 and u z,0 represent the horizontal and vertical displacements of the central node respectively, whereas u x,j and u z,j represent the horizontal and vertical displacements of the surrounding nodes, respectively.The above equation can be expressed in matrix form.
C is the impedance matrix with size 2N × 2N (N represents the total number of discrete nodes), which is obtained by equation( 15); and F represents the source term.

C. DISPERSION ANALYSIS
To derive the dispersion relation of the GFDM, we consider a plane wave in the frequency domain, expressed as GFDM approximations for the homogeneous elastic wave equations in matrix form are as follows. And where ω is the angular frequency.The wavenumber k = ϕ is the incident angle, h j = h j , l j is a vector of the distance between the node and the center node, and h j = x j − x 0 , l j = z j − z 0 .For a non-trivial solution to (18), the determinant of the coefficient matrix must be zero.This leads to the following equation.
The roots of the above equation are Then the angular frequency can be expressed as According to the reference [31], the distance vector h j = h j , l j can be replaced by h j = (L x , L z ).L x = h j and L z = l j are fixed when node positions are determined.The term 1 2 can then be extracted from A, B and C, and the original angular frequency ω = 2πf = vk.By dividing the calculated angular frequency by the original angular frequency, we can obtain the dispersion function with respect to the phase velocities of the compressional and shear waves.
G = 2π k represents the number of nodes per shear wavelength, and and v s are the phase velocities for the compressional and shear waves.To compare with the conventional FDM, the dispersion curve below is calculated using the consistent rectangular node distribution of the FDM with = dx = dz.
Fig. 2 and Fig. 3 display the dispersion curves of GFDM, MF-FDM and FDM when Poisson's ratio is 0.25 and 0.4 respectively.It can be seen from Fig. 3 and Fig. 4 that dispersion exists in the results obtained using different calculation methods.P wave dispersion is isotropic,  while shear wave dispersion exhibits noticeable anisotropic.It can be observed in the figure that FDM has the smallest dispersion curve fluctuation among the three methods, requiring only 2.5 sampling points per shear wavelength to achieve a 5% error margin.However, for a frequency dispersion limit of 1%, FDM necessitates 8.3 sampling points per shear wavelength when Poisson's ratio is 0.25, while GFDM and MF-FDM require 5.6 and 8.3 points, respectively.At Poisson's ratio of 0.4, FDM requires 12.5 sampling points per shear wavelength, while GFDM and MF-FDM requires 7.1 and 8.3 points, respectively.These findings indicate that GFDM outperforms other methods in media with high Poisson's ratios.

D. BOUNDARY CONDITION
In practice, seismic waves propagate in an infinite underground medium after excitation.Numerical simulations of elastic wave equation are challenging due to the difficulty of simulating an infinite medium.Therefore, it is necessary to limit the model size and create an artificial boundary.If the boundary is not treated numerically, elastic waves will reflect upon reaching it, disrupting the wavefield inside the model.In the frequency domain, boundary layer thickness directly impacts memory usage and computational efficiency.To achieve optimal absorption effect in a thinner absorbing layer, we incorporated MCPML boundary conditions into meshless GFDM for elastic wave modeling.The complex stretching functions for CPML boundary conditions are: κ x and κ z ≥ 1, complex frequency shift factors α x and α z can avoid the singular value of the attenuation factor when the frequency is too small.The incorporation of complex frequency shift factors can enhance the attenuation capability of traditional PML boundaries.The MCPML boundary conditions combine the MPML and CPML boundary conditions by retaining the complex frequency shift coordinate transformation and incorporating mutually orthogonal attenuation factors within the attenuation region.The complex stretching functions are: x and d p z/x is the scale coefficient, which can be obtained through experience or trial calculations (p z/x is set as 1 in this study).
The main direction damping factors were determined by the following equation: is the dominant frequency of the wavelet, l is the distance from a node in the PML layer to the inner boundary of the PML region, L is the thickness of the PML boundary, N is set as 2, q is 1, a is 1.79 in this paper.In the frequency domain, according to the complex stretching functions, we obtain Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
To assess the effectiveness of MCPML boundary conditions, a 1000m×1000m dimension homogeneous medium model was designed with compressional wave velocity of 4000m/s, shear wave velocity of 2309m/s, and density of 2000kg/m 3 .A Ricker wavelet with a dominant frequency of 20Hz served as the source wavelet.Numerical simulations were conducted using 20 layers of PML, CPML, MPML and MCPML boundary conditions for comparison.The X and Z components of the 600ms wavefields are presented in Fig. 4.
It is evident from Fig. 4 that without PML boundaries, the wavefield is obscured by reflections, with no discernible wavefronts.In contrast, with PML, CPML, MPML and MCPML absorption boundaries, clear wavefronts are visible, and the energy of boundary reflections is controlled with different boundary conditions.However, it is still difficult to discern the advantages of MCPML boundary conditions from 600ms wavefields.For further comparison, the shot records in Figure 5 were obtained using 40 PML layers.The proposed MCPML boundary exhibits the best absorption effect, followed by the MPML boundary.
As shown in Figure 6, extracted results from the 80 th channel indicate that the MCPML boundary condition has the highest capability of suppressing boundary reflections.Consequently, the MCPML boundary condition proposed can effectively reduce boundary reflections with fewer absorbing layers and higher computational efficiency.

III. EXAMPLES A. HOMOGENEOUS MODEL
To validate the GFDM's accuracy in solving the elastic wave equation in the frequency domain, a 2D homogeneous medium model was created for forward modeling with dimensions 1000m×1000m.The P-wave velocity was 4000m/s, S-wave velocity was 2309m/s (Poisson's ratio was 0.25) and density was 2000kg/m 3 .A Ricker wavelet with a dominant frequency of 20Hz was used as the source, located at the model center.For comparison, all sources were applied to the vertical displacement component.The GFDM was compared with the analytical solution and 25-points FDM.The node distribution aligned with the FDM's rectangular grid, with a grid (node) spacing of 8m.
The GFDM results with 40 layers of MCPML are presented in Fig. 7 (c) and (d)).For comparison, FDM with 60 layers of MCPML (Fig. 7 (e) and (f)) and 40 layers of MCPML (Fig. 7 (g) and (h)) were also used for forward modeling.The GFDM and FDM results align well with the analytical solution overall, with variations in amplitude near the source location.FDM results with 60 layers of MCPML exhibit smoother events than those with 40 layers, suggesting incomplete absorption of boundary reflections under 40 layers (visible in Fig. 7(g) and (h)).In contrast, GFDM fully absorbs boundary reflections under 40 layers.This is because, in the area close to the outer boundary, FDM uses the reduced order form for difference stencil computation, while GFDM still uses a 22-points difference stencil near the boundary.GFDM has significantly involved more nodes in the attenuation calculation than FDM, and the absorption of boundary reflections can be accomplished in a thicker boundary, thereby reducing the size of the impedance matrix.
The GFDM impedance matrix size was 84872 × 84872, while the FDM's was 121032 × 121032.The number of non-zero elements in the GFDM was 3734368, and that in the FDM was 2745168.However, the difference in calculation time for solving the equation was not significant (Table 1).The solution time for GFDM was 3.0s (single-frequency), and the solution time for FDM was 3.2s (single-frequency).GFDM is simpler in assigning an impedance matrix than FDM, with a total calculation time of 5.6s for GFDM, and 5.9s for FDM (single-frequency).The single-frequency calculation time of the GFDM method was shorter than that of the FDM method.It is worth noting that the numerical simulation in this paper was conducted using C language, Intel icc compiler, CentOs Linux computing platform, and an Intel Xeon Gold 6248 CPU with 2.5GHz frequency, 48 cores and 384GB memory.
For further compare, the GFDM and FDM results were transformed into the time domain to obtain a 300ms wavefields (Fig. 8) and shot records (Fig. 9).The GFDM achieved good boundary reflection absorption using 40 layers of MCPML, reducing calculation nodes.Fig. 10 displays comparison results from the 80 th traces in shot records, revealing close alignment with analytical solutions for both methods, highlighting GFDM's high precision.Additionally, GFDM demonstrated greater calculation efficiency than FDM.
To assess GFDM's performance on media with high Poisson's ratio (S-wave velocity: 1633m/s, Poisson's ratio: 0.4), we extracted the 80 th traces from shot records for comparison (Fig. 11).Before 0.3 seconds, both GFDM and FDM  results closely align with the analytical solution.However, after 0.3 seconds, FDM shows a phase lag, while GFDM remains consistent with the analytical solution, highlighting its precision in such media.

B. CORNER EDGE MODEL
The corner Edge model (Fig. 12) is commonly used to validate frequency-domain forward modeling.The model dimensions were 2000m × 2000m.Regular rectangular nodes (the local node distribution in the green square in Fig. 12(a) is shown in Fig. 12(b)) were initially set with a grid spacing of 10m and a time sampling interval of 0.5ms, maximizing the simulation time to 1.0s.The Ricker wavelet's dominant frequency was 20Hz, and the source was placed at (1000m, 500m).The forward modeling frequency range was 1-60Hz, with a frequency sampling rate of 1Hz, and all ground nodes acted as receivers.
Fig. 13 displays wavefields at 10Hz, 20Hz and 40Hz.The wavefronts radiate outward from the central source, forming circular arcs.Wavefields at a single frequency are time-harmonic and stable.Velocity regions yield distinct wavelengths; phase changes are evident at velocity interfaces with abrupt changes.
All simulated single-frequency wavefields were converted to the time domain using a Fourier transformation.Fig. 14 displays the wavefields computed by GFDM at 300ms,  600ms and 900ms, while Fig. 15 shows the wavefields computed by FDM at the corresponding times.The shot records obtained by GFDM and FDM are displayed in Fig. 16.Comparing the shot records and wavefields at different times, we observed consistent event positions and energy distri-    GFDM achieved a single-shot simulation time of 730.8s, which was more efficient than the 981.4s achieved by FDM.
4696 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
To highlight the benefits of meshless GFDM, a variable density node distribution (Fig. 12(c)) was employed to enhance the efficiency of forward modeling by reducing the total number of nodes.In the low-speed region, the node spacing remained 10m, while in the high-speed region, it was increased to 15m.As a result, the total number of nodes decreased from 103041 to 78541, eliminating 24500 equations.The shot records obtained using both node distributions are depicted in Fig. 17

C. UNDULATING INTERFACE MODEL
The use of FDM in dealing with undulating interfaces can result in false diffraction due to the stepped grid, which in turn affects the accuracy of elastic wave modeling.However, GFDM overcomes this limitation by adapting nodes to changes in the undulating interface, thereby eliminating stepped grid diffraction.The undulating interface model depicted in Fig. 18(a) was employed, with dimensions of 2000m × 2000m.A time-sampling interval of 0.5ms and a total computation time of 2.0s were utilized.The source was positioned at (1000m, 0m), and the Ricker wavelet had a dominant frequency of 20Hz.The forward modeling frequency range spanned 1-60Hz, with a frequency sampling rate of 0.5Hz, and all ground nodes served as receivers.
FDM and GFDM were used for forward modeling.The local grid distribution of the FDM is shown in Fig. 18(b), where the stepped grid can be seen clearly, while the local node distribution of the GFDM is shown in Fig. 18(c), where the nodes are directly distributed on the interface accurately.The seismic records in Fig. 19(c) and (d) show that there are many diffracted waves behind the primary reflected wave by the FDM, whereas the diffraction waves are not seen in the GFDM record (Fig. 19(a) and (b)).The results demonstrate that the forward modeling of the GFDM significantly suppressed the effect of stepped grid diffraction and was suitable for the forward modeling of the formation with undulating interfaces.

D. OVERTHRUST MODEL
The 2D Overthrust model in Fig. 20 was used for elastic waves simulation with dimensions of 10km × 4.64km.The time sampling interval was 0.5ms, with source located at (5000m, 0m), and dominant Ricker wavelet frequency of 10Hz.The forward frequency range was 1-60Hz, with a frequency sampling rate of 0.25Hz.Forward simulations were conducted using FDM and 22 points GFDM simulations were  performed using a rectangular node distribution with a 10m node spacing.We tested that differences in seismic records between 100 and 90 layers of MCPML are negligible, hence we opted for 100-layer boundary results as references (Fig. 21(a) and (f)).Fig. 21(b) and (g) depict the results of 40 layers of MCPML, while Fig. 21(c) and (h) show those of 60 layers of PML, Fig. 21(d) and (i) display the results of 40 layers of PML.Comparing the green boxed records, we observed minimal differences with 40MCPML and 60PML compared to the references.However, clear discrepancies are evident in Fig. 21(d) and (i).These variations are clearer in Fig. 21(e) and (j), which display the 100 th reference traces and differences between various boundaries at the red line (0.8s∼1.5s).The figures suggest that MCPML generally achieves equivalent absorption on boundary reflections in fewer layers than PML.This indicates that MCPML can reduce the impedance matrix size by absorbing boundary reflections in fewer absorption boundary layers than PML.It can be observed that the two methods yield consistent results in terms of the event position and relative amplitude, with minor differences in some details (in the green box), and these differences are relatively small.This result indicates that GFDM forward modeling is suitable for simulating elastic wave propagation in complex media such as FDM.
In addition, Table 3 compares the number of nodes, calculation time, and memory usage of the two methods.It shows that GFDM outperforms FDM in terms of computation time and memory usage.Due to the relatively large scale of the model, the computation occupies a large amount of memory; at the same time, to simulate a wave propagation of 4.0 seconds, it is necessary to calculate many frequencies, so the calculation time is very long.Here, we choose to convert the frequency-domain wavefield to the time domain solely for comparison purposes to validate the method's accuracy.In general applications, such as frequency domain full waveform inversion, we only need to calculate the wavefield in the frequency domain.

IV. CONCLUSION
In this paper, we proposed a GFDM for simulating elastic wave propagation in the frequency domain.We evaluated its effectiveness through dispersion analysis and numerical experiments, introducing the MCPML boundary condition.The following conclusions can be drawn from theoretical analysis and numerical calculations.
(1) GFDM requires fewer sampling points per wavelength to achieve the same accuracy as MF-FDM due to the distance function.It outperforms FDM in media with high Poisson's ratio.
(2) MCPML can absorb the boundary reflection more effectively than traditional PML, CPML and MPML.With the same type of absorption boundary, GFDM achieves good absorption of boundary reflected waves using a thinner absorption boundary than FDM, reducing the number of nodes involved in the calculation and enhancing computational efficiency.
(3) In inhomogeneous medium, adaptive node distribution based on the model parameters improves computational efficiency, making meshless GFDM suitable for complex velocity models with large velocity contrast.
(4) When using GFDM for forward modeling, the nodes can be set in consistent with the actual velocity interface by generating an appropriate distribution of nodes.When the velocity interfaces are accurately depicted, diffractions due to the step grid can be avoided.
(5) For full waveform inversion is performed, only a few frequencies are required for forward modeling in the frequency domain, unlike the time domain where the entire time record is simulated before the next inversion step.Therefore, GFDM is efficient and flexible for frequency domain full waveform inversion.

FIGURE 1 .
FIGURE 1.An irregular cloud of points and the selection of star via distance criterion.

FIGURE 2 .
FIGURE 2. Normalized phase velocities for a Poisson's ratio of 0.25, (a) P waves by GFDM; (b) P waves by MF-FDM; (c) P waves by FDM; (d) S waves by GFDM; (e) S waves by MF-FDM; (f) S waves by FDM.Dispersion curves are plotted for propagation angles of 0 • , 15 • , 30 • , and 45 • with respect to the x-axis.

FIGURE 3 .
FIGURE 3. Normalized phase velocities for a Poisson's ratio of 0.4, (a) P waves by GFDM; (b) P waves by MF-FDM; (c) P waves by FDM; (d) S waves by GFDM; (e) S waves by MF-FDM; (f) S waves by FDM.Dispersion curves are plotted for propagation angles of 0 • , 15 • , 30 • , and 45 • with respect to the x-axis.

z
are the damping factors of the main direction, d (z) x and d (x) z are the damping factors of the secondary direction.α (x) x and α (z) z are the complex frequency shift factors of the main direction, α (z) x and α (x) z are the complex frequency shift factors of the secondary direction.The relationship between the damping factors of the main and secondary directions, and the relationship between the complex frequency shift factors of the main and secondary directions are as follows.

FIGURE 6 .
FIGURE 6.Comparison between the 80 th trace of shot records, (a) X component; (b) Z component.

FIGURE 12 .
FIGURE 12. Corner edge model and local node distribution, (a) Corner edge model; (b) Regular node distribution; (c) Variable density node distribution.

FIGURE 17 .
FIGURE 17. Shot records of Corner edge model obtained by GFDM with different node distribution, (a) Regular node X component; (b) Variable density node X component; (c) Comparison between the 80th trace of X component; (d) Regular node Z component; (e) Variable density node Z component; (f) Comparison between the 80th trace of Z component.

17
(c) and (f)) suggests that a reduction in nodes does not compromise calculation accuracy.Additionally, a comparison of calculation time and memory usage (Table2) reveals a decrease in calculation time by 188.4s (25.8%) and a reduction in memory usage by 191.75Mb (31.0%).

Fig. 21
Fig.21shows the X and Z components of seismic shot records obtained by the GFDM with different absorption boundaries, highlighting the discrepancies between methods.We tested that differences in seismic records between 100 and 90 layers of MCPML are negligible, hence we opted for 100-layer boundary results as references (Fig.21(a) and (f)).Fig.21(b) and (g) depict the results of 40 layers of MCPML, while Fig.21(c) and (h) show those of 60 layers of PML, Fig.21(d) and (i) display the results of 40 layers of PML.Comparing the green boxed records, we observed minimal differences with 40MCPML and 60PML compared to the references.However, clear discrepancies are evident in Fig.21(d) and (i).These variations are clearer in Fig.21(e) and (j), which display the 100 th reference traces and differences between various boundaries at the red line (0.8s∼1.5s).The figures suggest that MCPML generally achieves equivalent absorption on boundary reflections in fewer layers than PML.This indicates that MCPML can reduce the impedance matrix size by absorbing boundary reflections in fewer absorption boundary layers than PML.Fig. 22 displays X and Z components of the seismic shot records obtained by the GFDM and FDM with 40 layers of

Fig. 22
Fig.21shows the X and Z components of seismic shot records obtained by the GFDM with different absorption boundaries, highlighting the discrepancies between methods.We tested that differences in seismic records between 100 and 90 layers of MCPML are negligible, hence we opted for 100-layer boundary results as references (Fig.21(a) and (f)).Fig.21(b) and (g) depict the results of 40 layers of MCPML, while Fig.21(c) and (h) show those of 60 layers of PML, Fig.21(d) and (i) display the results of 40 layers of PML.Comparing the green boxed records, we observed minimal differences with 40MCPML and 60PML compared to the references.However, clear discrepancies are evident in Fig.21(d) and (i).These variations are clearer in Fig.21(e) and (j), which display the 100 th reference traces and differences between various boundaries at the red line (0.8s∼1.5s).The figures suggest that MCPML generally achieves equivalent absorption on boundary reflections in fewer layers than PML.This indicates that MCPML can reduce the impedance matrix size by absorbing boundary reflections in fewer absorption boundary layers than PML.Fig. 22 displays X and Z components of the seismic shot records obtained by the GFDM and FDM with 40 layers of

TABLE 2 .
Comparison table of regular node and variable density node for GFDM simulation.

TABLE 3 .
Comparison table of GFDM and FDM shot records computation.