Damage, softening and healing

MILEVA has two independent ways of weakening material that has deformed:

mechanism	state variable	weakens	heals	switch
grid damage	$D\in [0,1]$ on the cells	the whole viscosity, by up to damage_max_weakening	yes, over damage_heal_time	cfg.useDamage
strain weakening	accumulated strain eps_acc on the markers	the cohesion of the yield stress	no	cfg.useFullPlasticity with cfg.strainWeakeningFactor > 1

Both default off.

Physical background

Why deformation history matters

A viscosity law of the kind described on the rheology page depends on the present state of a parcel, that is its temperature, pressure, composition and strain rate. Real rocks also remember how they have deformed. A shear zone that has accommodated large strain is weaker than the same rock next to it, because deformation has reduced its grain size, aligned its minerals, ground its fault surfaces smooth or let fluids and melt in. Once formed, such a zone stays weak for some time after the deformation that created it has stopped, and can be reactivated later.

This memory matters most at the top of a convecting planet. A cold lithosphere following an Arrhenius law is many orders of magnitude stiffer than the mantle beneath it, and on its own it forms a single immobile lid, the stagnant-lid regime that is thought to describe Venus at present (Solomatov & Moresi 1996). Earth's lithosphere instead is broken into plates that are wide and rigid in their interiors but separated by narrow, weak boundaries. Plate velocities are a few centimetres per year (Schubert, Turcotte & Olson 2001), and a plate boundary that takes 5 cm/yr of relative motion across a zone of order 50 km width shears at roughly $3\times {10}^{-14}$ s${}^{-1}$ (velocity over width, so a strain-rate invariant of about half that), while the plate interior next to it barely deforms at all. Producing that contrast requires a rheology in which deformation makes rock weaker, so that a zone that starts to deform faster concentrates more deformation, a positive feedback called strain localisation (Bercovici, Ricard & Richards 2000). Plastic yielding of the lithosphere alone can break the lid (Moresi & Solomatov 1998), but the weak zones it creates vanish as soon as the stress drops. Weakening that depends on accumulated deformation, and recovers only slowly, lets weak zones persist and be inherited from one tectonic episode to the next (Bercovici & Ricard 2014).

Grid damage on a row of cells across a 50 km wide zone of a ${10}^{23}$ Pa s lid, deforming at a peak strain-rate invariant of $1.5\times {10}^{-14}$ s${}^{-1}$ (about 5 cm/yr across the zone) for 20 Myr and then not at all, stepped by evolve_damage! and apply_damage_to_viscosity! every 50 kyr. A: the imposed strain rate and the viscosity after 20 Myr (dashed) for three accumulation rates; the default weakens every cell above the threshold by the full factor 100. B: the viscosity at the zone centre. When the deformation stops the zone stays weak and recovers over the 100 Myr healing time, so it outlives the flow that made it. Script: docs/scripts/fig_damage_sections.jl.

Weakening mechanisms in rocks

Several microphysical processes weaken deforming rock, and they act in different parts of the lithosphere.

• Loss of cohesion and friction in the brittle regime. At shallow depth rock strength is limited by frictional sliding on fractures, described to first order by Byerlee's law, a shear stress of about $0.85{\sigma }_n$ for normal stresses below 200 MPa and of about $50\text{MPa}+0.6{\sigma }_n$ above it (Byerlee 1978). As a fault accumulates offset its cohesion and effective friction can fall. Continuum models represent this as cohesion or friction that decreases with accumulated plastic strain (Lavier, Buck & Poliakov 2000; Huismans & Beaumont 2003).

• Grain-size reduction in the ductile regime. Dynamic recrystallisation during dislocation creep reduces the grain size, which in turn speeds up grain-size-sensitive diffusion creep and lowers the viscosity. This is the basis of the grain-damage theory of Bercovici & Ricard (2012) and of the damage equation for grain size of Rozel, Ricard & Bercovici (2011). MILEVA models it explicitly on the grain size page.

• Fabric and other effects. The development of a fabric, in which weak phases become aligned and interconnected, can weaken ductile shear zones (Montési 2013). Shear heating, fluids and melt act in the same direction.

