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) 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
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
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)
Institute of Earth Sciences Coders 