Continental collision
Convergent orogeny — crustal thickening and shortening.
Run it with:
bash
julia -t auto --project=. main.jl cases/collision.tomlFull case file — cases/collision.toml:
toml
# ══════════════════════════════════════════════════════════════════════════════
# Continental collision — two blocks converging at 6 cm/yr total.
#
# Domain: 2000 × 400 km.
# Left plate (x 0–900 km) and right plate (x 1100–2000 km) each push
# inward at 3 cm/yr; a 200 km wide suture zone in the centre localises
# deformation and decouples the two plates.
#
# Rocks:
# 1 = upper crust (felsic, 25 km, C=10 MPa, φ=30°)
# 2 = lower crust (mafic, 15 km, C=20 MPa, φ=30°)
# 3 = lith. mantle (olivine, 80 km, C=40 MPa, φ=37°)
# 4 = asthenosphere (hot, η=5×10¹⁹ Pa·s, weak plastic floor)
# 5 = suture zone (very weak, C=2 MPa, φ=6°, decoupler)
#
# Expected physics:
# • Upper crust forms fold-and-thrust belt + pop-up structure.
# • Lower crust channels into the orogen core.
# • Lithospheric mantle underthrusts / possibly delaminates.
# • Radiogenic + shear heating warm the thickened crust.
# ══════════════════════════════════════════════════════════════════════════════
[model]
case = "collision"
dim = 2
geometry = "box"
xsize = 2000000.0 # 2000 km
ysize = 400000.0 # 400 km (120 km lithosphere + 280 km asthenosphere)
nx = 101
ny = 21
mnx = 3
mny = 3
gx = 0.0
gy = 9.81
pinit = 0.0
eta_min = 1e18
eta_max = 1e24
str_min = 1e-20
maxxystep = 0.5
maxtkstep = 50.0
maxtmstep = 1.5778e13 # 500 kyr cap per step (keeps each step mildly nonlinear)
ntot = 200 # ~20 Myr
do_stokes = true
do_heat = true
do_advect = true
do_adiabatic = false
do_shear_heating = true
do_friction = true
do_phase_change = false
do_melting = false
do_free_surface = false
do_radiogenic = true
do_ale_free_surface = false
heat_stabilisation = "entropy"
fe_order_v = 2
fe_order_p = 1
fe_order_T = 1
n_picard = 10
picard_tol = 1e-4 # viscoplastic collision needs looser tol than purely viscous cases
nonlinear_method = "picard"
newton_damping = 1.0
solver_type = "mumps"
krylov_tol = 1e-7
krylov_maxit = 200
gmg_levels = 3
output_dir = "output/collision"
output_every = 10
T_top = 273.0
T_bot = 1600.0
T_ref_density = 273.0
marker_T_mode = "pic_average"
# ── Collision kinematics ───────────────────────────────────────────────────────
# Each plate pushes inward at 3 cm/yr (6 cm/yr total convergence rate).
# plate_thickness / continent_thickness = depth of the Dirichlet BC window.
# bc_deep = "prescribed" auto-computes mass-conserving bottom outflow.
bc_left = "prescribed"
bc_right = "prescribed"
bc_surface = "free_slip"
bc_deep = "prescribed"
vx_plate = 9.5129e-10 # m/s = +3 cm/yr (left plate → rightward)
vx_plate_right = -9.5129e-10 # m/s = -3 cm/yr (right plate → leftward)
plate_thickness = 120000.0 # m — full lithosphere on left side
continent_thickness = 120000.0 # m — full lithosphere on right side
trench_x = 1000000.0 # domain centre (Config compat.; unused by BCs)
x_plate_left = 0.0
x_plate_right = 2000000.0
plate_age = 0.0
air_thick = 0.0
air_rock_id = 1
surface_smooth = 0.0
# ── Rocks ──────────────────────────────────────────────────────────────────────
# Yield criterion (Drucker-Prager): τ_yield = C(ε) + sin(φ(ε)) · P
# Strain-softening: C and sin(φ) weaken linearly from e0 → e1 accumulated strain.
# Field mapping: a0=cohesion_init, a1=cohesion_weak, b0=friction_init (=sin φ),
# b1=friction_weak, e0=strain_thresh_init, e1=strain_thresh_weak.
# ── Rock 1: upper continental crust (felsic, 0–25 km) ─────────────────────────
[[rock]]
nu0 = 1e22
ro0 = 2700.0
alpha = 3e-5
beta = 0.0
cp = 1000.0
kt = 2.5
ht = 1.5e-6 # W/m³ — high radiogenic heat (K-rich felsic crust)
mm = 1.0
dh = 0.0
dv = 0.0
ss = 0.0
a0 = 10e6 # cohesion_init 10 MPa
a1 = 2e6 # cohesion_weak 2 MPa
b0 = 0.5 # sin(φ_init) = 0.5 → φ ≈ 30°
b1 = 0.1 # sin(φ_weak) = 0.1 → φ ≈ 6° (strain-softened)
e0 = 0.0
e1 = 0.5
nu_min = 1e18
nu_max = 1e23
immobile = false
G = 3e10
# ── Rock 2: lower continental crust (mafic / granulite, 25–40 km) ─────────────
[[rock]]
nu0 = 1e23
ro0 = 2950.0
alpha = 3e-5
beta = 0.0
cp = 1000.0
kt = 2.5
ht = 3e-7 # W/m³ — lower radiogenic (mafic)
mm = 1.0
dh = 0.0
dv = 0.0
ss = 0.0
a0 = 20e6 # cohesion_init 20 MPa
a1 = 5e6 # cohesion_weak 5 MPa
b0 = 0.5
b1 = 0.1
e0 = 0.0
e1 = 0.5
nu_min = 1e19
nu_max = 1e24
immobile = false
G = 4e10
# ── Rock 3: lithospheric mantle (cold, strong olivine, 40–120 km) ─────────────
[[rock]]
nu0 = 1e23
ro0 = 3300.0
alpha = 3e-5
beta = 1e-11
cp = 1200.0
kt = 3.0
ht = 2e-8 # W/m³ — trace radiogenic
mm = 1.0
dh = 0.0
dv = 0.0
ss = 0.0
a0 = 40e6 # cohesion_init 40 MPa
a1 = 10e6 # cohesion_weak 10 MPa
b0 = 0.6 # sin(φ_init) = 0.6 → φ ≈ 37° (Byerlee)
b1 = 0.1
e0 = 0.0
e1 = 0.5
nu_min = 1e20
nu_max = 1e24
immobile = false
G = 7e10
# ── Rock 4: asthenosphere (hot convecting mantle, >120 km) ────────────────────
[[rock]]
nu0 = 5e19
ro0 = 3300.0
alpha = 3e-5
beta = 1e-11
cp = 1200.0
kt = 3.0
ht = 2e-8
mm = 1.0
dh = 0.0
dv = 0.0
ss = 0.0
a0 = 1e6
a1 = 1e6
b0 = 0.1
b1 = 0.05
e0 = 0.0
e1 = 1.0
nu_min = 1e18
nu_max = 1e21
immobile = false
G = 7e10
# ── Rock 5: suture / collision zone (very weak, inherited oceanic relict) ──────
[[rock]]
nu0 = 1e20
ro0 = 2800.0
alpha = 3e-5
beta = 0.0
cp = 1000.0
kt = 2.5
ht = 5e-7
mm = 1.0
dh = 0.0
dv = 0.0
ss = 0.0
a0 = 2e6 # cohesion_init 2 MPa — yields immediately at the suture
a1 = 0.5e6 # cohesion_weak 0.5 MPa
b0 = 0.1 # sin(φ_init) = 0.1 → φ ≈ 6°
b1 = 0.02
e0 = 0.0
e1 = 0.1
nu_min = 1e18
nu_max = 1e21
immobile = false
G = 2e10
# ── Initial geometry ───────────────────────────────────────────────────────────
# y = 0 at TOP (surface), y = 400 km at BOTTOM.
# Evaluated in order — LATER entries override EARLIER ones.
[[layer]]
name = "asthenosphere"
rock_id = 4
shape = "background"
[[layer]]
name = "left_lith_mantle"
rock_id = 3
shape = "rectangle"
x_min = 0.0
x_max = 900000.0
y_min = 40000.0
y_max = 120000.0
[[layer]]
name = "right_lith_mantle"
rock_id = 3
shape = "rectangle"
x_min = 1100000.0
x_max = 2000000.0
y_min = 40000.0
y_max = 120000.0
[[layer]]
name = "left_lower_crust"
rock_id = 2
shape = "rectangle"
x_min = 0.0
x_max = 900000.0
y_min = 25000.0
y_max = 40000.0
[[layer]]
name = "right_lower_crust"
rock_id = 2
shape = "rectangle"
x_min = 1100000.0
x_max = 2000000.0
y_min = 25000.0
y_max = 40000.0
[[layer]]
name = "left_upper_crust"
rock_id = 1
shape = "rectangle"
x_min = 0.0
x_max = 900000.0
y_min = 0.0
y_max = 25000.0
[[layer]]
name = "right_upper_crust"
rock_id = 1
shape = "rectangle"
x_min = 1100000.0
x_max = 2000000.0
y_min = 0.0
y_max = 25000.0
[[layer]]
name = "suture_zone"
rock_id = 5
shape = "rectangle"
x_min = 900000.0
x_max = 1100000.0
y_min = 0.0
y_max = 40000.0 # suture replaces crust in the 200-km collision front
# ── Initial temperature ────────────────────────────────────────────────────────
# Background: warm adiabat (asthenosphere, T=1300 K at y=0 to 1600 K at y=400 km).
# Continental plates: linear geotherm 273 K (surface) → 1573 K (LAB at 120 km).
# Suture zone: same geotherm as surrounding crust (linear 273→1573 K over 120 km).
# Later [[init_T]] entries override earlier ones at overlapping points.
[[init_T]]
name = "background_adiabat"
shape = "background"
profile = "adiabatic"
T_surf = 1300.0
T_mantle = 1600.0
[[init_T]]
name = "left_continental_geotherm"
shape = "rectangle"
profile = "continental_geotherm"
x_min = 0.0
x_max = 900000.0
y_min = 0.0
y_max = 120000.0
T_surf = 273.0
T_mantle = 1573.0
[[init_T]]
name = "right_continental_geotherm"
shape = "rectangle"
profile = "continental_geotherm"
x_min = 1100000.0
x_max = 2000000.0
y_min = 0.0
y_max = 120000.0
T_surf = 273.0
T_mantle = 1573.0
[[init_T]]
name = "suture_geotherm"
shape = "rectangle"
profile = "linear"
x_min = 900000.0
x_max = 1100000.0
y_min = 0.0
y_max = 120000.0
T_top = 273.0
T_bot = 1573.0