Practice: Stress and Strain - EPS 164 Seismology Course Materials

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

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:

# 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!")

Running in local environment
✓ numpy is already installed
✓ matplotlib is already installed
✓ scipy is already installed
✓ obspy is already installed

✓ All required packages are installed!

Part (a): Calculate Lamé Parameters¶

The University of California, San Diego, operates the Piñon Flat Observatory (PFO) in the mountains northeast of San Diego (near Anza). The observatory serves as a “testbed” for the global geophysical community, housing some of the world’s most sensitive instrumentation. It is perhaps best known for its long-base laser strainmeters, which can detect shifts in the earth as small as the width of an atom. These instruments, along with deep borehole sensors and advanced seismic arrays, allow researchers to monitor the gradual buildup of tectonic strain that precedes major earthquakes. Instruments include high-quality strain meters for measuring crustal deformation. 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:

Calculate the shear modulus $\mu$ using $\mu = \rho v_s^2$
Calculate the first Lamé parameter $\lambda$ using $\lambda = \rho v_p^2 - 2\rho v_s^2$
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')

  Cell In[3], line 2
    vp =     # m/s
             ^
SyntaxError: invalid syntax

Part (b): Stress from Landers Earthquake¶

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).

$\sigma_{11} = \lambda(\epsilon_{11} + \epsilon_{22}) + 2\mu\epsilon_{11}$
$\sigma_{22} = \lambda(\epsilon_{11} + \epsilon_{22}) + 2\mu\epsilon_{22}$
$\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¶

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

Tasks:

Construct the 2D stress tensor as a NumPy array
Calculate eigenvalues and eigenvectors using numpy.linalg.eig
Report the principal stress magnitudes (eigenvalues)
Plot the principal stress orientations by either hand-sketching on paper or by using the demo notebook. Make sure your plot:
- 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}

(1)

Hint: The formula for azimuth from eigenvector [v₁, v₂] is $\text{azimuth} = \arctan2(v_1, v_2) \times \frac{180°}{\pi}$ . In python, 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¶

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:

Calculate the accumulated total strain: $\epsilon_{ij} = \dot{\epsilon}_{ij} \times 1000$ yr
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¶

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)

(2)

where $L_0 = 1000$ m.

Tasks:

Calculate area changes for:

(e1) Annual change from tectonic deformation¶

Use the annual strain rates (not 1000-year accumulation)
Construct strain tensor and compute eigenvalues
Calculate new area and report the change in m²

(e2) Change from Landers earthquake¶

Use the Landers strain values from Part (b)
Construct strain tensor and compute eigenvalues
Calculate new area and report the change in m²
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¶

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}

(3)

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

Tasks:

Calculate $\tau_{\text{max}}$ for the Landers earthquake stress field
Compare to typical fault friction stress (~10–30 MPa): Is the Landers stress change significant?
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.

# 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)
✅ 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.