Opposing all of these is healing. Grain growth, crack sealing and the re-equilibration of fabrics return the rock towards its undeformed strength. The balance between damage and healing decides whether weak zones are persistent features or short-lived ones. In convection models with damage the occurrence of plate-like behaviour depends on that balance, and a hotter surface, which speeds up healing, favours a stagnant lid (Landuyt & Bercovici 2009; Foley, Bercovici & Landuyt 2012). This is one proposed reason why Venus, with a surface several hundred kelvin hotter than Earth's, lacks plate tectonics.

A: Byerlee's law with the normal stress taken as ${\rho }_{0}gz$, against the yield stress of composite_rheology (creep branches switched off, default cohesion 20 MPa and friction 0.2, cap 300 MPa), intact and after an accumulated strain of 2 with strainWeakeningFactor = 5. Cohesion loss lowers the strength near the surface by 16 MPa and matters little where friction dominates. B: the diffusion-creep viscosity of composite_rheology against grain-size reduction at 1400 K and 60 km: halving the grains weakens the rock eight times, a tenfold reduction a thousand times. Script: docs/scripts/fig_damage_sections.jl.

How MILEVA represents them

MILEVA does not model these microphysical processes one by one in the two mechanisms on this page. Each is a phenomenological state variable that stands for their combined effect.

• Grid damage is a single number $D$ between 0 (intact) and 1 (fully damaged) per cell. It grows where the strain rate exceeds a threshold, heals exponentially with a fixed time constant, and multiplies the whole viscosity. It has the structure of the damage models used to generate plates in convection calculations (Bercovici & Ricard 2012; Foley, Bercovici & Landuyt 2012), but is driven by the strain rate rather than by the deformational work, and has no grain-size or pinning physics behind it.

• Strain weakening is the brittle-regime softening of the continuum lithosphere models, in which the accumulated strain carried by each marker lowers the cohesion of the yield stress (Lavier, Buck & Poliakov 2000; Huismans & Beaumont 2003). Strain-dependent weakening has also been used in global convection models (Tackley 2000, part 2). It never heals.

One cell deformed at ${10}^{-14}$ s${}^{-1}$ for 10 Myr and then left alone. A: grid damage from evolve_damage! ($k_{\text{acc}}=10$) saturates within a million years and heals over ${\tau }_{\text{heal}}$; the strain-weakening fraction $W$ of the integrated ${\epsilon }_{\text{acc}}$ reaches 1 at a strain of 1.5, after 4.8 Myr, and stays there. B: what each one weakens, the viscosity through $1-w_{\text{max}}D$ and the cohesion of the yield branch of composite_rheology with strainWeakeningFactor = 5 at 10 km. Script: docs/scripts/fig_damage_sections.jl.

Scales and dimensionless numbers

Three ratios control what grid damage does in a run.

1. Damage number. The product $\mathrm{\Lambda }=k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }{\tau }_{\text{heal}}$ compares the rate of damage accumulation with the rate of healing. Its steady state is $D_{\text{ss}}=\mathrm{\Lambda }/(1+\mathrm{\Lambda })$ (derived below). For $\mathrm{\Lambda }\ll 1$ the field stays close to zero; for $\mathrm{\Lambda }\gg 1$ it saturates at 1.

2. Memory ratio. The healing time compared with the convective time scale, for example the time a parcel takes to cross the mantle, decides whether a weak zone outlives the flow that made it. With the default ${\tau }_{\text{heal}}\approx 100$ Myr a zone abandoned by the flow stays weak for about that long.

3. Threshold against the flow. The threshold ${\dot{\epsilon }}_{\text{thr}}$ sets which parts of the domain can accumulate damage at all. As an order-of-magnitude estimate, a mantle flowing at 1 cm/yr over a length of 1000 km deforms at about $3\times {10}^{-16}$ s${}^{-1}$, so the default threshold of ${10}^{-17}$ s${}^{-1}$ lies below typical convective strain rates and above those of a rigid lid.

