A seismic ray tracing method based on Fibonacci search

We design a fast ray tracing technique to simulate the response of seismic sources located at depth such as passive data due to fluid injection, hydraulic fracturing, volcanic tremors or earthquakes, with receivers located at the surface recording the displacement (and particle-velocity) field. The ray tracing is based on the Fibonacci search algorithm. We consider interfaces of arbitrary geometry and homogeneous layers defined by the seismic velocity, the quality factor and the mass density. Amplitude losses consider geometrical spreading, intrinsic attenuation (related to the quality factor) and transmission losses at the interfaces. The traces corresponding to time spikes are then convolved with a Ricker wavelet to obtain band-limited synthetic seismograms. We compare the results with full-wave seismograms computed with a direct algorithm based on the Fourier pseudospectral method. The ray tracing algorithm is much faster and from a practical point of view can be used as a forward modeling algorithm for the location of different sources of


Introduction
Location algorithms of different energy sources such as micro-seismicity due to fluid injection, hydraulic fracturing, volcanic tremors or earthquakes require a fast and practical forward modeling algorithm (e.g., Battaglia and Aki, 2003;Fischer et al., 2008, Kumagai et al., 2009;Shapiro, 2015].In particular Battaglia and Aki [2003] use a very simple analytical equation to compute the amplitude decay from the source location (a volcanic tremor) to the receiver.This is due to the fact that the medium is assumed homogeneous.However, in more realistic situations, the medium is heterogeneous and a numerical algorithm is required.
A ray tracing algorithm including the source radiation pattern, attenuation and transmission losses is sufficiently efficient for this purpose.
We propose a ray tracing technique based on the Fibonacci search method.The modeling intends to simulate the response of seismic sources located at depth and recorded at the surface.We model the source radiation pattern due to tensile sources, which may approximate fairly well the energy released in cases of micro-cracks due to fluid injection, fracking and magma flow inducing volcanic tremors.We focus on direct arrivals since multiple events have a lower amplitude and are difficult to pick and detect.We compare our results to fullwave synthetic seismograms computed with the Fourier pseudospectral method [Carcione et al., 1988, Carcione, 2014].

Methodology
We compute synthetic seismograms of a viscoacoustic 2D heterogeneous medium by using a ray tracing method based on a Fibonacci search algorithm, where the geological model is described in terms of vertically stacked layers, each uniform and isotropic.We consider a single source, located at depth, and a set of receivers on the surface.

Definition of the model
Each layer is isotropic, anelastic and homogeneous, and characterized by its acoustic velocity, quality factor and density.The interfaces can have any geometrical shape provided that they do not cross each other (e.g., pinchouts are not allowed) and extend from one side of the model to the other.They are defined by a small number (∼50) of equi-spaced points along the horizontal direction (x-variable).Spline interpolation between the starting points is performed at a more densely sampled set of equi-spaced points (∼5000).The larger set of points is connected by line segments.Receivers are located on the surface.

Initial value problem
For a given source point and starting ray angle, the determination of the ray path constitutes an initial value problem.Proper ray tracing requires the solution of a two-point problem, to find the ray (defined by its starting angle) going from the source location to a receiver.We find such rays by performing a Fibonacci search on the starting angle, which requires to iteratively solve the initial value problem as follows.
With the geological model defined in the (x, z)-plane, we consider a ray departing from source point S, making an angle θ with the vertical z-axis.Because each layer is isotropic and homogeneous, the ray follows a straight path until it intersects an interface.Each interface is composed of a large number of straight segments.To determine where an intersection between the ray and an interface occurs, the algorithm loops over all the segments of an interface to the left of the source point if θ < 0, and to the right of the source point if θ > 0, checking if an intersection exists.If the loop ends without positive results then there is no intersection within the model boundaries.
Starting in the layer where the source is located, the algorithm performs three tests: First, that no intersection with the underlying interfaces exists.Second, that an intersection with the overlying interface exists within the model boundaries.And third, that no total internal reflection occurs.If any of these checks fails, no ray path is computed.If all three tests are successful, the intersection point and transmission angle are saved, and the three tests are run again on the upper layer, considering the new ray starting point and direction.This procedure continues until an intersection with the surface layer is found, or until one of the tests fails.

