Alper Soyler

Sevil Cikrikci

Cagri Cavdaroglu

Dylan Bouillaud

Jonathan Farjon

Dr Patrick Giraudeau
email: patrick.giraudeau@univ-nantes.fr

Dr Mecit H Oztop

Multi-scale 1H benchtop NMR analysis for milk analysis

Keywords: Benchtop NMR spectroscopy, milk, lactose hydrolysis, classification, flow NMR, artificial neural network (ANN)

Credit Author Statement

-The study is a part of the PhD dissertation of Dr. Alper Soyler -Dr. Farjon and Mr. Bouillaud helped Dr. Soyler in the experiments.

-Dr. Cikrikci helped to arrange the data and Mr. Çavdaroğlu conducted and wrote the statistical analysis.

-Dr. Oztop and Dr. Giraudeau were the advisors of the study.

J o u r n a l P r e -p r o o f

INTRODUCTION

Milk is a unique product including significant amounts of nutrients such as calcium, proteins and vitamins [START_REF] Dutra Rosolen | Lactose Hydrolysis in Milk and Dairy Whey Using Microbial β-Galactosidases[END_REF]. It can be considered as the 'magic liquid' due to its rich and various composition and confirmed benefits on the newborns. In this context, accessible analytical tools are crucial to ensure the quality of milk, both to optimize the production processes and after production to ensure its characteristic properties.

The first focus of this paper is about the production stage, on the example of the potential production of lactose-free milk. One of the milk constituents, lactose (aka milk sugar) is a disaccharide composed of glucose and galactose that are linked by a β-1,4 glycosidic bond [START_REF] Dong | Dispersible Biopolymer Particles Loaded with Lactase as a Potential Delivery System to Control Lactose Hydrolysis in Milk[END_REF]. However, lactose found in the milk or milk products cause digestion problems in some part of the population. The absence of the digestive enzyme lactase which is also known as 'lactose intolerance' is the reason of this problem. Significant amount of world population suffers from lactose intolerance thus demand for lactose free products is getting increasing attention [START_REF] Churakova | Accurate analysis of residual lactose in low-lactose milk: Comparing a variety of analytical techniques[END_REF]. Lactose free milk is produced by the pre-digestion of the lactose in milk with the addition of the lactase enzyme. Controlling the production of lactose-free milk during its production has a major industrial impact. There exist many different techniques recognized by analytical standard agencies to determine the lactose content of milk, such as polarimetry, mid-infrared detection, gravimetry, differential pH, enzymatic methods detecting either the glucose or galactose moiety of lactose and HPLC [START_REF] Churakova | Accurate analysis of residual lactose in low-lactose milk: Comparing a variety of analytical techniques[END_REF]. Each of these methods has its own advantages and disadvantages, but to our knowledge none of them have the ability to monitor the hydrolysis process online in the continuous mode. Therefore, developing new online tools for controlling this reaction is of major importance for the milk industry.

On the other hand, the differentiation of milk products not only regarding the sugar type (lactose free or not) but also according to the other parameters such as process, origin, fat J o u r n a l P r e -p r o o f content, etc. requires considerable attention especially considering the regulatory specifications and quality of the milk. In recent years, several studies have focused on discriminating the products coming from different farming systems, animals or feeding strategies, origin etc. [START_REF] Bergamaschi | Chemometric authentication of farming systems of origin of food (milk and ripened cheese) using infrared spectra, fatty acid profiles, flavor fingerprints, and sensory descriptions[END_REF].

Traditional methods like chromatographic techniques or sensory analyses might require highly skilled operators and they could not be easily adapted to online monitoring. Moreover, they could be time consuming and expensive. Easy, quick to use, minimum chemical requiring methods have always been on the interest of manufacturers for ensuring the quality and authenticity.

Nuclear magnetic resonance (NMR) spectroscopy, is a powerful, fingerprint tool for the simultaneous identification and quantification of compounds in complex mixtures [START_REF] Monakhova | NMR spectroscopy as a screening tool to validate nutrition labeling of milk, lactose-free milk, and milk substitutes based on soy and grains[END_REF]. It is recognized as a valuable tool for food analysis enabling not only the characterization of food matrices in terms of authenticity and quality but also identification of counterfeits [START_REF] Sobolev | Use of NMR applications to tackle future food fraud issues[END_REF]. Also, studies of 1 H NMR-based metabolomics to get metabolite profiles during storage [START_REF] Jansson | Chemical and proteolysis-derived changes during long-term storage of lactose-hydrolyzed ultrahigh-temperature (UHT) milk[END_REF] or to analyze heritability of individual milk metabolites [START_REF] Buitenhuis | Estimation of genetic parameters and detection of quantitative trait loci for metabolites in Danish Holstein milk[END_REF] in animal based products have been explored [START_REF] Bertram | NMR-Based Metabolomics: Quality and Authenticity of Milk and Meat[END_REF]. However, most of the NMR spectroscopy instruments work at high field ( 1 H frequencies >300 MHz) and requires a massive dedicated place with specific installation and trained staff as well as a substantial investment cost and expensive cryogenic fluids for maintenance [START_REF] Bouillaud | Benchtop NMR for the monitoring of bioprocesses[END_REF]. To overcome such limitations, benchtop systems capable of conducting spectroscopy experiments have been developed [START_REF] Perlo | High-Resolution NMR Spectroscopy with a Portable Single-Sided Sensor[END_REF]. These systems have been successfully used to monitor authenticity in some kind of foods such as in coffee [START_REF] Defernez | Low-field 1H NMR spectroscopy for distinguishing between arabica and robusta ground roast coffees[END_REF] and in meat [START_REF] Jakes | Authentication of beef versus horse meat using 60 MHz 1H NMR spectroscopy[END_REF]. Furthermore, thanks to their compact character, NMR spectroscopy can easily be used for online monitoring under flow conditions [START_REF] Anderssen | Online monitoring of enzymatic hydrolysis of marine by-products using benchtop nuclear magnetic resonance spectroscopy[END_REF][START_REF] Bouillaud | Benchtop flow NMR spectroscopy as an online device for the in vivo monitoring of lipid accumulation in microalgae[END_REF][START_REF] Knox | Benchtop flow-NMR for rapid online monitoring of RAFT and free radical polymerisation in batch and continuous reactors[END_REF][START_REF] Rönnols | Lignin analysis with benchtop NMR spectroscopy[END_REF][START_REF] Soyler | Real-time benchtop NMR spectroscopy for the online monitoring of sucrose hydrolysis[END_REF] J o u r n a l P r e -p r o o f Giraudeau, & Oztop, 2020), which makes them ideal candidates for the online control of food production processes. For instance, benchtop NMR spectroscopy has been used to monitor continuous processes such as enzymatic hydrolysis of marine by-products [START_REF] Anderssen | Online monitoring of enzymatic hydrolysis of marine by-products using benchtop nuclear magnetic resonance spectroscopy[END_REF], sucrose hydrolysis [START_REF] Soyler | Real-time benchtop NMR spectroscopy for the online monitoring of sucrose hydrolysis[END_REF], in vivo monitoring of lipid accumulation in microalgae [START_REF] Bouillaud | Using benchtop NMR spectroscopy as an online non-invasive in vivo lipid sensor for microalgae cultivated in photobioreactors[END_REF] and lignin analysis [START_REF] Rönnols | Lignin analysis with benchtop NMR spectroscopy[END_REF].

However it is important to note that implementing NMR spectroscopy in a flow system is challenging and can introduce additional complexity relative to static samples [START_REF] Knox | Benchtop flow-NMR for rapid online monitoring of RAFT and free radical polymerisation in batch and continuous reactors[END_REF]. The impact of flow rate on signal intensities should be investigated carefully and the parameters should be set accordingly [START_REF] Soyler | Real-time benchtop NMR spectroscopy for the online monitoring of sucrose hydrolysis[END_REF].

To our knowledge, this will be 1 st study where milk will be investigated by using a bench top NMR spectroscopy system from different aspects. We have explored the versatility of benchtop NMR on milk samples for two different purposes. In the first part of the study, benchtop NMR spectroscopy was utilized in conditions mimicking the online control of lactose hydrolysis in milk and lactose hydrolysis was monitored first in a model system, then in a milk sample -with no need for deuterated solvents through the use of a tailored solvent suppression pulse sequence (Gouilleux, Charrier, Akoka, & Giraudeau, 2017a). Kinetic analysis of the reaction was also performed. In the second part, 36 milk samples from Turkey and France; having different fat and sugar contents; production styles (conventional vs organic) were analysed by NMR spectroscopy together with multivariate approaches to assess the potential of benchtop NMR as a differentiation tool.

MATERIALS AND METHODS

Materials

Lactose in the form of D-(+)-Lactose monohydrate was purchased from Fluka BioChemika, Switzerland. Lactase enzyme (Maxilact ® LGi 5000) for lactose hydrolysis was supplied from DSM, (Heerlen, The Netherlands). D-(+)-Glucose, D-(+)-Galactose, imidazole and 3-(Trimethylsilyl) propionic-2,2,3,3-d 4 acid (TSP) were provided from Sigma-Aldrich Co. (St.

J o u r n a l P r e -p r o o f

Louis, MO, USA). HPLC grade water from a water purification system (Nanopure Infinity, Barnstead International, IA) was used for the preparation of the solutions. For the lactose hydrolysis experiments in milk, a UHT whole milk (Lactel, France) bought from a grocery store in France was used.

For the second part of the study, 36 different samples of milk (different cow milks: whole milk, organic whole milk, low-fat milk, lactose-free milk and goat milk) were purchased from local Turkish and French supermarkets. Information on the milk samples are given in Table 1.

Sample Preparations

For online monitoring analysis, hydrolysis reactions were performed both for a model system and on a milk sample. As the model solution, a solution of glucose, galactose and lactose was used. 5% (w/w) glucose, galactose and lactose solutions were prepared separately by dissolving in non-deuterated HPLC grade water in a 50 mL tube. TSP was used as the internal standard.

