Rheology

MILET's constitutive model is a pointwise visco-elasto-plastic (VEP) composite evaluated per marker at $(T,P,{\dot{\epsilon }}_{II},{\epsilon }_{\mathrm{pl}},d)$ and harmonically averaged to cells. The core lives in src/materials/Rheology.jl (_compute_viscosity_core, compute_VEP_viscosity, viscous_mechanism_partition) with calibrated flow laws in src/materials/FlowLaws.jl.

The total deviatoric strain rate decomposes as a serial (Maxwell-type) arrangement sharing one deviatoric stress:

$$\dot{\epsilon }={\dot{\epsilon }}_{\mathrm{visc}}+{\dot{\epsilon }}_{\mathrm{plast}}+{\dot{\epsilon }}_{\mathrm{elast}}.$$

1. Viscous branch

Calibrated dislocation creep (creep_A / disl_prof)

The secondary-invariant power law (Ranalli 1995 convention)

$${\dot{\epsilon }}_{II}=A{\sigma }_{II}^n\exp (-\frac{E+PV}{RT}),A{\text{in MPa}}^{-n}{\mathrm{s}}^{-1}$$

is inverted to an effective viscosity (dislocation_viscosity in FlowLaws.jl):

$${\eta }_{\mathrm{disl}}={10}^6\cdot \frac{1}{2}{A}^{-1/n}{\dot{\epsilon }}_{II}^{(1-n)/n}\exp (\frac{E+PV}{nRT})[\mathrm{Pa}\mathrm{s}],$$

with $E$ in J/mol, $V$ in m³/mol, $P$ in Pa; the ${10}^6$ converts the MPa-based stress convention. A rock activates this path by setting creep_A > 0 explicitly or by naming a catalog profile:

[[rock]]
disl_prof = "Dry_Olivine-Ranalli_1995"

lookup_flow_law resolves the name from the FLOW_LAWS catalog and overrides activation_energy, activation_volume, stress_exponent, and creep_A. Shipped entries:

Profile	$A$ (MPa⁻ⁿ s⁻¹)	$n$	$E$ (kJ/mol)	$V$ (m³/mol)
Dry_Olivine-Ranalli_1995	2.5e4	3.5	532	0
Wet_Olivine-Ranalli_1995	2.0e3	4.0	471	0
Mafic_Granulite-Ranalli_1995	1.4e4	4.2	445	0
Plagioclase_An75-Ranalli_1995	3.3e-4	3.2	238	0
Quartzite-Ranalli_1995	6.7e-6	2.4	156	0
Wet_Quarzite-Ranalli_1995	3.2e-4	2.3	154	0
Dry_Olivine_disl-Hirth_Kohlstedt_2003	1.1e5	3.5	530	14e-6

Grain-size-sensitive diffusion creep (diff_A / diff_prof)

A parallel linear-viscous mechanism (Hirth & Kohlstedt 2003 form):

$${\dot{\epsilon }}_{II}=A{\sigma }_{II}{d}^{-m}\exp (-\frac{E+PV}{RT})⟹{\eta }_{\mathrm{diff}}={10}^6\frac{(d\cdot {10}^6{)}^m}{2A}\exp (\frac{E+PV}{RT}),$$

with lab units $A$ in MPa⁻¹ µm${}^m$ s⁻¹ and grain size $d$ in metres at the call site (diffusion_viscosity). Configure with explicit diff_A, diff_m, diff_E, diff_V keys or a catalog profile (DIFF_CREEP_LAWS): Dry_Olivine_diff-Hirth_Kohlstedt_2003 ($A=1.5\times {10}^9$, $m=3$, $E=375$ kJ/mol, $V=6\times {10}^{-6}$) or Wet_Olivine_diff-Hirth_Kohlstedt_2003 ($A={10}^9$, $m=3$, $E=335$ kJ/mol, $V=4\times {10}^{-6}$). diff_A = 0 (default) disables the mechanism. The grain size $d$ is a marker property evolved by the paleowattmeter — see Grain-size evolution.

Harmonic composite

Both mechanisms act in parallel at common stress, so compliances add (the standard LaMEM/ASPECT "harmonic" approximation, evaluated at the total strain-rate invariant):

$$\frac{1}{{\eta }_{\mathrm{visc}}}=\frac{1}{{\eta }_{\mathrm{disl}}}+\frac{1}{{\eta }_{\mathrm{diff}}(d)}.$$

Legacy Arrhenius path (eta0)

When creep_A = 0 and eta0 > 0, the pre-exponential viscosity form is used:

$${\eta }_{\mathrm{creep}}={\eta }_0\exp (\frac{Q+V{P}_{\mathrm{bar}}}{nRT}),$$

with $n = $ stress_exponent ($n=1$ → diffusion creep, $n>1$ → dislocation creep; the $1/n$ in the exponent is essential to avoid $\sim {10}^8\times$ overestimated viscosity at depth). Unit note for this path: activation_energy $Q$ is in J/mol but activation_volume $V$ is in J/(mol·bar) since $P$ enters in bar. Strain-rate dependence is supplied by the power-law blend when ref_yield_stress ${\sigma }_{\mathrm{ref}}>0$ and $n>1$:

$$x=(\frac{2{\eta }_{\mathrm{creep}}{\dot{\epsilon }}_{II}}{{\sigma }_{\mathrm{ref}}}{)}^{1-1/n},{\eta }_{\mathrm{visc}}=\frac{{\eta }_{\mathrm{creep}}}{1+x},$$

which interpolates continuously between the linear regime (${\dot{\epsilon }}_{II}\to 0$) and the power-law regime without an explicit switch. A legacy "air-like" branch (activation_energy < 0) reproduces the I3ELVIS depth-linear softening, and a legacy grain_exponent modifier scales eta0-path viscosity by $(d/d_{\mathrm{ref}}{)}^m$.

Viscous strain softening (visc_soft_factor)

Fabric/CPO weakening of ductile shear zones: the aggregate viscous viscosity is divided by a factor that ramps linearly from $1$ to visc_soft_factor over the same accumulated-strain interval $[{\epsilon }_0,{\epsilon }_f]$ = [strain_thresh_init, strain_thresh_weak] as the plastic weakening. visc_soft_factor = 1.0 (default) disables it; e.g. 0.3 gives a 3.3× weaker steady shear zone.

2. Plastic branch (Drucker–Prager + strain weakening)

Drucker–Prager yield with linear strain weakening of both cohesion and friction:

$${\tau }_{\mathrm{yield}}=min(C({\epsilon }_{\mathrm{pl}})+\sin \phi ({\epsilon }_{\mathrm{pl}})\cdot P,{\sigma }_{y,max}),$$

where $C$ and $\sin \phi$ interpolate linearly from (cohesion_init, friction_init) to (cohesion_weak, friction_weak) as the accumulated plastic strain ${\epsilon }_{\mathrm{pl}}$ grows across [strain_thresh_init, strain_thresh_weak]. Note friction_init / friction_weak are $\sin \phi$ directly (dimensionless), not angles.

The absolute yield cap yield_stress_max (default 5e8 Pa, after Kaus 2010 / Glerum et al. 2018) prevents $C+\sin \phi P$ from growing to many GPa at depth — without it cold slab cores accumulate unbounded elastic stress and trigger spurious break-off cascades; set <= 0 to disable.

Plasticity is realized rate-dependently as an effective viscosity

$${\eta }_{\mathrm{plast}}=\frac{{\tau }_{\mathrm{yield}}}{2{\dot{\epsilon }}_{II}},$$

so that $2{\eta }_{\mathrm{plast}}{\dot{\epsilon }}_{II}={\tau }_{\mathrm{yield}}$ exactly when plasticity is engaged.

Visco-plastic combination

The combination is stress-limiting, not harmonic — the weaker mechanism governs:

$$\frac{1}{{\eta }_{\mathrm{VP}}}=max(\frac{1}{{\eta }_{\mathrm{visc}}},\frac{1}{{\eta }_{\mathrm{plast}}})⟺{\eta }_{\mathrm{VP}}=min({\eta }_{\mathrm{visc}},{\eta }_{\mathrm{plast}}),$$

which caps the stress at yield once plasticity engages rather than summing compliances. _compute_viscosity_core also returns a yielded flag (inv_plast > inv_visc); its marker average forms the per-cell yield fraction state.yield_frac_c, which gates plastic-strain accumulation on markers (update_marker_strain! in src/particles/Projection.jl, toggled by strain_gate_plastic) so weakening strain builds only in shear zones that actually yield. Strain healing (heal_tau0, heal_E, heal_Tref) anneals ${\epsilon }_{\mathrm{pl}}$ back — see Softening & healing.

