Neil J Goulding

Neil M Ribe

Olivier Castelnau

Andrew M Walker

James Wookey

Analytical parametrization of self-consistent polycrystal mechanics: Fast calculation of upper mantle anisotropy

Keywords: Mantle processes, Creep and deformation, Seismic anisotropy

Progressive deformation of upper mantle rocks via dislocation creep causes their constituent crystals to take on a non-random orientation distribution (crystallographic preferred orientation or CPO) whose observable signatures include shear-wave splitting and azimuthal dependence of surface wave speeds. Comparison of these signatures with mantle flow models thus allows mantle dynamics to be unraveled on global and regional scales. However, existing self-consistent models of CPO evolution are computationally expensive when used with 3-D and/or time-dependent convection models. Here we propose a new method, called ANPAR, which is based on an analytical parametrization of the crystallographic spin predicted by the second-order (SO) self-consistent theory. Our parametrization runs ≈2-6 × 10 4 times faster than the SO model and fits its predictions for CPO and crystallographic spin with a variance reduction >99 per cent. We illustrate the ANPAR model predictions for the deformation of olivine with three dominant slip systems, (010)[100], (001)[100] and (010)[001], for three uniform deformations (uniaxial compression, pure shear and simple shear) and for a corner-flow model of a spreading mid-ocean ridge.

I N T RO D U C T I O N

Seismic anisotropy observed in Earth's upper mantle is typically explained by the partial alignment of the lattices of the constituent olivine and pyroxene crystals caused by deformation associated with mantle convection (e.g. [START_REF] Nicolas | Formation of anisotropy in upper mantle peridotites: a review, in composition, structure and dynamics of the lithosphere-asthenosphere system[END_REF][START_REF] Silver | Seismic anisotropy beneath the continents: probing the depths of geology[END_REF][START_REF] Long | Mantle dynamics and seismic anisotropy[END_REF]. Because each crystal is elastically anisotropic, this non-random distribution of crystallographic directions (called a crystallographic preferred orientation, or CPO) will impart elastic anisotropy to the bulk material. The seismically observable consequences of this anisotropy include shearwave birefringence or 'splitting' (e.g. [START_REF] Crampin | An introduction to wave-propagation in anisotropic media[END_REF][START_REF] Silver | Shear-wave splitting and subcontinental mantle deformation[END_REF] and the azimuthal dependence of surface-wave speeds (e.g. [START_REF] Montagner | Global upper mantle tomography of seismic velocities and anisotropies[END_REF]. Simulation of the development of CPO in models of mantle deformation, and comparison of this with seismic observations of the Earth, allow mantle dynamics to be unraveled on global (e.g. [START_REF] Becker | On the relationship between azimuthal anisotropy from shear wave splitting and surface wave tomography[END_REF]) and regional scales (e.g. [START_REF] Long | Constraints on subduction geodynamics from seismic anisotropy[END_REF]). However, these simulations are computationally challenging when performed for time-dependent models of mantle convection or at the high spatial resolution needed for finite frequency simulation of seismic wave propagation. Here we describe an accurate and computationally efficient alternative to existing methods for the simulation of CPO development in the upper mantle.

The principal cause of CPO and seismic anisotropy in the mantle is the progressive deformation experienced by mantle rocks as they participate in the global convective circulation. Under appropriate conditions of stress, temperature, and grain size, olivine and pyroxene crystals deform via dislocation creep, whereby internal dislocations move through the crystal to accommodate strain. The dislocations move on crystallographic planes and in directions set by the crystal structure, and the combination of a plane and direction define the limited number of slip systems available to allow the crystal to deform. Deformation of this type constrains the crystallographic axes to rotate relative to a fixed external reference frame, much as a tilted row of books on a shelf rotates when one pushes down on it. Because crystals with different orientations rotate at different rates, the overall distribution of orientations evolves with time in a way that reflects both the geometry of the slip systems and the character of the imposed deformation.

Because CPO and seismic anisotropy are so directly linked to progressive deformation, observations of seismic anisotropy have the potential to constrain the pattern of convective flow in the mantle. Realizing this potential, however, requires a reliable polycrystal mechanics model that can predict how the individual crystals in an aggregate deform and rotate in response to an imposed macroscopic stress or strain rate. Three broad classes of polycrystal models have been proposed to date.

The first class comprises the full-field models. In these, the polycrystal is treated explicitly as a spatially extended body, and the stress and strain within it are field variables that vary continuously as a function of position. Full-field models allow the stress and strain to vary both among and within individual grains in a physically realistic way. This approach can be implemented as a finite element problem (e.g. [START_REF] Sarma | Effects of interactions among crystals on the inhomogeneous deformation of polycrystals[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] or, more efficiently, using a method based on fast Fourier transforms [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Lebensohn | N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform[END_REF][START_REF] Suquet | Multiscale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]. Predictions from full-field models agree remarkably well with laboratory experiments [START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF] and analytical results available for simple cases [START_REF] Lebensohn | Fast Fourier Transformbased modelling for the determination of micromechanical fields in polycrystals[END_REF]). However, their great computational expense makes them too slow (by many orders of magnitude) for routine use in convection calculations.

This disadvantage is overcome to some extent by so-called 'homogenization' models, in which the detailed spatial distribution of the grains is ignored and the aggregate is treated as a finite number of grains with different orientations and material properties. In this mean-field approach compatibility of stress and strain equilibrium is not enforced between spatially contiguous grains, but rather between each grain and a 'homogeneous effective medium' defined by the average of all the other grains. For viscoplastic behaviours as considered here, a well-known member of this class makes use of the so-called 'tangent' anisotropic scheme of [START_REF] Molinari | A selfconsistent approach of the large deformation polycrystal viscoplasticity[END_REF] and Lebensohn & Tome (1993). In this model the local stress and strain rate tensors vary among the grains. In the geophysical literature this approach is generally known as the viscoplastic self-consistent (VPSC) model, and we use this name for this firstorder approximation. The VPSC model has been widely used in solid-earth geophysics including studies of CPO development in the upper mantle (e.g. [START_REF] Wenk | Modelling plastic deformation of peridotite with the self-consistent theory[END_REF][START_REF] Tommasi | Upper mantle tectonics: threedimensional 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][START_REF] Tommasi | Structural reactivation in plate tectonics controlled by olivine crystal anisotropy[END_REF][START_REF] Wenk | Modeling dynamic recrystallization of olivine aggregates deformed in simple shear[END_REF][START_REF] Mainprice | Pressure sensitivity of olivine slip systems and seismic anisotropy of Earths upper mantle[END_REF][START_REF] Bonnin | Numerical modelling of the upper-mantle anisotropy beneath a migrating strike-slip plate boundary: the San Andreas Fault system[END_REF][START_REF] Leo | Development of texture and seismic anisotropy during the onset of subduction[END_REF], the transition zone [START_REF] Tommasi | Strain-induced seismic anisotropy of wadsleyite polycrystals and flow patterns in the mantle transition zone[END_REF], in the lowermost mantle [START_REF] Wenk | Modeling lower mantle anisotropy development in a subducting slab[END_REF][START_REF] Merkel | Deformation of (Mg,Fe)SiO 3 postperovskite and D anisotropy[END_REF][START_REF] Wenk | Deformation in the lowermost mantle: from polycrystal plasticity to seismic anisotropy[END_REF][START_REF] Mainprice | Predicted glide systems and crystal preferred orientations of polycrystalline silicate Mg-Perovskite at high pressure: implications for the seismic anisotropy in the lower mantle[END_REF][START_REF] Walker | Elastic anisotropy of D predicted from global models of mantle flow[END_REF][START_REF] Dobson | Strong inheritance of texture between perovskite and post-perovskite in the D layer[END_REF][START_REF] Nowacki | Evaluating post-perovskite as a cause of D anisotropy in regions of palaeosubduction[END_REF]Amman et al. 2014;[START_REF] Cottaar | Synthetic seismic anisotropy models within a slab impinging on the core-mantle boundary[END_REF], in the inner core (e.g. [START_REF] Wenk | A convection model to explain anisotropy of the inner core[END_REF][START_REF] Deguen | Texturing in Earth's inner core due to preferential growth in its equatorial belt[END_REF][START_REF] Lincot | Is inner core seismic anisotropy a marker for plastic flow of cubic iron?[END_REF], or in ice [START_REF] Castelnau | Viscoplastic modelling of texture development in polycrystalline ice with a self-consistent approach: comparison with bound estimates[END_REF][START_REF] Castelnau | Modelling viscoplastic behavior of anisotropic polycrystalline ice with a selfconsistent approach[END_REF]. However, as noted by Masson et al. (2000), there is an inconsistency in the common first-order VPSC approach in the definition of the stress localization tensor, leading to an inaccurate estimation of the effective rheology for highly anisotropic viscoplastic polycrystals such as olivine [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF].

More recently, an improved 'second order' (SO) self-consistent homogenization scheme has been proposed by Ponte Castañeda (2002). In the SO model the stress and strain rate varies among grains with the same orientation and physical properties, and these fluctuations are used to derive the effective polycrystal behaviour. As a result, its predictions of quantities such as the effective average stress in the aggregate are much more accurate than those of simpler homogenization schemes [START_REF] Lebensohn | Selfconsistent modelling of the mechanical behaviour of viscoplastic polycrystals incorporating intragranular field fluctuations[END_REF][START_REF] Castelnau | Micromechanical modeling of the viscoplastic behavior of olivine[END_REF]. Recent examples of the application of the SO approach to olivine deformation are provided by [START_REF] Castelnau | Micromechanical modeling of the viscoplastic behavior of olivine[END_REF][START_REF] Castelnau | Numerical simulations of texture development and associated rheological anisotropy in regions of complex mantle flow[END_REF][START_REF] Castelnau | Microstructures and rheology of the Earth's upper mantle inferred from a multiscale approach[END_REF] and [START_REF] Raterron | Multiscale modeling of upper mantle plasticity: from single-crystal rheology to multiphase aggregate deformation[END_REF].

While the physical self-consistency of the SO and simpler VPSC models is appealing, both are computationally expensive when ap-plied to typical mantle minerals deforming by dislocation creep. The reason is the strongly nonlinear rheology of such minerals, which makes it necessary to use iterative methods to solve the equations of stress compatibility among the large number (∼10 3 -10 4 ) of grains required to represent the polycrystal. Moreover, the number of iterations required at each deformation step increases rapidly as the CPO becomes progressively more strongly anisotropic. These difficulties render the VPSC and SO models unsuitable for calculations of evolving CPO in complex time-dependent mantle flow fields, unless powerful computer capacity and elaborate computation strategies are used. Indeed, because of these computational constraints, none of the studies referenced above make use of the VPSC or SO approaches to directly compute the elasticity on a fine spatial scale (suitable for finite frequency forward modelling of the seismic wave field) from a time-varying description of mantle flow. Instead various approximations are used, such as limiting the calculation to selected ray-theoretical paths [START_REF] Blackman | Seismic anisotropy of the upper mantle: 1. Factors that affect mineral texture and effective elastic properties[END_REF][START_REF] Nowacki | Evaluating post-perovskite as a cause of D anisotropy in regions of palaeosubduction[END_REF][START_REF] Leo | Development of texture and seismic anisotropy during the onset of subduction[END_REF], interpolating the calculated elasticity [START_REF] Bonnin | Numerical modelling of the upper-mantle anisotropy beneath a migrating strike-slip plate boundary: the San Andreas Fault system[END_REF], or simplifying the model of mantle flow [START_REF] Raterron | Multiscale modeling of upper mantle plasticity: from single-crystal rheology to multiphase aggregate deformation[END_REF].

A final degree of physical simplicity and computational efficiency is reached in models of the 'kinematic' class, which are based on either an analytical expression for the deformation-induced rate of crystallographic rotation [START_REF] Ribe | A theory for plastic deformation and textural evolution of olivine polycrystals[END_REF][START_REF] Kaminski | A kinematic model for recrystallization and texture development in olivine polycrystals[END_REF][START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF] or on a simple relationship between finite strain and the expected CPO [START_REF] Muhlhaus | Emergent anisotropy and flow alignment in viscous rock[END_REF][START_REF] Lev | Prediction of anisotropy from flow models: A comparison of three methods[END_REF]. One example, the DRex model [START_REF] Kaminski | A kinematic model for recrystallization and texture development in olivine polycrystals[END_REF][START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF]) has been widely used to predict CPO and seismic anisotropy from flow models (e.g. [START_REF] Lassak | Seismic characterization of mantle flow in subduction systems: Can we resolve a hydrated mantle wedge?[END_REF][START_REF] Conder | Rapid mantle flow beneath the Tonga volcanic arc[END_REF][START_REF] Becker | Radial seismic anisotropy as a constraint for upper mantle rheology[END_REF][START_REF] Long | Mantle dynamics and seismic anisotropy[END_REF][START_REF] Faccenda | Development of mantle seismic anisotropy during subduction-induced 3-D flow[END_REF][START_REF] Long | Constraints on subduction geodynamics from seismic anisotropy[END_REF][START_REF] Faccenda | Mid mantle seismic anisotropy around subduction zones[END_REF]. Kinematic models are computationally 10-100 times faster than homogenization models, and predict very similar CPO. However, the physical principle underlying the expression for the spin is ad hoc, and has not yet been adequately justified. Moreover, because the kinematic approach does not account explicitly for stress compatibility among grains, it cannot be used to predict rheological properties of a deforming aggregate.

In view of the above limitations, it would clearly be desirable to have a polycrystal model that combines the physical rigour of the self-consistent approach with a much lower computational cost. The aim of this paper is to derive such a model. For purposes of illustration, we consider the case of a pure olivine polycrystal (dunite), a relevant (albeit simplified) representation of the mineralogy of the upper ≈400 km of Earth's mantle. Our approach is to examine in detail the predictions of the SO model for dunites subject to different kinds of deformation, and to extract from those predictions a simple parametrization that can be expressed analytically. The active slip systems of mantle olivine are believed to vary with temperature, pressure and hydrogen content (e.g. [START_REF] Jung | Water-induced fabric transitions in olivine[END_REF][START_REF] Cordier | Dislocations and slip systems of mantle minerals[END_REF][START_REF] Mainprice | Pressure sensitivity of olivine slip systems and seismic anisotropy of Earths upper mantle[END_REF]. In order to reduce the parameter space that we need to consider, we limit ourselves to the three dominant slip systems relevant to deformation under anhydrous upper mantle conditions: (010)[100], (001)[100] and (010) [001]. However, we allow the relative importance of these three slip systems to vary and, as discussed below, these slip systems include the most important ones under mantle conditions down to 410 km [START_REF] Castelnau | Microstructures and rheology of the Earth's upper mantle inferred from a multiscale approach[END_REF][START_REF] Raterron | Multiscale modeling of upper mantle plasticity: from single-crystal rheology to multiphase aggregate deformation[END_REF].

The most important prediction of the SO model is the crystallographic spin ġ as a function of the crystal's orientation g, which is what controls the evolution of CPO. Accordingly, this paper focusses on the task of deriving an analytical parametrization of ġ Figure 1. Schematic comparison of the analytical approach and homogenization methods for a single strain increment. Using SO to calculate the spin, ġ, for the nth grain and update its orientation, g, requires knowledge of the strain, and thus spin, of all other crystals in the aggregate necessitating an expensive self-consistent solution. Our analytical approach replaces this information with a record of previous deformation stored as an auxiliary finite strain ellipsoid. This, combined with a handful of other parameters, A, enables rapid calculation of the spin.

that agrees with the SO model predictions. We first note that the total spin ġ is the sum of spins ġ[s] due to the activities of each of the slip systems s = 1, 2, . . . , S within the crystals. We then derive an analytical expression for ġ[s] that depends on the crystal's orientation, the macroscopic strain rate tensor (SRT), the already existing texture and the parameters that characterize the rheology of the slip systems. This expression is then compared, for each slip system separately, with the spins ġ[s] (g) predicted by the SO model for an aggregate of crystals with several simultaneously active slip systems (S > 1). Remarkably, we find that the analytical expression for ġ[s] matches the SO prediction exactly for each slip system s, to within a set of amplitudes A ijkl that can be determined by leastsquares fitting. We uncover surprising symmetries that reduce the number of independent non-zero components of the 'spin' tensor A from 25 to just 2. Finally, we use full SO solutions to determine how these coefficients depend on the relative strengths of the slip systems and on the finite strain experienced by the aggregate.

For irrotational, time-independent deformation, the finite strain ellipsoid (FSE) has the same shape and orientation as the virtual ellipsoid generated by the instantaneous global SRT. In this simplified case, we show that we require only one amplitude. However, when the two ellipsoids are not aligned (see Fig. 1), an extra amplitude is required. We show that predictions of evolving CPO using these analytical parametrizations (which we call ANPAR) are indistinguishable from those of the SO model, and cost only ≈0.01 per cent as much time to compute.

T H E O R E T I C A L P R E L I M I N A R I E S

We begin by reviewing how the orientation and internal deformation of crystals in an aggregate are described mathematically, using the particular case of olivine as an example.

Crystal orientation and orientation distribution

Consider an aggregate comprising a large number N of olivine crystals deforming by dislocation creep. When the aggregate as a whole is subject to a given macroscopic deformation, its constituent crystals respond by deforming via internal shear on a small number S of 'slip systems'. Each slip system s = 1, 2, . . . , S is defined by a unit vector n [s] normal to the slip (glide) plane and a unit (Burgers) vector l [s] parallel to the slip direction. In this study we assume that olivine has three dominant slip systems (010)[100], (001)[100] and(010)[001], corresponding to the indices s = 1, 2, and3, respectively. The degree of anisotropy of an aggregate can be described by specifying for each crystal the three Eulerian angles (φ, θ, ψ) ≡ g that describe its orientation relative to fixed external axes. The definition of these angles that we use is shown in Fig. 2, following [START_REF] Bunge | Texture Analysis in Materials Science[END_REF] convention. The associated matrix of direction cosines a ij is

a i j (g) = ⎛ ⎜ ⎝ cφ cψ -sφ sψ cθ sφ cψ + cφ sψ cθ sψ sθ -cφ sψ -sφ cψ cθ -sφ sψ + cφ cψ cθ cψ sθ sφ sθ -cφ sθ cθ ⎞ ⎟ ⎠ , (1)
where c and s indicate the cosine and sine, respectively, of the angle immediately following. The quantity a ij is the cosine of the angle between the crystallographic axis i and the external axis j.

In the limit as the number of grains N → ∞, the distribution of their orientations can be described by a continuous 'orientation distribution function' (ODF) f (g, t), defined such that f (g, t) dg is the volume fraction of crystals with orientations between g and g + dg at time t. For crystals with triclinic symmetry, the volume of the space of Eulerian angles ('Euler space') required to include all possible orientations is φ ∈

[0, 2π], θ ∈ [0, π ], ψ ∈ [0, 2π]. For olivine, which is orthorhombic, it is sufficient to work in a reduced space (so-called 'irreducible space') φ ∈ [0, π ], θ ∈ [0, π ], ψ ∈ [0, π ],
as the entire Euler space can be retrieved from the irreducible space by symmetry operations. The condition that the total volume fraction of crystals with all possible orientations is unity is then

f (g, t)dg ≡ π 0 π 0 π 0 f (g, t) dφ dψ dθ sin θ = 1, (2)
which implies that f = (2π 2 ) -1 ≡ f 0 for an isotropic orientation distribution.

Kinematics of intracrystalline slip

The time evolution of the ODF is governed by the equation [START_REF] Clement | Prediction of deformation texture using a physical principle of conservation[END_REF])

0 = ∂ f ∂t + ∇ • ( ġ f ) ≡ ∂ f ∂t + ∂ ∂φ ( φ f ) + ∂ ∂ψ ( ψ f ) + 1 sin θ ∂ ∂θ ( θ sin θ f ), (3)
where ( φ, θ, ψ) ≡ ġ is the rate of change of the orientation ('spin') of an individual crystal with orientation g. Eq. ( 3) is a conservation law which states that the rate of change of the volume fraction of crystals having orientations in a small element dg of the Euler space is equal to the net flux of crystal orientations into that element. The [s] , θ [s] , ψ [s] ) such that both the slip direction l [s] and the normal n [s] to the slip plane are in the plane perpendicular to the ψ[s] -axis.

spins ( φ, θ, ψ) are related to the Cartesian components ω i of the spin by

φ = ω 3 + (ω 2 cos φ -ω 1 sin φ) cot θ, (4a) θ = ω 1 cos φ + ω 2 sin φ, (4b) ψ = (ω 1 sin φ -ω 2 cos φ) csc θ. ( 4c
)
Note also that the crystallographic spin ω is just the sum of the externally imposed rotation rate and a contribution ω (c) due to intracrystalline slip, or

ω i = i -i jk l j n k γ ≡ i + ω (c) i , (5)
with ijk the Levi-Civita symbol (component of the permutation tensor).

The spin ġ is the fundamental quantity that will concern us in this study. It depends on the instantaneous macroscopic velocity gradient tensor D, the components of which are

D i j = E i j -i jk k , ( 6
)
where E ij and k are the components of the SRT E and the macroscopic rotation rate of the polycrystal, respectively. When the aggregate is deformed, each crystal within it responds by deforming in simple shear on planes normal to n(g) at a rate γ (g). The local velocity gradient tensor inside the crystal is thus

d i j = γ l i n j . ( 7
)
The local SRT e ij is the symmetric part of d ij , or

e i j = γ 2 (l i n j + l j n i ) ≡ γ S i j , (8)
and in eq. ( 5) is associated with the antisymmetric part of d ij . Here S ij is the Schmid tensor, which resolves the SRT inside each crystal onto the natural frame of reference of the slip system. It is symmetric and traceless and therefore has only five independent components. These can be expressed in terms of generalized spherical harmonics (GSH) T mn l of degree l = 2, where