For milk differentiation experiments, 0.5% (w/w) imidazole was added to the 36 milk samples (Table 1) as the internal reference for NMR measurements since TSP was masked with the peaks of the fats present in the milk. Imidazole had 1 H chemical shifts of 7.26 ppm and 8.11 ppm. All the samples were stirred for 5 min for complete dissolution of imidazole in the milk.

Afterwards, the samples were transferred into 5 mm NMR tubes and stored in the dark and refrigerated before analysis to prevent spoilage.

NMR Spectroscopy

NMR spectroscopy experiments were performed on a benchtop spectrometer operating at a 43 MHz 1 H frequency with a compact permanent magnet based on the Halbach design (Magritek Spinsolve, Wellington , NZ) (E. [START_REF] Danieli | Highly stable and finely tuned magnetic fields generated by permanent magnet assemblies[END_REF][START_REF] J O U R N A L P R E -P R O O F Danieli | Small magnets for portable NMR spectrometers[END_REF]. The spectrometer was equipped with a gradient coil along the B 0 axis in the transverse plane of the NMR tube that could produce a maximum field gradient of 0.16 T.m -1 and also had an external lock system which allowed the use of non-J o u r n a l P r e -p r o o f deuterated solvents. The experiments were performed at 29 °C since this was the temperature at which the magnet stability was optimal. Since deuterated solvents were not used, an experiment to suppress the water signal was used, based on a WET-180-NOESY [START_REF] Gouilleux | Gradient-based solvent suppression methods on a benchtop spectrometer[END_REF]. This pulse sequence provides an optimal solvent suppression for small molecules on a benchtop spectrometer, leading to a lower and narrower water signal with a clean phase with a minimal impact on nearby peaks.

An inversion recovery sequence was used for measuring longitudinal relaxation times T 1 using an inversion time range of 0.1-10,000 ms with 15 points.

The 1D 1 H spectra were obtained with 64 scans for a total experiment time of 6 min for model lactose solution hydrolysis and 7 min for milk hydrolysis and 'milk differentiation' experiments. The 90 o pulse angle was achieved by a pulse length of 6.7 µs at 0 dB. The FIDs were recorded with 16 K points, a dwell time of 200 µs, and a repetition time of 6 s for model lactose solution and 7 s for milk, corresponding to 5 times the longest T 1 in the sample to ensure quantitative conditions. 1D data were processed with MestReNova. To align all the spectra correctly, TSP was used to calibrate the chemical shift axis at 0 ppm for the hydrolysis experiments. For milk spectra, the signal of imidazole at 7.26 ppm was used as the chemical shift reference. All spectra were processed with a 0.2 Hz exponential apodization, an automatic phase correction and an automatic baseline correction via a Whittaker smoother algorithm. Manual phase correction and manual baseline correction were also performed in addition to the automatic corrections. Relevant peak integrals were calculated by integration with the MestReNova software. Deconvolution tools were also evaluated but yielded a slightly lower performance than integration, probably due to the non-ideal line shapes.

Lactose hydrolysis in continuous flow system

For on-line monitoring, the system included a glass flow cell with 5 mm outer diameter, a peristaltic pump (Reglo Digital, Ismatec, Wertheim, Germany) and PEEK tubing. A heating plate (RCT Basic, IKA-Werke GmbH & Co. KG, Staufen, Germany) was added to the system J o u r n a l P r e -p r o o f to control the temperature of the hydrolysis reaction.

NMR signal is very sensitive to the flow rates as the signal becomes broader and the resolution decreases with increasing flow rates. The approach followed in a previous study [START_REF] Soyler | Real-time benchtop NMR spectroscopy for the online monitoring of sucrose hydrolysis[END_REF] was used to determine the optimal flow rate of 0.5 mL/min.

Differentiation milk samples

Experiments were performed in 5 mm NMR tubes in static mode. For the quantitative analysis, lactose, glycerol and fat contents were determined from the spectra as will be shown later.

Considering the factors given in Table 1 (country (TR/FR), milk type (Goat/Cow), production style (Conventional/Organic), processing (UHT/Pasteurized), fat and sugar content), Analysis of Variance (ANOVA), was conducted on the responses of glycerol; lactose and fat contents to see if any of the responses showed differences with milk type. If significant difference was detected, means were compared by the Tukey test (p < 0.05) using MINITAB (Version 19, Coventry, U.K) software.

In addition, to obtain some grouping info for the given data set, an Artificial Neural Networks (ANN) model was developed by using NMR spectra results. Artificial Neural Networks (ANN) model is a nonlinear mathematical model with the capability of developing meaningful relationships between input and output variables through a learning process [START_REF] Zheng | Application of Artificial Neural Network (ANN) and Partial Least-Squares Regression (PLSR) to predict the changes of anthocyanins, ascorbic acid, total phenols, flavonoids, and antioxidant activity during storage of red bayberry juice based on fractal ana[END_REF]. The function and organization of the human brain gave inspiration to its function and structure [START_REF] Bila | Accurate wavelet neural-network-based model for electromagnetic optimization J o u r n a l P r e -p r o o f of microwave circuits[END_REF]. The mathematical model designed to classify milk samples with respect to their fat and sugar contents, milk source, production and processing type was based on nonlinear ANNs. These intelligent algorithms have been used to link chemical information coming from 36 milk samples with milk type to achieve a potential tool for accurate differentiation and quality control.

ANN model was implemented with the 'neural net package' (Version 1.44. A value for R 2 between 0.50 and 0.65 can be used for discrimination between high and low concentrations. A value for R 2 between 0.66 and 0.81 shows approximate quantitative predictions., R 2 between 0.82 and 0.90 are considered as good prediction and above 0.91 are excellent [START_REF] Saeys | Potential for onsite and online analysis of J o u r n a l P r e -p r o o f pig manure using visible and near infrared reflectance spectroscopy[END_REF]. Categorical variables were used as binary data, while numeric data was used for normalization and scaling.

RESULTS and DISCUSSION

I. Lactose hydrolysis in continuous flow system

1D solvent suppression enhanced NMR Spectra

At first, 1D 1 H NMR experiments were carried out with a water signal suppression pulse sequence on model solutions of isolated compounds to identify the peaks which could be used to monitor the lactose hydrolysis reactions. Fig. 1 shows the stacked spectra of glucose, galactose and lactose alone. Due to low magnetic field, the spectra of the sugars were highly overlapped. The suppressed water peak was observed at 4.9 ppm. As can be seen in Fig. 1, alpha anomeric proton peaks of glucose, galactose and lactose are overlapped at 5.3 ppm.

Although beta anomeric proton peaks of glucose and galactose are masked by the water peak, the beta anomeric proton peak of lactose is seen as a distinctive peak at 4.4 ppm. Despite the fact that it is difficult to integrate only alpha lactose site due to alpha glucose and galactose signal overlap, it is still possible to monitor the lactose hydrolysis with the isolated beta anomeric signal. Moreover, anomeric ratio between alpha and beta anomers is constant during J o u r n a l P r e -p r o o f the monitoring therefore it is reliable to follow the beta anomeric site along the course of the hydrolysis.

Online Monitoring and Kinetic Modelling of Lactose Hydrolysis in the Model System

Lactose when hydrolyzed yields glucose and galactose as seen below.

For model lactose solution hydrolysis, a 5% lactose solution was hydrolyzed with 20 µL lactase enzyme. The NMR spectra were obtained at 6 min intervals for 240 min (Fig. 2a). The consumption of lactose can be seen from ß-lactose peaks (Fig 2b).

Concentrations of lactose for each time point were calculated using Equation 1. The number of protons contributing to the signal was 9 for TSP, and 3 for imidazole. Lactose data were fitted to Equation 2 (R 2 >0.99) and the average rate constant was found as k=1.66*10 -2 min -1 at 29°C.

= (Eq. 2)
J o u r n a l P r e -p r o o f

Lactose Hydrolysis in Milk

Based on the good repeatability of the model lactose solution hydrolysis, milk was also hydrolyzed with 20 µL of enzyme. Milk is a much more complex fluid than the model lactose solution. The major constituents of milk are water, fat, proteins and lactose. Milk also contains trace amounts of vitamins, minerals, organic acids, enzymes. Therefore, 1D 1 H NMR spectrum of milk was quite different than the model lactose solution (Fig. 4a). The NMR spectra were obtained at each 7 min for 280 min. As seen in Fig. 4b, the consumption of lactose was followed using the peak of ß-lactose. The average rate constant from the fractional conversion model yield a value of 1.52x10 -2 min -1 at 29°C which was consistent with the previous study of Panesar (2007) where it was found as 1.56x10 -2 min -1 at 30 o C (Fig. 5).

Results of this part showed that lactose hydrolysis can be controlled online in a continuous mode successfully by using benchtop NMR spectroscopy, thus the method has high potential to be adapted for industrial processes. For lactose free milk production, it should be ensured that residual lactose amounts are within the legal limits for that product category. This can further be confirmed by other high-resolution spectroscopic techniques.

II. Differentiation of Milk Samples by Benchtop 1 H-NMR Profiling and

Chemometrics

In this part of the study, 36 milk samples of which the properties have been listed in Table 1 were analyzed. A sample spectrum of milk was given in Figure 4a. As seen in the spectrum, there were broad fat peaks between 0.5 -2.5 ppm. For the differentiation experiment of the samples, integral area of glycerol, lactose and fat peak were used. The imidazole peak integral area was used for the normalization of other peaks. The normalized integral areas are given in Table 2. Lactose-free milks (LFM(T/F)) were found to have no lactose and half-fat milks (HFM(T/F)) had nearly half fat content compared to whole milk samples.

J o u r n a l P r e -p r o o f

Glycerol peak was observed next to ß-lactose peak between 4.0 -4.25 ppm. When the glycerol contents were compared, Turkish and French milks were found to be significantly different from each other; glycerol contents of French milks being higher (p<0.05).

Two reasons could suggest for the presence of glycerol in the milk. Firstly, glycerol is used for the treatment of ketosis which is a metabolic disease caused by the negative energy balance during lactation [START_REF] Johnson | The treatment of ketosis with glycerol and propylene glycol[END_REF]. And secondly, glycerol is used in the feed as an energy source. Expansion of the biodiesel industry decreased the availability of corn for animal feed and increased its prices. Therefore, glycerol which is also a by-product of the biodiesel industry has started to be used as an energy source for cows [START_REF] Carvalho | Replacing corn with glycerol in diets for transition dairy cows[END_REF]. The energy concentration of glycerol (1.98-2.29 Mcal/kg) is almost equal to corn starch [START_REF] Schröder | Glycerol as a by-product of biodiesel production in diets for ruminants[END_REF]. And inclusion of that in the diet was eventually observed in the milk.

According to ANOVA results, only the fat and sugar status of the milk were found to be significant for fat and lactose contents (p<0.05). Production type, country and processing were found to be insignificant on lactose and fat content (p>0.05).

Artificial Neural Network (ANN) Model

As stated before, ANN was used as the multivariate approach for classification of the different milk samples based on the outputs of the NMR spectra. The most commonly used quantitative prediction multivariate technique is Partial Least Square Regression (PLSR), which is mostly used to analyse spectroscopic data, especially FTIR. However, ANN was found to show much better differentiative prediction ability as compared to other classification models [START_REF] Efenberger-Szmechtyk | Implementation of chemometrics in quality evaluation of food and beverages[END_REF]. This is explained with several reasons. ANNs learn relationship between the input variables and the output values through successive trainings and the non-linear technique allows model to fit the data better. Analysis is noise insensitive which provides accurate prediction in the presence of noisy data [START_REF] Berrueta | Supervised pattern recognition in food analysis[END_REF]. There has been several studies on foods which show that ANN is superior to J o u r n a l P r e -p r o o f PLSR in prediction [START_REF] Panagou | A comparison of artificial neural networks and partial least squares modelling for the rapid detection of the microbial spoilage of beef fillets based on Fourier transform infrared spectral fingerprints[END_REF][START_REF] Perai | A comparison of artificial neural networks with other statistical approaches for the prediction of true metabolizable energy of meat and bone meal[END_REF][START_REF] Vasquez | Comparison between artificial neural network and partial least squares regression models for hardness modeling during the ripening process of Swiss-type cheese using spectral profiles[END_REF].

In this work, ANN model has been designed and optimized to discriminate milk samples according to their glycerol, lactose and fat peak areas obtained from the 1D NMR spectra. The mathematical tool employed six categorical parameters (country, fat content, sugar content, milk source, production type and processing type) as variables to accurately classify the milk samples.

RMSECV values of the model were used to determine ruggedness of the model. The breakpoint in the RMSECV vs. iteration number (step) plot were identified as the number of iteration (step) parameter for the analyses. Following such an approach prevented overtraining of the model. Since the RMSECV values are usually an integral part of the ANN analysis, they were not saved thus not reported in the model. Feedforward supervised method was used for training ANN model with 5 neurons and 2 neurons in the two hidden layers.

Multilayer feedforward neural network diagrams are also given in the supplement. RMSEC value was obtained as 0.017 and R 2 cal was 0.997. When external validation data was fitted to model, RMSEP was observed as 0.128 and R 2 pred was 0.819 and these values indicated good prediction. Each sample was correctly classified in the cross validation (24/24), while one of the 12 samples was misclassified in the external validation set (Fig. 6). These results showed that glycerol, fat and lactose contents, supplied with information such as production type, processing method, milk source and analysed with ANN, could be used to locate the origin (geography) of milk.

ANN was also used to differentiate organic milk samples from conventional milks samples, and pasteurized milk samples from UHT processed milk samples. For classification of organic and conventional samples, RMSEC value was obtained as 0.121 and R 2 cal was 0.893. When external validation data was fitted to model, RMSEP was observed as 0.232 and R 2 pred was 0.419. Although R 2 pred was relatively low, 0.42 is adequate for discrimination between distinct classes. Prediction ability of model according to samples' production type was 24 out of 24 for J o u r n a l P r e -p r o o f calibration set and 9 out of 12 for validation set. In the validation set, two organic samples were misclassified as conventional and one sample was misclassified as organic.

Same ANN parameters were used to classify pasteurized milk samples and UHT processed milk samples. RMSEC for calibration set was 0.018 and RMSEP for prediction was 0.153. R 2 for calibration was obtained as 0.999 while R 2 for prediction was 0.851. Each sample in calibration set and pasteurized samples in validation set were classified correctly while two of UHT processed milk samples in validation set were misclassified. R 2 for prediction value showed that model was good enough for prediction.

In summary, the ANN model was successfully used to predict origins of milk samples, production type and process type. For all parameters evaluated in the study, the maximum number of incorrect matches was 2 which was a good indicator for the strong classification ability of the model. These results confirmed the existence of a relationship between the six categorical parameters and three main components in the classification of milks, leading to potential design of an objective evaluating method for classification of diverse milk samples.

The model supported potential use of benchtop NMR spectra for differentiating milk samples and provided input for a more comprehensive authentication study.

CONCLUSION

Milk is a very valuable food for the human being. Ensuring its quality and authenticity is significant. In this study, multiple aspects of benchtop NMR spectroscopy have been evaluated

to investigate its potential as a multi-scale control tool at different stages of milk production.

The first part of the study focused on the online monitoring of lactose hydrolysis in milk.

Enzymatic hydrolysis of lactose was successfully followed and a kinetic model was obtained in agreement with previous studies. In the second part, milk samples that have different properties were examined. Glycerol, fat, sugar content of the milks was obtained from the NMR spectra and an Artificial Neural Network (ANN) model was developed to seek for classifications based on the production, processing type, milk source and geographic origin.

The developed model supported the potential use of benchtop NMR spectra for differentiating J o u r n a l P r e -p r o o f milk samples. Overall, it was shown that even for a complex fluid like milk, benchtop NMR spectroscopy provides valuable information and the technique after a full analytical validation that is out of the scope of this proof-of-concept study has the potential to be used in industry for monitoring production, quality control and authenticity determination (but on larger sample set). J o u r n a l P r e -p r o o f J o u r n a l P r e -p r o o f

  2) in R (Stefan Fritsch, Frauke Guenther, Marvin N. Wright, Marc Suling, & Sebastian M. Mueller, 2019). For training (calibration), 2/3 of the data were used and 1/3 was separated to test the model J o u r n a l P r e -p r o o f (external validation). In order to obtain the best score, ANN model was trained using different number of neurons between 2 to 15 in the hidden layers. ANN model yielded possible classifications of the samples and possibility ratios. R 2 cal was used to investigate the prediction ability of the trained ANN model and R 2 pred for external validation. Error values were stated as root mean square error of calibration (RMSEC) and root mean square error of prediction (RMSEP). R 2 values should be close to 1 while error values should be low to obtain a robust prediction model (Uncu & Ozen, 2019). Optimal number of learning steps was selected according to the minimum root means square error cross validation (RMSECV) value for the calibration dataset.

  MW were the concentrations in mg/g, the molecular weight. REF referred to the internal reference. A and N, denoted the integral area in a fully relaxed 1 H NMR spectrum and the number of hydrogens contributing to the signals. Hydrolysis experiments were repeated three times and the change on the average concentrations of lactose from the experiments were plotted as a function of time in Fig. 3.

Funding:

  Alper Soyler got funding from Short Term Scientific Mission (STSM) program under COST Action CA15209 EURELAX -European Network on NMR Relaxometry (http://www.cost.eu/COST_Actions/ca/CA15209) to conduct experiments at the University of Nantes. Authors from Nantes acknowledge support from the Corsaire metabolomics core facility (Biogenouest). J.F. acknowledges Sandrine Bouchet for a continuous essential assistance. Middle East Technical University, Scientific Research Projects Coordination Unit with grant # BAP-07-02-2017-004-077 is also acknowledged.

Figure Captions Figure 1 .Figure 3 .

 Captions13 Figure Captions

Figure 4 .Figure 5 .Figure 6

 456 Figure 4. 1D 1 H NMR spectrum of milk with the assignment for its main constituents (a)

Fig. 1 J

 1 Fig. 1

  Fig. 2

Table 2 .

 2 Integral areas of glycerol, lactose and fat NMR peaks of analysed milks, normalizedto the imidazole signal.

TABLES

Table

1

Milk samples used in the study J o u r n a l P r e -p r o o f

J o u r n a l P r e -p r o o f

J o u r n a l P r e -p r o o f


Jean Furstoss
email: furstoss@geoazur.unice.fr

Carole Petit

Andrea Tommasi

Clément Ganino

Daniel Pino Muñoz

Marc Bernacki

On the role of solute drag in reconciling laboratory and natural constraints on olivine grain growth kinetics

Keywords: Microstructure, numerical modelling, mantle processes

We investigate the effect of solute drag on the grain growth (GG) kinetics in olivine-rich rocks through full field and mean field modelling. Considering a drag force exerted by impurities on grain boundary migration allows reconciling laboratory and natural constraints on olivine GG kinetics. Solute drag is implemented in a full field level-set framework and in a mean field model that explicitly accounts for a grain size distribution. After calibration of the mean field model on full field results, both models are able to both reproduce laboratory GG kinetics and predict grain sizes consistent with observations in peridotite xenoliths from different geological contexts.

Introduction

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 its mechanical behavior through the grain size dependence of the diffusion creep 1 regime [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF]. In this regime, the rock strength is inversely proportional to the grain size to a power p, with p = 2 [START_REF] Karato | Rheology of synthetic olivine aggregates: influence of grain size and water[END_REF].

Moreover grain growth may lead the rocks into the grain-size independent, but stress-dependent dislocation creep regime. A variation in the mean grain size of olivine rocks may therefore produce marked changes in the upper mantle rheology, which may control strain localization [START_REF] Braun | A simple parameterization of strain localization in the ductile regime due to grain size reduction: A case study for olivine[END_REF]. Consistently, preservation of small grain sizes (mylonitic or ultramylonitic microstructures within ductile shear zones [START_REF] Vissers | Shear zones in the upper mantle: A case study in an alpine iherzolite massif[END_REF]) has been proposed as one of the mechanisms allowing to preserve weak plate boundaries over geological times [START_REF] Bercovici | Plate tectonics, damage and inheritance[END_REF].

However, numerical models of grain growth without stored energy [START_REF] Furstoss | 2d and 3d simulation of grain growth in olivine aggregates using a full field model based on the level set method[END_REF][START_REF] Chu | Olivine rheology, shear stress, and grain growth in the lithospheric mantle: geological constraints from the kaapvaal craton[END_REF] based on experimental data on olivine GG [START_REF] Karato | Grain growth kinetics in olivine aggregates[END_REF][START_REF] Ohuchi | Grain growth in the forsterite-diopside system[END_REF][START_REF] Hiraga | Grain growth systematics for forsterite±enstatite aggregates: Effect of lithology on grain size in the upper mantle[END_REF] speak against the persistence of small grain sizes through million years in pure olivine rocks. These models simulate well the experimental data, but they fail even to predict the commonly observed plurimillimetric grain sizes observed in dunites when run on geologically relevant timescales (i.e., over millions of years), for which they predict meter scale grain sizes [START_REF] Furstoss | 2d and 3d simulation of grain growth in olivine aggregates using a full field model based on the level set method[END_REF][START_REF] Chu | Olivine rheology, shear stress, and grain growth in the lithospheric mantle: geological constraints from the kaapvaal craton[END_REF]. To obtain a GG kinetics compatible with natural observations in upper mantle rocks, the presence of second phases is often considered [START_REF] Hiraga | Grain growth systematics for forsterite±enstatite aggregates: Effect of lithology on grain size in the upper mantle[END_REF][START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF]. Nevertheless, to preserve very small olivine grain sizes (<< 100µm) over millions of years, these models have to introduce some questionable features to increase the impediment of olivine grain boundary migration (GBM), such as very slow GG of the second phases (Nakakoji and Hiraga, 2018), grain size reduction of all phases and enhanced mixing between phases [START_REF] Bercovici | Mechanisms for the generation of plate tectonics by two-phase grain-damage and pinning[END_REF] or the presence of small fixed particles [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF] .

Our understanding of olivine GG in itself is limited by the inconsistency between laboratory and natural time and spatial scales. All GG models are calibrated on experimental data obtained for ultra-fine grains (1 to 30µm), for which the grain size evolution is measurable at experimental timescales (a few hours to several days), and then extrapolated to natural conditions. As a consequence, if a physical mechanism becomes prominent only for grain sizes larger than 50µm, its effect will not be captured in the experiments. An example of such a physical mechanism is the drag force exerted by impurities also called solute drag, which has been very seldom considered in studying GG in rocks (orthopyroxene [START_REF] Skemer | Effects of solute segregation on the graingrowth kinetics of orthopyroxene with implications for the deformation of the upper mantle[END_REF], halite [START_REF] Guillope | Dynamic recrystallization during creep of single-crystalline halite: An experimental study[END_REF] and never accounted for modeling GG of olivine.

In the present work we tested the effect of solute drag on the GG kinetics of olivine. We performed 2D full field GG simulations using a level-set (LS) framework [START_REF] Bernacki | Level set framework for the finite-element modelling of recrystallization and grain growth in polycrystalline materials[END_REF] in which we implemented the solute drag effect. These models show that considering this mechanism allows for consistent simulation of olivine GG kinetics at both experimental and upper mantle conditions. We then adjust a mean-field model on results of the full field simulations to propose an analytical expression allowing to compute efficiently the grain size evolution within geodynamic large scale models accounting for grain size-dependent rheologies.

Modeling solute drag

Impurities present within the crystal lattice and segregated at grain boundaries can have an impact on their migration kinetics through the so called solute drag effect. The grain boundary migration velocity (v) is generally expressed as [START_REF] Humphreys | Recrystallization and related annealing phenomena[END_REF] :

v = M P, (1)
where M (m 4 .J -1 .s -1 ) is the grain boundary mobility and P is the sum of the pressures exerted on the grain boundary. The solute drag effect can be described in terms of drag pressure exerted on the grain boundary, which is a function of the grain boundary velocity and impurities concentration and nature. This drag pressure has been quantified theoretically [START_REF] Lücke | A quantitative theory of grain-boundary motion and recrystallization in metals in the presence of impurities[END_REF] . Its intensity follows three main regimes (high, intermediate and low velocity), similarly to dislocation impediment by Cottrell atmospheres [START_REF] Cottrell | Dislocation theory of yielding and strain ageing of iron[END_REF]. Impurities segregated around the grain boundary interact with it and when the grain boundary migrates, the impurity cloud tends to accompany it. In the high velocity regime, the grain boundary moves so fast that the segregated impurities cannot follow the interface and the interactions between the impurities and the grain boundary are strongly reduced. In the low velocity regime, the impurity cloud stays segregated around the moving grain boundary, but the intrinsic drag of the grain boundary due to the intrinsic defects within the interface is generally higher than the drag exerted by the impurities (Fig. 1). Between these two regimes, an intermediate regime exists where the impact of solute drag on GBM is important.

Even if the drag effect is expected to follow different mathematical relationships depending on the velocity regime [START_REF] Lücke | A quantitative theory of grain-boundary motion and recrystallization in metals in the presence of impurities[END_REF], an unified expression for the drag pressure P i exerted by a c 0 concentration of impurities (within the bulk) valid for the three regimes has been proposed by [START_REF] Cahn | The impurity-drag effect in grain boundary motion[END_REF] :

P i = αvc 0 1 + β 2 v 2 , ( 2
)
where v is the grain boundary velocity norm, c 0 the concentration of im-purities within the bulk, and α(J.s.m -4 ) and β(s.m -1 ) are two parameters modulating the intensity of the solute drag. α controls the intensity of the solute drag pressure and β constrains the grain boundary velocity which is the most impacted by the presence of the impurities (Fig. 1). The solute drag effect is most effective for grain boundary velocities close to 1/β. The α and β parameters can be expressed as integrals of functions D(x) and E(x) representing the variation of the impurity diffusion coefficient and of the interaction energy, respectively, along the grain boundary normal (x). As these two functions are difficult to constrain experimentally, their mathematical expressions are generally hypothesized in an integrable way, which permits to obtain analytical expressions for α and β [START_REF] Cahn | The impurity-drag effect in grain boundary motion[END_REF] :

α = N v (k b T ) 2 E 0 D (sinh( E 0 k b T ) - E 0 k b T ), (3)
and,

β 2 = αk b T δ 2N v E 2 0 D , (4)
where N v (m -3 ) is the number of atoms per unit volum, k b the Boltzmann constant, E 0 (J) the interaction energy, D(m 2 .s -1 ) the diffusion coefficient of the impurity within the bulk and δ(m) is the characteristic segregation length of the impurities around the grain boundary.

Semi-explicit implementation of solute drag within the level-set framework

The LS framework [START_REF] Bernacki | Finite element model of primary recrystallization in polycrystalline aggregates using a level set framework[END_REF], already used to model olivine GG [START_REF] Furstoss | 2d and 3d simulation of grain growth in olivine aggregates using a full field model based on the level set method[END_REF][START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF], proposes an implicit description of the polycrystal through the use of LS functions representing the signed distance function to the grain boundaries surrounding the grain they represent (positive inside the grain and negative elsewhere) in a finite element (FE) context. The microstructural evolution is simulated by moving the LS functions according to physical laws describing GBM [START_REF] Cruz-Fabiano | Assessment of simplified 2d grain growth models from numerical experiments based on a level set framework[END_REF] or by creating LS functions to represent new grains nucleated during the recrystallization mechanisms [START_REF] Maire | Modeling of dynamic and post-dynamic recrystallization by coupling a full field approach to phenomenological laws[END_REF]. Theoretically, each grain of a polycrystal is represented by its own LS function. However to reduce the computation time and memory storage, several non-neighboring grains in the initial microstructure can be grouped to form Global Level Set (GLS) functions thanks to a graph coloration technique. Re-coloration is then used to avoid numerical grain coalescence during grain boundary motion [START_REF] Scholtes | New finite element developments for the full field modeling of microstructural evolutions using the level-set method[END_REF][START_REF] Scholtes | 3d level set modeling of static recrystallization considering stored energy fields[END_REF]. The GLS functions dis-placement is computed within an efficient FE framework [START_REF] Scholtes | New finite element developments for the full field modeling of microstructural evolutions using the level-set method[END_REF][START_REF] Shakoor | An efficient and parallel level set reinitialization method-application to micromechanics and microstructural evolutions[END_REF] using anisotropic mesh refinement around interfaces [START_REF] Resk | Adaptive mesh refinement and automatic remeshing in crystal plasticity finite element simulations[END_REF]. The initial microstructure is generated using a Voronoï-Laguerre Dense Sphere Packing algorithm [START_REF] Hitti | Precise generation of complex statistical representative volume elements (rves) in a finite element context[END_REF], which respects precisely an imposed initial grain size distribution.

If GG is solely controlled by capillarity (reduction in the grain boundary energy)

and solute drag, the grain boundary velocity is controlled by the difference between the capillarity pressure (P capi ) and solute drag pressure (P i described by Eq.2) as :

v = M (P capi -P i ) n = M (-γκ - αc 0 v 1 + β 2 v 2 ) n, (5)
where M (m 4 .J -1 .s -1 ) and γ(J.m -2 ) are the grain boundary mobility and energy respectively, κ is the grain boundary local curvature (in 2D) or the sum of the main local curvatures (in 3D) and n the outward unit normal to the grain boundary.

By introducing the GLS function Φ i to represent the grains i, we can obtain the following implicit and explicit first order time discretizations for the temporal derivative of Φ i :

∂Φ i ∂t = Φ t+∆t -Φ t ∆t , (6)
for the GLS function velocity :

v = Φ t+∆t -Φ t ∆t , (7)
and for the squared velocity :

v 2 = ( Φ t -Φ t-∆t ∆t ) 2 = v 2 old . (8)
The displacement of the GLS functions is computed using the convective LS equation [START_REF] Osher | Fronts propagating with curvaturedependent speed: algorithms based on hamilton-jacobi formulations[END_REF] :

∂Φ i ∂t + v • ∇Φ i = 0. (9)
Considering the geometrical properties of LS function κ = -∆Φ, ∇Φ. ∇Φ = 1, ∇Φ = -n, using the above time discretization and substituting Eq.5 within the above 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 0 1 + β 2 v 2 old . (11)
The fact that Eq.11 uses v 2 old instead of v 2 is a simplification that reduces the non-linear behaviour of the problem. In practical terms, the time marching scheme uses a small timestep, thus the error of replacing v 2 by v 2 old is small and this greatly simplifies the numerical scheme allowing its implementation in a generic FE code. This formulation is similar to the diffusion formulation classically used for capillarity driven GG in the LS formalism [START_REF] Cruz-Fabiano | Assessment of simplified 2d grain growth models from numerical experiments based on a level set framework[END_REF]. By analogy with the heat diffusion equation, M is equivalent to a mass term (i.e. the product between the specific heat capacity and density).

The main difference relative to the classic LS approach for capillarity-driven GG is that when solute drag is modelled the mass term differs from one and depends on the velocity of the LS function at time t -∆t. This heterogeneous mass term is computed on each node of the mesh and linearly interpolated.

A major drawback of the LS approach lies in the fact that during grain boundary migration, the GLS are no longer distance functions || ∇φ|| = 1. This is par-ticularly problematic when a remeshing technique depending on the distance property is used at grain interfaces. In addition, the new diffusive formulation proposed in Eq.10 requires a distance function as it is based on the respect of || ∇φ|| = 1, at least in a thin layer around the interface. For these reasons, the GLS functions need to be reinitialized at each time step in order to restore their metric property. Numerous approaches exist for this reinitialization procedure.

Here we use a fast and accurate approach usable in unstructured FE mesh proposed by [START_REF] Shakoor | An efficient and parallel level set reinitialization method-application to micromechanics and microstructural evolutions[END_REF]. The residual errors inherent to this approach are discussed in [START_REF] Florez | A new finite element strategy to simulate microstructural evolutions[END_REF].

The introduction of the solute drag pressure is expected to reduce the GG kinetics. Thus, its influence has to be accounted for in the adaptative time stepping scheme to allow larger steps when the grain size evolves slowly. The timestep is computed by imposing a maximal incremental displacement corresponding to a

given fraction (H) of the LS reinitialized width (E p ) :

∆t = HE p v m , (12)
where v m is the grain boundary mean velocity :

v m = M ( γ R - c 0 α| R| 1 + β 2 R2 ), (13)
where R and R are the mean grain radius and its temporal evolution, respectively. To avoid negative timestep values, we impose a lower bound for v m at M γ/R max , where R max is the maximum grain radius.

The explicit time discretization in Eq.8 may have an impact on the numerical resolution 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 case with solute drag. Thus, the strong formulation (Eq.10) for the displacement of LS functions accounting for solute drag can be validated. When accounting for solute drag by impurities, the GG kinetics follows the expected trend and is much slower than in cases without drag (Fig. 2).

Figure 2: Full field model of the mean grain radius evolution using M = 2.9 • 10 -2 mm 4 .J -1 .s -1 , γ = 10 -6 J.mm -2 (representing olivine material parameters at 1573K from [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF]) for a case without impurities (dashed lines) and for a case with c 0 = 1000ppm, α = 10 5 J.s.mm -4 and β = 10 5 s.mm -1 (solid lines). The H coefficient controls the timestep (Eq.12).

Mean field approach for GG with solute drag

To construct a mean field model describing GG kinetics including solute drag, it is important to account for the grain size distribution. Computing GG kinetics using only the mean grain size evolution will hide the dispersion of the solute drag pressures exerted within the microstructure due to variations in the capillarity force, which is a function of grain size. This will result in three clear stages of grain size evolution as illustrated in Fig. 1 : an initial mean grain size evolution unimpacted by solute drag (high velocity regime), a second phase, the most impacted by solute drag, where the mean grain size will be quasi-static and a final phase also unimpacted by solute drag (low velocity regime). To account for the coexistence of these three regimes within the microstructure and predict realistic mean grain size evolutions, the initial dispersion of individual grain sizes and their evolution rates have to be considered.

For this purpose, we adapted the [START_REF] Hillert | On the theory of normal and abnormal grain growth[END_REF] model, which proposes a discrete representation of the grain size distribution (GSD) in the microstructure. By considering R i the radius of the i-grain size bin, the Hillert's model allows computing the evolution of each bin by accounting for the capillarity pressure [START_REF] Maire | Improvement of 3d mean field models for capillarity-driven grain growth based on full field simulations[END_REF]. This procedure enables to follow the evolution of each bin of the initial GSD. This mean field approach is able to reproduce, in terms of GSD, the predictions of full field simulations in context of pure GG in 2D [START_REF] Cruz-Fabiano | Assessment of simplified 2d grain growth models from numerical experiments based on a level set framework[END_REF] and in 3D [START_REF] Maire | Improvement of 3d mean field models for capillarity-driven grain growth based on full field simulations[END_REF] even for complex initial GSD (like bimodal ones).

expressed as M γκ R 2 R 2 ( R R 2 -1 Ri ), with κ ≈ 1.6 from
To account for solute drag, we subtract from the capillarity pressure the solute drag pressure (Eq.2) replacing v by Ṙi and v 2 by Ṙ2 i,old (i.e. the grain size evolution rate at the precedent increment), which gives :

Ṙi = κ R 2 R 2 M (γ( R R 2 - 1 R i ) - αc 0 Ṙi 1 + β 2 Ṙ2 i,old ), (14)
which can be reformulated as :

Ṙi = γ( R R 2 - 1 R i )/( R 2 M κR 2 + αc 0 1 + β 2 Ṙ2 i,old ), (15)
To account for topological effects or a non-uniform mass term along grain bound-aries, equation 15 is generalized through the following expression :

Ṙi = γ( R R 2 - 1 R i )/( R 2 M κR 2 + C α αc 0 Ṙn i,old 1 + (C β β) 2 Ṙ2 i,old ), (16)
where C α , C β and n are mean field parameters which have to be calibrated on full field simulations. It can be noticed that the equation 16 is equivalent to the non-generalized form (Eq.15) for n = 0 and

C α = C β = 1.
In practice, each bin R i are equipropables in number. We begin with a list of bins generated from an imposed initial GSD and compute iteratively their evolutions. If the grain radius of a bin becomes lower than 0.1µm it is considered as consumed by the growth of neighboring grains and removed from the bin list.

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 number of bins after long annealing time is very large and the computational cost becomes prohibitive. Thus, we use a repopulation strategy allowing to do the calculation 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 bin list by adding ten new bins for each existing bin R i with radii ranging between 0.9R i and 1.1R i with a step of 0.02R i . Those range of values have been chosen in order to minimize the difference between GSD before and after repopulation step (see an example in Fig. 3). The minimal number of bins is fixed at 40 based on a convergence study.

Solute drag in olivine

In mantle rocks, several incompatible elements and impurities are present in concentrations ranging from few to hundreds ppm [START_REF] De Hoog | Trace-element geochemistry of mantle olivine and application to mantle petrogenesis and geothermobarometry[END_REF]. Some of them are more or less homogeneously dispersed through the bulk material, but others can be enriched at grain boundaries [START_REF] Suzuki | Grain-boundary enrichment of incompatible elements in some mantle peridotites[END_REF]. The elements that can influence grain growth through solute drag are those that exhibit a partitioning between grain interiors and boundaries. The major constituent of mantle rocks is olivine, an orthorhombic material of composition (F e, M g) 2 SiO 4 . In olivine-rich rocks, chromium (Cr), aluminum (Al) and calcium (Ca) [START_REF] Suzuki | Grain-boundary enrichment of incompatible elements in some mantle peridotites[END_REF][START_REF] Hiraga | Chemistry of grain boundaries in mantle rocks[END_REF][START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF] display a strong partitioning between olivine grain matrix and boundaries, the latter being often qualified as «incompatible element reservoirs ». The major parameters controlling their effect on GBM are the element concentration, diffusion coefficient and interaction energy with grain boundary. The latter will control at first order the intensity of the drag pressure exerted by the impurity. As interstitials, vacancies and even grain boundaries could have a non-null electrical charge in materials as olivine, the full formulation of this energy term should include, in addition to an elastic part, an electrostatic component. However for sake of simplicity we will assume a pure elastic interaction corresponding to the lattice distortion due to the misfit between the impurity size and the size of the typical host ion the impurity replaces. Close to and within olivine grain boundaries the impurities may replace Mg ions (a depletion in Mg is observed at grain boundaries) [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF]. Within the grain interiors, Mg ions are located in the octahedral M sites, and the interaction energy between impurity and grain boundary can be computed, as a first order approximation, using the characteristic length r 0 of those sites [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF]) :

E 0 = 4πE α ( r 0 2 (r i -r 0 ) 2 + 1 3 (r i -r 0 ) 3 ), (17)
with E α is the Young modulus of the M lattice site and r i is ionic radius of the impurity. E α is close to the bulk Young's modulus and r 0 can be computed from the length of the bonds between M and oxygen sites, and between oxygen and oxygen sites [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF]. The temperature dependency of these parameters is small and we will consider in the following E α = 159GP a and r 0 = 0.064nm [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF].

Using the above expression and Eqs.2, 3 and 4 one can compute the drag pres- 1.9 • 10 -15 , 1.9 • 10 -16 from [START_REF] Jollands | Substitution and diffusion of cr2+ and cr3+ in synthetic forsterite and natural olivine at 1200-1500 c and 1 bar[END_REF][START_REF] Zhukova | A subsidiary fast-diffusing substitution mechanism of al in forsterite investigated using diffusion experiments under controlled thermodynamic conditions[END_REF][START_REF] Spandler | Survival times of anomalous melt inclusions from element diffusion in olivine and chromite[END_REF])), we estimate that the characteristic drag pressure exerted by Ca is at least 1 order of magnitude higher than the ones exerted by Ni and Al. Thus in the following, the only impurity considered for solute drag will be Ca.

The Ca solute drag parameters α Ca and β Ca can be computed using equations 3 and 4 with a segregation length (δ) of 5nm [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF]) and an

Arrhenius-type temperature dependence of the diffusion coefficient (D) where the reference value and activation energy are D 0 = 3.16.10 -6 mm 2 .s -1 and [START_REF] Coogan | Experimental determination of the diffusion coefficient for calcium in olivine between 900 c and 1500 c[END_REF] respectively. We obtain α Ca = 3.10 7 J.S.mm -4 and β Ca = 9.6.10 5 s.mm -1 at 1573K and α Ca = 3.10 11 J.s.mm -4

Q D = 200kJ.mol -1 (
and β Ca = 2.8.10 9 s.mm -1 at 1073K.

4 Laboratory and natural constraints on grain sizes in olivine-rich rocks

The main goal of this study is to show that laboratory experiments and natural observation on olivine GG can be reconciled by accounting for Ca solute drag.

In the following section, we define a priori constraints on GG kinetics provided by natural and experimental data.

Laboratory constraints

The annealing experiments on ultrafine grained natural San Carlos olivine of [START_REF] Karato | Grain growth kinetics in olivine aggregates[END_REF] seem appropriate to define the experimental reference, particularly the high pressure (300M P a) dry runs at 1573K. This sample has probably as much impurities as a natural mantle rock because it was synthesized using crushed grains of natural olivine, so the solute drag was probably as active in these experiments as in nature. These experiments allow constraining grain growth kinetics for short times (≤ 10h) and small grain sizes (2 -25µm).

Natural constraints

Natural constraints on GG kinetics have to be considered cautiously because of the large number of uncertainties in the determination of physical conditions of natural systems (initial size distribution, precise thermal history, pure static conditions, etc.). Contraints on thermal history, strain and grain size evolution of peridotite xenoliths can be found in the literature (e.g. [START_REF] Pokhilenko | New data on the mineralogy of megacrystalline pyrope peridotite from the udachnaya kimberlite pipe, siberian craton, yakutian diamondiferous province[END_REF][START_REF] Baptiste | Petrophysical constraints on the seismic properties of the kaapvaal craton mantle root[END_REF]), but these rocks are a mixture of olivine, pyroxene and other secondary phases and their microstructural evolution is certainly affected by pinning effects due to their polymineralic nature [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF] in addition to the solute drag effect. Nevertheless, when solely controlled by pinning, the maximum grain size is dictated by the spacing between static pinning phases like spinels [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF] and the latter is always larger than the actual GS in natural polymineralic peridotites, like harzburgites or lherzolites [START_REF] Herwegh | The role of second phases for controlling microstructural evolution in polymineralic rocks: A review[END_REF]. An alternative could be to focus on dunites, which are composed by ≥ 90% of olivine. However, in such rocks, the initial crystal size distribution is clearly related to the mechanism of their formation, which usually involves extensive interaction with melts percolating the mantle [START_REF] Berger | Les dunites en enclaves dans les basaltes alcalins des îles océaniques: approche pétrologique[END_REF][START_REF] Kelemen | Reaction between ultramafic rock and fractionating basaltic magma i. phase relations, the origin of calc-alkaline magma series, and the formation of discordant dunite[END_REF] and may be very different from the GS distribution used as initial conditions for our models. Considering these two limitations, we chose to compare our models results with classical (i.e., polymineralic) peridotites, being aware that our model only captures a part of the processes that limit the maximum grain size of natural rocks.

To focus on natural peridotites with textures typical of thermal (without stored energy) annealing according to the terminology defined by [START_REF] Harte | Rock nomenclature with particular relation to deformation and recrystallisation textures in olivine-bearing xenoliths[END_REF] we selected examples which were described as coarse granular or protogranular.

Porphyroclastic, granoblastic and tabular textures were not considered because they reflect significant rock deformation that was not subsequently fully annealed, or grain growth in presence of melts or fluids. We briefly recall here the geological setting, age and temperature history of the selected examples.

The Udachnaya (Siberia) kimberlite xenoliths

The Udachnaya kimberlite pipe in Siberia is well known due to its diamond mine and because of the occurrence of megacrystalline harzburgite and dunite xenoliths, mainly composed by coarse olivine crystals (up to 10 cm wide; [START_REF] Pokhilenko | New data on the mineralogy of megacrystalline pyrope peridotite from the udachnaya kimberlite pipe, siberian craton, yakutian diamondiferous province[END_REF]). The age of Udachnaya pipe has been determined between 345 and 385 Ma depending on the dating method [START_REF] Dehvis | New data on the age of yakutian kimberlites obtained by uranium-lead study of zircons[END_REF][START_REF] Ilupin | Kimberlites[END_REF]. Equilibrium temperatures and depth of the megacrystalline peridotite xenoliths range between 1173 and 1373K and 150-200 km [START_REF] Pokhilenko | Megacrystalline pyrope peridotites in the lithosphere of the siberian platform: mineralogy, geochemical peculiarities and the problem of their origin[END_REF][START_REF] Griffin | Thermal state and composition of the lithospheric mantle beneath the daldyn kimberlite field, yakutia[END_REF]. Olivine in the megacrystalline harzburgites from Udachnaya have low Ca contents (mainly < 150ppm; [START_REF] Pokhilenko | New data on the mineralogy of megacrystalline pyrope peridotite from the udachnaya kimberlite pipe, siberian craton, yakutian diamondiferous province[END_REF]). It is impossible to determine precisely how much time they have spent at these temperatures before being extracted by the kimberlitic eruption in Late Devonian -Early Carboniferous times. Re/Os ages in diamond sulfides have provided an age for the formation of the cratonic lithosphere around 1.8 Ga [START_REF] Ionov | The age and history of the lithospheric mantle of the siberian craton: Re-os and pge study of peridotite xenoliths from the obnazhennaya kimberlite[END_REF], if we assume that these rocks represent samples from the oldest part of the cratonic lithosphere, we can estimate a maximum residence time of these rocks at 1173-1373K of 1450 Ma. A minimum residence time is difficult to estimate, since the cratonic mantle has been modified by metasomatism (interaction with percolating melts) after its stabilization [START_REF] Ionov | The age and history of the lithospheric mantle of the siberian craton: Re-os and pge study of peridotite xenoliths from the obnazhennaya kimberlite[END_REF].

The Kaapvaal (South Africa) kimberlite xenoliths

The Archean lithosphere of the Kaapvaal craton in South Africa stabilized around 3 Ga ago according to Re-Os isotopic data on peridotite xenoliths [START_REF] Pearson | Re/os, sm/nd, and rb/sr isotope evidence for thick archaean lithospheric mantle beneath the siberian craton modified by multistage metasomatism[END_REF] extracted by kimberlite magmatism in Late Jurassic to Cretaceous times, i.e. between 180 and 90 Ma [START_REF] Griffin | Emplacement ages and sources of kimberlites and related rocks in southern africa: U-pb ages and sr-nd isotopes of groundmass perovskite[END_REF]. Except of the deepest sections of the lithosphere (> 180km) where mylonitic microstructures and lherzolitic compositions are common, most peridotite xenoliths consist of coarse (5-8 mm) to cm-size grained harzburgites in which olivine has low intragranular deformation and rounded to polygonal shapes indicating significant annealing posterior to an early stage of deformation (e.g. [START_REF] Boullier | Classification of textures and fabrics of peridotite xenoliths from south african kimberlites[END_REF][START_REF] Boyd | Composition and structure of the kaapvaal lithosphere, southern africa[END_REF][START_REF] Baptiste | Petrophysical constraints on the seismic properties of the kaapvaal craton mantle root[END_REF]). Ca content in the rims of olivine grains reaches up to 1200 ppm [START_REF] Hervig | Lherzolite xenoliths in kimberlites and basalts: petrogenetic and crystallochemical significance of some minor and trace elements in olivine, pyroxenes, garnet and spinel[END_REF]) and grain interiors have in general lower Ca content (30-600 ppm, median 100 ppm, cf. Georoc database). PT estimates range along a low geothermal gradient with temperatures of 873 -1273K at depths between 80 and 150 km [START_REF] Boyd | Evidence for a 150-200-km thick archaean lithosphere from diamond inclusion thermobarometry[END_REF][START_REF] Chu | Olivine rheology, shear stress, and grain growth in the lithospheric mantle: geological constraints from the kaapvaal craton[END_REF][START_REF] Baptiste | Petrophysical constraints on the seismic properties of the kaapvaal craton mantle root[END_REF]. In the Jagersfontein pipe, very coarse-grained peridotite xenoliths exhibit grain sizes from 5 to 20 mm with temperature estimates ranging from 973 -1223K [START_REF] Winterburn | Peridotite xenoliths from the jagersfontein kimberlite pipe: I. primary and primary-metasomatic mineralogy[END_REF]. As precise constraints on the thermal evolution of the cratonic root are not available, in our model, we assume that the equilibrium cratonic conductive geotherm was reached rapidly after the stabilization of the cratonic lithosphere and that these rocks have spent 2 to 3 Ga at their equilibrium temperatures before being erupted. By using a mean temperature of 1173K, we obtain estimates for the olivine grain sizes at depths of 100-125km [START_REF] Baptiste | Deformation, hydration, and anisotropy of the lithospheric mantle in an active rift: constraints from mantle xenoliths from the north tanzanian divergence of the east african rift[END_REF].

The Kerguelen hotspot in the Indian Ocean

The Kerguelen archipelago is part of a Large Igneous Province, the Kerguelen Plateau, formed above the Kerguelen plume [START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF][START_REF] Mattielli | Kerguelen basic and ultrabasic xenoliths: evidence for long-lived kerguelen hotspot activity[END_REF][START_REF] Grégoire | The kerguelen archipelago: an hypothetic continental mafic protolith[END_REF]. Plume-related volcanism forming the Ker-guelen Islands started around 45 Ma ago and lasted until 0.1 Ma ago ( [START_REF] Cottin | The kerguelen volcanic plateau: the second largest oceanic igneous province (lip) on earth and a witness of the indian ocean opening[END_REF] and references therein). Ultramafic xenoliths brought at the surface in the Kerguelen Islands by the plume-related volcanism are harzburgites and dunites typical of a depleted mantle which has undergone a large degree of partial melting ( [START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF] and references therein). Equilibrium PT conditions determined on xenoliths close to the crust-mantle boundary are around 1 GPa and 1173-1273K [START_REF] Grégoire | The kerguelen archipelago: an hypothetic continental mafic protolith[END_REF]. Some protogranular harzburgites have mean grain sizes of 2-10 mm while equigranular dunites have a mean grain size between 0.5 and 1 mm [START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF]. Here again, it is difficult to estimate the annealing time of these xenoliths. Based on geochemical and petrological analyses the Kerguelen harzburgites were interpreted as residues from a partial melting episode linked with the Kerguelen plume, which were subsequently affected by melt percolation forming the dunites [START_REF] Mattielli | Kerguelen basic and ultrabasic xenoliths: evidence for long-lived kerguelen hotspot activity[END_REF][START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF]. Peridotite xenoliths have been sampled in lavas dated between 28 and 7 Ma ( [START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF] and references therein) so they might have spent up to 38 Ma at temperatures close to 1173-1273K.

The temperature, rock type, mean grain size and inferred residence times for the three contexts are summarized within table 1.

Table 1: Summary of constraints on temperature, grain size and geodynamic context for the four selected peridotite samples that were used to estimate our model performance. a (Baptiste and Tommasi, 2014); b [START_REF] Winterburn | Peridotite xenoliths from the jagersfontein kimberlite pipe: I. primary and primary-metasomatic mineralogy[END_REF]; c [START_REF] Pokhilenko | New data on the mineralogy of megacrystalline pyrope peridotite from the udachnaya kimberlite pipe, siberian craton, yakutian diamondiferous province[END_REF]; d [START_REF] Goncharov | Thermal state, oxygen fugacity and c o h fluid speciation in cratonic lithospheric mantle: New data on peridotite xenoliths from the udachnaya kimberlite, siberia[END_REF]; e [START_REF] Yudin | The age of udachnaya-east kimberlite: U/pb and 40 ar/39 ar data[END_REF]; f [START_REF] Grégoire | The kerguelen archipelago: an hypothetic continental mafic protolith[END_REF]; g [START_REF] Bascou | An integrated study of microstructural, geochemical, and seismic properties of the lithospheric mantle above the kerguelen plume (indian ocean)[END_REF].

The geological contexts presented above will be used in the following to test the performance of our mean-field solute drag model. In a first step, to calibrate the solute drag parameters, we run full-field model searching for solute drag parameters that both reproduce the experimental data and predict grain sizes ranging between 0.5 and 10mm for annealing times ranging between 0.1Ma and 1Ga at a constant temperature of 1073K.

