Frequency-Domain Finite Element Modeling of Seismic Wave Propagation Under Different Boundary Conditions

Abstract

Seismic wave propagation in complex terrains, especially in the presence of air layers, plays a crucial role in accurate subsurface imaging. However, the influence of different boundary conditions on seismic wave propagation characteristics has not been fully explored. This study employs the finite element method (FEM) to simulate and analyze seismic wavefields under different boundary conditions, including perfectly matched layer (PML), Neumann free boundary conditions, and air layer conditions. First, the finite element solution for the 2D frequency-domain acoustic wave equation is introduced, and the correctness of the algorithm is validated using a homogeneous model. Then, both horizontal and undulating terrain interfaces are designed to investigate the kinematic and dynamic characteristics of the wavefields under different boundary conditions. The results show that PML boundaries effectively absorb seismic waves, prevent reflections, and ensure stable wave propagation, making them an ideal choice for simulating open boundaries. In contrast, Neumann boundaries generate significant reflected waves, particularly in undulating terrains, complicating the wavefield characteristics. Introducing an air layer alters the dynamics of the wavefield, leading to energy leakage and multi-path effects, which are more consistent with real-world seismic-geophysical models. Finally, the computational results using the Overthrust model under different boundary conditions further demonstrate that different boundary conditions significantly affect wavefield morphology. It is essential to select appropriate boundary conditions based on the specific simulation requirements, and boundary conditions with an air layer are most consistent with real seismic geological models. This study provides new insights into the role of boundary conditions in seismic numerical simulations and offers theoretical guidance for improving the accuracy of wavefield simulations in realistic geological scenarios.

1. Introduction

In recent years, seismic wave propagation in complex terrains has attracted considerable attention due to its critical role in understanding wave scattering and diffraction phenomena. However, forward modeling of seismic waves in undulating terrains with air layers and various boundary conditions remains underexplored [1,2]. The presence of an air layer significantly influences seismic wave propagation paths, reflection and refraction behaviors, and energy attenuation characteristics in such terrains. At the same time, the choice of boundary conditions—such as perfectly matched layer (PML), Neumann, or combinations involving air layers—directly affects the accuracy of wavefield simulations. Neglecting the air layer’s impact or improperly applying boundary conditions, particularly in frequency-domain forward modeling, can lead to substantial inaccuracies in wavefield estimation [3,4,5,6,7].

The finite element method was first introduced into seismology by Lysmer and Drake in 1972 [8]. Seron et al. used the finite element method to numerically simulate the propagation process of seismic waves [9,10]. Marfurt pointed out that frequency-domain seismic wave simulation can efficiently simulate wave propagation and attenuation effects in anisotropic media and is very suitable for handling a large number of shot points simultaneously [11]. Sarma et al. proposed non-reflective boundary conditions for elastic wave propagation in finite element numerical simulations [12]. Shi Ruiqi et al. used unstructured grids and finite element methods in the time domain to achieve numerical simulation of seismic wavefields under complex structures [13]. Subsequently, Liu Shaolin et al. investigated the dispersion relations of acoustic and elastic waves in finite element methods using triangular meshes and analyzed the factors influencing these relations [14]. Liu et al. developed a mixed mesh finite element method to simulate the propagation of seismic waves in two-dimensional complex media [1]. Chen et al. proposed a method for generating a three-dimensional finite element mesh of layered geological bodies for finite element numerical simulation of intersecting fault zones, providing a new approach for geological modeling of seismic finite element numerical simulation of fault zones [15]. However, existing seismic simulation methods seldom consider the impact of the air layer, particularly in regions with undulating topography.

When using the wave equation method for seismic wave forward numerical simulation, it is usually necessary to set artificial boundaries to limit the originally infinite simulation area within a finite range. How to construct effective and stable absorbing boundary conditions has also become an important research topic in seismic wave numerical simulation. In seismic numerical simulation, perfectly matched layer or PML boundary conditions are often used. Berenger proposed them in 1994 [16]. By adding an absorbing medium outside the calculation area to absorb the energy of seismic waves transmitted to the boundary, the energy of seismic waves is exponentially attenuated until it disappears. A reasonable setting of the absorbing medium layer can ensure that seismic waves of any frequency and angle will not reflect after entering the absorbing layer. Compared with traditional absorption boundaries, PML absorption boundaries have superior absorption effects and have been widely used in seismic wave simulation [17,18,19]. Neumann free boundaries are commonly used to approximate the ground [20,21,22,23]. The real seismic, geological model has an air layer. Currently, there is limited research on the intrinsic laws of PML boundaries, Neumann boundaries, and air layer boundaries on the wavefield. Their impact on the wavefield should be further clarified.

Therefore, this paper begins with the two-dimensional frequency-domain acoustic wave equation and employs the finite element method to investigate the wavefield characteristics under undulating topography and the presence of an air layer. The study aims to provide theoretical and methodological guidance for seismic numerical simulation and inversion imaging.

2. Methods

2.1. Frequency-Domain 2D Acoustic Wave Equation

The two-dimensional frequency-domain acoustic equation is

(1)${\displaystyle \frac{{𝜕}^{2}p}{𝜕x^2}}+{\displaystyle \frac{{𝜕}^{2}p}{𝜕z^2}}+{\displaystyle \frac{{\omega }^2}{c^2}}p=-S,$

where $p\left(x,z,\omega \right)$ represents the frequency-domain acoustic wavefield, with units of Pa; $c$ represents the acoustic wave velocity, with units of $\mathrm{m}/\mathrm{s}$; $\omega$ represents the angular frequency, with units of rad/s, $\omega =2\pi f$; $f$ represents the frequency; $-S\left(x_0,z_0,\omega \right)$ represents the seismic source; and $\left(x_0,z_0\right)$ represents the coordinates of the seismic source. The seismic source is a Ricker wavelet [24,25], and its time-domain expression $s\left(t\right)$ is given by

(2)$s\left(t\right)=\left[1-2{\left(\pi f_{m}t\right)}^2\right]e^{-{\left(\pi f_{m}t\right)}^2},$

where t represents time, with units of seconds, and $f_m$ denotes the dominant frequency of the Ricker wavelet.

The frequency-domain expression of the Ricker wavelet is

(3)$s\left(\omega \right)={\displaystyle \frac{2f^2}{f_m^3\sqrt{\pi }}}{e}^{-{{\displaystyle \left({\displaystyle \frac{f}{f_m}}\right)}}^2},$

where f represents the frequency of the Ricker wavelet.

For solving Equation (1), the finite element method is employed. Equation (1) can be simplified as

(4)${\nabla }^{2}p+k^{2}p+S=0$

