Stress and Strain

Stress and Strain#

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

Learning Objectives:

  • Define stress and strain tensors

  • Rotate stress tensors to principal axes

  • Relate stress and strain via Hooke’s Law

  • Visualize tensor deformation

Prerequisites: Linear algebra (matrices, eigenvalues), Calculus

Reference: Shearer, Chapter 2 (Stress and Strain)

Notebook Outline:

  • Time:** 50 minutes

  • Origin**: This notebook came out Heiner Igel’s seismolive project.

Hide code cell source

 
# Install dependencies (for Google Colab or missing packages)
import sys

# Check if running in Colab
try:
    import google.colab
    IN_COLAB = True
    print("Running in Google Colab")
except:
    IN_COLAB = False
    print("Running in local environment")

# Install required packages if needed
required_packages = {
    'numpy': 'numpy',
    'matplotlib': 'matplotlib',
    'scipy': 'scipy',
    'obspy': 'obspy'
}

missing_packages = []
for package, pip_name in required_packages.items():
    try:
        __import__(package)
        print(f"✓ {package} is already installed")
    except ImportError:
        missing_packages.append(pip_name)
        print(f"✗ {package} not found")

if missing_packages:
    print(f"\nInstalling missing packages: {', '.join(missing_packages)}")
    import subprocess
    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q"] + missing_packages)
    print("✓ Installation complete!")
else:
    print("\n✓ All required packages are installed!")

Hide code cell source

 
# Imports
import matplotlib.pyplot as plt
import numpy as np
from numpy import linalg as la

Part (a): Calculate Lamé Parameters#

Time estimate: 5–7 minutes

At 5 km depth beneath PFO, seismic tomography gives:

  • P-wave velocity: \(v_p = 6.0\) km/s

  • S-wave velocity: \(v_s = 3.5\) km/s

  • Density: \(\rho = 2700\) kg/m³

Tasks:

  1. Calculate the shear modulus \(\mu\) using \(\mu = \rho v_s^2\)

  2. Calculate the first Lamé parameter \(\lambda\) using \(\lambda = \rho v_p^2 - 2\rho v_s^2\)

  3. Convert both values to GPa

Formulas provided:

  • \(v_p = \sqrt{\frac{\lambda + 2\mu}{\rho}}\)

  • \(v_s = \sqrt{\frac{\mu}{\rho}}\)

# Given values (CONVERT TO SI UNITS)
vp =     # m/s
vs =     # m/s
ro =     # kg/m**3

# Calculate Lamé parameters
mu =     # Pa
lam =    # Pa

print('Lamé Parameters:')
print(f'  μ (mu)    = {mu:.3e} Pa = {mu/1e9:.2f} GPa')
print(f'  λ (lambda) = {lam:.3e} Pa = {lam/1e9:.2f} GPa')

Part (b): Stress from Landers Earthquake#

Time estimate: 7–9 minutes

After the 1992 Landers earthquake (M7.3), GPS measurements 80 km north of PFO recorded the following horizontal strain components:

  • \(\epsilon_{11} = -0.26 \times 10^{-6}\) (East-East)

  • \(\epsilon_{12} = -0.69 \times 10^{-6}\) (East-North shear)

  • \(\epsilon_{22} = +0.92 \times 10^{-6}\) (North-North)

Indices: 1 = East, 2 = North

Tasks:

Calculate the stress changes at 5 km depth using the 2D stress-strain relationship and your Lamé parameters from Part (a).

  1. \(\sigma_{11} = \lambda(\epsilon_{11} + \epsilon_{22}) + 2\mu\epsilon_{11}\)

  2. \(\sigma_{22} = \lambda(\epsilon_{11} + \epsilon_{22}) + 2\mu\epsilon_{22}\)

  3. \(\sigma_{12} = 2\mu\epsilon_{12}\)

Formulas provided:

  • General form: \(\sigma_{ij} = \lambda \epsilon_{kk}\delta_{ij} + 2\mu\epsilon_{ij}\)

  • Trace: \(\epsilon_{kk} = \epsilon_{11} + \epsilon_{22}\) (in 2D)

# Strain components from Landers earthquake
e11_landers =     # dimensionless
e12_landers =     # dimensionless
e22_landers =     # dimensionless

# Calculate stress changes (use lam and mu from Part a)
s11_landers =     # Pa
s22_landers =     # Pa
s12_landers =     # Pa

print('Stress Changes from Landers Earthquake:')
print(f'  σ₁₁ = {s11_landers:.3e} Pa = {s11_landers/1e3:.2f} kPa')
print(f'  σ₂₂ = {s22_landers:.3e} Pa = {s22_landers/1e3:.2f} kPa')
print(f'  σ₁₂ = {s12_landers:.3e} Pa = {s12_landers/1e3:.2f} kPa')

Part (c): Principal Stress Directions#

Time estimate: 10–15 minutes

