domain: spi sde

text:
Emmanuil Veveakis

Jean Sulem

Ioannis Stefanou

Modeling of fault gouges with Cosserat Continuum Mechanics: Influence of thermal pressurization and chemical decomposition as coseismic weakening mechanisms

Keywords: Shear Heating, Cosserat Continuum, Reaction kinetics, Undrained adiabatic shearing, Micro-Inertia

come

Modeling of fault gouges with Cosserat Continuum Mechanics: Influence of thermal pressurization and chemical decomposition as coseismic weakening mechanisms

Introduction

Geomaterials exhibit failure in various patterns, like the propagation of single fractures in otherwise intact rock, crack growth in presence of multiple interacting fractures, initiation and percolation of damage zones under quasistatic and dynamic loading, strain localization into shear-or compaction bands.

The theory describing localized failure has been a direct extension of the classic Mohr's strength of materials theory, and is applied for geomaterials using the so-called Thomas -Hill -Mandel shear-band model [START_REF] Hill | Acceleration waves in solids[END_REF]; [START_REF] Mandel | Conditions de stabilite et postulate de drucker[END_REF]), which was introduced in the early 60s and it was widely publicized by the paper of [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF]. More recently, the mathematical formulation of bifurcation and post-bifurcation phenomena and related instabilities was summarized by Vardoulakis and Sulem (1995) to form the basis of a contemporary continuum theory of failure of geomaterials.

Apart from the theoretical point of view of failure, shear banding has been claimed as the primary failure pattern for faults [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]) and landslides (Vardoulakis (2002a)). However, field evidence from exhumed faults reveals that seismic events take place in even narrower zones within shear-bands, formed from post-failure evolution [START_REF] Chester | Ultracataclasite structure and friction processes of the punchbowl fault, san andreas system, california[END_REF]), dominated by weakening mechanisms that are considered to be mainly thermal of origin [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]). Indeed, field observations suggest that slip in individual events may then be extremely localized, and may occur primarily within a thin shear zone, < 1 -5 mm thick, called the principal slipping zone (PSZ) [START_REF] Sibson | Thickness of the Seismic Slip Zone[END_REF]). This localized shear zone lies within a finely granulated fault core of typically tens to hundreds millimeter thickness, as shown schematically in the idealized description of Fig. 1.

Evidence for that morphology can be found in the North Branch San Gabriel fault [START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF]) and Punchbowl [START_REF] Chester | Ultracataclasite structure and friction processes of the punchbowl fault, san andreas system, california[END_REF]) fault of the San Andreas system in southern California, the Median Tectonic Line fault in Japan [START_REF] Wibberley | Internal structure and perme-43 ability of major strike-slip fault zones: The median tectonic line in mid prefecture, southwest japan[END_REF]), and in the Hanaore fault in southwest Japan [START_REF] Noda | Thermal pressurization and slip-weakening distance of a fault: An example of the Hanaore fault, Southwest Japan[END_REF]). Note also that from drilling cores in the Aigion system, central Greece, the fault core of clay size particles consisting of finely crushed radiolarites, extended to about 1 m, was found to be intercepted by a "fresh"distinct slip surface of sub-millimiter size [START_REF] Sulem | Stress orientation evaluated from strain localisation analysis in Aigion fault[END_REF]; [START_REF] Cornet | Drilling through the active Aigion fault: the aig10 well observatory[END_REF]).

As perhaps the best characterized case, a thin principal slipping zone was identified along an exposure of the Punchbowl fault at 2 to 4 km depth, which has accommodated 44 km of slip [START_REF] Chester | Ultracataclasite structure and friction processes of the punchbowl fault, san andreas system, california[END_REF]). Also, a thin central slip zone was identified by [START_REF] Wibberley | Internal structure and perme-43 ability of major strike-slip fault zones: The median tectonic line in mid prefecture, southwest japan[END_REF] along the Median Tectonic Line fault. The nominal thickness of such a shear zone is usually of the order of sub-millimeter and, as described by [START_REF] Smith | Principal slip zones in limestone: Microstructural characterization and implications for the seismic cycle (Tre Monti fault, Central Apennines[END_REF], it could be met in calcite-bearing fault rocks even in the absence of pressure-solution and precipitation, providing thus field evidence of extreme localization during coseismic slip in natural faults. On the other hand the [START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF] and [START_REF] Wibberley | The Internal Structure of Fault Zones -implications for mechanical and fluid flow properties[END_REF]. The fault core, d, is of the order of centimeters, whereas the Principal Slipping Zone (PSZ) D is of the order of millimeter or less.

thickness of the cataclastic core may vary up to a few centimeters.

In geomechanics there are several studies focusing on the thickness of shear bands in soils. [START_REF] Morgenstern | Microscopic structures in kaolin subjected to direct shear[END_REF] calculated that shear bands for a clay, kaolin may have thicknesses of order 10 -100 times a clay platelet size (∼ 0.5 µm). For granular materials on the other hand usually use the mean particle diameter, d 50 , is used as a scaling measure, defined for nonuniform size distributions such that 50 % by weight of the particles have larger size.

Under purely mechanical effects, [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF] calculated for sands the thickness of the initial shear band (with thickness d in Figure 1) to be around d ∼ 16d 50 . In any case, this zone seems to be rather wide to be considered as the primary slip surface [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]). Indeed, [START_REF] Chester | Microscale characterization of natural and experimental slip surfaces relevant to earthquake mechanics[END_REF] showed that the nominal thickness of the Punchbowl's cataclastic core varies from 0.6 to 1.1 mm. However, as [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF] mentions: "within this ∼ 1 mm thick nominal shear zone, most of the shearing seems to have been accommodated within a zone of extreme shear localization having an apparent thickness of 100 -300 µm ". [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF] showed that d 50 is not always a characteristic length scale for cataclastic and ultracataclastic formations, since apparently d 20 -d 30 seems to be more appropriate for modeling the PSZ. Thus, in principle the PSZ represents a very fine shear localization, consisting of ultrafine particles that have usually undergone mechano-chemical degradation (gelification, decarbonation and dehydration reactions, melting, as thoroughly discussed by Di [START_REF] Toro | Fault lubrication during earthquakes[END_REF]). Therefore, other mechanisms that are responsible for the formation of these structures and the dynamic weakening of faults should be looked after.

Since fault zones usually exhibit the presence of fluid interacting with the rock, hence inducing coupled effects including shear heating and pore fluid pressurization, such mechanisms are expected to act as weakening mechanisms during fault slip [START_REF] Sibson | Interaction between temperature and pore-uid pressure during earthquake faulting a mechanism for partial or total stress relief[END_REF]; [START_REF] Lachenbruch | Frictional heating, uid pressure, and the resistance to fault motion[END_REF]; [START_REF] Wibberley | Earthquake slip weakening and asperities explained by thermal pressurization[END_REF]). These mechanisms have also been suggested for weakening in catastrophic landslides (Vardoulakis (2002a); [START_REF] Veveakis | Thermoporomechanics of creeping landslides: The 1963 Vaiont slide, Northern Italy[END_REF]; [START_REF] Goren | Long runout landslides: The role of frictional heating and hydraulic diffusivity[END_REF]; [START_REF] Pinyol | Criteria for rapid sliding II.thermo-hydromechanical and scale effects in Vaiont case[END_REF]; [START_REF] Goren | The long runout of the Heart Mountain landslide:heating, pressurization, and carbonate decomposition[END_REF]).

Thermo-poro-mechanical couplings due to shear heating [START_REF] Sulem | Thermoporo-mechanical properties of the aigion fault clayey gouge -application to the analysis of shear heating and fluid pressurization[END_REF]; [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]) can be also associated to chemical effects such as dehydration of minerals or decomposition of carbonates, theoretically studied recently by [START_REF] Brantut | Thermochemical pressurization of faults during coseismic slip[END_REF]Brantut et al. ( , 2011)); [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF]; [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF] and reported to take place in real faults [START_REF] Hirono | A chemical kinetic approach to estimate dynamic shear stress during the 1999 Taiwan Chi-Chi earthquake[END_REF]) and experiments at laboratory conditions [START_REF] Han | Ultralow friction of carbonate faults caused by thermal decomposition[END_REF] In a recent paper, [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF] have shown that by resorting to a Cosserat continuum (Vardoulakis and Sulem (1995)) for describing the kinetics of a granular material, a preferred wave length for the instability mode of a layer under shear can be selected. They showed that instability can occur even in the hardening regime of the stress-strain curve, due to thermal pressurization, confirming the claim that thermal in origin mechanisms play a major role in the formation and evolution of the PSZ. The localizing structures obtained were comparable in thickness with the ones observed in nature, rationalizing thus the viscous and second gradient regularization initially proposed by [START_REF] Vardoulakis | Stability and bifurcation of undrained, plane rectilinear deformations on water-saturated granular soils[END_REF][START_REF] Vardoulakis | Dynamic stability of undrained simple shear on watersaturated granular soils[END_REF][START_REF] Vardoulakis | Deformation of water-saturated sand: II. effect of pore water flow and shear banding[END_REF] for regularizing these type of problems. In this work we extend this analysis by including chemical reactions that may take place inside a rapidly deformed shear-band and study the impact of certain decomposition reactions to the dynamic weakening and localization of deformation in granular media.

Problem Formulation

In classical Continuum mechanics, where the continuum is allowed only for translational degrees of freedom (displacements), one may be able to predict the onset of a shear band [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF]). However it can be shown that the instability tends to localize in a strip of zero thickness.

Classical theories are thus unable to simulate the localization of the deformation in a shear band of finite thickness. This drawback can be traced back to the fact that conventional constitutive models do not contain material parameters with dimension of length which can scale the shear band thickness (internal length).

Thus, a fault is usually modelled as a discontinuity surface for the velocity fields. Such an oversimplifying assumption however would overestimate significantly the heat produced in a fault, and thus the temperature (as shown by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]), imposing difficulties in our effort to determine the processes taking place during dynamic weakening of faults. To raise this discrepancy, the introduction of a thin but finite-in-width slipping zone of intense shear is requested. [START_REF] Vardoulakis | Stability and bifurcation of undrained, plane rectilinear deformations on water-saturated granular soils[END_REF][START_REF] Vardoulakis | Dynamic stability of undrained simple shear on watersaturated granular soils[END_REF][START_REF] Vardoulakis | Deformation of water-saturated sand: II. effect of pore water flow and shear banding[END_REF] showed that introducing viscous (rate-dependent) considerations, combined with second gradient plasticity models, may alleviate the ill-posedness of the problem. Later, [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF] showed that indeed non-linear viscous responses provide the necessary counterbalancing effect to stabilize the problem of simple shear.