The strength of the weakening is set by $w_{\text{max}}$. A fully damaged cell is $1/(1-w_{\text{max}})$ times weaker, 100 by default. A cold lid can be many orders of magnitude stiffer than the mantle below it, so whether that factor, together with yielding, is enough to mobilise the lid is an outcome of the calculation and not something the damage parameters decide alone.

A: the steady damage of evolve_damage! run to equilibrium (1000 Myr in 50 kyr steps) at strain rates from just above the threshold to ${10}^{-14}$ s${}^{-1}$ and three accumulation rates, against $D_{\text{ss}}=\mathrm{\Lambda }/(1+\mathrm{\Lambda })$; all points collapse onto the one curve of $\mathrm{\Lambda }$. B: the time scale on which $D$ approaches $D_{\text{ss}}$, against the strain rate. With the default ${10}^{14}$ it is seconds to hours at any strain rate above the threshold, far below a 50 kyr time step, which is why the default field acts as a switch. Script: docs/scripts/fig_damage_sections.jl.

Grid damage

Every time step, each grid cell does a small piece of bookkeeping. It first asks how much faster it is deforming than a threshold strain rate. That excess drives damage up, while healing pulls damage back down in every cell, in proportion to the damage already there. Their difference is applied over one time step, the result is kept between 0 and 1, and the new damage then lowers the cell's viscosity before the flow is solved. The diagram follows one deforming cell and one quiet cell through these steps.

One time step of grid damage, stage by stage, with the default damage parameters and $\mathrm{\Delta }t$ = 50 kyr. Right: a shear-zone cell (${\dot{\epsilon }}_{II}={10}^{-15}$ s${}^{-1}$, $D=0.30$, $\eta ={10}^{23}$ Pa s) overshoots to $D^{\ast }\approx {10}^{11}$, is clamped to $D=1$ and drops to ${10}^{21}$ Pa s; a quiet cell below the threshold (${10}^{-18}$ s${}^{-1}$, $D=0.60$, $\eta ={10}^{22}$ Pa s) only heals, to $D=0.5997$, and is weakened by a factor 0.406 to $4.1\times {10}^{21}$ Pa s. The values are checked against evolve_damage! and apply_damage_to_viscosity!. Script: docs/scripts/fig_damage_pipeline.jl.

physics/damage.jl carries a scalar damage field $D$ with the dimensions of the cell arrays, a phenomenological damage variable in the spirit of the damage-rheology models of Bercovici & Ricard (2012). It evolves as

$$\frac{dD}{dt}=k_{\text{acc}}max({\dot{\epsilon }}_{II}-{\dot{\epsilon }}_{\text{thr}},0)(1-D)-\frac{D}{{\tau }_{\text{heal}}},$$

accumulating where the strain-rate invariant exceeds the threshold and healing everywhere at rate $1/{\tau }_{\text{heal}}$. evolve_damage! takes one explicit Euler step per time step and clamps $D$ to $[0,1]$. apply_damage_to_viscosity! then reduces the viscosity,

$$\eta \leftarrow \mathrm{clamp}[\eta (1-w_{\text{max}}D),{\eta }_{\text{min}},{\eta }_{\text{max}}],$$

so a fully damaged cell is $1/(1-w_{\text{max}})=100$ times weaker with the default $w_{\text{max}}=0.99$.

The symbols in these two equations are:

symbol	meaning	units	set by
$D$	damage, 0 intact and 1 fully damaged	none	evolved, cell field damage
${\dot{\epsilon }}_{II}$	second invariant of the strain-rate tensor, $\sqrt{\frac{1}{2}{\dot{\epsilon }}_{ij}{\dot{\epsilon }}_{ij}}$	s${}^{-1}$	grid field strainII
${\dot{\epsilon }}_{\text{thr}}$	threshold below which no damage accumulates	s${}^{-1}$	cfg.damage_strain_rate_threshold
$k_{\text{acc}}$	accumulation rate per unit strain above the threshold	none	cfg.damage_acc_rate
${\tau }_{\text{heal}}$	healing time constant	s	cfg.damage_heal_time
$w_{\text{max}}$	fractional viscosity reduction at $D=1$	none	cfg.damage_max_weakening
$\eta$	viscosity from the rheology, before and after the damage factor	Pa s	ETA
${\eta }_{\text{min}}$, ${\eta }_{\text{max}}$	viscosity bounds	Pa s	cfg.eta_min, cfg.eta_max

