Skip to content

Grain-size evolution ​

With cfg.evolveGrainSize = true every marker carries a grain size d that grows by normal grain growth and is reduced by dynamic recrystallisation, after the wattmeter of Austin & Evans (2007), with an optional piezometric lower bound for grain-boundary sliding. The grain size feeds back on the diffusion-creep and grain-boundary-sliding viscosities of the composite rheology. The update is 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,

ε˙=Aσnd−mexp⁡(−E+PVRT),

where ε˙ is the strain rate, σ the differential stress, d the grain size, E and V the activation energy and volume, P the pressure, T the temperature and R the gas constant. Diffusion creep has n=1 and a grain-size exponent m of 2 when the diffusion runs through the grain interiors (Nabarro-Herring creep) and 3 when it runs along the grain boundaries (Coble creep) (Karato 2008); the olivine laws MILEVA provides use m=2.5 (Karato & Wu 1993) and m=3 (Hirth & Kohlstedt 2003). Dislocation creep of olivine has n≈3.5 and m=0 (Hirth & Kohlstedt 2003). Since the viscosity is η=σ/2ε˙, the diffusion-creep viscosity scales as dm. With m=3, halving the grain size makes the rock eight times weaker, and a tenfold change in grain size moves the viscosity by three orders of magnitude.

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 10−15 s−1 the diffusion-creep viscosity rises as d3 while dislocation creep does not depend on d; the composite follows whichever is weaker, and the two cross near the 1 mm reference size. B: the grain size at which the two branches are equally strong, against strain rate at three temperatures; finer rock deforms by diffusion creep, coarser rock by dislocation creep. Script: 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

dp−d0p=k0exp⁡(−EgRT)t,

with d0 the initial size, k0 a rate constant, Eg the activation energy of grain-boundary migration and p the growth exponent. Growth driven by boundary curvature alone gives p=2. Experiments on olivine and other rocks often give larger exponents, which is usually attributed to pinning of the boundaries by second phases, pores and impurities (Karato 1989; Evans et al. 2001). Growth always slows as the grains coarsen, because the curvature that drives it falls as 1/d, and it is strongly temperature dependent.

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 1/d, the reduction rate is proportional to the work rate σε˙ times d2 divided by the grain-boundary energy. The balance between this reduction and normal growth then sets a steady grain size that depends on temperature and on the work rate, and with published growth and creep laws it reproduces the measured relations between stress and recrystallised grain size of quartz and olivine rocks without fitting (Austin & Evans 2007). Rozel et al. (2011) derived a thermodynamically self-consistent evolution equation of the same kind for grain size during dynamic recrystallisation.

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 d∝t1/2, the default exponent p=2, until the 10 cm clamp. B: the growth rate and the reduction rate of the code at 1300 K and 100 km, recovered exactly from single steps. Growth falls as 1/d, reduction CWε˙II does not depend on d, and where they cross is the steady size, 4.59 mm at 10−13 s−1; at 10−15 s−1 the crossing lies beyond the upper clamp. Script: 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 10−15 and 10−13 s−1, reduction at 10−12 s−1, where the grain size falls to 0.46 mm. B: the diffusion-creep viscosity follows d3, ten times weaker for the shear-zone marker and four orders of magnitude stiffer for the slow one, but with the default constants dislocation creep carries more than 95 per cent of the strain at these conditions, so the composite viscosity moves by a few per cent at most. The feedback of the text acts only where diffusion creep carries a large share. Script: docs/scripts/fig_grain_size_sections.jl.

Scales and typical values ​

Two time scales compete. For growth with the exponent p the time to reach a size d from a much smaller one is

τgrow≈dppkT,kT=kge−Eg/RT,