Early approaches by [START_REF] Hill | Acceleration waves in solids[END_REF] and [START_REF] Mandel | Conditions de stabilite et postulate de drucker[END_REF] envisioned such a structure with width d, which however could only be artificially imposed, since -as said before-in classical continua d is undetermined. It seemed as if the solutions obtained from classical continua correspond to the limit d → 0 of a corresponding higher grade extension of it (see Vardoulakis and Sulem (1995), p. 268 for a comprehensive discussion). Following this conclusion, as well as experimental observation of the influence of the effect of particle irregularities and rotations inside the shear bands, [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF] resorted to concepts from Cosserat Continuum mechanics to naturally determine d. Thus, by allowing both particle displacements and particle rotations, they accounted for the angularities of the grains, and enriched the classical continuum with additional kinematic and static fields, and calculated for sands that d ∼ 16d 50 . Other mechanisms leading to localization, and initially neglected by [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF], are reported to be flash heating and powder lubrication [START_REF] Toro | Fault lubrication during earthquakes[END_REF][START_REF] De Paola | Fault lubrication and earthquake propagation in thermally unstable rocks[END_REF]; [START_REF] Han | Strong velocity weakening and powder lubrication of simulated carbonate faults at seismic slip rate[END_REF][START_REF] Han | Granular nanoparticles lubricate faults during seismic slip[END_REF][START_REF] Tisato | Experimental investigation of flash weakening in limestone[END_REF] or grain breakage [START_REF] Nguyen | Nonlocal regularisation of a model based on breakage mechanics for granular materials[END_REF]).

In this study we follow [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF] and consider a layer of saturated granular material of thickness D that is sheared in plane strain (Fig. 2) in such a way that there is no extensional strain in the x 1 -direction and displacements u 1 and u 2 of material points vary only with x 2 (and time t). The additional rotational degree of freedom of the considered 2D Cosserat continuum is ω c . The grains themselves are not allowed for breakage or wear, excluding thus breakage and flash heating, since the scope of this study is to provide insights on the thickness of the PSZ at different temperature regimes.

Constitutive equations of elasto-plasticity

The incremental constitutive equations for the considered 2D Cosserat continuum are derived from the Muhlhaus -Vardoulakis plasticity model [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF]). The main feature of the model is that a 2D flow theory of plasticity for granular media with Cosserat microstructure can be derived by keeping the same definitions for the yield surface and the plastic potential as in the classical theory and by generalizing appropriately the stress and strain invariants involved in these definitions. The details of the model for isothermal deformation are given in Vardoulakis and Sulem (1995) and in [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]. They are briefly recalled in the following.

In a two-dimensional Cosserat Continuum each material point has two translational degree of freedom (v 1 , v 2 ) and one rotational degree of freedom ωc . The four components of the rate of the non-symmetric deformation tensor are given as

ǫ11 = ∂v 1 ∂x 1 ; ǫ12 = ∂v 1 ∂x 2 + ωc (1) ǫ21 = ∂v 2 ∂x 1 -ωc ; ǫ22 = ∂v 2 ∂x 2 ,
and the two components of the curvature of the deformation rate (gradient of the Cosserat rotation rate)

κ1 = ∂ ωc ∂x 1 ; κ2 = ∂ ωc ∂x 2 . ( 2
)
Moreover we decompose the stresses and strain rates into spherical and deviatoric parts, σ ij = s ij + σ kk δij/2 and ǫij = ėij + ǫkk δ ij /2, where δ ij is the Kronecker's delta. The following generalized stress and strain invariants are used:

σ = σ ii 2 ; τ = h 1 s ij s ij + h 2 s ij s ji + h 3 m k m k /R (3) ǫp = ǫp kk ; γp = g 1 ėp ij ėp ij + g 2 ėp ij ėp ji + R 2 g 3 κk κk ,
where

h i = [3/4, -1/4, 1], g i = [3/2, 1/2, 1]
for the so-called static Cosserat model, and R is the internal length (Vardoulakis and Sulem (1995)).

Besides the four components of the non-symmetric stress-tensor σ ij , there are two couple stresses m k (k = 1, 2). As in classical small-strain plasticity theory, the deformation is decomposed into elastic and plastic parts

ǫij = ǫe ij + ǫp ij ; κk = κe k + κp k . (4)
Following [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF] we assume a Terzaghi decomposition of the stress tensor in effective stress and pore pressure (σ ij = σ ′ ij -pδ ij ) and assume a Coulomb yield stress and plastic potential for the effective stresses:

F = τ + µ(σ + p), (5)
Q = τ + β(σ + p),
where the mobilized friction coefficient µ and the dilatancy coefficient β are functions only of the accumulated plastic strain γ p . Under these considerations, the rate thermo-poro-elasto-plastic relationships are expressed as

follows γ = τ G + γp (6) ǫ = 1 K ( σ + ṗ) + α s Ṫ + ǫp ,
where G and K are the elastic shear and bulk modulus respectively, α s the thermal dilation coefficient of the solid skeleton and T the temperature. The rate of plastic deformation is written as

γp = 1 H [ τ + µ( σ + ṗ)] ; ǫp = β γp , (7)
where H = H(γ p ) = h(σ + p) (with h = dµ/dγ p ) is the plastic hardening modulus, which is related to the tangent modulus H tan of the τ versus γ curve through the relationship H tan = H 1+H/G and is either positive or negative according to whether the τ -γ curve is increasing (hardening) or decreasing (softening). Notice that, in principle H tan is varying with strain, as thoroughly presented by [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF]. When H tan is non positive then the material is enterring an unstable regime, where softening induces localization of deformation in a shear band [START_REF] Holcomb | Inelastic constitutive properties and shear localization in Tennessee marble[END_REF]).

However, when the material is fluid saturated [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]), thermal pressurization may cause localization of deformation even at hardening regimes, before reaching H tan = 0, where the material is supposed to be stable. In this work we study the effect of a chemical reaction at high temperatures, near the activation temperature of the reaction, to assess whether the reaction itself may induce softening and localization at even higher values of the hardening modulus.

Stress Equilibrium

For the considered 2D Cosserat continuum and the configuration of Fig.

2, the local equilibrium equations are

∂σ 12 ∂x 2 = 0, ( 8
)
∂σ 22 ∂x 2 = 0, ∂m ∂x 2 + σ 21 -σ 12 = 0,
where m = m 2 and m 1 = 0. The stresses applied at the boundary are a shear stress τ and a normal stress σ n in the x 2 direction. It is assumed that no couple stress is imposed at the boundary. Prior to localization, the state of stress and strain is uniform and considering the couple free boundary condition, the couple stress is identically zero in the sheared layer. Therefore the rock behaves as a classical continuum before localization.

Reaction kinetics

We assume the presence of a decomposition reaction that decomposes the solid skeleton and produces excess pore fluid. In particular, as an example we consider the reaction of calcite decomposition (CaCO 3 → CaO + CO 2 ), producing CO 2 which at high pressures and temperatures (conditions similar to the ones met in faults at depth) is at supercritical, liquid state [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF]). Evidence of the presence of this reaction in faults is found in many active crustal faults [START_REF] Famin | Earthquakes produce carbon dioxide in crustal faults[END_REF]), and experimentally reproduced for example by [START_REF] Han | Ultralow friction of carbonate faults caused by thermal decomposition[END_REF] from experiments on simulated faults in Carrara marble. Other types of chemical reactions such as mineral dehydration (Brantut et al. (2011)) could be studied following the same methodology.

Following [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF] and [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF], we assume that the reaction is taking place with a rate r, which on its turn can be expressed from first order reaction kinetics as [START_REF] Law | Combustion Physics[END_REF])

r ≈ ρ s M CaCO 3 A 0 e -Tc/T , (9)
where T c = E/R a the activation temperature of the reaction, E the activation energy and R a the universal gas constant. In this expression A 0 is the preexponential factor, a parameter that may vary significantly, obtaining for calcite decomposition values between 10 2 s -1 and 10 20 s -1 (Lvov et al. ( 2002)).

The Arrhenius model assumed in Eq. ( 9), is the simplest one, but has received much criticism for its validity, mainly due to the strongly varying magnitudes of both A 0 and T c for a given reaction (Lvov et al. (2002)). In the pertinent literature there are numerous models suggested to replace it [START_REF] Law | Combustion Physics[END_REF] The internal plastic deformation induced by milling during frictional sliding might further reduce the activation energy of the reaction (Fisher. (1988)).

Hence, emphasis should be given on thorough experimental detemination of the reaction parameters at various temperature and pressure regimes.

Fluid-mass Conservation

Following the decomposition of the plastic strain rates into elastic and plastic parts (Eq. 4), we consider the volumetric plastic strain rate to consist of a mechanical part ǫp m subject to the constitutive behavior of Eq. 7 and a chemical part ǫp c , influenced by the reaction rate r

ǫp = ǫp m + ǫp c . ( 10
)
Under the above considerations and by assuming that the solid matrix is plastically incompressible, the diffusion-reaction equation for the pore pressure is obtained from the fluid mass equation as [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF])

∂p ∂t = c hy ∂ 2 p ∂x 2 2 + Λ ∂T ∂t - β β ⋆ ∂ γp ∂t + ρ s -ρ f ζ ρ f β ⋆ A 0 M CaCO 3 e -Tc/T , (11)
where

ζ = M CaCO 3 M CO 2 - ρ CaCO 3 ρ CaO M CaO M CO 2
, ρ i the density of the i-th constituent and ρ s the density of the solid skeleton, c hy the hydraulic diffusivity, β ⋆ the storage capacity and Λ = (λ f -λ n )/β ⋆ the pressurization coefficient, expressed as the difference of the thermal expansion coefficimobilizedents of the fluid (λ f ) and the pore volume (λ n ).

Heat Equation

By assuming that all the mechanical work input is converted into heat [START_REF] Chester | Fracture surface energy of the Punchbowl Fault, San Andreas System[END_REF]; Pittarello et al. ( 2008)), and that the heat is expressed through Fourier's law, we may obtain the diffusion-reaction temperature equation [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF]; [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF])

∂T ∂t = c th ∂ 2 T ∂x 2 2 + 1 ρC τ γp - |∆H| ρC ρ s M CaCO 3 A 0 e -Tc/T , (12)
where c th is the thermal diffusivity, ρC the specific heat capacity of the mixture and |∆H| ≈ E the specific enthalpy of the reaction, expressing the energy consumed during an endothermic reaction.

Undrained-adiabatic limit

