# Travel-time curve calculation for the Spherical Earth (PREM Model)

I have used the Preliminary Reference Earth Model of Dziewonski and Anderson (1981) to calculate the travel time curve. This model is very robust in the sense that it is designed to fit a wide variety of data including free-oscillation, surface wave dispersion observations, , travel time data for a number of body wave phases etc. To download the model, click here.

The travel time curve is the plot of first arrival time of a ray from source to receiver vs the epicentral distance (surface to surface) traveled. The sum of the time taken by the ray in each layer from the source to the receiver gives the total travel time. The ray path depends on the ray-parameter (which is constant for a laterally homogeneous, flat layers). Ray with different ray-parameter will propagate along different paths for a given earth model. For a normal velocity model (velocity increases with depth), smaller ray parameters will turn deeper in Earth and hence travel farther. The ray-parameter (p) varies along the travel time curve and is give by its slope at an instant.

I have taken two approach to calculate the travel time; first by modeling the earth using laterally homogeneous layers, and second by modeling using the velocity gradient. For details see Chapter 4, Shearer, 2009.

In the first case, I stack all the layers and the distance traveled by the ray is calculated by summing the travel time in each layer until the ray turns (which is decided by comparing slowness and ray-parameter value). When the slowness (reciprocal of velocity) of the layer is equal to the ray-parameter, the ray turns in that medium and then start traveling towards the surface. There is an inherent uncertainty associated with this method depending on the thickness and large number of the layer.

In different approach, we can parametrize the velocity model at a number of discrete points in depth and evaluate the distance and travel time by assuming an appropriate interpolation function between the model points. Here, for my calculation, I have considered the linear velocity gradient in the layer.

We have used the Earth-flattening transformation to calculate the travel-time for the spherical Earth. The curvature of the travel time curve in a spherical Earth can be simulated in a flat Earth model if a special velocity gradient is introduced in the half-space (Mu¨ller, 1971; Shearer, 2009).

Fortran Program for layered earth model:

```program comp_travel_time
!Calculates the travel time data for a layered, flat earth for given layer thicknesses and velocities
implicit none
real, parameter :: pi = 2*asin(1.0), a=6371
integer, parameter:: np=401
real(kind=16), allocatable, dimension(:) :: delztmp, delz, slow, vel, eta, slowS, velS, etaS, utop, ubot, ustop, usbot
real(kind=16) :: xx0, xx, tt0, tt, delz0, rr, pinit,pend, dp, xx0S, xxS, tt0S, ttS
integer:: i,j,numl, nlyr
real(kind=16), dimension(1:np) :: p
character*11 :: file1 !name of the model file
CHARACTER(len=32) :: arg1 !number of layers as argument
CALL get_command_argument(1, arg1)
allocate(delztmp(1:numl), delz(1:numl), slow(1:numl), vel(1:numl), eta(1:numl), &
slowS(1:numl), velS(1:numl), etaS(1:numl),utop(1:numl), ubot(1:numl), ustop(1:numl), usbot(1:numl))

file1='ps_prem.dat' !model file name
i=1
open(1,file=file1, status='old')
do while (i .le. numl)
read(1,*) delztmp(i), vel(i), velS(i) !layer thickness and vel from file
i=i+1
end do
close(1)
i=1
j=1
do while(i .le. numl) !For flat Earth, calculate the layer thickness, and average velocity
delz(j)=delztmp(i+1) - delztmp(i)
utop(j)=vel(i)
ustop(j)=velS(i)
ubot(j)=vel(i+1)
usbot(j)=velS(i+1)
vel(j)=(utop(j)+ubot(j))/2
velS(j)=(ustop(j)+usbot(j))/2
!write(*,'(i2,3f12.2)') j, delz(j), vel(i), vel(i+1)
nlyr=j
i=i+1
j=j+1
end do
nlyr=nlyr-1 !total number of layers

pinit=0.0017 !lower bound of ray parameter range
pend=0.1128 !upper bound of ray parameter range
dp=(pend-pinit)/np
!write(*,'(a,1x,a,1x,a,1x,a)') "p in s/km","X(p) in km","X(p) in degrees","Time in mins"
do j=1,np
p(j)=pinit + (j-1)*dp !ray paramter increasing by dp
xx0=0
xx0S=0
tt0=0
tt0S=0
!P-wave calculation
do i=1,nlyr
slow(i)=1/vel(i) !in s/km
if (slow(i) .gt. p(j) .and. slow(i) .le. 100000) then
call travel_time_calc(xx0,tt0,slow(i),p(j),delz(i),eta(i),xx,tt)
end if
!S- wave calculation
slowS(i)=1/velS(i)
if (slowS(i) .gt. p(j) .and. slowS(i) .le. 100000) then
call travel_time_calc(xx0S,tt0S,slowS(i),p(j),delz(i),etaS(i),xxS,ttS)
end if
end do

open(2, file='p_travel_time.dat')
write(2,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60
!write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60

open(3, file='s_travel_time.dat')
write(3,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
!write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
end do

deallocate(delztmp, delz, slow, vel, eta, slowS, velS, etaS)
end program comp_travel_time

!!SUBROUTINES

!reading the p, and s wave velocity, and depth from the model and convert it from spherical to flat earth
real, parameter :: pi = 2*asin(1.0), a=6371
real(kind=16):: delztmp, vel, velS
real(kind=16) :: rr
rr=a-delztmp !distance from the center of earth
delztmp=-a*log(rr/a) !earth flattening transformation
vel=(a/rr)*vel !earth flattening transformation for velocity
velS=(a/rr)*velS !earth flattening transformation for shear velocity
return
end

!Calculation of surface to surface distance of source and receiver and the travel time of ray
subroutine travel_time_calc(xx0,tt0,slow,p,delz,eta,xx,tt)
real(kind=16) :: eta,slow,p,xx,xx0,delz,tt,tt0
eta=sqrt(slow**2 - p**2)
xx=xx0+(2*p*delz*(1/eta))
tt=tt0+(2*slow**2*delz*(1/eta))
xx0=xx
tt0=tt
end```

