%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
Background¶
Videos:¶
http://www.iris.edu/hq/inclass/animation/4station_seismograph_network_records_a_single_earthquake
http://www.iris.edu/hq/inclass/animation/seismic_wave_motions
http://www.iris.edu/hq/inclass/video/travel_time_curves_described
http://www.iris.edu/hq/inclass/video/types_of_seismic_wave_paths_through_the_earth
Specifically, we are going to explore the seismic waves produced by the August 24, 2016 earthquake in Italy.
http://www.iris.edu/hq/files/programs/education_and_outreach/retm/tm_160824_italy/160824italy.pdf
Reference Earth Model¶
There are several reconstructions of the interior of the Earth. These are called Earth Models. http://ds.iris.edu/spud/earthmodel
To first order, the Earth is radially symmetric and so we can use a 1D Reference Earth Model. These models give the density and velocites of P- and S-Waves as a function of radius.
In this lab, we will use the IASP91 model (Kennett B.L.N. and Engdahl E.R., 1991. "Travel times for global earthquake location and phase association." Geophysical Journal International, 105:429-465.)
Our Earth Reference Model:
model_z = np.array([-2, -1, 0,19,20,34,35,77.5,120,165,209,210,260,310,360,409,410,460,510,560,610,659,660,710,760,809.5,859,908.5,958,1007.5,1057,1106.5,1156,1205.5,1255,1304.5,1354,1403.5,1453,1502.5,1552,1601.5,1651,1700.5,1750,1799.5,1849,1898.5,1948,1997.5,2047,2096.5,2146,2195.5,2245,2294.5,2344,2393.5,2443,2492.5,2542,2591.5,2641,2690.5,2740,2740,2789.67,2839.33,2888,2889,2939.33,2989.66,3039.99,3090.32,3140.66,3190.99,3241.32,3291.65,3341.98,3392.31,3442.64,3492.97,3543.3,3593.64,3643.97,3694.3,3744.63,3794.96,3845.29,3895.62,3945.95,3996.28,4046.62,4096.95,4147.28,4197.61,4247.94,4298.27,4348.6,4398.93,4449.26,4499.6,4549.93,4600.26,4650.59,4700.92,4751.25,4801.58,4851.91,4902.24,4952.58,5002.91,5053.24,5103.57,5153.9,5153.9,5204.61,5255.32,5306.04,5356.75,5407.46,5458.17,5508.89,5559.6,5610.31,5661.02,5711.74,5762.45,5813.16,5863.87,5914.59,5965.3,6016.01,6066.72,6117.44,6168.15,6218.86,6269.57,6320.29,6371])
model_Vp = np.array([0, 0, 5.8,5.8,6.5,6.5,8.04,8.045,8.05,8.175,8.3,8.3,8.4825,8.665,8.8475,9.03,9.36,9.528,9.696,9.864,10.032,10.2,10.79,10.9229,11.0558,11.144,11.23,11.314,11.396,11.4761,11.5543,11.6308,11.7056,11.7787,11.8504,11.9205,11.9893,12.0568,12.1231,12.1881,12.2521,12.3151,12.3772,12.4383,12.4987,12.5584,12.6174,12.6759,12.7339,12.7915,12.8487,12.9057,12.9625,13.0192,13.0758,13.1325,13.1892,13.2462,13.3034,13.361,13.419,13.4774,13.5364,13.5961,13.6564,13.6564,13.6679,13.6793,13.6908,8.0088,8.0963,8.1821,8.2662,8.3486,8.4293,8.5083,8.5856,8.6611,8.735,8.8072,8.8776,8.9464,9.0134,9.0787,9.1424,9.2043,9.2645,9.323,9.3798,9.4349,9.4883,9.54,9.59,9.6383,9.6848,9.7297,9.7728,9.8143,9.854,9.892,9.9284,9.963,9.9959,10.0271,10.0566,10.0844,10.1105,10.1349,10.1576,10.1785,10.1978,10.2154,10.2312,10.2454,10.2578,11.0914,11.1036,11.1153,11.1265,11.1371,11.1472,11.1568,11.1659,11.1745,11.1825,11.1901,11.1971,11.2036,11.2095,11.215,11.2199,11.2243,11.2282,11.2316,11.2345,11.2368,11.2386,11.2399,11.2407,11.2409])
model_Vs = np.array([0, 0, 3.36,3.36,3.75,3.75,4.47,4.485,4.5,4.509,4.518,4.522,4.609,4.696,4.783,4.87,5.07,5.176,5.282,5.388,5.494,5.6,5.95,6.0797,6.2095,6.2474,6.2841,6.3199,6.3546,6.3883,6.4211,6.453,6.4841,6.5143,6.5438,6.5725,6.6006,6.628,6.6547,6.6809,6.7066,6.7317,6.7564,6.7807,6.8046,6.8282,6.8514,6.8745,6.8972,6.9199,6.9423,6.9647,6.987,7.0093,7.0316,7.054,7.0765,7.0991,7.1218,7.1449,7.1681,7.1917,7.2156,7.2398,7.2645,7.2645,7.2768,7.2892,7.3015,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.4385,3.4488,3.4587,3.4681,3.477,3.4856,3.4937,3.5013,3.5085,3.5153,3.5217,3.5276,3.533,3.5381,3.5427,3.5468,3.5505,3.5538,3.5567,3.5591,3.561,3.5626,3.5637,3.5643,3.5645])
model_ρ = np.array([0, 0, 2.72,2.72,2.92,2.92,3.3198,3.3455,3.3713,3.3985,3.4258,3.4258,3.4561,3.4864,3.5167,3.547,3.7557,3.8175,3.8793,3.941,4.0028,4.0646,4.3714,4.401,4.4305,4.4596,4.4885,4.5173,4.5459,4.5744,4.6028,4.631,4.6591,4.687,4.7148,4.7424,4.7699,4.7973,4.8245,4.8515,4.8785,4.9052,4.9319,4.9584,4.9847,5.0109,5.037,5.0629,5.0887,5.1143,5.1398,5.1652,5.1904,5.2154,5.2403,5.2651,5.2898,5.3142,5.3386,5.3628,5.3869,5.4108,5.4345,5.4582,5.4817,5.4817,5.5051,5.5284,5.5515,9.9145,9.9942,10.0722,10.1485,10.2233,10.2964,10.3679,10.4378,10.5062,10.5731,10.6385,10.7023,10.7647,10.8257,10.8852,10.9434,11.0001,11.0555,11.1095,11.1623,11.2137,11.2639,11.3127,11.3604,11.4069,11.4521,11.4962,11.5391,11.5809,11.6216,11.6612,11.6998,11.7373,11.7737,11.8092,11.8437,11.8772,11.9098,11.9414,11.9722,12.0001,12.0311,12.0593,12.0867,12.1133,12.1391,12.7037,12.7289,12.753,12.776,12.798,12.8188,12.8387,12.8574,12.8751,12.8917,12.9072,12.9217,12.9351,12.9474,12.9586,12.9688,12.9779,12.9859,12.9929,12.9988,13.0036,13.0074,13.01,13.0117,13.0122])
This model is given use the velocity (in km/s) for P-waves and S-waves at different depths (in km).
Re = 6372 # km
Notices that this Reference Earth Model gives the values in terms of depth (called $z$) (in km) from the Earth's surface. We will need the radius (called $r$) from the centre of the Earth.
$$ r = R_E - z$$and
$$ z = R_E - r $$To use this Referece Earth Model, we can use the scipy module to calculate interpolations functions so that given a radius, $r$, the function Vp(r) and Vs(r) will give use the velocity of the P-wave and S-wave, respectively.
from scipy import interpolate
Vp = interpolate.interp1d(Re - model_z, model_Vp, fill_value='extrapolate')
Vs = interpolate.interp1d(Re - model_z, model_Vs, fill_value='extrapolate')
ρ = interpolate.interp1d(Re - model_z, model_ρ, fill_value='extrapolate')
To understand how this interpolation works, consider the following plot of velocity vs radius.
fig, axes = plt.subplots()
r = np.arange(0, Re, 10)
plt.plot(Vp(r), r, label = 'Vp')
plt.plot(Vs(r), r, label = 'Vs')
plt.xlabel('velocity (km/s)')
plt.ylabel('radius (km)')
plt.legend(loc='upper left')
plt.xlim(-1, 15)
plt.title('Seismic Velocities Depending on Depth')
For S-waves, the velocity is zero between $r = 1200$ and $r = 3400$. Why?
Polar plots¶
While these 1D plots are a good way of representing the velocity structure of the Earth, it can also be helpful to draw it a polar plot.
def plot_background():
# set up a polar projection plot
fig, ax = plt.subplots(figsize=(8,8),
subplot_kw=dict(projection='polar'))
# define a grid for all angles then (60 divisions around the circle)
ntheta = 60
θ = np.radians(np.linspace(0, 360, ntheta))
r = np.arange(0, Re, 10)
# define a 2D grid in radius and theta
r_plt, theta_plt = np.meshgrid(r, θ)
# the density needs to the same for all values of theta
ρ_plt = ρ(r)*np.ones((ntheta, 1))
# choose a colour map
plt.set_cmap('viridis_r')
# make a contour plot of 200 levels
cax = ax.contourf(theta_plt, r_plt, ρ_plt, 200, vmax=20)
# add a colorbar
cb = fig.colorbar(cax)
cb.set_label("Density (10$^3$ kg/m$^3$)")
return ax
plot_background()
plt.show()
Ray tracing equations¶
These the ray tracing equations introduced in Lecture 15. We need to make one change in notation: let now $\phi$ will represent the direction of the ray (as an angle to the horizontal). We do this so that we can plot our results in polar coordinates: $r$ and $\theta$.
For simplicity, let's define $x$ and $y$ to have an origin at the centre of the Earth. Then the coordinate transforms needed are
\begin{align} r &= \sqrt{x^2 + y^2} \\ \theta &= \arctan(y/x) \end{align}and
\begin{align} x &= r\cos(\theta) \\ y &= r\sin(\theta) \\ \end{align}# to make the code easier to read
sin = np.sin
cos = np.cos
arctan2 = np.arctan2
π = np.pi
# Centered difference scheme
def NDc(f, x0, dx):
return (f(x0+dx) - f(x0-dx)) / (2*dx)
# define the right hand side, F(q)
def F( q, t, velocity_model):
"""
velocity_model is the velocity from the Reference Earth Model to use
"""
# separate out the variables
x, y, ϕ = q
# need a function which gives v from x, y
# we need to first get the radius
v = lambda x, y: velocity_model(np.sqrt(x**2 + y**2))
# evaluate the right hand side of the ODE
dxdt = v(x,y)*cos(ϕ)
dydt = v(x,y)*sin(ϕ)
dvdx = NDc(lambda x: v(x,y), x, 1)
dvdy = NDc(lambda x: v(x,y), y, 1)
dϕdt = sin(ϕ)*dvdx - cos(ϕ)*dvdy
dqdt = np.array([dxdt, dydt, dϕdt])
return dqdt
Solve our ray equations.
# define initial values (x0, y0, ϕ0)
x0 = 0
y0 = Re
# ray travel into the ground at an angle to the horizontal
ϕ0 = -π/4
q0 = [x0, y0, ϕ0]
tmax = 800
dt = 10
# create an array for the time variables
t = np.arange(0, tmax, dt)
sol = integrate.odeint(F, q0, t, args = (Vp,))
x, y, ϕ = sol.T
fig, axes = plt.subplots()
plt.plot(x, y)
plt.xlabel('x (km)')
plt.ylabel('y (km)')
plt.title('Ray Trace')
It helps to understand what this really means by looking at the time series; plot each variable against time.
fig, axes = plt.subplots()
r = np.sqrt(x**2 + y**2)
plt.plot(t, x, label='x')
plt.plot(t, y, label='y')
plt.plot(t, r, label='r')
plt.ylabel('x, y, r(km)')
plt.xlabel('t (s)')
plt.legend(loc='lower right')
plt.title('Ray Trace')
This particular ray, heading a direction 45$^\circ$ to the vertical into the ground, reaches the Earth surface again just before $t=600$ s (about 10 minutes).
It also good to visualize the rays not in only Cartesian coordinates ($x$ and $y$) but also in polar coordinates ($\theta$ and $r$). Let's write a helper that computes the ray path and then converts the results to polar coordinates.
def compute_ray(dϕ = 0.1, tmax = 1200, dt = 1, velocity_model=Vp):
"""
Propagate a ray starting at θ = 90 deg and the surface r=Re,
into the ground at an angle dϕ to the vertical.
returns the ray path in polar coordinates.
"""
x0 = 0
y0 = Re
ϕ0 = -π/2 + dϕ # always send waves into the ground
q0 = [x0, y0, ϕ0]
# create an array for the time variable
t = np.arange(0, tmax, dt)
sol = integrate.odeint(F, q0, t, args=(velocity_model,))
x, y, ϕ = sol.T
r = np.sqrt(x**2 + y**2)
θ = arctan2(y, x)
return θ, r, t
A note on arctan2¶
When converting from Cartesian to polar coordinates, special attention has to be given to which quadrant we require. As well, if $x = 0$ and $y$ is non-zero, we want $\theta$ to take on a value of $\theta = -\pi/2$ even though $y/x$ would technical be undefined.
There is a special variation of the arctan function, called arctan2, that has exactly this behaviour. Notice this is not the square of arctan but just a different version of arctan. See https://en.wikipedia.org/wiki/Atan2 for more information.
Example¶
Send a P-wave into the Earth.
θ, r, t = compute_ray(dϕ = π/4, tmax=800, dt=10)
We can superimpose the ray path on our background.
plot_background()
plt.plot(θ, r, color='r')
plt.show()
The other way to visualize this wave is by its time series.
def plot_timeseries(θ, r, t, axes = None):
if axes is None:
fig, axes = plt.subplots(2,1, figsize=(8,6))
plt.sca(axes[0])
plt.plot(t, r)
plt.xlabel('t (s)')
plt.ylabel('r (km)')
plt.sca(axes[1])
plt.plot(t, np.degrees(θ))
plt.xlabel('t (s)')
plt.ylabel(r'$\theta$')
plot_timeseries(θ, r-Re, t)
Travel time¶
From the time series plot we can read off that the ray took about 560 s to travel from $\theta=90^\circ$ to $\theta=35^\circ$ or about $55^\circ$. Rather than just reading this off the graph, let's try and write a function to find these values automatically.
def travel_time(θ, r, t):
"""
Determine the travel time and distance (in degrees) for a ray to again reach the surface
"""
# start search beyond the initial condition
for i in range(1,len(r)):
if r[i] >= Re:
# if we beyond the radius of the Earth, stop.
break
if r[i] < Re:
print ('Warning: ray did not reach the surface of the Earth')
time = t[i]
distance = np.degrees(θ[0] - θ[i])
return time, distance
travel_time(θ, r, t)
Notice that the results are consistent with our qualitative estimates.
Travel Time Curves¶
Now we can put all the pieces together to create a travel time curve.
Start by generating rays for many different initial angles. Since it takes some time to calcuate each ray path, an information message is printed at the start of every calculation.
plot_background()
for dϕ in np.arange(0.1, 1.0, 0.2):
print("Calculating for dϕ = ", dϕ)
θ, r, t = compute_ray(dϕ=dϕ, dt=10)
plt.plot(θ, r, color='r')
Notice that it appears waves never get seem to make it to between $θ = 315^\circ$ to $350^\circ$. This is the P-Wave shadow zone. Compare with this figure here (from Wikipedia - Shadow zone):
To compute the travel times we need consider the time series
fig, axes = plt.subplots(2,1, figsize=(8,6))
for dϕ in np.arange(0.1, 1.0, 0.2):
print("Calculating for dϕ = ", dϕ)
θ, r, t = compute_ray(dϕ=dϕ, tmax=1200)
plot_timeseries(θ, r, t, axes=axes)
It is important that each one of the above rays actually reaches the Earth's surface for the travel_time routine to work. If the ray stops short of the reaching the Earth's surface, you need to use a larger value of tmax.
times = []
distances = []
for dϕ in np.arange(0.1, 1.0, 0.2):
print("Calculating for dϕ = ", dϕ)
θ, r, t = compute_ray(dϕ=dϕ, tmax=1300)
time, distance = travel_time(θ, r, t)
times.append(time)
distances.append(distance)
times = np.array(times)
distances = np.array(distances)
The final step would be to plot a travel time curve.
fig, axes = plt.subplots()
# plot travel time for Vp
plt.plot(times / 60, distances, 'o-', label='Vp')
plt.ylabel('Distance (degrees)')
plt.xlabel('Travel Time (minutes)')
plt.title('Travel Time Curve')
plt.legend(loc='upper left')
Application to Earthquake Detection¶
In the slides on the August 24, 2016 earthquake in Italy, there was a seismograph showing "the record of the earthquake in Bend, Oregon (BNOR) is illustrated below. Bend is 9365 km (5819 miles, 84.4°) from the location of this earthquake."
http://www.iris.edu/hq/files/programs/education_and_outreach/retm/tm_160824_italy/160824italy.pdf
Task 1: P Wave Travel Time Curve¶
Improve the resolution of the Travel Time Curve for the P wave. Specifically, can you use this curve to justify the statement from the slides on the August 24, 2016 earthquake in Italy that said: "Following the earthquake, it took 12 minutes and 31 seconds for the compressional P waves to travel a curved path through the mantle to Bend". Adjust the line of code that says
for dϕ in np.arange(0.1, 1.0, 0.2):
so that you are able to get a better estimate of the travel time from Italy to Oregon USING THE REFERENCE EARTH MODEL IN THIS LAB.
plot_background()
times = []
distances = []
for dϕ in np.arange(0.1, 1.0, 0.2):
print("Calculating for dϕ = ", dϕ)
θ, r, t = compute_ray(dϕ=dϕ, tmax = 1300)
plt.plot(θ, r, color='r')
time, distance = travel_time(θ, r, t)
times.append(time)
distances.append(distance)
times = np.array(times)
distances = np.array(distances)
fig, axes = plt.subplots()
plt.plot(times / 60, distances, 'o', label='Vp')
plt.axhline(84.4, color='k')
plt.ylabel('Distance (degrees)')
plt.xlabel('Travel Time (minutes)')
plt.title('Travel Time Curve')
plt.legend(loc='upper left')
plt.show()
Our model predicts that a P wave would travel the distance to Bend, Oregon (84.4$^\circ$) in 12.66 minutes (about 12 minutes and 40 seconds). This is about 9 seconds longer than the actual observation.
Answer: (0.56, 0.61, 0.01), tmax=14*60
Task 2: S Wave Travel Time Curve¶
Construct a travel time curve for an S-wave. Notice that the routine compute_ray is already set up to use a velocity model for S waves as well as for P waves.
The important line of code that needs to change is
θ, r, t = compute_ray(dϕ=dϕ, tmax=1300)
to
θ, r, t = compute_ray(dϕ=dϕ, tmax=1300, velocity_model=Vs)
You should then be able compare your result with the statement "S waves are shear waves that follow the same path through the mantle as P waves. S waves took 22 minutes and 55 seconds to travel from the earthquake to Bend."
plot_background()
times = []
distances = []
for dϕ in np.arange(0.1, 1.0, 0.2):
print("Calculating for dϕ = ", dϕ)
θ, r, t = compute_ray(dϕ=dϕ, tmax = 1300, velocity_model=Vp)
plt.plot(θ, r, color='r')
time, distance = travel_time(θ, r, t)
times.append(time)
distances.append(distance)
times = np.array(times)
distances = np.array(distances)
fig, axes = plt.subplots()
plt.plot(times / 60, distances, 'o', label='Vs')
plt.axhline(84.4, color='k')
plt.ylabel('Distance (degrees)')
plt.xlabel('Travel Time (minutes)')
plt.title('Travel Time Curve')
plt.legend(loc='upper left')
Our model predicts about 23 minutes for an S Wave to travel to Bend, Oregon from Italy which is consistent with the observation from the seismograph.
Answer: (0.50, 0.64, 0.03), tmax=27*60, Vs