The first approximation in our effort to assess the influence of these weakening mechanisms during earthquake slip is to prohibit drainage and heat flux at the boundaries of the layer (q f = q h = 0). This case corresponds to the well known undrained-adiabatic limit, acknowledged to be established at or near seismic rates, i.e. when the heat production rate and pore pressure increase in the slipping zone is significantly larger than the heat and pore pressure diffusion times towards the fault core and in the damage zone (Gara- 2011)).

In addition to the above, we assume that the normal stress σ n acting on the sheared layer is constant. From Eqs. ( 7),( 11),( 12) we may obtain the following relationship between the shear stress rate, the strain rate and the temperature rate:

τ = H 1 + H/G γ + B 1 + H/G Ṫ , (13)
where

H = H + µ β β ⋆ -µτ 1 β ⋆ |∆H| ρ s -ζρ f ρ f , (14)
B = µ ρC β ⋆ |∆H| ρ s -ζρ f ρ f -Λ .
Obviously, when the thermal effects are neglected (i.e. in the absence of thermal pressurization and chemical reaction) and incompressible solid and fluid are assumed (so that β ⋆ = 1/K) the classical expression for dilatant hardening effect, as proposed by [START_REF] Rice | On the stability of dilatant hardening for saturated rock masses[END_REF] is retrieved

τ = H + µβK 1 + (H + µβK)/G γ. ( 15
)
For dilatant material (β > 0) we may recognize in (13) a hardening effect due to dilatancy with the term µβ/β ⋆ and a strong softening effect from the reaction due to the term µτ

β ⋆ ρs-ζρ f ρ f |∆H|
(positive for calcite decomposition, as

shown by [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF]).

As shown by [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF], in the absence of a chemical reaction instability may occur even in the hardening regime of the τ -γ curve (H > 0)

if the dilatancy effect is very weak (as it is the case at great depth) and if the thermal pressurization effect is significant. Here we verify this claim and indicate that the presence of a chemical reaction may lead to instability at even higher hardening regimes.

Linear Stability Analysis of undrained adiabatic shearing of a

Cosserat layer

At any given state, the evolution time of the kinematic quantities, the pore fluid pressure and the temperature determine the stability of the shearing of the Cosserat layer. The small perturbations considered herein are defined as follows:

u 1 (x 2 , t) = u 0 1 (x 2 , t) + ũ1 (x 2 , t), (16)
u 2 (x 2 , t) = u 0 2 (x 2 , t) + ũ2 (x 2 , t), ω c (x 2 , t) = ω 0 c (x 2 , t) + ωc (x 2 , t), p(x 2 , t) = p 0 (x 2 , t) + p(x 2 , t), T (x 2 , t) = T 0 (x 2 , t) + T (x 2 , t).
The governing equations for the perturbed fields (ũ 1 , ũ2 , ωc , p , T ) are of the same form as the ones derived in the previous section. Similarly, the fluid and heat fluxes are prevented at the boundaries of the layer. However, the internal fluid and heat flows are permitted inside the layer.

Normalized equations

The system of equations can be brought in a dimensionless form by introducing the dimensionless quantities:

x = x 1 R , z = x 2 R , ūi = u i R , ( 17
) p = p σ n , σij = σ ij σ n , t = c th R 2 t , T = Λ σ n T,
where R is the internal length of the Cosserat model, which can be related to the mean grain size R = d 50 /2. The linearized forms of the pore pressure and heat equation (linearized around T 0 = T c , i.e. at higher temperature regimes, near the activation temperature of the reaction) may be reduced in their dimensionless forms as (the superimposed bars are omitted for convenience):

∂p ∂t = η ∂ 2 p ∂z 2 + ∂T ∂t + δ ∂ γp ∂t + χ T -T c T c , (18)
∂T ∂t = ∂ 2 T ∂z 2 + α ∂ γp ∂t -ψ T -T c T c , where η = c hy c th , δ = β β ⋆ σ n , α = Λτ 0 ρCσ n (19) χ = ρ s -ζρ f ρ f β ⋆ R 2 Λc th A 0 ψ = ρ s ∆H ρCM CaCO 3 R 2 c th Λ σ n A 0 , and
∂γ p ∂t = 1 H ⋆ ∂σ 12 ∂t + σ 21 ∂t + µ σ ′ 11 ∂t + σ ′ 22 ∂t , (20)
H ⋆ = H σ n = h σ ′ σ n .
The corresponding dimensionless form of the linear and angular momentum with inertia and micro-inertia terms is as follows:

∂σ 12 ∂z -I ∂ 2 u 1 ∂t 2 = 0, ( 21
) ∂σ ′ 22 ∂z - ∂p ∂z -I ∂ 2 u 2 ∂t 2 = 0, ∂m ∂z + σ 21 -σ 12 - I 2 ∂ 2 ω c ∂t 2 = 0,
where I =

Numerical example

We consider the case of a 7-km deep fault, of a calcite-made gouge material having the parameters depicted in Table 4.2. With these values, and by performing the method presented in the Appendix A, we obtain a critical value for the hardening modulus h cr ≈ 6, such that for h ≥ h cr all the roots of the polynomial equation have negative real part, and thus the system is stable. For h < h cr roots with positive real parts exist, bringing the system to instability.

We observe that the destabilizing effect of the chemical reaction by the fact that instability can occur in the hardening regime, for values of the hardening modulus significantly larger than in the case of thermal pressurization, h ′ cr = 0.014 [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]). Hence, as temperature increases reaching the activation temperature of the reaction, the material becomes unstable even at regimes mechanically considered as stable in the absence of chemical reaction.

In order to assess the localizing effect of the chemical reaction, we plot in the Appendix A the greatest root (i.e. the one corresponding to the insta- Parameter Value

Units

Parameter Value Units bility mode with the fastest growth in time) as a function of the wavelength λ (Figure 6b) and select the wavelength number λ = λ m corresponding to the maximum value of the growth coeffcient s, so that a wavelength with the fastest growth in time appears. In Figure 3 this selected wavelength is plotted against the hardening modulus. As in [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF] the selected wavelength tends to infinity for h → h - cr denoting that at initial instability no localization takes place. As the hardening modulus decreases the wavelength decreases as well, reaching a minimum at h ≈ 0.16, where λ m ≈ 3.95. [START_REF] Muhlhaus | Thickness of shear bands in granular materials[END_REF] similarly obtained that at the bifurcation state the shear-band thickness for dry granular materials is infinite and decreases beyond the bifurcation state as the strain localization process evolves.

φ 0 0.03 - c hy 10 -5 m 2 /s σ ′ n 200 MP a c th 10 -6 m 2 /s M CaCO 3 0.1 kg/mol R 0.1 mm M CaO 0.056 kg/mol G 10 4 MP a M CO 2 0.044 kg/mol K 2 • 10 4 MP a ρ s 2.5 • 10 3 kg/m 3 µ 0.5 - ρ f 10 3 kg/m 3 β 0 - ρC 2.8 MP a/ • C A 0 10 6 1/s Λ 0.5 MP a/ • C |∆H| 200 kJ/mol
We notice that the selected wavelength (λ m ≈ 3.95) in this case is significantly smaller than in the absence of the chemical reaction (λ ′ m ≈ 187), revealing the crucial, localizing role of the chemical reaction. Since for this example we have used indicative values for the reaction parameters, like the pre-exponential factor and the activation temperature, we proceed with an analysis of these parameters in order to reveal their influence on the width of the PSZ and the type of instability.

Parametric Analysis. Impact of the chemical reaction and the pressurization mechanism

From the parameters appearing in 4.2, those being less constrained in values are the parameters of the reaction A 0 and T c . Since they always appear in the parameter groups as a ratio A 0 /T c , we emphasize on the influence of this ratio on the stability of the system by treating it as a parameter.

Effect of the chemical reaction

To assess the influence of the chemical reaction, we plot in Figure 4 the values of h cr with the ratio of the reaction parameters A 0 Tc varying from 10 -4 to 10 12 . We observe that as A 0 Tc → 0, then h cr → h ′ cr = 0.014, thus it tends to the case of pressurization alone as weakening mechanism [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]). As A 0 Tc increases, h cr increases as well, reaching an upper threshold at h m cr ≈ 6.

A 0 /T c ).
It is to be noted that as soon as h cr reaches h m cr , then the maximum positive root changes from purely real to complex (Figure 4). This effect takes place at a critical A 0 Tc (in this case A cr 0 = 10 7 s -1 for T c = 800 • C) and signifies the onset of oscillatory instabilities due to the chemical reaction.

The existence of oscillatory instabilities is questionable because when the temperature decreases below the given threshold the reaction stops. In other words, the linear stability analysis is not valid for non monotonous processes in case of non-linear behaviour. For this reason the exact type of these oscillatory instabilities, as well as their time evolution cannot be determined by the present linear stability analysis. To this end nonlinear analyses and time integration of the full system of equations is required, to determine whether imaginary eigenvalues may correspond to stick-slip type of instability (appearing as limit cycles on a phase diagram, as for example in Alevizos and Veveakis (2011))

Effect of thermal pressurization mechanism

As already mentioned and anticipated, when A 0 Tc → 0, then h cr → h ′ cr = 0.014, thus the system is influenced by pressurization alone as weakening mechanism [START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]). However, when the reaction becomes strong and A 0 Tc increases so that h cr → h m cr , the effect of thermal pressurization, compared to the chemical effect (i.e. the production of excess pore fluid due to CaCO 3 lattice breakdown and CO 2 production), is negligible. This fact becomes apparent in the Appendix A, where we compare the s -λ curves of Figures 7b and6b with the corresponding ones in the absence of thermal pressurization, appearing in Figure 8.

We notice that the same response is obtained, both for the magnitudes of h, s and λ and for the type of instability (nodal instability for A 0 = 10 6 s -1 and oscillatory for A 0 = 10 7 s -1 , T c = 800 • C). The pre-exponential factor for the considered depths cannot be realistically estimated, since it is pressure dependent (Lvov et al. (2002)) and there is a lack of experimental data at high pressures, however usually it is accepted to be more than 10 10 s -1

(for example [START_REF] Dollimore | The kinetic interpretation of the decomposition of calcium carbonate by use of relationships other than the Arrhenius equation[END_REF] suggest a value of the order of 10 15 s -1 , used also by [START_REF] Sulem | Thermal decomposition of carbonates in fault zones: slip-weakening and temperature-limiting effects[END_REF]). Thus, we may conclude that at high temperatures where reaction is triggered (around T c ), and for the decomposition reaction considered, the main weakening mechanism of the two is the reaction, verifying the recent experimental findings of [START_REF] Han | Ultralow friction of carbonate faults caused by thermal decomposition[END_REF] and and localization is driven by thermal pressurization, a fact that could not be analyzed in this study due to the choice of the reference temperature T c , around which we linearized in Section 4.1, in order to emphasize on the reaction itself. As temperature increases the process is progressively localizing towards an ultimate shear-zone whose thickness is determined by the reaction itself. To address this claim, we proceed with a shear band thickness analysis.