Source-receiver ray tracing
We define f s,i (θ) as the distance between the i-th receiver and the ray's intersection point with the surface measured along the x-axis.In order to find the source-receiver ray paths, a Fibonacci search (outlined in Appendix A) is carried out on θ to minimize f s,i (θ) for each receiver.The method requires the specification of an initial angular interval [a 1 , b 1 ] for each receiver, and the desired final interval length l.It then prescribes a fixed number (q) of f s,i (θ) evaluations at sequentially determined θ values.The integer q is chosen in accordance with equation ( 12), so that b q − a q < l.All f s,i (θ) evaluations are carried out by solving the initial value problem as described in the preceding section.
The starting interval for each receiver is obtained by shooting a large number of rays with starting angles uniformly distributed in every direction and selecting the pair of angles which produce the closest arrivals on each side.While keeping the source point fxed, the computational cost of this procedure becomes relatively smaller as the number of receivers increases, because the same set of bracketing intervals is used for each receiver.
The method is only feasible if f s,i (θ) is unimodal in the starting angular interval.This is clearly not the general case.For instance, if the receiver lies in a shadow zone, the prescribed bracketing intervals could contain the shadow zone limits on both sides, which would correspond to two different minima.Additionally, even if the receiver is not in a shadow zone, it can occur that more than one ray reaches it.In these cases the algorithm will likely fail to converge to the receiver positions.
If the starting interval is adequate, the Fibonacci method reduces the starting interval length for each receiver by a factor 2F q , where F q is the q -th Fibonacci number.Let θ * i be the central value of the final interval corresponding to the i-th receiver.Rather than having accuracy in the angle parameter we are interested in accurate surface arrival positions.Thus, for each receiver we check that f s,i (θ * i ) < d with d the arrival distance tolerance.
The complete ray tracing results are the ray paths, as defined by the intersection points and the incidence/ transmission angles across each interface, for all receivers where acceptable rays are found.Receivers for which the tolerance condition is not satisfied, and those for which a premature termination condition is triggered at any point during the ray tracing procedure are completely discarded from additional computations.

Amplitude losses and travel time
Once we have determined the ray path from the source to a given receiver, we can compute the amplitude.The intrinsic (physical) attenuation along the ray is described by where r is the distance, a is the attenuation factor, f is the frequency of the signal, c is the wave velocity and Q is the quality factor.This equation holds for Q >> 1 [Carcione, 2014;eq. (2.123)].Velocity dispersion is neglected.
Energy is also lost at interfaces, where transmission rays do not reach the critical angle.We neglect intrinsic attenuation effects on the transmission coefficient, assumed as a 2nd-order effect.Let's assume that we record the vertical displacement at the surface.
The displacement transmission coeffi cient is (2) where Z = ρc is the impedance, ρ is the mass density and θ is the ray angle with respect to the line perpendicular to the interface.Primed quantities correspond to the transmission medium.Equation ( 2) is complemented with Snell's law: sin θ/c = sin θʹ/cʹ.The loss by transmission involves the energy fl ux and is given by (3) [Carcione, 2014].Equations ( 2) and ( 3) are demonstrated in Appendix B.Moreover, the signal decays by geometrical spreading.In 2D media, the decay is proportional to √r, i.e., a pulse with amplitude 1 at the source location will decay to 1/√r at a distance r from the source.
The total amplitude loss from the source to the surface is given by (4) where N is the number of interfaces and r i is the length of the ray in medium i.In addition to the amplitude loss, we compute the travel times of the ray inside each layer as r i /c i and add them over all the traversed layers to obtain the total travel time t * .A pulse at time t * , with amplitude A given by equation ( 4), represents the system impulse response if the source has amplitude 1 at time t = 0.Only fi rst arrivals are considered.

