MILETModeling in Lithospheric Evolution Tectonics

A thermomechanical finite-element framework for lithosphere and mantle dynamics, developed entirely in Julia on top of Gridap and named after Miletus, a historic center of early scientific thought.

Getting started

Governing equations

Parameter reference

M1 subduction benchmark — evolving thermo-mechanical model

Stokes + energy with markers-in-cell

Inf-sup stable Taylor–Hood (Q2/Q1) Stokes coupled to heat transport; RK4 marker advection carries composition, temperature, stress, plastic strain, grain size and porosity.

Visco-elasto-plastic rheology

Composite diffusion + dislocation creep with experimentally calibrated flow laws, Drucker–Prager yield, strain softening, Arrhenius healing, and reaction softening.

Geometric multigrid in 2-D and 3-D

FGMRES on the full saddle-point system, preconditioned by Galerkin-coarsened geometric multigrid V-cycles with Chebyshev smoothing — variable viscosity inherited exactly, no rediscretisation.

Two-phase McKenzie flow

Porous melt migration through a compacting viscous matrix — a 3-field (u, p_f, p_c) extension of the Stokes solve that reduces exactly to single-phase Stokes at zero porosity.

Grain-size evolution

Arrhenius normal grain growth balanced against Austin–Evans (2007) paleowattmeter reduction, two-way coupled to the diffusion/dislocation mechanism partitioning.

Phase change, melting, free surface

Divariant rate-limited metamorphic reactions with reaction-progress softening, melting parameterisations, sticky-air and ALE free surfaces, AMR, checkpoint/restart.

What is MILET?

Modeling in Lithospheric Evolution Tectonics is a thermo-mechanical finite-element code for lithospheric and mantle dynamics. It solves the coupled Stokes and energy equations with a mixed, inf-sup stable Taylor–Hood discretisation ( $Q_{2}$ velocity, $Q_{1}$ pressure, built on Gridap.jl) combined with a marker-in-cell representation of material. It is designed for problems in which composition, rheology, and metamorphic phase change evolve together — lithospheric foundering, delamination, drips, and subduction.

MILET is a flat-include research code: main.jl includes the sources under src/ in dependency order and drives one simulation per run, fully specified by a single TOML parameter file (deep-merged onto cases/_defaults.toml).

Capabilities

Incompressible Boussinesq Stokes flow (anelastic ALA/TALA forms scaffolded), 2-D and 3-D Cartesian plus spherical-chunk geometry — Governing equations
Visco-elasto-plastic composite rheology with calibrated dislocation-creep flow laws (Ranalli 1995 catalogue), grain-size-sensitive diffusion creep, and Drucker–Prager yield — Rheology
Marker-in-cell transport: RK4 advection, reseeding, and subgrid-scale temperature diffusion after Gerya & Yuen (2003) — Markers
Divariant, rate-limited (Arrhenius) phase transitions with reaction-progress-dependent softening — Phase changes
Melting parameterisations and melt weakening — Melting
Heat transport with shear, adiabatic, and radiogenic heating, plus SUPG / entropy-viscosity stabilisation options — Heat transport
Sticky-air free surface with FSSA stabilisation and an ALE surface-tracking option — Free surface
Direct (MUMPS) and iterative (GMG, PETSc) solvers with cached factorisations and warm starts — Linear solvers
Picard and Newton nonlinear loops, CFL-limited time stepping, checkpoint/restart, VTK output — Time stepping
Adaptive mesh refinement and ALE — AMR & ALE

New capabilities

Geometric multigrid 3-D Stokes solver — Geometric multigrid. Right-preconditioned FGMRES on the monolithic saddle-point system; the velocity block is preconditioned by V-cycles of a Galerkin geometric multigrid ( $A_{ℓ + 1} = P_{ℓ}^{T} A_{ℓ} P_{ℓ}$ , so variable viscosity and Dirichlet elimination are inherited exactly from the fine operator), with Chebyshev smoothing and a cached-LU coarse solve; the pressure Schur complement uses the $η$ -weighted lumped pressure mass (Elman–Silvester–Wathen). Select with [solver] type = "gmg" (gmg_levels, gmg_cheb_degree, gmg_cycles, gmg_coarse_dofs); shipped 3-D cases: cases/plume3d_cube_gmg.toml, cases/plume3d_slab_gmg.toml. Implementation: src/solvers/GMG.jl, src/solvers/BlockMG.jl.

