# Lithospheric Strength Envelope Calculator

## Overview
This script calculates and visualizes the **strength envelope** and **viscosity** of the Earth's lithosphere and mantle at different depths. It models how rocks deform under stress and heat, showing which deformation mechanism dominates at different depths.


---
## Key Material Parameters

### **Plastic (Brittle) Material Parameters**

For brittle deformation (shallow crust), the strength is calculated as:

$$\tau = C + \mu \times \sigma_n = C + \tan(\text{friction angle}) \times \rho g \text{depth}$$

| Parameter | Symbol | Typical Range | Meaning |
|-----------|--------|---|---------|
| Density | ρ | 2700-3500 kg/m³ | Rock density; used to calculate pressure σ_n = ρ·g·depth |
| Cohesion | C | 10-100 MPa | Material strength at surface (when depth = 0); represents pre-existing structure |
| Friction angle | μ | 20-60° | Angle of internal friction; converted to coefficient: μ = tan(angle) |

:::success
**Parameter effects:**
- **Higher cohesion** → Stronger at shallow depths, weaker rocks are "stickier"
- **Higher friction angle** → Stronger at all depths, especially deep
- **Higher density** → Stronger at depth (due to higher normal stress)

:::

### **Viscous Material Parameters (Dislocation & Diffusion)**

#### **Viscosity:**

$$\eta_i = \frac{1}{2}A_i^{-1/n_i} d_i^{m_i/n_i} \dot{\epsilon}^{(1-n_i)/n_i} \exp\left(\frac{E_i + PV_i}{n_i RT}\right)$$

#### **Viscous Strength:**

$$\sigma_{viscous} = 2 \eta \times \dot{\epsilon}$$

| Parameter | Symbol | For Dislocation | For Diffusion | Meaning |
|-----------|--------|---|---|---------|
| Creep coefficient | A | A_disl | A_diff | Material constant; larger = weaker |
| Stress exponent | n | n_disl (~3) | n_diff (~1) | Controls stress dependence |
| Activation energy | E | E_disl | E_diff | Temperature sensitivity (J/mol) |
| Activation volume | V | V_disl | V_diff | Pressure sensitivity (m³/mol) |
| Grain size | d | d | d | Larger = stronger in dislocation |
| Grain size exponent | m | m_disl | m_diff | Controls grain size effect |

:::info
**Parameter Effects:**
**A (Creep Coefficient):**
- **Larger A** = weaker material (deforms more easily)
- **A_diff >> A_disl**: Diffusion is much weaker than dislocation at comparable stresses
- Example: If A_disl = 1e-23 and A_diff = 1e-15, then diffusion is 100 million times weaker!

**n (Stress Exponent):**
- **Higher n** = more stress-sensitive (strength changes rapidly with stress)
- Dislocation (n~3): Very nonlinear, highly stress-sensitive
- Diffusion (n~1): Nearly linear (Newtonian-like)

**E (Activation Energy):**
- **Higher E** = stronger temperature dependence
- At 1000°C increase in temperature:
- If E = 300 kJ/mol: viscosity changes by factor of ~10
- If E = 500 kJ/mol: viscosity changes by factor of ~100
- Larger E means deep, hot rocks are much weaker

**V (Activation Volume):**
- **Larger V** = more pressure-sensitive
- At 100 km depth increase:
- If V = 1e-5: significant strength increase
- If V = 5e-6: moderate strength increase

**d (Grain Size):**
- **Larger grains** in dislocation creep → Stronger material
- **Smaller grains** in diffusion creep → Weaker material
- Grain size evolution is important in real rocks (not constant here)

**m (Grain Size Exponent):**
- **Higher m** = grain size has stronger effect
- Typical m = 3: stress ∝ d³

:::

### Combining Multiple Mechanisms: Harmonic Mean

When both dislocation and diffusion operate:

$$\eta_{eff} = \left(\frac{1}{\eta_{diff}} + \frac{1}{\eta_{dis}}\right)^{-1}$$

**Physical meaning:**
- Both mechanisms can deform the material in parallel
- The weaker mechanism dominates
- But both contribute to total deformation

**With plastic:**
$$\eta_{total} = \left(\frac{1}{\eta_{eff}} + \frac{1}{\eta_{plast}}\right)^{-1}$$

### The Strain Rate Exponent: Critical Term

The key difference in the viscosity formula is: $\dot{\epsilon}^{(1-n_i)/n_i}$

#### **For Dislocation (n = 3):**
$$\eta_{disl} \propto \dot{\epsilon}^{-2/3}$$

- **Higher strain rate** → Lower viscosity → Weaker (easier to deform)
- **Lower strain rate** → Higher viscosity → Stronger
- This is **physically correct**: faster deformation is easier for dislocations

