On the role of solute drag in reconciling laboratory and natural constraints on olivine grain growth kinetics
Jean Furstoss, Carole Petit, Andrea Tommasi, Clément Ganino, Daniel Pino Muñoz, Marc Bernacki
OOn the role of solute drag in reconcilinglaboratory and natural constraints on olivinegrain growth kinetics
J. Furstoss , C. Petit , A. Tommasi , C. Ganino , D. PinoMuñoz , and M. Bernacki Université Nice Côte d’Azur, CNRS, OCA, IRD, Géozur, France MINES ParisTech, PSL Research University, CEMEF-Centre demise en forme des matériaux, CNRS UMR 7635, France Géosciences Montpellier, Univ. Montpellier, CNRS, Montpellier,FranceJuly 10, 2020
Abstract
We investigate the effect of solute drag on grain growth (GG) kinetics in olivine-rich rocks through full field and mean field modelling. Considering a drag forceexerted by impurities on grain boundary migration allows reconciling laboratoryand natural constraints on olivine GG kinetics. Solute drag is implemented ina full field level-set framework and on a mean field model, which explicitlyaccounts for a grain size distribution. After calibration of the mean field modelon full field results, both models are able to both reproduce laboratory GGkinetics and to predict grain sizes consistent with observations in peridotitexenoliths from different geological contexts.
Olivine is the major constituent of the Earth upper mantle, and its grain growth(GG) kinetics is of major importance in several geodynamic processes. In fact,a variation in the mean grain size of a mantle rock may drastically change itsmechanical behavior through the grain size dependence of the diffusion creepregime [1]. In this regime, grain size reduction produces significant weakening,whereas grain growth may lead the rocks into the grain-size independent, butstress-dependent dislocation creep regime. A variation in the mean grain sizeof olivine rocks may therefore produce marked changes in the upper mantle1 a r X i v : . [ phy s i c s . g e o - ph ] J u l heology, which may control strain localization [2]. Consistently, preservationof small grain sizes (mylonitic or ultramylonitic microstructures within ductileshear zones [3]) has been proposed as one of the mechanisms allowing to preserveweak plate boundaries over geological times [4].However, numerical models [5, 6] based on experimental data on olivineGG [7, 8, 9] cannot explain the persistence of small grain sizes through millionyears in pure olivine rocks. Although they simulate well the experimental data,these models also fail to predict the commonly observed plurimillimetric grainsizes observed in dunites when run on geologically relevant timescales (i.e., overmillions of years), for which they predict meter scale grain sizes [5, 6]. Toobtain a GG kinetics compatible with natural observations in upper mantlerocks, the presence of second phases is often considered [9, 10]. Nevertheless,to preserve small olivine grain sizes over millions of years, these models have tointroduce some questionable features, such as very slow second phase GG [11],an enhanced mixing between phases [12] or the presence of small fixed particles[10] to increase the impediment of olivine grain boundary migration (GBM).Our understanding of olivine GG in itself is limited by the inconsistencybetween laboratory and natural time and spatial scales. All GG models arecalibrated on experimental data obtained for ultra-fine grains ( to µm ), forwhich the grain size evolution is detectable at the experimental timescale (afew hours to several days), and then extrapolated to natural conditions. As aconsequence, if a physical mechanism becomes prominent only for grain sizeslarger than µm , its effect will not be captured in the experiments. Such aphysical mechanism may be the drag force exerted by impurities also calledsolute drag, which has been very seldom considered in studying GG in rocks(orthopyroxene [13], halite [14]) and never accounted for in modeling GG ofolivine.In the present work we tested the effect of solute drag on the GG kinetics ofolivine. We performed 2D full field GG simulations using a level-set (LS) frame-work [15] in which we implemented the solute drag effect. These models showthat considering this mechanism allows for consistent simulation of olivine GGkinetics at both experimental and upper mantle conditions. Finally we adjusta mean-field model on results of the full field simulations to propose an ana-lytical expression allowing to compute efficiently the grain size evolution withingeodynamic large scale models accounting for grain size-dependent rheologies. Impurities present within the crystal lattice and segregated at grain boundariescan have an impact on their migration kinetics through the so called solute drageffect. The grain boundary migration velocity ( v ) is generally expressed by [16]: v = M P, (1)where M ( m .J − .s − ) is the grain boundary mobility and P is the sum of2he pressures exerted on the grain boundary. The solute drag effect can bedescribed in terms of dragging pressure exerted on the grain boundary, which is afunction of the grain boundary velocity and impurities concentration and nature.The quantification of this drag pressure has been studied theoretically [17] .Its intensity follows three main regimes (high, intermediate and low velocity),similarly to well known dislocation impediment by Cottrell atmospheres [18]. Infact, impurities segregated around grain boundary interact with it and when thegrain boundary migrates the impurity cloud tends to accompany it. In the highvelocity regime, the grain boundary moves so fast that the segregated impuritiescannot follow the interface and the interactions between the impurities and thegrain boundary are much reduced. In the low velocity regime, the impuritycloud can stay segregated around the grain boundary, but the intrinsic dragof the grain boundary (due to drag exerted by the intrinsic defects within theinterface) is generally higher than the drag exerted by the impurities (Fig.1).Between these two velocities, a third regime exists where the impact of solutedrag on GBM is the most important.Even if the drag effect is expected to follow different mathematical relationshipsdepending on velocity regime [17], an unified expression for the drag pressure P i exerted by a c concentration of impurities (within the grain matrix) describingthe two velocity regimes has been proposed by [19] : P i = αvc β v , (2)where v is the grain boundary velocity, c the concentration of impuritieswithin the bulk, α ( J.s.m − ) and β ( s.m − ) are two parameters modulating theintensity of the drag. α controls the intensity of the solute drag pressure and β constrains the grain boundary velocity which is the most impacted by thepresence of the impurities (Fig.1). The drag effect is most effective for grainboundary velocities close to of /β . 3igure 1: Solute drag pressure intensity as a function of the grain boundaryvelocity calculated using eq.2 for c = 100 ppm and different values of α and β .The α and β parameters can be expressed as integrals of functions D ( x ) and E ( x ) representing the variation of the impurity diffusion coefficient and inter-action energy, respectively, along the grain boundary normal ( x ). As these twofunctions are difficult to constrain experimentally, their mathematical expres-sions are generally hypothesized in an integrable way, which permits to obtainanalytical expressions for α and β [19] : α = N v ( k b T ) E D ( sinh ( E k b T ) − E k b T ) , (3)and, β = αk b T δ N v E D , (4)where N v ( m − ) is the number of atoms per unit volum, k b the Boltzmannconstant, E ( J ) the interaction energy, D ( m .s − ) the diffusion coefficient ofthe impurity within the grain matrix and δ ( m ) is the characteristic segregationlength of the impurities around the grain boundary.4 .1 Semi-explicit implementation of solute drag withinthe level-set framework The LS framework [20], already used to model olivine GG [5, 10], proposesan implicit description of the polycrystal through the use of LS functions rep-resenting the signed distance function to the grain boundaries surrounding thegrain they represent (positive inside the grain and negative elsewhere) in a finiteelement (FE) context. The microstructural evolution is simulated by movingLS functions according to physical laws describing GBM [21] or by creatingLS functions to represent new grains nucleated during the recrystallization pro-cesses [22]. Theoretically, each grain of a polycrystal is represented by its ownLS function. In order to reduce the computation time and memory storage, sev-eral non-neighboring grains in the initial microstructure can be grouped to formGlobal Level Set (GLS) functions thanks to a graph coloration technique. Re-coloration technique is used to avoid numerical grain coalescence during grainboundary motion [23, 24]. The GLS functions displacement is computed withinan efficient FE [23, 25] framework using anisotropic mesh refinement around in-terfaces [26]. The initial microstructure is generated using a Voronoï-LaguerreDense Sphere Packing algorithm [27] allowing to respect precisely an imposedinitial grain size distribution.If GG is solely controlled by capillarity (reduction in the grain boundary energy)and solute drag, the grain boundary velocity can be expressed as the sum of thecapillarity pressure ( P capi ) and impurity drag pressure ( P i described by Eq.2)as : (cid:126)v = M ( P capi − P i ) (cid:126)n = M ( − γκ − αc (cid:126)v · (cid:126)n β v ) (cid:126)n, (5)where v = || (cid:126)v || , M ( m .J − .s − ) and γ ( J.m − ) are the grain boundary mobil-ity and energy respectively, κ is the grain boundary local curvature (in 2D) orthe sum of the main local curvatures (in 3D) and (cid:126)n the outward unit normal tothe grain boundary.By introducing the GLS function Φ i to represent the grains i , we can do thefollowing implicit and explicit first order time discretizations, for the temporalderivative of Φ i : ∂ Φ i ∂t = Φ t +∆ t − Φ t ∆ t , (6)for the GLS function velocity : (cid:126)v = Φ t +∆ t − Φ t ∆ t (cid:126)n = − Φ t +∆ t − Φ t ∆ t (cid:126) ∇ Φ , (7)and for the velocity squared : (cid:126)v = ( Φ t − Φ t − ∆ t ∆ t ) = v old . (8)5he displacement of the GLS functions are then computed using the convectiveLS equation [28] : ∂ Φ i ∂t + (cid:126)v · (cid:126) ∇ Φ i = 0 . (9)Considering the geometrical properties of LS function κ = − ∆Φ , (cid:126) ∇ Φ .(cid:126) ∇ Φ = 1 , (cid:126) ∇ Φ = − (cid:126)n , using the above time discretization and substituting Eq.5 within theabove convective LS equation, we obtain the FE strong formulation : M Φ t +∆ t ∆ t − M γ ∆Φ t +∆ t = M Φ t ∆ t , (10)where : M = 1 + M αc β v old . (11)The fact that Eq.11 uses v old instead of v is a simplification that reduces thenon-linear behaviour of the problem. In practical terms, the time marchingscheme uses a small timestep, thus the error of replacing v by v old is small.However it greatly simplify the numerical scheme allowing its implementationon a generic FE code. This formulation is similar to the diffusion formulationclassically used for capillarity driven GG in the LS formalism [21]. By analogywith the heat equation, M is equivalent to a mass term (i.e. the product betweenthe specific heat capacity and density). The main difference with the classic LSapproach for capillarity-driven GG when solute drag is modelled is that themass term differs from one and depends on the velocity of the LS function attime t − ∆ t . This heterogeneous mass term is computed on each node of themesh and linearly interpolated.A major drawback of the LS approach lies in the fact that during grain boundarymigration, the GLS are no longer distance functions || (cid:126) ∇ φ || (cid:54) = 1 . This is par-ticularly problematic when a remeshing technique depending on the distanceproperty is used at grain interfaces. In addition, the new diffusive formulationproposed in Eq.10 requires a distance function as it is based on the respect of || (cid:126) ∇ φ || = 1 , at least in a thin layer around the interface. For these reasons, theGLS functions need to be reinitialized at each time step in order to restore theirmetric property. Numerous approaches exist for this reinitialization procedure.Here we use a new direct fast and accurate approach usable in unstructuredFE mesh proposed by Shakoor et al. [25]. The residual errors inherent to thisapproach are discussed in [29].The introduction of the solute drag pressure is expected to reduce the GGkinetics. Thus its influence has to be accounted for in the adaptative timestepping scheme in order to allow larger steps when the grain size evolves slowly.The timestep is computed by imposing a maximal incremental displacementcorresponding to a given fraction ( H ) of the LS reinitialized width ( E p ) : ∆ t = HE p v m , (12)6here v m is the grain boundary mean velocity : v m = M ( γ ¯ R − c α | ¯˙ R | β ¯˙ R ) , (13)where ¯ R and ¯˙ R are the mean grain radius and its temporal evolution respectively.In order to exclude a negative timestep value, we impose a lower bound for v m to M γ/R max where R max is the maximum grain radius.The explicit time discretization in Eq.8 may have an impact on the numericalresolution if the non-linearity of the grain boundary velocity (Eq.5) is strong.This has been evaluated by performing computations with different timesteps(Fig.2), which show that mean grain size evolution does not depend on H ,neither for the reference case without solute drag ( M = 1 ), nor for the casewith solute drag.Figure 2: Full field model of the mean grain radius evolution using M = 2 . · − mm .J − .s − , γ = 10 − J.mm − for a case without impurities (dashedlines) and for a case with c = 1000 ppm , α = 10 J.s.mm − and β = 10 s.mm − (solid lines). The H coefficient controls the timestep (Eq.12).The GG kinetics when solute drag accounts for impurities exhibit the expectedtrend, being slower than in cases without drag (Fig.2). Thus, the strong formu-lation (Eq.10) for the displacement of LS functions accounting for solute dragcan be validated. 7 .2 Mean field approach for GG with solute drag To construct a mean field model describing GG kinetics accounting for solutedrag, it is important to account for the initial grain size distribution. In fact,computing GG kinetics accounting only for mean grain size evolution will hidethe dispersion of the solute drag pressures exerted within the microstructure,due to variations in the capillarity force, which is a function of grain size. Thiswill result in three clear stages of grain size evolution, an initial mean grain sizeevolution unimpacted by solute drag (high velocity regime), a second phase, themost impacted by solute drag, where mean grain size will be quasi-static and afinal phase also unimpacted by solute drag (low velocity regime). To accountfor the coexistence of these three phases within the microstructure and predictrealistic mean grain size evolutions, the initial dispersion of individual grainsizes and their evolution rates have to be considered.For this purpose, we adapted the Hillert’s model [30], which proposes a discreterepresentation of the grain size distribution (GSD) in the microstructure. Byconsidering R i the radius of the i-grain size bin, the Hillert’s model allowscomputing the evolution of each bin by accounting for the capillarity pressureexpressed as M γ ( RR − R i ) multiplied by a coefficient ( = κ R R with κ ≈ . ,[31]), which is adjusted based on experimental and/or full field results. Thisprocedure enables to follow the evolution of each bin of the initial GSD. Thismean field approach is able to reproduce, in terms of GSD, the predictions offull field simulations in context of pure GG (in 2D [21] and in 3D [31]) even forcomplex initial GSD (like bimodal ones).To account for solute drag, we subtract from the capillarity pressure the solutedrag pressure (Eq.2) replacing v by ˙ R i and v by ˙ R i,old (i.e. the grain sizeevolution rate at the precedent increment), which gives : ˙ R i = κ R R M ( γ ( RR − R i ) − αc ˙ R i β ˙ R i,old ) , (14)which can be reformulated as : ˙ R i = γ ( RR − R i ) / ( R M κR + αc β ˙ R i,old ) , (15)To account for topological effects or non-uniform mass term along grain bound-aries, equation 14 is generalized through the following expression : ˙ R i = γ ( RR − R i ) / ( R M κR + C α αc ˙ R ni,old C β β ) ˙ R i,old ) , (16)where C α , C β and n are mean field parameters which have to be calibrated onfull field simulations. It can be noticed that the equation 16 is equivalent to thenon-generalized form (Eq.15) for n = 0 and C α = C β = 1 .In practice, each bin contains an unique grain of radius R i . We begin with a listof grain radii generated from an imposed initial GSD and compute iteratively8heir evolutions. If a grain radius becomes lower than . µm it is consideredas consumed by the growth of neighboring grains and removed from the grainlist. As this mean field model is intended to be used for long term calculations(Myr), the initial number of bins needed to preserve a representative numberof grains after long annealing time is very large and the computational costbecomes prohibitive. Thus, we use a repopulation strategy allowing to do thecalculation with a reasonable number of bins all along the simulation. To do so,when the number of bins is less than a minimal number, we repopulate the binlist by adding ten new bins for each existing bin R i with radii ranging between . R i and . R i with a step of . R i . Those range of values have been chosenin order minimize the difference between GSD before and after repopulationstep (see an example in Fig.3). The minimal number of bins is fixed at 40 basedon a convergence study.Figure 3: Example of GSD (in number) before and after a repopulation stepwith a minimal bin number of 40. Both histograms are constructed using 20classes ranging between the radii of the bins with the minimal and maximalgrain radius. In mantle rocks, several incompatible elements and impurities are present inconcentrations ranging from few to hundreds ppm [32]. Some of them are moreor less homogeneously dispersed through the bulk material, but others can beenriched at grain boundaries [33]. The elements that can influence grain growththrough solute drag are those that exhibit a partitioning between grain interiors9nd boundaries. In olivine-rich rocks, nickel (Ni), aluminum (Al) and calcium(Ca) [33, 34, 35] display a strong partitioning between olivine grain matrix andboundaries, the latter being often qualified as «incompatible element reservoirs». The major parameters controlling their effect on GBM are the element con-centration, diffusion coefficient and interaction energy with grain boundary. Thelatter will control at first order the intensity of the drag pressure exerted by theimpurity. As interstitials, vacancies and even grain boundaries could have a non-null electrical charge in materials as olivine, the full formulation of this energyterm should include, in addition to an elastic part, an electrostatic component.However for sake of simplicity we will assume a pure elastic interaction corre-sponding to the lattice distortion due to the misfit between the impurity sizeand the size of the typical host ion the impurity replaces. Close to and withinolivine grain boundaries the impurities may replace Mg or Si ions (a depletionin Mg is observed at grain boundaries) [35]. Within the grain interiors, Mgions are located in the octahedral M sites, and the interaction energy betweenimpurity and grain boundary can be computed, as a first order approximation,using the characteristic length r of those sites [35] : E = 4 πE α ( r r i − r ) + 13 ( r i − r ) ) , (17)with E α is the Young’s modulus of the M lattice site and r i is ionic radius ofthe impurity. E α is close to the bulk Young’s modulus and r can be computedfrom the length of the bonds between M and oxygen sites, and between oxygenand oxygen sites [35]. The temperature dependency of these parameters is smalland we will consider in the following E α = 159 GP a and r = 0 . nm [35].Using the above expression and Eqs.2, 3 and 4 one can compute the drag pres-sure exerted by Ni ( r i = 0 . nm ), Al ( r i = 0 . nm ) and Ca ( r i = 0 . nm )for different grain boundary velocities and for natural characteristic c concen-trations of , and ppm respectively [32]. Considering the diffusioncoefficients of the three elements at K [36, 37], we estimate that the char-acteristic drag pressure exerted by Ca is at least 1 order of magnitude higherthan the ones exerted by Ni and Al. Thus in the following, the only impurityconsidered for solute drag will be Ca.The Ca solute drag parameters α Ca and β Ca can be computed using equations 3and 4 with a segregation length ( δ ) of nm [35] and a calcium diffusion coefficient( D ) taken as Arrhenius’s laws where the reference value and activation energyare equal to D = 3 . . − mm .s − , Q D = 200 kJ.mol − [37] respectively.We obtain α Ca = 3 . J.S.mm − and β Ca = 9 . . s.mm − at 1573K and α Ca = 3 . J.S.mm − and β Ca = 2 . . s.mm − at 1073K. The main goal of this study is to show that laboratory experiments and naturalobservation on olivine GG can be reconciled by accounting for Ca solute drag.10n the following section, we define a priori constraints on GG kinetics providedby natural and experimental data.
The annealing experiments for ultrafine grained natural San Carlos olivine [7]seem appropriate to define the experimental reference, particularly the dry runsat K . This sample has probably as much impurities as a natural mantlerock (because it was synthesized using crushed grains of natural olivine), so thesolute drag model should fit these data. These experiments allow constraininggrain growth kinetics for short times ( ≤ h ) and small grain sizes ( − µm ). Natural constraints on GG kinetics have to be considered cautiously becauseof the large number of uncertainties in the determination of physical conditionsof natural systems (initial size distribution, precise thermal history, pure staticconditions, etc.). Contraints on thermal history, strain and grain size evolutionof peridotite xenoliths can be found in the literature (e.g. [38, 39]), but theserocks are a mixture of olivine, pyroxene and other secondary phases and theirmicrostructural evolution may be controlled by their polymineralic nature [10]in addition to the impurity effect. However, in natural polymineralic peridotiteslike harzburgites or lherzolites, the maximum grain size is dictated by the spac-ing between static, pinning phases like spinels, and is always larger than theiractual GS [40]. An alternative could be to focus on dunites, a coarse-grainedrock, which mineral assemblage is greater than 90% olivine with only minoramounts of secondary minerals. However, in such rocks, the initial crystal sizedistribution is clearly related to the mechanism of their formation, which usu-ally involves extensive interaction with melts percolating the mantle and maybe very different from the GS distribution used as initial conditions for ourmodels [41], [42]. Considering these two limitations, we chose to compare ourmodels results with classical (i.e., polymineralic) peridotites, being aware thatour model only captures a part of the processes that limit the maximum grainsize of natural rocks.Considering only rock samples with textures typical of thermal annealing andaccording to the terminology defined by [43] we selected samples which weredescribed as “coarse”, which is equivalent to “protogranular” [44], and equant(or granular). Porphyroclastic, granuloblastic and tabular textures were notconsidered because they reflect significant rock deformation that was not subse-quently fully annealed by grain growth, or grain growth in presence of melts orfluids. We briefly recall here the geological setting, age and temperature historyof the selected samples.
The Udachnaya (Siberia) kimberlite xenoliths
The Udachnaya kimberlite pipe in Siberia is well known due to its diamondmine and because of the occurrence of megacrystalline olivine ( ≈ The Kaapvaal (South Africa) kimberlite xenoliths
The Archean lithosphere of the Kaapvaal craton in South Africa stabilizedaround 3 Ga ago according to Re-Os isotope studies [51]. Mantle xenolithshave been sampled by magmatism in Late Jurassic to Cretaceous times, i.e.between 180 and 90 Ma [52]. Most peridotite xenoliths consist of coarse (5-8mm) to cm-size grained harzburgites indicating significant annealing posteriorto an early stage of deformation [53]. PT estimates range along a low geothermalgradient with temperatures of − K at depths between 80 and 150 km[54, 55, 53]. In the Jagersfontein pipe, very coarse-grained peridotite xenolithsexhibit grain sizes from 5 to 20 mm with temperature estimates ranging from − K [56]. Olivine CaO content reaches up to 1200 ppm [57]. For ourmodel, we assume that a mean temperature of 1073K was reached rapidly afterthe stabilization of the cratonic lithosphere, and we consider that these rockshave spent 2 to 3 Ga at this temperature before being erupted. The Kerguelen hotspot in the Indian Ocean
The Kerguelen archipelago is part of a Large Igneous Province, the KerguelenPlateau, formed above the Kerguelen plume [58, 59, 60]. Plume-related volcan-ism forming the Kerguelen Islands started around 45 Ma ago and lasted until0.1 Ma ago ([61] and references therein). Ultramafic xenoliths brought at thesurface in the Kerguelen Islands by the plume-related volcanism are harzbur-gites and dunites typical of a depleted mantle which has undergone a largedegree of partial melting ([58] and references therein). Equilibrium PT condi-tions determined on xenoliths close to the crust-mantle boundary are around 1GPa and 1173-1273K [60]. Some protogranular harzburgites have mean grainsizes of 2-10 mm while equigranular dunites have a mean grain size between0.5 and 1 mm [58]. Here again, it is difficult to estimate the annealing time ofthese xenoliths. Based on geochemical and petrological analyses the Kerguelenharzburgites were interpreted as residues from a partial melting episode linkedwith the Kerguelen plume, that were subsequently affected by melt percolationforming the dunites [59, 58]. Peridotite xenoliths have been sampled in lavasdated between 28 and 7 Ma ([58] and references therein) so they might havespent up to 38 Ma at temperatures close to 1173-1273K.12he temperature, rock type and mean grain size for different contexts raiseabove are summarized within table 1.Table 1: Summary of constraints on temperature, grain size and geodynamiccontext for the four selected peridotite samples that were used to estimate ourmodel performance. a [53]; b [56]; c [45]; d [62]; e [63]; f [60]; g [58].The geological contexts presented above will be used in the following to testthe performance of our mean-field solute drag model. In a first step, in orderto calibrate the solute drag parameters on natural constraints for the full-fieldmodel, we will consider that the model should produce grain sizes ranging be-tween 0.5 and 10mm for annealing times ranging between 0.1Ma and 1Ga. Forthe full field models being supposed to fit those natural constraints, we tooka temperature of 1073K because this temperature corresponds to a typical up-per mantle temperature. Moreover, it corresponds to the temperature at whichdiffusion becomes so slow that textural rebalancing under static conditions be-comes negligible. In this section, we present the full field and mean field results concerning thelong term annealing of olivine aggregates. The olivine grain boundary energy( γ ), calcium concentration within grain matrix ( c ) are taken as J.mm − [64]and ppm respectively while the grain boundary mobility ( M ) is taken as anArrhenius’s law where the reference value and activation energy are equal to M = 4 . mm .J − .s − , Q M = 185 kJ.mol − [10] respectively.To perform long term annealing full field simulations, passing from micrometerto millimeter scale grain sizes, we need to define a model chaining strategy. Thefull field simulations begin with approximately 2000 grains respecting an initialgrain size distribution corresponding to the one used in laboratory experiments[7]. When the number of grains within the simulation domain is less than 200,the simulation is stopped and a new set of 2000 grains and a larger domain isgenerated for the next job following the final grain size distribution of the latestrun. The error due to this chained calculation is minimized by sampling veryprecisely the final GSD and by imposing it as the new initial GSD using theVoronoï-Laguerre Dense Sphere Packing algorithm [27].In the following, we first propose an adjustment of the solute drag materialparameters, based on full field simulations, which permits to reconcile laboratoryand natural observations by producing the adequate GG kinetics. Then, the13ean field model presented within section 2.2 is calibrated on full field results.Finally we show an application of this framework within the geological contextspresented in section 4.2. We begin our full field simulations using the solute drag parameters α Ca and β Ca presented in section 3. However at 1573K the computed GG kinetics doesnot match the experimental results [7] (Fig.4). The GG curve quickly deviatesfrom the laboratory data and thus we stopped this simulation at the second run.Figure 4: Olivine aggregates GG kinetics at 1573K and 1073K. Thin solid linesrepresent the usual extrapolations (not accounting for solute drag) at 1573K(red) and 1073K (blue) of the laboratory results [7] (yellow). The bold purpleline present the predictions from full field simulation using calcium solute dragparameters α Ca and β Ca calculated in section 3, while the bold red and bluesolid lines represent full field simulations with α = α Ca and β = 100 β Ca at1573K and 1073K, respectively . The dashed lines represent the GG kineticscomputed using the mean field model with adjusted parameters based on thefull field simulations at 1573K (red) and 1073K (blue). The green rectanglerepresents olivine mean grain sizes and annealing times typical of lithosphericmantle rocks (see section 4.2).We tested different values for the solute drag parameters and found that α = α Ca and β = 100 β Ca produce results consistent with laboratory results14t 1573K and have GG kinetics compatible with natural observations at 1073K(Fig.4). Using those values, the mean grain size evolution begins to deviatefrom the classical extrapolation of laboratory results (which does not accountfor solute drag) for a mean grain size near 30 µm at 1573K and near 50 µm at1073K (Fig.4). The GG kinetics is slowed down by solute drag for mean grainsizes of up to ca. 1mm. For coarser mean grain sizes, GG is weakly influencedby solute drag (low velocity regime in Fig.1). The distribution of the mass term M (see section 2.1, Eq.11) within the simulated microstructure (Fig.5) clearlyshows an increase of its value from the small grains to the large ones.Figure 5: Distribution of mass term M along grain boundaries for the first threeruns at 1573K using α = 10 α Ca and β = 500 β Ca . The first run (left), a . mm aside square, is also represented in the right bottom corner of the second run(center), a . mm aside square, which is also represented at the right bottomcorner of the third run (right), a mm aside square.At first, for ultra-fine grained aggregates as the ones used in laboratory ex-periments, the mass term is nearly equal to 1 within all the microstructure(Fig.5). This first phase correspond to the high velocity regime in which thegrain growth kinetics is nearly unimpacted by the solute drag (Fig.1). Af-terwards, as the grain sizes increase the mass term becomes heterogeneouslydistributed between one and its maximum value depending on the local grainboundary velocities, which are controlled by the local curvature. Within thisregime, the long straight grain boundary segments have an higher mass termand their velocities are near the velocity the most impacted by solute drag.Those segments are thus slowed down by solute drag. They also exert a pinningforce on other grain boundaries, which enhances the deceleration of the graingrowth kinetics. Finally, when grain sizes are sufficiently large, the grain bound-ary velocities become very small and one enters the low velocity regime in whichthe grain growth kinetics is again weakly impacted by the solute drag (Fig.1).The mass term almost reaches its maximum value in the whole microstructure15Fig.5).The analysis of the GG kinetics for simulations with different initial mean grainsizes highlights a first phase of very slow or even null mean grain size evolution(Fig.6).Figure 6: Mean grain size evolution at 1573K for α = 10 α Ca and β = 500 β Ca for the first four runs (without applying the time offset), the full lines representthe full field model and the dashed lines represent the mean field computations.The duration of this phase of slow GG is longer when the initial mean grain sizeis coarse (Fig.6). This initial slow GG kinetics can be explained by analyzingthe GSD (Fig.7). In fact, even when the initial microstructure is composed ofcoarse grains, some grains have to shrink to let the other ones grow. This canbe seen by comparing in Fig.7a the GSD of the initial microstructure and at theend of the initial slow GG phase for the experiment with the coarser initial meangrain size in Fig.6. This comparison highlights an enrichment in small grains(Fig.7a). Before disappearing, the GBM velocities of the shrinking grains willnecessarily pass through a velocity regime highly impacted by the presence ofimpurities. This will affect the grain growth kinetics directly (by slowing downthe shrinking) and indirectly by impeding the movements of the other grainboundaries through pinning mechanisms.From these observations, we infer that initial grain size does not affect GGkinetics apart from delaying the start of the rapid GG stage.16 a) (b) Figure 7: GSD (in number) for the fourth run presented in figure 6 (at 1573Kfor α = 10 α Ca and β = 500 β Ca ), 7a : initial GSD (blue), mean field (orange)and full field (yellow) modeled GSD after s , 7b : mean field (orange) andfull field (yellow) modeled GSD after . . s .After calibration, the mean field model (Eq.16) reproduces well the full fieldmodeled GG kinetics for both temperatures (Figs.4 and 6). The best fittingmean field parameters C α , C β and n are equal to . , . and . respectively.The GSD predicted by the mean field approach is consistent with the full fieldmodeled ones (Fig.7) even if some differences can be seen after some annealingtime (Fig.7b). This mean field approach allows to be predictive on olivineaggregates GG kinetics for a much lower computational cost, which permits totest our formalism on different geological contexts. In order to estimate how well our mean field model predicts the average grain sizeof natural samples with reasonable c (Ca concentration within bulk) values, wetest it against different contexts in terms of annealing temperature and measuredgrain sizes. For the mean field GG model, the temperature is kept constant and grain sizegrows indefinitely; we compute the GG curve for residence times up to 2 Ga.GG curves are computed from four temperatures of 1073, 1173, 1273 and 1373Kand c values of 600, 800, 1000 and 1200 ppm.All models were run initially with an initial grain size of 20 µm and a standarddeviation of 2 µm ; then, isothermal models were run again with an initial grainsize of 0.5mm and a standard deviation of 50 µm and an intial grain size of 2mmand a standard deviation of 200 µm . In fact, the initial grain size has no effect17f the grain growth curve or final grain size, apart from shifting the beginningof the beginning of the positive slope on the grain growth curve (Fig.8).Figure 8: Grain size evolution for isothermal model with T = 1173 K , c =800 ppm and three initial grain sizes of 20, 500 and 2000 µm . The c parameter, in the range of tested values, has only a moderate effectof the final or intermediate grain sizes (table 2 and figure 9). However, thesegrain sizes are considerably smaller than the ones predicted without impurities.Indeed, for isothermal models, grain sizes reach 2 to 8 meters after 1 Ga attemperatures between 1173 and 1373K, respectively, while they do not exceed47 cm for the scenario with the largest temperature (1373K) and the lowestimpurity concentration (600 ppm). 18able 2: Temperature and geochemical ( c ) conditions for mean field modelsand their results in terms of grain sizes at two different steps (10 and 1000 Ma).Mean field models predict a GG rate which decreases exponentially with time(Fig.9). Extremely large grain sizes of several cm can be reached for condi-tions (temperature and annealing times) consistent with those inferred for theUdachnaya kimberlite xenoliths (i.e., temperatures of 1173 to 1273K and resi-dence times between 100 and 1450 Ma). The Kaapvaal peridotites, which havesmaller grain sizes and potentially larger residence times than the Udachnayaones, range slightly below the 1073K temperature curve.19igure 9: Grain size evolution for isothermal models with 4 different temper-atures and 4 different impurity ( c ) concentrations. Solid squares indicate theapproximate location of Siberia and Kerguelen and Kaapvaal xenoliths on thiscurve according to the literature (see references in Tab.1).For annealing times of 1 to 40 Ma, predicted mean grain sizes range between2 mm (for 1173K and c = 1200 ppm) and 4 cm (for 1273K and c = 600ppm). These predictions match the grain sizes of the harzburgites from theKerguelen Islands. The grain size measured in Kerguelen dunites is smallerthan these predictions; it suggests either rapid exhumation (hence very limitedgrain growth) or a lower residence temperature ( ≈ The GG kinetics simulated using β Ca cannot adequately reproduce both labo-ratory experiments and natural observations on grain size evolution of olivineaggregates (Fig.4). However we find that a value of β two orders of magnitudelarger than β Ca allows our models to be consistent with those two constraints.This difference could arise from two main reasons :First, the two functions D ( x ) and E ( x ) needed for the calculation of solute dragparameters α and β (see section 2) are poorly known even for metallic mate-20ials in which the interaction between an impurity ion and a grain boundaryis mostly elastic. For ceramics-like materials such as olivine, this interactionshould also account for electrostatic interaction between the solute ion and thegrain boundary and the interaction between solute-vacancy dipoles and electricalfield around grain boundaries [65]. The quantification of all of those interac-tions should be done by atomistic calculations using systematic approaches fordescribing grain boundaries and their interactions with solute [66].Secondly, it is well accepted that the impurity concentration at grain boundaryevolves with the grain size, increasing when grain size increases [67] which couldalso be expressed by a grain size dependent partition coefficient [35]. Even if thesolute drag pressure (Eq.2) depends on grain matrix impurity concentration,one can find expressions for α and β as functions of the partition coefficient[19]. The expected effect of indirectly introducing this grain boundary impurityconcentration may be similar to an increase of β , as we did here. In fact,increasing the β will increase the space of the high velocity regime for which thegrain boundary migration will be poorly affected by impurities. An increase ofthe impurity concentration with grain size will also result in a lower impact ofthe solute drag for small grains.Our full field and mean field simulations show an initial phase of very slowgrain size evolution (Figs.6 and 8) due to the direct slowed down of the grainshrinkage and growth by solute drag, and indirectly by the impediment of GBMby the impurity slowed grain boundaries. This quasi-static phase, which islonger for coarser initial grain sizes, implies a very weak dependency of thegrain growth kinetics on initial grain size. This small dependency frees us fromthe need of precise constrains on the initial grain size, which is very difficultto know in geological contexts. Taking advantage of this, we can apply oursolute drag model for different stable geological contexts. For the majority ofthe geological contexts presented here, our formalism shows consistent grain size/ time predictions (Fig.9). It is difficult to evaluate precisely the performanceof our models with respect to the evolution of natural samples, given the lack ofdata concerning their temperature and grain size evolution. However, impuritydrag due to Ca concentration in olivine, in the range of commonly observedvalue, explains a grain size reduction of several orders of magnitude at geologicaltimescales, compared to models without impurities. In this work we have demonstrated that accounting for the presence of impuritieswithin olivine rich-rocks permits to explain GG kinetics of both experimentallyand naturally annealed olivine aggregates. The solute drag parameters neededfor this approach are however quite different from the ones expected for thecalcium, which is the mantle rock impurity expected to be the most impact-ing on olivine grain boundary migration. Atomistic calculations or a strongexperimental framework could be needed to explain this gap.We have developed a new mean field model accounting for the presence of those21mpurities and showed this approach could be used to predict grain size ofolivine-rich rocks in geological contexts such as oceanic cooling or isothermalevolution. As this approach successfully reproduced natural grain size in an-nealed peridotite, healing kinetics may be implemented in large scale numericalgeodynamic models based on this framework. In order to have a useful grainsize evolution law for geodynamic model, a model accounting for the compe-tition between grain size reduction due to deformation and growth should bedeveloped.In order to be consistent with the real multiphase nature of natural rocks, theinfluence of second phases, such as pyroxenes and alumina phases, on GG ki-netics should also be considered. Such a framework may open the door to apaleo-chronometer based on grain size evolution within geological contexts inwhich deformation has stopped and GG predominates.
This work was supported by CNRS INSU 2018-programme TelluS-SYSTER.The support of the French Agence Nationale de la Recherche (ANR), Arcelor-Mittal, FRAMATOME, ASCOMETAL, AUBERT&DUVAL, CEA, SAFRANthrough the DIGIMU Industrial Chair and consortium are gratefully acknowl-edged.
The data for supporting all figures of the paper are avalaible upon request tothe authors as well as the mean field code used in this work.
References [1] Shun-Ichiro Karato, Mervyn S Paterson, and John D FitzGerald. Rheologyof synthetic olivine aggregates: influence of grain size and water.
Journalof Geophysical Research: Solid Earth , 91(B8):8151–8176, 1986.[2] Jean Braun, Jean Chéry, Alexei Poliakov, David Mainprice, Alain Vauchez,Andrea Tomassi, and Marc Daignières. A simple parameterization of strainlocalization in the ductile regime due to grain size reduction: A case studyfor olivine.
Journal of Geophysical Research: Solid Earth , 104(B11):25167–25181, 1999.[3] RLM Vissers, MR Drury, EH Hoogerduijn Strating, and D van der Wal.Shear zones in the upper mantle: A case study in an alpine iherzolite massif.
Geology , 19(10):990–993, 1991.[4] David Bercovici and Yanick Ricard. Plate tectonics, damage and inheri-tance.
Nature , 508(7497):513, 2014.225] Jean Furstoss, Marc Bernacki, Clément Ganino, Carole Petit, and DanielPino-Muñoz. 2d and 3d simulation of grain growth in olivine aggregatesusing a full field model based on the level set method.
Physics of the Earthand Planetary Interiors , 283:98–109, 2018.[6] Xu Chu and Jun Korenaga. Olivine rheology, shear stress, and grain growthin the lithospheric mantle: geological constraints from the kaapvaal craton.
Earth and Planetary Science Letters , 333:52–62, 2012.[7] SI Karato. Grain growth kinetics in olivine aggregates.
Tectonophysics ,168(4):255–273, 1989.[8] Tomohiro Ohuchi and Michihiko Nakamura. Grain growth in the forsterite–diopside system.
Physics of the Earth and Planetary Interiors , 160(1):1–21,2007.[9] Takehiko Hiraga, Chihiro Tachibana, Naoki Ohashi, and Satoru Sano.Grain growth systematics for forsterite ± enstatite aggregates: Effect oflithology on grain size in the upper mantle. Earth and Planetary ScienceLetters , 291(1-4):10–20, 2010.[10] Jean Furstoss, Marc Bernacki, Carole Petit, Julien Fausty, Daniel Pino-Muñoz, and Clément Ganino. Full field and mean field modeling of graingrowth in a multiphase material under dry conditions: application to peri-dotites.
Journal of Geophysical Research: Solid Earth , page e53942, 2020.[11] T Nakakoji and T Hiraga. Diffusion creep and grain growth in forsterite+20 vol% enstatite aggregates: 2. their common diffusional mechanism andits consequence for weak-temperature-dependent viscosity.
Journal of Geo-physical Research: Solid Earth , 123(11):9513–9527, 2018.[12] David Bercovici and Yanick Ricard. Mechanisms for the generation of platetectonics by two-phase grain-damage and pinning.
Physics of the Earth andPlanetary Interiors , 202:27–55, 2012.[13] Philip Skemer and Shun-ichiro Karato. Effects of solute segregation on thegrain-growth kinetics of orthopyroxene with implications for the deforma-tion of the upper mantle.
Physics of the Earth and Planetary Interiors ,164(3-4):186–196, 2007.[14] M Guillope and JP Poirier. Dynamic recrystallization during creep ofsingle-crystalline halite: An experimental study.
Journal of GeophysicalResearch: Solid Earth , 84(B10):5557–5567, 1979.[15] Marc Bernacki, Roland E Logé, and Thierry Coupez. Level set frameworkfor the finite-element modelling of recrystallization and grain growth inpolycrystalline materials.
Scripta Materialia , 64(6):525–528, 2011.[16] Frederick John Humphreys and Max Hatherly.
Recrystallization and relatedannealing phenomena . Elsevier, 2012.2317] K Lücke and K Detert. A quantitative theory of grain-boundary motion andrecrystallization in metals in the presence of impurities.
Acta Metallurgica ,5(11):628–637, 1957.[18] Ao H Cottrell and BA Bilby. Dislocation theory of yielding and strainageing of iron.
Proceedings of the Physical Society. Section A , 62(1):49,1949.[19] John W Cahn. The impurity-drag effect in grain boundary motion.
Actametallurgica , 10(9):789–798, 1962.[20] Marc Bernacki, Héba Resk, Thierry Coupez, and Roland E Logé. Finiteelement model of primary recrystallization in polycrystalline aggregatesusing a level set framework.
Modelling and Simulation in Materials Scienceand Engineering , 17(6):064006, 2009.[21] Ana Laura Cruz-Fabiano, R Logé, and Marc Bernacki. Assessment of sim-plified 2d grain growth models from numerical experiments based on a levelset framework.
Computational Materials Science , 92:305–312, 2014.[22] Ludovic Maire, Benjamin Scholtes, Charbel Moussa, Nathalie Bozzolo,Daniel Pino Muñoz, Amico Settefrati, and Marc Bernacki. Modeling of dy-namic and post-dynamic recrystallization by coupling a full field approachto phenomenological laws.
Materials & Design , 133:498–519, 2017.[23] Benjamin Scholtes, Modesar Shakoor, Amico Settefrati, Pierre-OlivierBouchard, Nathalie Bozzolo, and Marc Bernacki. New finite element devel-opments for the full field modeling of microstructural evolutions using thelevel-set method.
Computational Materials Science , 109:388–398, 2015.[24] Benjamin Scholtes, Romain Boulais-Sinou, Amico Settefrati, Daniel PinoMuñoz, Isabelle Poitrault, Aurore Montouchet, Nathalie Bozzolo, and MarcBernacki. 3d level set modeling of static recrystallization considering storedenergy fields.
Computational Materials Science , 122:57–71, 2016.[25] Modesar Shakoor, Benjamin Scholtes, Pierre-Olivier Bouchard, and MarcBernacki. An efficient and parallel level set reinitialization method–application to micromechanics and microstructural evolutions.
AppliedMathematical Modelling , 39(23-24):7291–7302, 2015.[26] Héba Resk, Laurent Delannay, Marc Bernacki, Thierry Coupez, andR Logé. Adaptive mesh refinement and automatic remeshing in crystalplasticity finite element simulations.
Modelling and Simulation in Materi-als Science and Engineering , 17(7):075012, 2009.[27] Karim Hitti, Patrice Laure, Thierry Coupez, Luisa Silva, and MarcBernacki. Precise generation of complex statistical representative volumeelements (rves) in a finite element context.
Computational Materials Sci-ence , 61:224–238, 2012. 2428] Stanley Osher and James A Sethian. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations.
Jour-nal of computational physics , 79(1):12–49, 1988.[29] Sebastian Florez, Modesar Shakoor, Thomas Toulorge, and Marc Bernacki.A new finite element strategy to simulate microstructural evolutions.
Com-putational Materials Science , 172:109335, 2020.[30] M Hillert. On the theory of normal and abnormal grain growth.
Actametallurgica , 13(3):227–238, 1965.[31] Ludovic Maire, Benjamin Scholtes, Charbel Moussa, Nathalie Bozzolo,D Pino Muñoz, and Marc Bernacki. Improvement of 3d mean field modelsfor capillarity-driven grain growth based on full field simulations.
Journalof materials science , 51(24):10970–10981, 2016.[32] Jan CM De Hoog, Louise Gall, and David H Cornell. Trace-element geo-chemistry of mantle olivine and application to mantle petrogenesis andgeothermobarometry.
Chemical Geology , 270(1-4):196–215, 2010.[33] Kazuhiro Suzuki. Grain-boundary enrichment of incompatible elements insome mantle peridotites.
Chemical Geology , 63(3-4):319–334, 1987.[34] Takehiko Hiraga, Ian M Anderson, and David L Kohlstedt. Chemistry ofgrain boundaries in mantle rocks.
American Mineralogist , 88(7):1015–1019,2003.[35] Takehiko Hiraga, Ian M Anderson, and David L Kohlstedt. Grain bound-aries as reservoirs of incompatible elements in the earth’s mantle.
Nature ,427(6976):699, 2004.[36] Carl Spandler, H St C O’Neill, and Vadim S Kamenetsky. Survival times ofanomalous melt inclusions from element diffusion in olivine and chromite.
Nature , 447(7142):303, 2007.[37] LA Coogan, A Hain, S Stahl, and S Chakraborty. Experimental determi-nation of the diffusion coefficient for calcium in olivine between 900 c and1500 c.
Geochimica et Cosmochimica Acta , 69(14):3683–3694, 2005.[38] Shiran Liu, Andréa Tommasi, Alain Vauchez, and Maurizio Mazzucchelli.Deformation, annealing, melt-rock interaction, and seismic properties ofan old domain of the equatorial atlantic lithospheric mantle.
Tectonics ,38(4):1164–1188, 2019.[39] HG Ave Lallemant, JC C Mercier, NL Carter, and JV Ross. Rheology ofthe upper mantle: inferences from peridotite xenoliths.
Tectonophysics ,70(1-2):85–113, 1980. 2540] Marco Herwegh, Jolien Linckens, Andreas Ebert, Alfons Berger, andSH Brodhag. The role of second phases for controlling microstructuralevolution in polymineralic rocks: A review.
Journal of Structural Geology ,33(12):1728–1750, 2011.[41] Emmanuel Berger and Michel Vannier. Les dunites en enclaves dans lesbasaltes alcalins des îles océaniques: approche pétrologique.
Bulletin deminéralogie , 107(5):649–663, 1984.[42] P. B. Kelemen. Reaction between ultramafic rock and fractionating basalticmagma i. phase relations, the origin of calc-alkaline magma series, and theformation of discordant dunite.
Journal of petrology , 31(1):51–98, 1990.[43] Ben Harte. Rock nomenclature with particular relation to deformationand recrystallisation textures in olivine-bearing xenoliths.
The Journal ofGeology , 85(3):279–288, 1977.[44] JC C Mercier and Adolphe Nicolas. Textures and fabrics of upper-mantleperidotites as illustrated by xenoliths from basalts.
Journal of Petrology ,16(1):454–487, 1975.[45] LN Pokhilenko, VG Mal’Kovets, DV Kuz’Min, and NP Pokhilenko.New data on the mineralogy of megacrystalline pyrope peridotite fromthe udachnaya kimberlite pipe, siberian craton, yakutian diamondiferousprovince. In
Doklady Earth Sciences , volume 454, page 179. Springer Sci-ence & Business Media, 2014.[46] GL Dehvis, NV Sobolev, and AD Khar’kiv. New data on the age of yakutiankimberlites obtained by uranium-lead study of zircons.
Doklady AkademiiNauk SSSR , 254(1):175–179, 1980.[47] IP Ilupin, VI Vaganov, and BI Prokopchuk. Kimberlites.
Moscow, Ne-dra.[In Russian.] , 1990.[48] NP Pokhilenko, NV Sobolev, FR Boyd, DG Pearson, and N Shimizu.Megacrystalline pyrope peridotites in the lithosphere of the siberian plat-form: mineralogy, geochemical peculiarities and the problem of their origin.
Russian Geology and Geophysics , 34:56–67, 1993.[49] WL Griffin, FV Kaminsky, CG Ryan, SY O’Reilly, TT Win, and IP Ilupin.Thermal state and composition of the lithospheric mantle beneath the dal-dyn kimberlite field, yakutia.
Tectonophysics , 262(1-4):19–33, 1996.[50] Dmitri A Ionov, Richard W Carlson, Luc S Doucet, Alexander V Golovin,and Oleg B Oleinikov. The age and history of the lithospheric mantle ofthe siberian craton: Re–os and pge study of peridotite xenoliths from theobnazhennaya kimberlite.
Earth and Planetary Science Letters , 428:108–119, 2015. 2651] DG Pearson, SB Shirey, RW Carlson, F Rl Boyd, NP Pokhilenko, andN Shimizu. Re/os, sm/nd, and rb/sr isotope evidence for thick archaeanlithospheric mantle beneath the siberian craton modified by multistagemetasomatism.
Geochimica et Cosmochimica Acta , 59(5):959–977, 1995.[52] WL Griffin, JM Batumike, Y Greau, NJ Pearson, SR Shee, and Suzanne YO’Reilly. Emplacement ages and sources of kimberlites and related rocks insouthern africa: U–pb ages and sr–nd isotopes of groundmass perovskite.
Contributions to Mineralogy and Petrology , 168(1):1032, 2014.[53] Virginie Baptiste and Andréa Tommasi. Petrophysical constraints on theseismic properties of the kaapvaal craton mantle root.
Solid Earth , 5(1):45,2014.[54] Xu Chu and Jun Korenaga. Olivine rheology, shear stress, and grain growthin the lithospheric mantle: geological constraints from the kaapvaal craton.
Earth and Planetary Science Letters , 333:52–62, 2012.[55] RL Saltzer, N Chatterjee, and TL Grove. The spatial distribution of garnetsand pyroxenes in mantle peridotites: pressure–temperature history of peri-dotites from the kaapvaal craton.
Journal of Petrology , 42(12):2215–2229,2001.[56] Peter A Winterburn, Ben Harte, and John J Gurney. Peridotite xeno-liths from the jagersfontein kimberlite pipe: I. primary and primary-metasomatic mineralogy.
Geochimica et Cosmochimica Acta , 54(2):329–341, 1990.[57] RL Hervig, JV Smith, and JB Dawson. Lherzolite xenoliths in kimber-lites and basalts: petrogenetic and crystallochemical significance of someminor and trace elements in olivine, pyroxenes, garnet and spinel.
Earthand Environmental Science Transactions of the Royal Society of Edinburgh ,77(3):181–201, 1986.[58] Jerôme Bascou, Guillaune Delpech, A Vauchez, BN Moine, Jean-Yves Cot-tin, and Guilhem Barruol. An integrated study of microstructural, geo-chemical, and seismic properties of the lithospheric mantle above the ker-guelen plume (indian ocean).
Geochemistry, Geophysics, Geosystems , 9(4),2008.[59] Nadine Mattielli, Dominique Weis, M Grégoire, Jean Paul Mennessier,Jean-Yves Cottin, and A Giret. Kerguelen basic and ultrabasic xenoliths:evidence for long-lived kerguelen hotspot activity.
Lithos , 37(2-3):261–280,1996.[60] M Grégoire, Jean-Yves Cottin, Nadine Mattielli, C Nicollet, DominiqueWeis, and A Giret. The kerguelen archipelago: an hypothetic continentalmafic protolith.
Terra Antarctica , 2(1):1ą6, 1995.2761] Jean-Yves Cottin, Gilbert MiChon, and Guillaume DelpeCh. The kerguelenvolcanic plateau: the second largest oceanic igneous province (lip) on earthand a witness of the indian ocean opening.
The Kerguelen Plateau: MarineEcosystems and Fisheries , 2011:29–42, 2011.[62] AG Goncharov, DA Ionov, Luc Serge Doucet, and LN Pokhilenko. Thermalstate, oxygen fugacity and c o h fluid speciation in cratonic lithosphericmantle: New data on peridotite xenoliths from the udachnaya kimberlite,siberia.
Earth and Planetary Science Letters , 357:99–110, 2012.[63] DS Yudin, AA Tomilenko, AV Travin, AM Agashev, NP Pokhilenko, andYu Orihashi. The age of udachnaya-east kimberlite: U/pb and 40 ar/39ar data. In
Doklady Earth Sciences , volume 455, pages 288–290. PleiadesPublishing, 2014.[64] RF Cooper and DL Kohlstedt. Rheology and structure of olivine-basaltpartial melts.
Journal of Geophysical Research: Solid Earth , 91(B9):9315–9323, 1986.[65] MF Yan, RM Cannon, and HK Bowen. Space charge, elastic field, anddipole contributions to equilibrium solute segregation at interfaces.
Journalof Applied Physics , 54(2):764–778, 1983.[66] Pierre Hirel, Gabriel Franck Bouobda Moladje, Philippe Carrez, andPatrick Cordier. Systematic theoretical study of [001] symmetric tilt grainboundaries in mgo from 0 to 120 gpa.
Physics and Chemistry of Minerals ,46(1):37–49, 2019.[67] Katharina Marquardt and Ulrich H Faul. The structure and compositionof olivine grain boundaries: 40 years of studies, status and current devel-opments.