SeisTimes

SeisTimes is a toolkit for computing first-arrival seismic traveltimes in heterogeneous 2D and 3D anisotropic media.

Implemented are: - Ray tracing using the bending method for horizontally layered media. - Wavefield construction via the Lax-Friedrichs approximation of the static Hamilton-Jacobi (Eikonal) equations for weakly anisotropic media. - Supports Fast Sweeping numerical schemes for general heterogeneous media on regular grids. - Includes 1st-, 3rd-, and 5th-order Lax-Friedrichs schemes for flexible accuracy.

Installation

julia using Pkg Pkg.add(url="https://github.com/wtegtow/SeisTimes.jl")

Quick Start - Wavefront Construction (LxFS)

Wavefront exports 3 functionalities: - Solid2D - Solid3D - fast sweep

The solid object constructors take x,z or x,y,z coordinates and nx,nz or nx,ny,nz arrays of P-wave and S-wave velocities for 2D and 3D problems, respectivly.

Optional keyword arguments can be used to define anisotropic media:

```julia

iso2D = Solid2D(xcoords, zcoords, vp, vs) vti2D = Solid2D(xcoords, zcoords, vp, vs; eps=eps, gam=gam, del=del)

iso3D = Solid3D(xcoords, ycoords, zcoords, vp, vs) ort3D = Solid3D(xcoords, ycoords, zcoords, vp, vs; eps1=eps1, eps2=eps2, gam1=gam1, gam2=gam2, del1=del1, del2=del2, del3=del3)

```

With:

• 2D Thomsen parameters:
  eps, gam, del

• 3D Tsvankin parameters:
  eps1, eps2, gam1, gam2, del1, del2, del3

Once defined, the solid objects can be passed to the fast_sweep function to compute traveltimes. The following example shows a basic 2D use case. The 3D version works analogously.

```julia

using SeisTimes using GLMakie Makie.inline!(true)

2d example

h = 5 xcoords = 0:h:500 zcoords = 0:h:500

vp = zeros(length(xcoords), length(zcoords)) .+ 2000 vs = zeros(length(xcoords), length(zcoords)) .+ 1000 eps = zeros(length(xcoords), length(zcoords)) .+ 0.25 del = zeros(length(xcoords), length(zcoords)) .- 0.1

add some heterogeneity

vp[:, 18:34] .= 2200.0
vs[:, 18:34] .= 750.0

source location

source = [(100, 400)]

create solid objects

iso = Solid2D(xcoords, zcoords, vp, vs) vti = Solid2D(xcoords, zcoords, vp, vs; eps=eps, del=del)

algorithm parameter

wavemode = :S # :P, :S scheme = :LxFS5 # LxFS1 -> 1st order, LxFS3 -> 3rd order, LxFS5 -> 5th order Lax-Friedrich schemes verbose = false
maxiter=200 maxerrortol = 1e-5 # convergence criterium viscositybuffer = 2.5 # stabilizer. If too small, solution diverge, if too large, computations take longer

compute travel times

ttiso = fastsweep(iso, source, wavemode, scheme; maxiter=maxiter, maxerrortol=maxerrortol, viscositybuffer=viscositybuffer)

ttvti = fastsweep(vti, source, wavemode, scheme; maxiter=maxiter, maxerrortol=maxerrortol, viscositybuffer=viscositybuffer)

visualize

name = ["S", "qS"] imgs = [ttiso, ttvti] fig = Figure(size=(700,400)) for (i, img) in enumerate(imgs) ax = Axis(fig[1,i], title=name[i]) im = contourf!(ax, xcoords, zcoords, img, levels=100, colormap=:glasbeybwminc20n256) Colorbar(fig[2,i], im, vertical=false, height=5, label="sec") end save("docs/assets/img1.png", fig; pxperunit=2) display(fig)

```

Quick Start - 2-Point Ray Tracing

Wavefront exports 1 functionality: - ray_bending

The ray_bending function expects the following inputs:

• Matrix with layer properties:

  • For 2D: a [nlayer × 4] matrix with columns representing vp, vs, epsilon, and delta.
  • For 3D: a [nlayer × 9] matrix with columns representing vp, vs, epsilon1, epsilon2, gamma1, gamma2, delta1, delta2, and delta3.

• Depths of layer interfaces

Once defined, these objects can be passed to the ray_bending function to compute traveltimes.

Note: This ray tracing works only for horizontally layered media.

For reference, the following 2D example computes the same travel time grids shown above.

```julia

using SeisTimes using GLMakie Makie.inline!(true)

2d example

h = 5 xcoords = 0:h:500 zcoords = 0:h:500 nx, nz = length(xcoords), length(zcoords)

specify layer properties

vplayers = [2000, 2200, 2000] vslayers = [1000, 750, 1000] epslayers = [0.25, 0.25, 0.25] dellayers = [-0.1, -0.1, -0.1]

assemble in [nlayer x 4] matrix

Miso = hcat(vplayers, vslayers, zeros(3), zeros(3)) Mvti = hcat(vplayers, vslayers, epslayers, dellayers)

specify interface depths

interfacedepths = [zcoords[18], z_coords[34]]

wavemode = :S # :P, :S source = (100, 400)

compute travel times with 2-point ray tracing

ttiso = zeros(nx,nz) ttvti = zeros(nx,nz)

for x in 1:nx, z in 1:nz rcv = (xcoords[x], zcoords[z]) if src != rcv ttiso[x,z] = raybending(wavemode, Miso, interfacedepths, src, rcv; verbose=false).t ttvti[x,z] = raybending(wavemode, Mvti, interfacedepths, src, rcv; verbose=false).t end end

visualize

name = ["S", "qS"] imgs = [ttiso, ttvti] fig = Figure(size=(700,400)) for (i, img) in enumerate(imgs) ax = Axis(fig[1,i], title=name[i]) im = contourf!(ax, xcoords, zcoords, img, levels=100, colormap=:glasbeybwminc20n256) Colorbar(fig[2,i], im, vertical=false, height=5, label="sec") end display(fig)

```

As shown, the Lax-Friedrichs (LxFS) method selects the slowest branch of the multivalued qS wavefront, rather than the faster cuspidal branches. In contrast, ray tracing is able to capture the singularities of the S-wave.

The 3D variants work analogously. See the examples/ folder for more complex use cases.

References

• Grechka, V., Anisotropy and Microseismics: Theory and Practice. Society of Exploration Geophysicists, 2020, Chapter 6.
• Jiang, G. S., and D. Peng, Weighted ENO schemes for Hamilton-Jacobi equations, 2000.
• Kao, C. Y., S. Osher, and J. Qian, Lax-Friedrichs sweeping schemes for static Hamilton-Jacobi equations, 2004.