## ABSTRACT

Seismic modeling in heterogeneous media is accomplished by using either approximate or fully numerical methods. A popular approximate method is ray-Born modeling, which requires the computation of 3D integrals. We have developed an integration technique for accurate and, under certain circumstances, efficient evaluation of the ray-Born integrals in the time domain. The 3D integrals are split into several 2D integrals, each of which gives the wavefield at a certain time, so that the waveform at each time step is computed independently of all other times. We compute seismograms for 3D heterogeneous acoustic media using this technique and compare these seismograms with seismograms computed using two other modeling methods: frequency-domain ray-Born modeling and finite-difference modeling of the acoustic wave equation. Our method can also be applied to elastic ray-Born modeling. Velocity models with smooth scatterers and the SEG/EAGE overthrust model are used for comparison. The ray-Born seismograms computed using the time- and frequency-domain ray-Born modeling methods are identical, as expected. The comparison between the ray-Born modeling and the finite-difference-modeling method indicates that the waveforms are similar for both types of velocity models. We evaluate the discrepancies in terms of multiple scattering and multipathing.

## INTRODUCTION

Modeling of seismic waves is very important for studying the Earth’s structure. For example, seismic algorithms such as full-waveform inversion require it. In practice, this modeling is often done using the acoustic wave equation, solved either by using fully numerical methods or by using approximate methods (e.g., Carcione et al., 2002; Virieux and Lambaré [2015], chapter 1.05). Fully numerical methods are popular and very useful, but they are not computationally efficient at higher frequencies, especially in large 3D models. They also have limitations in the case of complicated structural models (e.g., finite-difference seismograms contain artifacts in the case of strongly undulating interfaces). Approximate methods do not necessarily suffer from these disadvantages; for example, ray theory is particularly valid at high frequencies in smooth models. Therefore, it is sometimes beneficial to compute wavefields using approximate methods.

A useful technique for modeling singly scattered waves is the (first-order) ray-Born approximation. Several papers have investigated the theoretical aspects of the Born (e.g., Knopoff and Hudson, 1964; Hudson, 1977; Hudson and Heritage, 1981) and ray-Born approximations (see references listed by Chapman [2004], section 10.3 and Červený [2005], section 2.6.2). In studies focusing on the numerical aspects of ray-Born modeling, Červený and Coppoli (1992) generate synthetic seismograms in the case of a background medium with two curved interfaces and a perturbation consisting of isolated scatterers. For a model with similar characteristics, Gibson et al. (1993) use ray-Born modeling to examine reflections from vertical seismic profile data. Coates and Charrette (1993) compare the 2D isotropic elastic generalized Born approximation (Coates and Chapman, 1991) and the conventional Born approximation to finite-difference modeling. In a study focusing on 3D ray-Born inversion, Thierry et al. (1999) compute synthetic traces from a migrated image and compare the results with real data. In a review paper, Moser (2012) pays special attention to the advantages of acoustic ray-Born modeling for modeling diffractions from specific structural discontinuities, such as edges and tips. Šachl (2013) investigates discretization effects of the elastic ray-Born scattering integral and suggests a cosine window to attenuate boundary diffractions. Tengesdal (2013) and Minakov et al. (2017) discuss 2D acoustic ray-Born and finite-difference modeling and show synthetic- and real-data examples as well as an application to full-waveform inversion. All of the above papers compute the ray-Born integrals, in either the time or frequency domain, by summing over individual volume or area elements (e.g., Coates and Charrette, 1993). An alternative method for the computation of ray-Born seismograms in the time domain was presented by Spencer et al. (1997) and Keers et al. (2002). That method used a specific integration over tetrahedrons.

We present a new method for 3D ray-Born modeling in the time domain using isochron surfaces. Integrals over isochrons are used to accurately compute seismograms for a band-limited signal. We compare the results with seismograms computed using frequency-domain ray-Born and finite-difference modeling. The method can be seen as an extension to 3D integrals of what was done for 1D Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) seismograms by Chapman (1978), where the band limiting is done before the discretization and for 2D Kirchhoff seismograms by Haddon and Buchen (1981), where the integrals are computed over isochrons.

The results in this paper are based on the acoustic wave equation, but the concepts can be extended to the elastic case. We give a brief description of this in the “Discussion” section.