and the time for reduction to remove a grain size d is τred≈d/d˙red. The grain size settles where the two are comparable, so hot, slowly deforming regions coarsen and cool, rapidly deforming regions refine. With MILEVA's default constants (below) and p=2, kT≈9.6×10−18 m2 s−1 at 1600 K and 200 km depth, so a grain grows to 1 mm in about 1.7 kyr, to 1 cm in about 0.17 Myr and to 10 cm in about 17 Myr. At 1300 K and 100 km depth growth to 1 mm takes about 34 kyr. These numbers follow from the defaults of the code, not from measured olivine kinetics, and they show why the grain size of a quiet mantle region at 1600 K reaches the upper clamp of 10 cm within about 20 Myr unless the constants are set for the material.

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, dp/(pkT), with the growth constant kT read from evolve_grain_size! at 200 km (and at 100 km for 1 mm). At 1600 K, kT=9.6×10−18 m² s⁻¹, as quoted above. B: the steady grain size kT/(CWε˙II) at 200 km, clamped to [dmin,dmax]; with the default constants most of the plane sits at the 10 cm clamp, and a steady size of 1 mm or less needs strain rates above about 10−14 s−1 even at 1100 K, and faster ones in hotter rock. Script: 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 Δt=1011 s from d=1 mm, default constants, cfg.evolveGrainSize = true. Shear-zone marker: 1300 K, 100 km, 10−12 s−1; reduction wins and d falls to 0.951 mm (steady size 0.459 mm), but its viscosity, 1.5×1020 Pa s, is set by yielding and does not change. Quiet mantle marker: 1600 K, 200 km, 10−15 s−1; growth wins and d rises to 1.956 mm (steady size capped at dmax), which stiffens diffusion creep 7.5 times and the composite viscosity from 5.3×1021 to 7.0×1021 Pa s. The grid viscosity each marker reads is the composite viscosity of its cell at 1 mm. Every value is checked against evolve_grain_size! and composite_rheology. Script: docs/scripts/fig_grain_size_pipeline.jl.

At each marker, with the viscosity η and strain-rate invariant ε˙II of the step interpolated from the grid (bilinear in 2-D, trilinear in 3-D) and T the marker temperature floored at 300 K:

d˙grow=kgexp⁡(−Eg/RT)dp−1,Eg=(200+0.075min(zkm,2000)) kJ mol−1,d˙red=CWΦσ,Φ=2ηε˙II2,σ=max(2ηε˙II, 1 kPa),

with zkm the depth below mat.R_surf in kilometres. The step is semi-implicit in the reduction term:

dn+1=dn+Δtd˙grow(dn)1+Δtd˙red/dn,

and the result is clamped to [dmin,dmax], mat.grain_d_min = 1 µm and mat.grain_d_max = 10 cm. To first order in Δt this is d˙=d˙grow−d˙red.

The symbols are

symbolmeaningcode
d, dnmarker grain size, at step n (m)M.d
d˙grow, d˙redgrowth and reduction rates (m s−1)d_dot_grow, d_dot_red
kggrowth prefactor (mp s−1)mat.grain_k_g
pgrowth exponentmat.grain_p
Egactivation energy of growth, rising linearly from 200 kJ mol−1 at the surface to 350 kJ mol−1 at 2000 km depth and constant belowE_g
Rgas constant, 8.314 J mol−1 K−1cfg.R_GAS
Tmarker temperature, floored at 300 K (K)M.T
η, ε˙IIgrid viscosity and second invariant of the strain rate, interpolated to the markereta_m, strain_m
Φ2ηε˙II2, half the viscous dissipation (W m−3)Phi_m
σstress 2ηε˙II, floored at 1 kPa (Pa)sigma_m
CWreduction constantmat.grain_C_W
Δttime step (s)dt

The growth term is the differential form of the classical growth law above. At fixed temperature, integrating d˙=kT/dp−1 gives dp=d0p+pkTt. The step treats reduction implicitly, dividing by 1+Δtd˙red/dn, so a strong reduction can never drive the grain size negative, and treats growth explicitly, evaluated at the old grain size.

Steady state. Wherever the stress exceeds 1 kPa, Φ/σ=ε˙II and the reduction rate is CWε˙II (CW then has units of length). The steady grain size is

dss=(kge−Eg/RTCWε˙II)1/(p−1),

a function of temperature, depth and strain rate only; the viscosity and the stress cancel. Φ here is 2ηε˙II2, half the dissipation 4ηε˙II2 the driver reports, a factor that is absorbed in CW.

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 σε˙d2, so the steady size depends on the work rate and falls as dssp+1∝kT/(σε˙). In MILEVA the ratio Φ/σ reduces to ε˙II and the reduction rate CWε˙II does not depend on d, so dss depends on the strain rate and not on the stress. The work rate is also the total one, from all deformation mechanisms, where the original model counts only the work done by dislocation creep. The two forms predict different stress and temperature dependences of the steady grain size, which is one more reason to calibrate kg and CW before a quantitative use.

