Subduction

Plate-driven trench with slab descent (oceanic–continental / oceanic–oceanic).

Run it with:

julia -t auto --project=. main.jl cases/subduction.toml

Full case file — cases/subduction.toml:

# ══════════════════════════════════════════════════════════════════════════════
# Oceanic subduction case — two-plate symmetric convergence (two-plate symmetric setup).
#
# Domain: 4000 km × 1400 km (y = 0 surface, +y depth).
# Q2-Q1 Taylor-Hood FEM, MUMPS direct solver.
# ══════════════════════════════════════════════════════════════════════════════

[model]
case = "subduction"

# ── domain & resolution ──────────────────────────────────────────────────────
dim   = 2
xsize = 4000000.0      # 4000 km
ysize = 1400000.0      # 1400 km
nx    = 201            # 200 cells, 20 km
ny    = 71             # 70 cells, 20 km
mnx   = 4
mny   = 4

# ── gravity ──────────────────────────────────────────────────────────────────
gx = 0.0
gy = 9.81
pinit = 0.0

# ── rheology limits ──────────────────────────────────────────────────────────
eta_min = 1e19            # asthenosphere floor (was 1e18 — gave unphysical 1000 cm/yr flows)
eta_max = 1e25

# ── time control ─────────────────────────────────────────────────────────────
maxxystep = 0.5
maxtkstep = 20.0
maxtmstep = 1.5768e13   # ≈ 500 kyr
ntot      = 100

# ── physics toggles ──────────────────────────────────────────────────────────
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  = true
do_radiogenic    = true

# ── FEM discretization ───────────────────────────────────────────────────────
fe_order_v = 2          # Q2 velocity
fe_order_p = 1          # Q1 pressure (Taylor-Hood)
fe_order_T = 1
n_picard   = 3
picard_tol = 1e-3

# ── solver ───────────────────────────────────────────────────────────────────
solver_type = "mumps"
krylov_tol  = 1e-7
krylov_maxit = 200
gmg_levels  = 3

# ── IO ───────────────────────────────────────────────────────────────────────
output_dir   = "output/subduction"
output_every = 10
str_min      = 1e-18
eii_seed        = 1e-15   # mild seed — engages power-law without overheating
body_force_full = false   # Boussinesq perturbation (ρ−ρ̄)·g: stable with our pressure null-space
marker_T_mode   = "overwrite"   # subduction convection-light — PIC not needed

# ── thermal BCs ──────────────────────────────────────────────────────────────
T_top = 273.0
T_bot = 1623.0

# ── subduction geometry ──────────────────────────────────────────────────────
trench_x            = 2000000.0
plate_thickness     = 80000.0
plate_age           = 4e7
continent_thickness = 120000.0
air_thick           = 40000.0
x_plate_left        = 294000.0
x_plate_right       = 3706000.0
vx_plate            =  9.51173113506658e-10
vx_plate_right      = -9.51173113506658e-10

# ── free surface ─────────────────────────────────────────────────────────────
air_rock_id    = 1
surface_smooth = 0.0

# ══════════════════════════════════════════════════════════════════════════════
# Rocks (1=air, 2=oceanic crust, 3=continental crust, 4=lith mantle, 5=asthen)
# ══════════════════════════════════════════════════════════════════════════════

[[rock]]
nu0=1e18
ro0=1000.0
alpha=0.0
beta=0.0
cp=1000.0
kt=50.0
ht=0.0
mm=1.0
dh=0.0
dv=0.0
ss=0.0
a0=0.0
b0=0.0
a1=0.0
b1=0.0
e0=0.0
e1=1.0
nu_min=1e18
nu_max=1e19
immobile=false
G = 0.0

[[rock]]
nu0=1.97e17
dh=154000.0
dv=0.3
ss=3e4
mm=2.3
ro0=2900.0
alpha=3e-5
beta=0.0
cp=1000.0
kt=2.5
ht=2.5e-7
a0=3e6
b0=0.6
a1=3e5
b1=0.0
e0=0.0
e1=0.1
nu_min=1e18
nu_max=5e18
immobile=false
G = 4e10

[[rock]]
nu0=1e23
dh=0.0
dv=0.0
ss=0.0
mm=1.0
ro0=2700.0
alpha=3e-5
beta=0.0
cp=1000.0
kt=2.5
ht=1e-6
a0=1e7
b0=0.6
a1=1e6
b1=0.06
e0=0.0
e1=0.1
nu_min=1e18
nu_max=1e25
immobile=false
G = 3e10

[[rock]]
nu0=1.97e17
dh=154000.0
dv=0.3
ss=0.0
mm=1.0
ro0=3300.0
alpha=3e-5
beta=1e-11
cp=1000.0
kt=3.0
ht=2.5e-8
a0=3e7
b0=0.1
a1=3e6
b1=0.0
e0=0.0
e1=0.1
nu_min=1e22                # lithospheric mantle: strong (was 1e18; never softer than 10²²)
nu_max=1e25                # ceiling raised to engage MC plasticity at hinge
immobile=false
G = 7e10

[[rock]]
nu0=1.97e17
dh=154000.0
dv=0.3
ss=3e4
mm=2.3
ro0=3300.0
alpha=3e-5
beta=1e-11
cp=1000.0
kt=3.0
ht=2.5e-8
a0=1e6
b0=0.03
a1=1e5
b1=0.01
e0=0.0
e1=0.1
nu_min=1e19                # asthenosphere floor (raised from 1e18 — too weak)
nu_max=1e22                # max — allows some stiffness in deep regions
immobile=false
G = 7e10

# ── Phase transitions (Olivine → Wadsleyite at 410 km; eclogitization of basalt at 70 km) ──
[[phase]]
rock_from = 5     # asthenosphere
rock_to   = 5
transition = 1
T_ref = 1700.0
P_ref = 14e9      # ~410 km
clapeyron = 2.5e-9   # +2.5 MPa/K (exothermic)
drho = 270           # +270 kg/m³
latent = -90000      # -90 kJ/kg (released)

[[phase]]
rock_from = 2     # oceanic crust
rock_to   = 4     # transforms to eclogite (heavier — joins lith mantle)
transition = 2    # depth-triggered
depth = 70e3
drho = 200
latent = 0.0

# ── Melting (Katz-style for oceanic crust + asthenosphere) ──
[[melt]]
rock_id           = 2     # oceanic crust
T_solidus_ref     = 1085.0
T_liquidus_ref    = 1780.0
solidus_dPdT      = 132.9e-9
liquidus_dPdT     = 89.0e-9
H_latent_melt     = 4e5
extract_threshold = 0.3
rock_after_extract = 2

[[melt]]
rock_id           = 5     # asthenosphere
T_solidus_ref     = 1085.0
T_liquidus_ref    = 1780.0
solidus_dPdT      = 132.9e-9
liquidus_dPdT     = 89.0e-9
H_latent_melt     = 4e5
extract_threshold = 0.3
rock_after_extract = 5