In the “Theory” section, we review the acoustic ray-Born modeling and give the theoretical background of the integration technique. Next, we discuss the numerical implementation of the ray-Born integrals. For the sake of completeness, this section also gives a brief description of the two-point ray-tracing algorithm as well as the frequency-domain ray-Born integrals and finite-difference modeling. In the “Results” section, we present seismograms computed using two types of velocity models with different scattering characteristics. The first type has smooth velocity perturbations, analogous to that used in ray-based tomographic model building (or full-waveform inversion with low-frequency velocity updates). The second type, the 3D SEG/EAGE overthrust model of Aminzadeh et al. (1997), contains rough and rapidly varying perturbations.

## THEORY

### Background Green’s function

### The ray-Born approximation

Later on in this paper, we briefly discuss multiple scattering and multipathing, and it is important to distinguish the two effects. Multipathing occurs when, given a velocity model $c$ and source $s$, there are two or more rays between $s$ and the receiver $r$. Multiply scattered waves form that part of the wavefield that is modeled using the higher order terms ($L2g0$, $L3g0$, etc.) of the scattering series 9. Therefore, multipathing depends on the source $s$, the velocity model $c$, but also the background velocity $c0$. Singly scattered waves are the waves modeled by the term $g1=Lg0$. Note that multipathing and multiple scattering are not mutually exclusive, so that a part of the wavefield can be interpreted in terms of multipathing and in terms of multiple scattering.

### Discrete ray-Born seismograms in the time domain

For $t=t1$, it is not possible to integrate along an isochron because the isochron for $t=t1\u2212\Delta t$ does not exist. In this case, there are two different cases that must be distinguished: (1) $t1\u2212\Delta t<Trs<t1\u2212\Delta t/2$ and (2) $t1\u2212\Delta t/2<Trs<t1$. In both cases, the ray-Born integral is given in terms of a volume integral as in equation 25. In our experience, these contributions in practice are not significant. However, this should be included here for the sake of completeness. The numerical implementation of equation 36 and the volume integrals are discussed in the next section.

## NUMERICAL IMPLEMENTATION

### Ray tracing

In preparation for ray-Born modeling using equation 36, the quantities $Trs$, $\u2207Trs$, and $Ars$ are computed using ray tracing. They are found by first computing $T(s,x)$, $\u2207T(s,x)$, and $A(s,x)$ from the source $s$ to all points $x$ in the scattering region $D$. Likewise, $T(r,x)$, $\u2207T(r,x)$, and $A(r,x)$ are computed from the receiver $r$ to the scattering points $x$.

If the velocity model is slowly varying, raypaths, traveltimes $T$, and wavefront normals $\u2207T$ are found by solving the kinematic ray equations given the initial conditions (position and take-off angles). In a similar manner, the amplitude is computed by solving the dynamic ray equations. Both sets of equations are solved numerically with a fourth-order Runge-Kutta method (Press et al., 2007) with Snell’s law applied when the ray encounters an interface. The result is two 3D fans of rays emanating from the source and the receiver.

For a given source $s$ and scattering point $x$, the two-point ray-tracing problem is to find a ray from $s$ to $x$. This is done in two steps. First, all points on the ray fan are used to construct a 3D Delaunay tessellation (de Berg et al., 2000), which covers the scattering region $D$. For each $x$ in $D$, we then interpolate values for $T$, $\u2207T$, and $A$ using the tessellation and linear interpolation on a tetrahedron (e.g., Farin, 2002). Likewise, the same procedure is done for a receiver $r$ and scattering point $x$.

### Ray-Born seismograms: Time domain

The isochron can be computed in various ways. In a heterogeneous background with the traveltime function $Trs$ given on a discrete grid, one of several isosurface extraction techniques can be used (e.g., Sutton et al., 2000). We have computed $Trs$ using the ray-tracing procedure described above, and we have used the marching cubes algorithm (Lorensen and Cline, 1987) to extract surfaces. In this case, $S\tau $ is defined on a triangular grid with vertex positions $x1=x1(n)$, $x2=x2(n)$, and $x3=x3(n)$, where $n=1,\u2026,ntr$ and $ntr$ is the number of triangles.

The volume integral 25 for the two special cases mentioned in the previous section can also be evaluated numerically. For this, the vertices of the first isochron as well as the points on the background ray, which is completely enclosed by the isochron, can be used. In this case, the integral 25 becomes a sum of the volume contributions of the tetrahedrons $\u2211n=1ntetI^n\Delta Vn$, where $ntet$ is the total number of tetrahedrons, $\Delta Vn$ is the volume of tetrahedron $n$, and $I^n$ is the integrand of tetrahedron $n$ (which may be computed as the average of the integrand on the vertices of each tetrahedron).

### Ray-Born seismograms: Frequency-domain comparison

