Wave Equations - Exercises¶

Run in Google Colab

Complete the following exercises to practice the concepts from the wave equation lecture.

In [ ]:

# Install dependencies (run this cell first on Google Colab)
!pip install obspy -q

# Install dependencies (run this cell first on Google Colab) !pip install obspy -q

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.geodetics import gps2dist_azimuth

import numpy as np import matplotlib.pyplot as plt from obspy import UTCDateTime from obspy.clients.fdsn import Client from obspy.geodetics import gps2dist_azimuth

---

Exercise 1: Calculate P and S Wave Velocities¶

Given the following rock properties, calculate the P-wave and S-wave velocities.

Recall:

• $V_P = \sqrt{\frac{\lambda + 2\mu}{\rho}}$
• $V_S = \sqrt{\frac{\mu}{\rho}}$

In [ ]:

# Rock properties for basalt
rho = 2900      # density in kg/m^3
mu = 35e9       # shear modulus in Pa
lam = 40e9      # Lamé's first parameter in Pa

# TODO: Calculate Vp and Vs (in km/s)
Vp = ???  # Vp = sqrt((lambda + 2*mu) / rho)
Vs = ???  # Vs = sqrt(mu / rho)

print(f"P-wave velocity: Vp = {Vp:.2f} km/s")
print(f"S-wave velocity: Vs = {Vs:.2f} km/s")
print(f"Vp/Vs ratio: {Vp/Vs:.2f}")

# Rock properties for basalt rho = 2900 # density in kg/m^3 mu = 35e9 # shear modulus in Pa lam = 40e9 # Lamé's first parameter in Pa # TODO: Calculate Vp and Vs (in km/s) Vp = ??? # Vp = sqrt((lambda + 2*mu) / rho) Vs = ??? # Vs = sqrt(mu / rho) print(f"P-wave velocity: Vp = {Vp:.2f} km/s") print(f"S-wave velocity: Vs = {Vs:.2f} km/s") print(f"Vp/Vs ratio: {Vp/Vs:.2f}")

Question 1a: What is the Vp/Vs ratio you calculated? How does it compare to the typical value of ~1.73?

Your answer here:

Question 1b: If a rock has $\mu = 0$ (like water), what happens to Vs? What does this mean physically?

Your answer here:

---

Exercise 2¶

Recap: A harmonic wave $u(x,t) = A\cos(kx - \omega t)$ is described by frequency $f$, period $T=1/f$, wavelength $\lambda = c/f$, angular frequency $\omega = 2\pi f$, and wavenumber $k = 2\pi/\lambda$. The wave speed satisfies $c = f\lambda = \omega/k$.

Wave Properties

A seismic wave has the following properties:

• Frequency: $f = 2$ Hz
• Wave speed: $c = 4$ km/s

Calculate the other wave properties.

Recall:

• $T = 1/f$ (period)
• $\omega = 2\pi f$ (angular frequency)
• $\lambda = c/f$ (wavelength)
• $k = 2\pi/\lambda$ (wave number)

In [ ]:

f = 2.0    # Hz
c = 4.0    # km/s

# TODO: Calculate the wave properties
T = ???  # T = 1/freq
omega = ???  # omega = 2*pi*freq
wavelength = ???  # wavelength = c/freq
k = ???  # k = 2*pi/wavelength

print(f"Period: T = {T} s")
print(f"Angular frequency: ω = {omega:.3f} rad/s")
print(f"Wavelength: λ = {wavelength} km")
print(f"Wave number: k = {k:.3f} rad/km")

f = 2.0 # Hz c = 4.0 # km/s # TODO: Calculate the wave properties T = ??? # T = 1/freq omega = ??? # omega = 2*pi*freq wavelength = ??? # wavelength = c/freq k = ??? # k = 2*pi/wavelength print(f"Period: T = {T} s") print(f"Angular frequency: ω = {omega:.3f} rad/s") print(f"Wavelength: λ = {wavelength} km") print(f"Wave number: k = {k:.3f} rad/km")

Question 2: Verify that $c = \omega/k$. Show your calculation:

In [ ]:

# TODO: Verify c = omega/k

# TODO: Verify c = omega/k

---

Exercise 3: Plot a Harmonic Wave¶

Plot the wave $u(x,t) = A \cos(kx - \omega t)$ at $t = 0$ and $t = T/4$.

Use the values you calculated in Exercise 2, with amplitude $A = 1$.

In [ ]:

A = 1.0
x = np.linspace(0, 10, 500)  # distance in km

