Stagnant lid that fails: Arrhenius viscosity and a yield stress (2-D)

When viscosity falls steeply with temperature, the cold top boundary layer becomes too stiff to take part in the flow. It stops moving and thickens into a conductive lid, and convection continues underneath it in a nearly isoviscous interior, driven by the small part of the temperature drop the mobile sublayer at the base of the lid can carry. This is the one-plate regime of terrestrial planets without plate tectonics. A real lid is not infinitely strong, though: once the stresses in it reach the strength of rock it can break, sink into the mantle and be replaced, which is the episodic regime. The example sets up both ingredients, an Arrhenius viscosity and a yield stress, in the spherical annulus, at a strength where the lid holds for hundreds of Myr and then fails. The viscosity-contrast example shows how the same shell moves from a mobile to a stagnant lid as the contrast grows, without yielding.

The model

Boussinesq, heated uniformly from within, with the core held at a fixed temperature, free slip and isothermal walls, in the common shell. The internal heating is what lets the lid reach a steady thickness: the lid grows until it conducts away what the mantle produces. Heated from below alone, the interior sits just under the core temperature, the bottom boundary layer carries a Nusselt number of about one at this Rayleigh number and activation energy (Solomatov & Moresi 2000), and the lid thickens by conduction without levelling off. The viscosity comes from the simple path, compute_viscosity in physics/compute_viscosity.jl, which main_SphCon uses whenever cfg.useCompositeRheology is false:

$$\eta (T)={\eta }_{\mathrm{ref}}\exp \left[\frac{E}{R}(\frac{1}{T}-\frac{1}{T_{\mathrm{ref}}})\right],{\eta }_{\min }\le \eta \le {\eta }_{\max },$$

with every layer reference viscosity (mat.*_eta0) set to ${\eta }_{\mathrm{ref}}$, mat.E_act_scale $=E/R$ in kelvin and mat.V_act_scale $=0$ (no pressure dependence). ${\eta }_{\mathrm{ref}}$ is set from the interior Rayleigh number

$${\mathrm{Ra}}_i=\frac{{\rho }_{0}g\alpha \mathrm{\Delta }T{D}^3}{\kappa {\eta }_{\mathrm{ref}}}.$$

Between the 300 K surface and the 1600 K interior the unclamped law spans about forty orders of magnitude, so the lid sits on the clamp and the effective viscosity contrast is ${\eta }_{\max }/{\eta }_{\mathrm{ref}}$, the third argument.

On top of that, cfg.useSimpleYield caps the stress at a yield stress ${\tau }_y$ (the fifth argument), the pseudo-plastic lid of Moresi & Solomatov (1998):

$${\eta }_{\mathrm{eff}}=min(\eta (T),\frac{{\tau }_y}{2{\dot{\epsilon }}_{II}}),$$

with ${\dot{\epsilon }}_{II}$ the second invariant of the strain rate. The cap is applied by simple_yield! in physics/compute_viscosity.jl, which reads the cohesion, friction and mat.maxYieldStress; the driver sets the cohesion and the cap to ${\tau }_y$ and the friction to zero, so the yield stress is the same at every depth. Where the lid stress stays below ${\tau }_y$ nothing changes. Where it would exceed it, the viscosity drops until the stress equals ${\tau }_y$, and since that viscosity depends on the strain rate it produces, each step iterates Stokes solve, strain rate and viscosity (Picard) until the velocity settles.

Which Arrhenius law

The time loop applies the Arrhenius law in compute_viscosity, driven by mat.E_act_scale and mat.T_ref, which is what this example sets. The separate cfg.useArrhenius helper (utils/arrhenius_viscosity.jl) is not used by main_SphCon.