Each term has a simple reading. The source term is proportional to the excess strain rate, so damage grows faster in faster-deforming cells, and to $1-D$, so it slows down as the cell approaches full damage and cannot exceed 1 in the continuous equation. The healing term is a linear relaxation towards the intact state with a single time constant, independent of temperature or strain rate. The viscosity factor $1-w_{\text{max}}D$ is linear in $D$, so half-damaged rock is weakened by about a factor 2 with the default $w_{\text{max}}$, and almost all of the factor 100 is reached only as $D$ approaches 1. The bounds ${\eta }_{\text{min}}$ and ${\eta }_{\text{max}}$ are the same as for the rest of the rheology, so the damaged viscosity stays inside them.

In the discrete update the source and sink are both evaluated with the damage at the start of the step, $D^{n+1}=\mathrm{clamp}(D^n+\mathrm{\Delta }t[k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }(1-D^n)-D^n/{\tau }_{\text{heal}}],0,1)$, which is why a large accumulation rate overshoots far above 1 before the clamp, as in the figure above. While the time step is still zero, as on the very first step, the driver passes $\mathrm{\Delta }t={10}^{10}$ s to evolve_damage! instead.

In the time loop. Damage is evolved and applied in step 2, after the viscosity law and the visco-elasto-plastic correction and before the Stokes solve, with the strain rate of the previous step. On the composite path the nonlinear iteration rebuilds the viscosity in every sweep and reapplies the damage factor to it, but does not evolve $D$ again: damage advances once per time step. $D$ is an Eulerian field: it stays at its grid position and is not carried with the flow. It is saved in checkpoints and restored on resume, and is not written to the VTK output.

Steady state and healing. Under a constant strain rate the damage approaches

$$D_{\text{ss}}=\frac{k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }{\tau }_{\text{heal}}}{1+k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }{\tau }_{\text{heal}}},\mathrm{\Delta }\dot{\epsilon }={\dot{\epsilon }}_{II}-{\dot{\epsilon }}_{\text{thr}},$$

on the time scale $1/(k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }+1/{\tau }_{\text{heal}})$, and once deformation stops it decays as $e^{-t/{\tau }_{\text{heal}}}$.

Default accumulation rate. $k_{\text{acc}}$ is a rate per unit strain. With the default ${10}^{14}$, $k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }\mathrm{\Delta }t$ is far above one for any mantle time step, so $D$ reaches 1 in the first step in which a cell exceeds the threshold and the field acts as a switch: a factor $1-w_{\text{max}}$ on the viscosity wherever ${\dot{\epsilon }}_{II}>{\dot{\epsilon }}_{\text{thr}}$, released over ${\tau }_{\text{heal}}$ after the deformation stops. Values of $k_{\text{acc}}$ of order 1 to 10 give the gradual accumulation of the loading figure below.

A worked example. Take a cell deforming steadily at ${\dot{\epsilon }}_{II}={10}^{-15}$ s${}^{-1}$ with the default threshold and healing time, so $\mathrm{\Delta }\dot{\epsilon }=9.9\times {10}^{-16}$ s${}^{-1}$ and ${\tau }_{\text{heal}}=3.15\times {10}^{15}$ s.

$k_{\text{acc}}$	$\mathrm{\Lambda }=k_{\text{acc}}\mathrm{\Delta }\dot{\epsilon }{\tau }_{\text{heal}}$	$D_{\text{ss}}$	time to approach $D_{\text{ss}}$	viscosity factor $1-w_{\text{max}}{D}_{\text{ss}}$
1	3.1	0.76	about 24 Myr	0.25
10	31	0.97	about 3 Myr	0.04
${10}^{14}$ (default)	$3\times {10}^{14}$	1	one time step	0.01

