E Moreau

Survey of roughness by stochastic oscillations

Keywords: numbers: 05.40.-a, 02.50.Ey, 68.35.Ja

In this paper, connection between surface roughness and directed polymers in random medium are studied, when the surface is considered as a directed line undergoing stochastic oscillations. This is performed by studying the influence of a stochastic elastic forcing term -κy +η(s), on a particle moving along a rough surface. Two models are proposed and analysed in this way: the random-walk process (RW) in its discrete and continuous form, and a Markovian process via the Ornstein-Uhlenbeck (O-U) process. It is shown that the continuous RW leads to an oscillator equation, via an effective action obeying a KPZ equation which is solved analytically. The O-U process allows to obtain information on the profile of surface for a long size substrate. The analogy with the roughness is achieved by introducing a quantity suited to directed line formalism: the height velocity variation ∂h/∂s.

I. INTRODUCTION

Characterization of surface roughness or interface growth is of first importance in many applications: microelectronics, preparation of catalysts, magnetic materials (see e.g. [1]- [4]). And numerous models based on numerical as well as analytical approaches have been proposed over past decades to treat interface properties. These works concern in particular the KPZ equation (scale laws, roughness, growth) [5]- [15], and models treating connections between surface growth and directed polymers in random medium [16]- [24]. The present paper focus on the second class of these models, considering the roughness of a surface as the result of specific stochastic processes, and based on the directed line formalism. This formalism, much used in polymer physics [25], allows to analyse the motion of a system by the way of transversal and curvilinear coordinates [26]- [28] (respectively noted y and s in this paper). And we propose models showing the analogies between the profile of a rough surface and the study of the transversal position equation. Note that, in this formalism, s plays the role of the time. For example dy/ds designs the transversal velocity. Moreover, this formalism is also linked to the statistical physics: the partition function of the system and its links with the KPZ or the Burgers equation, which follow from this formalism, represent useful tools [29]- [31]. Concerning the roughness analysis, different works have treated this problem [29], [32]- [37], but the roughness is nevertheless poor studied from a purely stochastic point of view. Among recent works treating stochastic surface properties, we can cite [38] and [39] which both have a Markovian approach, using a Langevin equation with a stochastic noise.

In the models presented here, we assumed that the profile of surface is described by the trajectory of a particle undergoing a directing stochastic strength in addition to diffusion. We will use for this purpose the elastic forcing term f (y, s) = -κy + η(s) inducing stochastic oscillations, or equivalently, a potential Φ, such as f = -∂Φ/∂y; η having for example the properties of a Langevin force. For deposition processes, we suppose that the surface start at the height h 0 , and we are interested in the height velocity variation (hvv), ∂h/∂s, which gives local information on the structure of the surface to analyse: for a fast variation of the height, we will have a high value of the hhv, (Fig. (1)). Hence, this quantity is the reflect of local irregularities. But, the hvv do not have to be confused with the quantity ∂h/∂x, where x stands for the horizontal coordinate: locally, the slope ∂h/∂x may takes an infinite value (i.e parallel to the h axis), whereas The outline of the paper will be as follow: we first derive properties related to profile of surface from pure and forced Random Walk (RW). In what follows, we will employed the term forced RW to design a RW undergoing a forcing term in addition to the diffusion. Then, we justify the use of an oscillator model comparing the directed line to a string vibrating. This will allow to use the restricted partition function to determine the transversal equation of motion, in presence of a quadratic stochastic potential. In the last part, we shall concentrate on the O-U equation of motion for the transversal position, and particularly on the asymptotic solutions for s → ∞.

(h(s 1 ) -h(s 2 ))/(s 1 -s 2 )) is weak.

II. TOPOLOGY OF SURFACE VIEW AS A RANDOM WALK PROCESS

The models presented here concern the use of the directed line formalism associated to a RW process. In its discrete form, we based our approach on the paper of Vilgis [25], treating directed polymers in random media. Our derivation is somewhat different of Vilgis, but leads to the wellknown partition function of polymer physics as showed later on.

A. Discrete random walk process

We start from a set of vectors {b i ≡ b i e s ; i = 1, . . . , N}, which are successively connected, and recovering the whole of the surface (Fig. (2)). We suppose they are statistically independent.

The probability to find such a set is given by:

P ({b i ; i = 1, . . . N}) = N i=1 p(b i ), (1)
where the probability p(b i ) corresponds to a step of the random walk of b i , with the mean value properties:

b i = 0 b i b j = b 2 δ ij . ( 2
)
δ ij is the Kronecker symbol and b is a characteristic length of the shift of position of the surface from one to the other. For the step b i , a normalized gaussian probability law is taken:

p(b i ) = (1/2πb 2 ) 1/2 exp(-b 2 i /2b 2 ), (3)
and the associated distribution of the surface length reads

P (R) = db 1 . . . db N δ R - N i=1 b i N i=1 p(b i ); (4)
the vector R = N i=1 b i characterizing the line. Then, using the Fourier transform representation of the Dirac delta, a Gaussian law for the length distribution arises quite naturally

P (R) = (1/2πb 2 N) 1/2 exp(-R 2 /2b 2 N). (5)
And from this expression, the quadratic mean "length" of the surface can be evaluated:

R 2 = R 2 P (R 2 ) d R = b 2 N. (6)
Note that the quantity R 2 preserves information on the length of the trajectory, contrary to the geometrical quantity R. However, this result is a limiting case of a pure RW trajectory since it implies that the entanglements, and so the length of the surface, could escape to infinity. Thereby, for this case we label N such as

R 2 rw = b 2 N rw , (7)
with the meaning that N rw steps are necessary in average for a pure RW to describe the profile of surface. The other limiting case, i.e. the straight line (sl), is given by

R 2 sl = b 2 N 2 sl ≃ L 2 , (8)
in such a way that the general case (g), representing a forced RW, and included between the cases ( 7) and ( 8) , may be written

R 2 g = b 2 N 2γ g ; N rw ≤ N g ≤ N sl , (9)
where γ is a critical exponent representative of a strength acting on the RW. In addition, we have the conditions

b 2 N 2γ g → b 2 N 2 sl for N g → N sl ; (10)
which implies γ = 1, and

b 2 N 2γ g → b 2 N rw for N g → N rw ; (11)
implying γ = 1/2. Thus, 1/2 < γ < 1, available for a line moving in the plane (h, x). This result may also be expressed in term of the length l ≃ R 2 g of the surface:

l = b N g γ . ( 12
)
In the case of a two-dimensional RW, γ may takes the form γ = 1/d f , where d f designs the fractal dimension of the directed line. As a matter of fact, d f describes the degree of irregularity of the surface, and for γ = 1 and γ = 1/2, d f becomes identical to the topological dimension of the system. Clearly, the exponent γ is linked to the number of steps N g . A relation between the two quantities may be obtained from the probability of the general case, P g (R), when the line undergoes both effects of diffusion and of an external force. Anticipating the continuum approach described afterwards, we write this probability

P g (R) = N exp -R 2 /2b 2 N g -βΦ(R) , (13)
where Φ(R) stands for an external potential. For example, in the case of a harmonic potential,

Φ = (κ/2)R 2 , the quadratic mean length is R 2 g = R 2 P g (R 2 ) d R = b 2 N g 1 + b 2 βκN g ( 14
)
which leads, with Eq. ( 9), to the expression

γ = ln (N g /(1 + b 2 βκN g )) 2 ln (N g ) . (15)
One can then obtain information on the topology of a surface, in the sense that the ratio l/L is the reflection of its irregularities and roughness. The fractal dimension may be evaluated experimentally, for example, by the box-counting method with boxes of length b (in the same way as for electric discharges [START_REF] Niemeyer | Scaling properties of growing zone and capacity of Laplacian fractal[END_REF]). In simulations, the varying parameter would be the number of steps of the RW.

B. Equations of evolution

If we express now the step as b i = R i -R i-1 , and writing shortly P = P ({b i ; i = 1, . . . N}), the probability (Eq. ( 1)) may be rewritten

P = 1 2πb 2 N/2 exp - 1 2b 2 N i=1 (R i -R i-1 ) 2 . ( 16
)
One can then associated to the exponential argument a symbolic Hamiltonian

βH 0 = 1 2b 2 N i=1 (R i -R i-1 ) 2 . ( 17
)
This is the energy of a mass-spring chain of a one dimensional lattice. This suggest to take the continuous limit of this discrete chain: R i -R i-1 → ∂R/∂s, and N i=1 → N 0 , where s is the dimensionless curvilinear coordinate. In this limit, the probability of the pure RW reads

P → N exp - 1 2b 2 N 0 ∂R ∂s 2 ds , (18)
where N is a normalization constant. From this result, one can construct the partition function as a sum over all the possible paths starting from s = 0 and reaching the linear size substrate at

s = N: Z = N All paths R(s) exp - 1 2b 2 N 0 ∂R ∂s 2 ds . ( 19
)
The argument of the exponential is equivalent to the transversal kinetic energy, or more precisely to the effective action of the system since we integrate over the "time" s. Furthermore, this path integral obeys the classical diffusion equation which arises naturally from a pure RW process:

∂ ∂s Z(R, s) = b 2 2 ∂ 2 ∂R 2 Z(R, s), (20)
with the boundary condition Z(R, 0) = δ(R). For a forced RW, we introduce a potential term Φ(R(s)) standing for an external (virtual) strength field, in such a way that the total effective action takes the form

β S = 1 2b 2 N 0 ∂R(s) ∂s 2 ds + N 0 Φ(R(s))ds. (21)
We note S to emphasize the fact that this quantity represents a pseudo action (the integration is made over a Hamiltonian and not over a Lagrangian). This yields to the total partition function of the system, as the Feynman-Kac path integral:

Z(R, s) = N exp -β S DR(s). (22)
It is then convenient to introduce the transversal coordinate y ≡ R/b and so to write the action in the form [26]

S = s 0 ds ′ 1 2 dy ds ′ 2 + Φ(y, s ′ ) = s 0 ds ′ H(y, s ′ ). (23)
H standing for the Hamiltonian of the system. By this way, we obtain the well-know expression of the partition function of directed polymers

Z(y, s) = N y(s) y(0) Dy exp -β s 0 ds ′ 1 2 dy ds ′ 2 + Φ(y, s ′ ) ; (24)
expression directly issued from a RW process with a directed line formalism. For a one dimensional system, Z(y, s) obeys the following PDE [START_REF] Kleinert | Path integrals in Quantum Mechanics[END_REF]:

∂ ∂s Z(y, s) = ν ∂ 2 Z(y, s) ∂y 2 -βΦ(y, s) Z(y, s), (25)
where ν = 1/2β is a viscosity coefficient.

Let us underline that there exists several versions of the Hamiltonian depending on system to study [16]- [18], [20]- [24].

In view of these results, several analogies may be drawn between directed lines and vibrating strings. In its discrete form, the directed line is composed of a succession of rigid chains and the whole evolving with a resulting motion depending on the form of the noisy potential. For the continuous form, the term of inertia in polymer physics is due to the thermal agitation of the medium, and the potential to an external field (e.g. electric field). In our case, the external potential is necessary because of the strength which have to be imposed to the line, in order to obtain realistic profiles of surface. So, the noisy potential of our model will be composed of a harmonic part and of a stochastic part which shall make the oscillations stochastic: Φ(y, s) = 1 2 ky 2 + η(y, s), where κ is the elastic constant of the line. The next step is to determine the equation of motion with the frequency of vibration of the line. In order to do this, let us proceed as follow:

In one hand, starting from Eq. ( 25), it is well-known that we may obtain successively a KPZ and a Burgers equation (approach (A)):

∂ S ∂s = ν ∂ 2 S ∂y 2 + 1 2 ∂ S ∂y 2 -Φ; (26)
and

∂u ∂s + u ∂u ∂y = ν ∂ 2 u ∂y 2 + ∂Φ ∂y . (27)
In this case the transformations are

u = - ∂ S ∂y and Z = exp S/2ν . ( 28
)
Usually one introduce the free energy. But the free energy tends to the total energy of the system S when the entropy tends to zero.

In the other hand, starting from a general Burgers equation for a velocity field u, we have (approach (B))

∂u ∂s + u ∂u ∂y = ν ∂ 2 u ∂y 2 - ∂V ∂y , (29)
and

∂F ∂s = ν ∂ 2 F ∂y 2 + 1 2 ∂F ∂y 2 + V, (30)
and finally, with Z = exp (F/2ν),

∂ ∂s Z(x, s) = ν∇ 2 Z(x, s) + V 2ν . (31)
Whether u represents the transversal velocity, we have in this case

u = dy ds = - ∂F ∂y . (32)
From these two approaches, a question arises: "May the quantity u of Eq. ( 27) be interpreted as the transversal velocity of Eq. ( 29) ?" In other words, is there identity between S and F ?

To remove this ambiguity, let us start from the Burgers equation ( 29), u given by (32)

The first and last integrals may be rewritten

In such a way that the "free energy" will read

F (y, s) = -ν - s 0 1 2 dy ds ′ 2 -V (y, s ′ ) ds ′ . (38)
Finally, the Hopf-Cole transformation Z = exp (F/2ν), with 1/2ν = β, yields explicitly to

Z = e -1/2 × exp -β s 0 ds ′ 1 2 dy ds ′ 2 -V (y, s ′ ) . ( 39
)
The heat equation ( 25) being invariant under the shift Z → Z = c × Z, c ∈ R, we wiil have equivalence between the two approaches if Φ = -V . And, therefore, one obtain the classical Hamilton equation of motion:

dy ds = - ∂F ∂y = - ∂ S ∂y ⇐⇒ d 2 y ds 2 = - ∂H ∂y (40)
The determination of an equation of motion for y is henceforth equivalent to solve the KPZ equation (26). Moreover, we have a relation between the height of the growing surface and the transversal position: |h(s)h 0 | ≡ |y(s)y(0)|, and so, by extension,

dy ds ≃ ∂h ∂s . ( 41
)
For example ∂h/∂s = 0 ⇔ dy/ds = 0. Thereby, the action may as well be written [START_REF] Sarma | [END_REF]43]

S = s 0 1 2 ∂h ∂s ′ 2 + Φ(h, s ′ ) ds ′ , (42)
We shall now concentrate on the solution of the KPZ equation (26).

Analytic solutions of noisy KPZ equations have already been studied in the literature in specific cases [44]- [46]. But the non-linearity of this equation added to the presence of a stochastic noise make the resolution difficult. Nevertheless, it is known that this equation may be transformed either into a Burgers equation or a Schrödinger equation with an imaginary time. And, among works which treat exact solutions of theses kind of equations, we can cite [47]- [51]. This is our concern here to treat the case of a line undergoing the stochastic elastic forcing term -κy + η(s),

where the function η contains all the noise that affects the oscillations (e.g. a Langevin force). This amounts to saying that we solve the KPZ equation ( 26) with a potential Φ given by Φ(y, s) = κ 2 y 2yη(s). Moreover, we assume a constant initial condition for the transversal velocity: u(y, 0) = u 0 .

Following methods presented in Ref. [51], it is shown in appendix that the solution is of the form:

S(y, s) = ξ(s) 2 y 2 -yψ(s) + Θ(s), (43)
where ξ(s), ψ(s) and Θ(s) are functions of s and of η(s) (see the appendix). But, we have the

relation dy ds = - ∂ S ∂y = -yξ(s) + ψ(s). (44)
Consequently, we obtain an equation standing for the equation of a forced oscillator with a term of friction:

d 2 y ds 2 + ξ(s) dy ds + ω 2 (s) y = F (s), (45)
where the pulsation ω 2 and the forcing term F are respectively given by dξ/ds and dψ/ds. We can notice from the expressions of ξ and ψ that, while the frequency f ∝ dξ/ds depends only on s, the term of force depends also on the stochastic noise η. Thus, we expect profiles of the form drawn in Fig.

Once y determined, we can deduce the hvv from Eq. ( 44):

∂h ∂s ≃ -ξ(s) exp - s 0 ξ(t)dt 1 + s 0 ψ(t) exp t 0 ξ(t ′ )dt ′ dt + ψ(s). (46)
Then, the degree of irregularity of a surface, and so its global roughness, may be estimated by adding up all the local slopes |∂h/∂s|, in such a way that the interesting quantity to determine is |∂h/∂s| , where the mean is taken over all the length of the surface.