where ${\nabla }^2$ represents the Laplace operator, and $k$ is the wavenumber, $k={\displaystyle \frac{\omega }{v}}$. Taking the variational form of the above equation yields

(5)$F\left(p\right)={\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathsf{\Omega }}{\int }\left[{\left(\nabla p\right)}^2-k^2{p}^2-2Sp\right]}d\mathsf{\Omega }$

Let $\delta F\left(p\right)=0$; the equation can be solved using the finite element method. The detailed process has been extensively described in previous studies and will not be repeated here [11,26,27]. At different frequencies, the finite element equations can all be expressed as the following complex linear system:

(6)$A\left(\omega \right)p\left(\omega \right)=S\left(\omega \right)$

where $A\left(\omega \right)$ is a large complex linear sparse matrix dependent on frequency, while $p\left(\omega \right)$ and $S\left(\omega \right)$ represent the wavefield and source terms, respectively. The finite element method uses triangles as basic elements, with linear interpolation functions. The mesh is adaptively refined, and its schematic diagram is shown in Figure 1.

2.2. Boundary Conditions

In acoustic wave numerical simulation, PML (perfectly matched layer) boundary conditions are commonly used to simulate open boundaries, while Neumann free boundary conditions are applied to simulate the ground surface. These two boundary conditions differ in their principles and application scenarios.

2.2.1. PML Boundary Conditions

The design concept of the PML (perfectly matched layer) absorbing boundary is to introduce a complex coordinate transformation, enabling rapid energy attenuation within the PML region, thereby preventing reflection at the boundaries.

The coordinate transformation is as follows:

(7)$\begin{array}{l}{\displaystyle \frac{𝜕}{𝜕x}}\to {\displaystyle \frac{1}{1+i{\sigma }_x/\omega }}{\displaystyle \frac{𝜕}{𝜕x}}\\ {\displaystyle \frac{𝜕}{𝜕z}}\to {\displaystyle \frac{1}{1+i{\sigma }_z/\omega }}{\displaystyle \frac{𝜕}{𝜕z}}\end{array},$

where ${\sigma }_x$ and ${\sigma }_z$ are the absorption coefficients within the PML region. The absorption coefficient $\sigma$ is a function of the distance L from the inner boundary of the PML to the current point, as shown in the following formula:

(8)$\sigma \left(d\right)={\sigma }_{\max }{\left({\displaystyle \frac{d}{L}}\right)}^m$

where L is the total thickness of the PML layer, set to 20 grid units in this study; m is the polynomial grading order, set to 3 in this study; ${\sigma }_{\max }$ is the maximum absorption coefficient, calculated using the following formula:

(9)${\sigma }_{\max }={\displaystyle \frac{\left(m+1\right)v_{\max }}{2L}}\ln \left({\displaystyle \frac{1}{R}}\right)$

where $v_{\max }$ is the maximum wave velocity in the model, and R is the target theoretical reflection coefficient, set to 10⁻⁶.

This complex coordinate transformation replaces each derivative in the acoustic wave equation with its complex form, enabling energy dissipation. In finite element discretization, the PML boundary condition is implemented by modifying the relevant terms in the matrix A(ω) at the boundary. The specific steps include the following:

1. Mark the PML Region: Add a layer of PML outside the computational domain boundary, with a thickness of 20 grid cells.

2. Adjust the Stiffness Matrix: Within the PML region, modify the relevant entries of the stiffness matrix by incorporating the absorption coefficient to ensure that wave amplitudes gradually decay in the PML region.

3. Ensure Zero Reflection: Choose an appropriate absorption coefficient to ensure that wave energy is completely dissipated at the boundary of the absorbing layer.

The selection of the absorption coefficient is crucial in PML boundary processing. If the absorption coefficient is too large, it can cause the wavefield to attenuate prematurely, leading to wavefield distortion near the boundary. In this case, it is advisable to reduce ${\sigma }_{\max }$ and increase the thickness of the PML layer to smooth the spatial variation of the absorption coefficient and avoid sharp wavefield attenuation. Conversely, if the absorption coefficient is too small, it may result in insufficient absorption, generating residual reflected waves and interfering with the wavefield. To address this, ${\sigma }_{\max }$ and the value of m can be appropriately increased to enhance the absorption effect.

2.2.2. Neumann Free Boundary Conditions

The Neumann boundary condition is a type of natural boundary condition, characterized by specifying the derivative of a variable on the boundary rather than the value of the variable itself. Under this boundary condition, the normal gradient of the wavefield on the boundary is zero,

(10)${\displaystyle \frac{𝜕\tilde{p}}{𝜕n}}=0,$

where $\tilde{p}$ represents the wavefield, and n is the direction normal to the boundary. The Neumann free boundary condition is typically used in wave equations to describe reflective boundaries, making it suitable for simulating scenarios with total reflection.

This paper analyzes the response characteristics of seismic wavefields under undulating topography with an air layer in both the frequency and time domains, focusing particularly on the similarities and differences in the kinematic and dynamic properties of waves under different boundary conditions. The kinematic properties primarily describe the wavefield propagation mode, wave velocity, wavefield morphology, and propagation characteristics, without involving energy transfer or mechanical effects. In contrast, the dynamic properties mainly describe the mechanical effects induced in the medium during wave propagation, including stress, strain, pressure changes, and energy transfer.

3. Results

In this section, a homogeneous model is first designed to verify the correctness of the algorithm. Subsequently, a horizontal interface is introduced, with the seismic source positioned both on the surface and underground. This setup is used to study the similarities and differences in the kinematic and dynamic characteristics of the seismic wavefield in both the frequency and time domains under various boundary conditions. The objective is to investigate the impact of different boundary conditions on the computational results.

3.1. Homogeneous Model

A homogeneous model was designed to verify the accuracy of the frequency-domain finite element algorithm, as shown in Figure 2. The medium velocity is 3000 m/s, with an acoustic point source located at the model’s center. A Ricker wavelet with a dominant frequency of 20 Hz is used as the seismic source. PML absorbing boundary conditions are applied. The frequency-domain wavefields at 20 Hz, 40 Hz, and 60 Hz are shown in Figure 3, demonstrating that the wavefields are essentially unaffected by boundary effects. The accuracy is evaluated using the relative root-mean-square error (Rrms) [28], which is less than 1%, confirming the correctness of the algorithm.

Calculations were performed for frequencies ranging from 1 Hz to 60 Hz to analyze the time-domain wavefield, with 60 frequency components evenly spaced. The time-domain wavefield snapshots were obtained through an inverse FFT transformation, as shown in Figure 4. The seismic record at z = 0 is shown in Figure 5.