#### **For Diffusion (n = 1):**
$$\eta_{diff} \propto \dot{\epsilon}^{0} = 1$$

- **Independent of strain rate** → Newtonian (linear) viscosity
- Only temperature and pressure affect it

---

### **Global Temperature Profile**

| Parameter | Value | Unit | Meaning |
|-----------|-------|------|---------|
| `GEOTHERM_CRUST` | 30 | K/km | Temperature gradient in crust (steep) |
| `MOHO_DEPTH` | 35 | km | Crust-mantle boundary depth |
| `GEOTHERM_MANTLE` | 2.5 | K/km | Gradient in lithospheric mantle (moderate) |
| `ASTHENOSPHERE_TOP` | 100 | km | Base of lithosphere |
| `GEOTHERM_ASTHENOSPHERE` | 0.5 | K/km | Gradient in asthenosphere (adiabatic) |

## Python Script

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt

# ============================================================================
# CONTROL FLAGS
# ============================================================================
INCLUDE_DISLOCATION = True #False  # Set to True to include dislocation creep, False for diffusion only
# ============================================================================

# ============================================================================
# TEMPERATURE PROFILE PARAMETERS
# ============================================================================
SURFACE_TEMP = 273      # K

# Lithospheric (crustal) gradient
GEOTHERM_CRUST = 30     # K/km - steep gradient in crust

# Mohorovičić discontinuity depth (boundary between crust and mantle)
MOHO_DEPTH = 35e3       # m - typically 30-70 km depending on continental/oceanic

# Temperature at Moho
TEMP_AT_MOHO = SURFACE_TEMP + GEOTHERM_CRUST * (MOHO_DEPTH / 1000)  # K

# Mantle gradient (below Moho) - more gradual
GEOTHERM_MANTLE = 2.5   # K/km - much lower than crust; can be adiabatic (~0.3-0.5) or conductive
# For asthenosphere: can use adiabatic gradient ~0.3-0.5 K/km
# For lithospheric mantle: typically 0.5-2 K/km

# Optional: Define asthenosphere boundary (where deformation becomes dominant)
ASTHENOSPHERE_TOP = 100e3      # m - base of lithospheric mantle
GEOTHERM_ASTHENOSPHERE = 0.5    # K/km - adiabatic gradient (very small)

# ============================================================================

# Other parameters
STRAIN_RATE = 1e-13     # s^-1
R = 8.314
g = 9.81


def calculate_temperature(depth):
    """
    Calculate temperature at a given depth with boundary-dependent gradients.
    
    Regions:
    - 0 to MOHO_DEPTH: Crust (steep gradient)
    - MOHO_DEPTH to ASTHENOSPHERE_TOP: Lithospheric mantle (moderate gradient)
    - > ASTHENOSPHERE_TOP: Asthenosphere/deeper mantle (adiabatic/low gradient)
    """
    if depth <= MOHO_DEPTH:
        # Crustal region: simple linear gradient
        T = SURFACE_TEMP + GEOTHERM_CRUST * (depth / 1000)
    elif depth <= ASTHENOSPHERE_TOP:
        # Lithospheric mantle: from Moho to asthenosphere top
        T = TEMP_AT_MOHO + GEOTHERM_MANTLE * ((depth - MOHO_DEPTH) / 1000)
    else:
        # Asthenosphere and deeper mantle: adiabatic or very low gradient
        temp_at_asthenosphere_top = TEMP_AT_MOHO + GEOTHERM_MANTLE * ((ASTHENOSPHERE_TOP - MOHO_DEPTH) / 1000)
        T = temp_at_asthenosphere_top + GEOTHERM_ASTHENOSPHERE * ((depth - ASTHENOSPHERE_TOP) / 1000)
    
    return T