# TODO: Calculate u(x, t) at t=0 and t=T/4
t1 = 0
t2 = ???  # T/4

u1 = ???  # A * cos(k*x - omega*0)
u2 = ???  # A * cos(k*x - omega*t2)

# Plot
plt.figure(figsize=(12, 5))
plt.plot(x, u1, 'b-', linewidth=2, label='t = 0')
plt.plot(x, u2, 'r-', linewidth=2, label='t = T/4')
plt.xlabel('Distance x (km)', fontsize=12)
plt.ylabel('Displacement u', fontsize=12)
plt.title('Harmonic Wave at Two Different Times', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

A = 1.0 x = np.linspace(0, 10, 500) # distance in km # TODO: Calculate u(x, t) at t=0 and t=T/4 t1 = 0 t2 = ??? # T/4 u1 = ??? # A * cos(k*x - omega*0) u2 = ??? # A * cos(k*x - omega*t2) # Plot plt.figure(figsize=(12, 5)) plt.plot(x, u1, 'b-', linewidth=2, label='t = 0') plt.plot(x, u2, 'r-', linewidth=2, label='t = T/4') plt.xlabel('Distance x (km)', fontsize=12) plt.ylabel('Displacement u', fontsize=12) plt.title('Harmonic Wave at Two Different Times', fontsize=14) plt.legend() plt.grid(True, alpha=0.3) plt.show()

Question 3: In which direction is the wave traveling? How can you tell from the plot?

Your answer here:

---

Exercise 4: S-P Time and Distance¶

A seismologist observes that the S-wave arrives 12 seconds after the P-wave at a station.

Using typical crustal velocities ($V_P = 6$ km/s, $V_S = 3.5$ km/s), estimate the distance to the earthquake.

In [ ]:

Vp = 6.0    # km/s
Vs = 3.5    # km/s
sp_time = 12  # seconds

# The S-P time is: t_S - t_P = d/Vs - d/Vp = d(1/Vs - 1/Vp)
# Therefore: d = sp_time / (1/Vs - 1/Vp)

# TODO: Calculate the distance
distance = ???  # distance from S-P time

print(f"Estimated distance to earthquake: {distance:.1f} km")

Vp = 6.0 # km/s Vs = 3.5 # km/s sp_time = 12 # seconds # The S-P time is: t_S - t_P = d/Vs - d/Vp = d(1/Vs - 1/Vp) # Therefore: d = sp_time / (1/Vs - 1/Vp) # TODO: Calculate the distance distance = ??? # distance from S-P time print(f"Estimated distance to earthquake: {distance:.1f} km")

Question 4: If the S-P time were 24 seconds instead, what would the distance be? Is there a simple relationship?

Your answer here:

---

Exercise 5: Explore Real Earthquake Data¶

Use ObsPy to download and analyze waveforms from a real earthquake.

Task: Select a different earthquake from the catalog and download waveforms from a nearby station.

In [ ]:

# Initialize the NCEDC client
client = Client("NCEDC")

# Search for recent Northern California earthquakes
starttime = UTCDateTime() - 365 * 24 * 3600
endtime = UTCDateTime()

catalog = client.get_events(
    starttime=starttime,
    endtime=endtime,
    minlatitude=35,
    maxlatitude=42,
    minlongitude=-125,
    maxlongitude=-118,
    minmagnitude=4.0,
    orderby="magnitude"
)

print("Available earthquakes:")
for i, event in enumerate(catalog[:10]):
    origin = event.origins[0]
    mag = event.magnitudes[0]
    print(f"  [{i}] M{mag.mag:.1f} - {origin.time.strftime('%Y-%m-%d')} - "
          f"Lat: {origin.latitude:.2f}, Lon: {origin.longitude:.2f}")

# Initialize the NCEDC client client = Client("NCEDC") # Search for recent Northern California earthquakes starttime = UTCDateTime() - 365 * 24 * 3600 endtime = UTCDateTime() catalog = client.get_events( starttime=starttime, endtime=endtime, minlatitude=35, maxlatitude=42, minlongitude=-125, maxlongitude=-118, minmagnitude=4.0, orderby="magnitude" ) print("Available earthquakes:") for i, event in enumerate(catalog[:10]): origin = event.origins[0] mag = event.magnitudes[0] print(f" [{i}] M{mag.mag:.1f} - {origin.time.strftime('%Y-%m-%d')} - " f"Lat: {origin.latitude:.2f}, Lon: {origin.longitude:.2f}")

In [ ]:

# TODO: Choose an earthquake by changing event_index
event_index = 0  # Choose an earthquake index from the list above

event = catalog[event_index]
origin = event.origins[0]
magnitude = event.magnitudes[0]

print(f"Selected: M{magnitude.mag:.1f} on {origin.time.strftime('%Y-%m-%d %H:%M:%S')}")

eq_lat = origin.latitude
eq_lon = origin.longitude
eq_time = origin.time

# TODO: Choose an earthquake by changing event_index event_index = 0 # Choose an earthquake index from the list above event = catalog[event_index] origin = event.origins[0] magnitude = event.magnitudes[0] print(f"Selected: M{magnitude.mag:.1f} on {origin.time.strftime('%Y-%m-%d %H:%M:%S')}") eq_lat = origin.latitude eq_lon = origin.longitude eq_time = origin.time

In [ ]:

# Find nearby stations
inventory = client.get_stations(
    network="BK,NC",
    channel="BHZ,HHZ",
    starttime=eq_time,
    endtime=eq_time + 300,
    latitude=eq_lat,
    longitude=eq_lon,
    maxradius=2,
    level="channel"
)

# List stations with distances
stations = []
for net in inventory:
    for sta in net:
        dist_m, _, _ = gps2dist_azimuth(eq_lat, eq_lon, sta.latitude, sta.longitude)
        stations.append((net.code, sta.code, sta[0].code, dist_m/1000))

# Remove duplicates and sort
seen = set()
unique = []
for s in sorted(stations, key=lambda x: x[3]):
    if s[1] not in seen:
        seen.add(s[1])
        unique.append(s)

print("Available stations:")
for i, (net, sta, cha, dist) in enumerate(unique[:10]):
    print(f"  [{i}] {net}.{sta:5s} - {dist:6.1f} km")

# Find nearby stations inventory = client.get_stations( network="BK,NC", channel="BHZ,HHZ", starttime=eq_time, endtime=eq_time + 300, latitude=eq_lat, longitude=eq_lon, maxradius=2, level="channel" ) # List stations with distances stations = [] for net in inventory: for sta in net: dist_m, _, _ = gps2dist_azimuth(eq_lat, eq_lon, sta.latitude, sta.longitude) stations.append((net.code, sta.code, sta[0].code, dist_m/1000)) # Remove duplicates and sort seen = set() unique = [] for s in sorted(stations, key=lambda x: x[3]): if s[1] not in seen: seen.add(s[1]) unique.append(s) print("Available stations:") for i, (net, sta, cha, dist) in enumerate(unique[:10]): print(f" [{i}] {net}.{sta:5s} - {dist:6.1f} km")

In [ ]:

# TODO: Choose a station by changing station_index
station_index = 0  # Choose a station index from the list above

net, sta, cha, dist_km = unique[station_index]

# Download waveform
st = client.get_waveforms(
    network=net,
    station=sta,
    location="*",
    channel=cha,
    starttime=eq_time - 10,
    endtime=eq_time + 120
)
st.merge(fill_value=0)

print(f"Downloaded waveform from {net}.{sta} ({dist_km:.1f} km)")

# TODO: Choose a station by changing station_index station_index = 0 # Choose a station index from the list above net, sta, cha, dist_km = unique[station_index] # Download waveform st = client.get_waveforms( network=net, station=sta, location="*", channel=cha, starttime=eq_time - 10, endtime=eq_time + 120 ) st.merge(fill_value=0) print(f"Downloaded waveform from {net}.{sta} ({dist_km:.1f} km)")

In [ ]:

# Plot the waveform with P and S arrival times
Vp = 6.0
Vs = 3.5

tr = st[0].copy()
tr.detrend('demean')
tr.filter('bandpass', freqmin=1, freqmax=10)

times = tr.times() + (tr.stats.starttime - eq_time)
t_P = dist_km / Vp
t_S = dist_km / Vs

plt.figure(figsize=(14, 5))
plt.plot(times, tr.data, 'k-', linewidth=0.7)
plt.axvline(x=t_P, color='blue', linestyle='--', linewidth=2, label=f'P ({t_P:.1f} s)')
plt.axvline(x=t_S, color='red', linestyle='--', linewidth=2, label=f'S ({t_S:.1f} s)')
plt.xlabel('Time since origin (s)', fontsize=12)
plt.ylabel('Counts', fontsize=12)
plt.title(f'M{magnitude.mag:.1f} Earthquake - Station {net}.{sta} ({dist_km:.1f} km)', fontsize=13)
plt.legend()
plt.grid(True, alpha=0.3)
plt.xlim(-10, 100)
plt.show()

print(f"\nExpected S-P time: {t_S - t_P:.1f} seconds")

# Plot the waveform with P and S arrival times Vp = 6.0 Vs = 3.5 tr = st[0].copy() tr.detrend('demean') tr.filter('bandpass', freqmin=1, freqmax=10) times = tr.times() + (tr.stats.starttime - eq_time) t_P = dist_km / Vp t_S = dist_km / Vs plt.figure(figsize=(14, 5)) plt.plot(times, tr.data, 'k-', linewidth=0.7) plt.axvline(x=t_P, color='blue', linestyle='--', linewidth=2, label=f'P ({t_P:.1f} s)') plt.axvline(x=t_S, color='red', linestyle='--', linewidth=2, label=f'S ({t_S:.1f} s)') plt.xlabel('Time since origin (s)', fontsize=12) plt.ylabel('Counts', fontsize=12) plt.title(f'M{magnitude.mag:.1f} Earthquake - Station {net}.{sta} ({dist_km:.1f} km)', fontsize=13) plt.legend() plt.grid(True, alpha=0.3) plt.xlim(-10, 100) plt.show() print(f"\nExpected S-P time: {t_S - t_P:.1f} seconds")

Question 5a: Does the P-wave arrival match the predicted time (blue dashed line)? If not, why might there be a difference?

Your answer here:

Question 5b: Try different bandpass filter values. What frequency range gives the clearest P and S arrivals?

Your answer here:

---

Exercise 6: Compare Waveforms at Different Distances¶

Download waveforms from 3 stations at different distances and compare the S-P times.

In [ ]:

# TODO: Choose 3 stations at different distances
station_indices = [0, 3, 6]  # Choose 3 station indices from the list above

fig, axes = plt.subplots(3, 1, figsize=(14, 9), sharex=True)

for ax, idx in zip(axes, station_indices):
    net, sta_name, cha, dist = unique[idx]
    
    try:
        st = client.get_waveforms(net, sta_name, "*", cha, eq_time - 10, eq_time + 120)
        st.merge(fill_value=0)
        
        tr = st[0].copy()
        tr.detrend('demean')
        tr.filter('bandpass', freqmin=1, freqmax=10)
        
        times = tr.times() + (tr.stats.starttime - eq_time)
        t_P = dist / Vp
        t_S = dist / Vs
        
        ax.plot(times, tr.data, 'k-', linewidth=0.6)
        ax.axvline(x=t_P, color='blue', linestyle='--', linewidth=2)
        ax.axvline(x=t_S, color='red', linestyle='--', linewidth=2)
        ax.set_title(f'{sta_name} - {dist:.1f} km (S-P = {t_S-t_P:.1f} s)', fontsize=11)
        ax.set_ylabel('Counts')
        ax.grid(True, alpha=0.3)
    except:
        ax.text(0.5, 0.5, f'{sta_name}: No data', transform=ax.transAxes, ha='center')

axes[-1].set_xlabel('Time since origin (s)', fontsize=12)
plt.xlim(-10, 100)
plt.tight_layout()
plt.show()

# TODO: Choose 3 stations at different distances station_indices = [0, 3, 6] # Choose 3 station indices from the list above fig, axes = plt.subplots(3, 1, figsize=(14, 9), sharex=True) for ax, idx in zip(axes, station_indices): net, sta_name, cha, dist = unique[idx] try: st = client.get_waveforms(net, sta_name, "*", cha, eq_time - 10, eq_time + 120) st.merge(fill_value=0) tr = st[0].copy() tr.detrend('demean') tr.filter('bandpass', freqmin=1, freqmax=10) times = tr.times() + (tr.stats.starttime - eq_time) t_P = dist / Vp t_S = dist / Vs ax.plot(times, tr.data, 'k-', linewidth=0.6) ax.axvline(x=t_P, color='blue', linestyle='--', linewidth=2) ax.axvline(x=t_S, color='red', linestyle='--', linewidth=2) ax.set_title(f'{sta_name} - {dist:.1f} km (S-P = {t_S-t_P:.1f} s)', fontsize=11) ax.set_ylabel('Counts') ax.grid(True, alpha=0.3) except: ax.text(0.5, 0.5, f'{sta_name}: No data', transform=ax.transAxes, ha='center') axes[-1].set_xlabel('Time since origin (s)', fontsize=12) plt.xlim(-10, 100) plt.tight_layout() plt.show()

Question 6: How does the S-P time change as distance increases? Is the relationship linear?

Your answer here: