## Abstract

At high pressure prevailing in the lower mantle, lattice friction opposed to dislocation glide becomes very high, as reported in recent experimental and theoretical studies. We examine the consequences of this high resistance to plastic shear exhibited by ringwoodite and bridgmanite on creep mechanisms under mantle conditions. To evaluate the consequences of this effect, we model dislocation creep by dislocation dynamics. The calculation yields to an original dominant creep behavior for lower mantle silicates where strain is produced by dislocation climb, which is very different from what can be activated under high stresses under laboratory conditions. This mechanism, named pure climb creep, is grain-size–insensitive and produces no crystal preferred orientation. In comparison to the previous considered diffusion creep mechanism, it is also a more efficient strain-producing mechanism for grain sizes larger than ca. 0.1 mm. The specificities of pure climb creep well match the seismic anisotropy observed of Earth’s lower mantle.

- Mantle convection Creep mechanism bridgmanite

## INTRODUCTION

Earth dissipates its internal heat through large-scale mantle convection. Although made of crystalline rocks, the mantle can flow like a viscous fluid at geological time scale. This behavior is governed by creep mechanisms, which involve the motion of crystal defects at the microscopic scale. The nature of the defects and the creep mechanisms involved have profound implications on the rheology and hence, on the dynamics of our planet. Recent progress in high-pressure experiments has expanded our capability to perform deformation experiments at high pressure and at high temperature. Recently, a few studies have shown that ringwoodite and bridgmanite can be deformed experimentally provided that extremely high stresses (of the order of the GPa) are applied (*1*–*3*). However, these data were acquired at laboratory strain rates (of the order of 10^{−5} s^{−1}) because reproducing creep deformation at high *P*,*T* at the extremely slow strain rates of the mantle (of the order of 10^{−14} s^{−1}) is still out of reach. Here, we report results from numerical simulations of creep, which reveal the importance of a mechanism where dislocations are sources and sinks for vacancy diffusion and produce plastic shear by climb. By introducing a characteristic distance for diffusion smaller than the grain size, this pure climb creep mechanism is found to be more efficient than the standard Nabarro-Herring (NH) (or Coble) creep features. Pure climb creep, which is grain size–independent, is then a very efficient mechanism that accounts for the mantle flow in planetary interiors, even under very high temperatures and high pressures, and does not produce crystallographic preferred orientation.

To describe plastic flow in crystalline rocks, one commonly considers two distinct potential processes classified as diffusion creep and dislocation creep (Fig. 1). In diffusion creep, plastic strain results directly from the motion of crystal point defects. Vacancy concentration close to a grain boundary under tension, being greater than that close to a grain boundary under compression, leads to a net flux of matter between sources and sinks (Fig. 1A). At high temperatures, vacancies can diffuse through the bulk of the grain as considered in the NH creep mechanism (*4*, *5*) or along the grain boundaries as proposed by Coble (*6*). The grain size defines the characteristic distance between sources and sinks and strongly limits the efficiency of both mechanisms with potential implications on the convection of terrestrial planets (*7*, *8*). Dislocation creep (Fig. 1B), as commonly observed in metallurgy, takes part in crystal recovery processes associated with heat treatments (*9*–*11*). Plastic strain is produced by the glide of a fraction of dislocations that are made free to move owing to two thermally activated mechanisms (*12*), that is, cross-slip (screw dislocations deviating from their initial glide plane) and climb (motion out of the glide planes of nonscrew dislocations after absorbing point defects) (*13*). Contrary to diffusion creep (NH or Coble), dislocation creep can produce crystal preferred orientations potentially strong enough to yield seismic anisotropy (*14*). Characterized by a higher stress sensitivity, dislocation creep is most importantly grain size–insensitive. For this reason, the average grain size is considered the key parameter in determining whether diffusion creep or dislocation creep controls crystal plasticity in the deep Earth. Unfortunately, rocks’ grain size is not known in Earth’s lower mantle. Here, we show that new developments in dislocation dynamics (DD) and creep modeling in high-pressure minerals shed new light on dislocation creep mechanisms and their relevance in mantle conditions.

## RESULTS AND DISCUSSION

First, using DD simulation (*12*, *15*), we investigated the complex interplay between dislocation glide and dislocation climb. These simulations are instructive in understanding how, during creep processes, the glide mobility of dislocations competes with their climb mobility. Figure 2 shows the ratio (*v*_{g}/*v*_{c}) of the glide over the climb velocity as a function of the temperature and the applied stress. Both mobilities strongly depend on the crystal structure and on specific atomic configurations that build the dislocation cores. The climb velocity also depends on the diffusivity of point defects. In olivine, which controls the rheology of upper mantle rocks, glide is always much faster than climb (see Fig. 2A), leading to a creep behavior very close to the one originally proposed by Weertman [see the study of Boioli *et al*. (*12*) and Keralvarma *et al*. (*16*)]. However, it has been already emphasized that the contribution of climb to high-temperature creep of olivine cannot be ignored (*17*). For the high-pressure phases existing in the deeper mantle (wadsleyite, ringwoodite, and bridgmanite), experiments available so far do not provide information on dislocation mobility. However, recent modeling based on atomic-scale computations has successfully yielded dislocation glide velocity (*v*_{g}) models for wadsleyite (*18*), ringwoodite (*19*), and bridgmanite (*20*), which are able to reproduce the rare experimental data available well. This shows that the high stress levels observed experimentally (*1*–*3*) result from the very high lattice friction exhibited by wadsleyite, ringwoodite, and bridgmanite at high pressure. Besides particularities of each crystal structure and of its defects, this general trend reflects the strong influence of pressure on bond strengths, especially the ionocovalent Si-O bond (a fact that is already noticeable in the pressure dependence of elastic constants). It is only at high stresses (see Fig. 2, B and C) that dislocation glide can be activated to produce strain and that standard dislocation creep is expected. However, at lower stresses (which one could expect in the convecting mantle), the situation is drastically different. Below ca. 1.2 GPa for ringwoodite (Fig. 2B) and ca. 1.8 GPa for bridgmanite (Fig. 2C), the dislocation glide velocity is much slower than climb velocity. Hence, the efficiency of dislocation glide as a strain-producing mechanism becomes negligible compared to climb. To the best of our knowledge, there is only one case in materials science where this situation is faced: quasicrystals. Quasicrystals are model metallic alloys with atomic arrangements that induce symmetries that are long assumed to be forbidden in periodic crystals. In quasicrystals, the dense planes are very corrugated, making shear along these planes necessarily difficult. Thus, lattice friction is very high and plastic shear can only be achieved at high temperatures by edge dislocations moving by pure climb (*21*–*23*). The outcome that climb could be a large strain-producing mechanism (and not only a recovery process) has long been predicted theoretically (*24*). It has also been established experimentally in materials where dislocation glide can be hindered geometrically, that is, for instance, in some hexagonal single crystals compressed perpendicular to the basal plane: Mg (*25*, *26*) and Be (*27*).

High lattice friction in high-pressure phases leads to reconsideration of the processes operating during dislocation creep because climb becomes the dominant strain-producing mechanism. The next question then is whether pure climb is also a dominant creep mechanism for high-pressure silicates under mantle conditions. To address this new question and evaluate the potential contribution of climb to mantle rheology, we calculate pure climb creep rates for bridgmanite and compare them with diffusion creep rates. To do that, the two-dimensional (2D) DD model previously used to investigate creep in olivine is modified here to consider only pure climb. Two orthogonal slip systems, [100](010) and [010](100), are introduced in the simulation. Dislocations are parallel straight lines of pure edge characters, perpendicular to the reference 2D plane of the simulation. Their Burgers vector lies in the reference plane and defines the slip direction. The climb direction is in the reference plane and is orthogonal to both the Burgers vector and the line direction. The force acting on each dislocation depends on the stress fields at the dislocation positions, which results from the action of both the external loading applied on the simulation cell and from the stress fields of all other dislocations. It is given by the so-called Peach-Koehler equation (see detailed description in Materials and Methods). In pure climb creep (*24*), vacancies migrate from edge dislocations with their Burgers vectors roughly parallel to the tensile axis (case 1 in Fig. 3A) to edge dislocations with their Burgers vectors roughly perpendicular to the tensile axis (case 2 in Fig. 3A). This vacancy exchange allows dislocations to move with a climb velocity, which not only depends on the mechanical force acting on the dislocation but also on the “chemical force” arising from the gradient in the vacancy concentration and produce plastic strain. At the steady state, dislocation multiplication is counterbalanced by dislocation density reduction due to annihilation events, that is, the destruction of pairs of dislocations with opposite Burgers vectors coming across each other. We find that the steady state is reached, characterized by an average linear increase of the plastic strain with time, as shown in Fig. 3B and a constant dislocation density. An example of the resulting dislocation microstructure is shown in Fig. 3 (C and D). When comparing the steady-state strain rates obtained from pure climb DD models with diffusion creep rates (Fig. 4), dislocation climb is always a more efficient plastic deformation mechanism (provided that grain size is in excess of 0.1 mm) than diffusion creep under lower mantle conditions.

## CONCLUSION

Here, we show that in bridgmanite, creep dominated by dislocation pure climb is a viable mechanism that is able to promote plastic flow at rates compatible with those expected in the mantle. Pure climb creep involves vacancy diffusion as NH creep does; however, the characteristic diffusion length (between vacancy sources and sinks) is controlled by the dislocation density instead of the grain size. The former being necessarily smaller than the latter, pure climb creep becomes more efficient than NH creep when grains are larger than 0.1 mm. One of the major implications of this study is to show that grain size, which is a major unknown in the mantle, is not a controlling parameter of the mantle viscosity as assumed for many years. Moreover, climb involves no lattice rotation. Hence, although strain is produced by dislocation motion, pure climb creep does not produce crystal preferred orientation. Thus, our model is compatible with the absence of strong seismic anisotropy in Earth’s lower mantle.

## MATERIALS AND METHODS

### Glide and climb velocity *v*_{g} and *v*_{c}

Silicates under high pressure are characterized by high lattice friction. To move, a dislocation must overcome the intrinsic resistance of the lattice (which is also related to the specific atomic arrangements that build the dislocation core). This is quantified by the Peierls potential *V*_{P}. The derivative of *V*_{P} defines the Peierls stress τ_{P}, which is commonly viewed as the critical resolved shear stress at 0 K, and as such, gives the mechanical measure of the lattice friction experienced by dislocations. At finite temperature, dislocation motion over the Peierls potential is assisted by the conjugate effect of stress and thermal activation. During this process, the dislocation does not move as a straight line but through the nucleation and propagation of kink pairs. The kink-pair nucleation process corresponds to a small segment of the dislocation line that bulges over the Peierls potential. The further propagation along the dislocation line of the kinks is responsible for the glide of the whole dislocation into the next stable position in the crystal lattice. The kink-pair nucleation process is usually associated with a critical change in enthalpy that has to be supplied by thermal activation under a given stress. Following the kink-pair mechanism, the resulting glide velocity of a dislocation undergoing a uniform resolved shear stress τ* is described with the following formulation (Eq. 1)(1)where *k*_{B} is the Boltzmann constant and *T* is the temperature. Δ*H*_{0}, parameters *p* and *q*, and pre-exponential term *v*_{0} are dislocation intrinsic quantities describing the glide velocities according to a kink-pair mechanism of the dislocation. The theoretical description of dislocation motion involving the kink-pair mechanism has been recently successfully applied to the understanding of elementary deformation processes in ringwoodite (*19*) and bridgmanite (*20*). Glide velocities for the easiest slip systems, that is, <110>{110} dislocations in ringwoodite and [100](010) in bridgmanite, were thus computed in this study using the set of parameters listed in table S1. It should be recalled here that all these parameters come from atomistic calculations based on either first-principles calculations (ringwoodite) or empirical potential parameterizations (bridgmanite). In the case of bridgmanite, the relevance of the empirical parameterization to large (nonelastic) displacements was further validated with a comparison to ab initio calculations of the so-called generalized stacking fault energies at each pressure of interest in the present study. Moreover, as shown by Ritterbex *et al*. (*19*) or Kraych *et al*. (*20*), it should be pointed out that atomistic calculations led to glide velocities for dislocations, which correspond to resolved shear stresses that agree with flow stresses extracted from recent experimental measurements. In olivine, the parameters associated with the glide mobility were obtained by fitting experimental data at ambient pressure from various sources [see for instance the work of Durinck *et al*. (*28*) and references therein]. In table S1, we report the parameters used to describe the velocity of [100] dislocations, which have been found to govern the plastic behavior at high temperatures (*T* > 1300 K).

With the usual assumption that climb is controlled by vacancy diffusion and that the dislocation line is saturated with jogs (steps along the dislocation line that act as sinks/sources of vacancies), the climb velocity under steady-state conditions is given by the following analytical expression (*29*)(2)where *D*^{sd} is the self-diffusion coefficient, Ω is the vacancy formation volume, and η is a geometrical factor that depends on the geometry of the flux field. *X*_{∞} is the vacancy concentration far from the dislocation, whereas *X*_{v} is the equilibrium vacancy concentration in a bulk at a given temperature *T*.

According to Eq. 2, climb velocities (as plotted in Fig. 2) were computed for the three minerals as a function of τ*_{c}, the effective climb stress, at the equilibrium bulk vacancy concentration (*X*_{v} = *X*_{∞}). Relevant parameters corresponding to vacancy self-diffusion coefficient are given in table S2. To explicitly introduce the dependence of the diffusion coefficient on the vacancy concentration, it is useful to rewrite it as follows(3)where *D*^{v} is the vacancy diffusion coefficient, given by the product of the attempt jump frequency, the square of the jump distance, and the jump probability. The latter is given by the exponential term , where the barrier height Δ*H*^{m} is the vacancy migration enthalpy. From diffusion experiments, it is usually possible to extract two distinct parameters: the exponential prefactor of the self-diffusion coefficient *D*^{sd,0} and the vacancy migration enthalpy Δ*H*^{m}. The climb velocities plotted in Fig. 2 were obtained by inserting the experimental values of *D*^{sd,0} and Δ*H*^{m} for the slowest diffusing species in Eq. 3 for the three minerals from different sources in the literature. The vacancy formation volume was calculated from the unit cell volume.

### Climb creep modeling using 2.5D DD

2D DD simulations were performed to model creep produced by climb in MgSiO_{3} bridgmanite. Two orthogonal slip systems, [100](010) and [010](100), were considered. The Burgers vector amplitude *b* was set to the value of the [100] lattice parameter (*b* = 4.65 Å). Within this model, the dislocations were materialized as parallel straight lines of pure edge characters, perpendicular to the reference 2D plane of the simulation. Their Burgers vector * b* lies in the reference plane and defines the slip direction. The climb direction

*is then identified by the direction in the reference plane orthogonal to both the Burgers vector and the line direction*

**n****(**

*l***). The positive orientation of**

*n*=*l*×*b***is taken along the vacancy emission direction, that is, from the dislocation core, the direction pointing to the extra half-plane that characterizes an edge dislocation. By analogy with NH loading conditions (**

*n**4*,

*5*), a constant loading was applied with a tensile stress along the [010] direction and a compressive stress along the [100] direction. Under these conditions, maximum resolved stresses on climb planes were ensured with [010] dislocations climbing by emitting vacancies, while [100] dislocations climb by absorbing them (see Fig. 3A). As a result, the average vacancy concentration in the bulk remained constant because the excess of vacancies created by one slip system was absorbed by the other slip system. The initial vacancy concentration

*X*

_{v}was assumed to be homogeneous. Initial dislocation microstructures were built for different values of dislocation density (from 10

^{8}to 10

^{13}m

^{−2}) by varying the size

*L*of the simulation box from 3 to 1000 μm, keeping the initial dislocation number equal to 100. The simulation area corresponds therefore to a square of size

*L*on which periodic boundary conditions were applied in the reference plane.

Once an initial random microstructure was set, the five following steps were iteratively repeated during the course of the simulation.

Force calculations. The force acting on each dislocation *i* depends on the stress fields at the dislocation position *x*_{i} and it is given by the so-called Peach-Koehler equation(4)where *b*_{i} and *l*_{i} are the Burgers vector and line direction of the *i*th dislocation, respectively, **σ**_{app} is the external applied load, and **σ**_{int} is the internal stress field induced by the dislocation microstructure. The component of the force contributing to climb, along the climb direction *n*_{i,} is then given by(5)

Velocity calculation and displacement predictions. According to Eq. 2, the climb velocity (Eq. 6) of each dislocation *i* not only depends on the mechanical force acting on the dislocation itself but also on the chemical force arising from the gradient in the vacancy concentration. Here, the effective stress in the climb direction is equivalent to the component of the Peach-Koehler force contributing to climb *F*_{c}_{,i} divided by the amplitude of the Burgers vector *b*(6)

Again, *X*_{v} is the equilibrium vacancy concentration and *X*_{∞} is the vacancy concentration at a distance *R* from the dislocation, which we allowed to vary in time and space during the simulations. *R* is taken as the half of the average dislocation distance [], where ρ is the instantaneous dislocation density and *r*_{c} is the core radius that we took equal to 5*b*.

To explicitly introduce the dependence of the climb velocity on the vacancy concentration *X*_{v}, we calculated *X*_{v}*D*^{v} by using Eq. 3, where [following the work of Ammann *et al*. (*30*)] we put the attempt frequency equal to the Debye frequency (*ν*_{D} = 10^{13} Hz), the average jump distance *ℓ* equal to the nearest-neighbor distance in bridgmanite (*ℓ* = 2.5 Å) and *ΔH*^{m} equal to the migration enthalpy of Mg measured in bridgmanite at *P* = 24 GPa and *T* = 1973 to 2273 K (*31*).

To evaluate *X*_{∞}, we divided the simulation box in smaller boxes, as sketched in fig. S1. *X*_{∞} is taken as the average vacancy concentration in the box *j*, where the dislocation is located, and in the first layer of the neighboring boxes (colored area in fig. S1). The displacements were calculated by using the explicit Euler forward algorithm, which is the standard integration method in DD, so that the dislocation position can be written as(7)where *dt* is the simulation time step.

Update of the local and average vacancy concentration. Each dislocation represents a source/sink of vacancies because it needs to emit or absorb vacancies to climb. Here, we considered the idea that each dislocation exchanges vacancies within the dislocations inside the same box. In particular, the rate of change of the local vacancy concentration in the box *j* is given by(8)leading to(9)where the summation is performed over the *N*_{j} dislocations belonging to the box *j*; *L*_{jx} and *L*_{jy} are the linear dimensions of the box *j*. Similarly, the rate of change of the average vacancy concentration is given by(10)and(11)

where *N* is the total number of dislocations and *L*_{x} and *L*_{y} are simulation box dimensions. We noticed that, because we defined the positive climb direction **n**_{i} as the direction of vacancy emission for the dislocation *i*, *v*_{c,i} is positive when the dislocation emits vacancies. This resulted in an increase of average vacancy concentration and on the *j*th box vacancy concentration *X*_{j}. On the contrary, when *v*_{c,i} is negative, the dislocation absorbs vacancies leading to a decrease in and *X*_{j}.

Plastic strain calculations. Starting from the dislocation displacements *dx*_{i} *= v*_{c,i} *dt*, we can evaluate the contribution of climb to the plastic strain tensor ε. Each dislocation moving by climb during the time step *dt* produced an increment of the plastic strain: *d*β_{i} *= b dx*_{i}/*L*_{x}*L*_{y} *= b v*_{c,i} *dt*/*L*_{x}*L*_{y}, so that the climb plastic strain rate produced by the *i*th dislocation is *= b v*_{c,i} *dt*/*L*_{x}*L*_{y}. From the increment of the plastic strain, we can calculate climb contribution to the plastic strain rate and the increments of the plastic strain tensor ε per time step(12)

In general, for each slip system, we can define the climb projection tensor (*Q*^{α} *=* *m*^{α}⊗*m*^{α}), where *m*^{α} is the slip direction for the slip system α. Because in our model we assumed that all the dislocations are edge in character, the slip direction *m*^{α} for the slip system α coincides with the Burgers vector direction *b*^{α}.

Test for dislocation multiplication and annihilation. DD simulations make use of local or constitutive rules to deal with dislocation reactions. Here, we included the possibility for dislocation with the same Burgers vector but opposite sign to annihilate when they are at a distance smaller than a critical radius *r*_{annihil} = 20*b*. This distance is much smaller than the linear box size *L* = *L*_{x} *= L*_{y} that we considered in our simulations: 6 × 10^{3} *b* < *L* < 2 × 10^{6} *b.* Additional rules were included in the 2.5D DD model to reproduce relevant 3D dislocation properties. In particular, a multiplication rule was used to reproduce the general 3D observation that the dislocation density ρ increases linearly with the plastic strain ε: *d*ρ/*d*ε *= m*. In analogy with the values adopted in previously published simulations, we imposed *m* = 2 × 10^{15} (*12*).

### Diffusion creep: NH and Coble

In Fig. 4, strain rates corresponding to NH (*4*, *5*) and Coble creep (*6*) were calculated. Diffusion creep is the result of plastic strain produced by the motion of point defects. When a deviatoric stress is applied to a polycrystalline material, a heterogeneous stress state is built in the material. The vacancy concentration at grain boundaries under tension is larger than that at grain boundaries under compression, producing a flux of vacancies through the grains. This mechanism is referred to as NH creep. Under this loading condition (see sketch in Fig. 3A), the strain rate induced by migration of point defect through the bulk can be expressed by the following equation(13)

NH creep takes into account mass transport through the bulk. Migration of point defects may also occur at grain boundaries. When migration at the interfaces between grains becomes the most effective diffusion path, the strain rate produced can be expressed by the following equation(14)

This type of creep mechanism is referred to as Coble creep. In Eqs. 13 and 14, σ is the applied stress, Ω is the vacancy formation volume of bridgmanite, *k*_{B} is the Boltzmann constant, *T* is the temperature, α is a geometrical factor, here assumed to be equal to (*6*), and *d* is the average grain size. Diffusion strain rates reported in Fig. 4 were calculated by substituting *d* = 0.1 mm and *d* = 10 mm in Eq. 13 and Eq. 14, respectively. We notice that the grain size dependence of the strain rate is 1/*d*^{3} and 1/*d*^{2} for Coble and NH creep, respectively. Thus, Coble creep is expected to dominate at low grain size, but it is less favorable than NH creep at larger grain size.

The vacancy diffusion coefficient *D*^{v} used to calculate the NH creep rate (Eq. 13) was taken equal to the value used to calculate the climb velocity (Eq. 6). To compute the Coble creep strain rate, we wrote the grain boundary diffusion coefficient as , where we consider *ν*_{a} = *ν*_{D} = 10^{13} Hz and *ℓ* equal to the nearest-neighbor distance in bridgmanite *ℓ* = 2.5 Å (as in Eq. 13), where Δ*H*^{m,gb} = 3.22 eV is the migration enthalpy at the grain boundary measured in bridgmanite at *P* = 25 GPa and *T* = 1673 to 2073 K (*32*) and δ is the effective grain boundary thickness taken equal to 0.1 nm.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/3/e1601958/DC1

fig. S1. Sketch of the simulation box.

table S1. Parameters used to compute the glide mobility laws for dislocation in olivine, ringwoodite, and bridgmanite.

table S2. Parameters used to compute the climb mobility laws for dislocation in olivine, ringwoodite, and bridgmanite.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**We acknowledge financial support from the European Research Council under the Seventh Framework Programme (FP7 ERC grant no. 290424–RheoMan).

**Author contributions:**The author list is set in alphabetic order. P. Cordier designed the study and supervised it with P. Carrez. F.B. modelled the creep with B.D. S.R. performed calculations on ringwoodite with K.G. and P. Carrez. A.K. modelled the dislocation glide in bridgmanite with P. Carrez and P.H. P. Cordier wrote the paper with feedback and contributions from all co-authors. All authors discussed and interpreted the results.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017, The Authors