Piezometric bound. With cfg.useGBS = true the grain size is kept above the grain-boundary-sliding piezometer,

d≥clamp(Apzσ−qpz, dmin, dmax),

with σ in pascals (mat.grain_A_pz = 10−2, 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 Apz.

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 d≡dref):

ηdiff∝(ddref)m,ηGBS∝(ddref)mg/ng,

with m = mat.diff_m (3), dref = mat.diff_d_ref (1 mm), mg = cfg.gbs_m and ng = cfg.gbs_n (both 2). Dislocation creep and the yield stress do not depend on d. The simple viscosity path ignores grain size. New markers created by reseeding take the cell grain size, and the field is written to the VTK output as 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, η∝dm/n. That is d3 for the Hirth & Kohlstedt (2003) and d2.5 for the Karato & Wu (1993) olivine diffusion laws. With cfg.useGBS = true and a laboratory dislocation law selected, grain-boundary sliding follows the Hirth & Kohlstedt (2003) dislocation-accommodated law (m=2, n=3.5, hence η∝d4/7). The crustal laws read the same grain size. With 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 d, lowers ηdiff and ηGBS, and lets that region deform faster; growth in a quiet region raises them. Because the composite viscosity is a harmonic sum, the grain size changes the viscosity only where diffusion creep or grain-boundary sliding carries a significant share of the deformation. In the shear-zone example of the figure above yielding sets the viscosity, so the smaller grains change nothing.

With the default constants the steady size reaches dmax at 1600 K for strain rates below about 10−13 s−1 (right panel of the figure below), where diffusion creep is (dmax/dref)3=106 times stiffer than at the reference size. Set kg, CW and the clamps for the material before relying on the coupling.

Left: a marker starting at 0.01 mm at 1600 K, 200 km depth and 10−15 s−1, from evolve_grain_size! with time steps of 100 yr and 10 kyr, against parabolic growth d2=d02+2kTt (the reduction term is negligible here). The growth rate is taken at the old grain size, so a step long compared with d2/kT overshoots; the 10 kyr step reaches dmax in one step. Right: grain size after 6000 steps of 1012 s at 200 km depth, for three temperatures and two viscosities (symbols overlap), against kT/(CWε˙II) clamped to [dmin,dmax].

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 10−15 s−1 (growth to the clamp) and at 1300 K and 10−13 s−1 (a steady size). The two methods differ only in how they interpolate the grid fields to the marker, so on uniform fields they agree to the last bit. Script: 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 p and kg. 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_g is not read.

  • Simplified reduction. The reduction rate is CWε˙II, 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 dp/kT overshoots (left panel of the figure above).

A: the growth constant kT at 1600 K against depth, read from 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 kg, CW, p and the clamps should be set for the material and the question before grain-size results are interpreted.

Configuration ​

fielddefaultmeaning
cfg.evolveGrainSizefalseevolve marker grain size and use it in the composite rheology
mat.grain_k_g10−10 mp s−1growth prefactor kg
mat.grain_p2growth exponent p
mat.grain_C_W10−3reduction constant CW
mat.grain_d_min, mat.grain_d_max1 µm, 0.1 mclamps
cfg.useGBSfalsegrain-boundary sliding in the rheology and the piezometric bound
mat.grain_A_pz, mat.grain_q_pz10−2, 1.5piezometer
mat.diff_m, mat.diff_d_ref3, 1 mmgrain-size exponent and reference size of diffusion creep
cfg.gbs_A, cfg.gbs_n, cfg.gbs_m6.5×10−16, 2, 2grain-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 ​

WhatFile, function
evolution law, 2-D and 3-Dphysics/evolve_grain_size.jl: evolve_grain_size!, _bilinear_g2m_gs, _trilinear_g2m_gs
use in the viscosityphysics/composite_rheology.jl: composite_rheology (d_ratio)
marker field d, projection, reseedingmarkers/init_markers.jl, markers/markers_to_grid.jl, markers/reseed_markers.jl
call in the time loopmain_SphCon.jl, step 9
material defaultsconfig/MaterialProperties.jl

Tests ​

  • tests/test_rheology_switches.jl: the viscosity does not depend on the grain size passed in when cfg.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