Two-phase (McKenzie) flow — Two-phase flow. Porous melt migration through a compacting viscous matrix, as a 3-field extension of the Stokes solve (Katz 2008; Rhebergen, Wells & Katz 2014):

- \nabla \cdot (2 η_{ϕ} \dot{ε} (u_{s})) + \nabla p_{f} + \nabla p_{c} = \bar{ρ} g, \nabla \cdot u_{s} + \frac{p_{c}}{ζ_{ϕ}} = 0, - \nabla \cdot (\frac{k_{ϕ}}{μ_{f}} (\nabla p_{f} - ρ_{f} g)) - \nabla \cdot u_{s} = 0,

discretised as $Q_{2}$ solid velocity / $Q_{1}$ fluid pressure / $Q_{1}$ compaction pressure, with porosity carried on markers. At $ϕ = 0$ the system reduces exactly to single-phase Stokes (a regression test enforces this). Enabled with [twophase] do_two_phase = true (k0, phi0, n_perm, mu_fluid, rho_fluid, alpha_phi, phi_min, pf_surface0, …). Implementation: src/physics/TwoPhase.jl.

Grain-size evolution (paleowattmeter) — Grain-size evolution. Per-marker grain size $d$ evolves by Arrhenius normal grain growth balanced against Austin & Evans (2007) reduction driven by the dislocation-creep work rate $Ψ_{disl} = σ_{I I} {\dot{ε}}_{disl}$ :

\frac{d d}{d t} = \frac{G_{0}}{p d^{p - 1}} \exp (- \frac{E_{g} + P V_{g}}{R T}) - \frac{λ}{c γ} d^{2} Ψ_{disl},

integrated with exact operator-split sub-flows (update_grain_size, src/materials/GrainSize.jl) and two-way coupled to the composite diffusion + dislocation rheology. Per-rock keys grain_G0, grain_p, grain_Eg, grain_Vg, grain_lambda, grain_gamma, grain_c, grain_init (grain_G0 = 0 freezes the grain size).

Strain softening & healing — Softening & healing. Linear cohesion/friction weakening over [strain_thresh_init, strain_thresh_weak], now with plastic-strain gating ([rheology] strain_gate_plastic = true: weakening strain accumulates only where the cell actually yields), temperature-dependent Arrhenius healing ( $τ_{heal} = τ_{0} \exp [\frac{E}{R} (\frac{1}{T} - \frac{1}{T_{ref}})]$ ; per-rock heal_tau0, heal_E, heal_Tref), optional viscous strain softening (visc_soft_factor), and a configurable yield-stress cap (yield_stress_max).

Quick start

bash

# install dependencies (once)
julia --project=. -e 'using Pkg; Pkg.instantiate()'

# verification benchmark: SolCx (variable-viscosity Stokes, analytic solution)
julia --project=. main.jl cases/solcx_fixed.toml

# convection benchmark: Blankenbach case 1a
julia --project=. main.jl cases/blankenbach_1a.toml

# a worked geodynamic example: thermo-mechanical subduction
julia --project=. main.jl cases/subduction.toml

# 3-D plume with the geometric-multigrid solver
julia -t auto --project=. main.jl cases/plume3d_cube_gmg.toml

Results are written as VTK files (solution fields, markers, topography) to output/<case>/, viewable in ParaView — see Getting started for output formats and the depth-positive-down camera setup. Every TOML key is documented in the parameter-file reference.

Verification

SolCx (Duretz et al. 2011): grid convergence of the variable-viscosity Stokes solution (cases/solcx_fixed_n*.toml)
Blankenbach et al. (1989) case 1a: steady convection diagnostics ( $N u$ , $V_{rms}$ ) at three resolutions (cases/blankenbach_*.toml)
Smoke and benchmark tests under test/ (julia --project=. test/runtests.jl) — see Benchmarks and Testing

Notes

The multi-species volatile transport module is part of a separate study in preparation and is not included in this release; shipped cases run with volatiles disabled.
3-D and spherical-chunk support are functional but less exercised than the 2-D Cartesian path.

Citing

A methods paper describing MILET and its application to eclogite-driven lithospheric foundering on Venus is under review at Earth and Planetary Science Letters. Until it is published, please cite the repository: github.com/okaragoz/MILET.