The stress tensor from Part (b) has principal axes where only normal stresses act (no shear). Find these directions using eigenvalue decomposition.

Tasks:

  1. Construct the 2D stress tensor as a NumPy array

  2. Calculate eigenvalues and eigenvectors using numpy.linalg.eig

  3. Report the principal stress magnitudes (eigenvalues)

  4. Hand-sketch the principal stress orientations on paper:

    • Draw East-North coordinate axes

    • Plot both eigenvectors as arrows from the origin

    • Label each arrow with its corresponding eigenvalue

    • Calculate and label the azimuth (angle clockwise from North) of the maximum principal stress

Formula for azimuth from eigenvector [v₁, v₂]: $\(\text{azimuth} = \arctan2(v_1, v_2) \times \frac{180°}{\pi}\)$

Hint: Use np.degrees(np.arctan2(v[0,i], v[1,i])) where i is the eigenvector index

# Initialize 2D stress tensor
sigma = np.zeros((2,2))
sigma[0,0] =     # σ₁₁
sigma[1,1] =     # σ₂₂
sigma[0,1] =     # σ₁₂
sigma[1,0] =     # σ₂₁ (must equal σ₁₂ for symmetric tensor)

print('Stress Tensor:')
print(sigma)

# Calculate eigenvalues and eigenvectors
w, v =     # Use numpy.linalg.eig

# Report results
print('\nPrincipal Stresses (Eigenvalues):')
print(f'  σ₁ = {w[0]:.3e} Pa = {w[0]/1e3:.2f} kPa')
print(f'  σ₂ = {w[1]:.3e} Pa = {w[1]/1e3:.2f} kPa')

print('\nPrincipal Stress Directions (Eigenvectors):')
print(f'  Direction 1: [{v[0,0]:.4f}, {v[1,0]:.4f}]')
print(f'  Direction 2: [{v[0,1]:.4f}, {v[1,1]:.4f}]')

# Calculate azimuths (angle clockwise from North)
azimuth_1 =     # degrees
azimuth_2 =     # degrees

print(f'\nAzimuths (clockwise from North):')
print(f'  σ₁ direction: {azimuth_1:.1f}°')
print(f'  σ₂ direction: {azimuth_2:.1f}°')

HAND-SKETCH REQUIRED: On paper, draw the coordinate system and principal stress arrows. Attach this sketch to your submission.

Part (d): Long-Term Tectonic Stress Accumulation#

Time estimate: 6–8 minutes

Continuous GPS at PFO measures annual strain rates:

  • \(\dot{\epsilon}_{11} = +0.101 \times 10^{-6}\) /yr

  • \(\dot{\epsilon}_{12} = +0.005 \times 10^{-6}\) /yr

  • \(\dot{\epsilon}_{22} = -0.02 \times 10^{-6}\) /yr

Tasks:

Assuming these rates continue for 1000 years without earthquakes:

  1. Calculate the accumulated total strain: \(\epsilon_{ij} = \dot{\epsilon}_{ij} \times 1000\) yr

  2. Calculate the resulting stress changes using the same formulas as Part (b)

    • \(\sigma_{11}\), \(\sigma_{22}\), \(\sigma_{12}\)

# Annual strain rates (per year)
e11_rate =     # /yr
e12_rate =     # /yr
e22_rate =     # /yr

# Time period
time_years =     # years

# Calculate accumulated strain
e11_1000yr =     # dimensionless
e12_1000yr =     # dimensionless
e22_1000yr =     # dimensionless

# Calculate stress changes
s11_1000yr =     # Pa
s22_1000yr =     # Pa
s12_1000yr =     # Pa

print('Stress Changes after 1000 years:')
print(f'  σ₁₁ = {s11_1000yr:.3e} Pa = {s11_1000yr/1e6:.2f} MPa')
print(f'  σ₂₂ = {s22_1000yr:.3e} Pa = {s22_1000yr/1e6:.2f} MPa')
print(f'  σ₁₂ = {s12_1000yr:.3e} Pa = {s12_1000yr/1e6:.2f} MPa')

Part (e): Land Area Change#

Time estimate: 8–12 minutes

A farmer owns 1 km² of land (1000 m × 1000 m) near PFO. Ground deformation changes the area!

Background: The eigenvalues of the strain tensor represent fractional length changes along principal axes. If \(\lambda_1\) and \(\lambda_2\) are the strain eigenvalues:

\[\text{New Area} = (L_0 + L_0\lambda_1)(L_0 + L_0\lambda_2) = L_0^2(1 + \lambda_1)(1 + \lambda_2)\]

where \(L_0 = 1000\) m.

Tasks:

Calculate area changes for:

(e1) Annual change from tectonic deformation#

  1. Use the annual strain rates (not 1000-year accumulation)

  2. Construct strain tensor and compute eigenvalues

  3. Calculate new area and report the change in m²