Shear band thickness, Wavelength selection

Indeed, it is expected that at temperatures near T c the width D of the PSZ of Fig. 1 is influenced by the chemical reaction. To verify this claim we may plot the minimum selected wavelength against the ratio A 0 /T c (Figure 5) to conclude that with increasing A 0 /T c , λ m and thus D decrease. In the absence of chemical reaction, λ m was calculated to be equal to about 187.

The ultimate value for λ m in the presence of a chemical reaction is around 3, denoting that the PSZ may be an extremely localized structure, obtaining λ m = 300 µm when d 50 = 0.1 mm as is the case for fine sands, or even less if we consider possible average grain sizes from principal slipping zones, up to 10 µm, in the lines of the field evidence provided by [START_REF] Smith | Principal slip zones in limestone: Microstructural characterization and implications for the seismic cycle (Tre Monti fault, Central Apennines[END_REF] and thoroughly discussed in the introduction.

Thus, indeed the onset of the reaction is taking place in an ultra localized zone, that could be significantly thinner than the initial one defined by thermal pressurization alone (λ m ∼ 187). This means that in a "nonreactive"fault thermal pressurization of pore fluids would result in a slipping zone of 18.7 mm thick (for a fault material with d 50 = 0.1 mm), while under the same ambient conditions and slip rate, in a "reactive" gouge (e.g., calcite) chemical reaction would result in a slipping zone of 0.3 mm . What have to be noted at this point, is that localizing mechanisms like flash heating and breakage have been neglected in this study, as discussed thoroughly in Sections 2 and 5.2. Including these mechanisms would provide a more realistic behaviour of the fault even at low temperatures where pressurization or the chemical reaction are absent. In addition, the thickness of this ultimate structure (the PSZ) seems to push continuum theories to its limits, since the width of the structure is comparable to the size of the average grain. However, we emphasize that in this study the average grain size used (d 50 ) is the one of the initial grain size distribution and that during the evolution of a chemical reaction, grains undergo chemical degradation and phase transitions that could significantly alter the average grain size, or even make the determination of grain size ambiguous if amorphization [START_REF] Yund | Amorphous material in high strain experimental fault gouges[END_REF]; [START_REF] Brantut | High velocity frictional properties of a clay-bearing fault gouge and im-plications for earthquake mechanics[END_REF]) is the chemical reaction. As shown by [START_REF] Veveakis | Chemical reaction capping of thermal instabilities during shear of frictional faults[END_REF]; [START_REF] Alevizos | Thermal and chemical mechanisms of frictional faults: Static and dynamic instabilities[END_REF] during shearing of faults we may identify a regime of lower temperatures, where mechanical effects determine the evolution of the system and a high temperature regime, where the reaction is triggered and dominates the response of the fault, irrespective of the mechanical behavior. Thus, the present Cosserat approach is valid up to the point that chemical reaction is fully set and determines all the mechanical effects of the fault, like grain size evolutions and rheological response. It is used in order to obtain an estimate of the thickness of the layer at which chemical reaction will be triggered. Pass this point, either discrete analyses or even higher order continuum theories [START_REF] Stefanou | Homogenization of interlocking masonry structures using a generalized differential expansion technique[END_REF]) must be applied to account for the chemical reaction effects (for example by introducing varying d 50 ).

Discussion

Thermal instabilities were shown to be important at elevated temperatures, and to determine the stability of a fault even at regimes that are considered to be stable under purely mechanical conditions. The formalism presented here however revealed the importance of parameters that are not well-constrained. In particular, the hardening modulus h of rock materials is in principle varying with shear strain, while the parameters of the reaction A 0 and T c strongly depend on the applied pressures at which the reaction takes place (Lvov (2007), pp. 79) and several values can be found in the literature (For example A 0 may vary from 10 2 s -1 to 10 20 s -1 for calcite decomposition, as discused by Lvov et al. ( 2002)).

Therefore, in order to apply the presented theory in real faults, there is a need of extensive laboratory testing of the strain-stress response of the fault gouge materials and on the effect of pressure and temperature on this response.

The results presented in this work could provide insight for the thermomechano-chemical coupling occurring in soil materials. Given the stressstrain response of the material, the ambient thermal conditions of the fault determine whether the fault may weaken thermally or not. At lower temperatures thermal pressurization may be triggered and cause softening at positive, but near zero, hardening modulus. At higher temperatures, near the activation of chemical reactions thermal softening induces a strong weakening effect since it may take place at extremely high (positive) values of the hardening modulus. Therefore at high temperatures even if the material is mechanically stable (steep hardening regime) the triggering of a chemical reaction would lead the system to instability and softening.

The model presented here contains parameters that need to be evaluated from interdisciplinary experiments at conditions that match the ones met in real faults at seismic depths. Parameters like those appearing in the reaction rate of the chemical reaction (A 0 , T c ) or in friction laws of geomaterials that would provide parameters like G c should be evaluated before applying such models in real cases.

Conclusions

In this work, we have used higher order continua considerations and allowed to the fault material additional, rotational degree of freedoms, along with the presence of a decomposition reaction. We have shown that the illposed problem of simple shear of a saturated fault can be remedied when accounting for the micro-inertia of the rotations and that, at temperatures near its activation energy, the impact of the reaction is significant to the stability of the mechanical problem and the localization of deformation.

We have shown that the width of the localizing shear zone at which dynamic weakening of faults take place depends strongly on the parameters of the reaction, i.e. the pre-exponential factor A 0 and the activation temperature T c . Small values of the ratio A 0 Tc correspond to broader zones, of the order of few centimeters.

A. Linear Stability Analysis

By substituting the perturbation fields ( 22) into the linearized governing equations (18, 21), we obtain a homogeneous algebraic system for the coefficients U 1 ,U 2 ,Ω,Π,Θ:

           sδσ(d+ǫ+2cµ) 2hλσ ′ sδσ(2f +(a+b)µ) 2hλσ ′ 0 -aαssTcδµσ+αsbsTcδµσ+2hσ ′ (sTc+χ) 2hTcσ ′ s + η λ 2 -sασ(d+ǫ+2cµ) 2hλσ ′ sασ(2f +(a+b)µ) 2hλσ ′ 0 1 λ 2 + s(aαsαµ]σ+αsbαµ]σ+2hσ ′ ) 2hσ ′ + ψ Tc 0 -s 2 I -d λ 2 -f λ 2 -d+ǫ λ 0 0 -c λ 2 -s 2 I -a λ 2 0 αs(a+b) 2λ 1 λ -d+ǫ λ 0 -2d + 2ǫ -s 2 I 2 -g λ 2 0 0                       U 1 U 2 Ω Π Θ            = 0, (23)
where where only negative roots appear, rendering the system stable, (b) for h < h cr . The dashed curves represent the real part of complex roots whereas the solid curves depict the purely real roots (zero imaginary part). The wavelength at which the growth coefficient s obtains its maximum value is the one corresponding to the selected wavelength λ m .