From the frequency-domain wavefield of the homogeneous model and the time-domain wavefield obtained through the inverse Fourier transform, it can be concluded that the frequency-domain finite element method for acoustic wave calculation is accurate and reliable. The seismic record further confirms that the time-domain gathers exhibit hyperbolic characteristics, validating the correctness of the method. Therefore, applying this method to study the wavefield distribution characteristics under conditions of undulating topography with an air layer is both practical and feasible.

3.2. Analysis of the Impact of Air Layer on Horizontal Interface

In seismic numerical simulation, different treatments of boundary conditions lead to varying computational results. The seismic source is placed on the surface and within the subsurface medium to study the impact of different boundary conditions on the seismic wavefield. The characteristics of wavefields at different frequencies and time steps under various boundary conditions are analyzed.

3.2.1. The Seismic Source Is Located on the Ground

The model shown in Figure 6 was designed, and the influence of different boundary conditions on the calculation results when the seismic source is located on the ground in seismic exploration was analyzed. In the model, a and b are models without considering the air layer, with a size of 1000 m × 500 m, and c is a model considering the air layer, with a thickness of 100 m. The underground medium has a sound velocity of 3000 m/s, while the air layer has a velocity of 340 m/s. Rayleigh waves with a main frequency of 20 Hz are used. The ground boundary conditions of model a are PML boundary conditions, model b adopts Neumann free boundary conditions, and model c introduces an air layer, which is consistent with the real model; PML absorbing boundary conditions are used above the air layer.

The real and imaginary parts of the frequency-domain wavefield for a horizontal interface with the seismic source positioned on the surface under different boundary conditions are shown in Figure 7 and Figure 8. The results show that under PML boundary conditions, the wavefield exhibits a regular pattern of concentric circles, with energy propagating smoothly outward without significant reflective interference. This is because the PML absorbing boundary effectively absorbs the energy incident on the boundary, making it well-suited for simulating open boundaries and closely approximating a reflection-free scenario. Under Neumann boundary conditions, however, the wavefield amplitude near the surface is significantly larger, displaying reflective effects. This occurs because Neumann boundaries implement free-surface boundary conditions, which do not absorb seismic wave energy, resulting in higher energy concentrations near the surface. As a result, Neumann boundaries are not appropriate for simulating open boundaries. When an air layer is introduced, and PML boundary conditions are applied above the air layer to simulate an open boundary, the wavefield gradually propagates outward within the air layer. The kinematic characteristics of the wavefield remain consistent with those observed under PML and Neumann boundary conditions, but the dynamic characteristics, specifically the amplitude, are altered. Introducing the air layer enables the simulation of coupling effects between the surface and the air. The PML boundary at the top ensures the stability of wave propagation, thereby enhancing the model’s consistency with realistic seismic and geophysical conditions. Consequently, the results better align with actual scenarios.

To clearly illustrate the wavefield morphology near the source, the real and imaginary parts of the 20 Hz wavefield (within x: 400∼600 m, z: 0∼200 m) under different boundary conditions are extracted, as shown in Figure 9. From the red arrows in the figure, it can be observed that the wavefield amplitude near the source under the Neumann boundary condition is higher than that under the PML boundary condition. This is because the Neumann boundary condition causes surface reflections, leading to an increase in wavefield amplitude. The yellow arrows indicate that the presence of an air layer has a significant impact on both the amplitude and morphology of the surface wavefield, suggesting that approximating a free boundary using the Neumann boundary condition introduces certain errors.

To further compare the wavefield patterns under different boundary conditions, the real and imaginary parts at 20 Hz, 40 Hz, and 60 Hz were analyzed, as shown in Figure 10, Figure 11, and Figure 12, respectively. The results indicate that while the kinematic characteristics of the wavefield show similarities under different boundary conditions, there are significant differences in their dynamic characteristics. PML boundary conditions effectively absorb wave energy in the frequency domain and suppress boundary reflections, making them well-suited for simulating open boundaries. The Neumann boundary condition is equivalent to a free boundary, which can cause strong reflection of waves at the boundary. When the incident wave encounters the reflected wave, a standing wave interference phenomenon is formed, especially near the boundary and at specific positions (such as x = 500) where resonance effects occur, resulting in an abnormal increase in local pressure amplitude. In seismic numerical simulations, this boundary condition is often used to approximate the ground surface without introducing an air layer, thereby improving computational efficiency and reducing memory requirements. The introduction of an air layer, however, accounts for the coupling effect between the surface and the air. With the top PML boundary effectively absorbing reflected waves, this setup allows for a more realistic simulation of open boundaries. Consequently, combining the air layer and PML boundary conditions can better replicate real-world surface scenarios, producing frequency-domain wavefield characteristics that align more closely with actual physical conditions.

By comparing the numerical simulation results under different boundary conditions, it is evident that while the kinematic characteristics of seismic waves in the frequency domain remain similar, their dynamic characteristics differ significantly. Therefore, if consistency with real seismic-geophysical models is required in seismic numerical simulations, the combination of an air layer and PML boundary conditions should be adopted. Although this approach demands greater computational resources, it provides results more consistent with the real situation.

After analyzing the frequency-domain wavefield results, the frequency domain results are transformed into the time domain through FFT to obtain time-domain wavefield snapshots under different boundary conditions, as shown in Figure 13 and Figure 14. Figure 13 shows the wavefield snapshot including the air layer, while Figure 14 compares the wavefield snapshots after removing the air layer, which can more intuitively display the wavefield propagation of underground media under different boundary conditions. Combining Figure 13 and Figure 14, it can be seen that under the PML boundary conditions, the simulation of open boundaries has a better effect, with no reflection at the boundary, but there are differences from the real seismic geophysical model. Under Neumann boundary conditions, due to the seismic source being set at z = 0 m, the reflection phenomenon on the ground is not obvious, and there is no obvious reflection around the PML absorption boundary. After adding an air layer and setting a PML absorption boundary above the air layer to simulate an open boundary, the wavefield has a high degree of similarity with the propagation pattern of real seismic waves. The wave’s energy enters the air layer and propagates upwards, forming a multi-path effect. Therefore, under this boundary treatment method, the wavefield characteristics have wave leakage and multi-path effects in kinematics, and some energy is transferred to the air layer in dynamics, which is suitable for simulating the actual situation of surface air coupling.