Results

In this section, we present the full field and mean field results for long term annealing of olivine aggregates. The olivine grain boundary energy (γ), calcium concentration within the bulk (c 0 ) are taken as 1J.m -2 [START_REF] Cooper | Rheology and structure of olivine-basalt partial melts[END_REF]) and 100ppm respectively while the grain boundary mobility (M ) is taken as an Arrhenius's law where the reference value and activation energy are [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF] respectively.

equal to M 0 = 4.10 4 mm 4 .J -1 .s -1 , Q M = 185kJ.mol -1 (
To perform long term annealing full field simulations, passing from micrometer to millimeter scale grain sizes, we need to define a model chaining strategy. The full field simulations begin with approximately 2000 grains respecting an initial grain size distribution corresponding to the one used in laboratory experiments [START_REF] Karato | Grain growth kinetics in olivine aggregates[END_REF]. 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 is generated for the next simulation, which takes as input the final grain size distribution of the latest one. The error due to this chained calculation is minimized by sampling very precisely the final GSD and by imposing it using the Voronoï-Laguerre Dense Sphere Packing algorithm [START_REF] Hitti | Precise generation of complex statistical representative volume elements (rves) in a finite element context[END_REF].

In the following, we first propose an adjustment of the solute drag material parameters, based on full field simulations, which permits to reconcile laboratory and natural observations by producing the adequate GG kinetics. Then, the mean field model presented within section 2.2 is calibrated on the results of the heterogeneously distributed between one and its maximum value depending on the local grain boundary velocities, which are controlled by the local curvature.

Within this regime, the long straight grain boundary segments have higher mass terms and 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 pinning force on other grain boundaries, which enhances the deceleration of the grain growth kinetics. Finally, when grain sizes are sufficiently large, the grain boundary velocities become very small and one enters the low velocity regime in which the grain growth kinetics is again weakly impacted by the solute drag (Fig. 1).

The mass term reaches almost its maximum value in the whole microstructure (Fig. 5). However, locally, smaller grains, which may correspond to those being consumed in the annealing process, have lower mass terms and intermediate velocities, closer to the one most impacted by solute drag. Migration of these boundaries is slowed down by solute drag and this affects the whole system.

Thus, in the low velocity regime, solute drag has a weaker, but non-negligible effect on the average grain growth velocity. The analysis of the GG kinetics for simulations with different initial mean grain sizes highlights a first phase of very slow or even null mean grain size evolution (Fig. 6). For the equivalent full field models without accounting for the solute drag, this first phase of very slow GG also exists but is much shorter than for the models accounting for solute drag (Fig. 6). The full field simulations, which explicitly consider solute drag predict that the initial grain size does not affect GG kinetics apart from delaying the start of the rapid GG stage. The duration of this phase of slow GG is longer when the initial mean grain size is coarse (Fig. 6). This initial slow GG kinetics can be explained by analyzing the GSD (Fig. 7). In fact, when the initial microstructure is composed of coarse grains, some grains have to shrink to let the other ones grow. This can be seen by comparing in Fig. 7a the GSD of the initial microstructure and that at the end of the initial slow GG phase for the experiment with the coarser initial mean grain size in Fig. 6. This comparison highlights an enrichment in small grains (Fig. 7a). Before disappearing, the GBM velocities of the shrinking grains will necessarily pass through a velocity regime highly impacted by the presence of impurities. This will affect the grain growth kinetics directly (by slowing down the shrinking) and indirectly by impeding the movements of the other grain boundaries through pinning mechanisms.

After calibration, the mean field model (Eq.16) reproduces correctly the full field modeled GG kinetics for both temperatures (Figs. 4 and6). The best fitting mean field parameters C α , C β and n are equal to 0.7, 0.25 and 0.15 respectively. The GSD predicted by the mean field approach is consistent with the full field modeled ones (Fig. 7) but some differences appear for long annealing times (Fig. 7b). However these differences between full field and mean field predicted GSD after long annealing times seem to do not strongly impact the mean grain size predictions since the mean field simulations are in good agreement with the full field calculations, in terms of mean grain size, even for long annealing times (Figs. 4 and6). This mean field approach allows to be predictive on olivine aggregates GG kinetics for a much lower computational cost, which permits to test our formalism on different geological contexts.

Implication for the microstructural evolutions in ultramafic rocks

To estimate how well our mean field model predicts the average grain size of natural samples with reasonable c 0 (Ca concentration within bulk) values, we test it against different contexts in terms of annealing temperature and measured grain sizes. Figure 7: GSD (in number) for the fourth run presented in figure 6 (the green one, at 1573K for α = 10α Ca and β = 500β Ca ), 7a : initial GSD (blue), mean field (orange) and full field (yellow) modeled GSD after 10 7 s, 7b : mean field (orange) and full field (yellow) modeled GSD after 3.6.10 8 s.

Mean field GG models

For the mean field GG model, the temperature is kept constant and grain size grows 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 1373K

and c 0 values of 600, 800, 1000 and 1200 ppm.

All models were run initially with an initial grain size of 20µm and a standard deviation of 2µm; then, isothermal models were run again with an initial grain size of 0.5mm and a standard deviation of 50µm and an intial grain size of 2mm and a standard deviation of 200µm. In fact, the initial grain size has no effect of the grain growth curve or final grain size, apart from shifting the beginning of the positive slope on the grain growth curve (Fig. 8).

Results

The c 0 parameter, in the range of tested values, has only a moderate effect of the final or intermediate grain sizes (table 2 and figure 9). In all simulations in which solute drag is considered grain sizes are considerably smaller than the ones predicted in simulations without impurities. In such simulations, for isothermal conditions, grain sizes reach 2 and 8 meters after 1 Ga at temperatures of 1173 and 1373K, respectively, while they do not exceed 47 cm for the solute drag simulation with the largest temperature (1373K) and the lowest impurity concentration (600 ppm).

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 conditions (temperature and annealing times) consistent with those inferred for the Udachnaya kimberlite xenoliths (i.e., temperatures between 1173 to 1273K and residence times between 100 and 1450 Ma). The Kaapvaal peridotites, which have smaller grain sizes and potentially larger residence times than the Udachnaya ones, range below the 1073K temperature curve. Thus, if the annealing times of these peridotites are, as inferred, > 1Ga their grain sizes cannot therefore be explained by solute drag only, even if the predictions of the simulations with solute drag are closer to the observed values. Table 2: Temperature and impurity concentrations (c 0 ) conditions for the mean field simulations and their results in terms of grain sizes at two different steps (10 and 1000 Ma).

Olivine grain sizes measured in the Kerguelen harzburgites (2-10 mm) may be simulated if annealing times at temperatures >1373K are smaller than 1Ma.

They are consistent with up to 38Ma of annealing at temperature ≤ 1173K.

The smaller olivine grain sizes in the Kerguelen dunites may only be explained by lower equilibration temperatures and annealing times (very short annealing at 1173K or up to 10Ma of annealing at 1073K).

Discussion

GG simulations, which explicitly consider solute drag by Ca impurities in olivine, adequately reproduce both grain size evolution in laboratory experiments and observed grain size ranges in mantle peridotites subjected to various annealing times in the lithospheric mantle, if solute drag parameters derived directly from theoretical considerations on the interactions between Ca impurities and olivine grain boundaries, in particular β, are increased by two orders of magnitude.

This difference could arise from two main reasons :

First, the two functions D(x) and E(x) (diffusion coefficient and interaction energy as function of distance to grain boundary x) needed for the calculation of solute drag parameters α and β (see section 2) are poorly known even for metallic materials in which the interaction between an impurity ion and a grain boundary is mostly elastic. For ceramics-like materials such as olivine, this interaction should also account for electrostatic interaction between the solute ion and the grain boundary and the interaction between solute-vacancy dipoles and electrical field around grain boundaries [START_REF] Yan | Space charge, elastic field, and dipole contributions to equilibrium solute segregation at interfaces[END_REF]. The quantification of all of those interactions should be done by atomistic calculations using systematic approaches for describing grain boundaries (e.g. [START_REF] Hirel | Systematic theoretical study of [001] symmetric tilt grain boundaries in mgo from 0 to 120 gpa[END_REF]) and their interactions with solute.

Secondly, it is well accepted that the impurity concentration at grain boundary evolves with the grain size, increasing when grain size increases [START_REF] Marquardt | The structure and composition of olivine grain boundaries: 40 years of studies, status and current developments[END_REF] which could also be expressed by a grain size dependent partition coefficient [START_REF] Hiraga | Grain boundaries as reservoirs of incompatible elements in the earth's mantle[END_REF]. Even if the solute drag pressure (Eq.2) depends on grain matrix impurity concentration, one can find expressions for α and β as functions of the partition coefficient [START_REF] Cahn | The impurity-drag effect in grain boundary motion[END_REF]. The effect of indirectly introducing an increase in the grain boundary impurity concentration with increasing grain size may be similar to an increase of β, as we did here. In fact, increasing the β value will increase the space of the high velocity regime for which the grain boundary migration will be poorly affected by impurities. An increase of the impurity concentration with grain size will also result in a lower impact of the solute drag for small grains. Moreover, it will also increase the domain of effect of solute drag and decrease GG rates of coarse grained polycrystals.

Our full field and mean field simulations show an initial phase of very slow grain size evolution (Figs. 6 and8) due to the direct slow down of the grain shrinkage and growth by solute drag, and indirectly by the impediment of GBM by the impurity slowed grain boundaries. This quasi-static phase, which is longer for coarser initial grain sizes, implies a very weak dependency of the grain growth kinetics on initial grain size. This small dependency frees us from the need of precise constrains on the initial grain size, which is very difficult to know in geological contexts. Taking advantage of this, we can apply our solute drag model for different stable geological contexts. For the majority of the geological contexts presented here, our formalism shows consistent grain size / time predictions (Fig. 9). It is difficult to evaluate precisely the performance of our models with respect to the evolution of natural samples, given the lack of data concerning their temperature and grain size evolution. However, impurity drag due to Ca concentration in olivine, in the range of commonly observed value, explains a grain size reduction of several orders of magnitude at geological timescales, compared to models without impurities.

In the simulations presented in this work, secondary phases are not taken into account, which could explain why some of our simulations applied to geological contexts overestimated grain sizes (particularly for Kaapvaal harzburgites, cf Fig. 9). Indeed, the presence of secondary phases such as pyroxenes and alumina phases in natural rocks is known to both impose, through Smith-Zener pinning mechanisms, a maximal grain size [START_REF] Tasaka | Influence of mineral fraction on the rheological properties of forsterite+ enstatite during grain size sensitive creep: 3. application of grain growth and flow laws on peridotite ultramylonite[END_REF]Nakakoji et al., 2018) (i.e. distance between second phase grains) and to reduce grain growth kinetics [START_REF] Furstoss | Full field and mean field modeling of grain growth in a multiphase material under dry conditions: application to peridotites[END_REF] of the rocks. In natural annealed peridotites, the distance and the size of second phase grains range between few mm to few cm.

