Grain-size evolution
With cfg.evolveGrainSize = true every marker carries a grain size evolve_grain_size! (physics/evolve_grain_size.jl), step 9 of the time loop, after the Stokes solve and the strain rate of the step.
Physical background
Why grain size controls viscosity
Mantle rocks are polycrystalline aggregates, in the upper mantle dominated by olivine. At mantle temperatures they flow by thermally activated creep, and two families of mechanisms carry most of the strain (Karato & Wu 1993; Hirth & Kohlstedt 2003). In diffusion creep atoms migrate, through the grain interiors or along the grain boundaries, from boundaries under compression to boundaries under tension. The diffusion distance is the grain size, and the number of boundaries that act as sources and sinks grows as the grains get smaller, so the strain rate rises steeply as the grain size falls. In dislocation creep the crystals deform by the glide and climb of lattice dislocations inside the grains, and the strain rate does not depend on grain size. Laboratory flow laws express both with one form,
where
The two mechanisms operate in parallel, so the faster one dominates. At given temperature, pressure and stress there is a transition grain size. Finer rock deforms mainly by diffusion creep, coarser rock mainly by dislocation creep. A deformation-mechanism map in the plane of stress and grain size shows the two fields and the boundary between them (Karato & Wu 1993; Karato 2008). Whether a region of the mantle has a linear (Newtonian) or a strongly stress-dependent (power-law) rheology therefore depends on its grain size as much as on its temperature.

From composite_rheology with cfg.evolveGrainSize = true, each branch isolated and the viscosity clamps lifted. A: at 1600 K, 100 km depth and docs/scripts/fig_grain_size_sections.jl.
How grain size changes
Grain size is a state variable of the rock and not a material constant. Two processes change it at the same time.
Normal grain growth. Grain boundaries carry a surface energy. They migrate towards their centres of curvature, so large grains consume small ones and the mean grain size increases while the total boundary area falls. The classical kinetic law is
with
Dynamic recrystallisation. During dislocation creep new grain boundaries form, by the progressive rotation of subgrains and by the migration of boundaries driven by differences in stored dislocation energy, and they cut the grains into smaller ones. In experiments the recrystallised grain size depends mainly on the stress and falls as a power of it, which is the basis of the recrystallised grain-size piezometer (Karato et al. 1980). De Bresser et al. (1998) proposed that recrystallisation drives the grain size towards the boundary between the dislocation and diffusion creep fields, because grains that become fine enough to deform by diffusion creep no longer recrystallise.
The wattmeter. Austin & Evans (2007) related the rate of grain-size reduction to the rate of mechanical work. If a fixed fraction of the work done by dislocation creep is stored as new grain-boundary area, and the boundary area per unit volume scales as

A: evolve_grain_size! on one marker at 200 km with no deformation, from 0.01 mm, with steps a twentieth of the elapsed time. After a short start each curve follows docs/scripts/fig_grain_size_sections.jl.
Why it matters for mantle and lithosphere dynamics
Grain-size reduction couples deformation to weakening. Faster deformation shrinks the grains, smaller grains make diffusion creep easier, and the weaker rock deforms faster still. Karato & Wu (1993) concluded that near the transition between the two creep regimes this feedback weakens the upper mantle and localises deformation. Natural peridotite shear zones record the process. Their mylonites are much finer grained than the host peridotite and deform by grain-size-sensitive creep (Warren & Hirth 2006). Bercovici & Ricard (2012) built a theory of plate generation on grain damage, in which grain-size reduction in a two-phase rock whose boundaries are pinned by the second phase produces narrow, long-lived weak zones. Bercovici & Ricard (2014) suggested that the hot surface of Venus, by speeding up grain growth in the lithosphere, heals such weak zones too fast for plate boundaries to develop.
The feedback also works in the other direction. Where stresses are low, grains grow and diffusion creep stiffens, so grain-size evolution can add viscosity contrasts that temperature alone would not produce. Hall & Parmentier (2003) found that it affects the onset of convective instability beneath a cooling boundary layer and can create viscosity contrasts exceeding one order of magnitude in finite-amplitude convection.

A marker at 1300 K and 100 km held at a fixed strain rate for 20 Myr, with evolve_grain_size! reading the composite viscosity of its current grain size at every 1 kyr step (no yield law). A: growth wins at docs/scripts/fig_grain_size_sections.jl.
Scales and typical values
Two time scales compete. For growth with the exponent
and the time for reduction to remove a grain size
Olivine grain sizes in upper-mantle peridotite are of order millimetres (Karato & Wu 1993; Karato 2008), which is the scale of the 1 mm reference size of the diffusion-creep law in MILEVA. The mylonites of mantle shear zones are much finer grained than their host peridotites (Warren & Hirth 2006). No mantle rock has been sampled on Venus, so its grain size is not constrained by observation and enters a model as a parameter.