parameter	value
$E$	300 kJ/mol (mat.E_act_scale = 36 080 K)
$T_{\mathrm{ref}}$	1600 K
${\mathrm{Ra}}_i$	${10}^6$ (first argument); ${\eta }_{\mathrm{ref}}=1.36\times {10}^{21}$ Pa s
clamps	${\eta }_{\min }={\eta }_{\mathrm{ref}}/10$, ${\eta }_{\max }={\eta }_{\mathrm{ref}}\cdot \mathrm{\Delta }\eta =1.36\times {10}^{25}$ Pa s with $\mathrm{\Delta }\eta ={10}^4$ (third argument)
yield stress	${\tau }_y=8$ MPa, depth independent (fifth argument; cfg.useSimpleYield)
internal heating	$H={10}^{-11}$ W/kg, uniform (fourth argument; cfg.useInternalHeating, cfg.internalHeating_H), about 20 mW/m² of heat production per unit surface
$T_{\mathrm{cmb}}$, $T_{\mathrm{surf}}$	1700 K (held fixed), 300 K
initial state	interior at 1600 K, exponential boundary layers of 100 km (top) and 50 km (CMB), a degree-4 CMB perturbation of 20 K and 5 K noise
grid	Nr × 8 Nr (default 64 × 512)
run length	1500 Myr (cfg.maxTimeMyr)
solver	MUMPS direct, markers for temperature

The driver as documented (${\mathrm{Ra}}_i={10}^6$, $\mathrm{\Delta }\eta ={10}^4$, $H={10}^{-11}$ W/kg, ${\tau }_y=8$ MPa, 64 × 512), 8854 steps from 0 to 1500 Myr, 114 frames at 5 frames per second. VTK output is written every 8 steps, about 5 Myr apart while the lid is quiet; the animation keeps one frame per 20 Myr there and one in every 16 outputs while the lid fails, when the steps are much shorter, so each failure plays in slow motion. Left: temperature, shown between 1100 and 1800 K; everything colder than 1100 K, the lid, is the saturated dark blue band. Right: ${\log }_{10}$ viscosity, from about ${10}^{20}$ Pa s in the hottest interior to the ${\eta }_{\max }$ clamp of $1.4\times {10}^{25}$ Pa s in the lid. The lid thickens for 200 Myr, breaks at 214 Myr and sinks into the mantle within about 2 Myr, and a new lid grows over the cooled interior. It fails again, more slowly, between 841 and 856 Myr, with shorter bursts of surface motion at about 868 and 987 Myr, and is still regrowing at 1500 Myr.

The same run at 1, 214, 846 and 1497 Myr, same colour range: the lid before it fails, the first failure (a cold slab of lid sinking at the left), the second, and the regrown lid over an interior still cooler than before the first failure.

Running it

export GKSwstype=nul
julia --project=. -t 4 examples/annulus_stagnant_lid_2d.jl 1e6 64 1e4 1e-11 8

Arguments are Ra_i Nr Δη H τ_y, with $H$ in W/kg and ${\tau }_y$ in MPa; H = 0 heats the shell from below only and τ_y = 0 switches yielding off. These are also the defaults, so the bare command runs the documented case. At 64 × 512 the run above took about 14 h of wall time on two threads (5.7 s per step, Picard sweeps included) with another long run sharing the machine; the smoke test below gives an uncontended timing. At the end the driver prints the means of lid_mobility, lid_thick_km and Nu_top over the last third of the rows of diagnostics.csv.

Diagnostics

The lid is read from diagnostics.csv (utils/diagnostics_csv.jl):

column	meaning
lid_mobility	rms surface velocity over rms velocity of the whole mantle, in the no-net-rotation frame
lid_thick_km	depth at which the horizontally averaged temperature crosses 1300 K, interpolated between cell centres
Vsurf_rms_cmyr, Vrms_cmyr	the two velocities behind the mobility
Nu_top, Nu_bot	wall heat flux over the conductive flux of the same shell

What to expect

Whether the lid survives depends on how its stresses compare with ${\tau }_y$. The same case was run at seven yield stresses (all at $\mathrm{\Delta }\eta ={10}^4$, $H={10}^{-11}$ W/kg, 64 × 512, to 1500 Myr or beyond):

${\tau }_y$	lid failures	behaviour
10 MPa	none	stagnant; the lid levels off at 252 km after about 1300 Myr, lid_mobility $3\times {10}^{-4}$
8 MPa (documented)	214 Myr, 841 to 856 Myr, shorter bursts at 868 and 987 Myr	episodic; the lid regrows to 290 km by 1500 Myr
7 MPa	164 Myr, 1017 Myr, and a smaller burst at 309 Myr	episodic
5 MPa	131 Myr, repeatedly between 283 and 454 Myr, 1150 Myr	episodic, failures more frequent
3 MPa	112 Myr, 1003 Myr, with bursts until about 1075 Myr	episodic; the lid regrows to 257 km by 1500 Myr
2 MPa	from 99 Myr on, without pause	mobile; surface velocity 2 to 10 cm/yr over 1000 to 1500 Myr
1 MPa	from 73 Myr on, without pause	mobile; surface velocity 4 to 9 cm/yr over 1000 to 1500 Myr