The computation of the ray-Born approximation in the frequency domain is done using equation 17. The discretization of the integral is accomplished using a regular grid over the scattering region $D$. By approximating the integrand in each grid element by a constant, the computation of $g1(r,s,\omega )$ becomes a sum over these integrand values, multiplied by the volume element.

The comparison with the time-domain ray-Born modeling method and finite-difference modeling is done in the time domain. Therefore, after computation of $g1(r,s,\omega )$ for all frequencies, an inverse fast Fourier transform is used.

### Finite-difference comparison

The finite-difference method (e.g., Etgen and O’Brien, 2007) is used to solve the wave equation 1. For our implementation, we do not include any absorbing boundary conditions. To avoid edge effects within the time window of interest, all velocity models are extended beyond the area of interest.

## RESULTS

In this section, seismograms computed using the three modeling methods (time-domain ray-Born using isochrons, frequency-domain ray-Born, and finite-difference) discussed in the previous section, are presented for two different velocity models. The first of these velocity models is a smooth model with random reflectivity. The second one is the 3D EAGE/SEG overthrust model (Aminzadeh et al., 1997). This is a complicated velocity model with a rapidly varying velocity. Small differences between the computed seismograms are assessed using single-trace displays. Larger scale differences, with significance for structural interpretation, are assessed using shot records and zero-offset sections.

### 3D random Gaussian model

First, we present the modeling results for the velocity models with random variations. These random velocities are generated with a Gaussian autocorrelation function characterized by a perturbation strength $\u03f5$ (root-mean-square slowness variation) and correlation distance $a$ (Baig and Dahlen, 2004). We use this to create a random realization with weak, medium-strength, and strong perturbations (details are given in Table 1). The correlation distance in all cases is $a=50\u2009\u2009m$ and the velocity model covers a volume of $1.8\xd72.4\xd72\u2009\u2009km$. A constant background velocity of $2000\u2009\u2009m/s$ is used. Figure 3 shows two cross sections of the particular realization that was used.

The model grid spacing is 10 m, and this is used for the ray-Born-modeling methods as well as the finite-difference-modeling method. This way, amplitudes, traveltimes, and velocities are all sampled on the same grid. The models are extended in all directions to avoid interference from artificial boundary reflections within the recording length of 2 s. The sampling interval used was 4 ms. Receivers are positioned at a 50 m interval from 0 to 2350 m in the $y$-direction and at depth of 100 m (Figure 3). Synthetics are generated for a shot located at $s=(800,1200,200)\u2009\u2009m$, using a source time function in the form of a causal Ricker wavelet, with a center frequency of 15 Hz.

For three discrete sample points in time, $ti=(0.65,0.85,1.05)\u2009\u2009s$ and a source-receiver distance of 800 m, Figure 4 shows the corresponding isochron surfaces $Trs=ti$. The values on the surfaces are the values of the integrand of equation 36 for the weak perturbation model. As the background medium is homogeneous, the isochron surfaces are part of a spheroid.

Figure 5a (0.5% perturbation) shows a trace offset by 800 m relative to the source position (see the dashed blue line in Figure 6 for the location of the receiver) for the weak perturbation model (Table 1). We show the finite-difference seismogram and ray-Born seismograms computed in the time and frequency domains. Next, we show traces for the same random realization, but with medium-strength and strong perturbation models in Figure 5b (10% perturbation) and 5c (15% perturbation), respectively. For the models with 0.5%, 10%, and 15% perturbation, the ray-Born seismograms all appear with the same phase, but with a relative difference in the amplitudes. The results computed in the time and frequency domain are basically identical, as expected. In the comparison with finite-difference modeling, the seismograms appear with changes in the waveforms for the medium-strength (Figure 5b) and more so for the strong perturbation model (Figure 5c). This is because finite-difference modeling accounts for complicated wave-propagation effects such as multiple scattering and multipathing that become more predominant for stronger perturbations.

Figure 6 shows synthetic shot gathers for the medium-strength perturbation model. The synthetics generated by the ray-Born method and finite differences are shown in Figure 6a and 6b, respectively. Across the offset positions, we can recognize events from the singly scattered wavefield in both seismograms.

### 3D SEG/EAGE overthrust model

In this subsection, time-domain ray-Born, frequency-domain ray-Born, and finite-difference synthetics are computed using the acoustic 3D SEG/EAGE overthrust model (Aminzadeh et al., 1997; Operto et al., 2003). The overthrust model represents a succession of sedimentary deposits in a compressional tectonic setting, resulting in strong vertical and lateral velocity variations. It was created to study 3D wave propagation and to test seismic processing algorithms.

