Stagnant lid

Temperature-dependent-viscosity convection with an immobile lid (mobility M ≪ 1).

Run it with:

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

Full case file — cases/solomatov_stagnant.toml:

# ══════════════════════════════════════════════════════════════════════════════
# Solomatov (1995) — temperature-dependent-viscosity convection, stagnant-lid.
# Frank-Kamenetskii viscosity η = η₀·exp(−θ·T*), T* ∈ [0,1] (cold/stiff top,
# hot/weak base); viscosity contrast Δη = exp(θ). Implemented via the engine's
# FK branch (activation_energy = −θ in nondimensional T).
#
# Baseline: θ = 9.21 (Δη = 1e4), Ra_surf (cold ref) = 100 → Ra_interior ≈ 1e6.
# Expected: thick stagnant lid, surface mobility M ≪ 1. This is the LAB
# substrate — vary θ and Ra to map stagnant ↔ episodic ↔ mobile lid.
# SECTIONED layout (required — see blankenbach_1a.toml note).
# ══════════════════════════════════════════════════════════════════════════════

[model]
case = "blankenbach"        # reuse the box-convection rock/T closures

[mesh]
dim   = 2
xsize = 1.0
ysize = 1.0
nx    = 49
ny    = 49
mnx   = 4
mny   = 4

[gravity]
gx = 0.0
gy = 1.0

[time]
maxxystep = 0.5
maxtkstep = 0.005
maxtmstep = 0.02
ntot      = 2000

[rheology]
eta_min  = 1e-6
eta_max  = 1e2
str_min  = 1e-18
eii_seed = 0.0

[temperature]
T_top         = 0.0
T_bot         = 1.0
T_ref_density = 0.0
rho_ref_heat  = 1.0
marker_T_mode = "pic_increment"

[physics]
do_stokes        = true
do_heat          = true
do_advect        = true
do_adiabatic     = false
do_shear_heating = false
do_friction      = false
do_phase_change  = false
do_melting       = false
do_free_surface  = false
do_radiogenic    = false
buoyancy_from_T  = true

[boundary_velocity]
bc_left    = "free_slip"
bc_right   = "free_slip"
bc_surface = "free_slip"
bc_deep    = "free_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   = 3
picard_tol = 1e-4

[solver]
type = "mumps"

[output]
dir   = "output/solomatov_stagnant"
every = 50

# Frank-Kamenetskii temperature-dependent viscosity.
#   nu0 = η₀ = cold (T*=0) reference viscosity = 1
#   dh  = −θ  (negative → FK branch; θ = 9.21 → Δη = exp(9.21) = 1e4)
#   dv  = FK reference temperature T*_ref = 0
#   alpha = Ra on the cold reference viscosity (Ra_interior = alpha·e^θ ≈ 1e6)
[[rock]]
nu0      = 1.0
dh       = -9.21
dv       = 0.0
ss       = 0.0
mm       = 1.0
ro0      = 1.0
alpha    = 100.0
beta     = 0.0
cp       = 1.0
kt       = 1.0
ht       = 0.0
a0       = 0.0
b0       = 0.0
a1       = 0.0
b1       = 0.0
e0       = 0.0
e1       = 1.0
nu_min   = 1e-5
nu_max   = 10.0
immobile = false
G        = 0.0