seismic_waves.tex

% Seismic Wave Propagation Template
% Topics: P-waves, S-waves, travel times, Earth structure, ray tracing
% Style: Technical report with observational data analysis

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}{Example}[section]
\newtheorem{remark}{Remark}[section]

\title{Seismic Wave Propagation: Earth Structure and Travel Time Analysis}
\author{Geophysical Observatory}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This technical report presents a comprehensive analysis of seismic wave propagation
through Earth's interior. We examine P-wave and S-wave velocities in different Earth
layers, compute travel times using ray theory, and analyze seismograms to determine
Earth structure. The analysis includes velocity-depth profiles based on the PREM model,
Snell's law for ray tracing, and the interpretation of seismic shadow zones that reveal
Earth's liquid outer core.
\end{abstract}

\section{Introduction}

Seismic waves generated by earthquakes provide the primary means of probing Earth's
deep interior. The velocity and attenuation of these waves depend on the elastic
properties and density of the materials through which they propagate.

\begin{definition}[Seismic Wave Types]
The two main body wave types are:
\begin{itemize}
\item \textbf{P-waves} (Primary): Compressional waves with particle motion parallel to propagation
\item \textbf{S-waves} (Secondary): Shear waves with particle motion perpendicular to propagation
\end{itemize}
S-waves cannot propagate through liquids (zero shear modulus).
\end{definition}

\section{Theoretical Framework}

\subsection{Wave Velocities}

\begin{theorem}[Seismic Velocities]
For an isotropic elastic medium with bulk modulus $K$, shear modulus $\mu$, and density $\rho$:
\begin{align}
V_P &= \sqrt{\frac{K + \frac{4}{3}\mu}{\rho}} \\
V_S &= \sqrt{\frac{\mu}{\rho}}
\end{align}
The ratio $V_P/V_S = \sqrt{(K/\mu + 4/3)} \approx 1.7$ for typical rocks.
\end{theorem}