# Materials: [name, depth_min(m), depth_max(m), density, cohesion(MPa), friction(°), 
#             A_disl, n_disl, E_disl(J/mol), V_disl, d(m), m_disl,
#             A_diff, n_diff, E_diff(J/mol), V_diff, m_diff, color]
MATERIALS = [
    ['background', 0, 800e3, 3300, 20, 30,
     1e-24, 3.0, 5.0e5, 1.5e-5, 1e-3, 3.0,
     1.5e-15, 1.0, 375e3, 6e-6, 3.0, 'gray'],
    
    ['oceanic_crust_A', 0, 30e3, 3000, 50, 30,
     5e-24, 3.0, 4.5e5, 1.6e-5, 1e-3, 3.0,
     2.0e-13, 1.0, 350e3, 8e-6, 3.0, 'dodgerblue'],
    
    ['litho_mantle_A', 30e3, 100e3, 3400, 10, 30,
     1e-23, 3.5, 5.2e5, 1.4e-5, 1e-3, 3.0,
     1.0e-21, 1.0, 400e3, 5e-6, 3.0, 'navy'],
    
    ['oceanic_crust_B', 0e3, 30e3, 3000, 50, 30,
     5e-24, 3.0, 4.5e5, 1.6e-5, 1e-3, 3.0,
     2.0e-16, 1.0, 350e3, 8e-6, 3.0, 'green'],
    
    ['litho_mantle_B', 30e3, 100e3, 3400, 10, 30,
     1e-23, 3.5, 5.2e5, 1.4e-5, 1e-3, 3.0,
     1.0e-20, 1.0, 400e3, 5e-6, 3.0, 'lightgreen'],
    
    ['continental_crust_C', 0, 50e3, 2800, 60, 30,
     2e-24, 2.8, 3.2e5, 1.8e-5, 1e-3, 3.0,
     5.0e-24, 1.0, 320e3, 10e-6, 3.0, 'salmon'],
    
    ['continental_mantle_C', 50e3, 100e3, 3350, 50, 60,
     8e-24, 3.2, 5.0e5, 1.5e-5, 1e-3, 3.0,
     1.2e-23, 1.0, 380e3, 6e-6, 3.0, 'red'],
    
    ['upper_mantle', 100e3, 660e3, 3300, 40, 30,
     5e-24, 3.0, 5.0e5, 1.5e-5, 1e-3, 3.0,
     1e-25, 1.0, 375e3, 6e-6, 3.0, 'darkgoldenrod'],
    
    ['lower_mantle', 660e3, 800e3, 3500, 30, 30,
     2e-24, 3.2, 5.3e5, 1.4e-5, 1e-3, 3.0,
     3.0e-15, 1.0, 450e3, 8e-6, 3.0, 'maroon']
]

def calculate_viscosity_mechanism(depth, A, n, E, V, d, m, density, strain_rate):
    """
    Calculate viscosity for a single deformation mechanism using:
    η_i = (1/2) * A^(-1/n) * d^(m/n) * ε̇^((1-n)/n) * exp((E + PV)/(n*R*T))
    
    This formula is from the paper and accounts for the stress exponent n properly.
    """
    T = calculate_temperature(depth)
    P = density * g * depth
    
    # Calculate viscosity with proper strain rate dependence
    viscosity = (1/2) * (A ** (-1/n)) * (d ** (m/n)) * (strain_rate ** ((1-n)/n)) * \
                np.exp((E + V*P) / (n*R*T))
    
    return viscosity

def plastic_strength(depth, cohesion, friction, density):
    """Calculate plastic (brittle) strength."""
    sigma_n = density * g * depth / 1e6  # MPa
    mu = np.tan(np.radians(friction))
    return cohesion + mu * sigma_n

def plastic_viscosity(depth, cohesion, friction, density, strain_rate):
    """Convert plastic strength to viscosity."""
    stress_MPa = plastic_strength(depth, cohesion, friction, density)
    stress_Pa = stress_MPa * 1e6
    if strain_rate > 0:
        viscosity = stress_Pa / (2 * strain_rate)
    else:
        viscosity = 1e99
    return viscosity

def harmonic_mean_viscosity(eta1, eta2):
    """Calculate harmonic mean: (1/η1 + 1/η2)^-1"""
    if eta1 > 0 and eta2 > 0:
        return 1.0 / (1.0/eta1 + 1.0/eta2)
    elif eta1 > 0:
        return eta1
    else:
        return eta2

fig, (ax2, ax3) = plt.subplots(1, 2, figsize=(12, 7))

# Store handles and labels for legend
handles_list = []
labels_list = []