To further compare the differences in wavefield at different boundaries at the same time, the curves of the wavefield at different depths at 0.05 s, 0.1 s, and 0.15 s were taken for comparison, as shown in Figure 15, Figure 16 and Figure 17. Based on the previous analysis, it is further proven that when the seismic source is located on the ground, its propagation path is basically the same under different boundary treatments, but there are significant differences in wave energy. The PML boundary absorbs the energy of the wavefield at the boundary, and its energy at the surface is lower than that of the Neumann boundary condition. The Neumann boundary does not absorb any energy, and its energy at the surface is the strongest. After adding to the air layer, the wave field energy is minimized due to the entry of wave energy into the air layer. Therefore, the combination of air layer and PML boundary conditions can effectively reduce the influence of reflected waves, but it will introduce certain energy leakage, resulting in a decrease in wave amplitude. This boundary condition is suitable for simulating the coupling effect between the actual surface and air.

When the seismic source is located on the ground, the seismic trace records for z = 0 m under three different boundaries are shown in Figure 18. From the figure, it can be seen that the seismic gathers of PML and Neumann boundaries are basically the same. The reason for this is that the source is located on the ground and there is no obvious reflection, so the shape of their gathers is basically the same. Under the boundary conditions of the air layer and PML, shallow high-frequency interference is introduced in the trace records due to the presence of the air layer. In addition to direct waves, the presence of air waves can also be clearly seen in the trace records. Therefore, waveform processing is required. However, this trace record has higher similarity with the trace records collected from real seismic exploration, which is conducive to simulating the coupling between the surface and the air layer in actual situations.

3.2.2. The Seismic Source Is Located Underground

After in-depth analysis of the wavefield under different boundary conditions when the seismic source is located on the ground, it is known that the boundary has a significant impact on the wavefield morphology. To further investigate its intrinsic mechanism, a model is designed as shown in Figure 19, with the seismic source located in the underground medium at coordinates (250 m, 250 m), and the remaining model parameters remain unchanged. The characteristics of the wavefield under different boundary conditions are studied.

The real and imaginary parts of the frequency-domain wavefield are shown in Figure 20 and Figure 21, respectively. Unlike the case where the seismic source is located on the ground, when the seismic source is located in an underground medium, there are significant differences in the wavefield of the three boundaries. Under the PML boundary conditions, the real and imaginary parts of the wavefield exhibit regular concentric circle diffusion at different frequencies without obvious boundary interference, indicating that the waves propagate uniformly from the source to the surrounding areas, and the PML boundary absorption effect remains good as the frequency increases. Under Neumann boundary conditions, there are reflected waves from the ground, interference, and superposition of wave fields, forming a complex wave field distribution. The higher the frequency, the more complex the interference of the wave field becomes. Therefore, the reflection caused by the Neumann boundary leads to the appearance of reflected waves at the boundary, changing the propagation law of the wave field. The higher the frequency, the more complex the wave field distribution becomes. Under the boundary between the air layer and PML, the underground wavefield is similar to the wavefield under the Neumann boundary, but there is mutual coupling in the wavefield at the interface between the surface and the air. The presence of the air layer causes some of the wavefield energy to leak into the air layer, and the PML layer at the top of the air layer still has a good absorption effect on the wavefield. Therefore, this boundary condition treatment method is consistent with the true seismic geophysical model, which can consider both the existence of the air layer and the simulation of open boundaries.

To further compare the wavefield shapes under different boundary conditions when the seismic source is located underground, the real and imaginary parts at 20 Hz, 40 Hz, and 60 Hz are taken, as shown in Figure 22, Figure 23, and Figure 24, respectively. From the analysis in the figure, it can be seen that on the surface, the Neumann boundary energy is still the strongest, and the wavefield energy is the weakest when there is an air layer present; differently, at depths close to the epicenter, the wavefield energy is much stronger compared to when the epicenter is on the ground. This is because the epicenter is located underground, resulting in less energy leakage into the air. Once again, the three boundary treatment methods have a significant impact on the calculation results of the wavefield, and it is necessary to select the appropriate boundary treatment method according to the specific problem.

The time-domain wavefield snapshots under different boundary conditions when the seismic source is located underground are shown in Figure 25. From the figure, it can be seen that under the PML boundary conditions in the time domain, the wavefield still uniformly diffuses in concentric circles, and the amplitude gradually decays during propagation. At different times, both the wavefield’s kinematic characteristics (propagation law) and dynamic characteristics (energy distribution) conform to the ideal open boundary conditions. They are suitable for simulating ideal open boundaries. Under Neumann boundary conditions, after 0.1 s, there is a clear reflection wave in the wavefield, which overlaps with the direct wave, especially in the vicinity of the boundary, causing the wavefield to lose its regularity. At longer time scales, the propagation of the wavefield becomes more complex. Therefore, in earthquake numerical simulation, to save computational resources, Neumann boundary simulation is often used to simulate the ground, which has a certain rationality for simulating the propagation state of seismic waves. Under the boundary between the air layer and PML, due to the presence of air, the wavefield propagates at a speed of 340 m/s after reaching the air layer, resulting in the attenuation of energy propagation underground. There is also strong reflection in the wavefield at the boundary between the underground medium and the air layer. However, due to a portion of the energy propagating into the air, the energy of the underground reflected waves is weak. Moreover, the vibration direction of the reflected waves in the presence of the air layer is opposite to that under Neumann boundary conditions. Therefore, using Neumann boundary approximation to the ground changes the propagation characteristics of seismic waves, which differs from the actual situation.

When the seismic source is located underground, the comparison curves of the wavefield at different depths under different boundary conditions are shown in Figure 26. From the analysis in the figure, it can be seen that at time 0.05 s, because the wavefield has not yet propagated to the boundary, the wavefield calculation results for the three boundaries are basically consistent; at 0.1 s, the wavefield propagates to the ground and reflects back into the underground medium. It can be inferred that at the ground level, the amplitude of the wavefield is strongest under the Neumann boundary and weakest under the presence of an air layer. This is due to the complete reflection of the Neumann boundary, where there is no energy leakage. Additionally, when the reflected wave meets the incident wave, it forms a standing wave phenomenon, resulting in an abnormal increase in amplitude. The presence of an air layer causes a portion of the wavefield energy to propagate to the air layer. At z = 100 m, there are no reflected waves at the PML boundary, only direct waves, while both the Neumann boundary and the air layer boundary have reflected waves with opposite vibration modes, and the amplitude of the Neumann boundary is stronger than that of the air layer; at the 0.2 s moment of the wavefield, under the PML boundary condition, there are still only direct waves, and the boundary between the Neumann boundary and the air layer still has reflected waves with opposite vibration directions. Once again, it is proven that Neumann boundaries cannot fully simulate real seismic geophysical models, and the air layer needs to be considered for calculation to obtain results that are consistent with the actual situation.