3. Maxwell elasticity {#maxwell-elasticity} {#3.-Maxwell-elasticity-{#maxwell-elasticity}}

A Maxwell body $\dot{\sigma }/G+\sigma /{\eta }_{\mathrm{VP}}=\dot{\epsilon }$ discretized with backward Euler over $\mathrm{\Delta }t$ gives (compute_VEP_viscosity):

$$Z=\frac{G\mathrm{\Delta }t}{{\eta }_{\mathrm{VP}}+G\mathrm{\Delta }t}\in [0,1],{\eta }_{\mathrm{VEP}}=Z{\eta }_{\mathrm{VP}},\chi =1-Z,$$$${\mathit{\sigma }}_{\mathrm{new}}=2{\eta }_{\mathrm{VEP}}\dot{\epsilon }+\chi {\mathit{\sigma }}_{\mathrm{old}}.$$

Limits: $G\mathrm{\Delta }t\gg {\eta }_{\mathrm{VP}}$ → $Z\to 1$, $\chi \to 0$ (pure viscous, old stress forgotten); $G\mathrm{\Delta }t\ll {\eta }_{\mathrm{VP}}$ → $Z\to 0$, ${\eta }_{\mathrm{VEP}}\to 0$, ${\mathit{\sigma }}_{\mathrm{new}}\to {\mathit{\sigma }}_{\mathrm{old}}$ (pure elastic). In between, stress relaxes over the Maxwell time ${\tau }_M=\eta /G$. The per-rock key is shear_modulus (Pa; typical silicate $\approx 7\times {10}^{10}$); shear_modulus <= 0 disables elasticity ($Z=1$, $\chi =0$). The $\chi {\mathit{\sigma }}_{\mathrm{old}}$ term feeds the Stokes RHS — see Governing equations.

4. Melt weakening

_compute_viscosity_core applies an exponential melt weakening to the total visco-plastic viscosity (Mei et al. 2002):

$$\eta \leftarrow \eta e^{-{\varphi }_w{F}_{\mathrm{melt}}},$$

via the melt_factor (${\varphi }_w$) and F_melt (marker melt fraction) keyword arguments — melt weakens both diffusion and dislocation creep, hence the application to the aggregate. The single-phase marker projection currently calls it with melt_factor = 0.0; active matrix melt weakening in the shipped physics goes through the two-phase module, which applies $\eta e^{-{\alpha }_{\varphi }\varphi }$ with the [twophase] key alpha_phi (default 27.0) in src/physics/TwoPhase.jl. See Melting and Two-phase flow.

5. Clamps and evaluation pipeline

The final viscosity is clamped to

$$\eta \in [max(\mathtt{\text{eta\_min}},\mathtt{\text{rock.eta\_min}}),min(\mathtt{\text{eta\_max}},\mathtt{\text{rock.eta\_max}})],$$

combining the global [rheology] bounds (defaults 1e19 / 1e25 Pa·s) with the per-rock ones.

project_markers_to_cells! (src/particles/Projection.jl) evaluates compute_VEP_viscosity per marker with marker temperature, lithostatic pressure $P={\rho }_{0}gy$ (converted to bar), the host cell's ${\dot{\epsilon }}_{II}$ from the previous iteration, the marker's accumulated strain, grain size, and melt fraction. Cell viscosity and thermal conductivity are harmonic marker means; $\rho$, $c_p$, $h_r$, $\chi$, ${\mathit{\sigma }}_{\mathrm{old}}$, and the yield fraction are arithmetic means.

Mechanism partition (wattmeter input)

viscous_mechanism_partition splits the viscous deformation at a point for the grain-size wattmeter. At the harmonic-composite approximation both mechanisms see the common stress ${\sigma }_{II}=2{\eta }_{\mathrm{visc}}{\dot{\epsilon }}_{II}$ and partition as ${\dot{\epsilon }}_{\mathrm{mech}}={\sigma }_{II}/(2{\eta }_{\mathrm{mech}})$. It returns $({\eta }_{\mathrm{visc}},{\eta }_{\mathrm{diff}},{\eta }_{\mathrm{disl}},{\dot{\epsilon }}_{\mathrm{diff}},{\dot{\epsilon }}_{\mathrm{disl}},{\mathrm{\Psi }}_{\mathrm{disl}})$ with the dislocation work rate ${\mathrm{\Psi }}_{\mathrm{disl}}={\sigma }_{II}{\dot{\epsilon }}_{\mathrm{disl}}$. Deliberately unclamped: applying the eta_min/eta_max numerical guards here would break the identity ${\dot{\epsilon }}_{\mathrm{diff}}+{\dot{\epsilon }}_{\mathrm{disl}}={\dot{\epsilon }}_{II}$ and inflate the wattmeter work rate. With diff_A = 0 the full work rate is attributed to dislocation creep (the conservative choice for grain-size reduction). See Grain-size evolution.

Key reference

Key (per [[rock]] unless noted)	Units	Meaning
eta0	Pa·s	legacy pre-exponential viscosity
activation_energy, activation_volume, stress_exponent	J/mol; J/(mol·bar) legacy / m³/mol calibrated; —	$Q$/$E$, $V$, $n$
ref_yield_stress	Pa	${\sigma }_{\mathrm{ref}}$ of the $(1+x)$ power-law blend
creep_A or disl_prof	MPa⁻ⁿ s⁻¹ / name	calibrated dislocation creep
diff_A, diff_m, diff_E, diff_V or diff_prof	lab units / name	grain-size-sensitive diffusion creep
cohesion_init, cohesion_weak	Pa	$C_0$, $C_f$
friction_init, friction_weak	—	$\sin {\phi }_0$, $\sin {\phi }_f$
strain_thresh_init, strain_thresh_weak	—	weakening strain bracket $[{\epsilon }_0,{\epsilon }_f]$
visc_soft_factor	—	viscous softening end-member (1 = off)
yield_stress_max	Pa	absolute yield cap (default 5e8; ≤ 0 off)
shear_modulus	Pa	Maxwell $G$ (0 = no elasticity)
eta_min, eta_max	Pa·s	per-rock clamp (also global in [rheology])
strain_gate_plastic ([rheology])	bool	gate strain accumulation by yield fraction
grain_exponent ([rheology])	—	legacy $d$-dependence of the eta0 path