Source spectrum and radiation pattern
To obtain a more realistic seismogram, we perform a convolution between the impulse response and a wavelet representing the source time history.We consider a Ricker wavelet of the form (5) where T is the period of the wave and we take t s = 1.4T.The peak frequency is f p = 1/T.The frequency involved in equation ( 1) is f = f p .The frequency spectrum is (6) The source has a radiation pattern, i.e., it does not emit the same amplitude at every angle.Let us de-note with δ, λ and φ the dip, rake and strike angles, respectively, and ϕ is the slope angle describing the tensility of the source, such that ϕ = 90° for pure extensive sources, ϕ = 0° for pure shear sources and ϕ = − 90° for pure compressive sources.According to Carcione et al. [2015], a tensile source in 2D space is obtained for ϕ = λ = φ = 90° , while shear sources are described by ϕ = 0 and λ = φ = 90° .Then, using equation (1) of Kwiatek and Ben-Zion [2013] we get the radiation pattern for P waves due to a tensile source: (7) [see also Ou, 2008], where θ is the angle between the ray and the vertical direction and ν is the Poisson ratio of the medium.The source strength versus θ is A 0 R(θ), where A 0 is a constant.

Examples
Figure 1 shows a 3-layer subsoil model, and the ray paths traced from a source at (1.5,−1) km to a set of equi-spaced receivers at the surface, in the 0.5-2.5 km range, separated by 50 m.Interfaces where initially defi ned by 61 points with 50 m spacing, and then interpolated with 1 m spacing, defi ning 3000 line segments.The initial bracketing angles were equi-spaced with π/20 radians, and a fi nal interval length shorter than 10 −5 radians was required.From equation ( 12) it follows that the total number of evaluations per traced ray is 20.The largest ray-source arrival distance is less than 3 cm.Having completed the ray tracing from the source point to the receivers, we compute the amplitude recorded at each receiver.Figure 2 shows the results with the surface receivers now separated by 10 m.We consider the radiation pattern of a pure tensile source, as given by equation ( 7), and compute the amplitudes for different values of the source dip angle δ (ϕ = λ = φ = 90°).The Poisson ratio of the lower medium is 0.25.Figure 3 shows the seismograms for three different dip angles.The source time history is a Ricker wavelet given by equation ( 5) with f p = 30 Hz.As can be seen in Figures 2 and 3, the dip angle highly affects the amplitude distribution.Figures 4, 5 and 6 compare the results of the ray tracing with full-wave simulations (vertical component of the particle velocity) for horizontal and dipping layers.In all the three cases we consider isotropic sources, with f p = 30 Hz.For the ray tracing, the initial bracketing angles, the number of points defi ning each interface, and the fi nal angular interval length were confi gured as in the example given in Figure 1.A velocity inversion is considered in Figure 5.
For the full-wave simulations, the peak frequency of the relaxation mechanism coincided with f p .The dipping interfaces of Figures 4 and 5 require a smaller spacing between the grid points of the full-wave simu-    lation than those of the horizontal layers of Figure 6.In turn, an smaller time step is needed in the Runge-Kutta method to satisfy the stability condition [Carcione, 2014;eq. 9.12].As a consequence, the run time of the full-wave simulations of Figures 4 and 5 is nearly 10 minutes, while the full-wave modelling of Figure 6 requires 140 s.The ray tracing is in good agreement with the full-wave simulations, requiring less than 2 s.The reported times correspond to an Intel Core i5 4590 CPU.

