Global Seismic Phases

Global Seismic Phases#

Colab note: This notebook is designed to run on Google Colab. The first code cell installs dependencies. Open In Colab

Learning Objectives:

  • Identify major global phases on seismograms

  • Understand phase naming conventions (P, S, c, K, i)

  • Analyze travel time curves for global phases

Prerequisites: Ray theory basics

Reference: Shearer, Chapter 4 (Ray Theory)

Notebook Outline:

This notebook will explore global propagation of rays as predicted from 1D Earth models and observed in data.

Hide code cell source

 
## Install dependencies (for Google Colab)
# !pip install obspy

Hide code cell source

 
## Import necessary libraries
import numpy as np
import matplotlib.pyplot as plt
import obspy

from obspy.taup import TauPyModel
model = TauPyModel(model="iasp91")
arrivals = model.get_travel_times(source_depth_in_km=55,
                                  distance_in_degree=67)
print(arrivals)

arrivals = model.get_travel_times(source_depth_in_km=100,
                                  distance_in_degree=45,
                                  phase_list=["P", "PSPSPS"])
print(arrivals)

arrivals = model.get_ray_paths(500, 130)
arrivals = model.get_pierce_points(500, 130)

arrivals = model.get_ray_paths(
    source_depth_in_km=500, distance_in_degree=130, phase_list=["ttbasic"])
ax = arrivals.plot_rays()

arrivals = model.get_ray_paths(source_depth_in_km=500,
                               distance_in_degree=130,
                               phase_list=["Pdiff", "Sdiff",
                                           "pPdiff", "sSdiff"])

from obspy.taup import plot_travel_times
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(9, 9))
ax = plot_travel_times(source_depth=10, phase_list=["P", "S", "PP"],
                       ax=ax, fig=fig)

Student exercise#

Download data from a large earthquakes.

You can find earthquakes on the USGS website: https://www.usgs.gov/programs/earthquake-hazards

Pick a large earthquakes, get his date. Another way to collect earthquake information data is by querying an obspy earthquake catalog

from obspy.clients.fdsn import Client
client = Client("IRIS")

eq = client.get_events(starttime=obspy.UTCDateTime("2023-01-01T00:00:00"),
                       endtime=obspy.UTCDateTime("2025-01-01T00:00:00"),
                       minmagnitude=7.5)

print(eq)
# find the largest magnitude
mags = [e.magnitudes[0].mag for e in eq]
idx = np.argmax(mags)
eqmax = eq[idx]
starttime = obspy.UTCDateTime(eqmax.origins[0].time)  # Earthquake origin time
endtime = starttime + 3*3600  # End time (3 hour after origin)

get station information#

inv = client.get_stations(network="II",station="*",channel="BHZ", starttime=starttime, endtime=endtime)
print(inv)

download waveform data#

Using Obspy, download the BHZ channel of the network of stations and append the data into a single array of traces.

st=[]#obspy.Trace()
for sta in inv[0]:
  print(sta.code)
  try:
    st.append(client.get_waveforms("II", sta.code, "*", "BHZ", starttime, endtime,attach_response=True))  # Get BHZ component data
  except:
    print(f"Could not get data for station {sta.code}")

calculate angular distance#

from obspy.geodetics import gps2dist_azimuth,locations2degrees


event_latitude = eqmax.origins[0].latitude  # Get event latitude from the first trace
event_longitude = eqmax.origins[0].longitude  # Get event longitude

distances = []
ang_distances = []
for tr in st:
    for network in inv:
        if network.code == tr[0].stats.network:
            for station in network:
                if station.code == tr[0].stats.station:
                  station_latitude = station.latitude
                  station_longitude = station.longitude
    dist, az, baz = gps2dist_azimuth(event_latitude, event_longitude, station_latitude, station_longitude)
    ang = locations2degrees(event_latitude, event_longitude, station_latitude, station_longitude)
    # also get the angular distance

    ang_distances.append(ang)
    distances.append(dist / 1000)  # Convert distance to kilometers
    print(f"Distance between event and station {station.code}: {dist / 1000} km and angular distances: {ang} deg")

calculate travel times#

model = TauPyModel(model="iasp91")
travel_times = []
for dist in ang_distances:
  arrivals = model.get_travel_times(source_depth_in_km=eqmax.origins[0].depth/1000, distance_in_degree=dist, phase_list=["P"])
  print(arrivals)
  if len(arrivals):
    travel_times.append(arrivals[0].time)
  else:
    travel_times.append(np.nan)

Plot vertical seismogram#

Make a plot of the vertical seismograms as a function of

from obspy.signal.trigger import classic_sta_lta
fig, ax = plt.subplots(figsize=(10, 6))


  # STA/LTA parameters
sta_window = 5 # Short-term window in seconds
lta_window = 1000  # Long-term window in seconds

for i, tr in enumerate(st):
  time_shift = 0#travel_times[i]  # Shift the time axis by the travel time
  newtr = tr.copy()[0] # some stations may have several sensors at the same site, just pick the first one.
  # filter and resample data
  newtr.filter("bandpass", freqmin=0.001, freqmax=1)
  newtr.resample(10)
  # normalize data
  newtr.data/=np.max(np.abs(newtr.data))

  # create sta/lta filter
  # newtr.taper(0.05, type='cosine')


  # Convert to sample points
  sta_samples = int(sta_window / newtr.stats.delta)
  lta_samples = int(lta_window / newtr.stats.delta)

  # Compute STA/LTA function
  impulsivity = classic_sta_lta(newtr.data, sta_samples, lta_samples)


  ax.plot(impulsivity+ang_distances[i],newtr.times()/3600, "k")  # Plot the trace with an offset for visualization
  # ax.plot(newtr.data+ang_distances[i],newtr.times()/3600, "k")  # Plot the trace with an offset for visualization

  # plot markers for the predicted arrival times of P waves
  ax.scatter(ang_distances[i], travel_times[i]/3600, marker="x", color="r")

ax.set_xlabel("Angular Distance (°)")
ax.set_ylabel("Time (hour)")
ax.set_title("Vertical Seismograms ")
plt.show()
Contents