Mid-ocean ridge spreading

Symmetric spreading with passive asthenospheric upwelling at the axis. Starter template.

Run it with:

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

Full case file — cases/examples/ridge.toml:

# ══════════════════════════════════════════════════════════════════════════════
# EXAMPLE TEMPLATE — Mid-Ocean Ridge spreading (divergent margin)
#
# Two oceanic plates pulled apart at a constant half-rate; passive asthenospheric
# upwelling fills the gap beneath the ridge axis (x = xsize/2). A hot Gaussian
# anomaly at the axis seeds the upwelling.
#
# This is a STARTER TEMPLATE built on the data-driven [[layer]]/[[init_T]] path
# (the [model] case name is unknown to select_case, so geometry comes from the
# layer/temperature tables below). Tune resolution, rates, and rheology before
# using it for production. Coordinate convention: y = depth, positive DOWN.
# ══════════════════════════════════════════════════════════════════════════════

[model]
case = "ridge"                     # unknown name → data-driven layers/init_T

[mesh]
dim   = 2
xsize = 600000.0                   # 600 km wide
ysize = 200000.0                   # 200 km deep (incl. 20 km sticky air)
nx    = 201
ny    = 67
mnx   = 4
mny   = 4

[gravity]
gy = 9.81

[time]
maxxystep = 0.5
maxtmstep = 3.1556e12              # ~100 kyr cap
ntot      = 200

[rheology]
eta_min = 1e18
eta_max = 1e24

[temperature]
T_top = 273.0
T_bot = 1623.0

[physics]
do_stokes        = true
do_heat          = true
do_advect        = true
do_shear_heating = true
do_friction      = true
do_free_surface  = true
do_radiogenic    = false

# Divergent kinematics: both side walls move outward at 3 cm/yr (half-rate),
# the bottom is left open (traction-free) so asthenosphere flows up to fill in.
[boundary_velocity]
bc_left    = "uniform"
bc_right   = "uniform"
bc_surface = "free_slip"
bc_deep    = "traction"            # open bottom: passive inflow

[boundary_velocity_values]
bc_left_vx  = -9.51173e-10         # −3 cm/yr (leftward)
bc_right_vx =  9.51173e-10         # +3 cm/yr (rightward)

[free_surface]
do_ale_free_surface = false
air_rock_id = 1

[fe]
n_picard   = 3
picard_tol = 1e-3

[solver]
type = "mumps"

[output]
dir   = "output/ridge"
every = 10

# ── Rocks: 1=air, 2=oceanic crust (basalt), 3=lithospheric mantle, 4=asthenosphere
[[rock]]                           # 1 — sticky air
nu0=1e18
ro0=1.0
alpha=0.0
cp=1000.0
kt=50.0
ht=0.0
mm=1.0
dh=0.0
dv=0.0
a0=0.0
b0=0.0
a1=0.0
b1=0.0
nu_min=1e18
nu_max=1e19
immobile=false
G=0.0

[[rock]]                           # 2 — oceanic crust (basalt/gabbro)
nu0=1e21
ro0=2900.0
alpha=3e-5
cp=1000.0
kt=2.5
ht=2e-7
mm=1.0
dh=0.0
dv=0.0
a0=1e7
b0=0.2
a1=1e6
b1=0.02
e0=0.0
e1=0.5
nu_min=1e19
nu_max=1e24
immobile=false
G=3e10

[[rock]]                           # 3 — lithospheric mantle (dry olivine)
nu0=1.97e17
ro0=3300.0
alpha=3e-5
beta=1e-11
cp=1000.0
kt=3.0
ht=2.5e-8
mm=1.0
dh=154000.0
dv=0.3
a0=3e7
b0=0.1
a1=3e6
b1=0.0
e0=0.0
e1=0.2
nu_min=1e19
nu_max=1e24
immobile=false
G=7e10

[[rock]]                           # 4 — asthenosphere
nu0=1.97e17
ro0=3300.0
alpha=3e-5
beta=1e-11
cp=1000.0
kt=3.0
ht=2.5e-8
mm=2.3
dh=154000.0
dv=0.3
ss=3e4
a0=1e6
b0=0.03
a1=1e5
b1=0.01
e0=0.0
e1=0.1
nu_min=1e18
nu_max=1e23
immobile=false
G=7e10

# ── Initial geometry (evaluated in order; later layers override earlier) ──────
[[layer]]
name    = "asthenosphere"
rock_id = 4
shape   = "background"

[[layer]]
name    = "lithospheric_mantle"
rock_id = 3
shape   = "rectangle"
x_min   = 0.0
x_max   = 600000.0
y_min   = 27000.0                  # 7 km below the crust base
y_max   = 60000.0

[[layer]]
name    = "oceanic_crust"
rock_id = 2
shape   = "rectangle"
x_min   = 0.0
x_max   = 600000.0
y_min   = 20000.0                  # top of crust = base of air
y_max   = 27000.0                  # 7 km thick

[[layer]]
name    = "air"
rock_id = 1
shape   = "rectangle"
x_min   = 0.0
x_max   = 600000.0
y_min   = 0.0
y_max   = 20000.0

# ── Initial temperature: linear geotherm in the mantle, cold air, hot axis ────
[[init_T]]
name    = "geotherm"
shape   = "background"
profile = "linear"
T_top   = 273.0
T_bot   = 1623.0

[[init_T]]
name    = "air"
shape   = "rectangle"
profile = "constant"
x_min   = 0.0
x_max   = 600000.0
y_min   = 0.0
y_max   = 20000.0
T_value = 273.0

[[init_T]]
name    = "ridge_upwelling"
shape   = "gaussian"
profile = "gaussian"
xc      = 300000.0                 # ridge axis (domain centre)
yc      = 60000.0
r0      = 40000.0
dT      = 250.0                    # +250 K hot anomaly drives upwelling
T_base  = 1623.0