When the deformation stops, a fully damaged cell returns to $D=e^{-1}\approx 0.37$ after one healing time, about 100 Myr, where it is still weaker than intact rock by a factor $1/(1-0.99\times 0.37)\approx 1.6$. Damage falls below 0.01 only after $\ln 100\approx 4.6$ healing times, about 460 Myr. A weak zone created by grid damage therefore outlives the flow that created it by a time comparable with a mantle overturn.

Strain weakening of the yield stress

With cfg.useFullPlasticity = true the driver adds ${\dot{\epsilon }}_{II}\mathrm{\Delta }t$ to every marker's eps_acc after each solve (step 8, interpolated from the grid). The increment is the total strain of the step, in every cell, not only where the yield branch is active, and eps_acc never decreases. Markers created by reseeding take the cell value. The field is projected to the cells and read by composite_rheology, which weakens the cohesion of the Drucker-Prager yield stress:

$$C_{\text{eff}}=C_0[1-W(1-\frac{1}{f_{sw}})],W=\mathrm{clamp}(\frac{{\epsilon }_{\text{acc}}-{\epsilon }_1}{{\epsilon }_2-{\epsilon }_1},0,1),{\tau }_Y=\mathrm{clamp}(C_{\text{eff}}+\mu P,0,{\tau }_{\text{max}}).$$

The cohesion falls linearly from $C_0$ at ${\epsilon }_1$ to $C_0/f_{sw}$ at ${\epsilon }_2$. The friction coefficient $\mu$ and the cap ${\tau }_{\text{max}}$ are not weakened. With the default $f_{sw}$ = cfg.strainWeakeningFactor = 1 the strain is accumulated but the cohesion is unchanged. Strain weakening acts only on the composite path; the simple viscosity path and the plastic cap of visco-elasto-plasticity do not read eps_acc.

The symbols are:

symbol	meaning	units	set by
${\epsilon }_{\text{acc}}$	accumulated strain, the sum of ${\dot{\epsilon }}_{II}\mathrm{\Delta }t$ over all steps	none	marker field eps_acc, projected to the cells
${\epsilon }_1$, ${\epsilon }_2$	strain at which weakening starts and at which it is complete	none	cfg.strainWeakeningE1, cfg.strainWeakeningE2
$W$	weakening fraction, 0 before ${\epsilon }_1$ and 1 after ${\epsilon }_2$	none	computed
$f_{sw}$	ratio of intact to fully weakened cohesion	none	cfg.strainWeakeningFactor
$C_0$, $C_{\text{eff}}$	intact and weakened cohesion	Pa	mat.cohesion_upper, mat.cohesion_lower
$\mu$	friction coefficient	none	mat.friction_upper, mat.friction_lower
$P$	pressure, lithostatic (mantle_pressure) unless the caller passes a field	Pa	computed
${\tau }_{\text{max}}$	upper cap on the yield stress	Pa	mat.maxYieldStress
${\tau }_Y$	yield stress	Pa	computed

The yield stress enters the rheology as the yield viscosity ${\eta }_Y={\tau }_Y/(2{\dot{\epsilon }}_{II})$, combined with the creep viscosity as described on the rheology page. Lowering the cohesion therefore lowers the stress at which a cell starts to yield, and the effect is largest near the surface, where $\mu P$ is small and the cohesion carries most of the strength. At depth the frictional term dominates and cohesion loss changes the yield stress little. Because ${\epsilon }_{\text{acc}}$ is a strain, the weakening is set by how much a parcel has deformed and not by how fast. A parcel deforming at ${10}^{-14}$ s${}^{-1}$ reaches the default onset strain of 0.5 in about 1.6 Myr, one deforming at ${10}^{-16}$ s${}^{-1}$ in about 160 Myr. The default window of 0.5 to 1.5 is a strain of order one, the range over which continuum lithosphere models commonly let a fault lose its cohesion (Lavier, Buck & Poliakov 2000; Huismans & Beaumont 2003).