This distance constrains the maximum mean grain size that may be attained in a polyphasic rock, making predictions of grain sizes of tens of cm and higher unrealistic. Other natural factors could also be taken into consideration for a holistic description of the mechanisms impacting grain growth. For example, it has been shown that the presence of melt and the water content in the material impact the GG within laboratory experiments by accelerating the growth kinetics [START_REF] Evans | A few remarks on the kinetics of static grain growth in rocks[END_REF]. This effect may be similar in a geological context in which metasomatism has been active and/or melt has percolated within the rock and/or water content is high.

Conclusions

In this work we have demonstrated that accounting for the presence of impurities within olivine rich-rocks permits to explain GG kinetics of both experimentally and naturally annealed olivine aggregates. The solute drag parameters in the simulations that allow for fitting both datasets are however quite different from those predicted for Ca, which is the impurity expected to have the highest impact on olivine grain boundary migration in mantle rocks, based on simple theoretical considerations. Atomistic calculations or a new experimental framework are needed to explain this gap.

We have developed a new mean field model accounting for the presence of those impurities and showed this approach could be used to predict grain size of olivine-rich rocks in geological contexts such as isothermal evolution. As this approach successfully reproduced natural grain size in annealed peridotite, healing kinetics may be implemented in large scale numerical geodynamic models based on this framework. In order to have a useful grain size evolution law for geodynamic model, a model accounting for the competition between grain size reduction due to deformation (through dynamic recrystallization phenomenon) and growth should be developed.