When the seismic source is located underground, the seismic trace records for z = 0 m under three different boundaries are shown in Figure 27. Unlike the case where the earthquake source is located on the ground when the earthquake is located underground, there are no reflected waves in the seismic trace under the PML boundary. However, under the boundary conditions of Neumann and air layer, there are obvious reflected wave records. Due to the short interval between direct waves and reflected waves, they interfere and stack with each other. The trace record is relatively complex compared to PML, and when an air layer exists, airwave records can be seen in the trace records.

3.3. Analysis of the Influence of Undulating Terrain and Air Layer

The previous text has proven that different boundary conditions have a significant impact on the uniformity of the wavefield when the seismic source is on the ground and underground, but they are all studied at the horizontal interface. However, in practical situations, there are undulating terrains, which make the wavefield more complex. To study the influence of three types of boundaries on the wavefield under undulating terrain conditions, a model is designed as shown in Figure 28. The epicenter is located on the ground (500 m, 0 m), with a 100 m high mountain peak on the left side and a 100 m deep mountain depression on the right side. The epicenter is still a Rayleigh wave with a main frequency of 20 Hz. The velocity of the air layer is 340 m/s, and the velocity of the underground medium is 3000 m/s. The wave field propagation characteristics of the model under different boundary conditions are studied separately.

Under undulating terrain conditions, when different boundary treatments are used, the real and imaginary parts of the frequency-domain wavefield are shown in Figure 29 and Figure 30, respectively. From the figure, it can be seen that under the PML boundary, due to its good absorption of the wave field at the boundary, the wave field is basically distributed in concentric circles, and the undulating terrain has little effect on the wave field. Under Neumann boundary conditions, reflected waves are particularly prominent in undulating terrain, where the reflected waves from the uneven terrain overlap with the direct waves, resulting in significant interference in the wavefield distribution. Under the combined conditions of air layer and PML, the undulating terrain and air layer simultaneously affect the wave field morphology, and the wave field energy is weak, which differs from the results of Neumann boundary conditions and is suitable for simulating the surface air coupled wave field.

To clearly illustrate the wavefield morphology near the source and the undulating terrain, the real and imaginary parts of the 20 Hz wavefield (within x: 400∼600 m, z: −150∼150 m) under different boundary conditions are extracted, as shown in Figure 31. Under the PML boundary condition, the wavefield exhibits approximately uniform diffusion. In contrast, the Neumann boundary condition results in strong reflections from the surface with larger amplitudes. When an air layer is present, the wavefields in the air layer and the subsurface medium are coupled, leading to more complex wavefield patterns near the surface.

To further analyze the differences in the wavefield under different boundary conditions under undulating terrain conditions, the real and imaginary parts at 20 Hz, 40 Hz, and 60 Hz were taken, as shown in Figure 32, Figure 33, and Figure 34, respectively. From the analysis in the figure, it can be seen that the amplitude of the wavefield under Neumann boundary conditions is the strongest, while the amplitude of the wavefield under the presence of an air layer is the weakest; Under PML conditions, terrain has less influence, while under Neumann boundary conditions, terrain has a significant impact, making it suitable for studying reflected and scattered waves. The boundary combination of air layer and PML boundary conditions is suitable for surface air coupling simulation, which is consistent with real seismic geological models.

Time-domain wavefield snapshots of different boundary conditions under undulating terrain conditions are shown in Figure 35. Analysis shows the following:

Under PML boundary conditions, the wavefield diffuses relatively uniformly at different times, with slight curvature in undulating terrain areas and no obvious reflected waves observed.

Under Neumann boundary conditions, at 0.05 s, reflected waves have already begun to form, especially in the concave-convex areas of undulating terrain (such as mountaintops and valleys), where the reflected waves overlap with the direct waves, resulting in local energy enhancement of the wavefield. The reflection effect of undulating terrain further amplifies the reflection wave influence of Neumann boundaries, especially in high-energy areas (such as areas with drastic terrain changes), where the amplitude of the wave field increases abnormally. The existence of reflected waves makes the wavefield no longer a single propagation direction wave, but a multi-directional interference wave. Therefore, under undulating terrain conditions, the reflected waves at Neumann boundaries significantly enhance the complexity of the wavefield, and the combined effect of reflected waves and terrain-scattered waves leads to local energy distortion.

When the air layer is combined with the PML boundary, at 0.05 s, the energy of the shallow wavefield is low, manifested as a significant decrease in wavefield amplitude, which is due to energy leakage into the air layer; at 0.1 s and 0.15 s, the shallow wavefield still exhibits leakage effects, with a significant decrease in amplitude. Although shallow wavefield leakage is obvious, the propagation regularity of the deep wavefield (such as areas below terrain) remains good, similar to the propagation characteristics of PML boundary conditions. The terrain scattering effect is not as significant as it is under Neumann boundary conditions in the presence of an air layer, but in the concave context due to shallow energy leakage, the influence of terrain on the wavefield is more pronounced, manifested as a significant curvature of the waveform. The boundary conditions of the air layer are suitable for simulating surface air coupling scenarios, but the amplitude of the wavefield is generally weak.

Therefore, under undulating terrain conditions, different boundary conditions have a more significant impact on the results. For different research objects, accurate boundary conditions need to be selected in order to obtain suitable results.

3.4. Overthrust Land Model Test

To evaluate the effects of different boundary conditions in complex models, the SEG/EAGE Overthrust land model was adopted, as shown in Figure 36. The model dimensions are 5 km × 1.5 km, with the source located at (2.5 km, 140 m). The source is a Ricker wavelet with a dominant frequency of 30 Hz, and adaptive triangular meshing was applied. The model with an air layer was created by adding a 20 m thick air layer above the original model, with PML boundary conditions applied above the air layer. Frequencies from 1 Hz to 80 Hz, with a 1 Hz interval, were transformed to the time domain using fast Fourier transform (FFT) to obtain time-domain wavefield snapshots and seismic gather data. Three boundary conditions were tested on the ground surface: PML, Neumann, and a combination of air layer and PML. The real part of the 30 Hz frequency-domain wavefield is compared in Figure 37, while the time-domain wavefield snapshots are compared in Figure 38. The receivers were placed on the ground with a spacing of 5 m, and the seismic gather comparison is shown in Figure 39.

The frequency-domain wavefield and its comparison curves reveal that different boundary conditions significantly affect the wavefield’s shape near the surface, with the discrepancies in dynamic characteristics increasing as the distance from the source grows. The wavefield amplitudes at the surface for the PML absorbing boundary and the air layer boundary are similar, while the amplitude for the Neumann boundary remains significantly larger. This is due to the strong reflections caused by the Neumann boundary, which lead to standing wave phenomena and result in abnormally large pressure amplitudes. The time-domain wavefield snapshots are consistent with previous test results, showing similar wavefronts at different times. PML boundaries exhibit no ground reflection, while Neumann and air-layer boundaries show reflections with opposite phases. From the seismic gather data, it is evident that PML exhibits minimal boundary reflections, making it the best choice for absorbing boundaries. In contrast, Neumann and air-layer boundaries show ground boundary reflections. Neumann boundaries, due to total energy reflection, produce the highest energy in the seismic gather, while air layers dissipate energy into the air layer, resulting in the lowest energy, which aligns more closely with real-world conditions.