\begin{definition}[Poisson's Ratio]
Poisson's ratio relates to seismic velocities:
\begin{equation}
\nu = \frac{V_P^2 - 2V_S^2}{2(V_P^2 - V_S^2)}
\end{equation}
For fluids, $\nu = 0.5$; for typical rocks, $\nu \approx 0.25$.
\end{definition}

\subsection{Ray Theory}

\begin{theorem}[Snell's Law for Seismic Rays]
At an interface between layers with velocities $V_1$ and $V_2$:
\begin{equation}
\frac{\sin i_1}{V_1} = \frac{\sin i_2}{V_2} = p
\end{equation}
where $p$ is the ray parameter (constant along a ray path).
\end{theorem}

\begin{remark}[Critical Angle]
When $V_2 > V_1$, a critical angle $i_c = \arcsin(V_1/V_2)$ exists. For incidence
angles greater than $i_c$, total internal reflection occurs.
\end{remark}

\subsection{Travel Time Equations}

\begin{theorem}[Travel Time for Constant Velocity Layer]
For a horizontal layer of thickness $h$ and velocity $V$:
\begin{equation}
T(x) = \frac{1}{V}\sqrt{x^2 + 4h^2}
\end{equation}
where $x$ is the horizontal distance (epicentral distance).
\end{theorem}

\section{Computational Analysis}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

np.random.seed(42)

# PREM-like Earth model (Preliminary Reference Earth Model)
def prem_velocities(depth):
    """Return Vp, Vs (km/s) for given depth (km)."""
    if depth < 15:  # Upper crust
        return 5.8, 3.2
    elif depth < 35:  # Lower crust
        return 6.8, 3.9
    elif depth < 220:  # Lithosphere/Upper mantle
        return 8.0, 4.4
    elif depth < 400:  # Upper mantle
        return 8.5, 4.6
    elif depth < 670:  # Transition zone
        return 9.9, 5.4
    elif depth < 2891:  # Lower mantle
        return 13.0 - (depth - 670) * 0.001, 7.0 - (depth - 670) * 0.0005
    elif depth < 5150:  # Outer core (liquid)
        return 10.0 + (depth - 2891) * 0.0004, 0.0
    else:  # Inner core
        return 11.0, 3.5

# Create velocity profile
depths = np.linspace(0, 6371, 1000)
Vp = np.array([prem_velocities(d)[0] for d in depths])
Vs = np.array([prem_velocities(d)[1] for d in depths])

# Density profile (simplified)
def prem_density(depth):
    if depth < 35:
        return 2.7
    elif depth < 670:
        return 3.4 + depth * 0.0005
    elif depth < 2891:
        return 4.4 + (depth - 670) * 0.001
    elif depth < 5150:
        return 10.0 + (depth - 2891) * 0.001
    else:
        return 13.0

rho = np.array([prem_density(d) for d in depths])

# Travel time calculation (simplified ray tracing)
def travel_time(distance_deg, wave_type='P'):
    """Calculate travel time for given epicentral distance."""
    # Simplified: average velocity approximation
    distance_km = distance_deg * 111.19
    if wave_type == 'P':
        # P-wave average velocity varies with distance
        if distance_deg < 100:
            v_avg = 8.5 + distance_deg * 0.02
        else:
            v_avg = 13.0  # Core phases
        return distance_km / v_avg
    else:  # S-wave
        if distance_deg < 100:
            v_avg = 4.8 + distance_deg * 0.01
            return distance_km / v_avg
        else:
            return np.nan  # S-wave shadow zone

# Calculate travel times
distances = np.linspace(0, 180, 181)
tt_P = np.array([travel_time(d, 'P') for d in distances])
tt_S = np.array([travel_time(d, 'S') for d in distances])

# Generate synthetic seismogram
def synthetic_seismogram(distance_deg, duration=600):
    t = np.linspace(0, duration, duration * 10)
    signal = np.zeros_like(t)

    # P-wave arrival
    t_p = travel_time(distance_deg, 'P')
    if not np.isnan(t_p):
        idx_p = int(t_p * 10)
        if idx_p < len(signal) - 50:
            signal[idx_p:idx_p+50] = np.sin(2*np.pi*np.arange(50)/10) * np.exp(-np.arange(50)/15)

    # S-wave arrival
    t_s = travel_time(distance_deg, 'S')
    if not np.isnan(t_s):
        idx_s = int(t_s * 10)
        if idx_s < len(signal) - 100:
            signal[idx_s:idx_s+100] = 1.5 * np.sin(2*np.pi*np.arange(100)/15) * np.exp(-np.arange(100)/25)

    # Add noise
    signal += 0.05 * np.random.randn(len(t))
    return t, signal

# Calculate Poisson's ratio
nu = np.zeros_like(depths)
for i in range(len(depths)):
    if Vs[i] > 0:
        nu[i] = (Vp[i]**2 - 2*Vs[i]**2) / (2*(Vp[i]**2 - Vs[i]**2))
    else:
        nu[i] = 0.5  # Liquid

# Vp/Vs ratio
Vp_Vs = np.where(Vs > 0, Vp / Vs, np.nan)

# Ray parameter calculation
def ray_parameter(takeoff_angle, v_surface=6.0):
    return np.sin(np.radians(takeoff_angle)) / v_surface

# Impedance contrast
impedance = rho * Vp

# Create figure
fig = plt.figure(figsize=(14, 12))

# Plot 1: Velocity profile
ax1 = fig.add_subplot(3, 3, 1)
ax1.plot(Vp, depths, 'b-', linewidth=2, label='$V_P$')
ax1.plot(Vs, depths, 'r-', linewidth=2, label='$V_S$')
ax1.axhline(y=2891, color='gray', linestyle='--', alpha=0.7, label='CMB')
ax1.axhline(y=5150, color='gray', linestyle=':', alpha=0.7, label='ICB')
ax1.set_xlabel('Velocity (km/s)')
ax1.set_ylabel('Depth (km)')
ax1.set_title('PREM Velocity Profile')
ax1.legend(fontsize=8)
ax1.invert_yaxis()

# Plot 2: Travel time curves
ax2 = fig.add_subplot(3, 3, 2)
ax2.plot(distances, tt_P, 'b-', linewidth=2, label='P-wave')
ax2.plot(distances, tt_S, 'r-', linewidth=2, label='S-wave')
ax2.axvline(x=103, color='gray', linestyle='--', alpha=0.7)
ax2.axvline(x=143, color='gray', linestyle=':', alpha=0.7)
ax2.set_xlabel('Epicentral Distance (deg)')
ax2.set_ylabel('Travel Time (s)')
ax2.set_title('Travel Time Curves')
ax2.legend(fontsize=8)
ax2.text(115, 800, 'Shadow\nZone', fontsize=8, ha='center')

# Plot 3: Synthetic seismogram
ax3 = fig.add_subplot(3, 3, 3)
t_seis, signal = synthetic_seismogram(60)
ax3.plot(t_seis, signal, 'k-', linewidth=0.5)
t_p_60 = travel_time(60, 'P')
t_s_60 = travel_time(60, 'S')
ax3.axvline(x=t_p_60, color='blue', linestyle='--', label=f'P: {t_p_60:.0f}s')
ax3.axvline(x=t_s_60, color='red', linestyle='--', label=f'S: {t_s_60:.0f}s')
ax3.set_xlabel('Time (s)')
ax3.set_ylabel('Amplitude')
ax3.set_title('Synthetic Seismogram (60 deg)')
ax3.legend(fontsize=8)
ax3.set_xlim([0, 600])

# Plot 4: Poisson's ratio
ax4 = fig.add_subplot(3, 3, 4)
ax4.plot(nu, depths, 'g-', linewidth=2)
ax4.axhline(y=2891, color='gray', linestyle='--', alpha=0.7)
ax4.axvline(x=0.5, color='red', linestyle=':', alpha=0.7, label='Liquid')
ax4.set_xlabel("Poisson's Ratio")
ax4.set_ylabel('Depth (km)')
ax4.set_title("Poisson's Ratio Profile")
ax4.invert_yaxis()
ax4.set_xlim([0.2, 0.52])
ax4.legend(fontsize=8)

# Plot 5: Vp/Vs ratio
ax5 = fig.add_subplot(3, 3, 5)
ax5.plot(Vp_Vs, depths, 'purple', linewidth=2)
ax5.axhline(y=2891, color='gray', linestyle='--', alpha=0.7)
ax5.set_xlabel('$V_P/V_S$')
ax5.set_ylabel('Depth (km)')
ax5.set_title('Velocity Ratio Profile')
ax5.invert_yaxis()
ax5.set_xlim([1.5, 2.5])

# Plot 6: Density profile
ax6 = fig.add_subplot(3, 3, 6)
ax6.plot(rho, depths, 'brown', linewidth=2)
ax6.axhline(y=2891, color='gray', linestyle='--', alpha=0.7)
ax6.set_xlabel('Density (g/cm$^3$)')
ax6.set_ylabel('Depth (km)')
ax6.set_title('Density Profile')
ax6.invert_yaxis()

# Plot 7: S-P time difference
ax7 = fig.add_subplot(3, 3, 7)
ts_tp = tt_S - tt_P
ax7.plot(distances[:100], ts_tp[:100], 'purple', linewidth=2)
ax7.set_xlabel('Epicentral Distance (deg)')
ax7.set_ylabel('S-P Time (s)')
ax7.set_title('S-P Time Difference')

# Plot 8: Record section
ax8 = fig.add_subplot(3, 3, 8)
for i, dist in enumerate(range(20, 100, 10)):
    t_rec, sig = synthetic_seismogram(dist, 500)
    ax8.plot(t_rec, sig * 0.3 + dist, 'k-', linewidth=0.5)
ax8.set_xlabel('Time (s)')
ax8.set_ylabel('Distance (deg)')
ax8.set_title('Record Section')
ax8.set_xlim([0, 500])

# Plot 9: Impedance contrast
ax9 = fig.add_subplot(3, 3, 9)
ax9.plot(impedance, depths, 'orange', linewidth=2)
ax9.axhline(y=2891, color='gray', linestyle='--', alpha=0.7)
ax9.axhline(y=35, color='blue', linestyle=':', alpha=0.7, label='Moho')
ax9.set_xlabel('Impedance (kg/m$^2$/s)')
ax9.set_ylabel('Depth (km)')
ax9.set_title('Acoustic Impedance')
ax9.invert_yaxis()
ax9.legend(fontsize=8)

plt.tight_layout()
plt.savefig('seismic_waves_analysis.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Key values
Vp_mantle = 8.0
Vs_mantle = 4.4
Vp_Vs_mantle = Vp_mantle / Vs_mantle
CMB_depth = 2891
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{seismic_waves_analysis.pdf}
\caption{Seismic wave analysis: (a) PREM velocity profile; (b) Travel time curves showing
shadow zone; (c) Synthetic seismogram at 60 degrees; (d) Poisson's ratio showing liquid
outer core; (e) $V_P/V_S$ ratio; (f) Density profile; (g) S-P time difference for
earthquake location; (h) Record section; (i) Acoustic impedance contrasts.}
\label{fig:seismic}
\end{figure}

\section{Results}

\subsection{Earth Structure}

\begin{pycode}
print(r"\begin{table}[htbp]")
print(r"\centering")
print(r"\caption{Earth Layer Properties from PREM}")
print(r"\begin{tabular}{lcccc}")
print(r"\toprule")
print(r"Layer & Depth (km) & $V_P$ (km/s) & $V_S$ (km/s) & $\rho$ (g/cm$^3$) \\")
print(r"\midrule")
print(r"Upper Crust & 0--15 & 5.8 & 3.2 & 2.7 \\")
print(r"Lower Crust & 15--35 & 6.8 & 3.9 & 2.9 \\")
print(r"Upper Mantle & 35--670 & 8.0--9.9 & 4.4--5.4 & 3.4--4.4 \\")
print(r"Lower Mantle & 670--2891 & 13.0 & 7.0 & 4.4--5.6 \\")
print(r"Outer Core & 2891--5150 & 10.0 & 0 & 10.0--12.2 \\")
print(r"Inner Core & 5150--6371 & 11.0 & 3.5 & 13.0 \\")
print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:layers}")
print(r"\end{table}")
\end{pycode}

\section{Discussion}

\begin{example}[S-Wave Shadow Zone]
The absence of S-waves between 103 and 180 degrees epicentral distance proves
that Earth's outer core is liquid:
\begin{itemize}
\item S-waves require shear strength to propagate
\item Liquids have zero shear modulus
\item P-waves are refracted through the core, creating a separate shadow zone
\end{itemize}
\end{example}

\begin{remark}[Earthquake Location]
The S-P time difference is used to estimate distance to an earthquake:
\begin{equation}
\Delta t_{S-P} = D \left(\frac{1}{V_S} - \frac{1}{V_P}\right)
\end{equation}
With three or more stations, the epicenter can be triangulated.
\end{remark}

\begin{example}[Moho Discontinuity]
The Mohorovicic discontinuity marks the crust-mantle boundary:
\begin{itemize}
\item Sharp velocity increase: $V_P$ from 6.8 to 8.0 km/s
\item Depth varies: 35 km (continents), 7 km (oceans)
\item Reflection coefficient depends on impedance contrast
\end{itemize}
\end{example}

\section{Conclusions}

This seismic wave analysis demonstrates:
\begin{enumerate}
\item Upper mantle velocities: $V_P = \py{f"{Vp_mantle}"}$ km/s, $V_S = \py{f"{Vs_mantle}"}$ km/s
\item $V_P/V_S$ ratio: $\py{f"{Vp_Vs_mantle:.2f}"}$ (typical of silicate rocks)
\item Core-mantle boundary at $\py{f"{CMB_depth}"}$ km depth
\item S-wave shadow zone beyond 103 degrees proves liquid outer core
\item Seismic tomography can map 3D velocity variations
\end{enumerate}

\section*{Further Reading}

\begin{itemize}
\item Stein, S. \& Wysession, M. \textit{An Introduction to Seismology, Earthquakes, and Earth Structure}. Blackwell, 2003.
\item Shearer, P.M. \textit{Introduction to Seismology}, 3rd ed. Cambridge, 2019.
\item Dziewonski, A.M. \& Anderson, D.L. Preliminary reference Earth model. \textit{Phys. Earth Planet. Inter.} 25, 297--356, 1981.
\end{itemize}

\end{document}