Eclogite foundering (Venus plateau)
A thick basaltic plateau eclogitises through the divariant mode-3 reaction and founders when the mantle beneath it is weak enough (Karagoz 2026, JGR: Planets, in revision). Commented lines switch on composite creep with a fixed grain size, wet olivine and the dry-eclogite rheology blend.
Run it with:
bash
julia -t auto --project=. main.jl cases/venus_eclogite_foundering.tomlFull case file — cases/venus_eclogite_foundering.toml:
toml
# ══════════════════════════════════════════════════════════════════════════════
# venus_eclogite_foundering — eclogite foundering beneath a thick-crust plateau
#
# Generic version of the reference set-up of the Venus eclogite-foundering
# study (stagnant lid, no imposed kinematics, no plume). A 65 km thick
# basaltic plateau crust sits between 20 km normal crust. Only the thick root
# crosses the divariant basalt→eclogite reaction, densifies, and can founder
# as a gravitational drip from beneath the plateau.
#
# Domain 1000 km × 420 km, 5 km cells (200 × 84 elements)
# Gravity Venus, g = 8.87 m/s², surface 735 K, mantle 1700 K
# Geotherm half-space cooling, plate age 150 Myr → Moho (65 km) ≈ 1214 K
# Rheology dislocation creep: dry olivine (mantle), mafic granulite (crust)
# η clamped to [1e19, 1e23] Pa·s
# Reaction divariant (transition = 3) basalt→eclogite, Δρ_max = 550 kg/m³,
# rate-limited (Arrhenius τ_rxn) with reaction-deformation
# softening of weaken_decades = 3 (η ← η·10^(−3·X))
#
# Run: julialauncher --project=. main.jl cases/venus_eclogite_foundering.toml
# The run stops at t_end_myr = 50 Myr, or earlier when the eclogite reaches the
# domain floor (eclogite_floor_stop). A 5 km run takes several hours on a
# laptop; set nx = 100, ny = 42 for a quick 10 km look.
# ══════════════════════════════════════════════════════════════════════════════
[model]
case = "venus_eclogite_foundering" # data-driven: [[layer]] + [[init_T]] below
# Optional granulite→eclogite rheology blend: as the crust eclogitises (X: 0→1)
# its dislocation viscosity is log-interpolated toward this flow law,
# η = exp((1−X)·ln η_granulite + X·ln η_eclogite).
# Off when absent.
# eclogite_flow_law = "Dry_Eclogite-Jin_2001"
eclogite_floor_stop = true # stop when the drip hits the floor (default false)
[mesh]
dim = 2
xsize = 1000000.0 # 1000 km
ysize = 420000.0 # 420 km
nx = 200 # 5 km cells
ny = 84
mnx = 4 # markers per cell (x)
mny = 4 # markers per cell (y)
[gravity]
gy = 8.87 # Venus
[time]
maxxystep = 0.5 # CFL: max marker displacement per step (cells)
maxtmstep = 6.0e13 # dt cap ≈ 1.9 Myr
ntot = 600 # step cap; t_end_myr normally stops the run first
[rheology]
eta_min = 1e19
eta_max = 1e23 # global cap on the cold lid / eclogite root
str_min = 1e-18
eii_seed = 1e-15
[temperature]
T_top = 735.0
T_bot = 1700.0
marker_T_mode = "pic_increment"
do_subgrid_diffusion = true
subgrid_d = 1.0
t_end_myr = 50.0 # self-stop at this model time (Myr)
[physics]
do_stokes = true
do_heat = true
do_advect = true
do_adiabatic = false
do_shear_heating = true
do_friction = true
do_phase_change = true # the basalt→eclogite reaction is the driver
do_melting = false
do_volatile = false
do_free_surface = false
do_radiogenic = false
[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
[amr]
do_amr = false
[fe]
fe_order_v = 2
fe_order_p = 1
fe_order_T = 1
n_picard = 8 # viscoplastic + reaction nonlinearity
picard_tol = 1e-4
[solver]
type = "mumps"
[output]
dir = "output/venus_eclogite_foundering"
every = 5
# ══════════════════════════════════════════════════════════════════════════════
# Rocks
# 1 = asthenosphere (dry olivine)
# 2 = basaltic crust (mafic granulite; eclogitises in place, X carried
# by phase_drho, rock type is NOT switched)
# 3 = lithospheric mantle (dry olivine, stiffer floor eta_min)
# 4 = eclogite (reserved rock slot; unused by transition 3)
#
# Composite creep (optional): add a grain-size-sensitive diffusion law in
# parallel with dislocation creep, 1/η = 1/η_disl + 1/η_diff(d). Grain size is
# constant unless grain_G0 > 0. For rocks 1 and 3:
# diff_prof = "Dry_Olivine_diff-Hirth_Kohlstedt_2003"
# grain_size_mm = 3 # same as grain_init = 3e-3 (m)
# Wet mantle (Hirth & Kohlstedt 2003, C_OH = 1000 H/10^6 Si):
# disl_prof = "Wet_Olivine_disl-Hirth_Kohlstedt_2003"
# diff_prof = "Wet_Olivine_diff-Hirth_Kohlstedt_2003"
# grain_size_mm = 3
# ══════════════════════════════════════════════════════════════════════════════
# Rock 1 — asthenosphere
[[rock]]
disl_prof = "Dry_Olivine-Ranalli_1995"
# diff_prof = "Dry_Olivine_diff-Hirth_Kohlstedt_2003"
# grain_size_mm = 3
reference_density = 3300.0
thermal_expansion = 3e-5
specific_heat = 1250.0
thermal_conductivity = 3.0
cohesion_init = 1e15 # effectively no plasticity
cohesion_weak = 1e15
eta_min = 1e19
eta_max = 1e23
shear_modulus = 0.0
# Rock 2 — basaltic crust (eclogitisable). ρ₀ = 3000; the eclogite density
# gain is carried by phase_drho (full eclogite = 3000 + 550 = 3550 kg/m³).
[[rock]]
disl_prof = "Mafic_Granulite-Ranalli_1995"
reference_density = 3000.0
thermal_expansion = 3e-5
specific_heat = 1250.0
thermal_conductivity = 3.0
cohesion_init = 20e6
cohesion_weak = 20e6
eta_min = 1e19
eta_max = 1e23
shear_modulus = 0.0
# Rock 3 — lithospheric mantle
[[rock]]
disl_prof = "Dry_Olivine-Ranalli_1995"
# diff_prof = "Dry_Olivine_diff-Hirth_Kohlstedt_2003"
# grain_size_mm = 3
reference_density = 3330.0
thermal_expansion = 3e-5
specific_heat = 1250.0
thermal_conductivity = 3.0
cohesion_init = 1e15
cohesion_weak = 1e15
eta_min = 1e21
eta_max = 1e23
shear_modulus = 0.0
# Rock 4 — eclogite (reserved; a transition-1/2 phase could switch into it)
[[rock]]
disl_prof = "Mafic_Granulite-Ranalli_1995"
reference_density = 3500.0
thermal_expansion = 3e-5
specific_heat = 1250.0
thermal_conductivity = 3.0
cohesion_init = 1e15
cohesion_weak = 1e15
eta_min = 1e21
eta_max = 1e23
shear_modulus = 0.0
# ══════════════════════════════════════════════════════════════════════════════
# Basalt→eclogite: divariant equilibrium (transition = 3), rate-limited.
# X_eq(P,T) = clamp((P − P_c(T))/W + ½, 0, 1), P_c(T) = P_ref + S·(T − T_ref)
# dX/dt toward X_eq with τ_rxn(T) = react_tau_ref·exp((react_Q/R)(1/T − 1/react_Tref))
# phase_drho = drho·X → eclogite fraction X = phase_drho/drho
# reaction softening: η ← max(η·10^(−weaken_decades·X), eta_min)
# Set react_tau_ref = 0 for instantaneous equilibrium, weaken_decades = 0 to
# switch the softening off.
# ══════════════════════════════════════════════════════════════════════════════
[[phase]]
rock_from = 2
rock_to = 2 # transition 3 does not switch rock type
transition = 3
drho = 550.0 # Δρ_max basalt → full eclogite (kg/m³)
band_width = 6.0e8 # W (Pa): divariant width (~20 km)
P_ref = 1.010517e9 # P_c at T_ref (Pa)
clapeyron = 1.520571e6 # S = dP/dT (Pa/K)
T_ref = 928.15 # K
latent = 0.0
react_tau_ref = 3.1557e13 # ~1 Myr at react_Tref
react_Q = 3.0e5 # J/mol
react_Tref = 973.15 # K (700 °C)
weaken_decades = 3.0
# ══════════════════════════════════════════════════════════════════════════════
# Initial temperature: half-space cooling of a 150 Myr plate (Moho ≈ 1214 K).
# Moho temperature is the main control on foundering; vary plate_age_yr.
# ══════════════════════════════════════════════════════════════════════════════
[[init_T]]
name = "hsc_geotherm"
shape = "background"
profile = "hsc"
T_surf = 735.0
T_mantle = 1700.0
plate_age_yr = 150.0e6
kappa = 1e-6
# ══════════════════════════════════════════════════════════════════════════════
# Initial geometry: flat layers; plateau crust between x = 250 and 650 km.
# ══════════════════════════════════════════════════════════════════════════════
[[layer]]
name = "asthenosphere"
rock_id = 1
shape = "background"
[[layer]]
name = "lithospheric_mantle"
rock_id = 3
shape = "rectangle"
x_min = 0.0
x_max = 1000000.0
y_min = 0.0
y_max = 110000.0
[[layer]]
name = "normal_crust"
rock_id = 2
shape = "rectangle"
x_min = 0.0
x_max = 1000000.0
y_min = 0.0
y_max = 20000.0
[[layer]]
name = "plateau_thick_crust"
rock_id = 2
shape = "rectangle"
x_min = 250000.0
x_max = 650000.0
y_min = 0.0
y_max = 65000.0