Lab 8: Earthquake Nucleation and Dynamic Rupture

Open In Colab

Note: If running on Colab, uncomment and run the pip install cell below.

Lab 8: Earthquake Nucleation and Dynamic Rupture#

Learning Objectives#

  • Simulate the spring–slider instability and identify the transition between stable and unstable slip.

  • Visualize slip-weakening friction and compute fracture energy \(G\).

  • Explore how friction parameters (\(\tau_s\), \(\tau_d\), \(D_c\)) control the critical nucleation length \(L_c\).

  • Use an energy-balance argument to distinguish rupture propagation from arrest.

Prerequisites#

  • Lecture: Earthquake Nucleation and Dynamic Rupture (Module 8)

  • Lecture: Earthquake Source Dynamics (Module 7d)

Estimated Completion Time#

~60 minutes

# Uncomment if running on Colab:
# !pip install numpy matplotlib scipy

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# ── figure defaults ──────────────────────────────────────────────────────────────────
plt.rcParams.update({
    "figure.dpi": 110,
    "axes.grid": True,
    "grid.alpha": 0.3,
    "font.size": 11,
    "axes.titlesize": 12,
})

1. Spring–Slider Instability#

The simplest model of fault instability is a block on a frictional surface, loaded through a spring by a moving plate.

The driving stress is:

\[\tau = k\,(v_{pl}\,t - \delta)\]

where \(k\) is the spring stiffness, \(v_{pl}\) is the plate velocity, and \(\delta\) is slip.

We use a velocity-weakening friction law:

\[\tau_f(v) = \tau_d + (\tau_s - \tau_d)\,\frac{v_0}{v + v_0}\]

where \(v\) is the slip velocity, \(v_0\) is a reference velocity, and \(\tau_s > \tau_d\) provides weakening at high velocity.

Predict#

Before running the code: If you decrease the spring stiffness \(k\) while keeping friction fixed, do you expect the system to become more or less unstable? Why?

def spring_slider_rhs(t, y, k, m, v_pl, tau_s, tau_d, v0):
    """
    Right-hand side for the spring–slider ODE system.

    State vector y = [delta, v] where delta is slip and v is slip velocity.
    """
    delta, v = y
    # Driving stress from the spring
    tau_spring = k * (v_pl * t - delta)
    # Velocity-weakening friction
    tau_friction = tau_d + (tau_s - tau_d) * v0 / (abs(v) + v0)
    # Newton's second law: m * dv/dt = tau_spring - tau_friction
    dvdt = (tau_spring - tau_friction) / m
    ddelta_dt = v
    return [ddelta_dt, dvdt]


# ── Physical parameters ─────────────────────────────────────────────────────────
m = 1.0          # block mass (normalized)
v_pl = 1e-3      # plate velocity
tau_s = 1.0      # static friction strength
tau_d = 0.4      # dynamic friction strength
v0 = 1e-3        # reference velocity for friction law
t_end = 5000     # simulation duration

# Initial conditions: block at rest, no initial slip
y0 = [0.0, 0.0]

# Compare two stiffness values
k_values = [0.3, 1.5]
labels = [f"k = {k} (weak spring)" for k in k_values]

fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

for k, label in zip(k_values, labels):
    sol = solve_ivp(
        spring_slider_rhs, [0, t_end], y0,
        args=(k, m, v_pl, tau_s, tau_d, v0),
        max_step=1.0, rtol=1e-8, atol=1e-10,
    )
    axes[0].plot(sol.t, sol.y[0], label=label)
    axes[1].plot(sol.t, sol.y[1], label=label)

axes[0].set_ylabel(r"Slip $\delta$")
axes[0].legend()
axes[0].set_title("Spring–slider simulation: slip vs. time")

axes[1].set_ylabel(r"Slip velocity $v$")
axes[1].set_xlabel("Time")
axes[1].legend()
axes[1].set_title("Slip velocity vs. time")

plt.tight_layout()
plt.show()

Explain#

Interpretation: Compare the two panels.

  • Which stiffness produces stick–slip behavior (sudden jumps in slip)?

  • Which produces smooth, stable sliding?

  • How does this relate to the instability criterion \(|d\tau_f/d\delta| > k\)?

