Seismic forward modeling of acoustic surface-related order-separated multiples

Seismic surface-related multiples have become a hot topic of great significance due to the buried geological information provided by broader illumination areas than primaries. In recent years, researchers attempt to extract the hidden hint of multiples rather than treating them as noise and eliminating them directly. The elimination methods, e.g., the surface-related multiple elimination (SRME) and the inverse scattering series free-surface multiple elimination (ISS-FSME), may be affected by the overlapping or proximity of primaries and multiples. Typical imaging methods, e.g., the reverse time migration (RTM) and the least-square reverse time migration (LSRTM), suffer severe crosstalk artifacts from multiples of inappropriate order and smooth migration velocities. To study the characteristics of primaries and surface-related multiples, whether for elimination or imaging, we propose a forward modeling method of acoustic surface-related order-separated multiples established on the areal/virtual source assumption. The free surface is replaced with an absorbing surface under the dipole source approximation and the ghost creation approach. We present two reflection operators to approximate the reflection at the free surface and apply them to the areal source to obtain ideal results. Numerical experiments on three models prove the effectiveness of the proposed forward modeling method of acoustic surface-related order-separated multiples.


Introduction
Traditionally, surface-related multiples are regarded as noise and should be eliminated in conventional primary-based data processing to prevent misleading the final geological interpretation. From another point of view, multiples are provided with broader illumination ranges than primaries, driving researchers to extract useful information from them (Lu 2021). Taking multiples as effective signals, one can map them into superior images with fewer shots.
Various multiple elimination approaches have been developed in the past few decades that can be classified into three classes in general. The first class is filtering-based methods, such as the predictive deconvolution method relying on the periodicity difference between primaries and multiples (Peacock and Treitel 1969;Taner 1980) and the transformation methods founded on the separability of primaries and multiples in a specific domain (Foster and Mosher 1992;Schonewille and Aaron 2007;Yilmaz 2001). Unlike the first class of techniques, the second class is wave-equation-based methods provided with higher precision and a wider range of applications. The surface-related multiple elimination, abbreviated to SRME (Berkhout and Verschuur 1997;Verschuur et al. 1992), has developed into a widely accepted method in the industry. SRME predicts surface-related multiples from the acquired seismic data without further underground information if the source estimation and construction of traces near the source are processed beforehand (Verschuur 2013). Weglein et al. (1997) presented an alternative multiple elimination method that is also data-driven, namely, the inverse 3 Page 2 of 17 scattering series free-surface multiple elimination (ISS-FSME). This method is applicable for both surface-related and internal multiples if the corresponding subseries of the inverse scattering series are used. The SRME predicts an accurate travel time and approximate amplitude and phase of multiples, which may cause problems when using adaptive subtraction methods, especially for the overlapping or proximity of primaries and multiples. Unlike the SRME, the ISS provides an accurate travel time, amplitude, and phase of surface-related and internal multiples. This gives an opportunity to deal with complex interfering situations of different events Zou et al. 2019). Van Der Neut et al. (2011) studied the application of seismic interferometry in the removal of ghost, multiples and blur effects. Meles et al. (2015) and da Costa  integrate the Marchenko equation and seismic interferometry to predict prestack internal multiples. The third class is the machine learning methods concentrating on mining features from seismic data through well-trained neural networks like the generative adversarial network (Tao et al. 2022), the U-Net network , the convolutional neural network (Li and Gao 2020;Liu et al. 2022) and the deep neural network . Once the networks are successfully trained, multiples can be removed efficiently.
In recent years, researchers have come around to viewing multiples as efficient signals and invested great effort to use the intrinsic geological information buried in multiples to image the underground structures, especially for the less illuminated regions. It should be noted that although multiple imaging becomes a hot topic, it does not mean that imaging methods play a more important role than elimination methods. In fact, the removal and use of multiples are of equal significance and both attempt to image primaries. The former is to get a less contaminated image of recorded primaries, and the latter is to get a more illuminated image of unrecorded primaries extracted from recorded multiples. So far, there are several imaging strategies of multiples according to the different theoretical backgrounds.
Firstly, the WRW-based strategies. Berkhout and Vershuur (1994) show that the migration/imaging of multiples can be achieved by specifying the total response as the source wavefield. Later, the focal transform was developed using a multichannel weighted cross-correlation to convert multiples into primaries (Berkhout and Verschuur 2003;Verschuur and Berkhout 2005). In other words, it transforms multiples of a specific order into multiples of the previous adjacent order.
Secondly, the Green-function-based strategies. Seismic interferometry converts multiples into quasi primaries for imaging (Claerbout 1968;He et al. 2007;Schuster et al. 2004;Shan 2003;Wapenaar 2004). To overcome the limitation of seismic interferometry, researchers developed the Marchenko imaging that can retrieve the Green's function from reflection data acquired only at the surface without further details about the medium through an iterative or direct inverse solution of the Marchenko equation (Broggini et al. 2012;Thorbecke et al. 2017;Wapenaar et al. 2013. Slob et al. (2014) offered an imaging method of internal multiples through solving the Marchenko-type equations, assuming a planar wave propagating in a 1D model. Behura et al. (2014) discussed the advantages and drawbacks of the autofocus imaging algorithm for primaries and internal multiples. After then, Marchenko imaging methods were extended to seismic data containing surface-related multiples, allowing us to retrieve the Green's function under a free surface condition (Singh et al. 2015(Singh et al. , 2017. Meanwhile, the elastic Marchenko methods also appeared (da Costa Filho et al. 2014Wapenaar 2014).
Thirdly, wave-equation-based strategies. The methods of this kind vary according to the different treatments of the source and receiver wavefield. Guitton (2002) and Liu et al. (2011) specified the total response, including primaries and multiples, as the source wavefield and the predicted multiples as the receiver wavefield. Similarly, to skip the prediction of multiples, Wang et al. (2014) replaced the source wavefield with the total response and a point wavelet and then replaced the receiver wavefield only with the total response. These methods deal with multiples of all orders but produce plenty of non-negligible crosstalk artifacts, blurring the final image of substructures ). For such a reason, imaging methods using order-separated multiples appeared. Liu et al. (2016) realized the leastsquare reverse time migration of order-separated multiples divided by a modified SRME algorithm. This method can remarkably reduce order-related artifacts but is ineffective in event-related multiples. Zhang et al. (2019) then realized the reverse time migration using order-separated water-bottomrelated multiples, where only the water-column multiples of two adjacent orders are cross-correlated, effectively reducing the event-related multiples. After then, the RTM and LSRTM methods of order-separated multiples are pushed into viscoacoustic media (Qu et al. 2020a, b). The multiple imaging methods of this kind, meaning the imaging of surface-related multiples of adjacent orders, is consistent with the viewpoint of extracting unrecorded primaries from recorded multiples (Weglein 2016(Weglein , 2018b. In other words, taking the (n-1)th order multiples as areal sources and the nth order multiples as receiver data in imaging methods, e.g., RTM, imply subtracting the travel time the former experienced and then obtaining the unrecorded primaries excited from the areal sources (receivers) and recorded at receivers.
As a new hot topic, multiple imaging remains problems. One is the crosstalk of nonadjacent order multiples, which is also known as artifacts caused by unrecorded multiples. Another is the limitation of smooth migration velocity models (Weglein 2022). When using a smooth velocity model, recorded multiples should be removed for a better image of recorded primaries. It is similar for the imaging of unrecorded primaries.
To build a solid theoretical base for the elimination or the imaging methods discussed before, researchers should pay more attention to the modeling or construction of multiples. In marine geophysical exploration, the water surface is appeared as a strong reflector generating plenty of waterbottom multiples or water reverberations. As the water layer is often known, water-bottom multiples can be modeled and predicted through a combination of wavefield extrapolation and bottom reflectivity estimation (Wiggins 1988). Innanen (2009) extended the Born series to satisfy the heterogeneous media and proved that the non-linear forward modeling of multiples containing an infinite number of scattering terms could be simplified through the inverse scattering method if primaries have been constructed. Unlike the wavefield extrapolation methods and the Born series methods, Berkhout (2014a, b, c) presented a novel operator-driven modeling and imaging framework, among which the full wavefield modeling (FWMod) uses local propagation and scattering operators and a reflectivity model (operator-driven modeling) instead of the property description of velocity and density (property-driven modeling). Constructing orderseparated multiples from seismic data directly without any simulation is an alternative way established via modifying the SRME or focal transformation (Li et al. 2018;Liu et al. 2016).
The four kinds of methods have limitations: the wavefield extrapolation method only considers the water reverberations, the Born modeling method is limited to the first order approximation, FWMod suffers great computational cost from the round-trip scheme, and the SRME-like methods share the same deficiency as their prototype. Hence, we propose an iterative finite-difference modeling method to simulate the acoustic surface-related order-separated multiples. The method is established on the areal/virtual source assumption. The surface-related multiples of the current order are injected into the modeling procedure to generate the next order multiples, among which a Ricker wavelet is injected to generate primaries. The perfectly matched layers (PML) boundary is used to avoid waves reflecting into the media underground, making it possible to simulate multiples of individual order. To replace a free surface with an absorbing surface, we need to add the source and receiver ghost in the iterative modeling procedure of surface-related orderseparated multiples. Furthermore, two reflection operators are presented and acted on the areal source as an approximation of the reflection at the free surface. Using the ghost creation approach and the reflection operators, we can simulate ideal surface-related order-separated multiples. At last, we test the validity of the two proposed reflection operators and the forward modeling method of surface-related orderseparated multiples in numerical experiments.

Free surface
The free surface boundary is the key factor in generating surface-related multiples, and its accurate representation determines the simulation accuracy of surface-related multiples. Here, we follow the stress image method introduced by Levander (1988) and then discussed by Graves (1996). This high-precision stable approach can be implemented straightforwardly in 2D acoustic media. The acoustic wave equation traveling in a 2D isotropic media with constant density is where p, c, f represent the acoustic pressure, velocity, and source function, respectively. The symbols and ∇ 2 denote the acoustic wave operator and the Laplace operator. , , t indicate the position of the observation and source point, and time, respectively. The pressure p vanishes at the free surface and is imaged from downward to upward of the surface: where k is the grid index of the free surface. i is the grid number away from the surface, which is specified as half of the order of the spatial finite-difference operator.

Dipole source approximation and ghost wave
The seismic wave from a source underneath the surface propagates downward to form a direct wave, propagates upward and reflects off the free surface to form a ghost wave. Suppose the source locates near the free surface. In that case, we are able to use a dipole source comprised of a direct wave and a ghost wave to convert the free surface boundary condition into a transparent boundary condition, which enables the forward modeling of surface-related order-separated multiples. As shown in Fig. 1, the direct wave can be viewed as being excited by the positive monopole (the original source). The ghost wave can be viewed as being excited by the symmetrical negative monopole. The dipole source approximation of the free surface is also valid for the areal source, as depicted in Fig. 2. The incident wave from the positive line and the reflected wave from the negative 3 Page 4 of 17 line join to form a dipole wave as an approximation of the free surface.
Seismic waves reflected off the free surface generate a source and receiver ghost. Typically, the ghost wave is removed as noise, but here we prefer to create the ghost wave in the proposed algorithms so that the simulated surfacerelated order-separated multiples are more approximative to those in the direct simulation with the free surface. The ghost creation approach is implemented by subtracting a specific piece of time in the ghost wave under the vertical propagation assumption (Cocher 2017).
where f is the source without a source ghost, f ghost is the source with a source ghost, p is the synthetic data obtained at receivers without a receiver ghost, p ghost is the synthetic data with a receiver ghost. The variables t, z s , z r , c represent time, the source position, the receiver position, and the velocity near the surface, respectively. The terms −f (t − 2z s c) and −p(t − 2z r ∕ c) are the source and receiver ghost propagating from the negative side to the positive side in the normal direction. Since the acquisition system is close to the free surface, we can neglect the propagation effect on the amplitude. According to Figs. 1, 2, and Eq. (3), the ghost wave is created through two approximations: the first is the dipole source approximation transforming the free surface into a transparent boundary using a dipole source/wave, and the second is the normal incident approximation shifting the ghost wave from the negative side to the positive side.

Reflection at the free surface
Reflection coefficients of seismic waves at a free surface are related to the incidence angle. Under the assumption of normal incidence, the coefficient is equal to − 1. Here we present two reflection operators, R 1 and R 2 , to describe the reflection effect at the free surface. By applying the operators to the incident wave, we can get a new areal source and then reinject it in the forward modeling process to obtain an accurate simulation of the surface-related order-separated multiples.
(1). The reflection operator R 1 is defined as a scaled mirror of the incident wave at the free surface, which is a critical modification of the source term frequently used in waveequation-based imaging methods of surface-related multiples, for example, Liu et al. (2011) and . We argue that in order to get accurate multiples, i.e., accurately simulating the reflection at a free surface, the additional scalar z r must be added to the source term and the ghost creation approach must be incorporated in the forward modeling process. The derivation of the scalar z r will be further explained in the following part of R 2 .
where the acoustic wave equation acts on the pressure wavefield p , and the Dirac function indicates the excitation of the areal source. , , t denote the areal source position, the receiver position, and time, respectively. The subscripts surf and src stand for the free surface and the areal source. z r is the depth of the areal source or receivers.
(2). The reflection operator R 2 is defined as a scaled time derivative of the mirrored incident wave. One can refer to the work accomplished by Cocher (2017) for a further derivation of this operator based on Born theory, in which the reflectivity located at the free surface is described by a Dirac function rather than a Heaviside function.
Diagram of the dipole source approximation for reflection at the free surface. Note that the direct wave can be regarded as being excited by the positive monopole (the original source), and the ghost wave can be regarded as being excited by the negative monopole symmetrical about the free surface. The direct wave + ghost wave forms a reflection at the free surface where the symbols are the same as Eq. (4). Here, we give another straightforward explanation of the reflection operator R 2 . We discrete the additional term 2 ∕ c t as follows: where the symbols are the same as Eq. (3). This is the origin of the scalar z r . Compared to Eq.
(3), we can conclude that the reflection operator R 2 converts the incident wave into a dipole wave which is composed of the incident wave from the positive side and the time-shifted reflected wave from the negative side as depicted in Fig. 2. In this situation, the mirror dipole wave is reinjected into the modeling process as a new areal source. Hence, it proves the validity of the dipole source approximation as well as the equivalence between the two operators R 1 and R 2 in theory.

Forward modeling of surface-related order-separated multiples
Both the reflection operators are proposed in the view of the Huygens-Fresnel principle or areal shot. It means that each point of the spherical wave front can be regarded as the secondary source, which will be superimposed at an observation point to form a total wavefield according to their respective contributions. The reflection operators and the ghost creation approach are used in a finite-difference forward modeling process to obtain the surface-related order-separated multiples.
With the dipole source approximation and the normal incidence approximation mentioned before, we can replace the free surface with an absorbing boundary so that the incident wave arrived at the free surface are avoided propagating back to the subsurface again. The (n-1)th order multiples recorded at the receivers are mapped to a new areal source by the two reflection operators, then the areal source participates in an additional round trip of the finite-difference forward modeling generating the nth order multiples. Primaries and the surface-related order-separated multiples are calculated using Eq. (7) for the reflection operator R 1 : where represents the acoustic operator. f , f ghost represent the source function without a ghost (monopole source) and with a ghost (dipole source). p i=0,1,2 , p i=0,1,2,ghost are the ith order multiples without and with a receiver ghost, respectively. Other symbols are the same as the previous equations. The wave p i,ghost in Eq. (7) is the desired surfacerelated order-separated multiples. It combines the incident wave p i (t) and its receiver ghost p i t − 2z r ∕ c to form a new dipole wave required by the reflection operator R 1 to approximate the free surface. That is to say, the receiver ghost is a necessary part of the approximation of the free surface and must be included in the areal source. Hence, we can simulate the surface-related order-separated multiples iteratively by using an absorbing boundary. Similar to the equation for R 1 , primaries and the surfacerelated order-separated multiples can be calculated using Eq. (8) or Eq. (9) for the reflection operator R 2 : Alternatively, we have derived a novel equivalent form: Because the reflection operator R 2 can convert the incident wave into the dipole wave, there is no need to be provided with the receiver ghost once again. Equation (9) has the same form as Eq. (7), but the areal source is selected as the incident wave with only the source ghost. Thus, we need to supply the receiver ghost individually. In summary, only the dipole wave, including the source and receiver ghost, is valid for the reflection operator R 1 , while both the dipole wave and the incident wave with only the source ghost are valid for the reflection operator R 2 . Meanwhile, higher order multiples are available straightforward according to the above principle.

Algorithm implementation
We only consider the algorithms described in Eq. (7, 8), which will be tested in the numerical experiment section. The implementation of Eq. (9) is quite similar to that of Eq. (7) and has been verified in private. It will not be presented in this article for simplicity. The basic steps of the surface-related order-separated multiples forward modeling method adopting the reflection operator R 1 are listed as follows: (1) Simulating the entire synthetic data under a free surface condition, including multiples of all orders.
(2) Simulating primaries containing the source ghost only. Replace the free surface with the absorbing boundary, construct the dipole source using Eq. (3), and inject the (9) dipole source in the forward modeling procedure to obtain primaries with only the source ghost.
(3) Receiver ghost creation. Supply the receiver ghost to primaries that contain only the source ghost using Eq. (3).
(4) Constructing the areal source and simulating the first order multiples. Apply the reflection operator R 1 to primaries with the source and receiver ghost and reinject it to the forward modeling procedure according to Eqs. (4, 7). Note that two details need attention. First, the areal source produces great noise owing to the end effect, which can be attenuated by a specific damping window. Second, the discrete Dirac function is defined as 1∕ Δh where Δh is the grid interval in space (Mittet 1994). Then we get the discrete equation: the subscripts i, j and superscript n of the wavefield p stand for the grid index of depth, distance and time, respectively, which is discretized by (Δh, Δt) in space and time. ir is the grid index of the receiver depth where the areal source locates. The symbol ∇ 2 d denotes the discrete Laplace operator. A Blackman damping window w ir,j is applied to the dipole source s ir,j (primaries with the source and receiver ghost). Other symbols follow the same definitions as the equations before. We can obtain the first order multiples containing only the source ghost through Eq. (10).
(5) Time shift correction. As illustrated in the left part of Fig. 3, the incident wave traveling to and reflecting away from the free surface experiences double the dept 2z r between the free surface and the receiver line. In the R 1 case, as shown in the middle part of Fig. 3, the additional receiver ghost together with the incident wave, converts the incident wave (primaries with only source ghost) to a dipole source which can be viewed as being obtained/excited at the free surface. Once the new dipole areal source is reinjected at the receiver line, the travel path is z r less than that of the free surface situation. Hence, the travel time of the synthetic multiples must be corrected by z r ∕ c . This constant is valid for all-order multiples except for primaries.
(6) Simulating the second or higher order multiples. Equation (10) is valid for the forward modeling of allorder multiples. The surface-related multiples of the current order can be calculated by the repetition of steps (3), (4), and (5), taking the multiples of the previous adjacent order as the input.
For the reflection operator R 2 , the basic steps of the surface-related order-separated multiples forward modeling method are mostly the same as the R 1 case except that: ir,j = p n+1 ir,j − c 2 ir,j Δt 2 s ir,j ⋅ w ir,j z r Δh Fig. 4 The homogeneous model. The asterisk is the source, and the dashed horizontal line is the receiver line (1) There is no need to add the receiver ghost to multiples.
(2) The discrete equation for the forward modeling of surface-related order-separated multiples is defined as: where d is the discretized derivative and s ir,j stands for the incident wave with the source and receiver ghost.
(3) Time shift correction. In the R 2 case, the incident wave obtained and reinjected at the receiver line leads to (11) the travel path being 2z r less than that of the free surface situation, according to Fig. 3. Hence, the travel time of the synthetic multiples must be corrected by 2z r ∕ c , which is also valid for all-order multiples except for primaries.
Overall, the treatment of the ghost wave, the areal source, and the time shift are critical to getting an ideal result. Besides, the direct wave of primaries and multiples must be removed when realizing the algorithms.

Numerical examples
To verify its correctness and effectiveness, we test the surface-related order-separated multiples forward modeling method in three models (the homogeneous model, the twolayer model, and the Pluto model). The forward modeling method is a finite-difference scheme with second order in time and eighth order in space. First, the upper surface is set as the free surface. Then we obtain the entire synthetic data including all order multiples. Secondly, the upper surface is set as a perfectly matched layer (Berenger 1994;Komatitsch and Tromp 2003), then both the reflection operators are used to replace the free surface to simulate the surface-related Fig. 8 The two-layer model. The asterisk is the source, and the horizontal dashed line is the receiver line Fig. 9 The comparison of the synthetic shot data for the two-layer model. (a) The data simulated with the free surface; (b) The summation of primaries and surface-related multiples simulated with the reflection operator R 1 ; (c) The relative error; (d-f) Primaries, the first order multiples, and the second order multiples by operator R 1 , respectively Fig. 10 Same as Fig. 9, but for the comparison between the free surface and R 2 shot data Fig. 11 Synthetic seismogram comparison between the free surface and R 1 solutions selected at the dashed lines in Fig. 9a, b. a Trace at a distance of 650 m; b Trace at a distance of 800 m; c Trace at a distance of 950 m Fig. 12 Same as Fig. 11, but for the comparison between the free surface and R 2 solutions selected at the dashed lines in Fig. 10a, b order-separated multiples. Finally, the relative error is introduced to analyze the precision of the two operators quantitatively through Eq. (12) where p fs and p R i represent the entire recorded seismograms obtained with the free surface and the two reflection operators, respectively.

The homogeneous model test
As shown in Fig. 4, the homogeneous model (1500 m/s) aims to study the approximation precision of the two reflection operators, R 1 and R 2 , to the free surface. The model is a rectangular region of 1000 × 600 m discretized with a 5 m grid spacing. The source, situated at (500 m, 400 m) and shown by an asterisk, is a 30 Hz Ricker wavelet. The receivers distribute with a 5 m spacing along the dashed line at a depth of 200 m. The wavefield is sampled every 0.5 ms until 1000 ms.
As the source is buried deep away from the surface, the wavefield is not affected by the source ghost allowing full concentrations on the validation of the reflection operators themselves. Figure 5a-c exhibits the synthetic data simulated Fig. 13 The part Pluto model. The asterisk is the source, and the horizontal dashed line is the receiver line (d-f) Primaries, the first order multiples, the second order multiples using the free surface and reflection operators R 1 and R 2 , respectively. The three dashed lines represent three selected traces located at distances of 650 m, 800 m, and 950 m, which are exhibited in detail in Fig. 6 for R 1 and Fig. 7 for R 2 in a sequence of a, b, and c. The first event in Fig. 5 is the direct wave emitted by the source and received by the receivers. The second one in the same figure is the reflected wave from the surface to the receiver dashed line. Figure 5d, e shows the relative error for the reflection operators R 1 and R 2 , respectively. The reflected waves are excellently simulated with R 1 and R 2 at a quite low degree of relative error. As depicted in Fig. 5d, e, the simulation precision of the two operators decreases with the offset apart from the source in general due to the normal incidence assumption, while the direct wave acquired at receivers is subtracted entirely. Figures 6, 7 further illustrate the effectiveness of both the reflection operators R 1 and R 2 by overlapping the traces accordingly. The direct waves (the first event) match perfectly because the surface has no effect on the sources and receivers, and the reflected waves (the second event) match almost perfectly owing to the excellent approximation of the two reflection operators. The reflected waves generated by R 1 are slightly less accurate than those generated by R 2 , which can be attributed to the individual supplement of the receiver ghost.

The two-layer model test
Unlike the homogeneous model, the two-layer model is designed to generate waves with both a source and receiver ghost, as depicted in Fig. 8. The model distributes 400 m in depth and 1000 m in horizontal distance, which is discretized with a 5 m grid spacing. The source using a 30 Hz Ricker wavelet is indicated with an asterisk and positioned in the surface center. Both the source and receivers distribute Fig. 15 Same as Fig. 14, but for the comparison between the free surface and R 2 shot data 5 m underneath the surface. The wavefield is sampled every 0.5 ms until 1000 ms.
Primaries and the first two order multiples are simulated in an iterative finite-difference process. Higher order multiples follow the same scheme and are neglected here. The synthetic data produced with the free surface and reflection operator R 1 are depicted in Fig. 9, including the total wavefield (Fig. 9a, b for the free surface and the R 1 cases, respectively), the relative error (Fig. 9c), and the order-separated multiples (Fig. 9d-f). Figure 10 follows the same arrangement. Primaries and the first two order multiples distribute independently along the time axis. As shown in Figs. 9c, 10c, the relative error for the reflection operators R 1 and R 2 mainly concentrates on primaries and the offset far away from the source. It can be attributed to the fact that both the reflection operators are derived from the normal incidence assumption, which is not fully satisfied for primaries generated from a point source but is nearly perfectly satisfied for multiples generated from an approximately planar areal source (the previous order multiples).
We extract three traces at distances of 650 m, 800 m, and 950 m from the synthetic data in Figs. 9, 10 to show the travel time, phase, and amplitude in detail in Figs. 11,12. For the trace near the source, the data simulated with the free surface and the two reflection operators matched perfectly in the travel time, phase, and amplitude. At the same time, the amplitude misfit between the different data mainly exists in primaries and increases with the offset apart from the source. In contrast, the first and second order multiples in the simulated data match almost perfectly because they are more satisfied with the normal incidence approximation than primaries. We perform three iterations herein, and the results are shown in Figs. 14d-f, 15d-f for the reflection operators R 1 and R 2 , respectively. Figures 14b, 15b are the summations of primaries and the order-separated multiples to compare with the total record simulated by the FD algorithm with the free surface in Figs. 14a, 15a. It can be observed that the total record generated by the FD algorithm with the free surface and both the reflection operators keep remarkable consistency. As depicted in Figs. 14c, 15c, we can observe the same phenomenon as the results of the previous models. The relative error increases with the offset apart from the source and decreases with the increasing order of multiples (if taking primaries as the zeroth order multiples), which can be attributed to the normal incidence approximation relating to the offset and orders.
Then traces marked with the dashed lines at a distance of 2050 m in Figs. 14, 15 are selected and presented in Figs. 16, 17, respectively. The consistency between the records generated by the FD algorithm using the free surface and the two Fig. 17 Same as Fig. 16, but for comparison between the free surface and R 2 solutions. Note that the corresponding traces are selected in Fig. 15 3 Page 14 of 17 reflection operators is much clearer (Figs. 16a, b, and 17a, b). As the iteration continues, multiples become increasingly weak and complicated . Different orders of multiples blend to form the ultimate waveform of the record simulated by the FD algorithm with the free surface and the reflection operators, which will bring significant difficulties for the multiple elimination and imaging. For details, we intercept three different time segments and display them in Fig. 18 for the reflection operator R 1 and in Fig. 19 for the reflection operator R 2 . The records show an excellent agreement in travel time, phase, and amplitude. The slight amplitude mismatch in Fig. 18b, c is still attributed to the fact that the individual supplement of the receiver ghost for the reflection operator R 1 introduces an additional normal incidence approximation when compared to the reflection operator R 2 incorporating the receiver ghost inherently. As we can see, for complex models like Pluto, in the majority of records, primaries and multiples of different order always mix up to form the final waveform, which will bring more difficulties in multiples elimination and imaging. The numerical results of the part Pluto model demonstrate that both the proposed reflection operators are capable of approximating primaries and order-separated multiples nearly perfectly for complex understructures.

Conclusions and discussions
The three numerical tests conform to our theoretical expectations and prove that the two proposed reflection operators and the surface-related order-separated multiples forward modeling method are valid and feasible. The entire process shows a nice agreement with the data generated with the free surface. Both the reflection operators are excellent approximations of the free surface. The precision decreases with the offset away from the source and increases with the order of multiples owing to the normal incidence assumption. Since the reflection operator R 1 introduces an additional normal Fig. 18 Synthetic seismogram comparison between the free surface and R 1 solutions for three different segments in Fig. 16a, b. (a) 1.5-2.8 s; (b) 3.4-4.5 s; (c) 5.0-6.5 s. Different wave types are marked with 'P' for primaries, 'm 1 ' for the first order multiples, and 'm 2 ' for the second order multiples incidence approximation when adding the receiver ghost individually, it is slightly less accurate than the reflection operator R 2 incorporating the receiver ghost inherently. The proposed method affords an opportunity to model ghost-free synthetic data directly without any deghost procedure, which will be discussed in the future.
Meanwhile, researchers should notice two limitations caused by the ghost creation approach. The surface of models is required to be planar, and the velocity close to the surface must be constant. In addition, the depth of sources and receivers must not be too far from the surface.

Declarations
Competing interests The authors have no competing interests to declare that are relevant to the content of this article.
Ethical approval Not applicable.

Consent for publication
The publication of this manuscript has been approved by all co-authors and the responsible authorities. It has not been published before, nor is it under consideration for publication anywhere else.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. Fig. 19 Same as Fig. 18, but for comparison between the free surface and R 2 solutions in Fig. 17a, b