4. Conclusions

This article conducts a systematic comparison and analysis of the seismic wave field propagation characteristics in the frequency domain and time domain using PML absorption boundary conditions, Neumann boundary conditions, and air layer conditions on the ground under horizontal and undulating terrain interfaces through finite element numerical simulation. The summary is as follows:

1. The PML boundary has an excellent absorption effect on the wave field, without reflected wave interference, and can maintain the regularity of wave field propagation. The terrain scattering effect has a relatively small impact on the wave field under PML boundary conditions, and the propagation direction and energy distribution of the wave field are not significantly disrupted, making it suitable for wave field simulation under open boundary conditions.

2. The Neumann boundary has certain advantages in studying complex wave phenomena such as reflected waves and scattered waves, but it is not suitable for simulating wave field propagation under open boundary conditions. Under undulating terrain conditions, Neumann boundaries amplify the interaction between terrain-scattered waves and boundary-reflected waves, significantly enhancing the complexity of the wavefield, manifested as an abnormal increase in shallow local wave amplitude and waveform interference enhancement.

3. When the air layer is combined with PML boundary conditions, there is a significant shallow energy leakage phenomenon in the propagation of the wavefield, especially under high-frequency conditions, where the amplitude of the shallow wavefield decreases significantly. However, the propagation regularity of deep wavefield is good, which is close to the results obtained by using only PML boundary conditions. This indicates that the air layer mainly affects the propagation of shallow wave energy. When considering the air layer, the main impact of undulating terrain is concentrated on the leakage of shallow wavefield and waveform distortion. Therefore, the air layer conditions are suitable for simulating surface air coupling scenarios and have greater similarity with real seismic geophysical models.

4. In existing seismic numerical simulations, Neumann free boundary conditions are commonly used instead of surface conditions, but the results show that this simplification can lead to significant differences in the wave field propagation characteristics compared to real earthquake physical models. To simulate an earthquake wavefield that is closer to the real situation, it is necessary to consider the air layer conditions.

5. The SEG/EAGE Overthrust land model was used to further validate the impact of different boundary conditions on wavefield simulation in complex models.

5. Discussion

This paper mainly discusses the wavefield propagation characteristics under two conditions: horizontal interface and undulating terrain. In addition, for the energy leakage problem under air layer conditions, energy compensation methods or improved boundary conditions can be attempted to improve simulation accuracy. At the same time, the numerical simulation results of this article can be combined with actual earthquake observation data to further verify the applicability of numerical simulation results under different boundary and terrain conditions.

Author Contributions

Conceptualization, S.D., Y.Z. and H.L.; methodology, Y.Z.; validation, Y.Z. and H.L.; formal analysis, Y.Z.; investigation, Y.Z. and H.L.; writing—original draft preparation, Y.Z.; writing—review and editing, Y.Z., H.Z. and H.L.; visualization, Y.Z. and H.L.; supervision, S.D.; project administration S.D. and Y.Z.; funding acquisition, S.D. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Acknowledgments

All contributions have been included by the authors, and there are no other acknowledgments.

Conflicts of Interest

Author Haiyang Liu was employed by the company Changjiang Geotechnical Engineering Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Schematic diagram of model adaptive meshing.

Figure 2. Schematic diagram of the homogeneous model. The red star represents a point source.

Figure 3. Frequency-domain wavefield of the homogeneous model. ((a–c) are the real parts of the wavefield of 20 Hz, 40 Hz, and 60 Hz, and (d–f) are the imaginary parts of the wavefield of 20 Hz, 40 Hz, and 60 Hz).

Figure 4. Time-domain wavefield snapshot of homogeneous model. (a–f) represent the wavefield snapshots at 0.05 s, 0.1 s, 0.15 s, 0.2 s, 0.3 s, and 0.4 s, respectively.

Figure 5. Seismic trace record at z = 0 m.

Figure 6. Models under different boundary conditions when the seismic source is located on the ground. The red star represents a point source. ((a) is the PML absorption boundary condition used on the ground, (b) is the free boundary condition used on the ground, and (c) is the addition of an air layer and the PML absorption boundary used above the air layer).

Figure 7. When the seismic source is located on the ground, the real parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the real parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the real parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the real parts of the wavefield with different frequencies considering the air layer).

Figure 8. When the seismic source is located on the ground, the imaginary parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the imaginary parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the imaginary parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the imaginary parts of the wavefield with different frequencies considering the air layer).

Figure 9. Comparison of wavefield details captured under different boundary conditions at 20 Hz. The red arrow compares the wavefield energy near the source; the yellow arrow compares the wavefield shape near the surface. ((a–c) are the real parts of the 20 Hz wavefield under PML, Neumann, and air layer boundary conditions, respectively; (d–f) are the imaginary parts of the 20 Hz wavefield under PML, Neumann, and air layer boundary conditions, respectively).

Figure 10. When the seismic source is located on the ground, the real and imaginary parts of the 20 Hz wavefield under different boundary conditions. ((a) is the comparison of the real part of the 20 Hz wavefield, and (b) is the comparison of the imaginary part of the 20 Hz wavefield).

Figure 11. When the seismic source is located on the ground, the real and imaginary parts of the 40 Hz wavefield under different boundary conditions. ((a) is the comparison of the real part of the 40 Hz wavefield, and (b) is the comparison of the imaginary part of the 40 Hz wavefield).

Figure 12. When the seismic source is located on the ground, the real and imaginary parts of the 60 Hz wavefield under different boundary conditions. ((a) is the comparison of the real part of the 60 Hz wavefield, and (b) is the comparison of the imaginary part of the 60 Hz wavefield).

Figure 13. Comparison of wavefield snapshots under different boundary conditions when the seismic source is located on the ground. ((a–c) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under PML boundary conditions; (d–f) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under Neumann boundary conditions; and (g–i) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s considering the air layer).

Figure 14. Comparison of wavefield snapshots under different boundary conditions after removing the air layer. ((a–c) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under PML boundary conditions; (d–f) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under Neumann boundary conditions; and (g–i) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s considering the air layer).

Figure 15. Comparison curves of wavefield at different depths at 0.05 s under different boundary conditions when the seismic source is located on the ground.

Figure 16. Comparison curves of wavefield at different depths at 0.1 s under different boundary conditions when the seismic source is located on the ground.

