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).

vel_model
The Earth-Flattening Transformation applied to the PREM Model

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)
 read(arg1, '(I2)') numl
 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
 call read_data(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, 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
subroutine read_data(delztmp, vel, velS)
 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

travel_time.png

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)
 read(arg1, '(I2)') numl
 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) 
end program comp_travel_time_vel_grad

!!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
 b=(v2-v1)/h !slope of velocity gradient

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

travel_time.png

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).

xp_curvep
X(p) curve for P-wave

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)