Note that this problem is different of the one of a stochastic oscillator, which is a specific problem [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Poulin | [END_REF]

III. LONG SIZE SUBSTRATE

We treat now the case of a particle of mass unity, moving along a rough surface when we study its profile at a fixed time.

Let us consider the propagation of a Brownian particle through virtual scatterers. Its trajectory describing the surface to analyse. We assume that the particle obeys a general Langevin equation where an external force -∂Φ/∂y = -κy is considered:

d 2 y ds 2 = -ξ dy ds -κy + Γ(s). (47)
This means that the particle undergoes effects of the medium (term of friction) in addition with an oscillator term. Then, the Langevin force Γ(s) may be expressed in term of a Gaussian with noise η(s), such as

d 2 y ds 2 = -ξ dy ds -κ y + √ 2Dη(s); (48)
where D = ξ 2 ν. By nature, the solution of this equation follows a Markovian process. When the inertial term is negligible, one obtain an Ornstein-Uhlenbeck process for the transversal position y:

dy ds = - κ ξ y + √ 2ν η(s). (49)
This well-known process [54,55] is associated to the following Fokker-planck equation for the transition probability P (y, s|y ′ , s ′ ):

∂ P ∂s = κ ξ ∂y P ∂y + ν ∂ 2 P ∂y 2 . ( 50
)
With the standard initial condition P (y, s|y 0 , s 0 ) = δ(yy ′ ), the solution writes [START_REF] Risken | The Fokker-Planck equation[END_REF] P (y, y ′ , s) = κ 2πξ ν (1e -(2κ/ξ)s ) exp

  - κ y -e -(2κ/ξ)s y ′ 2 2ξ ν 1 -e -(2κ/ξ)s   . ( 51
)
This probability is asymptotically "stationary":

lim s→∞ P (y, y ′ , s) = κ 2πξ ν exp - κy 2 2ξ ν . ( 52
)
In practice, the condition s ≫ ξ/2κ is sufficient to obtain this solution. This implies that, under this condition, the fluctuations of the transversal position y becomes independent of the curvilinear abscissa. Moreover, the quantity Λ = ξ ν/κ, appearing in Eq. ( 52) with the dimension of a length, takes the sense of the mean range for the transversal fluctuations with respect to a given position:

+∞ -∞ y 2 P (y)dy = Λ 2 . ( 53
)
And the Einstein relation ν = 1/ξβ gives us the expression

Λ ∝ κ-1/2 (54)
Thereby, the higher the elastic constant, the lower the range of fluctuations. The connection between the O-U process and the roughness is then highlighted in Fig. (4), where we can notice that, for s ≫ 1, the maximal fluctuation values of the height h(y, s) are of the order of Λ. As a consequence, the hvv will be such as

    ∂h ∂s     ≤ Λ ∀s ≫ ξ/2κ. (55)
Consequently, in this limit,

    ∂h ∂s     ≤ Λ; (56)
available for models concerning processes of deposition executed from the left to the right of the figure, and on a "large" length of substrate.

Then, an interesting point lies in a connexion between the O-U process and KPZ equations.

Indeed, the transformation P (y, s) = Z(y, s)e -κy 2 /4ν , (57) applied to Eq. ( 50), leads to the heat equation and so to the KPZ equation

∂ Z ∂s = ν ∂ 2 Z ∂y 2 + κ 2 - κ2 y 2 4ν Z, (58)
∂ S ∂s = ν ∂ 2 S ∂y 2 + 1 2 ∂ S ∂y 2 + κ 2 - κ2 y 2 4ν . (59)
And once more, the introduction of a quantity u such as u = -∂ S/∂y, yields to a Burgers equation:

∂u ∂s + u ∂u ∂y = ν ∂ 2 u ∂y 2 - ∂ Φ ∂y . (60)
But the same interrogation as the one of the RW process may be done, namely, may u be considered as the velocity dy/ds of Eq. ( 49) ? As a matter of fact, these relations have been obtained with the assumption d 2 y/ds 2 ≃ 0. So, let us assume that u = dy/ds. The consequence is that the first term of Eq. (60) vanishes, and so we obtain

1 2 ∂ ũ2 ∂y = ν ∂ 2 ũ ∂y 2 - ∂ Φ ∂y ⇔ ∂ ∂y 1
Thus, u = dy/ds ⇔ u obeys to Eq. ( 62). Then, we refound Eq. ( 49) by putting c 1 = -κ/ξ and

c 2 (s) = √ 2νη(s).
Once again, a stochastic process (here Markovian) may be connected to the coupled KPZ/Burgers equations, and by extension to the Hamiltonian of the system. In addition, a partition function like the one of the RW may be obtained with appropriated ν and κ. Moreover, this

shown the relevant choice of a harmonic potential ∼ y 2 to treat these kind of fluctuations, since a Markovian process and a RW process yields both to the same class of equations.

IV. CONCLUSION

Two short stochastic models of directed line have been presented and applied to the description of rough surfaces. The background of our approach was that roughness of a one-dimensional surface is created by the motion of a virtual particle dived in a stochastic harmonic potential. It Although we have presented one-dimensional models, these approaches could open interesting ways of investigation for numerical and analytical studies, in particular if one needs to generate surfaces with specific profiles or properties.

FIG. 1 :

 1 FIG. 1: Profile of surface: result of the conflict between the random medium and a directing strength -∂Φ/∂y. The local slope (a) (resp. (b)) stands for a weak (resp. high) value of the hvv.

FIG. 2 :

 2 FIG. 2: Random Walk: N steps of length b are required in order to cover the surface on a linear size substrate L.

FIG. 3 :

 3 FIG. 3: Variations of the transversal position y for a forcing term inducing oscillations depending on the intensity of the noise. Each mode (a), (b), and (c) leads to very different roughness.

FIG. 4 :

 4 FIG. 4: Asymptotic profile of surface for an O-U process. Λ has to be interpreted as the maximal range of fluctuations.

  have that -∂H/∂y = d 2 y/ds 2 ≃ 0. So the last part of Eq. (61), leads to a solution of the form dy ds = c 1 (s)y + c 2 (s).

  brings out of the continuous RW that we can rigorously linked the partition function of the line to the Hamilton equation of motion for the transversal position. The survey of roughness amounts therefore to solve a KPZ equation with a particular potential. Then, interpreting the roughness as a result of a stochastic elastic forcing term acting on a surface, we have obtained a forced oscillator equation where the force behaves as a noise for the line. The analytical expressions of the pulsation and of the force have been determined too. Consequently, local analysis of the roughness may be performed thanks to the expression of the hvv (46) which, besides, obeys a Burgers equation.For a sufficient long substrate, the O-U model shows that the maximal values of roughness are asymptotically bounded in a range of length Λ ∝ κ-1/2 . The elastic constant being an adjustable parameter mastering the magnitude of the fluctuations. Moreover, we have shown that the two models, used with a same kind of force, had common properties since they lead both to the same class of equations. As a matter of fact, the two equations of motion are connected to the partition function of the system. Finally, these two models belong to stochastic processes described by the Burgers dynamics.

APPENDIX: SOLUTION OF THE KPZ EQUATION

We solve here the KPZ equation (26) or, equivalently, the heat equation (25) written in the following way:

with the initial condition Z(y, 0) = exp (-u 0 y/2ν). The resolution is based on the Time Space Transformation method (TST) of [51]. We give only the main points of resolution. We would like to underline for the comprehension that the notations x and y are reversed between the present paper and [51]:

The first step is to determine the initial condition. For this, we need to determine the expression for the variables a 1 , a 2 (s), p(s), q(s) and τ (s) of [51]. We obtain successively:

In such a way that the resolution of (A.1) amounts to solve the equation:

with the initial condition

The solution reads explicitly

Then, reversing the different transformations, we obtain as a solution of (A.1)

expression that we write for further analyse

where

√ κ e 2 √ κs tanh( √ κs) , (A.12)

The quantity a 1 , a 2 (s) and q(s) being respectively given by Eqs. (A.2), (A.3) and (A.5). Note that we recover the initial condition Z(y, 0). Then, since Z = exp S/2ν , we can write the solution on the final form S(y, s) = ξ(s) 2 y 2yψ(s) + Θ(s) (A.15)


Mickael Laumonier
email: m.laumoniermickael.laumonier@gmail.com

Bruno Scaillet

Laurent Arbaret

Joan Andújar

Rémi Champallier

Experimental mixing of hydrous magmas

Keywords: Magma mixing, mingling, hydrous, shearing, deformation, texture, enclave

Deformation experiments involving hydrous magmas of different compositions (basalt and haplotonalite) have been performed in a Paterson press at 300 MPa, in the temperature range 600°C-1020°C, with water-saturated melts, during 2-4 hours. Prior to deformation the two end-member magmas were annealed at either 950°C or 1000°C, yielding magmas with crystal contents in the range 31-53 wt% and 2 sets of viscosity contrasts. Under the experimental conditions investigated (i.e. moderate shear rates <10 -3 s -1 ), mixing/mingling textures appear at temperatures > 950°C. In the temperature range 950-985°C a few mixing and mingling textures occur, though both end-members essentially retain their physical integrity. It is only at, or above, 1000°C that a dramatic jump in mingling efficiency happens, corresponding to a crystal fraction of 45 vol%. Textures include entrainment of mafic crystals into the felsic magma, mafic-felsic banding, enclave formation, diffusion-induced interface, the latter only

A C C E P T E D

M A N U S C R I P T

INTRODUCTION

Magma mixing is commonly observed in nature, in particular in arc magmas, which are characterised by an abundance in volatiles, notably water (Anderson, 1976;Armienti et al., 1983;[START_REF] Blake | Magma mixing and hybridization processes at the alkalic, silicic, Torfajokull central volcano triggered by tholeiitic Veidivotn fissuring, south Iceland[END_REF][START_REF] Castro | Magma mixing in the subvolcanic environment: petrology of the Gerena interaction zone near Seville, Spain[END_REF][START_REF] Coombs | Experimental and textural constraints in mafic enclave formation in volcanic rocks[END_REF][START_REF] De Rosa | Relationships between deformation and mixing processes in lava flows: a case study from Salina (Aeolian Islands, Tyrrhenian Sea)[END_REF][START_REF] Druitt | Santorini Volcano[END_REF]Martin et al., 2006a,b;[START_REF] Pons | Development of Igneous layering during growth of pluton: the Tarçouate Laccolith (Morocco)[END_REF][START_REF] Pal | Dacite-andesite of Narcondam volcano in the Andaman sea -An imprint of magma mixing in the inner arc of the Andaman-Java subduction system[END_REF][START_REF] Woods | Magma mixing triggered during volcanic eruptions[END_REF]Davì et al., 2010;[START_REF] Eichelberger | Messy magma mixtures[END_REF][START_REF] Perugini | The mixing of magmas in plutonic and volcanic environment:Anologies and differences[END_REF]. The effects of mixing are particularly evident when evolved magma chambers stored at upper crustal pressures (P < 400 MPa) are replenished by batches of hot (up to 1250°C) and mafic magmas (e.g. [START_REF] Sparks | Magma mixing: A mechanism for triggering acid explosive eruptions[END_REF][START_REF] A C C E P T E D M A N U S C R I P T Accepted Manuscript | Evidence of magma mixing: petrological study of Shirouma-Oike calc-alkaline andesite volcano, Japan[END_REF]1981;Bacon, 1986;[START_REF] Civetta | Magma mixing and convective compositional layering within the Vesuvius magma chamber[END_REF][START_REF] Nakamura | Continuous mixing of crystal mush and replenished magma in the ongoing Unzen eruption[END_REF][START_REF] Pallister | Magma mixing at Mount Pinatubo: Petrographic and chemical evidence from 1991 deposits[END_REF][START_REF] Mandeville | Magma mixing, fractional crystallization and volatile degassing during the 1883 eruption of Krakatau volcano, Indonesia[END_REF][START_REF] Wiebe | Mafic-silicic layered intrusions: The role of basaltic injections on magmatic processes and the evolution of silicic chambers[END_REF][START_REF] Venezky | Preeruption conditions and timing of dacite-andesite magma mixing in the 2.2 ka eruption at Mount Rainier[END_REF][START_REF] Clynne | A Complex Magma Mixing Origin for Rocks Erupted in 1915, Lassen Peak, California[END_REF][START_REF] Miller | Mount Dutton Volcano, Alaska: Aleutian Arc analog to Unzen volcano, Japan[END_REF][START_REF] Browne | Magma mingling as indicated by texture and Sr / Ba ratios of plagioclase phenocrysts from Unzen volcano, SW Japan[END_REF][START_REF] Pal | Dacite-andesite of Narcondam volcano in the Andaman sea -An imprint of magma mixing in the inner arc of the Andaman-Java subduction system[END_REF][START_REF] Eichelberger | Messy magma mixtures[END_REF][START_REF] Ruprecht | Pre-eruptive reheating during magma mixing at Quizapu volcano and the implications for the explosiveness of silicic arc volcano[END_REF]. Magma mixing has an important role in the dynamics of the processes occurring in magmatic reservoirs and it has been proposed as a triggering mechanism of volcanic eruptions (e.g., [START_REF] Sparks | Magma mixing: A mechanism for triggering acid explosive eruptions[END_REF][START_REF] Pallister | Magma mixing at Mount Pinatubo: Petrographic and chemical evidence from 1991 deposits[END_REF][START_REF] Eichelberger | Messy magma mixtures[END_REF][START_REF] Kent | Preferential eruption of andesitic magmas through recharge filtering[END_REF][START_REF] Ruprecht | Pre-eruptive reheating during magma mixing at Quizapu volcano and the implications for the explosiveness of silicic arc volcano[END_REF][START_REF] Felice | The 2009 paroxysmal explosions at Strombolie (Italy): magma mixing and eruption dynamics[END_REF][START_REF] Druitt | Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano[END_REF]. Mixing depends on the rheological properties of magmas that are in turn strongly affected by volatiles. In particular, dissolved water strongly decreases melt viscosity [START_REF] Dingwell | The effect of water on the viscosity of a haplogranitic melt under P-T-X conditions relevant to silicic volcanism[END_REF], affects the fraction of crystals present at a given temperature, which, in turn, affects magma viscosity [START_REF] Champallier | Experimental investigation of magma rheology at 300MPa: From pure hydrous melt to 76 vol.% of crystals[END_REF][START_REF] Caricchi | Non-Newtonian rheology of crystal-bearing magmas and implications for magma ascent dynamics[END_REF][START_REF] Picard | Rheology and microstructures of experimentally plagioclase suspensions[END_REF].

Consequently water and volatiles, in general, strongly influence magma behavior and eruptive style, depending on whether the magma holds or loses its volatiles during its transfer to surface (e.g., [START_REF] Jaupart | Gas content, eruption rate and instabilities of eruption regime in silicic volcanoes[END_REF][START_REF] Martel | Magma storage conditions and control of eruption regime in silicic volcanoes: experimental evidence from Mt[END_REF][START_REF] Laumonier | Porosity redistribution enhanced by strain localization in crystal-rich magmas[END_REF]. However, dynamic experiments designed to explore magma mixing processes at the water-rich conditions found in arc settings in particular are absent [START_REF] Kouchi | Experimental study of mixing of basaltic and dacitic magmas[END_REF][START_REF] Kouchi | A model for mixing basaltic and dacitic magmas as deduced from experimental data[END_REF][START_REF] Watson | Behavior of alkalies during diffusive interaction of granitic xenoliths with basaltic magma[END_REF][START_REF] Wyllie | Interactions among magmas and rocks in subduction zone regions: experimental studies from slab to mantle to crust[END_REF]Carroll & Wyllie, 1989;[START_REF] Van Der Laan | Experimental interaction of granitic and basaltic magmas and implications for mafic enclaves[END_REF][START_REF] De Campos | Heterogeneities in magma chambers: Insights from the behavior of major and minor elements during mixing experiments with natural alkaline melts[END_REF]2011). So far, most works have been conducted at

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

atmospheric pressure (i.e. using dry magmas), high temperature (1200 to 1400°C), under static conditions or at high strain rates (~10 2 s -1 by [START_REF] Kouchi | Experimental study of mixing of basaltic and dacitic magmas[END_REF][START_REF] Kouchi | A model for mixing basaltic and dacitic magmas as deduced from experimental data[END_REF]; ~10 -1 s -1 by De [START_REF] De Campos | Heterogeneities in magma chambers: Insights from the behavior of major and minor elements during mixing experiments with natural alkaline melts[END_REF]2011; and ~10 -2 s -1 by Morgavi et al., 2013a;[START_REF] Morgavi | Time evolution of chemical exchanges during mixing of rhyolitic and basaltic melts[END_REF]. Such conditions, in particular strain rate, are consistent with volcanic eruption processes, but exceed those characteristic of magma reservoirs prior to eruption (10 -4 s -1 , Chadwick et al., 1988; Albertz et al., 2005; Hodge et al., 2012). Recently, Laumonier et al. (2014a) deformed two chemically distinct and dry magmas at conditions close to expected magmatic ones, i.e. 900 < T < 1200°C, 300 MPa and strain rate of 10 -5 to 10 -3 s -1 . The geometry used imposes simple shear at the interface between the two end-members, to simulate replenishment of a reservoir by a mafic dyke intrusion (Fig. 1A). These dry experiments reproduce most of the textures of magma mixing and mingling observed in rocks (enclaves, stretched filament, isolated crystal in disequilibrium with their host), clarifying the details of the mechanisms occurring during the incipient stages of mixing. The fraction and the arrangement of crystals appears to be a critical factor for magma mixing, in particular via the existence of crystal network, which controls the magma rheology (e.g., [START_REF] Philpotts | Role of plagioclase crytal chains in the differentiation of partly crystallized basaltic magma[END_REF]Martin et al., 2006a;Laumonier et al., 2014a;[START_REF] Caricchi | Insights into the mechanisms and timescales of pluton assembly from deformation patterns of mafic enclaves[END_REF]. However, the dry conditions explored by Laumonier et al. were aimed at limiting quench effects, to better document mixing textures and related mechanisms. Here, we extend our efforts toward hydrous conditions using a methodology similar to that in Laumonier et al. (2014a). To the best of our knowledge this is the first time that the effect of water on the mixing capacity of magmas is experimentally investigated under conditions of temperature, pressure, and strain rate relevant to subduction zone settings.

EXPERIMENTAL PROCEDURE

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

2.1.

Starting material

A natural basalt and a synthetic haplotonalite were selected as the mafic and felsic end-members for torsion experiments, allowing a direct comparison with previous dry experiments (Laumonier et al., 2014a). The basalt (composition in Table 1) was sampled from the massive part of the Cape Balos flow on Santorini volcano, Greece [START_REF] Nicholls | Petrology of Santorini Volcano, Cyclades, Greece[END_REF][START_REF] Druitt | Santorini Volcano[END_REF]. It contains phenocrysts of olivine (Fo75) that are set in a matrix of plagioclase, clinopyroxene, magnetite with rare ilmenite and orthopyroxene. Xenocrysts of olivine (Fo 78 )

and plagioclase (An 90 ) are also present, the latter sometimes displaying core sieved textures and inverse compositional zoning (cores: An55-60, rims An80) [START_REF] Nicholls | Petrology of Santorini Volcano, Cyclades, Greece[END_REF]Andújar et al., 2015). The phase equilibria of the mafic material have been investigated by Andujar et al.

( 2015) and are used as a guideline for our work. The haplotonalitic glass (Table 1) was produced by Schott A.G. (Germany). In the P-T-H 2 O range explored it crystallizes plagioclase only and the relationships between crystal fraction, water content, and temperature have been previously determined [START_REF] Picard | HP-HT deformation of silicic magmas: experimental constrains on the structural evolution and the rheological threshold at middle and high crystallinities[END_REF][START_REF] Picard | Rheology and microstructures of experimentally plagioclase suspensions[END_REF][START_REF] Laumonier | Porosity redistribution enhanced by strain localization in crystal-rich magmas[END_REF]. The chemical and rheological behaviors of such a plagioclase suspension in a felsic glass are similar to natural felsic magmas (trachyte to rhyolite, e.g. [START_REF] Calanchi | Silicic magma entering a basaltic magma chamber: eruptive dynamics and magma mixingan example from Salina (Aeolian islands, Southern Tyrrhenian Sea)[END_REF][START_REF] Ferla | Evidence of magma mixing in the 'Daly Gap' of alkaline suites: a case study from the enclaves of Pantelleria (Italy)[END_REF][START_REF] Davi | Mafic enclaves in the rhyolitic products of Lipari historical eruptions; relationships with the coeval Vulcano magmas (Aeolian Islands, Italy)[END_REF].

These two starting materials were first hydrated and annealed in an Internally Heat Pressure Vessel (IHPV) at 300 MPa and 950 or 1000°C, depending on the target crystal fraction. The felsic magmas synthesized at 950°C and 1000°C have crystal contents (Φs) of 38 and 31 vol.% respectively, with ~10 µm long plagioclase (An29 to An35) having an aspect ratio of 2-3, along with a few percent of bubbles (<2% in volume) and glass (Fig. 2A; Table 1). The mafic syntheses contain 53 vol% (at 950°C) and 45 vol % (at 1000°C) of amphibole + plagioclase + pyroxene + magnetite and glass (Fig. 2C & D, Table 1). More details about the preparation and syntheses can be found in the supplementary information.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

2.2.

Deformation experiments

2.2.1.

Experimental set up for deformation experiments

Cylinders from first step synthesis products were drilled out (13.78 to 14.96 mm diameter) and cut into thin disks (1.16 to 3.34 mm thick) to build the experimental torsion assembly while some pieces from different locations in the syntheses were selected to check the suspension homogeneity. Our deformation set up consists of 4 interleaved wafers of felsic/mafic magma synthesized at the same temperature and alternating in composition, always with a felsic disk located atop of the "stack" (Fig. 1B). Layers are named according to their position and composition, e.g. numbering from the top, and ρ and β symbols for felsic and mafic compositions, respectively (Fig. 1B). In all experiments except the run conducted at the lowest temperature (3 layers: ρ1, β2 and ρ3), the upper layer was ρ1 and the lower one is

ACCEPTED MANUSCRIPT

reservoir by a mafic dyke (Fig. 1; e.g. [START_REF] Castro | On granitoid emplacement and related structures. A review[END_REF]. During preliminary tests, we used a different sample geometry consisting of two half cylinders, one mafic and one felsic, sharing a vertical interface. However such a geometry was not mechanically sustainable, the rapid collapse of the less viscous material imparting an off-axis rotation of the entire column which could lead to severe damages of the apparatus: hence it was not used any further.

Experimental conditions

Seven torsion, and one static experiments were performed at 300 MPa, in the temperature interval 600-1020°C (Table 2). The prevailing oxygen fugacity was measured with a Shaw membrane during a calibration experiment and found to be 1 log unit above the Nickel-Nickel oxide solid buffer (fO 2 ≈ NNO+1±0.5) [START_REF] Scaillet | Improvements of the Shaw membrane technique for measurement and control of fH2 at high temperatures and pressures[END_REF]. Torsion was applied during 1 to 4 hours at strain rates 2.10 -4 <  < 9.10 -4 s -1 allowing to achieve relatively large finite strains (0.7 < γ < 5.1). Most torsion experiments were performed in the temperature range 950°C-1020°C, with an iron jacket (Table 2). Torsion experiments started by the pressurization of the vessel up to ~220 MPa, followed by heating which permits to reach the final pressure (300 MPa). Mixing experiments were performed using starting materials whose synthesis temperature was the closest to that of the deformation experiment (Table 2). The temperature difference relative to that of the synthesis may induce melting (PP262, PP293 & PP295, Table 2) or crystallization (PP258, PP265 & PP296), during the 1hour lap time allowed before torsion was applied. The lap time of low temperature experiments (600°C and 715°C) lasted 15 minutes so as to minimize crystallization. Torsion was applied until jacket failure was noticed by the drop of the load applied on the column holding the sample and experiments were ended by an isobaric quench (~60°C/mn).

Additional details on the experimental set-up and techniques can be found in [START_REF] Paterson | Rock deformation tests to large shear stress in torsion[END_REF].

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

After experiments, samples were recovered and embedded in epoxy resin for textural observation and chemical microanalyses. Generally, the initial geometry of each layer is well preserved, no significant compaction or vertical shrinkage was noticed (see lengths before and after experiments, determined from SEM pictures, Table 2). Observations were made on tangential sections far enough from the sample periphery to avoid interferences resulting from the jacket deformation (see details in supplementary information). The textural and chemical results of the four experiments characterizing the transition from non-mixing to mixing conditions are detailed (950-1000°C; Table 2).

RESULTS

3.1.

Strain distribution and macroscopic observations

During torsion experiments the jacket that moulds the sample is also affected by the deformation [START_REF] Paterson | Rock deformation tests to large shear stress in torsion[END_REF]. Thus, the local strain distribution along sample length (affecting each layer) can be directly inferred by photo reconstruction and image analysis of the outer jacket (Fig. 3). Whenever the strain is not homogeneously distributed across the stack, the deflection of vertical linear markers on the jacket, together with standard trigonometric relationships, are used to determine the local strain distribution (Fig. 3). The quantification of the local strain was confirmed by calculations based on the twisting angle and the sample dimension [START_REF] Paterson | Rock deformation tests to large shear stress in torsion[END_REF].

While experiments conducted at subsolidus temperature exhibit a strong partitioning of the deformation, the ones above 950°C display deformational textures generally affecting both layers (Fig. 3). In detail, SEM pictures show that the felsic material, also being less viscous, still preferentially accommodated the strain until 985°C (Fig. 4A-C). In contrast, at T > 1000°C, abundant mingling features such as irregular interfaces between end-members, boudinage and mafic enclaves (in the sense of aggregates of mafic minerals) in the felsic

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

magma were produced, despite a relatively low strain γ bulk < 1.3 (Fig. 4D &E). In those experiments conducted at T > 1000°C, both end-members were clearly deformed, in agreement with the near homogeneous strain distribution observed on the jacket (Fig. 3D).

Lastly, two of the seven experiments (PP295 and PP261) were affected by extrusions of material through holes in the jacket (Fig. 3B &D). Although such extrusions were uncontrolled they shed light on mingling and mixing features produced at higher strain (supplementary information).

3.2.

Microscopical observations

3.2.1.

Phase proportions

The crystal content was determined by SEM-image analysis of areas in mafic layers that were close to the felsic component and having a surface > 0.8 mm² to avoid local heterogeneities.

Textural analyses concerning the PP261 experiment were performed over a smaller area (0.35 mm²) for the enclave and farther from the contact because of the interface complexity (see below). While the experiments conducted at low temperature (600 and 715°C) have a phase assemblage and proportions similar to the starting material (Fig. 2C & 5A; Table 1), felsic and mafic layers in experiments conducted at T>950°C display considerable textural changes relative to starting products, in particular a higher crystal content in both magmas, as a result of quench crystallization. For instance, the melt fraction in mafic layers decreased from 48-58 to 24-32 vol.% (Fig. 5A). The experiment conducted at 985°C (grey squares on Figure 5) has the lowest melt fraction as a result of a large increase in plagioclase proportion, more than 4

times the starting values (Fig. 5B). Olivine occurs in all experiments at 950°C or higher temperature (up to 3% in volume). Oxides are generally absent, but occasionally observed (e.g. PP261, PP285). In felsic layers the crystal fraction is similar to that of the corresponding synthesis at low temperature (715°C), but increases in all run products annealed at higher

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
temperatures, up to 65 vol.% in the products of PP296. No dependence of the melt/crystal fraction with temperature, duration or finite strain of the experiment could be established.

Changes in crystal fraction in the felsic layers also result from diffusion at the interface, more or less affected by deformation.

Textures

In the run products of PP285 (950°C), PP295 (975°C) and PP296 (985°C), interfaces are no more straight but wavy (dashed line labeled 1, Fig. 6A to C). Despite their regular and nearly cylindrical shape, mafic layers underwent little deformation. A few single plagioclase and pyroxene crystals coming from the mafic layers are found in the felsic one (yellow arrow on Figure 6A). The SPO intensity defined by plagioclase in the felsic layers is low (<1.02) except in the experiment PP296 that suffered a high bulk strain. Crystal fractions of synthesis experiments are no longer preserved: the plagioclase population in the felsic layers is denser close to the mafic layer creating a halo 50 µm (PP285) to 270 µm (PP296) thick (dashed line labeled 2, Fig. 6), with crystal fractions up to 65%, whereas areas farther from mafic layers contain 47% of plagioclase.

Despite the relatively low amount of deformation reached in the experiments at 1000 and 1020°C (respectively =1.3 and =0.7), interfaces between end-members are also characterized by embayments and lobes at the scale of tens of µm (Fig. 6D). Wavy contacts are accompanied by isolated crystals and parcels of mafic magma, hereafter called "mafic enclaves", the latter having rounded shape and sizes ranging from tens of µm (aggregates constituted of a few crystals) up to 1 mm (Fig. 6D &E). Enclaves have a texture similar to the inner mafic layer one. Around the enclave presented in Figure 6E, a light SPO developed in the felsic layer made of 49% plagioclase, with intensities of 1.03 and 1.04, coinciding with the direction of the interface between the enclave and its host (Fig. 6E). The poorly developed

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

SPO of plagioclases around the enclave in the felsic layer supports the localization of the deformation instead of the development of a pervasive fabric.

3.3.

Compositions of phases

Unfortunately, the extensive quench crystallization affecting the run products severely complicates the interpretation of compositional data, in particular melt, hence obscuring their possible implications for natural contexts. For this reason we limited data acquisition essentially to the analyses of crystals in run products obtained at 1000°C considering the following domains or textures:

-">500": inner parts of mafic layers, relatively far (500 µm minimum) from the contact with the felsic component. Such parts are likely to have had little to no interactions with the felsic component;

-"<150": external parts of the mafic layers, at 150 µm maximum from the contact;

-"interface": such as the layer bridging the boudins (Fig. 6D) or that close to the lobe (Fig. 6D,F);

-"enclave": the enclave (in the sense of rounded crystals aggregate) produced in the PP261 experiment (Fig. 5E);

-"isolated": crystals located in the felsic layer with a composition and a texture different from the typical plagioclase population crystallizing from the haplotonalitic liquid.

-"ρ>300": concerns only the plagioclases in the felsic end-member at a distance >300 µm from any mafic layer and may be compared with the composition of plagioclase from the felsic synthesis (Htn). The phases from torsion experiments are compared with those from the 1000°C synthesis and the static experiment also conducted at 1000°C.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Plagioclase

Plagioclase cores are similar to the plagioclase in the starting suspension (An80-An92)

whereas rims are more sodic with An50-An76 (Fig. 7A). Plagioclases produced in the static experiment systematically developed a rim wherever they are located in the mafic layer. In general, plagioclases from torsion experiments are similar, with a large rim (An39-An77) to core (An79-An93) zonation. Analyzed plagioclases located within the mafic layers (>500)

have a rim that extends to An44, slightly less calcic than those from the mafic inner parts (>500) of the static experiment (An50). In contrast, rims of crystals close to felsic component (<150µm; An37-An77) and non zoned acicular crystals located in the interface layer (An40-An67, Fig. 7A) have compositions less calcic than those in the residual mafic melt. Such acicular plagioclases are larger when located close to the mafic layer. The smallest ones occur close to the felsic component and could not be analyzed because of their size.

Some plagioclases having similar size and aspect ratio to those from the mafic layer do also occur in the interface layer, with an anorthite content typical of the core composition (An82-An93). Theses plagioclases are frequently surrounded by acicular crystals, whose occurrence makes difficult the analysis of the thin rim overgrowth. Isolated large plagioclases in the haplotonalite, presumably from the mafic component, are also zoned with a calcic core composition (An79-An90) and a sodic rim (An40). Plagioclases with a typical felsic composition (~An35) were only observed within the felsic layer (ρ>300).

Olivine

Olivine, which is absent from the starting syntheses performed in conventional IHPV, is present in both static and torsion experiments conducted in the Paterson Press (Table 1). In the static experiment, olivine has a Fo content ranging between Fo70 and Fo77 (Fig. 7B). The olivines of torsion experiments display a larger Fo content variation, especially in the mafic

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

areas close to the felsic component (Fo64-Fo81 in <150 areas). Olivines within the interface layer, from enclaves or isolated, display a more restricted variability in Fo (Fo75-80), though Fo is not higher than in <150 areas, from where they might be sourced (Fig. 7B).

Pyroxene

The composition of pyroxene is shown in Fig. 7C. Pyroxenes from the static experiment have a composition similar to those from the starting mafic suspension when they are located within the mafic layer (>500), whereas close to the haplotonalite (<150) their composition varies along a trend with constant Ferrosilite ≈ 12%, and varying Ca and Mg contents (Fig. 7C) though changes in En and Wo contents remain small. In torsion experiments, most pyroxenes have a composition quite similar to the starting pyroxenes (Fig. 7C). Pyroxenes from the interface layers, enclaves, the inner (>500) and the outer (<150) parts of the mafic layers, also roughly follow the Fs = 0.12 trend (Ca concentration ranging between 32 and 47 molar %).

Glass

Glass analyses of static and torsion experiments reveal SiO 2 contents ranging between 55 and 69 wt% (Fig. 8), falling between the SiO 2 content of the starting glasses, though some analyses are more mafic. In contrast, the CaO concentration ranges between 2.6 and 8.2 wt.%, which plot slightly above the trend between both end-members (Fig. 8). Most melt analyses located in the external parts of the mafic layers (<150) and within enclaves have compositions slightly more evolved than the starting mafic one (higher SiO 2 and lower CaO content). In contrast, melt pockets between acicular crystals of the interface layer display heterogeneous melt compositions (SiO 2 ranging from 56 to 68 wt.%; Fig. 8).

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

DISCUSSION

4.1.

General considerations

The textures produced in the torsion experiments record, in part, the incipient stages of mafic-felsic interaction processes when the shearing deformation is parallel to the interface at relatively low strain rates, once the thermal equilibrium has been reached. As already stated, this scenario reproduces the local deformation pattern happening during the replenishment of a reservoir when a magma rises through dykes in a reservoir, though experimental strain rates might also be applicable to some volcanic processes [START_REF] Castro | On granitoid emplacement and related structures. A review[END_REF]Albertz et al., 2005;[START_REF] Chadwick | The mechanics ground deformation precursory to dome-building extrusions at Mt. St. Helens 1981-1982[END_REF]; Fig. 1A). Although we do not claim that our experiments reproduce faithfully all processes happening during magma mixing, we believe that they likely mimic the incipient stage of mixing at a microscale in natural hydrous magmas.

Due to faster kinetics resulting from the high water content (e.g. [START_REF] Dingwell | The effect of water on the viscosity of a haplogranitic melt under P-T-X conditions relevant to silicic volcanism[END_REF], the cooling rate of the vessel could not prevent the development of quench crystals whose occurrence partly obscures textural observations. Interpretation of experimental products is based on the first step syntheses that were ended by a drop quench, hence allowing better preservation of the HP/HT phase proportions. Although attention was paid to the development of fabric with applied strain, quench effects are too severe to extract any useful information in this field, as well as for the quantification of mixing intensity. The fact that plagioclase fabrics of the felsic layers are nearly isotropic (Fig. 5) should be viewed as a direct consequence of crystallisation during cooling in an almost isotropic stress field. Only in the lowest temperature experiments (600°C-700°C) does the fabric record the effect of deformation, showing that under the investigated conditions the mafic magma is unaffected, while plagioclases in the felsic magma are re-oriented at the vicinity of the contact with the mafic layer (Supplementary materials).

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

The presence of olivine and oxide in torsion products most likely reflects changes in the prevailing oxygen fugacity, as other parameters (T, P, H 2 O) were kept constant.

Experimental work performed on the Balos basalt from Santorini supports this observation: at NNO+3.5 and water-saturated conditions olivine is not stable whereas magnetite + clinopyroxene or magnetite + clinopyroxene + amphibole predominate at ≥1000ºC and 950ºC, respectively (Andújar et al., 2015). In contrast, at NNO, the stable mineral assemblage is olivine + clinopyroxene + plagioclase ± magnetite for temperatures > 985ºC; and amphibole co-crystallizes with these phases in the range 985-975ºC (Andújar et al., 2015). Thus, the occurrence of olivine in torsion experiments records the lower fO 2 conditions in the Paterson experiments relative to those prevailing in synthesis experiments, but without significantly changing the crystal fraction (Andújar et al., 2015).

4.2.

Interpretation of chemical results

Chemical data reflect interactions produced during the experiments, i.e. from the 1hour lap time before deformation to the cooling of the experiment. The diffusion distance of the main elements anticipated from the application of the diffusivity laws is consistent with what is observed in run products; for instance, the self diffusivity of Mg is about 10 -8 to 10 -9 cm²/s in hydrous felsic melts (van der Laan & [START_REF] Van Der Laan | Experimental interaction of granitic and basaltic magmas and implications for mafic enclaves[END_REF][START_REF] Mungall | Chemical diffusivities of 18 trace elements in granitoid melts[END_REF][START_REF] Zhang | Diffusion data in silicate melts[END_REF] and the expected diffusion distance after 2 hours is ~90 to ~360 µm (calculated using the relationships x = (Dt) 1/2 , x being distance in m, D diffusivity in m 2 /s and t time in s), in agreement with the chemical halo highlighted by elemental mapping (supplementary information) which extends over 200 µm at higher temperature (1000°C, PP261). Similarly, the expected distances by self diffusion of Ca (10 -8 to 2.10 -7 cm²/s) after 2 hours are 260 to 520 µm (van der Laan & [START_REF] Van Der Laan | Experimental interaction of granitic and basaltic magmas and implications for mafic enclaves[END_REF][START_REF] Watson | Diffusion in magmas at depth in the Earth: The effects of pressure and dissolved H 2 O[END_REF][START_REF] Zhang | Diffusion data in silicate melts[END_REF]. Therefore, assuming that the presence of crystals does not affect diffusion properties, chemical mixing occurring prior to deformation is likely to have been affected by strain.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

There are no significant compositional differences between clinopyroxenes crystallizing under either static or dynamic conditions whereas plagioclase and olivine define larger compositional intervals (Fig. 7), which may reflect the effect of mixing but also that of quench. Similar compositions of crystals whatever their localisation indicate efficient crystalmelt separation compared to the timescale of chemical equilibration. In static experiments olivine composition is Fo77-Fo70 and plagioclase is An80-An50 whilst in related torsion experiments, olivine and plagioclase compositions extend towards lower Mg (Fo64) and Ca (An40) values, suggesting that deformation enhances mixing process, probably via the entrainment of mafic minerals into felsic liquids. This is also supported by the wide compositional diversity of liquids which are sometimes more mafic than those obtained under synthesis conditions (Fig. 8). Clearly, however, it is also on liquid compositions that quench effects are maximised, which calls for caution in interpreting experimental trends solely in terms of mixing processes. Owing to profuse quench effects, chemical inferences from our experiments are necessarily very limited, though in agreement with observations made on dry compositions (Laumonier et al., 2014a), which show that the stacked, contrasted magmas produces a range of phase compositions, either solids or liquids. This is despite the relatively short run durations, < 4 hours. A large range of compositions has been also observed in previous experimental studies simulating mixing under static conditions (e.g., [START_REF] Johnston | Interaction of granitic and basic magmas: experimental observations on contamination processes at 10 kbar with H 2 O[END_REF].

4.3.

Effect of water on mixing

4.3.1.

Development of mingling/mixing features in hydrous magmas

The hydrous experiments described here document the progressive development of mingling and mixing textures over a 50°C temperature interval (under our experimental conditions). The first mingling feature appears at 950°C and corresponds to the mechanical

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

incorporation of isolated mafic crystals into the felsic magma (Fig. 6B), such as pyroxene and calcic plagioclase (An80-An92). When immersed in the felsic component, isolated plagioclases developed a sodic rim (An40 to An70), producing a normal zoning pattern (Fig.

& 7

). This contrasts with isolated plagioclases observed in dry conditions, characterized by resorption features and poorly developed rim (Laumonier et al., 2014a). At 950°C and higher temperatures, diffusion operates and can be observed in the felsic layers: the crystallization of plagioclase with more calcic composition around mafic layers results from the diffusion of calcium from the mafic component (Fig. 6). The effective diffusion of Fe and/or Mg is also marked by the halo of tiny Fe-Mg minerals in the PP296 experiment (Fig. 6C). The PP296 experiment can be compared to the PP156 run of Laumonier et al. (2014a): both lack mingling textures, despite being just 10-15°C below the temperature at which profuse mixing/mingling textures were produced. However, diffusion has clearly occurred in the hydrous experiment. Therefore, water allows mixing to operate even at temperatures slightly below the crystal rheological threshold (~0.50).

At 975 and 985°C, interface instabilities are generated (Fig. 6B &C) but the mafic layers still maintain their cylindrical shape in spite of a relatively high finite strain (γ = 5.1, PP296, Fig. 4). The first clear change in layer shape is observed in the runs performed at T>1000°C where both end-members are folded and stretched to boudin-like geometries, presumably because of the collapse of the touching crystal network. Hence, it appears that, at crystal fractions lower than 45%, mafic layers dislocate and produce enclaves, larger-scale lobes and embayments (e.g. [START_REF] Hodge | Experimental constraints on the deformation and breakup of injected magma[END_REF].

Origin of diktytaxitic textures

The presence of acicular microcrystals within the margins of enclaves is often interpreted as resulting from fast cooling of the hotter mafic enclave in its cold host (e.g. Bacon, 1986;[START_REF] Gerbe | Role of magma mixing in the petrogenesis of tephra erupted during the 1990-98 explosive activity of Nevado Sabancaya, southern Peru[END_REF][START_REF] Troll | Pre-eruptive magma mixing in ash-flow deposits of the tertiary Rum Igneous Centre, Scotland[END_REF][START_REF] Martin | Textural analysis of magmatic enclaves from the Kameni Islands, Santorini, Greece[END_REF]Davì et al., 2010). Such a fast crystallization leads to volatile exsolution and resulting bubbles which are entrapped within the crystals network, forming diktytaxitic texture (Fig. 6F & Fig. 9A).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
However, our experimental observations offer an alternative explanation: they suggest that the rim composed of crystals in many enclaves is not the result of a fast cooling only, but may come from the existence of a chemical gradient in the transition zone. As shown in dry experiments, the interface between chemically contrasted magmas is subject to chemical reactions (Laumonier et al., 2014a). Our hydrous experiments also produced a layer composed of plagioclase (fraction ~0.35), amphibole (~0.20), inherited pyroxenes (0.05) and dacitic melt (0.40). Most amphiboles and plagioclases in the layer are acicular but they are relatively large, approaching the size of crystals of the starting material (Fig. 6G), which suggest that they have started to grow before the quench, reflecting crystallization at P and T.

The absence of such a rim around the enclave produced in PP261 can be explained by the fact that chemical diffusion was erased by higher local strain, in spite of a Mg-enriched haplotonalite halo (Fig. S2). A chemical gradient between two hydrous magmas was also the cause of the development of a reaction zone such as reported by [START_REF] Johnston | Interaction of granitic and basic magmas: experimental observations on contamination processes at 10 kbar with H 2 O[END_REF]. In these experiments, a change in clinopyroxene composition with respect to its position to the interface was reported. Carroll & Wyllie (1989) also produced a reaction layer in their hydrous experiments between peridotite and granite at 1050 °C and 1.5 GPa, consisting of small pyroxenes. The weak chemical variation observed in clinopyroxenes in our experiments is possibly the result of the deformation, favoring the interaction zone and thus enhancing homogenization. Therefore, chemical gradients between enclaves and their magma host may play a role on the development of rims, which is not limited to the chilling effect of a hot mafic blob. Textural features from enclaves from Soufrière Hills supports such an hypothesis (Fig. 9); the external parts of a basaltic-andesitic enclave contain significantly more amphiboles, and a halo of plagioclases separates the enclave from the andesitic host,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

resembling the amphibole and plagioclase layers produced in our experiments (Fig. 6F & 11B). The orientation of the amphiboles at the periphery of the enclave is not perpendicular to the contact, but roughly parallel to it, indicating that some flow occurred after their crystallization. Therefore, amphibole is more likely to result from a chemical transition rather than from a quench process only. According to our results, such a texture is possible at a relatively low crystal fraction, which is consistent with the abundance of melt in both the enclave and its host.

Importance of the crystal fraction and viscosity contrast

The comparison with the dry experiments of Laumonier et al. (2014a) shows that the presence of water decreases the temperature at which magma mixing and mingling happens by 170°C; the widespread mafic magma remobilization indeed sets in suddenly at 1000°C, as compared to 1170°C under dry conditions. Such a difference is close to the one between the liquidus temperature between hydrous and dry basaltic magmas (Andùjar et al., 2015). Water shifts down the liquidus temperature, thus the crystallisation and consequently magma mixing capacities. Our experiments are thus in agreement with the concept of a rheological threshold, as proposed by e.g. Philppots et al. (1998) and Martin et al (2006a) and more recently by Laumonier et al. (2014a). The onset of dismembering of the whole mafic layer occurs between 985 and 1000°C, corresponding to a crystal fraction close to 45% and similar to mixing experiments at dry conditions (Laumonier et al., 2014a). At 985°C and lower temperatures, the crystal fraction is high enough to build a rigid (touching) crystal network controlling the mush rheology that even strain cannot break at our experimental conditions. A modest temperature rise of 15°C results in the diminution of the crystal fraction, allowing the dismembering of the crystal network and subsequent magma mingling. However, in detail, a few minor and discrete mingling textures (isolated crystals, crenulated contact) do occur before the onset of large scale dislocation of the mafic layers, such as in the charge at 950°C

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

(Fig. 6A), establishing a more progressive transition between mixed and unmixed magmas than that observed under dry conditions. Such a progressive transition might result in the lower melt viscosity under hydrous conditions.

At subsolidus tempertures (600 and 715°C), the viscosity contrast between the mafic and felsic melts does not exceed 0.3 log unit (Laumonier et al.;[START_REF] Laumonier | On the conditions of magma mixing and its bearing on andesite production in the crust[END_REF]. However, neither mingling nor mixing features were produced in those conditions. In fact, the viscosity contrast between magmas (melt plus crystals) is about 1 log unit, which seems large enough to prevent mixing (under our deformation conditions), as observed under dry conditions (Laumonier et al., 2014a). Therefore, our results support the important role of the crystal fraction, which controls bulk viscosity, on the remobilization and mixing of dry and hydrous magmas. Within the range of strain rates explored here, a difference in the viscosities of interacting magmas (either hydrous or dry) of about an order of magnitude is enough to prevent them from being mixed.

Lastly, the material extruded out of the samples during PP261 and PP295 experiments bears texturally complex mingling/mixing features (see supplementary materials). Although strain, strain rate and also the deformation geometry cannot be precisely determined, the other parameters (in particular T) are similar to those reigning in the sample. The PP295 experiment conducted at 975°C shows that the rheological transition has been upset in the case of extrusions: there, the crystal network hindering magma mingling in the interior sample was locally broken by a different regime of deformation when the magma was extruded through a hole in the jacket, creating filament and rounded shape textures. On this basis it can be concluded that mixing processes are considerably enhanced by elevated rates of deformations (compare Fig. 4 and Fig. S2), as found in other experimental studies or by a particular geometry of the interface forcing both magmas to be deformed (e.g. De [START_REF] De Campos | Enhancement of magma mixing efficiency by chaotic dynamics :an experimental study[END_REF]Morgavi et al., 2013a, b). The textures produced in the extrusions are similar to those

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
observed in the main samples at T>1000°C (PP261 & PP293) but better developed, presumably as a result of the high strain underwent by the magma (also including the production of stretched filaments; Fig. S2). Hence, we suggest that a relatively high crystal fraction (~0.50) may not be adequate to mixing and mingling in magma reservoir, but still feasible (until a certain fraction) in the volcanic conduit where higher strain and strain rate occur or in highly dynamic reservoirs where convective rates allow high strain rates to be attained [START_REF] Woods | Magma mixing triggered during volcanic eruptions[END_REF][START_REF] Pritchard | Eruption-triggered mixeing of extra-caldera basalt and rhyolite complexes along the East Gallatin-Washburn fault zone[END_REF].

4.4.

Implications for reservoir replenishments

In the light of these results we suggest that reservoirs being slowly replenished will have less possibility to produce well-mixed magmas, in accord with the results of analogical experiments (e.g. [START_REF] Turner | Convection and Mixing in Magma Chambers[END_REF]Jellinek et al., 1999;[START_REF] Hodge | Experimental constraints on the deformation and breakup of injected magma[END_REF]. Strain rates in convecting reservoirs are inferred to be quite low, <10 -6 s -1 (e.g., [START_REF] Spera | Rheology of melts and magmatic suspensions. 1. Design and calibration of a concentric cylinder viscometer for application to rhyolitic magma[END_REF]Albertz et al., 2005) which therefore suggests that mixing between rheologically contrasted magmas (i.e. with viscosities differing by more than 1 order of magnitude) are unlikely to be efficient at this location. In contrast, mixing during vigorous, or forcefull, injection may lead to thoroughly mixed material or hybrids, as illustrated by the extrusions and experiments involving natural basalt and rhyolite differing in viscosity by ~ 4 orders of magnitude, which nevertheless were mingled and mixed at higher temperature (T>1350°C) and strain rates (~10 - 2 s -1 ; De [START_REF] De Campos | Enhancement of magma mixing efficiency by chaotic dynamics :an experimental study[END_REF]Morgavi et al., 2013a[START_REF] Morgavi | Time evolution of chemical exchanges during mixing of rhyolitic and basaltic melts[END_REF]. Considering that the P-T-H 2 O conditions investigated here cover those likely to prevail in most arc magmas reservoirs [START_REF] Scaillet | Phase equilibrium constraints on the viscosity of silicic magmas, 1. Volcanic-plutonic comparison[END_REF]Andújar and Scaillet, 2012), our experiments constrain the temperature at which mixing can operate. As a general rule, our results suggest that injection of a mafic forerunner into a crystal-rich silicic host will not produce hybrid products if the equilibration temperature of the blended mixture remains below 950°C. The final temperature of equilibration depends, inter alia, on the relative masses of interacting end-members (Sparks [START_REF] Sparks | Thermal and mechanical constraints on mixing between mafic and silicic magmas[END_REF][START_REF] Frost | Field, chemical, and physical constraints on mafic-felsic magma interaction in the Lamarck Granodiorite, Sierra Nevada, California[END_REF][START_REF] Laumonier | On the conditions of magma mixing and its bearing on andesite production in the crust[END_REF] and so there is no unique solution to this problem. However, it can be anticipated that the slow replenishment via small mafic inputs of a large felsic reservoir, will not give rise to prolonged and widespread thermal excursions, and hence will lead to conditions unfavourable to mixing: slow replenishment will more likely end up producing stratified intrusions rather than a homogeneous whole body. Consequently, we suggest that mixing features encountered in highly crystallized magmas (> 50% crystals) may have been produced before their crystallization (closer to liquidus conditions, [START_REF] Scaillet | Phase equilibrium constraints on the viscosity of silicic magmas with implications for mafic-silicic mixing processes[END_REF][START_REF] Caricchi | Insights into the mechanisms and timescales of pluton assembly from deformation patterns of mafic enclaves[END_REF] or in extreme strain rate environments such as in volcanic conduits where the strain rates, during eruption, exceed the average plutonic ones (e.g., [START_REF] Koyaguchi | The dynamics of magma mixing in a rising magma batch[END_REF][START_REF] Pritchard | Eruption-triggered mixeing of extra-caldera basalt and rhyolite complexes along the East Gallatin-Washburn fault zone[END_REF]. That the mixing efficiency changes rapidly over a small temperature interval suggests in turn that small temperature variations may give rise to abrupt changes in the fluid dynamics of injected reservoirs. Eruptions triggered by mafic inputs may thus record the case of felsic reservoirs lying close to the rheological threshold prior to mafic injection.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

CONCLUSION

We present the results of the first experiments on mixing and mingling of hydrous

ACCEPTED MANUSCRIPT

(within one order of magnitude), and at temperatures shifted by ~200°C when compared to similar but dry magmas (Laumonier et al., 2014a). Finally, the transition between mixing and unmixing regimes appears more progressive in water-bearing magmas. However, mingling can occur at crystal fractions higher than the critical crystallinity, in particular under extreme strain rate conditions. Those conditions prevail in volcanic conduits, suggesting that mixing involving highly crystallized magmas may be the consequence of volcanic eruptions.

β4 providing 3

 3 interfaces between end-members. The faces of each disk were ground to have parallel sides (thickness variation lower than 0.02 mm) and polished to reduce interface irregularities. The stack of disks (5.33 to 10.93 mm thick) was wrapped in a platinum foil, which did not significantly chemically interact with the sample during the experiments. This assemblage was inserted in turn in a copper or iron jacket and sandwiched between pistons to be located in the isothermal zone of the furnace (± 2°C on 35 mm length determined during furnace calibration). Although in nature strong thermal gradients likely exist between mixing magmas, in this study we tried to avoid any temperature gradient to better constrain the role of melt fraction on mixing. Our experiments hence model the conditions of mixing once thermal equilibrium has been reached between the end-members. The column so prepared was inserted in a Paterson press (Paterson instrument, Australian Scientific Instruments) at ISTO to perform torsion experiments at constant strain rate. Such a deformation configuration reproduces simple shear conditions and allows the deformation of one material only without forcing the whole stack to deform, such as happening during the replenishment of a felsic A C C E P T E D M A N U S C R I P T

  magmas performed at P-T-H 2 O conditions relevant to arc contexts (P = 300 MPa, T < 1020°C and water contents up to 6.5 wt.% in the melt). Our low strain rate torsion experiments on hydrated wafers of felsic/mafic magmas fully reproduce the mingling textures observed in natural magmas (isolated crystals, boudin-like structures, enclaves) and mixing textures like the development of chemical halo and transition layers. Such layers consist of a dense population of plagioclase, and acicular crystals, melt and bubbles, altogether resembling the diktytaxitic texture often observed in natural enclave margins that may result from the rapid crystallization of a chemical transition between magmas. Mixing and mingling of hydrous magmas appears feasible at low crystal fraction (< 50 vol.%) and low viscosity contrast A C C E P T E D M A N U S C R I P T

Figure 1 :

 1 FIGURE CAPTIONS

Figure 2 :

 2 Figure 2: SEM picture of starting synthetic suspensions of haplotonalite (felsic) and basalt

Figure 3 :

 3 Figure 3: Reconstructed panorama of deformed samples still embedded in their jacket. Stack

Figure 4 :

 4 Figure 4: SEM pictures of large sections showing the relative position of layers after torsion

Figure 5 :

 5 Figure 5: (A) Phase proportions of the mafic layers. (B) Phase proportions for the area

Figure 6 :

 6 Figure 6: Details on particular textures produced in torsion experiments (see Fig. 4 for

Figure 7 :Figure 8 :

 78 Figure 7: Molar anorthite content of plagioclase in 1000°C and "reheated" Syntheses, static

Figure 9 :

 9 Figure 9: Natural examples of magma mixing from (A) Adamello Massif (Blundy and Sparks,

FUNDING

This work was supported by the French Agence National de Recherche [ANR-08CEA080 to B. S], and Equipex PLANEX.

ACKOWLEDGEMENTS

We greatly thank I. Di Carlo for analytical support and E. Lemoing and P. Teulat for technical assistance. The reviews of K. Russel, M. Pistone and an anonymous reviewer helped in focusing the topic of this manuscript. We are also grateful to K. Mezger and the guest editors of this special issue for editorial assistance.

Albertz, M., Paterson, S.R., Okaya, D. (2005). Fast strain rates pluton emplacement: Magmatically folded leucocratic dikes in aureoles of the Mount Stuart Batholith, Washington, and the Tuolumne Intrusive Suite, California, Geological Society of America Bulletin 117, 450-465. Anderson, A.T. (1976). Magma mixing: petrological process and volcanological tool, Journal of Volcanology and Geothermal Research 1, 3-33, doi:10.1016/0377-0273(76)90016-0.

Andújar, J., Scaillet, B. (2012). Relationships between pre-eruptive conditions and eruptive styles of phonolitestrachyte magmas. Lithos, doi:10.1013/j.lithos.2012.05.009.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Andújar, J., Scaillet, B., Druitt, T.H., Pichavant, M. (2015). Differentiation conditions of a basaltic magma from Santorini, and its bearing on the production of andesite in arc settings. Journal of Petrology 56, 765-794.

Armienti, P., Barberi, F., Bizouard, H., Clocchiatti, R., Innocenti, F., Metrich, N., Rosi, M., Sbrana, A. (1983).

The phlegraean fields: Magma evolution within a shallow chamber. Journal of Volcanology and Geothermal

Research 17 (1)(2)(3)(4), 289-311, DOI: 10.1016/0377-0273(83)90073-2. Bacon, C. R. (1986). Magmatic Inclusions in Silicic and Intermediate Volcanic Rocks. Journal of Geophysical Research 91 (B6), 6091-6112, doi:10.1029/JB091iB06p06091. Barnhoorn, A., Bystricky, M., Kunze, K., Burlini, L., Burg, J.P. (2005). Strain localization in bimineralic rocks: Experimental deformation of synthetic calcite-anhydrite aggregates. Earth and Planetary Science Letters 240, 748-763.

The felsic Dry* glass is the nominal composition given by Schott, A.G., Hydrated* glass is calculated. Standard deviations in terms of least unit cited are indicated by italic font.

Compositions were acquired using a Cameca SX 50 microprobe at 15kV and 6nA. Glass analyses were done with defocused beam (covering an area of ~ 25μm 2 ) whereas a focused beam (~ 1μm 2 ) was used for crystal analyses. Glass and minerals are normalized to 100% anhydrous. Water content was calculated by mass difference. Mg# is the molar ratio Mg/(Mg+Fe*) and An is the molar ratio Ca/(Ca+Na). (duration of the experiment, including equilibration lap time and deformation time). ρ1, β2, ρ3, β4 are the layers composing the sample (see Fig. 1). Dimensions (diameter, initial length, thickness of each layer and length after experiment) are given in mm.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Highlights (mandatory)

We performed the first dynamic mixing experiments between hydrous magmas Mixing textures are isolated crystals, boudins, enclaves, "diktytaxitic" layers Mixing of hydrous magmas occurs at T=170°C lower than between similar but dry magmas

Mingling and mixing appears feasible at crystal fraction ϕ < 50 vol%

Mingling may occur at ϕ>50 vol% under higher strain rate conditions (volcanic conduit)


Fanny Morin

Hadrien Courtecuisse

Matthieu Chabanas

Yohan Payan

F Morin Ab

Rest shape computation for highly deformable model of brain

Keywords: biomechanical simulation, rest shape, collisions, brain model

Rest shape computation for highly deformable model of brain.

Introduction

In medical simulations, finite element (FE) models are often extracted from medical images. This initial geometry usually corresponds to a state of the model at equilibrium between external and internal forces. Therefore, an important issue in biomechanics is to estimate the initial stress of FE model before running a simulation.

A typical example concerns gravity. Most biomechanical models ignore the gravity and directly run the simulation from the initial mesh with a null stress. However, this will obviously yield to different results, especially for living tissue that are usually soft and undergo large deformations. When the gravity force is applied to a model, it will shrink down unless other conditions are accounted for: boundary conditions, internal forces (e.g. muscle contraction) or other external forces or pressures. A rest shape must therefore be determined so that the model converges to the targeted equilibrium state when loads are applied.

Some methods propose to directly compute these loads whereas others prefer to estimate the rest shape of the model that will imply pre-stressed model at equilibrium. These inverse problems admit an exact formulation as a set of differential equations [START_REF] Govindjee | Computational methods for inverse finite elastostatics[END_REF] which has the benefit to find an exact and unique solution. However, the derivation of this formulation is cumbersome and depends on biomechanical laws of the model. On the other hand, [START_REF] Sellier | An iterative method for the inverse elastostatic problem[END_REF] proposed an iterative technic in order to compute rest position. That position may not be exact but its formulation doesn't rely on biomechanical laws.

In this paper, a method is presented in the context of a biomechanical model of the brain that will be used to model brain-shift during surgery. This phenomenon is mainly caused by a loss of cerebro-spinal fluid (CSF) surrounding brain. Our approach is based on [START_REF] Sellier | An iterative method for the inverse elastostatic problem[END_REF], with the addition of contact forces which enable the iterative method to converge towards an equilibrium state.

Methods

Biomechanical model

We used FE model with co-rotational approach [START_REF] Müller | Stable real-time deformations[END_REF] for the brain and collisions with skull are handled with Lagrangian Multipliers as explained in [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF].

Gravity and CSF pressure are applied on the model. CSF pressure is modeled as a uniform force field applied on the outer nodes of the mesh along their normal.

Computation of the rest shape

For this purpose, we extended the method proposed by [START_REF] Sellier | An iterative method for the inverse elastostatic problem[END_REF] In his method, [START_REF] Sellier | An iterative method for the inverse elastostatic problem[END_REF] proposed to compute ∆ in step 4) as the maximum of the ‖∆ 𝑖 ‖. But, this maximal value is often reached for few nodes set in highly constrained areas (such as contact regions or areas close to fixed points). Thus, computing the average value has the benefit to give a more global regularization to the model.

Moreover, as stated by [START_REF] Sellier | An iterative method for the inverse elastostatic problem[END_REF] in his proof of convergence, his method converges when deformations are small enough. In our case, due to the low stiffness of the brain, the biomechanical model undergoes large deformations when gravity and CSF pressure are applied. The model oscillates and finally diverges. Thus, no equilibrium state can be found.

Therefore, we propose in step 2) to account for collisions with the skull during this pre-stressing process. The normals of the skull contact surface are oriented toward the inside of the cavity, with collisions handled on direct orientation only. This means a node cannot move out of the skull boundaries but is allowed to cross it from the outside to the inside. This point is crucial since some of the rest positions computed with the above algorithm may be outside the skull.

*Corresponding author. Email: fanny.morin@imag.fr

Handling such contacts limits the deformations of the biomechanical model within the skull boundaries, which finally enables the iterative algorithm to converge towards an equilibrium state.

Results and discussion

Rest shape

The brain model was segmented from MR images using ITK-SNAP (www.itksnap.org), including main anatomical parts such as tumor area, brainstem, cerebellum and skull. Then, the meshing step was completed with CGAL (www.cgal.org) which allows creating labelled tetrahedral FE meshes from segmentation images. Finally, the FE model was developed with the simulation framework Sofa (www.sofa-framework.org). This is illustrated in Figure 1.

Figure 1: visual (a) and collision (b) models in Sofa

For the simulation, we used a FE mesh with 5266 elements and 1081 nodes. The Young modulus was fixed to 1.5 kPa for the brain and 10 kPa for the tumor area whereas Poisson ratio was uniformly set to 0.45 [START_REF] Schiavone | In vivo measurement of human brain elasticity using a light aspiration device[END_REF]). Nodes of brainstem and cerebellum were fixed.

Rest shape (fig. 2(a)) was found in 6 iterations, with average error ∆ equal to 0.746mm. Figure 2(b) shows that highest errors are located on the surface of the model, close to contact areas. Nevertheless, inside the volume, the errors are low and the shape is well interpolated. Finally, the heterogeneous region of the tumor also undergoes high stress because of its higher stiffness.

Impact of the pre-stressing step

In this section, we would like to show the importance of having pre-stressed models taking external forces into account.

On that purpose, we used two FE models. The first model was pre-stressed with our method, accounting for gravity and CSF pressure. Conversely, the second one was not pre-stressed, and thus, did not account for external forces. Then, we applied the same force field on each model.

Results show that both models deformed on the same direction but with different amplitude of displacement (

Conclusions

This paper presents a method to compute rest shape for highly deformable models. This technique is based on Sellier (2011) algorithm, with the addition of collisions with a surface surrounding the simulated mesh. This was presented on a specific biomechanical model of the brain, nevertheless the method could be generalised to other contexts. Moreover, we also presented results showing the importance of taking external forces into account for medical biomechanical simulations.

Figure 2 :

 2 Figure 2: (a) simulated model (grey) with rest positions (blue) and target positions (red) (b) error repartition on FE mesh from yellow (close to 0 mm) to dark blue (2.69 mm) (c) von Mises stress (in Pa) on a slice of the mesh at equilibrium Figure 2(c)shows von Mises stress repartition on a slice of the mesh at equilibrium. The blue area corresponds to the fixed cerebellum part: stress is quasi null here. Above this region, there is a red area with high stress due to the gravity acting on the mesh.

  . With 𝑋 𝑖 𝑖𝑛𝑖 the initial segmented position of node 𝑖, 𝑋 𝑖 𝑡𝑎𝑟 the target position, we look for the rest

	position 𝑋 𝑖 𝑟𝑒𝑠𝑡 so that the equilibrium state 𝑋 𝑖 𝑒𝑞 is as close as possible to 𝑋 𝑖 𝑡𝑎𝑟 . Our algorithm is:
	1) ∀𝑖, 𝑋 𝑖 𝑡𝑎𝑟 = 𝑋 𝑖 𝑟𝑒𝑠𝑡 = 𝑋 𝑖 𝑖𝑛𝑖
	2) Simulate until equilibrium is reached,
	accounting for contacts 3) ∀𝑖, ∆ 𝑖 = 𝑋 𝑖 𝑒𝑞 𝑡𝑎𝑟 -𝑋 𝑖 ∑‖∆ 𝑖 ‖ 4) ∆ = 𝑛𝑏 𝑁𝑜𝑑𝑒𝑠 5) If ∆ > 𝜀, ∀𝑖 𝑋 𝑖 𝑟𝑒𝑠𝑡 + ∆ 𝑖 and go to step (2)
	else accept (𝑋 𝑖

𝑟𝑒𝑠𝑡 ) as the final rest positions.

Table 1 )

 1 : as expected, the pre-stressed model appears much stiffer than the other one.

		Mean displ.	Max. displ.
	1. Pre-stressed	3.0347 mm	6.9218 mm
	2. Not pre-stressed	4.7843 mm	13.3968 mm
	Increase	57.7%	93.5%

Table 1 :

 1

impact of an force field applied on prestressed (

1

) and not pre-stressed model (

2

)

Acknowledgements

The authors would like to thank Pr. O. Palombi from the department of Neurosurgery, Grenoble Hospital. This work was partly supported by the French ANR within the Investissements d'Avenir program (Labex CAMI) under reference ANR-11-LABX-0004.