Left and middle: evolve_damage! and apply_damage_to_viscosity! on a single cell loaded at ${10}^{-15}$ s${}^{-1}$ for 100 Myr and then left undeformed, with $\mathrm{\Delta }t$ = 50 kyr and three accumulation rates. The default jumps to $D=1$ in the first step; all three heal on ${\tau }_{\text{heal}}\approx 100$ Myr (dashed: $e^{-t/{\tau }_{\text{heal}}}$). Right: the cohesion recovered from the yield branch of composite_rheology (creep and friction switched off, 10 km depth) against the accumulated strain, for three weakening factors.

Assumptions and simplifications

The two mechanisms are deliberately simple, and their simplifications are worth knowing when interpreting a run.

• Grid damage is driven by strain rate, not by work. Damage theories derive the damage rate from the fraction of the deformational work that goes into creating new grain boundaries or voids (Bercovici & Ricard 2012). MILEVA uses the excess strain rate, with no dependence on stress, temperature or grain size. Grain-size reduction as a physical process is covered separately by the grain size model.

• Healing is linear with one time constant. In rocks healing by grain growth is strongly temperature dependent, so hot regions heal faster than cold ones (Landuyt & Bercovici 2009). ${\tau }_{\text{heal}}$ is the same in every cell.

• Damage is Eulerian. $D$ stays at its cell and is not advected, so a damaged zone that is carried away by the flow leaves its damage behind. This is a reasonable approximation where the deformation and the material stay together, as in a plate boundary fixed relative to the grid, and a poor one where damaged material travels far before it heals. The accumulated strain eps_acc, in contrast, lives on the markers and moves with the material.

• The strain-rate history lags by one step. Damage is evolved once per step from the strain rate of the previous solve, so it does not take part in the nonlinear iteration of the current step.

• Strain weakening counts all strain and never heals. eps_acc sums the total strain, including viscous creep outside yielding regions, so material that has flowed far in the ductile mantle arrives at the surface already weakened. Brittle softening laws in lithosphere models usually count plastic strain only. Without healing a weakened marker stays weak for the rest of the run.

A: material moving at 2 cm/yr through a fixed 50 km zone for 30 Myr. The grid damage of evolve_damage! ($k_{\text{acc}}=10$) stays where the zone is, while the strain the markers accumulate, $\sum {\dot{\epsilon }}_{II}\mathrm{\Delta }t$ as in the driver, leaves with them, 0.7 for every parcel that has crossed. B: the cohesion of the yield branch of composite_rheology for a marker deforming steadily, strainWeakeningFactor = 5. At the ${10}^{-15}$ s${}^{-1}$ of ordinary ductile flow it starts to weaken after 16 Myr and is fully weakened after 48 Myr, and nothing restores it. Script: docs/scripts/fig_damage_sections.jl.

Scope

Both mechanisms are phenomenological. Their parameters are not calibrated against laboratory data, so results that depend on them should be presented together with a sensitivity test over $k_{\text{acc}}$, ${\tau }_{\text{heal}}$, $w_{\text{max}}$ or $f_{sw}$.

2-D and 3-D

evolve_damage! and apply_damage_to_viscosity! act element-wise on arrays of any dimension, and the marker strain accumulation and its projection have 2-D and 3-D methods. Both mechanisms behave identically in the annulus and on the Yin-Yang grid.

evolve_damage! on a 16 × 64 array and on a 16 × 12 × 36 array holding the same strain rates, $k_{\text{acc}}=10$, ${10}^{-15}$ s${}^{-1}$ for 100 Myr and then none. The update is element-wise, so a cell follows the same history to the last bit whatever the shape of the array. Script: docs/scripts/fig_damage_sections.jl.

Configuration