Figure 17. Comparison curves of wavefield at different depths at 0.15 s under different boundary conditions when the seismic source is located on the ground.

Figure 18. Seismic trace records under different boundary conditions when the source is on the ground. (a) PML. (b) Neumann. (c) With air layer.

Figure 19. Model under different boundary conditions when the seismic source is located underground. The red star represents a point source. ((a) is the ground using PML absorption boundary conditions, (b) is the ground using free boundary conditions, and (c) is the addition of an air layer and the use of PML absorption boundaries above the air layer).

Figure 20. When the seismic source is located underground, the real parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the real parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the real parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the real parts of the wavefield with different frequencies considering the air layer).

Figure 21. When the seismic source is located underground, the imaginary parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the imaginary parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the imaginary parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the imaginary parts of the wavefield with different frequencies considering the air layer).

Figure 22. When the seismic source is located underground, the real and imaginary parts of the 20 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 20 Hz wavefield, and (B) is the comparison of the imaginary part of the 20 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 23. When the seismic source is located underground, the real and imaginary parts of the 40 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 40 Hz wavefield, and (B) is the comparison of the imaginary part of the 40 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 24. When the seismic source is located underground, the real and imaginary parts of the 60 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 60 Hz wavefield, and (B) is the comparison of the imaginary part of the 60 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 25. Comparison of wavefield snapshots under different boundary conditions when the seismic source is located underground. ((a–c) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.2 s under PML boundary conditions, respectively; (d–f) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.2 s under Neumann boundary conditions, respectively; and (g–i) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.2 s, respectively, when considering the air layer).

Figure 26. Comparison curves of wavefield at different depths under different boundary conditions when the seismic source is located underground. ((a) is the waveform comparison curve at 0.05 s, (b) is the waveform comparison curve at 0.1 s, and (c) is the waveform comparison curve at 0.2 s).

Figure 27. Seismic trace records under different boundary conditions when the epicenter is located underground.

Figure 28. Model under different boundary conditions in undulating terrain conditions. The red star represents a point source. (The seismic source located on the ground). (a) The ground adopts PML absorption boundary conditions. (b) The ground adopts free boundary conditions. (c) Adds an air layer and uses PML absorption boundary above the air layer.

Figure 29. Under undulating terrain conditions, the real parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the real parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the real parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the real parts of the wavefield with different frequencies considering the air layer).

Figure 30. Under undulating terrain conditions, the imaginary parts of the wavefield with frequencies of 20 Hz, 40 Hz, and 60 Hz under different boundary conditions. ((a–c) are the imaginary parts of the wavefield with different frequencies under PML absorption boundary conditions, (d–f) are the imaginary parts of the wavefield with different frequencies under Neumann boundary conditions, and (g–i) are the imaginary parts of the wavefield with different frequencies considering the air layer).

Figure 31. Comparison of the details of the real and imaginary parts of the wave field captured at 20 Hz under different boundary conditions in undulating terrain conditions. The red arrow compares the wavefield energy near the ground. ((a–c) are the real parts of the 20 Hz wavefield under PML, Neumann, and air layer boundary conditions, respectively; (d–f) are the imaginary parts of the 20 Hz wavefield under PML, Neumann, and air layer boundary conditions, respectively).

Figure 32. Under undulating terrain conditions, the real and imaginary parts of the 20 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 20 Hz wavefield, and (B) is the comparison of the imaginary part of the 20 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 33. Under undulating terrain conditions, the real and imaginary parts of the 40 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 40 Hz wavefield, and (B) is the comparison of the imaginary part of the 40 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 34. Under undulating terrain conditions, the real and imaginary parts of the 60 Hz wavefield under different boundary conditions. ((A) is the comparison of the real part of the 60 Hz wavefield, and (B) is the comparison of the imaginary part of the 60 Hz wavefield. In (A), (a−f) represent the real part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively. In (B), (a−f) represent the imaginary part of the wavefield at depths of 0 m, 100 m, 200 m, 300 m, 400 m, and 500 m, respectively).

Figure 35. Comparison of snapshot images of different boundary wavefields under undulating terrain conditions. ((a–c) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under PML boundary conditions; (d–f) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s under Neumann boundary conditions; and (g–i) are snapshots of the wavefield at 0.05 s, 0.1 s, and 0.15 s when considering the air layer).

Figure 36. Overthrust onshore model.

Figure 37. Real part of 30 Hz seismic wavefield under different boundary conditions. ((a–c) are the real parts of the 30 Hz wavefield under PML, Neumann, and air layer boundary conditions, respectively; (d–f) represent the vertical variation curves of the wavefield at x = 1000 m, 2500 m, and 4000 m, respectively).

Figure 38. Comparison of wavefield snapshots under different boundary conditions. ((a–c) are snapshots of the wavefield at 0.1 s, 0.2 s, and 0.3 s under PML boundary conditions; (d–f) are snapshots of the wavefield at 0.1 s, 0.2 s, and 0.3 s under Neumann boundary conditions; and (g–i) are snapshots of the wavefield at 0.1 s, 0.2 s, and 0.3 s when considering the air layer).

Figure 39. Comparison of seismic traces under different boundary conditions (z = 0 m). (a) PML. (b) Neumann. (c) With air layer.

References

1. Liu, S.; Li, X.; Wang, W.; Liu, Y. A Mixed-Grid Finite Element Method with PML Absorbing Boundary Conditions for Seismic Wave Modelling. J. Geophys. Eng.; 2014; 11, 055009. [DOI: https://dx.doi.org/10.1088/1742-2132/11/5/055009]

2. Abubakar, A.; Habashy, T.M. Three-Dimensional Visco-Acoustic Modeling Using a Renormalized Integral Equation Iterative Solver. J. Comput. Phys.; 2013; 249, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.jcp.2013.04.008]