a = G + K -(µβK 2 )/(G + H + Kβµ) σ , (24
of the characteristic polynomial. To this end we consider the case of a 7km deep fault, of a gouge material having the parameters depicted in Table 4.2. With these values we obtain a critical value for the hardening modulus h cr ≈ 6, such that for h ≥ h cr all the roots of the polynomial equation have negative real part, and thus the system is stable (Figure 6a). For h < h cr roots with positive real parts exist, bringing the system to instability (Figure 6b).

In order to asses the localizing effect of the chemical reaction, we may plot the greatest root (i.e. the one corresponding to the instability mode with the fastest growth in time) as a function of the wavelength λ. This curve exhibits a maximum at λ = λ m so that a wavelength with the fastest growth in time appears (Figure 6b)

In the absence of thermal pressurization (i.e. when Λ = 0 in the original The maximum root is complex in this case, denoting the onset of oscillatory instability.

Figure 1 :

 1 Figure 1: Idealized description of a natural fault zone, as presented by[START_REF] Chester | Internal structure and weakening mechanisms of the san andreas fault[END_REF]

  ; Ferri et al. (2010); Brantut et al. (2011a); De Paola et al. (2011)).

Figure 2 :

 2 Figure 2: Problem formulation. Cosserat stresses, couple stresses and degrees of freedom.

  gashand Rudnicki. (2003a,b)). It provides a first insight on the impact of each mechanism on the dynamic weakening of faults. The actual value of heating rate required for the onset of the weakening mechanisms varies with the material of the fault gouge(Di Toro et al. (

  and the expressions of the stresses and couple stresses are derived from the elastoplastic model presented in Section 2.1[START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]).FollowingVardoulakis and Sulem (1995) and[START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF], we introduce the perturbed fields (16) for a given initial configuration, and the spatial dependency of the perturbations is decomposed into Fourier modes with wavelength λ. Hence we introduce perturbations of the form ũ1 (x 2 , t) = U 1 e st sin z λ = D/(πNR) and N equal to an integer satisfying the zero heat and fluid flux boundary conditions at z = x 2 R = 0. By substituting the perturbation fields (22) into the linearized governing equations (18, 21), we obtain a homogeneous algebraic system for the coefficients U 1 ,U 2 ,Ω,Π,Θ. A nonzero solution is possible only when the determinant of the linear system vanishes. This request results in a characteristic polynomial equation of degree 8 for the growth coefficient s. If a root has a positive real part, then the corresponding perturbation grows exponentially in time.

Figure 3 :

 3 Figure 3: (a) Wave length λ corresponding to the fastest growth coefficient of the instability versus the hardening modulus h. (b) The same as in (a), focused in the area 0 < h < 1, where the selected (minimum) wavelength, λ m , is located. The results are calculated for A 0 = 10 6 .

Table 1 :

 1 Indicative material parameters for a fault at 7 km depth, where the initial temperature is about 200 • C. For these values of the parameters and for T c = 800 • C the corresponding dimensionles groups are ζ = 1.24, α = 0.06, δ = 0, η = 10, I = 1.25 • 10 -9 , ᾱs = 10 -2 , χ ≈ 2A 0 = 2 • 10 6 , ψ ≈ 0.01A 0 = 10 4

Figure 4 :

 4 Figure 4: The critical hardening modulus h cr as a function of A 0 /T c (logarithmic scale for

Ferri

  et al. (2010). However, as also discussed by[START_REF] Han | Strong velocity weakening and powder lubrication of simulated carbonate faults at seismic slip rate[END_REF] and DePaola et al. (2011) important weakening mechanisms like powder lubrication and flash heating (see also[START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]) should not be excluded if one would like to model all the mechano-physical processes taking place in a fault, since they could have equal or even more dramatic role than chemical pressurization. In addition secondary chemo-mechanical effects like the breakdown of calcite grains expressed through the production of holes in the Calcite grains due to decarbonation reduce also the compressive strength, forcing the microstructure to collapse. In this study we have excluded these mechanisms to emphasize on the two pressurization mechanisms, i.e. thermal and chemical pressurization.However, we notice that the initial temperature for a 7-km fault is around 200 • C, far below the activation temperature T c ≈ 800 • C of calcite decomposition. Up to the temperature regime around T c the reaction term is negligible

Figure 5 :

 5 Figure 5: The selected wavelength λ m as a function of A 0 /T c

Figure 6 :

 6 Figure 6: Real parts of the roots of the characteristic equation (growth coefficient s) as a function of the wavelength number λ, for A 0 = 10 6 and T c = 800 (a) for h ≥ h cr = 6

Figure 7 :

 7 Figure 7: Real parts of the roots of the characteristic equation (growth coefficient s) as a function of the wavelength number λ, for A 0 = 10 7 and T c = 800 (a) for h ≥ h cr = 6 where only negative roots appear, (b) for h < h cr . The dashed curves represent the real part of complex roots, whereas the solid curves depict the purely real roots (zero imaginary part).

Figure 8 :

 8 Figure 8: Real parts of the roots of the characteristic equation (growth coefficient s) as a function of the wavelength number λ, in the absence of thermal pressurization for T c = 800 • C and (a) A 0 = 10 6 s -1 , (b) A 0 = 10 7 s -1 .

Acknowledgments

The authors would like to thank G. Di Toro, N. Brantut and an anonymous reviewer for their comprehensive reviews and fruitfull comments, from which this paper has greatly benefited. The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n 228051.


Sandra Kim Tiam
email: sandra.kim-tiam@cemagref.fr

Agnès Feurtet-Mazel

François Delmas

Nicolas Mazzella

Soizic Morin

Guillemine Daffe

Patrice Gonzalez

Development of q-PCR approaches to assess water quality: Effects of cadmium on gene expression of the diatom Eolimna minima

Keywords: Eolimna minima, Diatom, Cadmium, Metal, Quantitative real time PCR

This study was undertaken to develop molecular tools to assess water quality using diatoms as the biological model. Molecular approaches were designed following the development of a rapid and easy RNA extraction method suited to diatoms and the sequencing of genes involved in mitochondrial and photosystem metabolism. Secondly the impact of cadmium was evaluated at the genetic level by q-PCR on 9 genes of interest after exposure of Eolimna minima diatom populations cultured in suspension under controlled laboratory conditions.

Their growth kinetics and Cd bioaccumulation were followed.

Population growth rates revealed the high impact of Cd at 100 μg/L with total inhibition of growth. These results are linked to the high bioaccumulation values calculated after 14 days of exposure, 57.0 ± 6.3 μgCd/g dw and 734.1 ± 70 μgCd/g dw for exposures of 10 and 100 µgCd/L respectively. Genetic responses revealed the impact of Cd on the mitochondrial metabolism and the chloroplast photosystem of E. minima exposed to 10 and 100 µgCd/L with induction of cox1, 12S, d1 and psaA after 7 days of exposure for the concentration of 100 µgCd/L and of nad5, d1 and psaA after 14 days of exposure for both conditions. This is the first reported use of q-PCR for the assessment of toxic pollution on benthic river diatoms. The results obtained presage interesting perspectives, but the techniques developed need to be optimized before the design of new water quality diagnosis tools for use on natural biofilms.

Over the past decades, the ever increasing release of agricultural, industrial and domestic waste led to a significant contamination of the environment and particularly of the aquatic compartment.

In regard to the resulting degradation of freshwater quality, various diagnostic tools have been designed. In France, in collaboration with Water Agencies and Environment Ministry services, CEMAGREF has developed several diatom indices to estimate global water quality like the "Indice de Polluosensibilité Spécifique" (IPS, Cemagref 1982) or the Biological Diatom Index (BDI, Coste et al. 2009), both routinely used for monitoring applications in several European countries. Nevertheless, the indices currently used for general water quality assessment -including the IPS and BDI -were not designed to assess specific toxic alterations (metals, synthetic organic pollutants) and are not really suitable for a sensitive diagnosis of this type of pollution. Thus, heavy metals are of particular concern because of the different toxic effects they can produce [START_REF] Sauvant | Toxicity assessment of 16 inorganic environmental pollutants by six bioassays[END_REF][START_REF] Bucio | Cadmium and mercury toxicity in a human fetal hepatic cell line (WRL-68 cells)[END_REF].

Recent works have studied the responses of diatom communities to metal pollution, for example Morin et al. (2008a) observed structural impact at the community level and morphological abnormalities in a metal-polluted stream. Other studies assessed effects of metal-induced oxidative stress on functional descriptors such as photosynthetic [START_REF] Antal | Fluorescence study of the effect of mercuric chloride and methylmercury chloride on the photosynthetic activity of the diatom Thalassiosira weissflogii[END_REF] or antioxidant [START_REF] Branco | Sensitivity of biochemical markers to evaluate cadmium stress in the freshwater diatom Nitzschia palea (Kützing) W. Smith[END_REF]) enzyme activities.

However, studies assessing genetic responses of diatoms to contamination exposure are still at the early stage principally because of the lack of available diatom nucleotide sequences in genomic databases and difficulties to access genetic material. Only Thalassiosira pseudonana and Phaeodactylum tricornutum -two marine diatom species -have as yet been entirely sequenced, and genomic information about freshwater species is still extremely scarce.

Extracting the genetic material of organisms is the first step to obtain sequences of interest.

Diatoms require a particular approach for this because of their unique external structure.

Diatoms are unicellular algae and their most obvious distinguishing characteristic is their siliceous cell wall called the frustule. Various approaches have been used to break the frustule, but these are often time consuming [START_REF] Wawrik | Real-time PCR quantification of rbcL (ribulose-1,5-bisphosphate carboxylase/oxygenase) mRNA in diatoms and pelagophytes[END_REF] and/or need specialized laboratory material such as a French press [START_REF] Stabile | Molecular analysis of intraspecific variation in the marine diatom Skeletonema costatum[END_REF][START_REF] Davis | Characterization of a modular, cell-surface protein and identification of a new gene family in the diatom Thalassiosira pseudonana[END_REF][START_REF] Hildebrand | Characterization of a silicon transporter gene family in Cylindrotheca fusiformis: Sequences, expression analysis, and identification of homologs in other diatoms[END_REF] or a MiniBeadBeater [START_REF] Fawley | A simple and rapid technique for the isolation of DNA from microalgae[END_REF]). More rapid and straightforward extraction methods are needed, specifically for diatoms. The expression of many genes is involved in the response of organisms to toxicants and can be disturbed by them. Once inside the cell, metals can cause oxidative stress [START_REF] Wang | Cadmium inhibits the electron transfer chain and induces Reactive Oxygen Species[END_REF] affecting the way in which the mitochondria [START_REF] Stohs | Oxidative mechanisms in the toxicity of metal ions[END_REF] or photosynthesis function [START_REF] Knauert | The role of reactive oxygen species in copper toxicity to two freshwater green algae[END_REF]. The mitochondrial superoxide dismutase (sodMn) gene was selected because of its involvement in antioxidant defences. The mitochondrial metabolism was investigated using the cytochrome C oxidase (cox1) subunit and the NADH dehydrogenase subunit 5 (nad5) genes. In addition, the quantity of mitochondria in the cells was estimated using mitochondrial 12S ribosomal RNA. D1 protein (d1) and PsaA protein (psaA) are components of PS2 and PS1 respectively and were selected to investigate photosystem metabolism. Cytochrome P450 1A1 (cyp1A1) was selected as a biomarker for exposure to polynuclear aromatic hydrocarbons. Two genes commonly used as references in numerous genetic studies were selected as for study in our work: β actin (act) and 18S ribosomal RNA.

Consequently, the first objective of this work was to develop a diatom-specific extraction method and to shape molecular tools to select and sequence genes of interest to assess metal contamination. We used cultures of Eolimna minima, a metal-tolerant freshwater diatom commonly collected within periphytic biofilm samples in running water and especially in metal-contaminated areas (Morin et al, in press). Our attention was particularly focussed on genes involved in mitochondrial and photosystem metabolism and responses of E. minima to a gradient of cadmium contamination was then characterised using these selected molecular tools. Cultures of E. minima were exposed in the laboratory for 14 days to 10 µgCd/L or 100 µgCd/L. The moderate exposure corresponds to the concentrations reported in cadmium polluted rivers like the river Riou-Mort (South West France) [START_REF] Feurtet-Mazel | Study of periphytic diatoms communities exposed to metallic contamination through complementary field and laboratory experiments[END_REF]) and the highest is comparable to those reported in more highly polluted European rivers [START_REF] Ivorra | Translocation of microbenthic algal assemblages used for in situ analysis of metal pollution in rivers[END_REF]. Cell numeration, bioaccumulation and genetic responses were followed 1, 2, 7 and 14 days after contamination.

2.Materials and methods

2.1.Molecular methodological developments

2.1.1.RNA extraction method

In order to access the genetic material of diatoms, a novel RNA extraction method was developed. Designing an efficient method able to break the external silica wall of diatoms without degrading the genetic material was a prerequisite. For this purpose, cultures of the diatom Eolimna minima were grown in batches from axenic strains provided by the UTEX algal collection (Texas University, USA). The diatom cells were cultured in sterile Dauta medium [START_REF] Dauta | Conditions de développement du phytopancton. Etude comparative du comportement de huit espèces en culture. I. Détermination des paramètres de croissance en fonction de la lumière et de la température[END_REF] with silica added to a final concentration of 10mg/L in 3L Erlenmeyer flasks. The cultures were maintained at between 17 and 18°C in a thermostatic room with a photon flux density of 160µmol.m -2 .s -1 and a 12:12 h light:dark cycle. Fresh cultures were inoculated every seven days to promote optimal growth of the diatoms until sufficient cells were obtained. 30 mL (2.10 6 cells/mL) were removed from the E. minima cultures and the cells harvested after centrifugation at 3 863g for 5 min at room temperature. The cell pellet was transferred to 2 mL Eppendorf tubes and resuspended in 1 mL of diethyl pyrocarbonate (DEPC) treated water. Then, the diatoms were centrifuged again at 2 419g for 5 min at room temperature and the supernatant was removed. 1 mL of Trizol (invitrogen) and 300 μL of glass beads (0.10 -0.11mm-diameter, B. Braun Biotech International) were added to the cell pellet. Then the sample was vortexed 3 times for 30 seconds, and the supernatant transferred into a clean Eppendorf tube and placed in a water bath at 30°C for 5 min. After addition of 200 µL of Chloro RECTAPUR (VWR), the sample was vortexed for 10 s and placed in water bath at 30°C for 2 min. The tubes were centrifuged at 11 360g for 5 min at room temperature and the aqueous phase containing the genetic material was transferred into clean tubes. The end of the extraction was realized using Absolultely RNA Miniprep Kit (Stratagène) according to the manufacturer´s instructions: 550 µL of 75% ethanol were added to the aqueous phase, the sample was vortexed, transferred to affinity columns and centrifuged for 1 min at 11 360g at room temperature.