To compile and run this program:

```datafile=ps_prem.dat
rm -f \$datafile *.exe
awk 'NR>2{print \$1,\$3,\$4}' PREM_MODELorg.model > \$datafile
if [ ! -f "\$datafile" ]; then #When the velocity model is not provided
cat >\$datafile<<eof
3 4 3
6 6 5
9 8 7
eof
fi

numl=`wc -l \$datafile | awk '{print \$1}'`
prog="comp_travel_time.f95"

exect=`echo \$prog | cut -d"." -f1` #execuatable

gfortran \$prog -o \${exect}.exe
if [ -f "\${exect}.exe" ]; then
./\${exect}.exe \$numl #argument is the number of layers```

For the velocity gradient model, the fortran program is below. For this case, I have made use of the LAYERXT subroutine based on Chris Chapman’s WKBJ Program.

```program comp_travel_time_vel_grad
!Calculates the travel time data for a layered, flat earth for given layer thicknesses and velocities
implicit none
real, parameter :: pi = 2*asin(1.0), a=6371
integer, parameter:: np=501
real(kind=16), allocatable, dimension(:) :: delztmp, delz, vel, velS, utop, ubot, ustop, usbot
real(kind=16) :: rr, pinit,pend, dp, delz0, xx, tt, xxS, ttS, xx0,tt0, xx0S, tt0S
integer:: i,j,numl, nlyr,irtr, irtrS
real(kind=16), dimension(1:np) :: p
character*12 :: file1 !name of the model file
CHARACTER(len=32) :: arg1 !number of layers as argument
CALL get_command_argument(1, arg1)
allocate(delztmp(1:numl), delz(1:numl), vel(1:numl), &
utop(1:numl), ubot(1:numl), ustop(1:numl), usbot(1:numl), velS(1:numl))

file1='ps_prem.dat' !model file name
i=1
open(1,file=file1, status='old')
open(5,file='flat_earth.dat')
do while (i .le. numl)
read(1,*) delztmp(i), vel(i), velS(i) !layer thickness and vel from file
rr=a-delztmp(i) !distance from the center of earth
delztmp(i)=-a*log(rr/a) !earth flattening transformation
vel(i)=(a/rr)*vel(i) !earth flattening transformation for velocity
velS(i)=(a/rr)*velS(i) !earth flattening transformation for shear velocity
write(5,'(i2,3f12.2)') i, delztmp(i), vel(i), velS(i)
i=i+1
end do
close(1)
i=1
j=1
do while(i .le. numl) !For flat Earth
delz(j)=delztmp(i+1) - delztmp(i)
utop(j)=1/vel(i)
ustop(j)=1/velS(i)
ubot(j)=1/vel(i+1)
usbot(j)=1/velS(i+1)
!write(*,'(i2,3f12.2)') j, delz(j), vel(i), vel(i+1)
nlyr=j
i=i+1
j=j+1
end do
nlyr=nlyr-1 !because of singularity at earth's center
write(*,*) "size_delz", nlyr
pinit=0.0017 !lower bound of ray parameter range
pend=0.1128 !upper bound of ray parameter range
dp=(pend-pinit)/(np-1)
do j=1,np
p(j)=pinit + (j-1)*dp !ray paramter increasing by dp
xx0=0
xx0S=0
tt0=0
tt0S=0
irtr=100
irtrS=100
do i=1,nlyr
if (irtr .ne. 2) then
call layerxt(p(j),delz(i),utop(i),ubot(i),xx,tt,irtr)
xx=xx0+xx
tt=tt0+tt
xx0=xx
tt0=tt
else
exit
end if
end do
open(2, file='p_travel_time.dat')
write(2,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), 2*xx, 2*xx*(360/(2*pi*a)), 2*tt/60
!write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60
do i=1,nlyr
if (ustop(i) .le. 100000 .and. usbot(i) .le. 100000) then
if (irtrS .ne. 2) then
call layerxt(p(j),delz(i),ustop(i),usbot(i),xxS,ttS,irtrS)
xxS=xx0S+xxS
ttS=tt0S+ttS
xx0S=xxS
tt0S=ttS
else
exit
end if
end if
end do
open(3, file='s_travel_time.dat')
write(3,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), 2*xxS, 2*xxS*(360/(2*pi*a)), 2*ttS/60
!write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
end do
write(*,'(a,1x,a,1x,a,1x,a)') "p in s/km","X(p) in km","X(p) in degrees","Time in mins"
deallocate(delztmp, delz, vel, &
utop, ubot, ustop, usbot, velS)

!!SUBROUTINES

! LAYERXT calculates dx and dt for a ray in a layer with a linear
! velocity gradient. This is a highly modified version of a
! subroutine in Chris Chapman's WKBJ program.
!
! Inputs: p = horizontal slowness
! h = layer thickness
! utop = slowness at top of layer
! ubot = slowness at bottom of layer
! Returns: dx = range offset
! dt = travel time
! irtr = return code
! = -1, zero thickness layer
! = 0, ray turned above layer
! = 1, ray passed through layer
! = 2, ray turned in layer, 1 leg counted in dx,dt
!
subroutine LAYERXT(p,h,utop,ubot,dx,dt,irtr)
implicit none
real(kind=16) :: p, h, utop, ubot,dx, dt, u1,u2,v1,v2,b,eta1,x1,tau1, &
dtau,eta2, x2,tau2
integer :: irtr
if (p.ge.utop) then !ray turned above layer
dx=0.
dt=0.
irtr=0
return
else if (h.eq.0.) then !zero thickness layer
dx=0.
dt=0.
irtr=-1
return
end if

u1=utop
u2=ubot
v1=1./u1
v2=1./u2

eta1=sqrt(u1**2-p**2)

if (b.eq.0.) then !constant velocity layer
dx=h*p/eta1
dt=h*u1**2/eta1
irtr=1
return
end if

x1=eta1/(u1*b*p)
tau1=(log((u1+eta1)/p)-eta1/u1)/b

if (p.ge.ubot) then !ray turned within layer,
dx=x1 !no contribution to integral
dtau=tau1 !from bottom point
dt=dtau+p*dx
irtr=2
return
end if

irtr=1

eta2=sqrt(u2**2-p**2)
x2=eta2/(u2*b*p)
tau2=(log((u2+eta2)/p)-eta2/u2)/b

dx=x1-x2
dtau=tau1-tau2

dt=dtau+p*dx

return
end```