Two minor modifications are made to the original model: Edges are extrapolated to avoid artificial boundary reflections, and a seafloor interface is added at a depth of 500 m using a water layer velocity of $1480\u2009\u2009m/s$. A subvolume of the larger model is displayed as a vertical cross section in Figure 7a and at a depth of 2 km in Figure 7b.

We create the background model by 3D boxcar averaging of the model without the water layer. The background model is shown in Figure 7c and 7d. The difference $c\u2212c0$ then gives the perturbation model, which, in Figure 7e and 7f, is displayed in terms of percent relative to the background. Perturbation strengths have a maximum of 41.1% and a minimum of $\u221221.7%$. On average, the positive velocity perturbation is 6.8% and the average value of the negative velocity perturbation is $\u22125.1%$. In contrast to the velocity models of the previous section, the perturbation model contains sharp velocity variations, and, therefore, it is quite rough.

The grid spacing is 12.5 m along all spatial axes, and it is the same for frequency-domain ray-Born-modeling and finite differences. To achieve an exact match between time-domain ray-Born modeling and frequency-domain ray-Born modeling, we reduced the vertical grid spacing to a 6.25 m space for the time-domain implementation. This allows for a more precise isochron surface extraction. The source-time function used for the computation of seismograms is defined by a causal Ricker wavelet with a center frequency of 10 Hz. In the case of the ray-Born modeling methods, we include the direct wave traveling in the water layer, but we do not calculate the seafloor reflection.

With a water layer, take-off angles for ray tracing are limited by the critical angle at the water-sediment interface. Therefore, we do not have any real valued raypaths beyond these angles. This limited ray coverage creates an artificial boundary in the scattering volume, with an associated end-point signal (Bender and Orszag, 1999). To avoid these artifacts, traveltime and amplitude values are interpolated by including points just above the seafloor from ray tracing in the water layer, and adding these to the collection of spatial locations from one-point ray tracing before the Delaunay tessellation process. Figure 8 shows traveltimes and amplitudes for a source at $s=(2500,5000,200)m$ after interpolation on the uniform grid. Figure 8 also shows the interpolated ray shadow zone, the triangular wedges between the dashed line and the seafloor (the solid line).

For a source positioned at $s=(2500,5000,200)\u2009\u2009m$ and a receiver positioned at $r=(2500,5000,100)\u2009\u2009m$, Figure 9 shows isochron surfaces in the overthrust background model for times $Trs=1.8\u2009\u2009s$ (Figure 9a) and $Trs=2.36\u2009\u2009s$ (Figure 9b). Values on the surfaces are computed with the integrand in equation 36. Because the background model is heterogeneous, the isochron surfaces look like deformed spheroids as expected.

#### Ray-Born and finite-difference zero-offset sections

We show an example of zero-offset synthetics computed across a region with complicated faulting as well as a structurally simpler area with near-horizontal layers. Sources (at a depth of 100 m) and receivers (at a depth of 200 m) are positioned from 500 to 5000 m in the $x$-direction and at 5000 m in the $y$-direction. This is indicated by the dotted black line in Figure 7b. With an $x$, $y$ aperture of 4000 m, the wavefield is computed using the two ray-Born-modeling methods and the finite-difference-modeling method.

Seismograms, computed using finite-difference and ray-Born modeling, are shown in Figure 10, for three positions in the overthrust model. The three positions are at 1000, 2500, and 4000 m, as indicated with the dashed blue lines in Figure 11. The seismograms in Figure 10a (1000 m) are extracted from a region with less complicated wave-propagation effects (the velocity model has horizontal and dipping layers). Here, the comparison with the finite differences shows that the phase and amplitude appear to be well-matched, indicating the relative importance of singly scattered waves. We also compare the seismograms shown in Figure 10b (2500 m). In this case, the source and receiver are located in a region where the velocity model has steeply dipping layers and faulted structures. These structures cause complicated wave-propagation effects (such as multiply scattered events) that appear as differences in the modeling results after traveltimes of approximately 1.4 s. In Figure 10c (4000 m), the early-time waveform mismatch in the time interval 1−1.4 s is caused by a triplication in the wavefield that is formed by the shallow syncline. Perturbations at the syncline are relatively strong (up to 30%) and are therefore not accurately captured by the ray-Born approximation.

A zero-offset section for the ray-Born synthetics is shown in Figure 11a and for finite differences in Figure 11b. Strong reflections from the horizontal and dipping layers are indicated in the blue rectangle labeled by I, confirming what was discussed above for Figure 10a. The wavefield becomes more complicated below the anticline, indicated by the blue rectangle labeled II. In this area, ray-Born synthetics exhibit stronger reflections than the finite-difference synthetics. This is likely an effect of multiple scattering and multipathing: Internal reflections in the finite-difference synthetic interfere with the main reflections, and these internal reflections do not exist in the case of ray-Born modeling. In the blue rectangle labeled III, we have indicated the triplication effect that causes the mismatch at early times in Figure 10c. This particular effect is present in the finite-difference synthetic, but poorly modeled in the ray-Born synthetic.