T mn l (φ, θ, ψ) = e imψ P mn l (cos θ )e inφ , ( 9
)
and P mn l (cos θ) is the associated Legendre polynomial (Bunge 1982, eq. 14.2). Explicit expressions for the independent components of the (slip-system specific) Schmid tensor S

[s] i j for the slip systems s = 1, 2, 3, are given in Appendix A.

Another kinematical object that plays an important role in our theory is the FSE associated with the deformation history experienced by a polycrystal. It is well known in fluid mechanics that an arbitrary time-dependent flow field transforms an initially spherical fluid element of infinitesimal size into an ellipsoid, called the FSE. The shape of the FSE can be characterized by the logarithms of the ratios of the lengths c 1 , c 2 and c 3 of its axes, viz.

r 12 = ln c 1 c 2 ≡ (E 11 -E 22 )t, r 23 = ln c 2 c 3 ≡ (E 22 -E 33 )t, r 31 = ln c 3 c 1 ≡ (E 33 -E 11 )t, ( 10
)
where E 11 , E 22 and E 33 are the principal strain-rates of the SRT E that generates the background texture. Incompressibility of the fluid implies r 12 + r 23 + r 31 = 0, so that only two of the quantities r ij are independent. We also define an 'equivalent strain'

r 0 = √ 2 3 r 2 12 + r 2 23 + r 2 31 1/2 = 2 3 r 2 12 + r 12 r 23 + r 2 23 1/2 . ( 11
)

Slip-system rheology

Following standard practice, we assume that the slip rate γ [s] on each slip system s obeys a power-law rheology of the form

γ [s] ∝ τ τ [s] m [s] -1 τ τ [s] , ( 12
)
where τ is the resolved shear stress (i.e. the shear stress acting on the slip plane in the slip direction), τ [s] is a 'critical resolved shear stress' (CRSS) that measures the inherent resistance of the slip system to slip, and m [s] is a power-law exponent. Although the standard notation is to use n as the stress exponent, we have chosen m in this paper to avoid confusion with all the different occurrences of n. We assume m [s] = 3.5 for all slip systems, following [START_REF] Bai | High-temperature creep of olivine single crystals. 1. Mechanical results for buffered samples[END_REF]. Because the macroscopic deformation rate of the aggregate is specified in our SO calculations, only the ratios of the parameters τ [s] (and not their absolute values) are relevant.

In our calculations we assume τ [1] /τ [2] ∈ [0.25, 4.0] and τ [2] /τ [3] ∈ [0.25, 4.0] (see Table 1). Olivine deforming under upper mantle conditions is typically modelled with values of τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667 [START_REF] Kaminski | D-Rex, a program for calculation of seismic anisotropy in the convective upper mantle[END_REF]) but this range of parameters also captures deformation under the low temperature low pressure (τ [1] /τ [2] = 0.53 and τ [2] /τ [3] = 0.58) and high-pressure high-temperature (τ [1] /τ [2] = 0.73 and τ [2] /τ [3] = 0.39) conditions considered by [START_REF] Castelnau | Microstructures and rheology of the Earth's upper mantle inferred from a multiscale approach[END_REF]. We note that these calculations also permit movement on the (100)[001], {021}[100] and {110}[001] slip systems, which are suppressed in our model. We characterize the CRSS ratios of the dominant slip systems s = 1, 2 and 3 in terms of the variables

p 12 = ln τ [1] τ [2] , p 23 = ln τ [2] τ [3] . ( 13
)
Note also that

p 31 = ln τ [3] τ [1] = -p 12 -p 23 . ( 14
)
In . In our calculations we assume τ [4] /τ [2] = τ [5] /τ [2] = 100 (see Table 1). While these slip systems contribute significantly to the intracrystalline stress, they have a negligible (≈1 per cent) effect on the slip rates of the dominant systems. The model therefore gives valid predictions of the evolution of CPO.

A N A LY T I C A L PA R A M E T R I Z AT I O N

The considerations of the previous section imply that the instantaneous crystallographic rotation rate ġ depends on the crystal's orientation g; the macroscopic SRT E; the already existing texture f; and the parameters p 12 , p 23 and m that characterize the rheology of the slip systems:

ġ = ġ (g, E, f, p 12 , p 23 , m) . ( 15
)
Here ġ is understood as the slip-induced rotation rate, without the contribution due to the macroscopic vorticity which is the same for all crystals and can simply be added to ġ.

Next, we note that the spin components (4) take a particularly simple form when rewritten in terms of slip system-specific Eulerian angles (φ [s] , θ [s] , ψ [s] ) defined so that both the slip vector l [s] and the vector n [s] normal to the slip plane are perpendicular to the ψ[s] -axis (Fig. 2b). The crystallographic spin ġ[s] produced by slip in the direction l [s] on the plane n [s] then has only a single non-zero component ψ[s] , and φ[s] = θ[s] = 0 identically. Fig. 2(b) implies that l

[1] i = a 1i and n [1] i = a 2i , where a ij are given by ( 1) with (φ, θ, ψ) → (φ [s] , θ [s] , ψ [s] ). Ignoring the macroscopic vorticity as explained above, we find that eqs (4) and ( 5) simplify to [s] .

φ[s] = θ[s] = 0, ψ[s] = -γ [s] g
(16)

Thus, the crystallographic spin due to slip is simply the negative of the shear rate on the slip system in question.

To go further, we first note the obvious difficulty that the space of possible background textures f is infinite. To make progress, therefore, we need to restrict and parameterize this space in some way. Our choice is to consider the space of all textures produced by uniform triaxial straining of an initially isotropic aggregate, which can be parameterized by the axial ratios r 12 and r 23 of the associated FSE. Accordingly, the functional dependence we need to determine becomes γ [s] = γ [s] g [s] , E, r 12 , r 23 , p 12 , p 23 , m .

(17)

This still seems impossibly complex, so we now call the SO model to our aid. Consider the case of uniaxial compression along the x 3 -axis at a rate ˙ , for which the nonzero components of the SRT are

E 33 = -˙ , E 11 = E 22 = ˙ /2.
The shear rate γ [1] and the ODF f are then independent of the Eulerian angle φ [1] by symmetry. Fig. 3 shows the spin ψ[1] (θ [1] , ψ [1] ) for the slip system (010)[100] (s = 1) predicted by the SO model with p 12 = p 23 = 0 (τ [1] = τ [2] = τ [3] ) at two different equivalent strains r 0 ≡ | ˙ 3 |t = 0 (Fig. 3a) and r 0 = 0.4 (Fig. 3b). Remarkably, the images of Figs. 3(a) and (b) appear to be the same function with different amplitudes. A more detailed investigation shows that this impression is correct, and that the function in question is F = b sin 2ψ [1] sin 2 θ [1] , where b is an unknown amplitude. Least-squares fitting of this expression to the numerical predictions yields b = 1.25 for Fig. 3(a) and b = 1.71 for Fig. 3(b), with a nearly perfect fit (variance reduction = 99.9 per cent) in both cases.

Next, we note that the function sin 2ψ [1] sin 2 θ [1] can be written as sin 2ψ [1]

sin 2 θ [1] = - 2 √ 3 3 T 20 2 = -2 S [1] 11 + S [1] 22 , ( 18
)
where (Bunge 1982, eq. 14.37), where i = √ -1. This result has two surprising and far-reaching implications. First, the angular dependence of the spin (∝ sin 2ψ [1] sin 2 θ [1] in this case) remains the same regardless of the strength of the background texture; it is only the amplitude of the function that depends on the texture. Second, it suggests that the angular dependence of the spin is always a GSH of degree l = 2, without any contribution from higher-degree harmonics. Noting further that the shear rate γ [s] must depend linearly on the imposed macroscopic strain rate E, we are led to propose the following expression for γ [s] :

T 20 2 is a GSH, defined as T mn l = 2 -1/2 i 1+m-n T mn l -T -m-n l
γ [s] = -ψ[s] = A [s] i jkl (r 12 , r 23 , p 12 , p 23 , m)S [s] i j E kl , ( 19
)
where A is a fourth-order 'spin tensor'. The superscripts [s] denote the index of the slip system (s = 1, 2 or 3). In the next section we determine how the spin tensor A ijkl depends on its five arguments. (θ [1] , ψ [1] ) for the slip system (010)[100], predicted by the SO model with τ [1] /τ [2] = τ [2] /τ [3] = 1 for an initially isotropic olivine aggregate deformed in uniaxial compression. The Eulerian angles θ [1] and ψ [1] are defined as in Fig. 2(b).

PA R A M E T R I Z AT I O N O F T H E S P I N T E N S O R A

SO model calculations

We now use the SO model as a benchmark to determine the tensor components A ijkl (r 12 , r 23 , p 12 , p 23 , m). The procedure comprises two steps: (1) generation of the background texture and (2) calculation of the instantaneous spin induced by applying a given rate of strain to the background texture. Thus in step (1), we first select the number of crystals N in the model aggregate (= 2000 in all cases) and the values of the slip-system parameters p 12 , p 23 , and m (= 3.5 in all cases). We also choose the components of the SRT E that generates the background texture. We work in the reference frame of the FSE, which means that

E = ⎛ ⎜ ⎝ E 11 0 0 0 E 22 0 0 0 E 33 ⎞ ⎟ ⎠ . ( 20
)
This limits the model to orthotropic CPO, which is sound for practical purposes as natural CPO most often exhibit such a symmetry. The SO model is then run starting from an isotropic initial condition until target values of the FSE axial ratios r 12 and r 23 are reached. In step 2, we apply an instantaneous SRT

E = ⎛ ⎜ ⎝ E 11 E 12 E 13 E 12 E 22 E 23 E 13 E 23 E 33 ⎞ ⎟ ⎠ (21)
to the background texture. Note that E need not be the same as E, which allows us to obtain results for arbitrary orientations of the principal axes of E relative to those of E.

The final result of the procedure described above is a set of slip rates γ [s] n on each of the three slip systems (s = 1, 2 or 3 in Table 1) and for each of the n phases (n = 1, 2, . . . , N). The calculated values of γ [s] n are then substituted into eq. ( 5) to obtain the 'partial' spins ω

[s] i due to the actions of the individual slip systems, and which are related to the total spin ω i by

ω i = 3 s=1 ω [s] i . ( 22
)
Finally, by substituting the partial spins ω

[s] i into eq. ( 4) and expressing the results in terms of the slip system-specific Eulerian angles (φ [s] , θ [s] , ψ [s] ) ≡ g [s] , we obtain the rotation rates ψ [s] n for all grains n and slip systems s.

ANPAR model for crystallographic spin

At this point, we have a large library of numerical solutions, but little idea of what they imply about the structure of the function A

[s] i jkl (r 12 , r 23 , p 12 , p 23 , m). As a first simplification, we assume that the dependence on r 12 and r 23 is separable from the dependence on p 12 and p 23 , that is, A [s] i jkl = H [s] ( p 12 , p 23 , m)Q i jkl (r 12 , r 23 ). ( 23)