(e2) Change from Landers earthquake#

  1. Use the Landers strain values from Part (b)

  2. Construct strain tensor and compute eigenvalues

  3. Calculate new area and report the change in m²

  4. Compare: Which caused more area change?

# ===== (e1) Annual Tectonic Deformation =====

# Annual strain rates (same as Part d, but for 1 year only)
e11_annual =     # /yr × 1 yr = dimensionless
e12_annual =     # /yr × 1 yr = dimensionless
e22_annual =     # /yr × 1 yr = dimensionless

# Construct strain tensor
epsilon_annual = np.array([[e11_annual, e12_annual],
                           [e12_annual, e22_annual]])

print('Annual Strain Tensor:')
print(epsilon_annual)

# Compute eigenvalues
eigenvalues_annual, _ =     # Use numpy.linalg.eig

print(f'\nStrain Eigenvalues (annual):')
print(f'  λ₁ = {eigenvalues_annual[0]:.6e}')
print(f'  λ₂ = {eigenvalues_annual[1]:.6e}')

# Calculate area change
L0 = 1000  # m
old_area = L0 * L0
new_area_annual =     # m²
area_change_annual = new_area_annual - old_area

print(f'\nOld area: {old_area} m²')
print(f'New area: {new_area_annual:.2f} m²')
print(f'Area change: {area_change_annual:.4f} m² per year')

# ===== (e2) Landers Earthquake =====

# Landers strain components (from Part b)
e11_landers_strain =     # dimensionless
e12_landers_strain =     # dimensionless
e22_landers_strain =     # dimensionless

# Construct strain tensor
epsilon_landers = np.array([[e11_landers_strain, e12_landers_strain],
                            [e12_landers_strain, e22_landers_strain]])

print('Landers Strain Tensor:')
print(epsilon_landers)

# Compute eigenvalues
eigenvalues_landers, _ =     # Use numpy.linalg.eig

print(f'\nStrain Eigenvalues (Landers):')
print(f'  λ₁ = {eigenvalues_landers[0]:.6e}')
print(f'  λ₂ = {eigenvalues_landers[1]:.6e}')

# Calculate area change
new_area_landers =     # m²
area_change_landers = new_area_landers - old_area

print(f'\nOld area: {old_area} m²')
print(f'New area: {new_area_landers:.2f} m²')
print(f'Area change: {area_change_landers:.4f} m² (earthquake)')

Comparison: Which caused a larger area change—one year of tectonic deformation or the Landers earthquake?

Your answer here:

Part (f): Maximum Shear Stress — ESS 512 ONLY#

Time estimate: 10–12 minutes

Graduate students only: Fault slip is driven by shear stress on the fault plane. The maximum shear stress in 2D is:

\[\tau_{\text{max}} = \frac{|\sigma_1 - \sigma_2|}{2}\]

where \(\sigma_1\) and \(\sigma_2\) are the principal stresses (eigenvalues from Part c).

Tasks:

  1. Calculate \(\tau_{\text{max}}\) for the Landers earthquake stress field

  2. Compare to typical fault friction stress (~10–30 MPa): Is the Landers stress change significant?

  3. Bonus conceptual question: If we assume plane stress (\(\sigma_{33} = 0\), no vertical stress), would the 3D maximum shear stress be larger or smaller than your 2D result? Explain briefly.

# ===== ESS 512 ONLY =====

# Principal stresses from Part (c)
sigma_1 =     # Pa (larger eigenvalue)
sigma_2 =     # Pa (smaller eigenvalue)

# Calculate maximum shear stress
tau_max =     # Pa

print(f'Maximum Shear Stress:')
print(f'  τ_max = {tau_max:.3e} Pa = {tau_max/1e6:.3f} MPa')

# Comparison to fault friction
fault_friction_typical = 20e6  # Pa (20 MPa, mid-range estimate)
ratio = tau_max / fault_friction_typical

print(f'\nComparison to typical fault friction (~20 MPa):')
print(f'  Ratio: {ratio:.4f} (or {ratio*100:.2f}% of frictional strength)')

Interpretation: Is the Landers-induced shear stress significant for triggering aftershocks?

Your answer here:

Plane stress conceptual question: If \(\sigma_{33} = 0\), how would 3D \(\tau_{\text{max}}\) compare to 2D?

Your answer here:

Summary & Reflection#

You’ve now calculated:

  • ✅ Elastic constants from seismic velocities

  • ✅ Stress changes from measured strain

  • ✅ Principal stress orientations

  • ✅ Long-term vs. earthquake deformation

  • ✅ Practical applications (land area changes)

  • ✅ (ESS 512) Maximum shear stress and fault mechanics

Key Takeaway: GPS strain measurements (~\(10^{-6}\)) translate to stress changes of kPa to MPa, which accumulate over decades and can trigger earthquakes when frictional strength is exceeded.

End of Test

Contents