Lucan Mameri
email: lucan.mameri@gm.univ-montp2.fr

Andrea Tommasi
email: andrea.tommasi@umontpellier.fr

Javier Signorelli

Lars Hestbjerg Hansen

Andréa Tommasi

Lars N Hansen

Andréa Tommasi

Predicting viscoplastic anisotropy in the upper mantle: a

Keywords: olivine, crystal preferred orientation, mechanical anisotropy, dislocation creep, mantle deformation, rheology

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The lithospheric mantle is predominantly composed of olivine, a mineral that displays highly anisotropic elastic and viscoplastic mechanical behavior [START_REF] Kumazawa | Elastic moduli, pressure derivatives, and temperature derivatives of single-crystal olivine and single-crystal forsterite[END_REF][START_REF] Bai | High-temperature creep of olivine single crystals 1. Mechanical results for buffered samples[END_REF]. Upper mantle deformation by viscoplastic processes produces crystallographic preferred orientation (CPO or texture) of olivine [START_REF] Nicolas | Formation of anisotropy in upper mantle peridotites: a review, in Composition, Structure and Dynamics of the Lithosphere-Asthenosphere System[END_REF]. These CPO result in anisotropic behavior at the scale of the rock (cm-to m-scale) and larger, as indicated by observations of coherent seismic anisotropy patterns correlated with the tectonic fabric of continental plates at the scale of hundreds of km (cf. review in [START_REF] Tommasi | Heterogeneity and anisotropy in the lithospheric mantle[END_REF]. Large-scale olivine CPO in the upper mantle has therefore been proposed as a source of structural memory in the upper mantle, leading to strain localization along shear zones formed during previous tectonic events [START_REF] Vauchez | Rheological heterogeneity, mechanical anisotropy and deformation of the continental lithosphere[END_REF][START_REF] Tommasi | Continental rifting parallel to ancient collisional belts: an effect of the mechanical anisotropy of the lithospheric mantle[END_REF].

The consequences of upper mantle viscoplastic anisotropy on the deformation of both the plates and convective mantle have been analyzed in a variety of numerical models. Many of these models used a simplified description of this anisotropy by imposing a constant difference between 'shear' and 'normal' viscosities and reorienting the viscosity tensor in response to ow based on a director vector approach (e.g., [START_REF] Honda | Strong anisotropic flow in a finely layered asthenosphere[END_REF][START_REF] Christensen | Some geodynamical effects of anisotropic viscosity[END_REF]Lev andHager, 2008, 2011). These models, which tested normal to shear viscosity ratios of 10-100, concluded that viscous anisotropy does change mantle flow patterns by stabilizing the boundary layers and enhancing strain localization. Most recent models have coupled polycrystal plasticity simulations into finite-element models; this allows the follow-up of the evolution of the orientation and intensity of the full viscosity tensor (e.g., [START_REF] Chastel | Anisotropic convection with implications for the upper mantle[END_REF][START_REF] Knoll | A multiscale approach to model the anisotropic deformation of lithospheric plates[END_REF][START_REF] Tommasi | Structural reactivation in plate tectonics controlled by olivine crystal anisotropy[END_REF][START_REF] Castelnau | Numerical simulations of texture development and associated rheological anisotropy in regions of complex mantle flow[END_REF][START_REF] Blackman | Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy[END_REF]. These more sophisticated models support the early conclusions on the role of viscous anisotropy on the stability of convection [START_REF] Chastel | Anisotropic convection with implications for the upper mantle[END_REF], and corroborate that the olivine CPO formed during preexisting deformation episodes and subsequently frozen into the lithospheric mantle may explain the reactivation of ancient continental-scale structures to form new plate boundaries [START_REF] Knoll | A multiscale approach to model the anisotropic deformation of lithospheric plates[END_REF][START_REF] Tommasi | Structural reactivation in plate tectonics controlled by olivine crystal anisotropy[END_REF] and change the mantle flow pattern beneath spreading ridges [START_REF] Castelnau | Numerical simulations of texture development and associated rheological anisotropy in regions of complex mantle flow[END_REF][START_REF] Blackman | Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy[END_REF]. However, at the time most of those calculations were carried out, no experimental data existed for corroborating the predictions of the polycrystal plasticity models.

Recently, laboratory experiments demonstrated that olivine aggregates with textures formed in torsion are ∼2 times weaker for further deformation in torsion than for extension perpendicular to the previous shear plane (Hansen et al., 2012a). These experiments provided the first quantitative data for validating the predictions of viscoplastic anisotropy in olivine-rich rocks by polycrystal plasticity models. This validation is important because the use of geodynamical models to study the effects of an anisotropic rheology on Earth dynamics depends directly on our ability to predict the anisotropy of viscoplastic behavior associated with any combination of CPO and loading geometry. Only polycrystal plasticity models can test in a physically consistent way the full range of orientation relations between applied stresses and the broad range of texture types observed in upper-mantle rocks.

In this work, we test how well different polycrystal plasticity models simulate the viscoplastic anisotropy of olivine polycrystals measured in the laboratory by calculating the instantaneous mechanical response in simple shear and axial extension for polycrystals with the same CPO as the experimental samples and comparing these predictions to (1) the final stresses measured in the torsion experiments that produced the CPO and (2) to the peak stresses measured in axial extension experiments subsequent to torsion. The comparison between the polycrystal plasticity models and the experimental data allows us to quantify the contribution of other processes to the deformation of olivine polycrystals under laboratory conditions and propose modeling approaches that predict the possible range of texture-induced viscoplastic anisotropy in upper mantle domains deforming by dislocation creep.

Methods

Modeling viscoplastic behavior and texture evolution

We tested the capability of different polycrystal plasticity models to predict the viscoplastic behavior of textured olivine-rich rocks at high temperature and moderate pressure by simulating the torsion and ex-tension tests of Hansen et al. (2012a). Hansen et al. (2016a) presents complementary results (extension followed by torsion), but a heterogeneous distribution of strain in the initial extension experiments resulted in large variability among textures and stresses. For testing the different polycrystal plasticity models, we require higher precision data. Therefore in the present study, we focus solely on simulating experiments from Hansen et al. (2012a) in which torsion was followed by extension. To do so, we calculated the instantaneous mechanical response of four olivine polycrystals that have precisely the same CPO as those measured in the experimental samples. We also applied the same loading geometries as in the experiments: simple shear parallel to the previous shear deformation that created the texture or axial extension parallel to the previous torsional axis (Fig. 1). The employed models are viscoplastic and include no elasticity. This simplification is justified by fullfield simulations of polycrystalline materials performed in the elasto-viscoplastic regime that demonstrate that the viscoplastic behavior controls the mechanical response of the polycrystal even at very small strains (< 1%, [START_REF] Lebensohn | An elastoviscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF]. Hansen et al. (2012a). Instantaneous stresses predicted by the models in simple shear and extension are compared to (a) final stresses in torsion and (c) peak stresses in extension, respectively. The stresses in extension are averaged values normalized to a constant strain-rate of 2 × 10 -4 s -1 . The number of grains in each CPO is indicated by the value n, and the CPO intensity by the value J (J-index). Pole figures are contoured as multiples of uniform distribution; the same contour range limits are applied to all pole figures for easy evaluation of variations in CPO strength.

The viscoplastic mechanical response of the olivine polycrystals to the imposed extensional or simple shear deformation has been modeled using three polycrystal plasticity models: the uniform stress (lower bound) approach [START_REF] Sachs | Zur Ableitung einer Fliessbedingung[END_REF][START_REF] Chastel | Anisotropic convection with implications for the upper mantle[END_REF], the tangent VPSC model [START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys[END_REF], and the more recent second-order formulation of the VPSC model [START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: Itheory[END_REF][START_REF] Lebensohn | Full-field vs. homogenization methods to predict microstructure-property relations for polycrystalline materials[END_REF]. The three approaches differ in the way the crystal interactions are modeled. The lower bound assumes constant stress in all grains and disregards strain compatibility. In the VPSC approach, each grain is represented as a viscoplastic inclusion embedded in and interacting with a homogeneous viscoplastic medium that represents the average polycrystal behavior. Stress equilibrium and strain compatibility at the aggregate scale are obtained by approximating the macroscopic strain rate and stress to the average of strain rate and stress at the grain scale, respectively. The second-order (SO)-and tangent (TGT)-VPSC models differ in how the non-linear behavior of the material is estimated. Both formulations are based on the same linearization scheme, which considers heterogeneous stress and strain rate within the polycrystal. However, only the second-order formulation accounts for intragranular strain rate fluctuations to predict the effective polycrystal behavior [START_REF] Ponte Castañeda | Second-order homogenization estimates for nonlinear composites incorporating field fluctuations: Itheory[END_REF][START_REF] Lebensohn | Full-field vs. homogenization methods to predict microstructure-property relations for polycrystalline materials[END_REF].

In all three approaches, individual crystals deform only by dislocation glide on a finite number of slip systems, which are a function of the crystal structure. The relative strength of these slip systems, that is, the critical resolved shear stress (CRSS) needed to activate dislocation glide on each system, depends on the temperature, pressure, and stress. For technical reasons, the experiments of Hansen et al. (2012a) were performed using olivine crystals enriched in iron (Fo 50 ) relative to typical mantle olivines (Fo 90 ), for which mechanical data are not available for creep of single crystals. In the present study, we use CRSSs (Table 1) derived from deformation of natural olivine single crystals (Fo 90 ) at high temperatures and low pressure [START_REF] Bai | High-temperature creep of olivine single crystals 1. Mechanical results for buffered samples[END_REF]. The choice of these CRSSs, although arbitrary, is justified by the analysis of the CPO produced in the torsion experiments of Hansen et al. (2012a), which is consistent with the dominant activation of the (010)[100] slip system. We also examined the effect of setting equal the values of CRSS for the (001)[100] and (010)[100] slip systems, as used in previous VPSC models studying olivine CPO evolution at high temperature and low pressure (e.g., [START_REF] Tommasi | Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations: implications for the upper mantle seismic anisotropy[END_REF], and we observed negligible variations in terms of mechanical behavior compared to the present model setup. Laboratory experiments and atomistic modeling suggest nevertheless a change in the relative resistance between the different slip systems, favoring glide of dislocations with [001] Burgers vectors relative to those with [100] Burgers vectors at high pressure (e.g., [START_REF] Durinck | Application of the Peierls-Nabarro model to dislocations in forsterite[END_REF][START_REF] Raterron | Pressure-induced slip-system transition in forsterite: single-crystal rheological properties at mantle pressure and temperature[END_REF]. Comparison of the seismic anisotropy predicted by polycrystal plasticity models with high and low pressure CRSS sets with seismic anisotropy observations suggests that the CRSS data used in the present study are representative of deformation in the upper 200-250 km of the mantle (Mainprice et al., 2005). The shear strain rate in each slip system is proportional to the resolved shear stress on the slip plane in the slip direction normalized to the CRSS of the slip system raised to the power n. The iron-rich olivine polycrystals display a stress exponent n of 4.1, whereas typical mantle olivine has lower n values (e.g., 3.5, [START_REF] Bai | High-temperature creep of olivine single crystals 1. Mechanical results for buffered samples[END_REF]. We tested n equal to 3 or 4 and observed little difference in terms of magnitude of anisotropy (cf. supplementary material, Fig. S1).

In the following we focus on the results for n = 3 in order to agree with most laboratory and previous numerical experiments on olivine polycrystals, which usually find stress exponent n = 3 ± 0.5.

A major challenge in modeling viscoplastic anisotropy of a rock (i.e., an aggregate of crystals with different orientations) is that traditional polycrystal plasticity models only account for deformation accom-modated by dislocation glide. Since olivine has fewer than 4 independent slip systems, deformation processes other than glide are required to ensure strain compatibility. To simulate the contribution of these other processes, [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF] modified the SO-VPSC model by adding an isotropic deformation process. In the present study we build on this already validated development. To close the yield surface and simulate the contribution of isotropic deformation processes to plastic deformation, we add three fictitious pyramidal slip systems in the VPSC models:{111}〈110〉,{111}〈011〉and {111} 〈 101 〉 , here referred as dummy slip systems. These slip systems were defined such that they may accommodate any imposed plastic strain, but do not contribute to plastic spin. They mimic therefore isotropic strain accommodation mechanisms, which do not contribute to CPO evolution. The contribution of the dummy systems to the total viscoplastic deformation is varied by changing their CRSS by a factor β (Table 1). Thus, the present implementation differs from previous VPSC model applied to olivine polycrystals, referred thereafter as 'standard' VPSC model, since in the latter the dummy slip systems contributed to plastic spin and were not allowed to accommodate > 5% (therefore contributing only in a minor way to the CPO evolution) by using β values ≥50 (e.g., [START_REF] Wenk | Modelling plastic deformation of peridotite with the self-consistent theory[END_REF][START_REF] Tommasi | Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations: implications for the upper mantle seismic anisotropy[END_REF][START_REF] Castelnau | Numerical simulations of texture development and associated rheological anisotropy in regions of complex mantle flow[END_REF].

Hardening is not considered in the present models since in most high-temperature ( 1200 °C) deformation experiments on olivine aggregates, in the absence of significant microstructural changes, steady-state mechanical behavior is achieved after relatively small strains (e.g., [START_REF] Chopra | The role of water in the deformation of dunite[END_REF]Hansen et al., 2012b;[START_REF] Thieme | Stress evolution and associated microstructure during transient creep of olivine at 1000-1200 °C[END_REF]. In the present models, diffusive processes are therefore implicitly ac-counted for in two ways: via their contribution to avoid hardening by dislocation entanglements and via the activation of the dummy slip systems. Indeed, although we do not model any specific additional process, but solely their contribution to strain accommodation in the olivine polycrystal, most possible additional processes, such as recrystallization, grain boundary migration or sliding, involve diffusion.

Model setup

Two series of tests were run to simulate the stresses measured in the laboratory: (1) simple shear parallel to a previous shearing deformation that produced the olivine CPO, that is, shearing in the X direction and in a plane normal to the Z direction; (2) axial extension normal to the previous shear plane, parallel to the Z direction (Fig. 1). Lower bound simulations were run under kinematic boundary conditions by im-posing all components of the velocity gradient tensor L (s -1 ). For easy comparison between simple shear (ss) and axial extension (ext) simulations, L is defined with the condition that the norm of the strain rate tensor in both cases equals 1 (cf. supplementary material):

In the VPSC simulations, mixed boundary conditions were enforced because they better mimic the actual conditions in laboratory experiments: imposed torsion on the plane normal to the piston axis or imposed axial extension and homogeneous stress in the plane normal to the piston axis.

*=not imposed

We also ran VPSC simulations using fully kinematic boundary conditions. Differences between stresses predicted by TGT or SO-VPSC simulations using either kinematic or mixed boundary conditions are low (< 10%, cf. supplementary material, Fig. S2), by consequence we focus on the results for mixed boundary conditions.

The models are run for a single deformation step of 1% Von Mises equivalent strain and the mechanical behavior is analyzed before actualization of the texture, in other words, we calculated the instantaneous mechanical response of the olivine polycrystals. The stresses predicted by the viscoplastic models in simple shear and in axial extension were compared to final stresses measured in the torsion experiments and to peak stresses in the extensional ones, respectively (Fig. 1). For internal consistency, stresses measured in the laboratory were normalized for a homogeneous Von Mises equivalent strain rate of 2 × 10 -4 s -1 (Table 2; the normalization procedure is presented in the Supplementary Material). Stresses calculated in the models are normalized by the strength of the easiest slip system, (010) [100]. In addition, for comparing the viscoplastic anisotropy observed in the laboratory to the one predicted by the different numerical models, Von Mises equivalent stresses for simple shear and axial extension were normalized by the Von Mises equivalent stress of the isotropic aggregate in each data set (experiments or different modeling approaches). The effect of the texture on the viscoplastic mechanical behavior was tested using four different initial CPOs: an untextured and three olivine CPOs measured at the end of the torsional tests by Hansen et al. (2012a), performed at high temperature and moderate pressure conditions (Fig. 1). These initial textures are representative of common CPO symmetries observed in natural and experimental samples deformed in simple shear and cover the entire range of CPO intensities measured in peridotite massifs and mantle xenoliths [START_REF] Tommasi | Heterogeneity and anisotropy in the lithospheric mantle[END_REF]. In the models, grains are initially spherical. The untextured polycrystal is made of 1000 grains. The CPO intensity is characterized by the J-index, which is the integral of the squared orientation distribution function (ODF), quantified using the MTEX toolbox based on an ODF calculated using a "de la Vallée Poussin" kernel with a half width of 10° [START_REF] Mainprice | Descriptive tools for the analysis of texture projects with large datasets using MTEX-Strength, symmetry and components[END_REF]. Sample PT0655, which was initially deformed to a shear strain of 1, is represented by a texture composed of 870 grains. This texture is very weak (J-index of 1.5); it is characterized by a subtle double maximum of [100] at low and high angles to the shear direction and by a dispersed girdle of [010] at high angle to the maximum finite elongation with a weak maximum normal to the foliation. Samples PT0640 and PT0651, which underwent initial shear strains of 4 and 10, are represented by textures composed of 489 and 665 grains, respectively. Both have very strong olivine CPOs (J-index of 14.5 and 20.5, respec-tively) and are characterized by alignment of the [100] maximum with the shear direction and of the [010] maximum normal to the shear plane. However, the textures differ by the distribution of the [010] axes, which forms a girdle normal to the shear direction (X) in PT640 and a point maximum normal to the shear plane (XY plane) in PT655 (Fig. 1).

To verify if the modified SO-VPSC approach used in this work is capable of reproducing the CPO evolution recorded in simple shear deformation, we run simple shear tests up to shear strains of 20 for an initially untextured olivine polycrystal composed of 1000 grains. The simulations were conducted applying the same slip system parameters and (mixed) boundary conditions used to model the viscoplastic mechanical behavior. This study focuses on the SO-VPSC predictions because they have been shown to better reproduce the solutions of full-field simulations [START_REF] Castelnau | Micromechanical modeling of the viscoplastic behavior of olivine[END_REF]. The predictions of the TGT-VPSC model are nevertheless presented for continuity with previous work (e.g., [START_REF] Wenk | Large strain shearing of halite: experimental and theoretical evidence for dynamic texture changes[END_REF][START_REF] Tommasi | Upper mantle tectonics: Three-dimensional deformation, olivine crystallographic fabrics and seismic properties[END_REF][START_REF] Tommasi | Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations: implications for the upper mantle seismic anisotropy[END_REF].

Results

Mechanical behavior

Fig. 2a presents a comparison between the stresses normalized to the strength of the untextured polycrystal predicted by the different viscoplastic models (colored markers) at standard conditions (CRSS of the dummy slip systems in TGT and SO-VPSC models, β = 50) and the experimental measurements (black circles), as a function of the intensity of the initial CPO. Both laboratory and numerical experiments show that textured olivine polycrystals deformed in simple shear parallel to the preexisting fabric (empty markers) are weaker than the untextured aggregate, whereas those deformed in axial extension normal to the preexisting fabric (filled markers) are harder. There is a positive correlation between CPO intensity and magnitude of viscoplastic anisotropy, defined as the ratio between the Von Mises equivalent stresses in extension and in simple shear tests, but the relationship is not linear (Fig. 2b). In simple shear, the lower bound and TGT-VPSC models predict normalized stresses in reasonable agreement with the experimental data, but SO-VPSC simulations predict lower normalized stresses. In axial extension, the models overestimate the strength of the olivine polycrystals relative to the experimental measurements, except for the sample PT0655, which despite a very weak CPO, displays a stronger anisotropy than the one predicted by the models. The predictions of the lower bound approach are the closest to the experimental data and the TGT-VPSC models are the farthest. The discrepancy between numerical and laboratory experiments, which increases with increasing intensity of the CPO, results therefore from overprediction of the strength of olivine polycrystals deformed numerically in extension. In this loading geometry, most crystals are in hard orientations relative to extension axis, that is, low resolved shear stresses are expected on the three main slip systems of olivine: The lower bound model, which only enforces stress equilibrium and not strain compatibility, predicts the lowest polycrystal strengths for all textures and loading geometries. The better fit of the experimental data by the lower bound approach points to the major role of enforcing strain compatibility in modeling the mechanical behavior of olivine polycrystals. Further analysis of the stresses predicted by the different models (Fig. 3; Table 2) corroborates this point. The SO-VPSC, which most strictly enforces strain compatibility and better fits fullfield polycrystal plasticity simulations [START_REF] Castelnau | The effect of strain heterogeneity on the work hardening of polycrystals predicted by mean-field approaches[END_REF], predicts the highest stress values for both simple shear and extension. The high stresses predicted in axial extension by the SO-VPSC model for the strongly textured polycrystals PT0640 and PT0651 result from low activities of the (010)[100] and (001)[100] slip systems, which have the lowest CRSSs (Fig. 4). For these textures, resolved shear stresses on these two slip systems are low in most crystals (i.e., extension is applied parallel to the maximum concentration of [010] axes and normal to the maximum concentrations of [100] and [001] axes). To accommodate the imposed extensional deformation, the very hard (β = 50) dummy slip systems have to be activated in many crystals. Even though these systems contribute to < 5% of the total deformation (Fig. 4), this moderate amount of activation still requires high stresses. The (010)[001] slip system, which has intermediate CRSS values, contributes significantly to the total deformation (∼55%) in axial extension. In most crystals, resolved shear stresses on the (010)[001] slip system are also low, but the high stresses attained in the extensional deformation allow their activation.

In contrast, by imposing simple shear parallel to the maximum concentration of [100] olivine axes in a plane parallel to the maximum concentration of (010) planes in the strongly textured polycrystals PT0640 and PT0651, we obtain maximum resolved shear stresses on the easy slip system (010)[100], which accommodates > 50% of the deformation. In both simple shear and axial extension, the activity of the (010)[100] slip system is higher in the more strongly textured polycrystal (PT0651) than in the polycrystal with a CPO of intermediate strength. A closer analysis of the slip system activities highlights nevertheless some subtle differences between predictions for samples PT0640 and PT0651 (Fig. 4). This variation may be explained by differences in the symmetry of the CPO between the two polycrystals. PT0640 (J-index of 14.1) has a girdle distribution of [010] and [001], which allows in axial extension for high resolved shear stresses on (010)[001] in a larger number of crystals relative to PT0651 (J-index of 20), which has point distributions of the three main axes. The observed discrepancy between laboratory results and model predictions advocates in favor of a significant contribution of others processes in addition to dislocation glide during deformation of olivine aggregates at laboratory conditions. This contribution appears to be particularly important if most crystals in the aggregate are unfavorably oriented to deform by dislocation glide on the easy slip systems, such as the present extensional tests on the strongly textured polycrystals. To evaluate if an increase of the contribution of isotropic deformation processes allows reproducing the stresses measured experimentally, the CRSS of the dummy slip systems was systematically decreased. The results are presented in the following section.

Strain accommodation by additional isotropic deformation processes

We tested the effect of increasing the contribution of isotropic deformation processes on the aggregate strength by decreasing the CRSS of the dummy slip systems (β = 10, 5, 3, 2, 1) in the SO-VPSC simulations. These dummy systems represent any additional process that allows olivine crystals to deform even if oriented in a manner that results in very low resolved shear stresses on the three main slip systems: ( 010 (Hansen et al., 2012a). For reference, the mechanical behavior of an untextured polycrystal (isotropic) is displayed by a thin dashed (axial extension) or solid (simple shear) lines in all three figures; for the color legend of the pole figures and loading geometries, refer to Fig. 1.

A decrease in β results in a non-linear decrease in the polycrystal strength for all textures and both loading geometries (Fig. 5ac). The rate of decrease varies as a function of CPO intensity and of loading geometry. The weakly textured aggregate deforms at similar stresses in both loading geometries (i.e., very weak anisotropic behavior) even at high values of β. In contrast, for the strongly textured aggregate (Fig. 5b,c), decreasing β leads to distinct behaviors for the two loading geometries: the reduction in the polycrystal strength is more marked in the axial extension test than in the simple shear simulations. It results therefore in reduction of the magnitude of anisotropy with decreasing β (Fig. 6a). The modified SO-VPSC simulation in which the dummy systems (that allow for isotropic deformation) are 3 times harder to activate than the (010)[100] and (001)[100] easy slip systems (β=3) better reproduces the mechanical behavior measured in the laboratory; except for the sample with the weakest CPO (PT0655), for which all models underestimate the magnitude of anisotropy (Fig. 6a,b). Differences in mechanical predictions between TGT-and SO-VPSC models are only observed for β > 10 (cf. supplementary material, Fig. S3). The change in mechanical behavior with decreasing β may be explained by variations in the relative activities of the actual olivine slip systems and of the dummy ones, which represent additional deformation processes activated to ensure strain compatibility constraints (Fig. 7).

In the VPSC model with β = 50, the three main slip systems (010)[100], (001)[100], and (010)[001] account for at least 95% of the viscoplastic deformation in axial extension and 99% of the deformation in simple shear. As β decreases, the activity of the dummy slip systems increases at the expense of the activity of the actual ones. In a untextured aggregate (Fig. 7a), reduction of β from 50 to 10 results in an increase of the activity of the dummy slip systems from < 1% to ∼11%, for both simple shear and axial extension. For β = 5, their activity increases to ∼30%. A further decrease in β to 3 or 2, which is the range of values that better agrees with the mechanical data measured in the laboratory, leads to activities of the dummy slip systems of ∼45% or ∼60%, respectively. For β = 1, > 80% of the deformation is accommodated by the dummy slip systems. These high activities of the relaxation mechanism, which is represented in the present study by the dummy slip systems, are similar to those previously observed by [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF].

The presence of an olivine CPO results in different relative slip system activities in axial extension and simple shear (Fig. 7b-d). This difference increases with increasing CPO intensity, because in a strongly textured aggregate, a larger proportion of the crystals will be either well or poorly oriented to accommodate the imposed deformation. The three olivine aggregates used in the present simulations have textures in which most crystals are well oriented to accommodate the imposed simple shear, but poorly oriented to accommodate the imposed axial extension. For β equal to 5 or 3, the isotropic mechanism in simple shear respectively accommodates ∼25% or ∼40% (PT0655), ∼15% or ∼25% (PT0640), and ∼10% or ∼ 20% (PT0651) of the total viscoplastic deformation, whereas in axial extension it accommodates ∼30% or ∼50% (PT0655), ∼50% or ∼70% (PT0640), and ∼60% or ∼75% (PT0651), (Fig. 7b-d). It is noteworthy that for β ≤ 5, the dummy slip systems remain active even if most crystals in the aggregates are well oriented to deform by dislocation glide on the easy (010)[100] and (010)[001] slip systems, such as when simple shear is applied on the strongly textured aggregates (PT0640 and PT0651). Their isotropic character may explain this behavior; there is always at least one system among the 12 systems that compose the three dummy slip modes on which the resolved shear stresses are high.

Olivine CPO evolution

To test if the modified SO-VPSC model that better fits the viscoplastic anisotropy measured in the experiments can also reproduce the texture evolution in simple shear, we compare the CPO evolution predicted in simulations with different CRSSs for the dummy slip systems (β=50, 5, or 3) to CPO data for olivine polycrystals deformed experimentally in simple shear at high temperature conditions [START_REF] Zhang | Lattice preferred orientation of olivine aggregates deformed in simple shear[END_REF][START_REF] Bystricky | High shear strain of olivine aggregates: rheological and seismic consequences[END_REF][START_REF] Demouchy | Deformation of olivine in torsion under hydrous conditions[END_REF]Hansen et al., 2012aHansen et al., ,b,c, 2014)). Fig. 8 displays the evolution of the CPO intensity characterized by the J-index and the angle between the maximum eigenvector of the [100] axis distribution and the shear direction as a function of the shear strain. Note that the CPO evolution of the olivine aggregates presented by Hansen et al. (2012a), for which the mechanical responses were simulated in the present study, is also included in Fig. 8. In this figure, we also compare the CPO evolution predicted by the present modified SO-VPSC model to those of a VPSC model modified to account for dynamic recrystallization by subgrain rotation in olivine aggregates (CDRX-TG-VPSC; for details on the model cf. [START_REF] Signorelli | Modeling the effect of subgrain rotation recrystallization on the evolution of olivine crystal preferred orientations in simple shear[END_REF]. The CRSSs in this previous model are slightly different from those used in the present study, but the relative strengths between slip systems are maintained. Thus, we do not expect significant differences in CPO evolution. Pole figures displaying the evolution of the olivine CPO in the SO-VPSC model with β = 50, 5, or 3 are shown in the Supplementary material (Fig. S5).

The SO-VPSC simulations in which the strength of the dummy slip systems (β) is equal to 50 show a rotation of the [100] axes towards the shear direction that closely follows the rotation of the finite strain ellipsoid, poorly reproducing the experimental data. A decrease in β leads to faster rotation of [100] maximum towards the shear direction, which better reproduces the CPO reorientation trend observed in the experiments. However, for low strengths of the dummy slip systems (β ≤ 3), that is, when the additional deformation process are very active, and at high shear strains (≥5), the CPO evolution departs from a continuous trend: the maximum concentration of [100] oscillates between 0° and 10° from the shear direction (Fig. 8). Models with β = 5 also exhibit a similar behavior, but the oscillations start at higher shear strains (∼10) and are smaller (≥5°). Fig. 8. CPO evolution as a function of increasing shear strain predicted by SO-VPSC simulations with different dummy slip system strengths (β). a) CPO orientation defined by the angle between shear direction and the maximum eigenvector of the [100]axis orientation distribution. b) CPO intensity evolution characterized by the J-index. Experimental data from [START_REF] Zhang | Lattice preferred orientation of olivine aggregates deformed in simple shear[END_REF], [START_REF] Bystricky | High shear strain of olivine aggregates: rheological and seismic consequences[END_REF], [START_REF] Demouchy | Deformation of olivine in torsion under hydrous conditions[END_REF]), and Hansen et al. (2012a,b,c,2014) are displayed for comparison. The CPO evolution of the olivine aggregates in Hansen et al. (2012a), in which the viscoplastic response were simulated in the present study, is highlighted as larger left-pointing triangles. The color gradient in the experimental data symbols indicate the average grain sizes. The predictions of the CDRX-TG-VPSC model by [START_REF] Signorelli | Modeling the effect of subgrain rotation recrystallization on the evolution of olivine crystal preferred orientations in simple shear[END_REF], which simulates recrystallization by subgrain rotation, are shown by the dotted-dashed line. The pale gray rectangle in (b) indicates the range of J-index measured in > 600 natural samples [START_REF] Tommasi | Heterogeneity and anisotropy in the lithospheric mantle[END_REF]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Analysis of the evolution of the CPO intensity shows that the decrease in the strength of the dummy slip systems results in stabilization of the CPO at lower J-index values (Fig. 8b). This observation highlights the essential role of additional processes, such as diffusion or dynamic recrystallization, on the evolution of the CPO intensity. CPO intensities measured in the high shear strain experiments are significantly overestimated for VPSC simulations with β=50 and underestimated for VPSC simulations with β=3, but they are well reproduced by models with β = 5. As for the CPO orientation, periodic oscillation of the CPO intensity is observed for shear strains ≥5 for β = 3 (Fig. 8). In Fig. 8, the experimental data is colored as a function of the average grain size. CPO intensity saturates at lower values in the finer grain size aggregates, being better fitted by models with low β values. This observation provides indication that the additional deformation process might be grain-size sensitive.

In conclusion, the modified SO-VPSC model with β = 5 better reproduces the CPO evolution in simple shear or torsion experiments than standard SO-VPSC simulations in which most of the deformation is accommodated by dislocation glide on the actual olivine slip systems (β = 50). CPO evolution in the modified SO-VSPC model with β = 5 is similar, though less stable at shear strains ≥10, to the one predicted by the VPSC model modified to account for dynamic recrystallization by subgrain rotation (CDRX-TG-VPSC model; [START_REF] Signorelli | Modeling the effect of subgrain rotation recrystallization on the evolution of olivine crystal preferred orientations in simple shear[END_REF]. The modified SO-VPSC model with β = 3, which better reproduces the viscoplastic anisotropy measured in the experiments, predicts well the evolution of the CPO orientation observed in simple shear or torsion experiments up to shear strains of ∼5, but it underestimates the CPO strength and diverges from the experimental data at higher shear strains.

Discussion

Viscoplastic anisotropy under laboratory conditions

All polycrystal plasticity models analyzed in the present study reproduce qualitatively the mechanical results of Hansen et al. (2012a): simple shear in a direction parallel to the [100] maximum on a plane normal to the [010] maximum requires lower stresses than axial extension parallel to the [010] maximum (Fig. 2a). This behavior may be explained by simple geometrical considerations. In simple shear, re-solved shear stresses on the easy slip systems, (010)[100] and (001) [100], are high for most olivine crystals composing the aggregate, whereas if axial extension is imposed parallel to the [010] maximum resolved shear stresses are very low on these easy slip systems for most crystals.

Despite the qualitatively coherent results, the polycrystal plasticity models overestimate the magnitude of viscoplastic anisotropy relative to laboratory experiments for polycrystals with a well-developed CPO (Fig. 2b). Based on the better agreement between the lower bound model predictions and the experimental data, we infer that the discrepancy between numerical models and experiments arises from imposing strain compatibility in a system that does not have enough degrees of freedom (< 4 independent slip systems). In the experiments, other strain accommodation mechanisms are active (recrystallization, grain boundary sliding or migration, diffusional mass transfer, climb, formation of kinks...). Some of these processes are isotropic, others are not. The fact that the viscous anisotropy measured in the laboratory is lower than the one predicted by VPSC simulations in which most of the deformation is accommodated by dislocation glide (β = 50) implies nevertheless that these additional processes play an important role in strain accommodation in the experiments. To test this hypothesis, we reduced the critical resolved shear stress of the dummy slip systems, which simulate the contribution of isotropic deformation processes to strain accommodation in the VPSC simulations. For strongly textured aggregates, models in which the strength of the dummy slip systems is 3 times higher than the strength of the easiest (010)[100] slip system in olivine do reproduce the magnitude of anisotropy measured in the laboratory experiments (Fig. 6). In these models, the dummy systems contribute more to the total deformation when the polycrystal is in hard orientation than when it is in a soft one (Fig. 7), decreasing the magnitude of anisotropy. The present results are consistent with the previous analysis by Hansen et al. (2016b), in which the mechanical data from the experiments simulated here together with other experiments on olivine aggregates were fitted by a pseudo-Taylor model that relaxed the strain compatibility constraint.

The SO-VPSC model modified to account for additional isotropic deformation processes, which produce plastic strain but no plastic spin, also better reproduces the CPO evolution measured in torsion experiments for β = 5 compared to predictions of the 'standard' VSPC model (β = 50). Higher contribution of the isotropic processes, which is associated with low β values, results in faster rotation of [100] maximum relative to shear plane and in slower evolution of the CPO intensity (Fig. 8). These simulations reveal a CPO evolution in simple shear similar to the one predicted by the VPSC model modified to account for dynamic recrystallization by subgrain rotation, in which the strain compatibility was relaxed within subgrains [START_REF] Signorelli | Modeling the effect of subgrain rotation recrystallization on the evolution of olivine crystal preferred orientations in simple shear[END_REF]. Both models therefore hint at the essential role of additional processes, which accommodate strain components that cannot be produced by dislocation glide in the deformation of olivine polycrystals. The modified VPSC model presented here therefore represents an acceptable alternative to predict olivine CPO evolution in simple shear, with a significantly lower calculation cost than VPSC models explicitly incorporating recrystallization.

The modified SO-VPSC simulations with β = 3, which better fit the viscoplastic anisotropy measured in the experiments, predict correctly the texture evolution observed in simple shear experiments up to shear strains of 5, but diverge from the experimental data by displaying an oscillatory behavior at higher shear strains. This oscillation is consistent with the fact that ellipsoidal inclusions embedded in a deformable media do not have any stable position in simple shear, but rotate continuously [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]. Indeed, low β values result in a low strength of the equivalent homogeneous medium that represents the polycrystal in VPSC. This interpretation is corroborated by the analysis of the SO-VPSC model in which the critical (maximum) aspect ratio of the crystals is changed; higher critical aspect ratios result in higher oscillations of the CPO orientation (cf. supplementary material, Fig. S4).

Comparison of the VPSC predictions to the experimental data suggests that although isotropic, the additional strain accommodation processes are associated with a non-linear mechanical behavior. The study of [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF], which first introduced additional isotropic deformation processes into the VPSC model, showed that significant activation of a Newtonian isotropic process resulted in a significant decrease of the apparent stress exponent of the polycrystal relative to the stress exponent associated with dislocation creep. The present study suggests that the mechanical data of Hansen et al. (2012a) are better reproduced by simulations in which the isotropic mechanisms accommodate a substantial proportion of the total viscoplastic deformation (≥20% in simple shear and ≥50% in axial extension). Such a significant contribution from the isotropic processes would surely decrease the stress exponent if these processes were Newtonian. Yet, stress exponents measured in torsional experiments at high finite strains on polycrystals comprised by olivine Fo 50 (n = 4.1, Hansen et al., 2012b) and Fo 90 (n = 3.3, Bystricky et al., 2000) are high. We therefore suggest that the additional strain accommodation processes are probably non-Newtonian.

Viscoplastic anisotropy under upper mantle conditions

In the present simulations, we do not model specific additional deformation processes, but analyze their effect on the viscoplastic mechanical behavior of the textured rock (polycrystal) by modeling their contribution to strain accommodation. In other words, we calculate how the operation of these processes reduces the geometrical hardening. However, to extrapolate the conclusions of our analysis to the upper mantle, assumptions must be made about the nature of these processes. Increased amounts of strain from diffusive mass transfer or sliding on grain boundaries relative to strain from dislocation glide are favored by smaller grain sizes. This grain-size dependence has been clearly imaged in situ by the difference in strain distribution in deformation experiments on halite aggregates with different grain sizes [START_REF] Bourcier | Multiscale experimental investigation of crystal plasticity and grain boundary sliding in synthetic halite using digital image correlation[END_REF]. Thus, these grain-size sensitive processes are potentially more active in experiments that use fine-grained olivine polycrystals (5-20μm) than during deformation of coarse-grained peridotites in the upper mantle. Indeed, strain rates measured during axial compression experiments on fine-grained (< 10 μm) olivine aggregates were up to 4.6 times higher than strain rates predicted by a lower bound model for the same textures and loading conditions, but no enhancement in strain rate relative to the models was observed for experiments on aggregates with coarser (> 100 μm) grain sizes [START_REF] Tielke | Observations of grain size sensitive power law creep of olivine aggregates over a large range of latticepreferred orientation strength[END_REF]. In the present study, grain-size sensitivity of the additional processes is suggested by the comparison of the evolution of the CPO intensity in simple shear predicted by the modified SO-VPSC model and the CPO intensity measured experimentally (Fig. 8b): the observations on the finest-grained aggregates are better fitted by models with lower β values, in which the dummy slip systems representing these additional processes are more active. One should note that the present results do not imply that grain-size sensitive processes are isotropicthese processes might be anisotropic in the presence of a shape-preferred orientation [START_REF] Wheeler | Anisotropic rheology during grain boundary diffusion creep and its relation to grain rotation, grain boundary sliding and superplasticity[END_REF]. The present work solely indicates that the additional processes operating in the experiments allow avoiding strain compatibility problems arising from the viscoplastic anisotropy of the olivine crystal by accommodating strain components that cannot be achieved by the activation of the easy slip systems in olivine.

The partitioning of deformation among different processes (dislocation glide, diffusive mass transfer, grain-boundary sliding, and others) depends on temperature, stress, and chemical environment, but all other conditions being similar, an increase in grain size should favor dislocation glide. Olivine grain sizes in naturally deformed upper mantle rocks are typically plurimillimetric.

Finer-grained peridotites are limited to shear zones (< 1 km wide) localized in the shallow parts of the lithospheric mantle and, even within these shear zones, in the absence of synkinematic reactions, olivine grain sizes are most often > 100 μm (cf. [START_REF] Tommasi | Heterogeneity and anisotropy in the lithospheric mantle[END_REF][START_REF] Vauchez | Faults (shear zones) in the Earth's mantle[END_REF]. Thus, although strain rates and stresses are lower than in the experiments, one may expect a greater contribution of dislocation glide to plastic deformation in most of the upper mantle. If that is the case, our results suggest that the viscoplastic anisotropy of the upper mantle might be higher than the one measured in the experiments and better represented by the VPSC simulations in which the contribution of the dummy slip systems is reduced (β > 10).

By consequence, previous studies that coupled the standard VPSC model to finite element simulations, such as [START_REF] Knoll | A multiscale approach to model the anisotropic deformation of lithospheric plates[END_REF], [START_REF] Tommasi | Structural reactivation in plate tectonics controlled by olivine crystal anisotropy[END_REF][START_REF] Blackman | Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy[END_REF], have estimated the maximum effect of olivine viscoplastic anisotropy on the reactivation of preexisting plate-scale structures for creating new plate boundaries (Tommasi et al., 2009) and on the flow and seismic anisotropy patterns beneath oceanic plates [START_REF] Blackman | Effects of crystal preferred orientation on upper-mantle flow near plate boundaries: rheologic feedbacks and seismic anisotropy[END_REF]. On the other hand, the ratios of 10-100 between normal and shear viscosity used in convection models by [START_REF] Christensen | Some geodynamical effects of anisotropic viscosity[END_REF], [START_REF] Lev | Rayleigh-Taylor instabilities with anisotropic lithospheric viscosity[END_REF], 2011[START_REF] Becker | On the role of anisotropic viscosity for plate-scale flow[END_REF] moderately to strongly overestimate viscous anisotropy associated with strain-induced CPO in the upper mantle. Indeed, the CPO-induced viscoplastic anisotropy in olivine-rich rocks predicted by standard SO-VPSC models does not surpass a factor 6.

Conclusion

The viscoplastic anisotropy of olivine polycrystals deformed in simple shear and axial extension was modeled using three polycrystal plasticity models: the uniform stress (lower bound) approach and the tangent and secondorder VPSC models, which was modified to account for variable contribution of other deformation processes in addition to dislocation glide. The instantaneous mechanical responses predicted in the various simulations were compared to stresses measured in laboratory experiments for similar textures, loading geometries, and boundary conditions. Laboratory data and numerical experiments suggest that the intensity of the viscoplastic anisotropy depends non-linearly on the CPO intensity; it increases faster for weak to moderate CPO intensities than at high CPO strengths. VPSC simulations in which strain is essentially accommodated by dislocation glide (β = 50) slightly underestimate the strength of the strongly textured aggregates relative to the strengths measured experimentally in torsion, but strongly overestimates the strength relatively to the experimental data in axial extension. This leads to overprediction of the viscoplastic anisotropy. Lower bound simulations, in which strain compatibility constraints are relaxed, better fit the experimental data. This suggest that in the laboratory experiments additional deformation processes operate, allowing the accommodation of the strain components that cannot be achieved by the activation of the easy slip systems in olivine.

The experimental data for the strongly textured aggregates is indeed well reproduced by the SO-VPSC model that includes dummy slip systems simulating additional isotropic strain accommodation processes without producing grain rotation. The best fit to the experimentally measured viscoplastic anisotropy is obtained for simulations in which the CRSS of the dummy systems is 3 times higher than the CRSS of the easiest slip system in olivine, (010) [100]. In these models, the dummy slip systems accommodate up to 50% of the viscoplastic deformation in axial extension and up to 20% of the deformation in simple shear. This modified SO-VPSC model also better reproduces the texture evolution measured in torsion experiments than the standard VPSC model. However, the best fit is obtained for simulations with lower activity of additional isotropic processes (β = 5); simulations with β = 3 underestimate the CPO strength relatively to the experimental data and do not evolve into a stable orientation at high shear strains.

The present work implies that the viscoplastic anisotropy associated with deformation by dislocation creep in the upper mantle may be decreased by the operation of additional strain accommodation processes. If we assume that in most of the upper mantle the contribution of deformation processes such as grain boundary sliding or diffusive mass transfer is lower than in the experiments due to larger grain sizes, the viscoplastic anisotropy should be well predicted by the standard SO-VPSC model. In high-stress shear zones in the shallow part of the lithospheric mantle, in which grain sizes are small, the viscoplastic anisotropy is probably reduced and better predicted by the modified SO-VPSC model. The standard (β = 50) and modified (β = 3-5) SO-VPSC models therefore reflect the possible range of texture-induced viscoplastic anisotropy in the upper mantle, that is, the viscous anisotropy relevant to domains of the upper mantle deforming dominantly by dislocation creep. The present simulations, which use CRSS for the olivine slip systems derived from high temperature low pressure single crystal experiments, provide good predictions of anisotropic viscous deformation for the upper 200-250 km of the mantle. Different CRSS derived from high-pressure experiments should be used for deeper parts of the upper mantle.

Fig. 1 .

 1 Fig. 1. (a) Stress-strain curves for torsion experiments by Hansen et al. (2012a). (b) CPO produced by the torsion experiments and used as initial textures in the extension experiments and in the polycrystal plasticity models. (c) Stress-strain curves for extension experiments byHansen et al. (2012a). Instantaneous stresses predicted by the models in simple shear and extension are compared to (a) final stresses in torsion and (c) peak stresses in extension, respectively. The stresses in extension are averaged values normalized to a constant strain-rate of 2 × 10 -4 s -1 . The number of grains in each CPO is indicated by the value n, and the CPO intensity by the value J (J-index). Pole figures are contoured as multiples of uniform distribution; the same contour range limits are applied to all pole figures for easy evaluation of variations in CPO strength.

Fig. 2 .

 2 Fig. 2. (a) Instantaneous normalized Von Mises equivalent stresses (relative to behavior of an isotropic polycrystal) and (b) magnitude of viscoplastic anisotropy (ratio between Von Mises equivalent stresses in extension and simple shear) for lower bound and VPSC simulations (colored symbols), in which the CRSS of the dummy slip systems (β) is set to 50 for the three initial olivine CPO (recalled in Fig. 2b and presented as a function of the CPO intensity -J-index). Black circles indicate stresses measured in the laboratory experiments. The empty markers present the results for simple shear parallel to the [100] maximum in a plane normal to the [010] maximum of the initial CPO. The filled markers present the results for axial extension parallel to the [010] maximum of the initial CPO. For calculating the magnitude of anisotropy the stresses in extension were averaged; for the color legend and labels of the pole figures, as well as loading geometries, refer to Fig. 1.

(

  010)[100], (001)[100] and (010)[001]. By consequence, high stresses are required to deform the polycrystal by dislocation creep.

Fig. 3 .

 3 Fig. 3. Instantaneous Von Mises equivalent stresses (normalized to the strength of the easy (010)[100] slip system) predicted by different polycrystal plasticity models as a function of the initial CPO intensity for simple shear (empty markers) and axial extension (filled markers). The strength of the dummy slip systems (β) is set to 50. Pole figures show the CPO measured at the end of the torsion experiments (Hansen et al., 2012a), which are used as the initial CPO in the models; for the color legend and labels of the pole figures, as well as loading geometries, refer to Fig. 1.

Fig. 4 .

 4 Fig. 4. Slip system activity predicted by the SO-VPSC model for simple shear (empty markers) and axial extension (filled markers). The sum of the activity of the three dummy systems, that is, the total contribution of the strain accommodation processes, is shown as black asterisks. Strength of the dummy slip systems (β) is set to 50 in these simulations. The activity of (100)[001], (011) [100], and (110)[100] slips systems, which accommodate less than 1% of viscoplastic deformation, is not represented. Pole figures show the CPO measured at the end of the torsion experiments (Hansen et al., 2012a), for which the instantaneous viscoplastic response are calculated; for the color legend and labels of the pole figures, as well as loading geometries, refer to Fig. 1.

Fig. 5 .

 5 Fig. 5. Effect of the critical resolved shear stress (CRSS) of the dummy slip systems (β) on the Von Mises equivalent stress (normalized to the strength of the easy (010)[100] slip system) predicted by the SO-VPSC model for simple

Fig. 6 .

 6 Fig. 6. a) Instantaneous Von Mises equivalent stresses normalized by the behavior of an isotropic polycrystal and (b) magnitude of the viscoplastic anisotropy predicted by SO-VPSC simulations with CRSS of the dummy systems (β) between 50 and 1 for the three experimental olivine CPO as a function of increasing CPO intensity (J-index). Empty markers show results for simple shear parallel to the [100] maximum in a plane normal to the [010] maximum of the preexisting CPO and filled markers indicate the results for axial extension parallel to the [010] maximum of the preexisting CPO. The experimental data is represented by the black circles; for the color legend of the pole figures and loading geometries, refer to Fig. 1.

Fig. 7 .

 7 Fig. 7. Effect of the dummy slip system strength (β) on the activity of the slip systems predicted by SO-VPSC models for the polycrystals: (a) untextured, (b) PT0655, (c) PT0640, and (d) PT0651. Empty markers represent simple shear simulations and filled markers axial extension simulations. The activity of the (100)[001], (011)[100], and (110)[100] slips systems, which accommodate less than 1% of viscoplastic deformation, is not shown. The sum of the activity of the three dummy systems, that is, the total contribution of the strain accommodation processes, is shown as black asterisks; for the color legend of the pole figures and loading geometries, refer to Fig. 1.

Table 1 :

 1 Slip systems parameters used in the VPSC simulations*

	Slip Systems	Critical Resolved Shear Stress#	Stress exponent
	(010)[100]	1	3
	(001)[100]	1.5	3
	(010)[001]	2	3
	(100)[001]	3	3
	(011)[100]	4	3
	(110)[001]	6	3
	{111}<110> a		3
	{111}<011> a		3
	{111}<101> a		3
	* Lower bound simulations only consider the first four slip systems
	# Adimensional values; normalized by the flow stress of the (010)[100] slip system
	a dummy slip systems, which represent additional strain accommodation processes and do not
	contribute to plastic spin

Table 2 :

 2 Model predictions and laboratory measurements

Instantaneous stresses predicted by the models and final (simple shear) and peak (extension) stresses measured in the experiments Instantaneous stresses* predicted by SO-VPSC models for different ß values Mechanical Solicitation

		Initial texture	SO-VPSC* ß=50	TGT-VPSC* ß=50	Lower Bound*	Laboratory (MPa) Data #	ß=10	5	3	2	1
		UNTEXTURED	16.77	9.64	3.94	-	8.33	6.03 4.56 3.53 2.10
	Simple	PT0655	14.86	8.57	3.59	148.23	7.54	5.56 4.29 3.39 2.09
	Shear	PT0640	8.43	5.87	2.86	118.47	4.82	3.82 3.18 2.73 1.95
		PT0651	7.58	5.19	2.57	109.43	4.38	3.48 2.92 2.50 1.78
		UNTEXTURED	16.72	10.05	3.89	168.00	8.27	5.97 4.52 3.50 2.08
	Extension	PT0655 PT0640	18.00 28.04	11.20 20.82	3.94 6.00	201.24 199.44	8.64 12.31 8.03 5.53 3.98 2.14 6.12 4.56 3.50 2.06
		PT0651	37.13	33.32	5.69	209.85	14.88 9.08 6.00 4.22 2.23

* Von Mises equivalent stress normalized by the strength of the (010)[100] slip system # Von Mises equivalent stress recalculated for a homogeneous strain rate of 2×10 -4 s -1 ; extension values are average ones.

Acknowledgments

We thank Ricardo Lebensohn and Carlos Tomé for making the code VPSC7c, which includes most of the last developments in viscoplastic self-consistent modeling, freely available. We also gratefully acknowledge Alain Vauchez, Sylvie Demouchy, Catherine Thoraval for discussions and Manuel Thieme and Sonia Ouadahi for comments on the figures. Anonymous reviewers are thanked for insightful comments. This project received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 642029 (ITN-CREEP) and from the CNRS (France) and CONICET (Argentine) under the Project International de Cooperation Scienti que MicroTex (No. 067785).


Augustin Charvet
email: augustin.charvet@univ-lorraine.fr

Loïc Wingert

Nathalie Bardin-Monnier

Stéphanie Pacault

Frantz Fournier

Denis Bemer

Dominique Thomas

Multi-staged granular beds applied to the filtration of ultrafine particles: An optimization of collector diameters

Keywords: Granular bed, Optimization, Performances, Nanoparticles

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A large amount of ultrafine particles are generated by manufacturing processes such as surface processing (thermal spraying), machining (drilling, grinding…), combustion, welding and cutting of metals, ... The development of suitable protection systems therefore appears crucial in order to protect the persons and the environment. The most efficient and widely used de-dusting methods for separating the particles from a carrier fluid are currently pleated fiber filters. These fibrous media are initially very efficient but generate a strong increase in pressure drop because of a quick and sometimes irreversible clogging as well as an extra cost resulting from their replacement. Following these observations, we propose to investigate the potential of other methods that could be applied to ultrafine particle collection. Among them, granular beds exhibit interesting performances, as far as ultrafine particle collection efficiency is concerned, and a higher resistance to mechanical and thermal stresses comparatively to fibrous filters. Moreover, another advantage of granular bed filters is their cleanability and reusability in opposition to the difficulty of unclogging the particles deposited on fibrous filters. This specificity is even more important because the reusability of the granules reduces the variability of pressure drop and collection efficiency in research investigations and of process performances in industrial applications.

Aerosol collection in granular beds has been the subject of many theoretical and experimental studies since 1970 and a thorough review of those can be found in Tien's book [START_REF] Tien | Granular Filtration of Aerosols and Hydrosols[END_REF]. Experimental studies focused mainly on the overall performances (pressure drop and collection efficiency) of a granular filter by changing the filter construction parameters, superficial gas velocity [START_REF] Tsubaki | Gas filtration in granular moving beds -an experimental study[END_REF], temperature [START_REF] Peukert | Influence of temperature on particle separation in granular bed filters[END_REF], bed depth and porosity, granular collector characteristics [START_REF] Saxena | Particulate removal from high-temperature, high-pressure combustion gases[END_REF], ... Results highlight that an optimization of the granular bed parameters (collector diameter, granular bed thickness, face velocity of the effluent, magnetic stabilization, etc.) could bring this technique to become a viable filtration method for ultrafine particles. However, all of these studies were limited to the characterization of granular beds in their initial state and the evolution of their performances during the particle loading, i.e. during the accumulation of particles within the bed, is much less examined. Moreover, these studies were dedicated to the performances of granular beds towards micronic and submicronic particles and consequently results do not guarantee the performances of such granular bed filters towards ultrafine particles, nanostructured or not. Furthermore, Bémer et al. [START_REF] Bémer | Experimental study of granular bed filtration of ultrafine particles emitted by a thermal spraying process[END_REF] showed that ultrafine particle deposit is essentially located in the first layers of a granular bed. Thus, only a small depth is useful thereby reducing both operation time of granular beds between two uncloggings and the process yield. In order to optimize the performance of granular beds, Thivel et al. [START_REF] Thivel | Magnetically stabilized fluidization of a mixture of magnetic and non-magnetic particles in a transverse magnetic field[END_REF] emphasized that it may be interesting to use a magnetically stabilized bed because it combines high dust collection efficiencies with a reduced energy cost. In more recent studies, special attention has been devoted to the regeneration of the filter granules. In this way, moving granular bed filters, in countercurrent [START_REF] Brown | Similitude study of a moving bed granular bed[END_REF][START_REF] El-Hedok | The influence of granular flow rate on the performance of a moving bed granular filter[END_REF] or cross-currents configurations [START_REF] Chen | Filtration of dust particulates with a moving granular bed filter[END_REF][START_REF] Chen | Removal of dust particles from fuel gas using a moving granular bed filter[END_REF][START_REF] Hsiau | Flow of filter granules in moving granular beds with louvers and sublouvers[END_REF] use a continuous flow of fresh collectors to limit the formation of a filtration cake and maintain a low and constant pressure drop, without the need for a periodic unclogging. However, moving granular beds have some drawbacks related to the transport of collectors and present a lower collection efficiency as dust cake is continuously re-entrained. Moreover, operating principle of panel beds is somewhat different and based on a periodic renewal of the granular bed surface. Hence, these panel beds mimic fabric filters by taking advantage of the high collection efficiency of a particle cake. To decrease the pressure drop, a cyclical offline puff-back induces the fall of this filtration cake and of the first layers of collectors which are replaced with fresh ones [START_REF] Lee | Granular-bed filtration assisted by filter-cake formation 1. Exploiting a new mode of soil failure for renewal of filtration surfaces in a panel bed[END_REF][START_REF] Rodon | Granular-bed filtration assisted by filter-cake formation 2. The panel bed gas filter with puffback renewal of gas-entry surfaces[END_REF][START_REF] Wu | Granular-bed filtration assisted by filter-cake formation 3. Penetration of filter cakes by a monodisperse aerosol[END_REF][START_REF] Squires | Granular-bed filtration assisted by filter-cake formation 4. Advanced designs for panel-bed filtration and gas treating[END_REF].

Our optimization approach is quite different and consists in increasing the time between two unclogging operations by reducing the pressure drop increase while maintaining a high collection efficiency. Tian et al. [START_REF] Tian | Cascade filtration properties of a dual-layer granular bed filter[END_REF] applied this approach, using a dual-layer granular bed for the filtration of micron-sized ash particles. Nevertheless, the very thick, and consequently very efficient, upper stage implemented in this study reduces the interest of the lower stage. Indeed, during the clogging, the deposit is mainly located in the first layers of collectors which render the remaining thickness of the upper stage unnecessary. With this in mind, a granular bed with a gradient of collector sizes was proposed. From a theoretical standpoint, as the collection efficiency is a decreasing function of the collector size in a diffusional regime, a granular bed with decreasing collector sizes in the depth should permit to firstly clog the lower stages of the granular bed and consequently delay the pressure drop increase (cf. Figure 1). To limit the experimental constraints, we chose to study the simplest "gradient" granular bed, consisting in only three stages. Each stage is composed of beads satisfying the condition that a given stage must be filled with smaller beads than those of the previous upper stage. Moreover, it was decided to limit the thickness of each stage given that only the first layers of collectors are effective in a conventional granular bed, i.e. a granular bed composed of beads with given diameter.

Experimental approach

Experiments were achieved using a three-stage 40 mm-diameter column filled with stainless steel beads. The performances of a multi-staged granular bed composed of 1 mm-diameter beads in the first stage, 0.8 mm in the second and 0.5 mm in the third were compared to those of a conventional one, i.e. a granular bed with bead diameter of 0.5 mm in the three stages. For the reasons mentioned above, the depth of each granular bed stage was chosen to approximately 11 mm. To ensure that the porosity was as close as possible to the one of a random close packing of spheres, the granular beds in each stage were packed by applying vibrations until reaching a relative gap lower than 5% between the experimental pressure drop and the one predicted by the Kozeny-Carman's law. Static differential pressure sensors were mounted upstream and downstream of each stage allowing the monitoring of the pressure drop evolution during the clogging.

The particles, composed of 85% of zinc and 15% of aluminium, were generated by thermal spraying of metal using an electric arc gun (Margarido® M25). These nanoparticles were then directed towards a sharp cut cyclone (BGI® SCC 1.829) to collect the micron-sized fraction of the produced aerosol and clog the granular bed only with the fine fraction. The particle size distribution of the resulting particles of Zn-Al was measured with a Scanning Mobility Particle Sizer (SMPS 5.400 Grimm®) and presented a median number diameter, expressed in equivalent electric mobility, of 78.3 nm with a geometric standard deviation of 1.6. The total number concentration was approximately 2.10 8 particles/cm 3 . The total mass concentration, measured by a gravimetric method after sampling on glass fiber filters, was close to 72.5 mg/m 3 . The air flowrate was controlled by a sonic nozzle and measured by a mass flowmeter (TSI® 4100). It was set to 15 L/min corresponding to a superficial velocity of 0.2 m/s. The particle number size distributions were measured upstream and downstream of each stage using a Nanoscan analyzer (SMPS Model 3910, TSI®) after a dilution stage (VKL100, Palas®). These sequential measurements, called cycles, consisted of three scans of one minute upstream of each stage followed by three scans downstream. The three scans of a same sampling cycle were averaged to obtain the mean upstream and downstream fractional number concentrations (Cn,i) which were converted according to equation 1 into fractional mass concentrations (Cm,i). It was therefore possible to obtain the fractional mass collection efficiencies and finally the total collected mass collected at each time step for each stage.

C ",$ = π • ρ ) 6 • C +,$ • d -,$ . = π • ρ /,$ 6 • C +,$ • d ",$ . (1)
ρ /,$ = 40.24 • d ",$ 45.67 [START_REF] Tsubaki | Gas filtration in granular moving beds -an experimental study[END_REF] where ρp is the density of the raw material (5740 kg/m 3 for the Zn-Al), dv,i the volume equivalent diameter, dm,i the electric mobility equivalent diameter and Cn,i the measured fractional number concentration and ρe,i the effective density of the Zn-Al agglomerates. The effective density, defined as the ratio of aggregate mass to equivalent mobility volume, allows converting the volume equivalent diameter into an electric mobility equivalent diameter and consequently calculating the fractional mass concentrations. It should be noted that in the equation 2, providing the relationship between the effective density of the agglomerates and the electric mobility equivalent diameter, the diameter is expressed in nm and the effective density in g/cm 3 [START_REF] Charvet | On the Importance of Density in ELPI Data Post-Treatment[END_REF]. Furthermore, this equation is only applicable for electric mobility equivalent diameters higher than the diameter of the primary particles constituting the agglomerates, i.e. 9 nm for Zn-Al primary particles ; otherwise the effective density will become higher than the raw particle material.

Figure 2: Total mass efficiency and pressure drop evolutions of conventional and multi-staged granular beds

The performances of the two configurations were compared in terms of global mass efficiency and pressure drop as function of the collected mass (cf. Figure 2). As collection efficiency is a decreasing function of the collector size, the results show a lower initial efficiency for the multi-staged granular bed compared to the conventional one. But, this gap is rapidly reduced during the clogging process. At the same time, a lower pressure drop increase of the multi-staged granular bed is shown by the results. Indeed, the pressure drop is quite three times lower than the one of the conventional bed when the total mass efficiency reaches a value close to unity. The plot of the mass collected according to time and stage confirms that the deposition in a conventional granular bed mainly occurs in the first stage (cf. Figure 3). In the case of the multi-staged granular bed, the mass collected seems to be more homogenously distributed. Indeed, this figure clearly shows that the use of a gradient of bead diameter permitted to firstly clog the downstream layers of the granular bed, i.e. the third stage of the multi-staged column. According to these experimental results, the improvement of the granular bed filtration by means of a collector diameter gradient appears quite promising in order to extend the operation time between two uncloggings. Nevertheless, it was noticed that the first stage of this configuration becomes efficient sooner than expected and consequently that the second and third stages did not contribute enough to the particle retention. This combination of bead diameters (1 mm / 0.8 mm / 0.5 mm) is consequently not optimal.

Modeling approach

As a result of these experimental observations, it was decided to find an optimized configuration of bead diameters so that the whole granular bed depth contributes to the filtration. As an experimental approach would be too timeconsuming, we shifted towards a modeling approach using a previously developed model predicting the temporal evolution of pressure drop and collection efficiency and determining the location of the mass collected in the granular bed depth, for given operating conditions (filtration velocity, pressure, temperature …) and given granular bed parameters (bead size, thickness …).

This model, developed for conventional granular beds, has been modified to evaluate the performances of multistaged granular beds. These adjustments only consist of a repetition of the original model for each stage of the granular bed placed in series. In this way, the output data of the first stage become the input data of the second stage and so on. A brief description of this clogging model is given below. For more detailed information, readers can refer to [START_REF] Wingert | Modeling of the deep granular bed clogging by nanoparticles[END_REF].

In the model, the deposit of ultrafine particles is assimilated to a uniform and spherical shell all over the collector's surface. This assumption has been used in previous studies [START_REF] Tien | Simulation of the dynamic behavior of deep bed filters[END_REF][START_REF] Fichman | Effect of particle loading on granular bed filtration -The cluster enhanced filter model[END_REF] and validated by microscopic observations [START_REF] Charvet | Visualization of airborne nanoparticle deposits onto spherical collectors[END_REF]. According to this uniformity feature of the deposit, the principle is based on the calculation of equivalent collector diameters depending on the collected mass of particles and to implement them in initial pressure drop (equation 3) and collection efficiency (equation 4) models.

ΔP = 36 • h < • µ • U • ? (1 -ϵ DE ) G ϵ DE . H • ? z DE d J,/K G H (3)
E DE = 1 -exp P- 3 2 • ? 1 -ϵ DE d J,/K H • z DE • ηRd J,/K ST (4)
with hk the Kozeny-Carman constant, µ the air viscosity, U the superficial velocity, zGB the granular bed depth, dc,eq the equivalent collector diameter, ϵGB the granular bed porosity and η(dc,eq) the single collector collection efficiency depending on the equivalent collector diameter.

Results of clogging experiments have highlighted that the deep bed filtration can be split in two phases requiring the use of two different kinds of collector equivalent diameter. During the first phase, called phase A, the pressure drop and the efficiency remain approximately constant indicating an insignificant change in the granular beds structure. The collector diameter is assumed to increase slightly without any blocking of the pore spaces. Thus, during the phase A, the equivalent collector diameter is considered as the one of the sphere having the same volume as a collector and its deposit layer and can be expressed as:

d J,/K,U = ?d J . + 6 • m X,J,U π • ρ X • (1 -ϵ Y ) H 7 .
(

) 5
with dc,eq,A the instantaneous equivalent collector diameter during the phase A, dc the initial collector diameter, mp,c,A the instantaneous mass collected per collector during the phase A, ρp the density of the raw material and εd the deposit's porosity.

The second clogging phase, called phase B, is characterized by a sharp increase of the pressure drop and of the collection efficiency. This performance evolution is assigned to an increase of the collector specific area because of the deposit accumulation. Consequently, during the phase B, the equivalent collector diameter is defined as the one of the sphere having the same specific area as a collector and its deposit layer:

d J,/K,E = ? π • d J,/K,U . • ρ X • (1 -ϵ Y ) • d Z,[5 + 6 • d Z,[5 • m X,J,E π • d J,/K,U G • ρ X • (1 -ϵ Y ) • d Z,[5 + 4 • (1 -ϵ Y ) • m X,J,E H (6)
with dv,50 the median volume equivalent diameter of the collected particles and mp,c,B the instantaneous mass collected per collector during the phase B.

As the equivalent diameters are calculated at each increment of collected mass, this model allows determining the temporal performances of the granular bed. Concerning the spatial arrangement of the deposit in the granular bed depth, it is made possible by means of its discretization in layers.

Multi-criteria optimization

All the operating conditions previously described are maintained (filtration velocity, thickness of each stage …) and the model was implemented for a set of possible combinations of collector diameters. The collector size ranged from 0.2 mm to 1.6 mm by increment of 0.1 mm. The multi-staged granular beds simulated had to satisfy the condition that a given stage must be filled with smaller beads than those of the previous upper stage.

These settings led to 455 combinations of collector diameters. The performances of these 455 multi-staged granular beds were compared with each other based on three criteria: § The retention capacity, i.e. the total mass collected per volume unit during the time necessary to reach an efficiency of the first stage close to 99%; it is considered that above this value, the particles will no more be collected in the granular bed depth but at the surface of the first stage. In this way, industrially, the beginning of the surface filtration will correspond to a need of an unclogging operation. § The mean mass efficiency during this same clogging duration, i.e. the ratio of the total collected mass to the total generated mass; § The reciprocal of the pressure drop.

To obtain the best performances, each of these criteria should be maximized. As the optimum of each criterion is obtained for different combinations and as each criterion has not the same importance, it is difficult to identify the best trade-off solution. A classical approach consists in merging all criteria in a single one. This aggregation approach induces some subjectivity as it is necessary to attribute a relative weighting to each criterion. We chose another approach in order to reduce this subjectivity. It is based on an optimization in a first step and introducing preferences, i.e. subjectivity, as late as possible in the decision process in order to identify the whole set of trade-off solutions. To identify this set of trade-off solutions, the Pareto's domination concept has been applied [START_REF] Pareto | Manuale di Economia Politica[END_REF]. This analysis relies on a comparison of all the combinations by pairs in order to determine if a given combination is dominated by another one.

Considering two combinations X and Y, X is regarded as not dominated by Y if: § All the criteria of X are higher than or equal to all the equivalent criteria of Y; § At least one criterion of X is strictly better than the equivalent criterion of Y.

During this analysis by pairs, a range of indifference has been added in order to take into account the model uncertainty. This range is defined as the difference between two values of the same criterion that is not significant enough to consider a value as preferred over another one.

By selecting only the never dominated configurations, this Pareto's domination analysis allows considerably reducing the number of possible combinations (61/455), without introducing subjective weightings between the three chosen criteria (cf. Figure 4). To go further, subjective weightings have to be introduced to define a ranking between these 61 selected combinations. A multi-criteria optimization technique, called rough set method, has been applied [START_REF] Slowinski | Intelligent decision support: Handbook of Applications and Advances of the Rough Sets Theory[END_REF][START_REF] Thibault | Multicriteria optimization of a high yield pulping process with rough sets[END_REF]. This method consists in the selection of a small sample of seven combinations, being broadly representative of those retained after the multi-criteria analysis described previously. These seven chosen combinations have to be manually ranked, from the most to the least preferred, by a human expert who has a deep knowledge of the process.

Then, the combinations of the sample are compared by pairs for each criterion in order to establish a set of rules based on the expert's ranking. If one criterion of a preferred combination is better than the same criterion of a less preferred one, the value of 1 is affected to this criterion. Otherwise, it obtains the value of 0. From these comparisons result binary rules of preference. For example, the preference rule 1/0/1 is obtained if the first criterion of a combination is better, the second criterion worse and the third one better than the same criteria of a less preferred combination. Repeating this step, from the less to the most preferred combination, permits to get the binary rules of non-preference. The binary rules appearing as both a preference and a non-preference rule are eliminated.

The remaining ones are then applied to the 61 previously selected combinations by comparing them by pairs. If a comparison of two combinations corresponds to a rule of preference, the ranking score of the first combination is incremented by one unit and the score of the second combination decremented by one unit. If the comparison leads to a non-preference rule, one counts -1 point for the first combination and 1 for the second one. When the comparison matches neither a preference nor a non-preference rule, the ranking scores of the two combinations remain unaltered. This allows ranking all the combinations according to their final score, the highest representing the one that has to be finally chosen.

After the two optimization phases, consisting in the Pareto's domination analysis and the Rough set method, the optimal trade-off solution obtained is a combination with small bead diameters (0.6/0.4/0.2 mm). The results are quite similar whatever the expert chosen. The modeled performances of this optimal multi-staged granular bed have been compared to the ones of the two combinations experimentally tested and to the ones of the extremum combination, i.e. constituted of the smaller collectors.

Figure 5: Efficiency and pressure drop evolutions during the clogging of multi-staged granular beds with various combinations of collector diameters

Figure 5 (left) represents the efficiencies modeled for each combination as a function of the collected mass. As the collection efficiency is a decreasing function of the collector diameter, the modeled efficiencies are all the more high that the combination is constituted of small beads. Furthermore, for the combinations 1.0/0.8/0.5 mm and 0.5/0.5/0.5 mm, we can notice a good agreement between experimental data and model-based predictions even if there is a slight underestimation of the collection efficiency by the model. This difference is mainly caused by a lack of precision in the correlation used in the model for the calculation of initial efficiency by Brownian diffusion [START_REF] Wingert | Modeling of the deep granular bed clogging by nanoparticles[END_REF]. As the Neale and Nader hydrodynamic factor implemented in the model tends to slightly underestimate the efficiency of a conventional granular bed, this trend is even more marked for three granular beds in series.

At the same time, the model satisfactorily predicts the experimental pressure drop evolution (Figure 5, right). The combination 1.0/0.8/0.5 mm presents the lower energy expenditure among the tested multi-staged granular beds but also the lower mass efficiency. The trade-off solution appears to be the combination 0.6/0.4/0.2 mm as the mass efficiency is higher than 80% and the pressure drop increase relatively slow. The model also allows a fairly good prediction of the clogging profile as modeled and experimental mass fractions per stage follow the same trend (cf. Figure 6). As previously mentioned, the experimental fractions of mass collected are calculated from the upstream and downstream fractional mass concentrations deduced from the Nanoscan data. Figures 6A and6D highlight that a gradient of collector diameters allows clogging firstly the deepest stage contrary to a conventional granular bed in which the particle deposition mainly occurs in the upper stage (cf. Figure 6B). After a long period of particle loading, the particles are more homogeneously distributed in the depth of the bed. Contrary to a conventional granular bed in which particles are collected in the first layers, the particle mass is retained in a higher volume and consequently leads to a lower pressure drop increase. The results obtained for the combination 0.4/0.3/0.2 mm show that if collector diameters in two consecutive stages are too close to each other, the multistaged granular bed loses all of its advantages as particle deposition quickly becomes mainly located in this first stage rending the two other stages unnecessary, as in a conventional granular bed (Figure 6C). This is all the more emphasized since the upper stage is efficient. Thus, the difference between each initial collection efficiency should be sufficient to preferentially collect the particles into the lower stage.

Conclusion

According to experimental results, using a multi-staged granular bed with a collector diameter gradient appeared promising in order to obtain a more homogeneous deposit and extend the operation time between two uncloggings. The choice of collector diameters is essential to reach simultaneously a high initial efficiency and a more homogeneous deposit, i.e. a delayed pressure drop increase. A predictive clogging model was also developed in order to simulate the temporal performance evolution of a set of possible combinations of collector diameter. Thereafter, a multi-criteria optimization was led and permitted to find an optimized configuration of bead diameters (0.6/0.4/0.2 mm) so that the whole granular bed depth contributes to the filtration.

Figure 1 :

 1 Figure 1: Schematic view of the multi-staged granular bed

Figure 3 :

 3 Figure 3: Mass of particle collected in each of the three stages during the clogging of a multi-staged (left) and a conventional (right) granular bed

Figure 4 :

 4 Figure 4: Result of the Pareto domination analysis (the points represent the 61 trade-off solutions among the 455 possible combinations within the tetrahedron)

Figure 6 :

 6 Figure 6: Mass fraction collected in each stage during the clogging of granular beds with various combinations of collector diameters. Solid, broken and dotted lines correspond to the modelling results and diamonds, squares and circles to experimental results for the first, second and third stages, respectively