for mat in MATERIALS:
    name, d_min, d_max, rho, coh, fric, A_d, n_d, E_d, V_d, d, m_d, A_df, n_df, E_df, V_df, m_df, color = mat
    
    depths = np.linspace(d_min, d_max, 50)
    depths_km = depths / 1000
    
    # Calculate viscosities for each mechanism
    visc_diff = np.array([calculate_viscosity_mechanism(d_depth, A_df, n_df, E_df, V_df, d, m_df, rho, STRAIN_RATE) 
                          for d_depth in depths])
    visc_plast = np.array([plastic_viscosity(d_depth, coh, fric, rho, STRAIN_RATE) 
                          for d_depth in depths])
    
    if INCLUDE_DISLOCATION:
        # Dislocation creep
        visc_disl = np.array([calculate_viscosity_mechanism(d_depth, A_d, n_d, E_d, V_d, d, m_d, rho, STRAIN_RATE) 
                             for d_depth in depths])
        
        # Combined viscosity using harmonic mean (parallel flow)
        visc_creep_combined = np.array([harmonic_mean_viscosity(visc_diff[i], visc_disl[i]) 
                                        for i in range(len(depths))])
        
        # Effective viscosity: harmonic mean of creep and plastic
        visc_envelope = np.array([harmonic_mean_viscosity(visc_creep_combined[i], visc_plast[i]) 
                                 for i in range(len(depths))])
    else:
        # Effective viscosity: harmonic mean of diffusion and plastic
        visc_envelope = np.array([harmonic_mean_viscosity(visc_diff[i], visc_plast[i]) 
                                 for i in range(len(depths))])
    
    # Plot 1: Viscosity (log scale)
    line1 = ax2.semilogx(visc_diff, depths_km, color=color, linestyle=':', linewidth=2.5, marker='^',
                markersize=3, markevery=5, label=f'{name} - Diffusion', alpha=0.6)
    line2 = ax2.semilogx(visc_plast, depths_km, color=color, linestyle='--', linewidth=2, marker='s',
                markersize=3, markevery=5, label=f'{name} - Plastic', alpha=0.6)
    line3 = ax2.semilogx(visc_envelope, depths_km, color=color, linewidth=3.5, alpha=0.7,
                label=f'{name} - Envelope', linestyle='-')
    
    if INCLUDE_DISLOCATION:
        line4 = ax2.semilogx(visc_disl, depths_km, color=color, linestyle='-', linewidth=2, marker='o',
                    markersize=3, markevery=5, label=f'{name} - Dislocation', alpha=0.6)
    
    # Plot 2: Strength from stress calculation (for reference)
    def dislocation_strength(depth, A, n, E, V, d, m, density):
        T = calculate_temperature(depth)
        P = density * g * depth
        stress = (STRAIN_RATE / (A * d**m)) ** (1/n) * np.exp((E + V*P) / (n*R*T))
        return stress / 1e6  # MPa
    
    def diffusion_strength(depth, A, n, E, V, d, m, density):
        T = calculate_temperature(depth)
        P = density * g * depth
        stress = d**m * (STRAIN_RATE / A) * np.exp((E + V*P) / (n*R*T))
        return stress / 1e6  # MPa
    
    disl_stress = np.array([dislocation_strength(d_depth, A_d, n_d, E_d, V_d, d, m_d, rho) for d_depth in depths])
    diff_stress = np.array([diffusion_strength(d_depth, A_df, n_df, E_df, V_df, d, m_df, rho) for d_depth in depths])
    strength_plast = np.array([plastic_strength(d_depth, coh, fric, rho) for d_depth in depths])
    envelope_stress = np.minimum(np.minimum(disl_stress, diff_stress), strength_plast)
    
    ax3.plot(diff_stress, depths_km, color=color, linestyle=':', linewidth=2.5, marker='^',
            markersize=3, markevery=5, label=f'{name} - Diffusion', alpha=0.6)
    ax3.plot(strength_plast, depths_km, color=color, linestyle='--', linewidth=2, marker='s',
            markersize=3, markevery=5, label=f'{name} - Plastic', alpha=0.6)
    ax3.plot(envelope_stress, depths_km, color=color, linewidth=3.5, alpha=0.7,
            label=f'{name} - Envelope')
    
    if INCLUDE_DISLOCATION:
        ax3.plot(disl_stress, depths_km, color=color, linestyle='-', linewidth=2, marker='o',
                markersize=3, markevery=5, label=f'{name} - Dislocation', alpha=0.6)

# Add boundaries as horizontal lines
for ax in [ax2, ax3]:
    ax.axhline(y=MOHO_DEPTH/1000, color='black', linestyle=':', linewidth=1.5, alpha=0.5)
    ax.axhline(y=ASTHENOSPHERE_TOP/1000, color='purple', linestyle=':', linewidth=1.5, alpha=0.5)
    ax.invert_yaxis()
    ax.grid(True, alpha=0.3)
    ax.set_ylabel('Depth (km)', fontsize=12, fontweight='bold')
    ax.set_ylim(100, -5)

# Plot 1: Viscosity
ax2.set_xlabel('Viscosity (Pa·s)', fontsize=12, fontweight='bold')
ax2.set_title('Effective Viscosity (Log Scale)', fontsize=12, fontweight='bold')
ax2.set_xlim(1e18, 1e25)
ax2.grid(True, alpha=0.3, which='both')

# Plot 2: Stress Strength
ax3.set_xlabel('Shear Strength (MPa)', fontsize=12, fontweight='bold')
ax3.set_title('Strength Envelope', fontsize=12, fontweight='bold')
ax3.set_xlim(0, 1200)

# Get legend from ax2 and place it below the figure
handles, labels = ax2.get_legend_handles_labels()
fig.legend(handles, labels, loc='lower center', ncol=5, fontsize=8, 
          bbox_to_anchor=(0.5, -0.15), frameon=True, fancybox=True, shadow=False)

plt.tight_layout(rect=[0, 0.08, 1, 1])
plt.show()

---