field	default	meaning
cfg.useDamage	false	evolve and apply grid damage
cfg.damage_acc_rate	${10}^{14}$	$k_{\text{acc}}$, per unit strain
cfg.damage_strain_rate_threshold	${10}^{-17}$ s${}^{-1}$	${\dot{\epsilon }}_{\text{thr}}$
cfg.damage_heal_time	$3.15\times {10}^{15}$ s (about 100 Myr)	${\tau }_{\text{heal}}$
cfg.damage_max_weakening	0.99	$w_{\text{max}}$
cfg.useFullPlasticity	false	accumulate eps_acc on the markers and weaken the cohesion
cfg.strainWeakeningFactor	1	$f_{sw}$, cohesion ratio at full weakening
cfg.strainWeakeningE1, cfg.strainWeakeningE2	0.5, 1.5	onset and saturation strain
mat.cohesion_upper, mat.cohesion_lower	20, 60 MPa	$C_0$

Where in the code

What	File, function
damage evolution and viscosity factor	physics/damage.jl: evolve_damage!, apply_damage_to_viscosity!
damage in the time loop and the nonlinear sweeps	main_SphCon.jl, step 2 and the Picard loop
strain accumulation on markers	main_SphCon.jl, step 8: interp_grid_to_markers!(M, strainII .* dt, grid, :eps_acc, :add)
projection of eps_acc	markers/markers_to_grid.jl: markers_to_grid
cohesion weakening	physics/composite_rheology.jl: composite_rheology, step 3
checkpointing of $D$	utils/save_checkpoint.jl, utils/load_checkpoint.jl

Tests

• tests/test_full_feature_smoke.jl: a short 2-D run with damage, eclogite kinetics and the other optional physics switched on stays finite.

• tests/test_config_timeloop.jl: a run with plastic yielding and useFullPlasticity survives the time loop.

• tests/test_config_matrix.jl: damage and full plasticity are off by default.

Key references

• Bercovici, D. & Ricard, Y. (2012). Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning. PEPI 202-203, 27-55.

• Foley, B. J., Bercovici, D. & Landuyt, W. (2012). The conditions for plate tectonics on super-Earths: inferences from convection models with damage. EPSL 331-332, 281-290.

• Lavier, L. L., Buck, W. R. & Poliakov, A. N. B. (2000). Factors controlling normal fault offset in an ideal brittle layer. JGR 105, 23431-23442.

Further reading

• Bercovici, D., Ricard, Y. & Richards, M. A. (2000). The relation between mantle dynamics and plate tectonics: a primer. In The History and Dynamics of Global Plate Motions, AGU Geophysical Monograph 121, 5-46. doi:10.1029/GM121p0005

• Bercovici, D. & Ricard, Y. (2014). Plate tectonics, damage and inheritance. Nature 508, 513-516. doi:10.1038/nature13072

• Byerlee, J. (1978). Friction of rocks. Pure and Applied Geophysics 116, 615-626. doi:10.1007/BF00876528

• Huismans, R. S. & Beaumont, C. (2003). Symmetric and asymmetric lithospheric extension: relative effects of frictional-plastic and viscous strain softening. JGR 108. doi:10.1029/2002JB002026

• Landuyt, W. & Bercovici, D. (2009). Variations in planetary convection via the effect of climate on damage. EPSL 277, 29-37. doi:10.1016/j.epsl.2008.09.034

• Montési, L. G. J. (2013). Fabric development as the key for forming ductile shear zones and enabling plate tectonics. Journal of Structural Geology 50, 254-266. doi:10.1016/j.jsg.2012.12.011

• Moresi, L. & Solomatov, V. (1998). Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the Earth and Venus. GJI 133, 669-682. doi:10.1046/j.1365-246x.1998.00521.x

• Rozel, A., Ricard, Y. & Bercovici, D. (2011). A thermodynamically self-consistent damage equation for grain size evolution during dynamic recrystallization. GJI 184, 719-728. doi:10.1111/j.1365-246X.2010.04875.x

• Schubert, G., Turcotte, D. L. & Olson, P. (2001). Mantle Convection in the Earth and Planets. Cambridge University Press. doi:10.1017/CBO9780511612879

• Solomatov, V. S. & Moresi, L.-N. (1996). Stagnant lid convection on Venus. JGR Planets 101, 4737-4753. doi:10.1029/95JE03361

• Tackley, P. J. (2000). Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations. 2. Strain weakening and asthenosphere. G3 1. doi:10.1029/2000GC000043