The first factor H [s] in eq. ( 23) describes how the activities of the three slip systems depend on the CRSS ratios at the initial instant (r 12 = r 23 = 0) of the deformation, while the factor Q ijkl (r 12 , r 23 ) describes how the activities of slip systems with equal strengths (p 12 = p 23 = 0) vary as a function of strain for arbitrary deformations.

Consider first the factor Q ijkl (r 12 , r 23 ). Since E ij and S

[s] i j are both symmetric and traceless, there are at most 25 independent products of them, or equivalently 25 independent Q ijkl . However, we have found that Q ijkl (r 12 , r 23 ) obeys surprising symmetries that reduce the number of its independent non-zero components to just two. We began by fixing H [s] (0, 0, 3.5) = 1 and performing a least-squares fit of the model ( 19) to the spin predicted by the SO model. We did this for 217 different values of (r 12 , r 23 ), where both r 12 and r 23 were in the range [-0.9, 0.9], which is within the range of finite strain that can be successfully modelled by the SO approach. The sampled points were equally spaced along radial lines in the (r 12 , r 23 )-plane. For each of the sampled points we repeated the fit for five random instantaneous SRTs, giving a total of 1085 fits of our model to the numerical solutions of the SO model. This allowed us to discover numerically that 18 of the coefficients Q ijkl were identically zero. We also found at this stage that Q 1122 ≈ Q 2211 . These two numerical results implied that the tensor Q ijkl exhibits major symmetry, that is,

Q ijkl = Q klij . This leaves only six independent, non-zero components of Q ijkl , namely Q 1111 , Q 1122 , Q 2222 , Q 1212 , Q 1313 and Q 2323 .
Relationships among the six remaining non-zero Q ijkl arise from the fact that the labelling of the coordinate axes is arbitrary. We can have a cyclic permutation of the coordinate axes from (1, 2, 3) to (2, 3, 1) or (3, 1, 2), or a non-cyclic permutation from (1, 2, 3) to (1, 3, 2), (2, 1, 3) or (3, 2, 1). The spin ψ[s] has to be invariant under a relabelling of the coordinate axes. Equating the expressions for ψ[s] in the original and the transformed coordinate systems allows us to derive rigourous transformation rules (in the reference frame of the FSE). Setting B = Q 1111 and C = Q 1212 , we find that

Q 1212 (r 12 , r 23 ) = C (r 12 , r 23 ) , (24a) Q 1313 (r 12 , r 23 ) = C (r 31 , r 23 ) , ( 24b
)
Q 2323 (r 12 , r 23 ) = C (r 23 , r 12 ) , ( 24c
)
Q 1111 (r 12 , r 23 ) = B (r 12 , r 23 ) , ( 24d
)
Q 2222 (r 12 , r 23 ) = B (r 12 , r 31 ) , ( 24e
)
Q 1122 (r 12 , r 23 ) = 1 2 [B (r 23 , r 12 ) -B (r 12 , r 23 ) -B (r 12 , r 31 )] . ( 24f
)
These symmetries reduce the number of independent coefficients Q ijkl to just two, which we take to be Q 1111 and Q 1212 . Now consider the factor H [s] in eq. ( 23). We have discovered numerically that H [2] and H [3] can be obtained from H [1] by simple variable transformations:

H [2] ( p 12 , p 23 ) = H [1] (-p 12 , -p 31 ), H [3] ( p 12 , p 23 ) = H [1] (-p 23 , -p 12 ). ( 25
)
Combining this with the previous results for Q ijkl , we obtain the following general ANPAR model for the crystallographic spin, which is valid on each slip system s: + C (r 31, r 23 ) E 31 S [s] 31

ψ[s] = 1 2 H [s] (
. ( 26
)
Note that the coefficient C multiplies the off-diagonal components of E, and is therefore not needed for coaxial deformations where the principal axes of the SRT and the FSE are aligned.

The final part of the ANPAR procedure is to calculate the total spin for each crystal. Substituting the above slip-specific spin (26) into eq. ( 16) provides us with the slip rates γ [s] on each of the three slip systems (s = 1, 2 or 3 in Table 1). The calculated values of γ [s] are then substituted into eq. ( 5) to obtain the Cartesian components of the spin ω i due to the actions of the individual slip systems, and which are related to the total spin ω i by eq. ( 22). Finally, by substituting ω i into eq. ( 4) we obtain the total spin ( φ, θ, ψ) ≡ ġ of an individual crystal, in terms of the Eulerian angles.

Numerical determination of the parametrization coefficients

The symmetries outlined above indicate that we only require analytical expressions for the three functions B, C and H [1] . Full details of the expressions obtained and the methods used are given in Appendix B. Briefly, we first obtain models for B and C, by setting equal slip-system strengths (p 12 = p 23 = 0) and fixing H [1] (0, 0) = 1 in eq. ( 26). We then capture B and C data by a least-squares fit of the model ( 26) to the spin predicted by the SO model, for sampled values of r 12 and r 23 , and random instantaneous SRTs. In each case, the variance reduction of the fit R ≥ 99.7 per cent. Simple polynomials in r 12 and r 23 are fitted to the B and C data (using least squares) to obtain the analytical expressions (B3). The root mean square (RMS) errors of the fits are 0.039 and 0.0070, respectively. Figs 4 and 5 display contour plots of the models (B3) against the B and C data, respectively.

To obtain an analytical expression for H [1] , we capture H [1] data over the entire admissible range of (p 12 , p 23 ) values for olivine. Using eq. ( 13) and the CRSS ratios assumed in Table 1, this gives admissible values of p 12 and p 23 in the range [-ln 4, ln 4] = [-1.386, 1.386]. We fit eq. ( 26) to 81 instantaneous (t = 0, i.e. isotropic CPO) numerical solutions of the SO model for uniaxial compression, with equally spaced points in the (p 12 , p 23 )-plane with p 12 and p 23 in the above admissible range. Simple polynomials in p 12 and p 23 are fitted to the H [1] data using least squares, leading to the analytical expression (B6). The RMS error of the fit is 0.0068.

Finally, we test the assumption (23) that A ijkl can be written as a product of a scalar H [s] (that depends on p 12 and p 23 ) and a tensor Q ijkl (that depends on r 12 and r 23 ). We substituted the analytical expressions for B, C and H [1] (see eqs B3 and B6, respectively) into the full model for the spin on each slip system (26). We then fitted these models to the spin predicted by the SO model for random background textures (formed from various r 12 , r 23 , p 12 and p 23 values) and random instantaneous SRTs. Remarkably, in each case, the variance reduction R > 99.1 per cent and in most cases R > 99.7 per cent.

E V O L U T I O N O F C P O D U R I N G P RO G R E S S I V E D E F O R M AT I O N

The results in the previous sections imply that the ANPAR model provides an accurate and efficient substitute for the much more computationally expensive SO model. We now demonstrate this in more detail by comparing the textures predicted by the two models for olivine polycrystals subjected to various kinds of finite deformation. In the following three test cases (for uniform deformation) the strain increment used ( r 0 = 0.025) is the same for both the SO and ANPAR models. The different components of our method for ANPAR CPO calculation are summarized in Algorithm 1.

Irrotational deformations

Our first test case is a uniaxial compression to a strain r 23 = 0.9, r 12 = 0, with CRSS ratios τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667. Fig. 6 shows the (100), ( 010) and (001) pole figures predicted for this case by the SO model (Fig. 6a) and the ANPAR model (Fig. 6b). The two sets of pole figures are practically indistinguishable (variance reduction R > 99.9 per cent. Variance reduction for pole figures is defined in Appendix C). However, the ANPAR model is a remarkable 1.75 × 10 4 times quicker than the SO model (0.0344 s for ANPAR versus 603 s for SO).

To quantify the agreement in another way, we used the obtained CPO to calculate the effective elastic behaviour of the aggregate for the SO and ANPAR models using the MSAT software [START_REF] Walker | MSAT -A new toolkit for the analysis of elastic and seismic anisotropy[END_REF]). An element-by-element comparison of the two Voigt-Reuss-Hill average elastic stiffness tensors gives a maximum absolute difference of 0.19 GPa between SO and ANPAR, which is not significant for geophysical purposes.

As a second test, Fig. 7 shows the predicted pole figures for uniform deformation by pure shear in the x 1 -x 3 plane to a strain r 12 = r 23 = 0.563, again with CRSS ratios τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667. Again, the two sets of pole figures and the predicted elasticity are nearly identical (variance reduction 99.3 per cent, maximum absolute difference in the predicted elasticity 0.16 GPa). In this case, the speed of the ANPAR model is 3.1× 10 4 greater than that of the SO model (0.0348 s for ANPAR versus 1090 s for SO).

Rotational deformations

Rotational deformations are those in which the axes of the FSE do not remain aligned with the principal axes of the SRT as the deformation progresses.

As an example, consider the case of simple shear, for which the major axis of the FSE is initially aligned with the SRT but then rotates progressively away from it towards the shear plane. As a result, both functions B and C in eq. ( 26) come into play.

Let ˙ 1 be the maximum rate of extension along the x 1 axis. The elongation of the FSE at time t can then be described by the axial ratio R = exp(r 12 ) = exp(2 sinh -1 [˙ 1 t]) (Ribe & Yu 1991, eq. 16). If we denote χ (t) the angle between the two frames, then χ (0) = 0 and lim t → ∞ χ (t) = -π/4. Using the standard tensor transformation rule, we obtain a velocity gradient tensor of the form D = E + W, where

E = ˙ 1 ⎛ ⎜ ⎜ ⎝ cos 2χ sin 2χ 0 sin 2χ -cos 2χ 0 0 0 0 ⎞ ⎟ ⎟ ⎠ , W = ˙ 1 ⎛ ⎜ ⎜ ⎝ 0 -1 0 1 0 0 0 0 0 ⎞ ⎟ ⎟ ⎠ (27)
are the respective strain-rate and rotation-rate, tensors and χ = -1 2 tan -1 (˙ 1 t). We used CRSS ratios τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667. We updated the velocity gradient tensor at each time step to remain in the frame of reference of the FSE. (velocity gradient tensor, using eq. 6) (iv) Set r 0 [k] (strain increment) (v) Calculate r 12 [k], r 23 [k] are shown for critical resolved shear stress ratios τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667. Figure generated using MTEX [START_REF] Bachmann | Texture analysis with MTEX -free and open source software toolbox[END_REF].

[k] = (φ n [k], θ n [k], ψ n [k]) (integrating forward in time) end do (6) Plot texture g n [K ] = (φ n [K ], θ n [K ], ψ n [K ]) (using MTEX)
Fig. 8 shows the pole figures predicted by our theory together with those predicted by the SO model for ˙ 1 = 1 and r 0 = 0.5. Yet again, the two sets of pole figures are nearly indistinguishable, with a variance reduction R = 99.5 per cent and a maximum difference in the predicted elasticity of 0.77 GPa. The speed of the ANPAR model is 5.6 × 10 4 greater than that of the SO model (0.367 s for ANPAR versus 2062 s for SO).