## DISCUSSION

The main results of this paper are the derivation and computation of ray-Born integrals as integrals over isochron surfaces. First, we give the continuous ray-Born integral 21 and then derive the band-limited ray-Born integral 36, which is used in the numerical implementation. In the derivation of the band-limited integral, we use the normal distance from the isochron $St$ to $St+\Delta t/2$ to reduce a volume integral to a surface integral. Then, we described a numerical integration over this surface in equation 40, where the integration is carried out over triangles.

The concepts discussed above are derived for the acoustic case, but they can be applied to the elastic case and also to the generalized ray-Born integrals (Chapman, 2004). For an elastic isotropic model, a natural starting point is to carry out the derivations from equations 21–36 using the elastic ray-Born integral (e.g., Dahlen et al. [2000], equation 56) as a starting point.

Waveforms computed using this new method are compared with ray-Born integrals computed on a regular grid in the frequency domain. This latter method is accurate as long as sampling intervals for the integration grid are kept small; however, relative to the new time domain method, the frequency-domain method is not as computationally efficient. This is because a summation is carried out over all volume elements in the scattering region and the summation is repeated for every frequency. The numerical efficiency of the time-domain isochron ray-Born method largely depends on the efficiency of the method used to compute the isochrons. If the isochrons are readily available, as they could be in a 1D background model, the integration over triangles is very cheap and the computational cost is insignificant compared with the frequency-domain method.

In the ray-Born integration method by Spencer et al. (1997) and Keers et al. (2002) (see also Chapman, 2004), the computations are made accurate in the time domain by having an integration grid that closely follows isochrons. This is important for modeling at higher frequencies because the isochron surface tends to have irregular curvature. However, compared with the ray-Born integrals using the isochrons explicitly, as discussed in this paper, the integration for the tetrahedron method is carried out over volumes and with a more involved implementation.

It should be noted that isochrons play an important role in many seismic methods because they associate the traveltime of a point in the subsurface with a source-receiver pair. Some examples are Kirchhoff imaging (Bleistein, 1987), the computation of Fresnel volumes (Červený and Soares, 1992), and finite-frequency kernels (Dahlen et al., 2000). Isochrons have also been studied by Iversen (2004), who presents a system of differential equations to compute isochron rays in a similar manner to how conventional rays are computed using the ray equations.

We also compare the new ray-Born modeling method with finite-difference modeling to assess the singly scattered wavefield as an approximation to the complete wavefield, where multiple scattering and multipathing are present. We broadly conclude that, for weak perturbations, or for regions where the wavefield to a large extent consists of singly scattered waves, the ray-Born approximation is the most successful. Based on its ability to produce multiple-free data, it is natural to consider the ray-Born approximation as a modeling method for full-waveform inversion (e.g., during early iterations when the background model is slowly varying) or for the purpose of testing various inversion and imaging algorithms. Furthermore, there are applications in target-oriented studies, for example, in the modeling of the seismic responses caused by small and localized scatterers.

## CONCLUSION

We have presented a physically intuitive method for computing 3D ray-Born integrals over isochrons in the time domain. The integrals are simplified because volume integrals are split in a natural way to several surface integrals. Each surface integral gives the waveform in the ray-Born integral for a certain time. The integration technique was used to compute ray-Born seismograms for two velocity models with different scattering characteristics: slowly and rapidly varying velocities. For both model types, the seismograms computed over surface integrals in the time domain are identical to ray-Born seismograms computed over volumes in the frequency domain.

## ACKNOWLEDGMENTS

We thank our colleagues at Schlumberger for their constructive comments that led to improvements of the manuscript. In particular, we thank J. Ø. H. Bakke, M. Nickel, and L. Sønneland. The work was funded by the Norwegian Research Council through the Petromaks2 program, project number #245587. This research was partly done while H. Keers was on sabbatical at the University of Kiel and the University of Washington, and support for this work (by the L. Meltzer Høyskolefond, the Faculty of Mathematics and Natural Sciences of the University of Bergen, the Institute of Geoscience at the Christian Albrechts University of Kiel, the Applied Physics Laboratory at the University of Washington, and especially hosts T. Meier [CAU] and B. Odom [UW]) is gratefully acknowledged.