Do: Find the critical stiffness#

Sweep through a range of \(k\) values and classify each run as “stick–slip” (unstable) or “stable” based on whether the peak slip velocity exceeds a threshold.

k_range = np.linspace(0.1, 2.0, 30)
max_velocities = []

for k in k_range:
    sol = solve_ivp(
        spring_slider_rhs, [0, t_end], y0,
        args=(k, m, v_pl, tau_s, tau_d, v0),
        max_step=2.0, rtol=1e-7, atol=1e-9,
    )
    max_velocities.append(np.max(np.abs(sol.y[1])))

max_velocities = np.array(max_velocities)

fig, ax = plt.subplots(figsize=(8, 4))
ax.semilogy(k_range, max_velocities, "ko-", ms=4)
ax.axhline(10 * v_pl, color="C3", ls="--", label=r"Threshold ($10 \times v_{pl}$)")
ax.set_xlabel("Spring stiffness $k$")
ax.set_ylabel("Peak slip velocity")
ax.set_title("Identifying the critical stiffness")
ax.legend()
plt.tight_layout()
plt.show()

# Estimate critical k
threshold = 10 * v_pl
k_crit_idx = np.where(max_velocities > threshold)[0]
if len(k_crit_idx) > 0:
    k_crit_approx = k_range[k_crit_idx[-1]]
    print(f"Approximate critical stiffness: k_c \u2248 {k_crit_approx:.2f}")
    print(f"Friction weakening rate (\u0394\u03c4/D_c analogy): "
          f"({tau_s} - {tau_d}) = {tau_s - tau_d:.1f}")

Check#

Does your estimated \(k_c\) match the theoretical instability criterion? In the slip-weakening framework, the weakening rate is \((\tau_s - \tau_d)/D_c\). With the velocity-weakening law used here, the effective weakening slope scales as \((\tau_s - \tau_d)/v_0\).

2. Slip-Weakening Friction Visualization#

Now let’s look at the slip-weakening friction law directly and compute the fracture energy \(G\).

Predict#

Before running: If you double \(D_c\) while keeping \(\tau_s\) and \(\tau_d\) fixed, how does \(G\) change? What about nucleation difficulty?

def slip_weakening_friction(delta, tau_s, tau_d, Dc):
    """Linear slip-weakening friction: tau(delta)."""
    return np.where(delta < Dc, tau_s - (tau_s - tau_d) * delta / Dc, tau_d)


def fracture_energy(tau_s, tau_d, Dc):
    """Fracture energy for linear slip-weakening: G = 0.5 * (tau_s - tau_d) * Dc."""
    return 0.5 * (tau_s - tau_d) * Dc


# Parameters to compare
Dc_values = [0.5, 1.0, 2.0]
tau_s_val = 1.0
tau_d_val = 0.4
delta = np.linspace(0, 4, 500)

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

# Left: friction curves
for Dc in Dc_values:
    tau = slip_weakening_friction(delta, tau_s_val, tau_d_val, Dc)
    G = fracture_energy(tau_s_val, tau_d_val, Dc)
    axes[0].plot(delta, tau, lw=2, label=rf"$D_c = {Dc}$, $G = {G:.2f}$")

axes[0].set_xlabel(r"Slip $\delta$")
axes[0].set_ylabel(r"Shear stress $\tau$")
axes[0].set_title("Slip-weakening friction for different $D_c$")
axes[0].legend(fontsize=9)
axes[0].set_xlim(0, 4)

# Right: fracture energy as a function of Dc
Dc_sweep = np.linspace(0.01, 3, 100)
G_sweep = fracture_energy(tau_s_val, tau_d_val, Dc_sweep)
axes[1].plot(Dc_sweep, G_sweep, "C1-", lw=2.5)
axes[1].set_xlabel(r"Critical slip distance $D_c$")
axes[1].set_ylabel(r"Fracture energy $G$")
axes[1].set_title(r"$G = \frac{1}{2}(\tau_s - \tau_d)\,D_c$")