Non-Newtonian corner-flow model for a spreading ridge

Our final example is a more complex and non-uniform geophysical flow, namely the flow in the mantle beneath an ocean ridge. This flow can be simply modeled using the 'corner flow' similarity solution of the Stokes equation in polar coordinates (r, ϕ) [START_REF] Batchelor | An Introduction to Fluid Mechanics[END_REF]. Fig. 9 shows the geometry and boundary conditions appropriate for a ridge crest [START_REF] Lachenbruch | Rise of a variable viscosity fluid in a steadily spreading wedge-shaped conduit with accreting walls[END_REF]. Flow in the asthenosphere 0 < ϕ < α is driven by the horizontal motion of wedge-shaped surface plates at velocity U 0 . The solid lines with arrows show typical streamlines of the flow for an asthenosphere with a power law rheology with power law index n = 3 [START_REF] Tovish | Mantle flow pressure and the angle of subduction: non-Newtonian corner flows[END_REF]. The two streamlines are for ϕ 0 = 10 • and 20 • , and we use α = 60 • throughout this subsection.

The steady incompressible Stokes equations and the boundary conditions in Fig. 9 can be satisfied if the stream function has the self-similar form

= U 0 r F(ϕ), ( 28
)
which is valid for both Newtonian (n = 1) and non-Newtonian (n = 1) fluids. Here we use n = 3, corresponding to a rheology that is close to that of olivine at high stresses (n ≈ 3.5; [START_REF] Bai | High-temperature creep of olivine single crystals. 1. Mechanical results for buffered samples[END_REF].

The function F(ϕ) for n = 3 is of the form

F(ϕ) = A sin ϕ + Ch(ϕ, D), ( 29
)
where

h(ϕ, D) = 27 cos √ 5 3 (ϕ + D) -cos √ 5(ϕ + D) . ( 30
)
The constants A, C and D are chosen to satisfy the boundary conditions, yielding where h ϕ = dh/dϕ. The maximum strain rate ˙ is

D = 3π 2 √ 5 , ( 31
) C = -h(α, D) cos α -h ϕ (α, D) sin α -1 , ( 32
) A = -C h(α, D) sin α + h ϕ (α, D) cos α , ( 33
)
˙ = U 0 |F + F| 2r , ( 34
)
where F = d 2 dϕ 2 F(ϕ), and the local rotation rate (= one-half the vorticity) is

= -U 0 F + F 2r . ( 35
)
To proceed, we require knowledge of the FSE as we progress along a streamline. We obtain the axial ratio R = exp(r 12 ) of the FSE and the orientation χ of the FSE by solving the following evolution equations [START_REF] Kellogg | Mixing and the distribution of heterogeneities in a chaotically convecting mantle[END_REF][START_REF] Ribe | On the relation between seismic anisotropy and finite strain[END_REF]):

Ṙ = 2R (E 11 cos 2χ + E 12 sin 2χ ) , ( 36
) χ = + 1 + R 2 1 -R 2 (E 11 sin 2χ -E 12 cos 2χ ) . ( 37
)
The above equations can be simplified by transforming the Cartesian strain rate components E ij to polar coordinates, and then expressing the time derivatives in terms of a ϕ-derivative (McKenzie 1979, eq. 6):

D Dt = - U 0 F r d dϕ . ( 38
)
This leads to the following simplified form for the evolution equations:

dR dϕ = -R F + F F sin 2 (χ -ϕ) , ( 39
) dχ dϕ = F + F 2F 1 + 1 + R 2 1 -R 2 cos 2 (χ -ϕ) (40)
which must be solved subject to the following initial conditions at ϕ = ϕ 0 : R (ϕ 0 ) = 1, ( 41)

χ (ϕ 0 ) = ϕ 0 + π 4 . ( 42
)
These evolution equations were solved using a fourth-order Runge-Kutta method. In Fig. 9, the FSE is plotted at different points along two different streamlines.

In polar coordinates the velocity gradient tensor is of the form

(McKenzie 1979) D = ⎛ ⎜ ⎝ 0 -2 0 0 0 0 0 0 0 ⎞ ⎟ ⎠ , ( 43
)
where is given by eq. ( 35). However, to calculate CPO, we have to transform back into Cartesian coordinates. In doing this, we obtain

⎞ ⎟ ⎟ ⎠ , W = ˙ ⎛ ⎜ ⎝ 0 1 0 -1 0 0 0 0 0 ⎞ ⎟ ⎠ . ( 44
)
We then use the standard tensor transformation rule to transform the SRT into the reference frame of the FSE. This gives

E = ˙ ⎛ ⎜ ⎝ sin 2(χ -ϕ) cos 2(χ -ϕ) 0 cos 2(χ -ϕ) -sin 2(χ -ϕ) 0 0 0 0 ⎞ ⎟ ⎠ . ( 45
)
Fig. 10 shows the pole figures predicted by our theory together with those predicted by the SO model for an equivalent strain r 0 = 0.6 (r 12 = 1.047, r 23 = -0.523). This was for the first streamline ϕ 0 = 10 • in Fig. 10, with ϕ = 49 • and χ = 67 • . Again, the two sets of pole figures are almost identical, with a maximum difference in the predicted elasticity of 0.93 GPa. When comparing the two pole figures, the variance reduction is 99.0 per cent. In this case, the speed of the ANPAR model is 5.8 × 10 4 greater than that of the SO model (0.384 s for ANPAR versus 2240 s for SO).

C O N C L U D I N G R E M A R K S

The new ANPAR method we describe in this article is an accurate and computationally efficient alternative to existing methods for the simulation of CPO development in olivine. Benchmark tests against the SO self-consistent model (Ponte-Casteñeda 2002) show that ANPAR runs 2-6 × 10 4 times faster, yet predicts textures that are nearly indistinguishable from those predicted by SO. The proposed method is limited to CPO calculations; ANPAR does not tackle the viscoplastic mechanical behaviour of the polycrystal associated with the predicted CPO.

The ANPAR model has some similarities with the D-Rex model of [START_REF] Kaminski | A kinematic model for recrystallization and texture development in olivine polycrystals[END_REF]. In the latter model, the slip rates γ [s] are predicted by minimizing for each grain the misfit between the local and global SRTs. This ad hoc principle yields γ [s] = 2AS [s] i j E i j ,

where A = 1 if global strain compatibility is not enforced and A = 5 if it is. Since the quantities S

[s] i j are GSHs of degree 2, D-Rex agrees with ANPAR concerning the spectral content of the crystallographic spin. However, D-Rex assumes that the spin does not depend on the background texture, and so the amplitude A does not increase as strain accumulates. This is in contrast to the amplitudes B and C in ANPAR, both of which increase strongly with increasing strain in order to satisfy global strain compatibility. The ANPAR model for CPO depends on the deformation history, which agrees with the recent studies of [START_REF] Skemer | The influence of deformation history on the interpretation of seismic anisotropy[END_REF] and [START_REF] Boneh | The effect of deformation history on the evolution of olivine CPO[END_REF].

In constructing the ANPAR model we assumed that the spin tensor A ijkl can be written as the product of a tensor Q ijkl that depends only on the axial ratios r 12 and r 23 of the FSE and a scalar H that depends only on the relative slip system strengths p 12 and p 23 . Although this seems to be a major assumption, the near-perfection of the fits we obtain to the SO predictions appears to justify it.