Conclusions
We have developed a seismic ray tracing algorithm based on Fibonacci search.The method is fast enough to be applied as a forward modeling kernel to be used in inversion algorithms to locate seismic energy sources and model their seismic responses as well.The algorithm considers the effects of geometrical spreading, transmission losses, intrinsic attenuation and the source radiation pattern, and shows a good agreement with full-wave numerical simulations.sult in a smaller final uncertainty interval.
Let f(x) be the function to minimize, q the total number of iterations and [a m , b m ] the uncertainty interval after the m-th iteration.The algorithm prescribes the evaluation of the function at two interior points, c m and d m given by ( 9) . This is due to the function's unimodality.In both cases, it follows from equation ( 9) that the new interval length is (10) Additionally, it can be shown (e.g., Bazaraa et al., 2006) that for the next iteration, either c m+1 = d m or d m+1 = c m .Thus every iteration after the first requires only one function evaluation.After q − 1 iterations (and q function evaluations) the interval length is (11) For the very last evaluation, equation (9) yields c q −1 = d q −1, both in the center of interval [a q−1 , b q−1 ].In order to distinguish them one of them is slightly displaced.Taking d q −1 = c q −1 + ε, with ε very small, it is possible to reduce final interval length by approximately 1/2, obtaining (12) The total number of iterations must be selected a priori, taking equation ( 12) into account to achieve the desired accuracy level.

Appendix B: Transmission losses
To prove equations ( 2) and (3), we perform a mathematical analogy between SH and P waves.The SH wave equation is (13) where µ is the rigidity, ρ is the mass density, σ is stress, v is particle velocity, a dot above a variable denotes time derivative and ∂ i is a spatial derivative with respect to the variable x i .The indices 1, 2 and 3 correspond to the spatial variables x, y and z, respectively.

Appendix C: Full-wave modeling method
The full-wave synthetic seismograms are computed with a modeling code based on the viscoacoustic stress-strain relation corresponding to a single relaxation mechanism, based on the Zener mechanical model.The equations are given in Section 2.10.4 of Carcione [2014].The 2D particle velocity-stress formulation is (15) where v is particle velocity, σ is stress, f i are directional forces, s is the source (explosion), e is a memory variable, and τ denotes relaxation times.These are given by and ( 16) where Q 0 is the minimum quality factor and τ 0 is defined as follows.If f p is the central frequency of the source wavelet, we assume that the relaxation peak is located at ω 0 = 1/τ 0 = 2πf p .The velocity c in these equations corresponds to the unrelaxed or high-frequency limit velocity.The numerical algorithm is based on the Fourier pseudospectral method for computing the spatial derivatives and a 4th-order Runge-Kutta technique for calculating the wavefield recursively in time [e.g., Carcione, 2014].

Figure 1 .
Figure 1.Ray paths traced from the source point S to a set of surface receivers, represented by circles.The smooth horizontal curves indicate the layer interfaces, Q is the quality factor, c is the acoustic velocity and ρ is the density.

Figure 2 .
Figure 2. Surface amplitudes corresponding to the geological model of Figure 1 and a pure tensile source with amplitude 1 at S. Each curve corresponds to a different source dip angle δ.

Figure 3 .
Figure 3. Synthetic seismograms corresponding to the geological model of Figure 1.The source time history is a Ricker wavelet given by equation (5) with f p =30 Hz.The source dip angle δ is different for each example.(a) δ = 0°, (b) δ = 45° and (c) δ = 90°.

Figure 4 .
Figure 4. (a) Geological model and ray paths determined from source S at (1.498,−1.204)km to geophones on the surface indicated by circles.Actual ray tracing is performed from S to 200 geophones equispaced in the 0.5-2.5 km range.(b) Ray-tracing amplitudes compared to full-wave simulations (vertical component of the particle velocity) corresponding to the model in (a).The amplitude values are scaled with respect the maximum ones.The full-wave simulation uses a 429×221 mesh with 7 m spacing, and 0.4 ms as time step.

Figure 5 .
Figure 5. (a) Geological model and ray paths determined from source S at (2.1,−1.204)km to geophones on the surface indicated by circles.Actual ray tracing is performed from S to 200 geophones equispaced in the 0.5-2.5 km range.(b) Ray-tracing amplitudes compared to full-wave simulations (vertical component of the particle velocity) corresponding to the model in (a).The amplitude values are scaled with respect the maximum ones.The full-wave simulation uses a 429×221 mesh with 7 m spacing, and 0.3 ms as time step.