A: the time to grow to 1 mm, 1 cm and 10 cm, evolve_grain_size! at 200 km (and at 100 km for 1 mm). At 1600 K, docs/scripts/fig_grain_size_sections.jl.
Evolution law
Each marker's grain size changes in every step through a contest between two processes. Heat makes grains grow, faster when the mantle is hot and the grains are still small, while deformation breaks them up by recrystallisation, faster when the rock is strained quickly. The code works out both rates at the marker, combines the two in one step, keeps the result within fixed bounds, and passes the grain size on to diffusion creep, where smaller grains make the rock weaker. The figure follows one step for two markers.
One step of cfg.evolveGrainSize = true. Shear-zone marker: 1300 K, 100 km, evolve_grain_size! and composite_rheology. Script: docs/scripts/fig_grain_size_pipeline.jl.
At each marker, with the viscosity
with mat.R_surf in kilometres. The step is semi-implicit in the reduction term:
and the result is clamped to mat.grain_d_min = 1 µm and mat.grain_d_max = 10 cm. To first order in
The symbols are
| symbol | meaning | code |
|---|---|---|
| marker grain size, at step | M.d | |
| growth and reduction rates (m s | d_dot_grow, d_dot_red | |
| growth prefactor (m | mat.grain_k_g | |
| growth exponent | mat.grain_p | |
| activation energy of growth, rising linearly from 200 kJ mol | E_g | |
| gas constant, 8.314 J mol | cfg.R_GAS | |
| marker temperature, floored at 300 K (K) | M.T | |
| grid viscosity and second invariant of the strain rate, interpolated to the marker | eta_m, strain_m | |
Phi_m | ||
| stress | sigma_m | |
| reduction constant | mat.grain_C_W | |
| time step (s) | dt |
The growth term is the differential form of the classical growth law above. At fixed temperature, integrating
Steady state. Wherever the stress exceeds 1 kPa,
a function of temperature, depth and strain rate only; the viscosity and the stress cancel.
Relation to the wattmeter of Austin & Evans (2007). MILEVA takes the structure of the wattmeter, normal growth against a reduction driven by the deformational work, but not its grain-size dependence. In the original model the reduction rate scales with
Piezometric bound. With cfg.useGBS = true the grain size is kept above the grain-boundary-sliding piezometer,
with mat.grain_A_pz = mat.grain_q_pz = 1.5). With these constants the bound is below 1 µm for any stress above 1 kPa, so it becomes active only with a larger
Coupling to the viscosity
The marker grain size is projected to the cells as a weighted cell mean (dg from markers_to_grid) and passed to composite_rheology, which uses it only when cfg.evolveGrainSize = true (otherwise
with mat.diff_m (3), mat.diff_d_ref (1 mm), cfg.gbs_m and cfg.gbs_n (both 2). Dislocation creep and the yield stress do not depend on GrainSize [m].
When published laboratory laws replace the reference-viscosity laws (mat.flowlaw_diffusion or mat.flowlaw_dislocation set, see Flow laws), the cell grain size enters them in absolute terms, floored at 1 µm, with the exponents of the law, cfg.useGBS = true and a laboratory dislocation law selected, grain-boundary sliding follows the Hirth & Kohlstedt (2003) dislocation-accommodated law (cfg.evolveGrainSize = false every laboratory law uses the fixed mat.flowlaw_grain_size (1 mm by default).
Physically, the coupling closes the feedback loop described in the background section. Recrystallisation in a fast-deforming region lowers
With the default constants the steady size reaches

Left: a marker starting at 0.01 mm at 1600 K, 200 km depth and evolve_grain_size! with time steps of 100 yr and 10 kyr, against parabolic growth
2-D and 3-D
evolve_grain_size! has a 2-D method (bilinear interpolation on the annulus) and a 3-D method (trilinear interpolation on the Yin patch), and the time loop calls it in both. The physics is identical.

The same marker evolved by the 2-D method on an annulus and by the 3-D method on a Yin patch (8 radial cells each), with uniform viscosity and strain-rate fields and 10 kyr steps for 30 Myr, at 1600 K and docs/scripts/fig_grain_size_sections.jl.
Assumptions and simplifications
One grain size per marker. Each marker carries a single representative grain size rather than a grain-size distribution, and the cell value is the weighted mean of its markers. Rocks with a wide or bimodal distribution, such as partly recrystallised shear zones, are represented by their mean.
One mineral. The same growth and reduction constants apply at every depth and in every composition. Pinning by a second phase enters only implicitly, through the values of
and . Phase transitions, melting and melt extraction do not reset the grain size. Depth instead of pressure. The pressure dependence of growth is represented by an activation energy that rises linearly with depth, rather than by an activation volume multiplying the pressure;
mat.grain_E_gis not read.Simplified reduction. The reduction rate is
, driven by the total strain rate and independent of grain size (see the relation to the wattmeter above). Grid quantities at the marker. Viscosity and strain rate are interpolated from the grid, so the reduction term sees cell-scale values from the current step, after the Stokes solve.
Explicit growth. Growth is evaluated at the old grain size, so a time step long compared with
overshoots (left panel of the figure above).

A: the growth constant evolve_grain_size!. The activation energy rising linearly with depth slows growth by more than four orders of magnitude between the surface and 2000 km, and is constant below. B: a cell holding two marker populations, 0.05 mm and 5 mm grains, in diffusion creep at 1600 K and 100 km. The cell uses the viscosity of the mean grain size (black). If the two populations deformed side by side at equal stress the cell would be as weak as its fine fraction allows (blue); in series at equal strain rate, as stiff as its coarse fraction (dashed). A few per cent of fine markers changes the first by orders of magnitude and the cell mean hardly at all. Script: docs/scripts/fig_grain_size_sections.jl.
Calibration
The default constants are generic. With them the steady size reaches the upper clamp in quiet mantle, so
Configuration
| field | default | meaning |
|---|---|---|
cfg.evolveGrainSize | false | evolve marker grain size and use it in the composite rheology |
mat.grain_k_g | growth prefactor | |
mat.grain_p | 2 | growth exponent |
mat.grain_C_W | reduction constant | |
mat.grain_d_min, mat.grain_d_max | 1 µm, 0.1 m | clamps |
cfg.useGBS | false | grain-boundary sliding in the rheology and the piezometric bound |
mat.grain_A_pz, mat.grain_q_pz | piezometer | |
mat.diff_m, mat.diff_d_ref | 3, 1 mm | grain-size exponent and reference size of diffusion creep |
cfg.gbs_A, cfg.gbs_n, cfg.gbs_m | grain-boundary sliding law |
The activation energy of growth is the depth-dependent expression above; mat.grain_E_g is not read. The initial grain size is 1 mm on every marker.
Where in the code
| What | File, function |
|---|---|
| evolution law, 2-D and 3-D | physics/evolve_grain_size.jl: evolve_grain_size!, _bilinear_g2m_gs, _trilinear_g2m_gs |
| use in the viscosity | physics/composite_rheology.jl: composite_rheology (d_ratio) |
marker field d, projection, reseeding | markers/init_markers.jl, markers/markers_to_grid.jl, markers/reseed_markers.jl |
| call in the time loop | main_SphCon.jl, step 9 |
| material defaults | config/MaterialProperties.jl |
Tests
tests/test_rheology_switches.jl: the viscosity does not depend on the grain size passed in whencfg.evolveGrainSize = false.tests/test_config_matrix.jl: grain-size evolution is off by default.
Key references
Austin, N. J. & Evans, B. (2007). Paleowattmeters: a scaling relation for dynamically recrystallized grain size. Geology 35, 343-346. doi:10.1130/G23244A.1
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
Ohuchi, T., Kawazoe, T., Higo, Y., Funakoshi, K., Suzuki, A., Kikegawa, T. & Irifune, T. (2015). Dislocation-accommodated grain boundary sliding as the major deformation mechanism of olivine in the Earth's upper mantle. Science Advances 1, e1500360. doi:10.1126/sciadv.1500360
Further reading
Bercovici, D. & Ricard, Y. (2012). Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning. Physics of the Earth and Planetary Interiors 202-203, 27-55. doi:10.1016/j.pepi.2012.05.003
Bercovici, D. & Ricard, Y. (2014). Plate tectonics, damage and inheritance. Nature 508, 513-516. doi:10.1038/nature13072
De Bresser, J. H. P., Peach, C. J., Reijs, J. P. J. & Spiers, C. J. (1998). On dynamic recrystallization during solid state flow: effects of stress and temperature. Geophysical Research Letters 25, 3457-3460. doi:10.1029/98GL02690
Evans, B., Renner, J. & Hirth, G. (2001). A few remarks on the kinetics of static grain growth in rocks. International Journal of Earth Sciences 90, 88-103. doi:10.1007/s005310000150
Hall, C. E. & Parmentier, E. M. (2003). Influence of grain size evolution on convective instability. Geochemistry, Geophysics, Geosystems 4. doi:10.1029/2002GC000308
Hirth, G. & Kohlstedt, D. L. (2003). Rheology of the upper mantle and the mantle wedge: a view from the experimentalists. In Inside the Subduction Factory, Geophysical Monograph 138, 83-105. doi:10.1029/138GM06
Karato, S. (1989). Grain growth kinetics in olivine aggregates. Tectonophysics 168, 255-273. doi:10.1016/0040-1951(89)90221-7
Karato, S. (2008). Deformation of Earth Materials. Cambridge University Press. doi:10.1017/CBO9780511804892
Karato, S., Toriumi, M. & Fujii, T. (1980). Dynamic recrystallization of olivine single crystals during high-temperature creep. Geophysical Research Letters 7, 649-652. doi:10.1029/GL007i009p00649
Karato, S. & Wu, P. (1993). Rheology of the upper mantle: a synthesis. Science 260, 771-778. doi:10.1126/science.260.5109.771
Warren, J. M. & Hirth, G. (2006). Grain size sensitive deformation mechanisms in naturally deformed peridotites. Earth and Planetary Science Letters 248, 438-450. doi:10.1016/j.epsl.2006.06.006