The simplicity of the ANPAR model is due in part to the orthorhombic symmetry of olivine and the resulting orthogonality of the three dominant slip systems (010)[100], (001)[100] and (010) [001]. This is the reason why the spin consists of only degree-2 GSHs and why the parameter space can be reduced so dramatically. We are investigating the possibility of adapting ANPAR to minerals with non-orthogonal slip systems. We are currently looking at post-perovskite, which is assumed to have slip systems which are aligned with the crystallographic axes plus some slip systems that are oblique to these axes (e.g. Carrez et al. 2007a,b;[START_REF] Merkel | Deformation of (Mg,Fe)SiO 3 postperovskite and D anisotropy[END_REF][START_REF] Metsue | Numerical modelling of dislocations and deformation mechanisms in CaIrO 3 and MgGeO 3 postperovskites comparison with MgSiO 3 postperovskite[END_REF][START_REF] Walker | Elastic anisotropy of D predicted from global models of mantle flow[END_REF]. For example, {110}[001], (001) 110 and {110} 110 have been proposed as sets of slip systems with symmetry equivalents that are oblique to the crystal axes. In our preliminary investigation it appears that ANPAR can be extended in this manner but the non-orthogonal slip systems result in higher-degree harmonics being required in the model for the crystallographic spin. , χ = 67 • , r 12 = 1.047, r 23 = -0.523 (r 0 = 0.6), with CRSS ratios τ [1] /τ [2] = 0.5 and τ [2] /τ [3] = 0.667. The red and blue lines have the same meaning as in Fig. 8.

In its current form the ANPAR method does not easily lend itself to modelling polymineralic aggregates; this will be the subject of future development work. However, the model can easily be extended to other orthorhombic minerals with less than four independent slip systems, such as enstatite, and without having to satisfy the von Mises criterion. [START_REF] Detrez | Effective viscoplastic behavior of polycrystalline aggregates lacking four independent slip systems inferred from homogenization methods; application to olivine[END_REF] show that polycrystalline aggregates lacking four independent systems for dislocation glide can deform in a purely viscoplastic regime only if additional deformation mechanisms (e.g. grain boundary sliding, diffusion, disclinations) are activated, and they assume that the unknown accommodation mechanism can be represented by an isotropic potential. The robust character of our results leads us to suppose that our approach can be generalized to minerals with other symmetries and also to polyphase rocks. If this is the case, the fact that ANPAR is based on the SO model implies that such a generalization could be applicable for a range of materials where first-order homogenization schemes have had limited success and the more computationally taxing SO approach has been considered essential. One such example is modelling of deformation of materials, such as ice, with extreme plastic anisotropy (e.g. [START_REF] Lebensohn | Selfconsistent modelling of the mechanical behaviour of viscoplastic polycrystals incorporating intragranular field fluctuations[END_REF]. Finally, the speed advantage of ANPAR over the SO model holds out the possibility that it could be incorporated efficiently in 3-D and time-dependent simulations of mantle convection.

A C K N O W L E D G E M E N T S

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013)/ERC grant agreement 240473 CoMITAC. AMW is supported by a fellowship from the Natural Environment Research Council (Grant Number NE/K008803/1). We are grateful to two reviewers whose comments have helped us in clarifying the text and also to A. Nowacki for helpful discussions.

Figure 2 .

 2 Figure 2. (a) General definition of the Eulerian angles (φ, θ , ψ). (b) Definition of slip-system specific Eulerian angles (φ[s] , θ[s] , ψ[s] ) such that both the slip direction l[s] and the normal n[s] to the slip plane are in the plane perpendicular to the ψ[s] -axis.

Figure 3 .

 3 Figure 3. Instantaneous crystallographic spin ψ[1] (θ[1] , ψ[1] ) for the slip system (010)[100], predicted by the SO model with τ[1] /τ[2] = τ[2] /τ[3] = 1 for an initially isotropic olivine aggregate deformed in uniaxial compression. (a) | ˙ 3 |t = 0; (b) | ˙ 3 |t = 0.4. Colour scale is in units of the axial shortening rate ˙ 3 < 0. The Eulerian angles θ[1] and ψ[1] are defined as in Fig.2(b).

  Figure 3. Instantaneous crystallographic spin ψ[1] (θ[1] , ψ[1] ) for the slip system (010)[100], predicted by the SO model with τ[1] /τ[2] = τ[2] /τ[3] = 1 for an initially isotropic olivine aggregate deformed in uniaxial compression. (a) | ˙ 3 |t = 0; (b) | ˙ 3 |t = 0.4. Colour scale is in units of the axial shortening rate ˙ 3 < 0. The Eulerian angles θ[1] and ψ[1] are defined as in Fig.2(b).

Figure 4 .

 4 Figure 4. Spin amplitude B as a function of deformation when the strengths of the dominant slip systems are all equal, shown as a function of the axial ratios of the finite strain ellipsoid. Solid contours show the amplitude B that best fits the predictions of the SO model, and the dashed contours show the fitting function (B3).

Figure 5 .

 5 Figure 5. Same as Fig. 4, but for the spin amplitude C.

Algorithm 1 .

 1 ANPAR CPO calculations for Olivine.(1) Set N (Number of grains (n = 1, . . . , N)) (2) Set initial isotropic texture g n [0] = (φ n [0], θ n [0], ψ n [0]). (3) Set p 12 , p 23 (CRSS ratios) (4) Set r 12 [0] = r 23 [0] = 0.(FSE initially a sphere) (5) For k = 1, . . . , K do (Calculate texture at each time-step) (i) Set E[k] and[k] (SRT and macroscopic rotation rate) (ii) Transform SRT into reference frame of FSE (iii) Calculate D[k]

Figure 6 .

 6 Figure 6. Pole figures (equal-area projections) for an olivine polycrystal deformed by uniaxial compression to a strain r 23 = 0.9, r 12 = 0. The arrows indicate the compression (x 3 -) axis, which extends from the bottom to the top of each figure. The predictions of (a) the SO model and (b) the analytical model (ANPAR)are shown for critical resolved shear stress ratios τ[1] /τ[2] = 0.5 and τ[2] /τ[3] = 0.667. Figure generated using MTEX[START_REF] Bachmann | Texture analysis with MTEX -free and open source software toolbox[END_REF].

Figure 7 .

 7 Figure 7. Same as Fig. 6, but for deformation by pure shear to r 12 = r 23 = 0.563. The axes of maximum extension (x 1 ) and compression (x 3 ) are indicated by the arrows.

Figure 8 .

 8 Figure 8. Same as Fig.6, but for deformation by simple shear to r 12 = 0.841, r 23 = -0.421 (r 0 = 0.5), with CRSS ratios τ[1] /τ[2] = 0.5 and τ[2] /τ[3] = 0.667. The black, red and blue lines indicate the axis of maximum instantaneous extension (x 1 ), the shear plane and the long axis of the FSE, respectively. The direction of shear is indicated. a velocity gradient tensor of the form D = E + W, where

φFigure 9 .

 9 Figure 9. Geometry and boundary conditions for the corner flow model of a ridge crest (right half only). The half spreading rate is U 0 and the asthenospheric wedge has angular dimension α = 60 • . Solid lines with arrows show typical streamlines for a power-law rheology with index n = 3. The FSEs are plotted at different points along the two streamlines. The black ellipse relates to the point in the flow at which the texture is approximated (see Fig. 10).

  Figure 9. Geometry and boundary conditions for the corner flow model of a ridge crest (right half only). The half spreading rate is U 0 and the asthenospheric wedge has angular dimension α = 60 • . Solid lines with arrows show typical streamlines for a power-law rheology with index n = 3. The FSEs are plotted at different points along the two streamlines. The black ellipse relates to the point in the flow at which the texture is approximated (see Fig. 10).

  Figure 9. Geometry and boundary conditions for the corner flow model of a ridge crest (right half only). The half spreading rate is U 0 and the asthenospheric wedge has angular dimension α = 60 • . Solid lines with arrows show typical streamlines for a power-law rheology with index n = 3. The FSEs are plotted at different points along the two streamlines. The black ellipse relates to the point in the flow at which the texture is approximated (see Fig. 10).

Figure 10 .

 10 Figure10. Same as Fig.6, but for the corner-flow model of a spreading ridge shown in Fig.9. Textures estimated by the SO and ANPAR models at the point of the filled FSE in Fig.9. α = 60 • , ϕ 0 = 10 • , ϕ = 49 • , χ = 67 • , r 12 = 1.047, r 23 = -0.523 (r 0 = 0.6), with CRSS ratios τ[1] /τ[2] = 0.5 and τ[2] /τ[3] = 0.667. The red and blue lines have the same meaning as in Fig.8.

Table 1 .

 1 Slip systems.

	Index s	Slip plane	Slip direction	τ [s] /τ [2]	Exponent
	1	(010)	[100]	0.25-4.0	3.5
	2	(001)	[100]	1.0	3.5
	3	(010)	[001]	0.25-4.0	3.5
	4	(101)	[101]	100.0	3.5
	5	( 1 0 1)	[101]	100.0	3.5

A P P E N D I X A : S L I P S Y S T E M -S P E C I F I C S C H M I D T E N S O R S

The Schmid tensor S is slip system specific, and thus different for each system. Let S [s] i j be the Schmid tensor for slip system s. For the slip system (010) [100] (s = 1), the relationships between the independent components of the Schmid tensor and the Eulerian angles are S [1] 11 -S [1] 22 = -1 2 cos 2φ sin 2ψ cos 2 θ + 1 + 2 sin 2φ cos θ cos 2ψ , (A1a)

S [1] 12 = -1 4 sin 2φ sin 2ψ cos 2 θ + 1 -2 cos 2φ cos θ cos 2ψ , (A1b)

S [1] 13 = -1 2 sin θ (sin φ cos θ sin 2ψcos φ cos 2ψ) , (A1c)

S [1] 23 = 1 2 sin θ (cos φ cos θ sin 2ψ + sin φ cos 2ψ) , (A1d)

Similarly, for the slip system (001)[100] (s = 2) we obtain

11 -S [2] 22 = sin θ (cos 2φ cos θ sin ψ + sin 2φ cos ψ) , (A2a)

S [2] 13 = -1 2 (sin φ cos 2θ sin ψcos φ cos θ cos ψ) , (A2c)

)

Finally, for the slip system (010)[001] (s = 3) we find

11 -S [3] 22 = sin θ (cos 2φ cos θ cos ψsin 2φ sin ψ) , (A3a)

)

)

S [3] 11 + S [3] 22 = -1 2 sin 2θ cos ψ. (A3e)

A P P E N D I X B : A M P L I T U D E O F T H E C RY S TA L L O G R A P H I C RO TAT I O N R AT E

In this appendix, we quantify the dependence of the slip system amplitudes A [s] i jkl on the strains r 12 and r 23 and the CRSS ratios p 12 and p 23 . We assume that A can be expressed more compactly, as shown in eq. ( 23).

B1 Limit A [s]

i jkl (r 12 , r 23 , 0, 0, 3.5) We first consider how the activities of slip systems with equal strengths (p 12 = p 23 = 0) vary as a function of strain. This enables us to find out how the Q ijkl depend on the parameters r 12 and r 23 that characterize the FSE. As we pointed out in Section 2.3, the value of m = 3.5 is assumed for all slip systems.

We explained in Section 4.2 the symmetries that allow us to reduce the number of independent non-zero components Q ijkl from 25 to just 2. The transformation rules (24), derived by noting that the labelling of the coordinate axes is arbitrary, are given in their full form here:

We also discovered, numerically, the symmetry condition

This means, for example, that the values of Q ijkl for uniaxial extension are identical to those for uniaxial compression. Using eq. ( B2) and the complete symmetry transformations (B1), we can reduce by a factor of eight the size of the parameter space (r 12 , r 23 ) that we have to explore numerically to determine how B and C depend on r 12 and r 23 .

We determine values of B and C at sampled points in the (r 12 , r 23 )-plane by fitting eq. ( 26) (with H [s] (0, 0) = 1) to the spin predicted by the SO model, using a standard least-squares procedure. These calculations yield the curves shown by solid lines in Figs 4 and5. In each case the variance reduction R > 99.7 per cent. Then, we fitted B and C data, obtained by the above method, to fourth-order polynomials of the form where

The RMS errors of the fits to the B and C data are 0.039 and 0.0070, respectively. The values of the non-zero coefficients z mn are given in Table B1, for both fitting functions. The quality of these fits can be viewed in Figs 4 and5.

B2 Limit A [s]

i jkl (0, 0, p 12 , p 23 , 3.5) We next consider how the amplitudes depend on p 12 and p 23 at the initial instant (r 0 = r 12 = r 23 = 0) of the deformation. This allows us to discover how the functions H [s] in eq. ( 23) depend on the parameters that characterize the relative strength of the slip systems. Here, we consider values of p 12 and p 23 in the range [-ln 4, ln 4] ≈ [-1.386, 1.386]. We discovered numerically the following transformation rules:

H [2] ( p 12 , p 23 ) = H [1] (-p 12 ,p 31 ) = H [1] ( p 23 , p 31 ) (B5b)

H [3] ( p 12 , p 23 ) = H [1] (-p 23 ,p 12 ) = H [1] ( p 31 , p 12 ). (B5c)

These transformations enable us not only to express H [2] and H [3] in terms of H [1] , but also to reduce the size of the parameter space (p 12 , p 23 ) that we have to explore numerically to determine how H [1] depends on p 12 and p 23 .

We obtain H [1] data by calculating, via least-squares, the amplitudes H [1] that fit eq. ( 26) to 81 instantaneous numerical solutions of the SO model for uniaxial compression, with equally spaced points in the (p 12 , p 23 )-plane with both p 12 and p 23 in the range [-1.386, 1.386]. We then fit quadratic polynomials, satisfying the above relations (B5), to the collected H [1] data. The results are H [1] ( p 12 , p 23 ) = 1 -0.0295 p 12 -0.0130 p 23 -0.00743 p 2 12 -0.00347 p 12 p 23 -0.00333 p 2 23 .

(B6)

The RMS error of the fit is 0.0068. We apply the transformations (B5) to form analytical expressions for H [2] and H [3] , respectively.

B3 General case

We construct a general form for A [s] i jkl (r 12 , r 23 , p 12 , p 23 ) that is consistent with the above limiting cases. We first substitute the model (B6) into eq. ( 23). We then use the coefficients in eq. ( B4) and Table B1 to formulate the models (B3) for B and C. The B and C models are substituted into eqs (24a)-(24f) to form Q ijkl , which is subsequently used in eq. ( 23). The resulting expression is the one we use in all our numerical calculations.