This approach is fast enough for implementation in large scale numerical geodynamic models and allows for better predictions of olivine grain sizes than the direct extrapolation of the experimentally-derived flow laws. However, to be consistent with the real multiphase nature of natural rocks, the influence of second phases, such as pyroxenes and Al-rich phases, on GG kinetics should also be considered. The present work shows that the latter plays in particular a major role in limiting the maximum grain size that may be achieved, whereas solute drag plays a very important role in the initial stages of the growth process, when olivine grain sizes are smaller than the second phase spacing. Such

a framework may open the door to a paleo-chronometer based on grain size evolution in geological contexts where deformation has stopped and static GG predominates.

Figure 1 :

 1 Figure 1: Solute drag pressure intensity as a function of the grain boundary velocity calculated using eq.2 for c 0 = 100ppm and different values of α and β.

Figure 3 :

 3 Figure 3: Example of GSD (in number) before and after a repopulation step with a minimal bin number of 40. Both histograms are constructed using 20 classes ranging between the radii of the bins with the minimal and maximal grain radius.

  sure exerted by Cr (r i = 0.062nm), Al (r i = 0.054nm) and Ca (r i = 0.1nm) for different grain boundary velocities and for natural characteristic c 0 concentrations of 200, 50 and 100ppm respectively[START_REF] De Hoog | Trace-element geochemistry of mantle olivine and application to mantle petrogenesis and geothermobarometry[END_REF]. Considering the diffusion coefficients of the three elements at 1573K (respectively, 1.4•10 -15 ,

Figure 4 :

 4 Figure4: Olivine aggregates GG kinetics at 1573K (red) and 1073K (blue). Thin solid lines represent the usual simulations not accounting for solute drag based on the direct extrapolation of the laboratory results of[START_REF] Karato | Grain growth kinetics in olivine aggregates[END_REF]) (yellow). The bold purple curve present the predictions from a full field simulation using calcium solute drag parameters α Ca and β Ca calculated in section 3, while the bold red and blue solid lines represent full field simulations with α = 10α Ca and β = 500β Ca at 1573K and 1073K, respectively . The dashed lines represent the GG kinetics computed using the mean field model with adjusted parameters based on the full field simulations at 1573K (red) and 1073K (blue). The green rectangle represents olivine mean grain sizes and annealing times typical of lithospheric mantle rocks (see section 4.2).

Figure 5 :

 5 Figure 5: Distribution of mass term M along grain boundaries for the first three runs at 1573K using α = 10α Ca and β = 500β Ca . The first run (left), a 0.2mm aside square, is also represented in the right bottom corner of the second run (center), a 1.6mm aside square, which is also represented at the right bottom corner of the third run (right), a 12mm aside square.

Figure 6 :

 6 Figure 6: Mean grain size evolution at 1573K for α = 10α Ca and β = 500β Ca for the first four steps of the simulation (without the time offset normally applied to get the whole grain growth kinetics presented in Fig.4), the full large lines represent the full field model accounting for impurity, the full tiny lines represent the mean field computations and the dashed lines represent the equivalent full field models without impurity. The black stars represent the two sampling points for the GSD presented in Fig.7.

Figure 8 :

 8 Figure 8: Grain size evolution for isothermal model with T = 1173K, c 0 = 800ppm and three initial grain sizes of 20, 500 and 2000µm.

Figure 9 :

 9 Figure 9: Grain size evolution for isothermal models with 4 different temperatures and 4 different impurity (c 0 ) concentrations. Solid squares indicate the approximate location of Siberia and Kerguelen and Kaapvaal xenoliths on this curve according to the literature (see references in Tab.1).

Acknowledgments

We greatly acknowledge Editor Juan Carlos Afonso and the two anonymous reviewers for their constructive comments. 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, SAFRAN through the DIGIMU Industrial Chair and consortium are gratefully acknowledged.

full field simulations. Finally we show an application of this framework within the geological contexts presented in section 4.2.

Full field and mean field simulations : adjustment of solute drag parameters

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 does not match the experimental results of [START_REF] Karato | Grain growth kinetics in olivine aggregates[END_REF]) (Fig. 4). The predicted GG (purple curve in Fig. 4) quickly deviates from the laboratory data (yellow curve); thus we stopped this simulation at the end of the second step in the model chaining strategy.

We tested different values for the solute drag parameters and found that α = 10α Ca and β = 500β Ca produce results consistent with laboratory results at 1573K and have GG kinetics compatible with natural observations at 1073K (Fig. 4). Using those values, the mean grain size evolution begins to deviate from the classical extrapolation of laboratory results (which does not account for solute drag) for a mean grain size near 30µm at 1573K and near 50µm at 1073K (Fig. 4). The GG kinetics is slowed down by solute drag for mean grain sizes of up to ca. 1mm. For coarser mean grain sizes, GG is weakly influenced by 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) clearly shows an increase of its value from the small grains to the large ones. At first, for ultra-fine grained aggregates as the ones used in laboratory experiments, 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 the grain growth kinetics is nearly unimpacted by the solute drag (Fig. 1). Afterwards, as the grain sizes increase the mass term becomes

Data availability

The data for supporting all figures of the paper are avalaible upon request to the authors as well as the mean field code used in this work.