In the travel-time curve, we can notice the general increase of travel time with delta. We can notice the large shadow zone around 100 degrees. This is because of the P wave transition from the solid mantle to the liquid core; which leads to sharp decrease in velocity. This is the low-velocity zone, where the rays are bent downward. So, no rays originating at the surface can turn within the LVZ. If the seismic rays originate inside the LVZ, it can never get trapped within the zone. If the seismic attenuation is lower in the LVZ layer, then the wave originated in it can propagate a long distance (Optical fiber cables work on the same principle).

In general, the epicentral distance increases as the p decreases. When the X increases with decreasing p, the branch of the travel time curve is prograde otherwise called retrograde. The retrograde branch can be observed occasionally because of the rapid velocity transition in the Earth. The point of change from prograde to retrograde is called caustics; where X value does not change with p. At this point, the energy is focussed since the rays with different p values (or take off angles) arrive at same range(X). We can also notice in the above X(p) curve that there is sharp change in X around 0.02 and 0.04. Around 0.02, the slope is positive, which could be because of the rapid velocity transition in Earth. Since, for lower p value, the ray travels deeper; this could be the inner core boundary (ICB). This is because at ICB, the P-wave velocity increases drastically.  At point 0.04, the slope is large negative; which could be the core-mantle boundary.

• Utpal Kumar (IES, Academia Sinica)