Filtrates were removed, 600 µL of low-salt buffer were placed on the column which was centrifuged at 11 360g at room temperature for 2 min. 5 µL of DNAse 1 [1U] and 50 µL of activity buffer were added to the column and placed in a water bath at 37°C for 15 minutes. 600 μL of high-salt buffer were added to the column which was centrifuged for 1 min at 14 500 rpm. 300 µL of low-salt buffer were added and the column was centrifuged for 2 min at 14 500 rpm. The columns were transferred into clean tubes and 30 µL of elution buffer at 60°C were added and centrifuged for 1 min at 12 250g at room temperature. Total RNA was analysed on 1% (w/v) agarose gel with ethidium bromide by UV light on an illuminator.

2.1.2.Reverse transcription of RNA

The first strand of cDNA was synthesised from 14 µL of total RNA (3 to 5 µg) using the Stratascript first strand synthesis system (Agilent). After the addition of 1 µL of oligo(dT)

[1µM], 1 µL of random primers [1µM], 0.8 µL of dNTPs [10mM] and 2 µL of 10x firststrand buffer the reaction was incubated for 5 min at 65°C. Then 1 µL of Stratascript reverse

Author-produced version of the article published in Water Research, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005 transcriptase [1U/µL] and 0.5 µL of RNase inhibitor [0.5U] were added, the reaction was incubated for 1 h at 42°C in an Eppendorf Mastercycler. The cDNA mixture was conserved at -20°C until it was used in a real-time PCR reaction.

2.1.3.Cloning and molecular characterization of the target genes

Genetic research focussed on seven genes involved in responses to environmental contamination and/or in which expression can possibly be disturbed by environmental factors (sodMn,nad5,d1,cox 1,psaA,cyp 1A1,and 12S) and for two reference genes (18S and act).

Primers were designed by performing ClustalW analysis from protein or nucleic sequences (12S and 18S) from marine diatoms and other phylogenetically related aquatic organisms available in NCBI databases. From these alignments, primer pairs about 20-25 bp long were designed in the most conserved regions for amplification of fragments between 280 and 900 bp.

PCR reactions were performed using these primers following the manufacturer's instructions.

After the addition of 1 µL of dNTP [10 mM], 3 μL of MgCl 2 [25 mM], 0.2 μL of Taq [5U/μL], 10 μL of activity buffer 5X, 0.5 μL of each primer (upstream and forward primers)

[100 μM], 34 μL of DEPC treated water and 1 μL of cDNA, 40 PCR cycles at 94°C for 1 min, 55°C for 1 min and 72°C for 1 min were carried out in an Eppendorf AG thermocycler.

PCR products were analyzed on 1% (w/v) agarose gel with ethidium bromide. The ethidium bromide stained bands were revealed by UV light on an illuminator.

After a purification step using the PCR purification kit (Qiagen) according to the manufacturer's instructions, cDNA was cloned with pGEM ® -T, (PROMEGA).

Successful insertion of the fragments was checked by PCR using T7 and SP6 universal primers following the manufacturer's instructions (40 cycles at 94°C for 1 min, 50°C for 1 min and 72°C for 1 min.). All the partial gene sequences obtained were submitted to GenBank under the accession numbers reported in Table1.

2.1.4.q-RT-PCR primer design and real time q-PCR

For each gene, specific q-RT-PCR primer pairs (table 1) were determined using the LightCycler probe design software (version 1.0, Roche).

Real time PCR reactions were performed in a LightCycler (Roche) following the manufacturer's instructions (one cycle at 95°C for 1 min and 50 amplification cycles at 95°C for 5 s, 60°C for 5 s and 72°C for 20 s).

Each 20 µL reaction contained 1 µL of activity buffer (Syber Green I (Roche), Taq Polymerase, dNTP), 3.2 µL of MgCl 2 [25µM], 2 µL of the gene-specific primer pair at a final concentration of 300 nM for each primer, 12.8 µL of DEPC-treated water and 1 µL of cDNA.

Specificity was determined for each reaction from the dissociation curve of the PCR product.

This dissociation curve was obtained by following the SYBR Green fluorescence level during gradual heating of the PCR products from 60 to 95°C.

Relative quantification of each gene expression level was normalized according to β-actin gene expression. Relative mRNA expression was generated using the 2-Δ CT method [START_REF] Livak | Analysis of relative gene expression data using realtime quantitative PCR and the 2(-ΔΔC(T )) method[END_REF].

2.1.5.Determination of sensitivity limit

In order to optimize the RNA extraction method, the sensitivity limit was quantified.

Ten dilutions of E. minima cultures were prepared to obtain theoretical numbers of cells ranging from 1.8×10 7 cells to 1.8 cells (1.8×10 7 ; 9×10 6 ; 1.8×10 6 ; 9×10 5 ; 1.8×10 5 ; 1.8×10 4 ; 1.8×10 3 ; 1.8×10 2 , 18 and 1.8 cells). Extraction and reverse transcription of RNA were performed as described previously on these different samples. The resulting cDNA was used for q-PCR with cox1 and act specific primer pairs as described above.

2.2.Experimental protocol to assess cadmium effects on Eolimna minima 2.2.1.Exposure conditions 250 mL of E. minima culture in suspension in modified Dauta medium were placed in 500-ml Erlenmeyers to reach a final concentration of 10 6 cells/mL. The organisms were directly exposed to three Cd concentrations (C 0 =0, C 1 =10±3.2 and C 2 =96±34.2 µgCd/L, mean values ± standard deviations over 14 days, from CdCl 2 stock solution, prepared from a 1000 mg vial, Merck, Germany). The cultures were maintained at between 16 and 18°C along the 14 days of exposure with a photon flux density of 59.19 ± 7.75 µmol.m -2 .s -1 and a 12:12 h light:dark cycle. The Erlenmeyer flasks were closed with sterile cotton wool and kept on an orbital shaker (60 rpm). Triplicate flasks were collected at each of the four different sampling periods (after 1, 2, 7 and 14 days of exposure to Cd) leading to 36 different experimental conditions. Thus 7.5 mL of culture were sampled after the exposure durations assigned for water Cd analyses, intracellular and total Cd concentration analyses and cell counting. Simultaneously, three 100 mL replicates were sampled for genetic analysis and stored at -80 °C.

2.2.2.Cd analysis

Cadmium concentrations in water and in diatom cells were determined by atomic absorption spectrophotometry (AAS) using a graphite tube atomiser after filtration through a GH polypropylene membrane (0.45 µm pores, 25 mm diameter, PALL) for diatom cell Cd concentration and on a 0.45µm Teflon filter (from CAS) for Cd concentration in water. The analytical method was simultaneously validated for each sample series by analysing standard biological reference materials (TORT 2, Lobster hepatopancreas; DOLT 3 dogfish liver; NCR/CNRC, Ottawa, Canada). 3 mL of nitric acid (65% HNO 3 ) were added to reference the material before digestion in a pressurized medium (borosilicate glass tube) at 100°C for 3 h. 15 mL of ultra pure water were added and samples were stored at 4°C before analysis.

Samples used for intracellular cadmium determination within diatoms were treated first with EDTA (10 mM), a strong metal complexing ligand, to remove the metal ions adsorbed to the diatom cell walls [START_REF] Behra | Copper and zinc content of periphyton from two rivers as a function of dissolved metal concentration[END_REF]. In this case, a drop of EDTA was added to each 2 mL sample and filtered after 10 min. Samples used for total cadmium determination in diatoms skipped this step. After 24h at 44°C, all samples were digested by nitric acid (3 mL of 65% HNO 3 ) in a borosilicate glass tube at 100°C for 3 h. The samples were then analysed by atomic absorption spectrophotometry as described previously.

2.2.3.Diatom cell density

1.5 mL aliquots were immediately fixed in formalin (37% formaldehyde, Prolabo, France) for counting. Each sample was counted in triplicate using a Nageotte counting chamber (Marienfeld, Germany). After 10 min ultrasonication, 200 µL of sample were placed on the counting chamber. The total number of individuals and the number of dead cells were recorded in 4 fields of the gridded area (1.25 μL each, 0.5 mm depth) under light microscopy at 400x magnification (Leitz photomicroscope). Distinction between dead and live organisms was estimated by the observation of the turgescence and colour of the chloroplast.

2.2.4.Statistical analysis

Statistical analyses were performed with R 2.12.0 using the Linear Mixed Effect model (library nlme). Homogeneity of variances was verified by Leven test and Fisher's LSD was calculated.

To assess the effects cadmium on E. minima by the molecular tools described above, expression of the nine target genes sequenced in this study was examined. The expression levels of the nine genes involved in antioxidant defences, or mitochondrial or photosynthetic metabolism were investigated by quantitative real-time PCR as described previously using the q-PCR primers designed in this study.

3.Results and discussion

3.1.Methodological developments

3.1.1.Extraction method and optimization

Electrophoresis revealed the presence of large quantities of q-PCR amplification products, underlining the efficiency of this method for breaking open the silicate frustules while preserving the RNA. Moreover, the analysis of triplicates from the cadmium exposure experiment by q-PCR showed the regularity in the threshold cycle of most of the genes (apart from 18S and 12S see below). As shown in figure 1, seven of the nine genes studied showed a percentage of variation within triplicates of under 8%. This homogeneity reveals the high efficiency of the extraction method and suggests a constant basal level in RNA production from these genes in E. minima.