# Mark the three Dc values
for Dc in Dc_values:
    G = fracture_energy(tau_s_val, tau_d_val, Dc)
    axes[1].plot(Dc, G, "ko", ms=7)
    axes[1].annotate(rf"$D_c={Dc}$", (Dc, G), textcoords="offset points",
                     xytext=(8, 5), fontsize=9)

plt.tight_layout()
plt.show()

Explain#

Interpretation:

  • How does \(G\) scale with \(D_c\)? Is this linear, quadratic, or something else?

  • Connecting back to Module 7d: the fracture energy \(G\) you computed here is the same quantity that enters the energy budget \(\Delta E_\text{strain} = E_R + E_G + E_H\). What does a large \(G\) imply for radiation efficiency \(\eta_R\)?

3. Critical Nucleation Length#

For a 2-D fault with linear slip-weakening friction, the critical nucleation half-length scales as:

\[L_c = \frac{C\,\mu\,D_c}{(\tau_s - \tau_d)}\]

where \(\mu\) is the shear modulus and \(C\) is a geometric constant of order unity (depends on crack mode and geometry).

Predict#

Before running: Which parameter has the strongest effect on \(L_c\) — shear modulus \(\mu\), strength drop \((\tau_s - \tau_d)\), or critical slip distance \(D_c\)?

def nucleation_length(mu, Dc, delta_tau, C=1.0):
    """
    Critical nucleation half-length for linear slip-weakening friction.

    Parameters
    ----------
    mu : float
        Shear modulus (Pa).
    Dc : float
        Critical slip distance (m).
    delta_tau : float
        Strength drop tau_s - tau_d (Pa).
    C : float
        Geometric constant (~1 for mode II).

    Returns
    -------
    Lc : float
        Critical nucleation half-length (m).
    """
    return C * mu * Dc / delta_tau


# Typical crustal values
mu = 30e9  # shear modulus, 30 GPa

# Sweep Dc and delta_tau
Dc_range = np.logspace(-3, 0, 50)       # 1 mm to 1 m
dtau_values = [1e6, 5e6, 10e6, 20e6]    # 1–20 MPa

fig, ax = plt.subplots(figsize=(8, 5))

for dtau in dtau_values:
    Lc = nucleation_length(mu, Dc_range, dtau)
    ax.loglog(Dc_range, Lc, lw=2,
              label=rf"$\tau_s - \tau_d = {dtau/1e6:.0f}$ MPa")

ax.set_xlabel(r"Critical slip distance $D_c$ (m)")
ax.set_ylabel(r"Nucleation half-length $L_c$ (m)")
ax.set_title(rf"Critical nucleation length ($\mu = {mu/1e9:.0f}$ GPa)")
ax.legend()

# Reference lines
ax.axhspan(0.1, 10, alpha=0.05, color="C2")
ax.text(2e-3, 2, "Lab scale", fontsize=9, color="C2")
ax.axhspan(10, 1000, alpha=0.05, color="C1")
ax.text(2e-3, 100, "Natural faults", fontsize=9, color="C1")

plt.tight_layout()
plt.show()

# Print some representative values
print("\nRepresentative nucleation lengths:")
for Dc_ex in [0.001, 0.01, 0.1, 1.0]:
    for dtau in [5e6, 10e6]:
        Lc = nucleation_length(mu, Dc_ex, dtau)
        print(f"  Dc = {Dc_ex:.3f} m, \u0394\u03c4 = {dtau/1e6:.0f} MPa  =>  Lc = {Lc:.1f} m")

Explain#

Interpretation:

  • For a typical crustal strength drop of ~5 MPa and \(D_c \approx 0.01\) m (a common estimate), what is the nucleation length?

  • Why does a larger \(D_c\) make nucleation harder? Connect your answer to the energy argument from the lecture.

  • What range of \(D_c\) gives nucleation lengths comparable to laboratory samples (~0.1 m)? To natural faults (~100 m)?

4. Energy Balance: Propagation vs. Arrest#

For a mode-II shear crack of half-length \(a\) under a uniform stress drop \(\Delta\sigma\), the energy release rate is:

\[\mathcal{G}(a) = \frac{\pi\,\Delta\sigma^2\,a}{\mu}\]

Rupture propagates when \(\mathcal{G}(a) \geq G_c\), where \(G_c = \frac{1}{2}(\tau_s - \tau_d)D_c\).

The critical crack size is:

\[a_c = \frac{\mu\,G_c}{\pi\,\Delta\sigma^2}\]

Predict#

Before running: If you increase the stress drop \(\Delta\sigma\), does the critical crack size \(a_c\) increase or decrease? Why?

def energy_release_rate(a, delta_sigma, mu):
    """Energy release rate for a mode-II crack."""
    return np.pi * delta_sigma**2 * a / mu


def critical_crack_size(mu, Gc, delta_sigma):
    """Critical crack half-length."""
    return mu * Gc / (np.pi * delta_sigma**2)


# Parameters
mu = 30e9          # shear modulus, 30 GPa
tau_s = 15e6       # 15 MPa static
tau_d = 10e6       # 10 MPa dynamic
Dc = 0.05          # 5 cm critical slip distance
Gc = 0.5 * (tau_s - tau_d) * Dc  # fracture energy

# Stress drops to compare
dsigma_values = [1e6, 3e6, 5e6, 10e6]  # 1–10 MPa
a = np.linspace(0, 2000, 500)  # crack half-lengths in meters

fig, ax = plt.subplots(figsize=(9, 5))

# Fracture energy line
ax.axhline(Gc, color="k", ls="--", lw=2, label=rf"$G_c = {Gc:.0f}$ J/m$^2$")

for dsig in dsigma_values:
    G = energy_release_rate(a, dsig, mu)
    ac = critical_crack_size(mu, Gc, dsig)
    ax.plot(a, G, lw=2,
            label=rf"$\Delta\sigma = {dsig/1e6:.0f}$ MPa, $a_c = {ac:.0f}$ m")
    ax.plot(ac, Gc, "o", ms=8, color="k", zorder=5)

ax.set_xlabel("Crack half-length $a$ (m)")
ax.set_ylabel(r"Energy release rate $\mathcal{G}$ (J/m$^2$)")
ax.set_title("Propagation vs. arrest")
ax.legend(fontsize=9)
ax.set_ylim(0, 5 * Gc)
ax.set_xlim(0, 2000)

plt.tight_layout()
plt.show()

# Summary table
print(f"\nFracture energy: G_c = {Gc:.0f} J/m\u00b2")
print(f"Shear modulus:   \u03bc = {mu/1e9:.0f} GPa")
print(f"D_c = {Dc*100:.0f} cm\n")
print(f"{'\u0394\u03c3 (MPa)':>12} {'a_c (m)':>12}")
print("-" * 26)
for dsig in dsigma_values:
    ac = critical_crack_size(mu, Gc, dsig)
    print(f"{dsig/1e6:12.0f} {ac:12.0f}")

Explain#

Interpretation:

  • How does \(a_c\) depend on \(\Delta\sigma\)? Why does higher stress drop give a smaller critical size?

  • For \(\Delta\sigma = 3\) MPa (a typical value), what is the critical crack size? Is this consistent with observed nucleation zone sizes?

  • How does this framework explain why small earthquakes are far more common than large ones?

5. Wrap-Up: Synthesis Questions#

Answer the following questions based on your explorations above.

  1. Spring–slider connection: In the spring–slider model, you found a critical stiffness separating stable from unstable behavior. How does this map to the nucleation length concept on a spatially extended fault?

  2. Parameter sensitivity: You explored how \(D_c\), \(\Delta\tau\), and \(\Delta\sigma\) affect nucleation and propagation. Rank these three in terms of their influence on whether an earthquake occurs, and justify your ranking.

  3. Energy budget bridge: In Lab 7d you computed radiated energy \(E_R\) and fracture energy \(G\) for a given earthquake. Using what you learned today, explain how \(G\) controls not just the energy budget of the earthquake, but also whether it nucleates at all.

  4. Real-world complexity: The models in this lab use homogeneous parameters. Name two features of real faults that would complicate nucleation, and describe qualitatively how each would modify your results.

Contents