3. Caunt, E.; Nelson, R.; MacArthur, T.; Luporini, F.; Gorman, G. A Novel Immersed Boundary Approach for Irregular Topography with Acoustic Wave Equations. Geophysics; 2024; 89, pp. T207-T226. [DOI: https://dx.doi.org/10.1190/geo2023-0515.1]

4. Li, B.; Liu, Y.; Sen, M.K.; Ren, Z. Time-Space-Domain Mesh-Free Finite Difference Based on Least Squares for 2D Acoustic-Wave Modeling. Geophysics; 2017; 82, pp. T143-T157. [DOI: https://dx.doi.org/10.1190/geo2016-0464.1]

5. Sochacki, J.; Kubichek, R.; George, J.; Fletcher, W.R.; Smithson, S. Absorbing Boundary Conditions and Surface Waves. Geophysics; 1987; 52, pp. 60-71. [DOI: https://dx.doi.org/10.1190/1.1442241]

6. Li, X.; Yao, G.; Niu, F.; Wu, D. An Immersed Boundary Method with Iterative Symmetric Interpolation for Irregular Surface Topography in Seismic Wavefield Modelling. J. Geophys. Eng.; 2020; 17, gxaa019. [DOI: https://dx.doi.org/10.1093/jge/gxaa019]

7. Li, L.; Wen, X.; Tang, C.; Zhou, D.; Zhang, S. Numerical Simulation of Acoustic Wave Propagation by Finite Element Method Based on Optimized Matrices. J. Geophys. Eng.; 2024; 21, pp. 1027-1039. [DOI: https://dx.doi.org/10.1093/jge/gxae055]

8. Lysmer, J.; Drake, L.A. A Finite Element Method for Seismology. Methods in Computational Physics: Advances in Research and Applications; Elsevier: Amsterdam, The Netherlands, 1972; Volume 11, pp. 181-216. ISBN 978-0-12-460811-5

9. Serón, F.J.; Badal, J.; Sabadell, F.J. A Numerical Laboratory for Simulation and Visualization of Seismic Wavefields. Geophys. Prospect.; 1996; 44, pp. 603-642. [DOI: https://dx.doi.org/10.1111/j.1365-2478.1996.tb00168.x]

10. Serón, F.J.; Sanz, F.J.; Kindelán, M.; Badal, J.I. Finite-Element Method for Elastic Wave Propagation. Commun. Numer. Methods Eng.; 2010; 6, pp. 359-368. [DOI: https://dx.doi.org/10.1002/cnm.1630060505]

11. Marfurt, K.J. Accuracy of Finite-difference and Finite-element Modeling of the Scalar and Elastic Wave Equations. Geophysics; 1984; 49, pp. 533-549. [DOI: https://dx.doi.org/10.1190/1.1441689]

12. Sarma, G.S.; Mallick, K.; Gadhinglajkar, V.R. Nonreflecting Boundary Condition in Finite-Element Formulation for an Elastic Wave Equation. Geophysics; 1998; 63, pp. 1006-1016. [DOI: https://dx.doi.org/10.1190/1.1444378]

13. Shi, R.; Wang, S.; Guo, R.; Zhao, J. Finite-Element Numerical Modeling of Elastic Wave on Unstructured Meshes. Shiyou Diqiu Wuli Kantan/Oil Geophys. Prospect.; 2013; 48, pp. 915-923.

14. Shao-Lin, L.; Xiao-Fan, L.; You-Shan, L.; Tong, Z.; Mei-Gen, Z. Dispersion Analysis of Triangle-Based Finite Element Method for Acoustic and Elastic Wave Simulations. Chin. J. Geophys.; 2014; 57, pp. 712-725. [DOI: https://dx.doi.org/10.1002/cjg2.20135]

15. Chen, Y.; Yang, H.; Ye, Y.; Li, J. Generation of 3D Finite Element Mesh of Layered Geological Bodies in Intersecting Fault Zones. PLoS ONE; 2024; 19, e0293193. [DOI: https://dx.doi.org/10.1371/journal.pone.0293193]

16. Berenger, J.-P. A Perfectly Matched Layer for the Absorption of Electromagnetic Waves. J. Comput. Phys.; 1994; 114, pp. 185-200. [DOI: https://dx.doi.org/10.1006/jcph.1994.1159]

17. Wang, Y.; Xing, W.; Xie, W.; Zhu, Z. Study of Absorbing Boundary Condition by Perfectly Matched Layer. J. China Univ. Pet. Ed. Nat. Sci.; 2007; pp. 19-27. [DOI: https://dx.doi.org/10.3321/j.issn:1000-5870.2007.01.004]

18. Qitong, S.; Youxi, Y. Wavefield Simulation of 2-D Viscoelastic Medium in Perfectly Matched Layer Boundary. Geophys. Prospect. Pet.; 2007; 46, pp. 126-130. [DOI: https://dx.doi.org/10.3969/j.issn.1000-1441.2007.02.005]

19. Li, N.; Xie, L.L.; Zhai, C.H. Comparison of Perfectly Matched Layer and Multi-Transmitting Formula Artificial Boundary Condition Based on Hybrid Finite Element Formulation. Acta Seismol. Sin.; 2007; 20, 12. [DOI: https://dx.doi.org/10.1007/s11589-007-0684-9]

20. Haney, M.M. Generalization of von Neumann Analysis for a Model of Two Discrete Half-Spaces: The Acoustic Case. Geophysics; 2007; 72, pp. SM35-SM46. [DOI: https://dx.doi.org/10.1190/1.2750639]

21. Zampieri, E.; Pavarino, L.F. Implicit Spectral Element Methods and Neumann–Neumann Preconditioners for Acoustic Waves. Comput. Methods Appl. Mech. Eng.; 2006; 195, pp. 2649-2673. [DOI: https://dx.doi.org/10.1016/j.cma.2005.06.005]

22. Gimperlein, H.; Özdemir, C.; Stephan, E.P. Time Domain Boundary Element Methods for the Neumann Problem: Error Estimates and Acoustic Problems. J. Comput. Math.; 2018; 36, pp. 70-89. [DOI: https://dx.doi.org/10.4208/jcm.1610-m2016-0494]

23. Harwood, A.R.G.; Dupère, I.D.J. Calculation of Acoustic Green’s Functions Using BEM and Dirichlet-to-Neumann-Type Boundary Conditions. Appl. Math. Model.; 2015; 39, pp. 4134-4150. [DOI: https://dx.doi.org/10.1016/j.apm.2014.12.027]

24. Wang, Y. Frequencies of the Ricker Wavelet. Geophysics; 2015; 80, pp. A31-A37. [DOI: https://dx.doi.org/10.1190/geo2014-0441.1]

25. Ricker, N. The Form and Laws of Propagation of Seismic Wavelets. Geophysics; 1953; 18, pp. 10-40. [DOI: https://dx.doi.org/10.1190/1.1437843]

26. Boore, D. Finite Difference Methods for Seismic Wave Propagation in Heterogeneous Materials. Methods Comput. Phys. Adv. Res. Appl.; 1972; 11, pp. 1-37.

27. Xu, S. The Finite Element Method in Geophysics; Science Press: Beijing, China, 1994.

28. Zhang, Y.; Dai, S.; Zhang, Q.; Ling, J. Iterative Algorithm for Computing Magnetic Field Considering Remanent Magnetization and Demagnetization. Geophys. J. Int.; 2024; 240, pp. 362-385. [DOI: https://dx.doi.org/10.1093/gji/ggae385]