In order to optimize the extraction protocol, the sensitivity limit was determined by performing q-PCR with act and cox1 q-PCR primers. Analysis of dissociation curves of the different amplification products revealed the amplification of cDNA for an initial number of cells in the range 1.8×10 7 to 1.8×10 5 . Below this number no amplification was detected through lack of sufficient genetic material. As the minimal number of cells needed to obtain q-PCR amplification was determined at 1.8×10 5 cells for the diatom E. minima, the minimum number of cells used in subsequent q-PCR analyses was fixed at 10 7 cells. The RNA extraction method developed in this study showed its high efficiency in the extraction of total RNA for the diatom E. minima. Moreover the method is simple and only requires standard laboratory equipment, which is an important advantage compared to existing methods that can require expensive and/or specialised apparatuses to extract RNA like a French pressure cell press [START_REF] Hildebrand | Characterization of a silicon transporter gene family in Cylindrotheca fusiformis: Sequences, expression analysis, and identification of homologs in other diatoms[END_REF] or MiniBeadBeater [START_REF] Fawley | A simple and rapid technique for the isolation of DNA from microalgae[END_REF]. The RNA extraction method developed is also rapid compared to existing methods using a muffle furnace overnight [START_REF] Wawrik | Real-time PCR quantification of rbcL (ribulose-1,5-bisphosphate carboxylase/oxygenase) mRNA in diatoms and pelagophytes[END_REF]. Moreover, the method has been successfully tested on some quite different diatom species (Achnanthidium minutissimum and Nitzschia palea) suggesting it will be applicable to a large range of diatom species.

3.1.2.Nucleotide sequence accession number and q-PCR primers

The 9 genes of interest selected in our study were successfully cloned and sequenced, accession numbers of sequenced genes and q-PCR primers are shown in table 1. cDNA sequences of the 9 genes have been deposited in the GenBank database.

3.2.Effects of cadmium on Eolimna minima

3.2.1.Cd water concentrations in solution

Cadmium concentration in water for the three different treatments is plotted against time in Figure 2. During the 14 days of exposure, mean Cd concentrations in water were 10±3.2 μg/L and 96±34.2 μg/L for nominal concentrations of 10 and 100 µgCd/L respectively. The cadmium concentrations in water were very stable from days 0 to 7 in both C 1 (11.5±0.6 μgCd/L) and C 2 (111.3±3.8 μgCd/L) treatments. A strong decrease occurred in cadmium concentrations on day 14 with values falling to 4.3±4 μgCd/L and 35.2±14.2 μgCd/L for treatments C 1 and C 2 respectively. This is unlikely to be explained by metal adsorption onto the surface of the experimental units, as no decrease was observed during the first seven days of exposure: bioaccumulation of cadmium by diatoms is more likely responsible for this decrease, as shown by the following results.

3.2.2.Effect of cadmium on the growth of E. minima

The number density, over the 14 days of direct exposure to the three different concentrations of cadmium, is plotted in Figure 3. Statistical analysis revealed that growth was not significantly different in controls and with a Cd concentration of 10.0±3.2 µg/L over the whole duration of the experiment. Besides, [START_REF] Gold | Impacts of Cd and Zn on the development of periphytic diatom communities in artificial streams located along a river pollution gradient[END_REF] and [START_REF] Duong | Seasonal effects of cadmium accumulation in periphytic diatom communities of freshwater biofilms[END_REF] observed a high proportion of small, adnate species like Achnanthidium minutissimum, Encyonema minutum and E. minima in stations contaminated with around 6 µg/L cadmium. This suggests the tolerance of these species to high levels of cadmium and thus their ability to develop in cadmium-contaminated conditions with around 10 µgCd/L. Nevertheless, densities in contaminated samples were on average lower than in control treatments throughout the 14 days of exposure. It could be interesting to study population kinetics over a longer period to see if, after 14 days of exposure, significant differences are observed between the C 1 treatment and the control as observed by Morin et al. (2008a) who reported differences in biofilm diatom densities between control and 10 µgCd/L-contaminated treatments only after 6 weeks of exposure.

After 7 days of exposure a significant difference occurred between the controls and the cadmium-contaminated treatments, but no such difference appeared between moderate and high contamination (C 1 and C 2 ) or between intracellular and total bioaccumulation. After 14 days of exposure, total bioaccumulation in the controls was significantly different from that in the cadmium-contaminated treatments (p = 0.008 and p < 0.001 respectively).

Total bioaccumulation was significantly higher than intracellular bioaccumulation for contamination C 2 (430.1 ± 86.4 μgCd/g dw versus 734.1 ± 70 μgCd/g dw).

The high values of cadmium bioaccumulated and the large increase between days 7 and 14 for the higher Cd concentration could so explain the total inhibition of growth observed for this treatment.

The ratios of adsorbed to absorbed cadmium are difficult to estimate owing to the heterogeneity of the results except for exposure to 100 µgCd/L after 14 days. In order to evaluate this difference better for shorter periods and lower concentrations, sampling volumes would have to be increased. In our study however, increasing sample volumes from 2 to 40 mL still led to very low dry weight measurements (between 2 and 6 mg for day 7 and 14) increasing the incertitude of measurement and errors in bioaccumulation calculations.

Therefore particular attention must be paid to the volumes tested when diatom cultures are considered.

3.2.4.Gene expression levels

In our study, two genes were selected as potential references, β-actin and 18S ribosomal RNA. Both are widely used as references in q-PCR analyses. The average percentage variation in the threshold cycle (Ct) of 18S ribosomal RNA within triplicates was 13.7% against 4.9% for β-actin (figure 1). The Ct value was defined as the number of cycles needed for the amplification signal to reach a specific detection threshold, and so is inversely correlated with the amount of cDNA template present in the PCR amplification reaction. The low stability and the high abundance of 18S rRNA underlined by our results (figure 5) shows the unsuitability of this gene as a reference. Consequently, 18S was finally not kept as a reference gene in our experiments. On the other hand, β-actin showed a high regularity and high threshold cycle values and was thus defined as reference gene in our study. 18S rRNA has also been described as a poor reference gene in other species such as buffalo (Bubalus bubalis) [START_REF] Castigliego | Selecting reference genes in the white blood cells of buffalos treated with recombinant growth hormone[END_REF] or humans [START_REF] Radonić | Guideline to reference gene selection for quantitative real-time PCR[END_REF]. A recent study underlined this point; [START_REF] Siaut | Molecular toolbox for studying diatom biology in Phaeodactylum tricornutum[END_REF] reported that 18S rRNA was a very poor reference gene for Phaeodactylum tricornutum, a marine diatom. Here, we chose β-actin as reference gene

Author-produced version of the article published in Water Research, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j. watres.2011.11.005 owing to the stability of its expression, although [START_REF] Siaut | Molecular toolbox for studying diatom biology in Phaeodactylum tricornutum[END_REF] show that the β-actin level appears to be tightly correlated to light with an increase during the light period. In our study, however, the photoperiod was controlled and sampling performed at a fixed time, which avoided this source of error and allowed us to use β-actin as reference gene. Nevertheless it would be interesting to study other reference genes as proposed by [START_REF] Siaut | Molecular toolbox for studying diatom biology in Phaeodactylum tricornutum[END_REF] (H4, RPS, Cdka and TBP) for future analyses and particularly in field samples.

Expression of genetic levels during the 14 days of exposure to cadmium at 10.0±3.2 and 96.0±34.2 µg/L are shown in table 2. During the 14 days, an impact resulting in the amplification of 5 of the genes studied was observed (cox1, nad5, d1, psaA and 12S).

Different genetic responses are expressed as a function of time and concentration of exposure.

Differential genetic expression revealed an effect on mitochondrial and photosynthetic metabolism, observed after 7 days of exposure only at the highest contamination pressure, and for both contamination pressures on day 14.

After 7 days of exposure, cox1, d1, psaA and 12S were up-regulated for the exposure concentration of 96.0±34.2 µgCd/L. After 14 days of exposure, there was a strong upregulation of nad5, d1 and psaA for both concentrations. Moreover, a large difference occurred between the 2 treatments: expression of nad5, d1 and psaA were 4, 12 and 19 times higher respectively for C 2 than C 1 . We can also see that responses of d1 and psaA in the C 2 treatment increased with exposure time by a factor of 4.4 and 6.4 respectively between days 7 and 14.