We have fitted this model to the spin predicted by the SO model for random background textures (formed using various r 12 , r 23 , p 12 and p 23 values) and random instantaneous SRT's. Remarkably, in each case, the variance reduction R > 99.1 per cent and in most cases R > 99.7 per cent.

A P P E N D I X C : VA R I A N C E R E D U C T I O N B E T W E E N P O L E F I G U R E S

To calculate the variance reduction between the pole figures shown in Section 5, we use the transformations

to map the Euler angles (φ n , θ n , ψ n ) of each grain onto an 'Euler cube' φ * n , θ * n , ψ * n . The Euler cube has a uniform metric, and each of its sides is of length 2.0. Let the set of Euler angles for each grain for the SO and ANPAR approximations be denoted by (φ * (n,SO) , θ * (n,SO) , ψ * (n,SO) ) and (φ * (n,AN) , θ * (n,AN) , ψ * (n,AN) ), respectively. For each grain the 'distance', d n , between these two sets of Euler angles is calculated by


S Marié

X Gloerfelt

Dynamical selective ltering for the Lattice Boltzmann Method

In this study, a new selective ltering technique is proposed for the Lattice Boltzmann Method. This technique is based on dynamical implementation of the selective lter coefcient σ. The proposed model makes the latter coecient dependent on the shear stress in order to restrict the use of the spatial ltering technique in the sheared stress region where numerical instabilities may occur. Dierent parameters are tested on a 3D decaying Taylor Green Vortex and compared to the classical static ltering technique and to the use of a standard subgrid-scale model. I.

Introduction

The Lattice Boltzmann Method [START_REF] Skordos | Initial and Boundary Conditions for the lattice Boltzmann Method[END_REF] (LBM) is used to simulate the behavior of weakly compressible ows. Some studies [START_REF] Marié | Comparison between Lattice Boltzmann Method and Navier-Stokes high order schemes for Computational Aeroacoustics[END_REF] have shown the capabilities of LBM to perform aeroacoustic simulations due to its low dissipation error. However, as a counterpart, LBM suers from numerical instabilities when Reynolds number becomes high. Moreover, numerical instability waves are often generated by unadapted initial conditions, geometric singularities or in region where large numerical approximations are done. In industrial applications several of these numerical defaults are present, thus computations often become catastrophically unstable.

Various studies have proposed stabilization techniques. 47 Lots of stabilizing strategies have a global eect on the viscosity thus modifying the eective Reynolds number. Another eect of stabilizing methods is to increase the global dissipation then impacting pressure uctuation from which it could be important to keep an accurate estimation (e.g. for aeroacoustic simulations). Ricot et al. [START_REF] Ricot | Lattice Boltzmann Method with selective viscosity lters[END_REF] proposed to use selective spatial lters [START_REF] Bogey | A family of low dissipative explicit schemes for ow and noise computations[END_REF] to stabilize the method by increasing the dissipation for high wavenumbers where the LBM instabilities occur and then keeping a low dissipation at small wavenumbers. This approach can be applied to aeroacoustic simulations by maintaining an acceptable level of dissipation error at low wavenumbers. However, this method basically applies to the whole domain and becomes useless outside of sheared region where numerical instabilities have less chances to develop. Furthermore, the use of high-order lters have the eect of increasing the stencil of LBM which is low by nature and lead to a loss of locality of the method which is known to be one of the main advantages of the method. Finally, from a dynamical point of view, the selective spatial lters have never been tested on transition situation where the dynamical evolution prediction is of major importance for accuracy of the results. Therefore, the need for improved stabilization procedure is justied and should be carried out, in particular in the framework of the Lattice Boltzmann Method.

The idea of the present study is to propose an improved ltering strategy restricted to highly sheared region keeping low-sheared ones free of articial dissipation. The proposed strategy is developed in the framework of the Lattice Boltzmann Method and applied to a dynamical test-case highly sensitive to the numerical dissipation. After a brief presentation of the Lattice Boltzmann model, the new ltering strategy is described and validated on a 3D decaying Taylor Green Vortex and then compared to the usual ltering techniques.

II.

Model

II.A. Lattice Boltzmann Method

The Lattice Boltzmann method 1 is described by the following algorithm:

g α (x + c α , t + 1) = g α (x, t) - 1 τ g [g α (x, t) -g eq α (x, t)] (1)
where g α are distribution functions computed on a regular velocity lattice c α , colliding and relaxing to a local equilibrium g eq α with a relaxing parameter

τ g = ν c2 0 + 1 2
where ν and c0 being the adimensional viscosity and speed of sound respectively. In this study, we use the D3Q19 model described by the following parameters:

g eq α (x, t) = ρω α (1 + u.c α c2 0 + (u.c α ) 2 2c 4 0 - |u| 2 2c 2 0 ) (2)
ω α =              1 3 , α = 0 1 18 , α = 1 → 6 1 36 , α = 7 → 18 c2 0 = 1 3 (3)
It should be noticed that the construction of Eq.( 1) enforces the space and time discretization to be linked by the relation:

∆t = c0 ∆x c 0 (4)
This enforces the CFL number (c 0 ∆t/∆x) to be the same for each simulations (CFL= c0 ). The macroscopic quantities ρ and u can be computed from the distribution functions with the discrete momenta:

ρ = α f α (5) ρu = α c α f α (6)
The pressure is recovered by the relation:

p = c2 0 ρ (7)
Based on these parameters, it can be shown 2 that LBM simulates the 3D compressible Navier-Stokes equations in the limit of low Mach numbers with a 2nd order accuracy in space and time.

II.B. Dynamical selective spatial lters

As proposed by Ricot et al., [START_REF] Ricot | Lattice Boltzmann Method with selective viscosity lters[END_REF] the stability of LBM can be enhanced by space-ltering the momenta of eq.( 5) and (6). Spatial ltering of a quantity Q is dened by subtracting a weighted combination of the symmetric neighboring points:

Q(x) = Q(x) -σ D j=1 N n=-N d n Q(x + n∆x j ) (8)
where σ is a coecient between 0 and 1 often taken to 0.1, d n are coecients depending on the lter order and D is the number of spatial dimension. In this study, classical 3-point, 5-point stencil and optimized 9-point stencil lters are used. [START_REF] Bogey | A family of low dissipative explicit schemes for ow and noise computations[END_REF] The coecients of the lters are summarized in table 1.

The idea of this study is to make the coecient σ dependent on the shear stress. For instance, let's consider σ d (x) to be of the form:

σ d (x) = σ 0 1 -e -(|S(x)|/S0) 2 2 (9) 2 of 8 d 0 d 1 d 2 d 3 d 4 SF-3 1/2 -1/4
SF-5 where |S| = 2S ij S ij , σ 0 is the static lter amplitude and S 0 is the shear stress amount from which the dynamical lter start to be active. Thus, when the shear stress is low, the lter has no eect (σ d ∼ 0) and when shear stress rises to higher values, the lter acts normally (σ d ∼ σ 0 ). Then S 0 should be seen as a sensitivity parameter and its inuence will be further discussed in section III.B.2. The quantity |S| involves spatial derivatives that must be computed during the simulations. From a Lattice Boltzmann point of view, this quantity could be computed from the second order momentum:

τ ij = 2ρνS ij = - α c α,i c α,j (g α -g eq α ) (10)
which gives:

|S| = Q f 2ρν (11)
with

Q f = 2P ij P ij and P ij = α c α,i c α,j (g α -g eq α ).
Relation ( 11) is often use for the implementation of subgrid scale models in the Lattice Boltzmann Method. [START_REF] Dong | Inertial consistent subgrid model for Large Eddy Simulation based on Lattice Boltzmann Method[END_REF] In the following, the present model is tested on a 3D decaying Taylor-Green vortex and compared to the classical static lters and to the Smagorinsky subgrid-scale model.

III.

Taylor Green Vortex

III.A. Test case implementation

In order to study the eect of the present ltering technique, the decaying Taylor-Green vortex is used.

