Mantle plumes in a spherical chunk (2-D)

Whole-mantle plumes in a 92° annular wedge with radial gravity: a hot core-mantle boundary layer, thickened at a few uneven angles, feeds plumes that rise toward a cold surface boundary layer.

Run it with:

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

Full case file — cases/chunk_plumes.toml:

# ══════════════════════════════════════════════════════════════════════════════
# Rising mantle plumes in a 2D spherical chunk (annular wedge, radial gravity).
# A hot CMB boundary layer, thickened at a few uneven angles, feeds plumes that
# rise through a well-mixed mantle toward a cold surface boundary layer.
# Whole-mantle depth (2890 km), 92° wedge, isoviscous η = 5e22 Pa s (Ra ≈ 3e6).
# ══════════════════════════════════════════════════════════════════════════════

[model]
case = "chunk_plumes"   # chunk IC with seeded CMB plumes (initial_temperature_chunk_plumes)

dim = 2

# Chunk dimensions: full-mantle-thickness wedge spanning 1.6 rad (about 92°).
geometry            = "chunk"
r_inner             = 3.481e6   # core-mantle boundary (m)
r_outer             = 6.371e6   # surface (m)
chunk_theta_extent  = 1.6       # rad

# nx = number of radial gridlines, ny = angular gridlines.
nx    = 49
ny    = 121
mnx   = 2
mny   = 2

# Radial gravity instead of Cartesian.
gravity_model = "radial"
g_mag         = 9.81
gx            = 0.0
gy            = 0.0
pinit         = 0.0

# Free-slip on all walls.
bc_left    = "free_slip"
bc_right   = "free_slip"
bc_surface = "free_slip"
bc_deep    = "free_slip"

eta_min = 1e18
eta_max = 1e24

maxxystep = 0.5
maxtkstep = 50.0
maxtmstep = 3.1556e13   # 1 Myr cap; the Courant limit sets dt once plumes rise
ntot      = 300

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

fe_order_v = 2
fe_order_p = 1
fe_order_T = 1
n_picard   = 1
picard_tol = 1e-6

nonlinear_method = "picard"
newton_damping   = 1.0

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

output_dir   = "output/chunk_plumes"
output_every = 5
str_min      = 1e-18

T_top         = 273.0
T_bot         = 3500.0     # CMB temperature
T_ref_density = 273.0      # Boussinesq reference
buoyancy_from_T = true     # buoyancy from the FE temperature (single composition)

# placeholder values for unused subduction fields
trench_x = 0.0
plate_thickness = 0.0
plate_age = 1.0
continent_thickness = 0.0
air_thick = 0.0
x_plate_left = 0.0
x_plate_right = 0.0
vx_plate = 0.0
vx_plate_right = 0.0
air_rock_id = 1
surface_smooth = 0.0

[[rock]]
# Constant isoviscous mantle (dh=0, ss=0, n=1), η = 5e22 Pa·s (Ra ≈ 3e6, so
# the thermal boundary layers span a few cells of this mesh).
nu0      = 5e22
dh       = 0.0
dv       = 0.0
ss       = 0.0
mm       = 1.0
ro0      = 4500.0     # mean mantle ρ
alpha    = 3e-5
beta     = 0.0
cp       = 1250.0
kt       = 4.0
ht       = 0.0
a0       = 0.0
b0       = 0.0
a1       = 0.0
b1       = 0.0
e0       = 0.0
e1       = 0.0
nu_min   = 5e22
nu_max   = 5e22
immobile = false
G        = 7e10