The disturbance in the expression of the cytochrome C oxidase (cox1) subunit and the NADH dehydrogenase subunit 5 (nad5) observed in our study suggests an impact on mitochondrial metabolism. This is in accordance with data in the literature. Indeed, after entering the cell, cadmium can inhibit electron transfer in mitochondrial respiration, and also lead to the formation of reactive oxygen species (ROS) [START_REF] Wang | Cadmium inhibits the electron transfer chain and induces Reactive Oxygen Species[END_REF][START_REF] Gonzalez | Comparative effects of direct cadmium contamination on gene expression in gills, liver, skeletal muscles and brain of the zebrafish (Danio rerio)[END_REF]). The oxidative stress generated then brings about DNA breaks and induction of apoptotic mechanisms leading to cell death (Bagchi et al. 2000, Chan and[START_REF] Chan | Cadmium-induced ectopic apoptosis in zebrafish embryos[END_REF]. Nevertheless, the ratio cox1/12S was constant throughout the experiment indicating that the number of mitochondria increases to maintain the level of ATP production. This increase of cox1 expression has already been reported in zebrafish Danio rerio [START_REF] Gonzalez | Comparative effects of direct cadmium contamination on gene expression in gills, liver, skeletal muscles and brain of the zebrafish (Danio rerio)[END_REF]. Electron microscopy studies also in the zebrafish have shown extensive disturbance of the ultrastructure of the mitochondria after metal exposure to methyl mercury [START_REF] Cambier | At environmental doses, dietary methylmercury inhibits mitochondrial energy metabolism in skeletal muscles of the zebra fish (Danio rerio)[END_REF]. In future work it will be of particular interest to use electron microscopy to view E. minima after Cd exposure so as to be able to evaluate damage at the mitochondrial level.

Secondly, with regard to the up-regulation of psaA and d1 after 7 days of exposure to 100 µgCd/L and after 14 days for the two contamination levels, it is clearly shown that the photosynthetic system of E. minima is strongly affected by exposure to cadmium. The impact of metals on the photosystem is well known and has been studied in numerous works.

Photosystem inhibition has been reported by [START_REF] Sudo | Gene expression and sensitivity in response to copper stress in rice leaves[END_REF] at physiological and genetic levels with, in particular, repression of d1 and psaA in rice leaves after copper stress. [START_REF] Hörcsik | Effect of chromium on photosystem 2 in the unicellular green alga, Chlorella pyrenoidosa[END_REF] studied the effect of chromium on PS2 activities in Chlorella pyrenoidosa cells: they showed the inhibition of PS2 by photodestruction of reaction centres in cells cultivated in the presence of chromium. At the ultrastructural level, [START_REF] Ouzounidou | Physiological and ultrastructural effects of cadmium on wheat (Triticum aestivum L.) leaves[END_REF] showed damage to chloroplast structures with changes in the shape and thickness of the thylakoid membranes in wheat (Triticum aestivum L.) after exposure to 265 μgCd/L for 7 days.

Moreover, the impact at the genetic level underlined by the up-regulation of nad5, d1 and psaA after 14 days of exposure revealed effects at a Cd concentration of 10 µg/L. This result is of particular interest because, at a more integrative level of population kinetics, this effect was not noted. Overall, the above findings stress the importance of considering molecular markers when studying pollution because of the potentially higher sensitivity and the earlier response of genomic markers to toxic pollution compared to more global endpoints like mortality or population kinetics.

No significant difference in genetic expression levels was observed for the genes cyp1A1 or sodMn during the experiment, at either of the concentrations tested. Mitochondrial superoxide dismutase is known to be involved in the oxidative stress response and several studies have shown the up-regulation of this gene after metal contamination. For example [START_REF] Gonzalez | Comparative effects of direct cadmium contamination on gene expression in gills, liver, skeletal muscles and brain of the zebrafish (Danio rerio)[END_REF][START_REF] Gonzalez | Comparative effects of dietary methylmercury on gene expression in liver, skeletal muscle, and brain of the zebrafish (Danio rerio)[END_REF] demonstrated the up-regulation of sodMn in the gills of zebrafish after 7 days of direct exposure to 86 ± 0.26 µgCd/L and in another study in skeletal muscle and liver after 21 and 63 days of dietary exposure to methyl mercury for concentrations of 5 to 13.5 µgHg/g dw.

In our work, this up-regulation was not observed, several hypotheses can be proposed: the sampling rate may be unsuitable and/or other enzymes of the oxidative stress response may be involved such as catalase, glutathione peroxidise or other types of superoxide dismutases (sodFe, sodCu/Zn or sodNi). Other authors have studied antioxidant enzyme activities, for example [START_REF] Branco | Sensitivity of biochemical markers to evaluate cadmium stress in the freshwater diatom Nitzschia palea (Kützing) W. Smith[END_REF] observed an increase in catalase activity compared to controls after 5 days of exposure to Cd from 100 to 300 µgCd/L on the diatom Nitzschia palea and in the global superoxide dismutase activity for concentrations from 200 to 300 µgCd/L. In a 13

Author-produced version of the article published in Water Research, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005 field study, [START_REF] Bonet | Antioxidant enzyme activities as biomarkers of metal contamination in biofilms: field validation using[END_REF] showed a clear inhibition of glutathione-S-transferase at high metal exposure and no difference in superoxide dismutase activity between sites impacted by metals at different levels.

For future analyses it will be interesting to sequence genes coding for catalase, gluthatione peroxidise, glutathione-S-transferase and the other superoxide dismutases (sodFe, sodCu/Zn or sodNi) in order to study their responses to metal exposure together with enzymatic activities. Cytochrome P450 1A1 plays a role in the metabolism of organic compounds such as PAH or pesticides, which explains the null genetic response observed in our study after cadmium contamination. This gene has been studied and q-PCR primers have been developed for use in future analyses for diatom exposure to compounds like PAH or PCB.

4.Conclusions and perspectives

In this study, a new glass-bead RNA extraction technique for diatoms was successfully developed and optimized. Nine genes of interest were sequenced for Eolimna minima allowing the application of q-PCR tools to this species.

Our results underlined the toxicity of Cd towards E. minima population kinetics only for the highest concentration, while q-PCR analyses revealed an impact on mitochondrial metabolism and the chloroplast photosystem for both Cd exposure concentrations. Genetic expression of nad5, cox1, 12S, d1 and psaA by q-PCR could thus constitute an early warning biomarker of metal pollution. Future studies should investigate sequences of genes coding for catalase or glutathion peroxidase in order to study the response to oxidative stress.

The present study is the first reported use of q-PCR on river benthic diatoms and the results obtained are extremely promising. The techniques developed were successfully tested using simplified mixtures of diatom species. Further interesting steps for the early and sensitive assessment of metal pollution would firstly involve validating the results obtained by examining sensitive vs tolerant diatom species response levels, then finding or confirming genetic biomarkers for use on natural multispecific biofilms for impact assessment of toxic pollution.

Author-produced version of the article published in Water Research, 2012, 46 (4)

Figure 2 .Figure 3 .

 23 Figure2. Cadmium concentration in water ± standard deviation versus exposure time for the 3 treatments (C 0 =0, C 1 =10 and C 2 =100 μgCd/L). Full and empty triangles represent contaminated conditions. Full squares represent control treatments.

Figure 4 .

 4 Figure 4. Intracellular and total Cd bioaccumulation in E. minima cultures for the two exposure times and three treatments (C 0 =0, C 1 =10 and C 2 = 100μgCd/L) in µg Cd/g dw. Significant differences between a, b, c, d and e (p value <0.05). White bars represent control treatments, grey bars represent C 1 treatments and black bars represent C 2 treatments. Full and hatched bars represent total and intracellular cadmium bioaccumulation respectively.

Table 1 .

 1 , 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005 Accession numbers and specific primer pairs for the 9 E. minima genes used in our study. -produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005

	Gene name	Accession number	Primer (5'-3')
	sodMn	HM 449706	GGTAGTAGGCGTGCTCCC a CCAGGACAACCCGCTC b
	cox1	HM 449704	CAGTAATTCTCACTGCCCAGC a CCGTGTACCCACCGTTG b
	nad5	HM 449708	TCAACTTGGTTTGCATACATGGC a TTGAACTAATCCTGTTGTGGAAGC b
	d1	HM 449711	ACCACCAAATACACCAGCAAC a GCGTCCTTGGATTTTCGTAGC b
	psaA	HM 449705	CATAAAGCGGCACCCAAAC a CTTGGATATAACTGACTCATTAACTCAGG b
	cyp1A1	HM 449709	ACGAAGGAGTCTCCCATCA a CTTGGCCCGAACGATCAG b
	act	HM 449707	GGCTCCACAAAACCCCAAG a GGCGTACCCCTCGTAGAT b
	12s	HM 449710	CGCGGTAATACGGAGGATGC a AGTGCCTTCGCCATCGG b
	18s	HM 449712	CATTGTCAGAGGTGAAATTCTTGGA a CCCCGGAACCCAAAAGT b

Abbreviations: sodMn-mitochondrial superoxide dismutase; cox1-cytochrome C oxidase subunit I; nad5-NADH dehydrogenase subunit 5; d1-D1 protein; psaA-PsaA protein; cyp1A1-cytochrome P450 1A1; act-βactin, 12s-mitochondrial ribosomal RNA 12S; 18s-ribosomal RNA 18S. a Upstream primer. b Forward primer.

Author

Table 2 .

 2 Differential gene expression as compared to actin from E. minima after 1, 2, 7 and 14 days of cadmium exposure to 10 and 100 μgCd/L by direct route. Average variation of threshold cycle within triplicates in percentage over the cadmium exposure experiment for the nine genes studied.Author-produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005

	Functions	Genes Cadmiu contaminate x erimenta units m d e p l
		C1 (10.0 3.2 μg ) ± /l		C2 (96.0	.2 μ l) ± 34 g/
		1	2	7	14	1	2	7	14
	Mitochondrial	cox1

  Author-produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005

  Author-produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005 Figure 5. Expression of the two potential reference genes 18S and act in all treatments. mRNA levels of each reference gene is shown in absolute Ct values during the 14 days of the experiment. Closed and open circle represent 18S and act respectively.Author-produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005

								18S
								act
		30
		25
	Ct value	15 20
		10
		5
		0
		0	2	4	6	8	1 0	1 2	1 4
			a a		a a Exposure time (days)

Author-produced version of the article published in WaterResearch, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j.watres.2011.11.005

Author-produced version of the article published in Water Research, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j. watres.2011.11.005 Significant differences in cell densities were observable from day 7 between the controls and the series with the highest contamination of cadmium (2 900 000 ± 250 000 cells/mL for control versus 1 300 000 ± 500 000 cells/mL). After 14 days of exposure, the diatom cell density was 3.2 times higher than with the control treatment (4 200 000 ± 390 000 cells/mL versus 1 300 000 ± 280 000 cells/mL). Moreover, diatom growth was null with treatment C 2 . Our results are in accordance with those of [START_REF] Tourtelot | Contamination expérimentale par le cadmium et le zinc de diatomées périphytiques d'eau douce, Achnanthidium minutissimum, en conditions contrôlées de laboratoire[END_REF] who studied the acute effects of Cd on Achnanthidium minutissimum. For this species, referenced as pollution resistant, the rate of growth was null after 14 days of exposure to 100 µgCd/L. Similar results were observed at the community level, Morin et al. (2008b) exposed bioflms in microcosms to 10 and 100 µgCd/L for 6 weeks and observed a significant reduction in diatom density for biofilms exposed to the higher concentration compared to the control and the low cadmium contamination treatment.

3.2.3.Cd bioaccumulation in diatoms

Intracellular and total bioaccumulation are not given for days 1 and 2 after exposure because of the insufficient quantity of cells in the samples. Intracellular and total Cd concentration in diatoms after 7 and 14 days of exposure to the different treatments (0, 10±3.2 and 96±34.2 µgCd/L) are reported in Figure 4.

Intracellular and total bioaccumulation increased significantly between days 7 and 14 for the 2 contaminated treatments. For the lower contamination concentration (C 1 ), the intracellular and total concentration increased from 3.6±1.8 to 51.4±5.3 μgCd/g dw and from 6.4±3.4 to 57±6.3 μgCd/g dw respectively. For the higher contamination pressure (C 2 ), there was an increase from 5.6 ± 2.7 μgCd/g dw to 430.1 ± 86.4 μgCd/g dw and from 11.4 ± 3.0 to 734.1 ± 70 μgCd/g dw respectively, this latter value is in the range of literature values (Morin et al. 2008a, Duong et al. 2008). The increase in cadmium bioaccumulation with exposure time was also shown in environmental samples by [START_REF] Duong | Seasonal effects of cadmium accumulation in periphytic diatom communities of freshwater biofilms[END_REF]. The cadmium content in the biofilm increased gradually for days 14 and then remained constant. This could be due to the saturation of binding sites in the biofilm leading to a limitation of Cd accumulation. In our experiment it is difficult to estimate if the saturation phase is reached after 14 days of exposure. In order to determine the kinetics of saturation and the saturation value for E.

minima, experiments will have to cover a longer period with a constant contamination pressure.

Author-produced version of the article published in Water Research, 2012, 46 (4), 934-942. The original publication is available at http://www.sciencedirect.com/science/article/pii/S0043135411006774 doi : 10.1016/j. watres.2011.11.005