In the documented case each failure shows as a burst in every column of diagnostics.csv: Vrms_cmyr rises from about 0.01 to 100 cm/yr at the first failure, the surface moves with the mantle (lid_mobility near 1), lid_thick_km collapses from 191 to 15 km, and q_surf_mWm2 peaks at about 410 mW/m². The sinking lid then fills the mantle with cold rock: the mean temperature falls from 1380 K at the first failure to 1160 K after the second, below the 1700 K core, so heat now flows out of the core. While the new lid is still thin, between about 350 and 420 Myr, the surface keeps moving at a few cm/yr before the lid stiffens again. In the quiet stretches the whole mantle is nearly still, which makes lid_mobility, a ratio of two tiny velocities, read 0.5 or so; read Vsurf_rms_cmyr (0.02 cm/yr or less) to see that the surface is not moving. Over 1000 to 1500 Myr the documented case gives:

quantity	value
lid_thick_km (1300 K isotherm)	243 km, growing
Vsurf_rms_cmyr	0.0075 cm/yr
Vrms_cmyr	0.013 cm/yr
q_surf_mWm2	25.4 mW/m²
q_cmb_mWm2	17.4 mW/m² (out of the core)

The energy budget closes to better than one part in a thousand per sample between failures (budget_resid); during a failure the steps are so short that the sampled interval misses most of the change and the residual rises to about 0.2 for a few samples. On 178 steps, during the failures and the recoveries after them, the Picard iteration stopped at its 50-sweep cap with the residual above its tolerance of ${10}^{-4}$ (37 such steps by step 2400, where the log reports $1.2\times {10}^{-4}$, and 178 by step 7200, late in the second failure, at $1.6\times {10}^{-3}$); picard_n and picard_rel in diagnostics.csv record it. With no yield stress (τ_y = 0) the lid can only thin by convective erosion, and the run sits in the stagnant column of the table at any contrast of ${10}^4$ or more.

Lower yield stresses: from episodic to mobile

${\tau }_y=3$ MPa at 125, 167, 1034 and 1175 Myr, temperature between 700 and 1800 K: just after the first failure, the thin new lid, the second failure, and the lid regrowing (it reaches 257 km by 1500 Myr) over the cold rock of the failed lids, which collects above the core. VTK output of this run is written every 1000 steps, so the strip shows the frames nearest to the failures rather than the failures themselves.

Below 3 MPa the lid never re-forms. Once it first fails, at 99 Myr for 2 MPa and 73 Myr for 1 MPa, the surface keeps moving with the mantle for the rest of the run (lid_mobility 0.8 to 1.1) and the shell becomes a mobile lid that recycles its cold top boundary layer continuously:

τ_y = 2 MPa

τ_y = 1 MPa

Temperature between 700 and 1800 K over 1500 Myr. The sinking lid cools the whole mantle: over 1000 to 1500 Myr the mean temperature is 843 K at 2 MPa and 771 K at 1 MPa, far below the 1700 K core, which therefore loses 240 and 290 mW/m². Since the mantle is colder than 1300 K almost everywhere, lid_thick_km has no isotherm to find and reports the whole shell; read Vsurf_rms_cmyr and lid_mobility instead. Wall time was 31 h (1 MPa) and 54 h (2 MPa) on two threads with other runs sharing the machine.

Smoke test

A quick check that the driver runs:

MILEVA_MAXSTEPS=6 julia --project=. -t 2 examples/annulus_stagnant_lid_2d.jl 1e6 16

6 steps at 16 × 128 with the default contrast, heating and yield stress in 37 s of wall time (mostly compilation), 2.2 GB; the driver reports a lid of 153 km and a lid mobility of ${10}^{-3}$. At 64 × 512, 40 steps took 0.70 s per step and 2.7 GB on two threads. Steps while the lid fails cost more, because the Picard iteration then needs up to tens of sweeps instead of two.