It is a fundamental case used as prototype for vortex stretching and production of small-scale eddies and therefore allows the study of the dynamics of transition to turbulence. This test-case has been widely used to study the dissipation errors of numerical schemes. [START_REF] Fauconnier | On the performance of relaxation ltering for large-eddy simulation[END_REF] In particular, Aubard et al. [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective ltering strategy for large-eddy simulations[END_REF] have recently used this test-case to confront the selective ltering techniques to the use of subgrid-scale models.

The initialization of the Taylor-Green vortex is done by setting velocity and pressure variables as follows:

                   p = p ∞ + ρ ∞ U 2 ∞ 16 [cos(2z) + 2][cos(2x) + cos(2y)] u = U ∞ sin(x) cos(y) cos(z) v = -U ∞ cos(x) sin(y) cos(z) w = 0 (12)
In order to reduce numerical oscillations at the beginning of the simulation, the distribution functions g α are initialized to their equilibrium state with an additional non-equilibrium part based on the Chapman-Enskog micro-scale expansion. More discussions about LBM initialization are detailed in Skordos. [START_REF] Skordos | Initial and Boundary Conditions for the lattice Boltzmann Method[END_REF] For this study, the Mach number is taken to M ∞ = 0.085 and ρ ∞ = 1 xing the other parameters to U ∞ = 0.049 and p ∞ = 1/3 in lattice unit. The relaxation parameter is set with the Reynolds number Re to

τ g = U ∞ ∆xc 2 0 Re + 1 2
. For the validation of this test case, the spectral data from Brachet et al. [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF] are used and compared to our reference simulation on a 256 3 grid.

III.B. Results and discussion

All the simulations of this study are performed on 2π-periodic cubic domain with a Reynolds number of Re = 1600. A reference simulation is performed on a 256 3 grid without any ltering technique. Figure 1 displays the evolution of the sheared quantities in the domain. The classical behaviour of the Taylor-Green vortex is observed, the initial eld gives rise to large vortices which are then stretched and lead to the production of small-scale eddies and decaying turbulence. In the following, a series of tests are performed to characterize the present ltering strategy. First, the inuence of the grid resolution on the non-ltered scheme is presented, then the inuence of parameters σ 0 and S 0 in equation ( 9) are scrutinized and nally, a comparison with the static ltering strategy and the subgrid-scale model is performed.

III.B.1. Inuence of the resolution

The lattice Boltzmann method is a 2nd-order accuracy scheme in space and time. In the literature, few studies have been published on the validation of the 3D Taylor-Green vortex with LBM simulations. Then it could be interesting to see the capabilities of the classical scheme to simulates the dynamics of a 3D decaying vortex. For that purpose, the evolution of the kinetic energy dissipation rate = -∂ t k is scrutinized for various resolution from 64 3 to 256 3 and compared with the spectral simulation of Brachet et al. [START_REF] Brachet | Small-scale structure of the Taylor-Green vortex[END_REF] The TGV dynamical evolution is characterized by three main steps visible in the time trace of . First, the initial laminar state is transitioning to turbulence until the stretched vortex tubes break down into small-scales around t = 5. Then the dissipation rate is rising to a sharp peak near t = 9 corresponding to the full turbulent state which is then decaying similarly to an isotropic and homogeneous turbulence.

The results of the classical LBM simulations without any technique are displayed in gure 2. The reference simulation on a 256 3 grid is observed to be in very good agreement with the spectral results. The 128 3 grid gives satisfactory results for the transition region but the peak in the dissipation rate is not fully reached and the decaying phase starts too early. Then the 64 3 and 96 3 simulations give rise to basic instability at dierent time. The 96 3 grid simulation exhibits a relatively good transition to turbulence but collapses just before the peak of dissipation around t = 8. Finally, the 64 3 grid simulation collapses around t = 5 when the stretched vortex tubes break down into small-scales.

Thus, the LBM scheme shows a good ability to simulate the Taylor-Green vortex dynamics when using ne grids but is limited by its inherent instability for coarser grids. Then the use of stabilizing ltering strategy becomes a crucial point to investigate turbulent simulation with the Lattice Boltzmann Method.

III.B.2. Tuning parameters

As described in section II.B, the dynamical ltering strategy is sensitive to two main parameters which are σ 0 and S 0 . The rst one determines the ltering amplitude when the shear stress is high and could be seen as the eciency of the lter. The latter one plays the role of shear stress selectivity and is essential for controlling the dissipation amount of the ltering technique.

The sensitivity to these parameters is then investigated on a 96 3 grid for various values of σ 0 , S 0 and the dierent selective lters presented in table 1. The nomenclature used to distinguish the simulation is based on a ssf f dd name where s = 100σ 0 , f is the number of points of the selective lter and d = 100S 0 . For example, the s5f 3d100 simulation refers to the use of the 3-points lter with σ 0 = 0.05 and S 0 = 1 in lattice unit. A rst sensitivity analysis is performed on the σ 0 parameter with a xed value of S 0 = 1. Results are reported in gure 3-left and exhibit relatively similar behaviour. A high value of σ 0 (thick dotted lines) induces more damped results, in particular after the vortex breakdown near t = 5 when the lter order is low (red curves). On the other hand, better results are obtained for low σ 0 values (thin lines). This improvement is particularly important for the 3-point lter which can handle the vortex breakdown around t = 5 when σ 0 is set to 0.01. Lower values of σ 0 have been tested but lead to unstable simulations. This suggest that the dynamical ltering strategy should be applied for values of σ 0 as low as possible in the limit of stability. By looking at gure 3-right where the time evolution of σ d is represented, one can see that the value of σ 0 has a direct impact on the time evolution of σ d . Indeed, when σ 0 changes, the ltering amount is modied and so is the shear stress which modies the local σ d value. Moreover, for higher values of σ 0 the maximum value of σ d is never reached suggesting that the ltering amount induces lower shear stress and then reducing σ d . As a counterpart, when σ 0 is low, the ltering amount allows larger shear stress values then increasing σ d to its maximum value xed by σ 0 . However, these considerations depends on the shear stress selectivity which is driven by the value of S 0 .

The second sensitivity analysis is performed on the S 0 parameter with a xed value of σ 0 = 0.05. Results are reported in gure 4-left. It shows that the selected S 0 values have an important impact on the results. A lower sensitivity (dashed lines) induces an earlier ltering which damp the initial laminar state and lead to a wrong dissipation rate evolution. Indeed, results of gure 4-right indicates that σ d starts to increase from a not zero value for low sensitivity lters. In contrast, when the sensitivity is high (solid lines), the lter does not act in low sheared region and is activated only when large vortices have broken up and small-scales structures are rising to turbulence. Moreover it could be seen from gure 4 that the impact of the lter order is reduced when increasing the sensitivity. Then, 3-points lter results are close to these of the 9-point lter when sensitivity is high (thick solid lines). This is an important consequence in terms of computational cost suggesting that the lter order could be reduced when increasing sensitivity.

These tuning tests have emphasized the important role played by the couple (σ 0 , S 0 ) on the simulation of a 3D decaying Taylor-Green vortex. In the following simulations, this couple will be set to (0.05, 1.).

III.B.3. Comparison with static lters

The static ltering technique could be seen as a particular case of the dynamical one when the sensitivity is set to zero (S 0 = 0). Then, static lters act everywhere with the same amount, with potential damping of important structures in the dynamical evolution. With the above discussion, the static ltering technique is then expected to give dissipated results. This is conrmed by looking at gure 5 which displays the comparison of static and dynamical ltering on a 96 3 grid (left) and a 128 3 grid (right). An important result is observed for 3-point stencil lter which gives completely wrong behaviour with a static strategy and exhibits results close to higher order lters for the dynamical strategy. This is particularly observed in the transition region where 3-point dynamical lter gives better results than the static 5-point lter and similar results with the static 9-point lter. Moreover, the dynamical ltering appears to better predict the dissipation peak around t = 9 than the static lters. The above discussions thus indicates that the results could be improved by increasing the ltering sensitivity and decreasing the σ 0 . Then the 96 3 grid simulation with σ 0 = 0.01 and S 0 2 is added for indication on gure 5-left and gives satisfactory results. This behaviour suggests that the wavenumber selectivity of the spatial lters plays a minor role compared to the shear stress selectivity. Indeed, the results indicates that applying a high-order wavenumber selectivity on non-sheared quantities (i.e. with static lters) could lead to over-dissipation in the transition region. As a counterpart, by introducing a shear-stress dependency on the lters, the wavenumber selectivity is forced to act on turbulent quantities with an important shear amount which corresponds to large-scale structures on which the wavenumber selectivity has no impact and could be reduced by decreasing the lter order. This result is of major importance in the framework of lattice Boltzmann method showing that relevant results could be obtained with a 3-point stencil selective lter which correspond to the LBM stencil, then conserving the locality of the scheme particularly important when dealing with heavy parallel implementations.

III.B.4. Comparison with SGS models

The Lattice Boltzmann implementation of subgrid-scale model is very close to the one of the present dynamical ltering. Relation ( 11) is used to estimate the eddy viscosity and thus the relaxation time τ g . The comparison between the implementation of the classical Smagorinsky model and the present methodology is presented on gure 6 where turb1 refers to the Smagorinsky model with a constant C s = 0.1 and turbSm with a constant 0.18. The results indicates that the SGS implementation appears to be too dissipative and describes poorly the laminar-turbulent transition region. Similar observations are detailed in Aubard et al. [START_REF] Aubard | Comparison of subgrid-scale viscosity models and selective ltering strategy for large-eddy simulations[END_REF] where various SGS strategy has been compared and found to be not adapted to this type of transition test case. The present dynamical strategy is found to better predict the transition region by ltering only the sheared region. The main dierence between the present implementation and the SGS methodology relies in the amount of ltering when shear stress is lower than the imposed sensitivity (S < S 0 ). Indeed for those values the lter have no impact and let the turbulent structures free of numerical dissipation. As a counterpart, the SGS implementation is directly proportional to the amount of shear stress and the eddy viscosity have a non zero value for shear stress amount close to S 0 thus imposing an over-dissipation for the concerned structures which are of primary importance for the dynamical evolution. Moreover, the SGS approach is not a bounded procedure and the eddy viscosity can reach arbitrary high value when highly sheared region are encountered whereas the presented strategy is limited by the σ 0 value restricting the numerical dissipation to a limited amount. Finally, from a general point of view, the SGS strategy could be seen as the present dynamical ltering technique with a low shear stress sensitivity and should be linked to the results of section III.B.2.

IV. Conclusion

A dynamical ltering technique for the lattice Boltzmann method has been presented and tested on a 3D decaying Taylor-Green vortex. It has been shown that the use of selective spatial lters with a coecient based on the amount of shear stress led to improved dynamical predictions of the 3D decaying Taylor-Green 7 of 8 vortex. In particular, it has been emphasized that the shear stress selectivity was the action of the lters to localized zones then reducing the global amount of numerical dissipation. The choice of shear stress selectivity has been motivated by the lattice Boltzmann framework for which the shear stress is a relevant quantity that can be accessed with a minimum amount of additional computational time. The results are particularly improved when the lter order is low, especially for the comparison with static selective lters suggesting that the wavenumber selectivity is dominated by the one of the shear stress. The comparison of the presented methodology to the classical subgrid-scale model methodology on a laminar-turbulent transition test case such as the Taylor-Green vortex have also led to promising results highlighting the importance of shear-stress selectivity in the prediction of the turbulence dynamical evolution. From a computational cost point of view, the presented strategy have shown interesting capabilities when using low-order lters thus reducing the eective stencil to the one of the LBM.

Then the dynamical ltering strategy should be seen as an enhanced stabilization procedure for the lattice Boltzmann method where the amount of numerical dissipation could be controlled in space. The presented methodology have been tested on a specic test-case and should be improved by making the sensitivity S 0 independent of the test case. Further studies must be carried out to overcome these limitations through a global stability analysis. Moreover, additional investigations must be carried out to characterize the dynamical ltering behaviour in the case of wall-bounded ows.

Finally it should be noticed that the presented methodology could be applied to a wide variety of numerical problems where only local dissipation is required. Moreover, the presented procedure should be extended to dierent physical problems by modifying the sensitive quantity and not restrict it to the shear stress. Furthermore, the presented dynamical ltering is not restricted to the Lattice Boltzmann method and could be applied to the classical nite dierences schemes.

Figure 1 .

 1 Figure 1. Isosurface of the Q-criterion colored by kinetic energy at time t = 0,t = 4,t = 10,t = 16 for Re = 1600 on a 128 3 grid.

Figure 2 .

 2 Figure 2. Non dimensional time evolution of the dissipation rate .

Figure 3 .

 3 Figure 3. Inuence of the σ0 parameter. Non dimensional time evolution of the dissipation rate (left) and maximum value of normalized dynamic coecient max[σ d (x)]/σ0 (right).

Figure 4 .

 4 Figure 4. Inuence of the S0 parameter. Non dimensional time evolution of the rate (left) and maximum value of normalized dynamic coecient max[σ d (x)]/σ0 (right).

Figure 5 .

 5 Figure 5. Comparison between static lters (solid lines) and dynamical lters (dashed lines) on a 96 3 domain (left) and a 128 3 domain (right).

Figure 6 .

 6 Figure 6. Comparison between dynamical ltering strategy and subgrid-scale model on a 96 3 domain (left) and a 128 3 domain (right).

Table 1 .

 1 Coecients of the selective lters: 8 dn = d-n.

	6/16	-4/16	1/16
	Bogey-9 (B9) 0.243527493120 -0.204788880640 0.120007591680 -0.045211119360 0.008228661760

of 8 American Institute of Aeronautics and Astronautics

of 8 American Institute of Aeronautics and Astronautics

of 8 American Institute of Aeronautics and Astronautics