Lithospheric drip & delamination
A dense lithospheric root founders / peels away (Holt & Condit 2021).
Run it with:
bash
julia -t auto --project=. main.jl cases/holt_condit_2021.tomlFull case file — cases/holt_condit_2021.toml:
toml
# ══════════════════════════════════════════════════════════════════════════════
# Holt & Condit (2021) "Slab temperature evolution over the lifetime of a
# subduction zone" — reference thermo-mechanical model.
# https://zenodo.org/records/4543414
#
# Full thermo-mechanical setup with the Hirth-Kohlstedt composite diffusion+
# dislocation Arrhenius creep, 90-Myr SP and 10-Myr OP, large lower-mantle
# domain. Scaled DOWN from Holt-Condit's 11600×2900 km (which needs AMR +
# 100 MPI ranks) to 5000×1500 km — still big but tractable on uniform mesh.
#
# Reference paper figures show slab tip reaching the upper-mantle base
# (660 km depth, in our 1500-km domain that's well above the bottom) in
# ~20-30 Myr, then continuing into the deeper mantle while folding /
# stagnating below 660 km.
# ══════════════════════════════════════════════════════════════════════════════
[model]
case = "subduction"
[mesh]
dim = 2
xsize = 5000000.0 # 5000 km — scaled down from H-C 11600 km
ysize = 1500000.0 # 1500 km — covers upper mantle + bit of lower
nx = 251 # 250 cells × 20 km
ny = 76 # 75 cells × 20 km
mnx = 4
mny = 4
[gravity]
gy = 9.81
[time]
maxxystep = 0.5
maxtkstep = 20.0
maxtmstep = 2e13 # ~650 kyr cap
ntot = 1000
[rheology]
eta_min = 2.5e18 # H-C Minimum Viscosity
eta_max = 2.5e24 # H-C Maximum Viscosity
str_min = 1e-20
eii_seed = 1e-15
[temperature]
T_top = 273.0
T_bot = 1694.5 # H-C reference / bottom T
T_ref_density = 273.0
marker_T_mode = "overwrite"
[physics]
do_stokes = true
do_heat = true
do_advect = true
do_adiabatic = false
do_shear_heating = false
do_friction = true
do_phase_change = false
do_melting = false
do_free_surface = false
do_radiogenic = false
[boundary_velocity]
bc_left = "free_slip" # H-C: tangential velocity on all four walls
bc_right = "free_slip"
bc_surface = "free_slip"
bc_deep = "free_slip" # H-C does free-slip bottom (not no-slip)
[free_surface]
do_ale_free_surface = false
air_rock_id = 1
[amr]
do_amr = false
[fe]
fe_order_v = 2
fe_order_p = 1
fe_order_T = 1
n_picard = 10
picard_tol = 1e-5
[solver]
type = "mumps"
[output]
dir = "output/holt_condit_2021"
every = 10
[subduction]
trench_x = 2000000.0
plate_thickness = 100000.0
continent_thickness = 50000.0
air_thick = 0.0
plate_age = 90e6
x_plate_left = 300000.0
x_plate_right = 4700000.0
vx_plate = 0.0
vx_plate_right = 0.0
# ══════════════════════════════════════════════════════════════════════════════
# Rocks — Holt & Condit 2021 reference (4 compositions)
# ══════════════════════════════════════════════════════════════════════════════
# Composite η = η_diff + η_disl, both Arrhenius. Dislocation params for
# dry olivine (Hirth & Kohlstedt 2003), diffusion params lower (typical
# diffusion creep activation energy).
# 1 = background mantle (dry olivine, with plasticity)
# 2 = SP crust (weak η=2e20, no plasticity — decoupling layer)
# 3 = OP (dry olivine, NO plasticity)
# 4 = SP mantle (dry olivine, with plasticity)
# All rocks have same ρ_ref = 3300; density CONTRAST comes purely from T.
[[rock]] # 1: background mantle
eta0 = 1e16 # reference viscosity (will be modulated by Arrhenius)
activation_energy = 540000.0 # Hirth-Kohlstedt dry olivine Q_disl
activation_volume = 1.2 # 12 cm³/mol × 1e5 (bar→Pa)
stress_exponent = 3.5
ref_yield_stress = 1e7
reference_density = 3300.0
thermal_expansion = 3e-5
specific_heat = 940.0 # H-C value
thermal_conductivity = 2.96 # κρcp = κ_th = 1e-6 with these values
cohesion_init = 5e7 # 50 MPa
cohesion_weak = 5e7
friction_init = 0.342 # sin(20°)
friction_weak = 0.342
strain_thresh_init = 0.0
strain_thresh_weak = 1e12
eta_min = 2.5e18
eta_max = 2.5e24
shear_modulus = 0.0
immobile = false
[[rock]] # 2: SP crust (weak decoupling layer)
eta0 = 2e20
activation_energy = 0.0
activation_volume = 0.0
stress_exponent = 1.0
ref_yield_stress = 0.0
reference_density = 3300.0
thermal_expansion = 3e-5
specific_heat = 940.0
thermal_conductivity = 2.96
cohesion_init = 1e15 # plasticity off
cohesion_weak = 1e15
friction_init = 0.0
friction_weak = 0.0
eta_min = 2.5e18
eta_max = 2e20 # stay weak
shear_modulus = 0.0
immobile = false
[[rock]] # 3: overriding plate (no plasticity)
eta0 = 1e16
activation_energy = 540000.0
activation_volume = 1.2
stress_exponent = 3.5
ref_yield_stress = 1e7
reference_density = 3300.0
thermal_expansion = 3e-5
specific_heat = 940.0
thermal_conductivity = 2.96
cohesion_init = 1e15
cohesion_weak = 1e15
friction_init = 0.0
friction_weak = 0.0
eta_min = 2.5e18
eta_max = 2.5e24
shear_modulus = 0.0
immobile = false
[[rock]] # 4: SP mantle (slab core)
eta0 = 1e16
activation_energy = 540000.0
activation_volume = 1.2
stress_exponent = 3.5
ref_yield_stress = 1e7
reference_density = 3300.0
thermal_expansion = 3e-5
specific_heat = 940.0
thermal_conductivity = 2.96
cohesion_init = 5e7
cohesion_weak = 5e7
friction_init = 0.342
friction_weak = 0.342
strain_thresh_init = 0.0
strain_thresh_weak = 1e12
eta_min = 2.5e18
eta_max = 2.5e24
shear_modulus = 0.0
immobile = false
# ══════════════════════════════════════════════════════════════════════════════
# Initial geometry — same generic [[layer]] schema used by Glerum cases
# ══════════════════════════════════════════════════════════════════════════════
[[layer]]
name = "background_mantle"
rock_id = 1
shape = "background"
[[layer]]
name = "overriding_plate"
rock_id = 3
shape = "rectangle"
x_min = 300000.0
x_max = 1800000.0
y_min = 0.0
y_max = 100000.0
[[layer]]
name = "bent_slab_mantle"
rock_id = 4
shape = "bent_band"
thickness = 100000.0
vertices = [[1750000.0, 250000.0], [2000000.0, 0.0]]
[[layer]]
name = "bent_slab_crust"
rock_id = 2
shape = "bent_band"
thickness = 10000.0
vertices = [[1750000.0, 250000.0], [2000000.0, 0.0]]
[[layer]]
name = "sp_crust"
rock_id = 2
shape = "rectangle"
x_min = 2000000.0
x_max = 4700000.0
y_min = 0.0
y_max = 10000.0
[[layer]]
name = "sp_mantle"
rock_id = 4
shape = "rectangle"
x_min = 2000000.0
x_max = 4700000.0
y_min = 10000.0
y_max = 100000.0
# ══════════════════════════════════════════════════════════════════════════════
# Initial temperature — data-driven [[init_T]]
# Holt & Condit: HSC per plate (SP 90 Myr, OP 10 Myr). Adiabat bottom.
# ══════════════════════════════════════════════════════════════════════════════
# Background = constant adiabat T. Then a thin surface BL rectangle
# overlays a linear T(0→1593) over the top 150 km so the lithosphere has
# a proper cold lid. Per-plate HSC entries below override the BL where
# applicable (younger or older HSC profiles). This prevents the
# pathology where a "linear-across-full-depth" background sets the deep
# mantle as colder than the slab (slab would float instead of sinking).
[[init_T]]
name = "background_adiabat"
shape = "background"
profile = "constant"
T_value = 1593.0
[[init_T]]
name = "surface_boundary_layer"
shape = "rectangle"
profile = "linear"
x_min = -1e9
x_max = +1e9
y_min = 0.0
y_max = 150000.0
T_top = 273.0
T_bot = 1593.0
# OP: HSC with 10 Myr age (H-C reference)
[[init_T]]
name = "OP_hsc"
shape = "rectangle"
profile = "hsc"
x_min = 300000.0
x_max = 1800000.0
y_min = 0.0
y_max = 100000.0
T_surf = 273.0
T_mantle = 1573.0
plate_age_yr = 10e6
kappa = 1e-6
# SP: HSC with 90 Myr age (H-C reference)
[[init_T]]
name = "SP_hsc"
shape = "rectangle"
profile = "hsc"
x_min = 2000000.0
x_max = 4700000.0
y_min = 0.0
y_max = 100000.0
T_surf = 273.0
T_mantle = 1573.0
plate_age_yr = 90e6
kappa = 1e-6
# Bent slab: HSC along band (slab core stays cold as it dips)
[[init_T]]
name = "bent_slab_hsc"
shape = "bent_band"
profile = "hsc"
thickness = 100000.0
vertices = [[1750000.0, 250000.0], [2000000.0, 0.0]]
T_surf = 273.0
T_mantle = 1573.0
plate_age_yr = 90e6
kappa = 1e-6