Françoise Ehrburger-Dolle
email: francoise.ehrburger-dolle@ujf-grenoble.fr

Isabelle Morfin

Françoise Bley

Frédéric Livet

Gert Heinrich

Luc Piché

Mark Sutton

Tensile stress relaxation and recovery behavior of a cross-linked EPDM rubber matrix loaded with different fillers

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1

Tensile stress relaxation and recovery behavior of a cross-linked EPDM rubber matrix loaded with different fillers

Françoise Ehrburger-Dolle* 1 , Isabelle Morfin 1 , Françoise Bley 2 , Frédéric Livet 2 , Gert Heinrich 3 , Luc Piché 4 , Mark Sutton 4 1 Univ. Grenoble 1 / CNRS, LIPhy UMR 5588, Grenoble, F-38041, France 2 SIMaP, UMR 5266 Grenoble INP / CNRS / UJF, 38402 Saint Martin d'Hères, France 3 Leibniz-Institut für Polymerforschung Dresden, 010169 Dresden, Germany 4 Physics Department, McGill University, Montreal, Quebec H3A 2T8, Canada ABSTRACT. Tensile stress-relaxation measurements have been performed on a series of crosslinked filled elastomers. The samples are constituted by the same Ethylene-Propylene-Diene-Monomer (EPDM) matrix. The different fillers are chosen in order to investigate the effect of the filler-filler and the filler-matrix interactions on the time dependence of the relaxation modulus ( ) at different strain values during about 4000 s. Each strain step is reached by a change of strain that is positive during the first (UP1) and the second (UP2) elongation or negative (DOWN) when the strain plateaus are decreased. Analysis of the experimental data indicates that most of the curves obtained after a positive strain change (UP1 or UP2) yielding a decrease of ( ), can be fitted by the empirical Chasset-Thirion equation ( ) [ ( ⁄ ) ]. For carbon black samples characterized by a strong filler-matrix interaction, the value of the exponent m may become very small and the data are better fitted by a logarithmic relation

( ) ( )[ ( )].
As for the increase of ( ) measured after a negative strain step (DOWN), it is shown, for the first time, that all data can be fitted by the following equation:

INTRODUCTION

The mechanical behavior of filled elastomers is an active research area since more than 50 years mainly because it is related to an ancient and important application which is the reinforcement of rubber. Several recent reviews in the literature provide information about new developments in the field of polymer nanocomposites. Basic principles and technological applications have been reviewed by Vilgis et al. [START_REF] Vilgis | Reinforcement of Polymer Nanocomposites: Theory, Experiments and Applications[END_REF] and Heinrich et al. [START_REF] Heinrich | Constitutive Models for Rubber VI[END_REF]. The latest advances in the viscoelastic behavior of rubbery materials were described and analyzed by C. M. Roland [START_REF] Roland | Viscoelastic Behavior of Rubbery Materials[END_REF]. A recent book [START_REF] Thomas | Rubber Nanocomposites: Preparation, Properties and Applications[END_REF] is dealing about preparation, properties and applications of rubber nanocomposites. A review dealing with structure, phase behavior and properties of nanocomposites has been published by Kumar and Krishnamoorti [START_REF] Kumar | [END_REF]. Oberdisse et al. [START_REF] Oberdisse | Structure Determination of Polymer Nanocomposites by Small Angle Scattering[END_REF] reported recent investigations devoted to the determination of the structure of polymer nanocomposites by means of small-angle scattering techniques. A survey of the present understanding of polymer-based nanocomposites as obtained from either theoretical or computational approaches was presented by Allegra et al. [START_REF] Allegra | [END_REF].

The mechanical behavior of filled elastomers is characterized by two major features. The first one concerns the non-linear viscoelastic behavior (or Payne effect) revealed by the amplitude dependence of the storage and the loss moduli in shear or tensile measurements [8][9][10]. The second one is stress-softening (or Mullins effect) yielding the lowering of the resulting stress for the same applied strain after the first load and, as a consequence, mechanical hysteresis. The Mullins effect was recently reviewed by Diani et al. [11] and Amin et al. [12]. Bhattacharya et al. [13] reviewed several "classical" models attempting to take into account the complexity of the problem. These authors also reported an extensive study of the nonlinear, time-dependent mechanical behavior of a cross-linked carbon black filled elastomer. The unusual non-linear and plastic behavior of nanofilled elastomers was recently attributed to kinetics of rupture and rebirth of glassy bridges [14][15][16]. Resting time and relaxation features are also involved in partial [17] or complete [18,19] recovery (or healing) of the Mullins effect at room temperature.

Early measurements of tensile stress relaxation in a series of carbon black N330 filled rubber vulcanizates were reported by MacKenzie and Scalan [20]. The experimental data could be fitted almost equivalently either by a logarithmic function of time or a power law with exponents smaller than 0.07. The authors produced evidences for the existence of a distinct relaxation process below 60 seconds that vanished above a given strain, depending on the amount of filler.

Similar results were obtained for vulcanized ethylene-propylene-diene-monomer (EPDM) rubber filled with carbon black N660 [21]. Geethamma et al. [22] showed that the slope of the lines obtained in semi-logarithmic coordinates was larger for plain natural rubber (NR) than for filled ones. The increase of the amount of filler consisting of short coir fibers lead to the decrease of the absolute value of the slope. Opposite features were reported for NR filled with titania or silica particles [23]. In the case of nanogel filled elastomers, Mitra et al. [24] fitted the experimental values of the relaxation modulus with power laws in which the exponents varied between 0.2 and 0.3. More precisely, these authors used the following equation:

( ) [ ( ⁄ ) ] (1)
which has been proposed in 1965 by Chasset and Thirion [START_REF] Chasset | Proceedings of the Conference on Physics of Non-Crystalline Solids[END_REF] for describing the isothermal tensile relaxation modulus data ( ) for many rubbers. In this equation, is the equilibrium modulus and  0 and m are parameters which depend on the system investigated. The same equation applies also to the shear modulus relaxation ( )

( ) [ ( ⁄ ) ] (2)
An analogous equation was used by Dickie and Ferry [START_REF] Dickie | [END_REF] for fitting the creep compliance relaxation ( ) ( ) in dicumyl peroxide vulcanizated natural rubbers. The Chasset-Thirion (CT) equation has been widely used to describe power law stress relaxation in unfilled polymers or elastomers [27,[28][29][30][31][32]. In several cases, the exponent m was shown to increase linearly with the cross-linking density. Theoretical models have been proposed by Curro et al [33,34] and Heinrich and Vilgis [35] for explaining this feature and for justifying a posteriori the CT equation. However, several authors have shown that the linearity of the relation between m and the cross-linking density fails, as in the case of end-linked polydimethylsiloxane (PDMS) elastomers with pendent chains [27]. More, a decrease of m when the cross-linking density increases was reported by Martin et al. [30] for EPDM networks cross-linked with resol under compressive strain at 100°C.

In fact, there are many classes of materials in which stress relaxation is not exponential but best represented as log(t) [36] or as a power law in time [37]. Examples of such materials include, among others, block copolymers [38] and, more generally foams, emulsions, soft colloidal glasses, food products and others [37]. Power-law stress relaxation implies also power-law behavior in the frequency-dependence of the viscoelastic storage G'() and loss G''() moduli in materials for which there is no single characteristic relaxation time. Similarities in the rheology of many soft materials led Sollich and coworkers [39][40][41] to propose the soft glassy rheology (SGR) model. Moreover, SGR features are also observed in many biological systems [42][43][44][45][46][47][48][49][50]. It appears, however, that the elementary biophysical and biochemical mechanisms which govern the cell mechanics may not be well described by the SGR model. Therefore, other models involving a broad distribution of relaxation times were proposed [48] to explain power law relaxation. Jaishankar and McKinley [37] presented recently a constitutive framework using fractional derivatives to model the power-law responses observed experimentally in very different complex systems. These authors also show that the "effective noise temperature" x appearing in the exponent in the SRG model is intimately related to the exponents involved in their fractional model.

This non-exhaustive bibliographic survey concerning stress relaxation or recovery processes suggests that the power law involved in the CT equation could not be limited to the particular situation of cross-linked elastomers. It follows also that the physical meaning of the exponent m could be different from that involved in the Curro model since the power law behavior is typical of structural disorder and metastability [37,[39][40][41] and implies a broad distribution of relaxation times. We have shown in a recent paper [51] that the decrease of the force ( ) requested to maintain an elongation of 60% of cross-linked silica filled elastomer can be fitted by a power law equation formally similar to the CT equation and an exponent equal to 0.17. In the case of a cross-linked carbon black filled elastomer, under the same conditions, the force decreases as log(t). Heterodyne X-ray photon correlation spectroscopy (HD-XPCS) measurements performed simultaneously revealed a significant difference in the value of the velocity of the filler particles.

For the carbon black filled sample, the logarithmic decrease of the tensile force was associated with a small velocity of the filler particles that decreased rapidly with time. These features would agree with the concept of glassy bridges in filled elastomers [14][15][16]. On the opposite, for the hydroxylated silica sample, the aggregate velocity was about 10 times larger than that of the carbon black ones and its decrease with time was slower. These measurements suggested (i), similarities between the behavior of the silica filled sample and that of colloidal gels or soft glasses [52] [53] and (ii) the absence of glassy domains as also observed in other silica filled elastomers [54][55][56]. This paper also emphasized the effect of the filler-filler (strong for hydroxylated silica fillers) and the filler-matrix (strong for carbon black fillers after cross-linking of the matrix [57]) interactions on the stress relaxation behavior of nanocomposites at the macroscopic and at the mesoscopic scale.

Progresses in fundamental knowledge on polymer-particle interaction and its link with macroscopic properties were recently reviewed by Pandey et al. [58]. Lin et al. [59] investigated the role of the nanofillers on the polymer dynamics in order to get new insights into understanding processing conditions and rheological behavior of polymer nanocomposites. The work described in the present paper was performed with similar objectives. Its aim is to show that the behavior of the relaxation and recovery of filled elastomers may provide a link between macroscale mechanical features and nanoscale polymer-filler interaction. This link could be the power law exponent m (or more precisely, the value of ). To this end a well-defined series of samples have been prepared. They consist of the same Ethylene Propylene Diene Monomer (EPDM) cross-linked rubber filled with an hydroxylated silica (strong filler-filler interaction), an hydrophobic silica having the same morphology (intermediate filler-filler and filler matrix interaction) and carbon black (strong filler matrix interaction) at two different concentrations. In order to determine the tensile relaxation and recovery (under strain), measurements were performed at different steps of strain reached after a positive step increase  where and being the initial strain or a negative step where . Tensile relaxation was also measured during the second increase of strain in order to get new information about the Mullins effect. Relevancy of the SRG model to strained filled elastomers will be discussed in the last paragraph. The thickness e of the samples used for the stress-strain measurements was about 2 mm. For the tensile stress relaxation, e was equal to1 mm. For all measurements, the one or two millimeter thick plates were punched out to the classical dumb-bell shape. The width l of the gauge ranged between 4 and 4.2 mm depending on the puncher used. Its length was equal to 25 mm.

MATERIALS AND EXPERIMENTAL METHODS

Tensile force measurements.

The tensile stress-elongation curves (monotonous stress-strain curves) were measured according ISO 527-2 with a universal testing machine Zwicki Z2.5. at the IPF in Dresden using S2 specimen geometry. The stretch ratio is defined as ⁄ where L and L 0 are the final and the initial lengths respectively. The strain  is defined as

( ) ⁄
(engineering strain). The strain rate ⁄ was equal to 0.133 s 1 corresponding to a crosshead speed of 200 mm/min and initial clamping length of 25 mm.

Tensile relaxation measurements performed at given steps of strain (incremental stress-strain curves) were achieved by means of the apparatus designed for simultaneous stress and XPCS study [51]. A complete description can be found in the master thesis of L.J.S. Halloran [START_REF] Halloran | A Rheological Study of Stress Relaxation in Elastomers for in situ X-Ray Diffraction Measurements[END_REF].

Each strain step was reached after a rapid extension (0.034 s 1 ) from the previous step to the next one. A strain gauge and a digitizing unit were used to measure the evolution of the force F every 1 or 2 s at each elongation step. The value of the stress  is deduced from F by ⁄ where is the section of the sample. The tensile relaxation modulus ( ) is defined as ( ) ( ) ⁄ . Fitting of the experimental curves is achieved by non-linear regression procedures by means of the Marquardt -Levenberg algorithm provided by SigmaPlot 10.0.

ANALYSIS OF THE STRESS-STRAIN CURVES

3.1. Monotonous stress-strain curves. Figure 1a shows the variation of the stress as a function of the strain  plotted in logarithmic coordinates for all samples. In the small elongation domain, all curves exhibit a power law behavior ( ) . The value of the exponent n is nearly the same for all curves ( ). The domain in which this power law is verified, however, is not exactly the same for all samples. For the plain EPDM elastomer the power law is observed for values ranging between 0.03 and 0.26.( ) For the CB filled samples, the limits are 0.015 and 0.14. In the case of the silica samples, the limits of the power law domain are shifted to smaller (0.004 to 0.06) values for Si-OH-40-c or the domain becomes narrower (0.02 to 0.06) for Si-OR-40-c. In the region of small strains, a linear variation of with  corresponding to an elastic deformation is expected. Thus, the slope n of the logarithmic plots shown in Figure 1a should be equal to 1. Because tensile stress-strain data obtained for filled elastomers are never plotted in logarithmic coordinates, it is difficult to determine whether or not such non-linear effect was already reported. Power law relationships can be found in the case of uniaxial tensile tests on metals and alloys (Hollomon equation) [START_REF] Wilson | [END_REF]68]. In the case of polymeric materials [69] or fiber composites under elongation [70] or under compression [71], high strain rates induce viscoplasticity. From this non-exhaustive list of publications it can be concluded that the stress-strain power law behavior observed at small deformations could reveal plastic or elastoplastic effects, i.e., structural damages or failures. The existence or not of a relation between a power law behavior and the hysteresis associated to strain-softening characteristic of the Mullins effect would require further investigations. Thus, the origin of the power-law behavior shown in Figure 1a can hardly be found without further experimental work. The fact that the value of the exponent is nearly the same for all samples suggests that the hot-pressing step common to all samples, could induce a particular internal structure that would be almost irreversibly destroyed during the first stretch. It is worth mentioning the occurrence of plastic events in model polymeric glasses at very small strains, in a regime where the material is generally described as elastic [72].

The existence of a power law domain at small elongations allows us to calculate the relative stress  r by dividing the experimental values by the prefactor  0 . The collapsed curves plotted in Figure 1b reveal differences in the mechanical behavior of the samples at larger values. The first comment concerns the two CB samples: up to 0.7 the increase of the stress is nearly similar to that of the un-filled cross-linked elastomer. Thus it may be assumed that the carbon black particles act as additional cross-links. At larger strains, the shape of the curves indicates that stiffening occurs without the foregoing of a well-defined yield region.

On the opposite, for the OH-silica filled sample Si-OH-40-c that appears to be the stiffest at small strains, there is a well-defined yield region corresponding to strain softening. It is likely that this effect results from a progressive opening of interparticle H-bonds in the silica network.

Figure 1b shows that the relative strain softening effect decreases when the amount of OH-silica decreases from 40 to 20 phr.

Incremental stress-strain curves.

The curves obtained for the different samples are shown in Figure 2. In all cases, the difference between each strain step is equal to 0.2. The maximum strain is 0.6 for Si-OR-20-c, Si-OH-40-c and CB-40-c. For the 3 other samples, the limit is 0.8.

Sample Si-OR-20-c broke during the second elongation (UP2) at 0.6. indicates that the value of the maximum of the stress at each step is systematically smaller than that measured at the same strain shown in Figure 1a. This result may be explained by the fact that the change from one strain step to the next one is slower (0.034 s 1 ) than the continuous increase of strain (0.133 s 1 ) [73]. It appears (Figure 3a) that the two series of data are simply related by a multiplying factor for strain values between 0.2 and 0.6. It follows that the mechanism of the increase of stress with strain is the same in both cases in the 0.2-0.6 strain range. Above 0.6, the stress seems to be slightly overestimated. Furthermore, the stress values measured for CB-40-c by strain steps equal to 0.2 and 0.1 (open and closed circles in Figure 3, respectively) are nearly the same.

Figure 3b shows the evolution with strain of the tensile modulus ( ) ( ) ⁄ . In the investigation of the tensile relaxation modulus ( ) at a given , presented in the next paragraphs, the data point shown in this figure correspond to the initial values named ( ) (after suppression of the shifting factor introduced in the plots). This figure indicates that the amplitude of the decrease of the tensile modulus is significantly larger for the hydroxylated silica filled sample (Si-OH-40-c) than for the carbon black filled one (CB-40-c). For the other samples, the amplitude of the decrease stays between these two limits, as it is also shown in Figure 1b. phr. Furthermore, for CB-40-c (UP1), the concavity is directed towards the time axis while for all other UP curves, it is the convexity. For the latter, the Chasset-Thirion (CT) equation (Equation 1) properly fits the experimental data for as shown in Figure 4. The values of the parameters deduced from the fits collected in Table 1 will be discussed in paragraph 5. We have observed that the experimental values E(0) are systematically very close to E(1) calculated by using the fitting equation (Equation 1).

As for the data obtained for CB-40-c (UP1), the empirical fitting equation becomes:

( ) ( )[ ( ⁄ ) ] (3)
with m < 0 as shown in Table 1. It follows that ( ) decreases from the initial value ( ) without reaching an asymptotic equilibrium value ( ) as it is the cases for the other UP curves shown in Figure 4.This feature could suggest the occurrence of plastic effects (damages). When the strain step is negative i. e., when the strain is decreased from  0.60 to 0.40 (DOWN curves), ( ) increases with time. It appears that all data (including CB-40-c) can be properly fitted by using the following equation:

( ) [ ( ⁄ ) ] with m > 0 (4)
To the best of our knowledge, such a modified CT equation has never been proposed for fitting the after-strain recovery. The recovery process was often described by kinetic aggregate breaking and reforming (flocculation) [74][75][76] described by exponential equations. In fact, figure 4 shows that the approach of the tensile relaxation modulus E(t) to a steady state differs depending upon the direction of the change in the strain amplitude as observed earlier by Wang and Robertson

[77]. These authors considered this asymmetry as one of the evidences of the existence of an analogy between strain-induced nonlinearity in the modulus of filled rubbers, the physics of the glass transition of glass-forming materials, and the jamming transition of vibrated granular materials. As a consequence, the asymmetric behavior cannot be interpreted by simple kinetic models. It seems now that the modified CT equation could capture the physical meaning of the unusual phenomena observed in filled rubbers.

These features which will be discussed in more details in paragraph 5, allow us to introduce an alternative equation that can be considered as a generalized Chasset-Thirion equation writing:

( ) [ ( ⁄ ) ] (5)
in which: for positive strain steps ("UP", )

for negative strain steps ("DOWN", ) involving recovery.

Table 1 indicates that the exponent m is always positive except for CB-40-c (Equation 3). In order to take into account this particular case in Equation ( 5), one will add the following condition: when m is negative. is the asymptotic equilibrium tensile modulus.

Interestingly, for all samples (except CB-40-c), the asymptotic values (DOWN) are very close from (UP2) determined during the second step elongation.

The characteristic times   (Table 1) resulting from the fit vary considerably from sample to sample. In the present paper,   will be considered as a characteristic time scale but it will not be further discussed. Concerning samples CB-20-c and CB-40-c (UP1), the best fit is obtained by a logarithmic equation:

Experimental results

( ) ( )[ ( )] (6)
The slope K is very similar and small (0.016 and 0.018, respectively) for both samples.

Samples Si-OR-20-c and CB-40-c (UP2) behave as does CB-40-c (UP1) for . m is negative and the data follow the generalized CT equation (Equation 5) with . However, the absolute values of the exponents m are smaller than it was for the smaller strain. Most of them display a minimum located between 100 and 150 s followed by a domain in which ( ) reaches a nearly constant value considered as being the asymptotic limit E   For sample Si-OH-40-c (UP1), the short time (t < 40 s) behavior can be described by a logarithmic decrease of ( ). The slope K is equal to 0.118 which is significantly larger than it is for the samples at larger strains characterized by a logarithmic decrease of ( ) over the whole time domain investigated ( ) is equal to 3.48 MPa. A logarithmic decrease of ( ) is also observed for the sample Si-OR-20-c (UP1) in a larger range of time (up to 440 s). The slope is equal to 0.036 and ( ) equals 1.27 MPa. For all other samples, the time dependence of the modulus at short times, up to the minimum, is nearly sigmoidal (in the semi-logarithmic representation). Existence of an initial fast decay of the relaxation modulus and a final slow decay at small strains was also reported in the case of carbon black filled natural rubbers or butadiene rubbers [20] or in carbon black filled EPDM elastomers [21]. The change of regime was often observed around 100 s as in the present case. At larger strains, the initial fast decay was no longer observed as also shown in Figures 4 to 6. Because this initial fast decay was not observed in un-filled rubbers, this feature was attributed by MacKenzie and Scalan [20] to a rapid re-formation of the filler network not completely destroyed at small strains.

Figure 7 show that the DOWN curves measured for  0.20 are all well fitted by Equation ( 5) (in which ) as for all DOWN curves obtained at larger strains.

Experimental results obtained for the uncross-linked

Si-OH-40-u sample. Figure 8 shows the curves ( ) obtained for this uncross-linked sample. It appears that their shape is very similar to that obtained for the corresponding cross-linked samples. For  = 0.10, the tensile relaxation modulus decreases logarithmically with time over the same time interval (1 to 40 s) as for Si-OH-40-c sample (Figure 7a). ( ) is equal to 3.33 MPa. The slope K = 0.189 is larger for this sample than for the cross-linked one but the strain is not the same ( = 0.10 for the former and 0.20 for the latter). For  = 0.20 and above, the curves are well fitted by the CT equation.

The values of the exponent m varying between 0.27 and 0.21, are slightly larger in the uncrosslinked sample than in the cross-linked one. In the same way, the characteristic times   comprised between 0.5 and 3 seconds are significantly larger in the uncross-linked sample than in the cross-linked one. A similar result was also reported for natural and styrene-butadiene rubbers [24].

DISCUSSION

Physical meaning of the exponent m in the Chasset-Thirion equation

We have shown in the previous paragraph that all ( ) data measured at a given strain  larger than 0.20 and reached by increasing the elongation (UP curves) can be fitted by one of the following equations:

-the generalized Chasset-Thirion equation (Equation 5): ( )

[ ( ⁄ ) ]
-a logarithmic equation (Equation 6): E( ) ( )[ ( )]

In both cases, the derivative which corresponds to the rate of change of the tensile modulus writes as:

| ( ) ⁄ | (7)
with for the CT equation and for the logarithmic decrease (Equation 6). When  is smaller than 1, m is negative. - is larger or close to 1.1 for the cross-linked hydroxylated silica samples (Si-OH-c) and larger than 1.2 for the uncross-linked one (Si-OH-u) - is smaller or close to 1.1 for the hydrophobic silica samples (Si-OR-c) and tends to 1 for the CB samples characterized by a log(t) decrease of tensile modulus - depends on the initial strain and on the filler-filler and filler-matrix interactions - is always strictly larger than 1 during recovery (DOWN).

As a consequence, the exponent and, therefore, the CT exponent may be considered as a fingerprint of the time-dependent properties of strained viscoelastic materials.

Within the approach of soft glass rheology (SGR) [39][40][41], the power law exponent generally labeled ( ) characterizing the time-dependent viscoelastic properties is equal to 0 ( ) as the glass state is approached and to 1 ( ) for a viscous fluid. Our results suggest that the Chasset-Thirion exponent m is equivalent to the SRG exponent ( ). In such conditions, in equation ( 7) would correspond to x considered as an "effective noise temperature" [40] (also named "mechanical noise temperature"). In the present work, the value of the exponent m is not a characteristic of the material only but depends also on the strain (intensity and rate) and on its direction (UP or DOWN). As an example, for sample CB-20-c m is small ( = 1.04) for (UP1). For (UP1), the tensile relaxation modulus decreases as log(t) which corresponds to  = 1 in Equation [START_REF] Allegra | [END_REF]. It is likely that this result can be interpreted as the strain induced jamming described by Wang and Robertson [77]. When the strain is reached in the decreasing step starting from , ( ) follows the generalized CT equation (

) with an exponent (Figure 5). This feature agrees with the asymmetrical behavior discussed by these authors [77] already mentioned in §4.1. It suggests jamming resulting from the increase of the strain between 0.40 and 0.60 and unjamming when the strain is decreased from 0.80 to 0.60. More generally, we think that jamming could be revealed by the logarithmic time dependence corresponding to in Equation 7. In the SGR approach indicates that the system becomes close to the glassy state which would be reached at the at the glass transition temperature T g of the plain unstrained elastomer. These features would agree with the concept of glassy bridges induced by strain in filled elastomers [14][15][16]. It agrees also with the increase of T g observed in filled elastomers [59]. It should be mentioned, however, that the physical meaning of the words jamming and glass transition is not exactly the same [78]. Both words are used here within the scope of reference [77].

On the opposite, none of the silica filled samples exhibits this behavior, at least for strains . The UP1 and the DOWN curves can all be fitted with the generalized CT equation.

Accordingly, as shown in figure 9, the exponents  never reach the limit . In the SRG model, it means that the system does not reach a glassy state under strain because x remains strictly above 1. This remark is in agreement with the absence of glassy domains observed in other silica filled elastomers [54][55][56]. More generally, the above comments agree with conclusions reached by Luo et al. [79] who combined NMR and mechanical measurements on silica-silane and carbon black filled styrene-butadiene rubber (SBR). Figure 9 also shows that the values of are decreasing for the following filler sequence: hydroxylated silica, hydrophobic silica and carbon black. Interestingly, is the largest for the uncross-linked sample filled with hydroxylated silica. In the SRG model, it means that the mechanical noise temperature x and, therefore, the chain mobility, decreases accordingly. This observation is consistent with the fact that cross-linking and/or filler matrix interaction decrease the chain mobility. It follows that one should have now the possibility to characterize the macroscopic dynamics in a filled system by means of an effective noise temperature x related to the polymer chain mobility that depends on the interaction with the filler particles and on the strain. Further experimental work would be necessary to examine the role of the morphology of the nanoparticles used as filler. hysteresis would be expected for   . It follows that the hysteresis actually observed   results from a modification of the silica particle network due to movements of the particles under strain and observed by XPCS. For CB-40-c (Figure 10c), the logarithmic decrease of the strain with time in the first loading does not yield an asymptotic limit. In this case it appears that the values of   (DOWN) are very close to that obtained during the second loading (UP2). This feature is also observed for CB-20-c (Figure 11b). The effect of cross-linking on the evolution of  0 and   with strain is reported in Figure 12 for the hydroxylated silica filler (Si-OH-40). For the uncross-linked sample Si-OH-40-u,  0 goes through a maximum at close to 0.20. Such a feature has already been observed for uncrosslinked elastomers filled with silica or carbon black. It was shown to be coincident with the appearance of a "butterfly" in the SAXS patterns and attributed to the breaking and the reorganization of the colloidal particle network in the direction of the strain [62]. As discussed in the previous paragraph, the value of are the largest for Si-OH-40-u which involves a larger chain mobility making possible the reorganization of the filler particle network. These features suggest that the generalized CT equation describes the breaking (or deflocculation) of the silica network and that this mechanism becomes independent of the strain above 0.20. The fact that the asymptotic limit   remains nearly constant suggests that the system tends to the same final deflocculated state.

Analysis of

SUMMARY AND CONCLUSION

In this study, we have first analyzed the stress-strain curves obtained according the standards procedure. The un-conventional plotting of the curves in logarithmic coordinates, allowed us to reveal the same non-linear stress-strain behavior for all samples. This feature suggests the occurrence of plastic events at small strain which could be similar to what was observed by simulations about the effect of tension in polymeric glasses. Therefore, it would be valuable to revisit, in the future, the analysis of previous experimental data.

The largest part of the study is dealing with the analysis of the behavior of the tensile relaxation modulus ( ) at a given strain level reached either by a positive or a negative strain step. The leading parameter of this investigation was the filler-filler or the filler matrix interaction. We show that for all silica filled samples strained above 0.20 by a positive strain step, ( ) follows the generalized Chasset-Thirion equation. For the carbon black samples characterized by a strong filler-matrix interaction the decrease of ( ) is mainly logarithmic. These features indicate that the rate of decrease of ( ) after a positive strain step follows a power law characterized by an exponent . A similar feature is also observed for the rate of increase (recovery) of ( ) after a negative strain step for all samples and at all strain plateaus including . The values of m are always positive and the limiting case ( ) is never observed.

We have then remarked that the above power law features could be signatures of soft glassy rheology (SGR) phenomena characterized by a broad distribution of relaxation times and multiscale dynamical processes. It follows that (i) the exponent m involved in the Chasset-Thirion equation would be equivalent to the exponent involved in the SGR model and, (ii) the exponent could be considered as an effective noise temperature. In other words, one should have now the possibility to characterize the macroscopic dynamics in a filled system by means of an effective noise temperature x related to the polymer chain mobility that depends on the interaction with the filler particles and on the strain. We present several indications arising from the literature that may support this hypothesis. We think that the coming analysis of the series of HD-XPCS data obtained simultaneously with the mechanical measurements, and giving information about the "movements" of the filler particles, will be able to provide additional evidences of the pertinence of the SGR models in the field of filled elastomers.

Finally, the analysis of the limiting values (asymptotic limits) deduced from the fit of the data with the generalized CT equation, brings new information about the relative role of relaxation and plastic effect on the hysteresis (Mullins effect) observed in stress-strain measurements.

  This result allows us to introduce a generalized Chasset-Thirion equation writing ( ) [ ( ⁄ ) ] in which when is positive and when is negative. It is suggested that the physical meaning of the exponent m stemming from the generalized Chasset-Thirion equation can be assimilated to the exponent in the soft glassy rheology (SGR) model.

2. 1

 1 Samples. The elastomer matrix investigated is an Ethylene Propylene Diene Monomer (EPDM) rubber (Buna EP G 6850, Lanxess, Leverkusen, Germany). Three different fillers were used: carbon black N330 (Evonik Carbon Black GmbH), hydroxylated pyrogenic silica AEROSIL® 200 (Degussa) and AEROSIL® R 974. The latter is a hydrophobic fumed silica aftertreated with dimethyldichlorosilane based on the hydrophilic AEROSIL® 200 with a specific surface area close to 200 m 2 /g. The morphological characteristics of the fillers have been previously described in details[60][61][62]. The filled elastomers were prepared by mixing 20 or 40 g of filler with 100 g of rubber (i.e., 20 or 40 phr) and 3 g of dicumyl peroxide (the cross-linking agent) in a two-roll mill for 5 minutes at 50°C. Afterwards, the mixtures were cured in molds at 160°C under pressure during 10 minutes. The volume fraction of filler corresponding to 40 phr is close to 0.16 for the carbon black filled sample and to 0.14 for the silica ones. In both cases, the volume fraction is expected to be above the "percolation" threshold which is close to 0.10 for carbon black N330[63] and 0.02 for AEROSIL® 200[64,65]. At 20 phr, the volume fraction is slightly below percolation for carbon black (0.09) but probably above percolation, for silica (0.07). In the following, the names Si-OH-40-c, Si-OR-40-c and CB-40-c will refer to crosslinked samples with 40 phr of AEROSIL® 200, AEROSIL® R 974 and carbon black, respectively. Samples with 20 phr will be named Si-OH-20-c, Si-OR-20-c and CB-20-c. The uncross-linked AEROSIL® 200 filled sample will be named Si-OH-40-u.

Figure 1 .

 1 Figure 1. Stress-strain curves measured at a constant strain rate (0.133 s 1 ) (a) and (b), relative

Figure 2 ( 13 3. 3 .

 2133 Figure 2 (a, b, c). Incremental stress-strain curves. The filled black symbols correspond to the

Figure 3 . 4 . 1 .

 341 Figure 3. Comparison between stress (a) and tensile modulus (b) data obtained during the

Figure 4 .

 4 Figure 4. Time dependence of the tensile relaxation modulus E for  = 0.40. (a) series Si-OH-c

obtained for  0.60 and 0. 80 .

 80 The experimental results obtained for are plotted in Figure5. In all cases, the spread of the curves UP1, DOWN and UP2 is narrower than it was for (Figure4). For samples Si-OH-40-c, Si-OR-20-c and CB-40-c there are no DOWN plots because these samples were not strained above 0.60. For the other ones, the modified CT equation provides a good fit for the experimental data measured during the DOWN step. The CT equation fits well the UP1 and UP2 curves obtained for the two hydroxylated silica samples (Si-OH-20-c and Si-OH-40-c) and one of the hydrophobic silica sample (Si-OR-40-c). The values of the exponent m are indicated in brackets on figure 5.

Figure 5 .

 5 Figure 5. Time dependence of the tensile relaxation modulus E for  = 0.60. The values of m are

Figure 6 . 4 . 3 .

 643 Figure 6. Time dependence of the tensile relaxation modulus E for  = 0.80. The values of m are

Figure 7 .

 7 Figure 7. Time dependence of the tensile relaxation modulus E for  = 0.20. The values of m

Figure 8 .

 8 Figure 8. Time dependence of the tensile relaxation modulus E at different strains for the

Figure 9 .

 9 Figure 9. Comparison of the values of the exponent obtained at different strains.

  the asymptotic limits. The determination of an asymptotic equilibrium by means of the CT equation allows us to examine the Mullins effect in these systems. To this end, we have plotted on Figure 10 and 11 the values of the stress ( ) measured at (Figure 2) and the values of determined by means of the generalized CT equation for the three different fillers (40 phr). As expected, the second loading curves (UP2) are located between the first ones (UP1) and the unloading curves (DOWN). These figures show that the area of the hysteresis loop consisting of UP1 and DOWN curves decreases when going from Si-OH to Si-OR and to CB fillers (for 40 phr), i.e., when the filler matrix interaction increases. The magnitude of the hysteresis loop area seems similar for 40 phr (Figure 10c) and 20 phr (Figure 11b) of carbon black. In the case of hydroxylated silica, the decrease of the filler content seems to decrease the area of the hysteresis loop. This qualitative observation agrees with the Mullins effect model [80] based on dissipative friction phenomena due to internal sliding of the macromolecular chains and of the filler particles. It also intuitively agrees with the very different behavior of the filler particle velocity and relaxation investigated by HD-XPCS for Si-OH-40-c and CB-40-c (0.60) [51]. For Si-OH-40-c, the area of the hysteresis loop obtained by considering the asymptotic limits   determined by means of the CT equation is significantly reduced. If the hysteresis observed for  0 would be related only to relaxation effects, no

Figure 10 .

 10 Figure 10. Evolution of the initial stress  0 and of the asymptotic limit   determined by the CT

Figure 11 .

 11 Figure 11. Evolution of the initial stress  0 and of the asymptotic limit   determined by the CT

Figure 12 .

 12 Figure 12. Comparison between the evolution of the initial stress  0 and the asymptotic limit  

Table 1 .

 1 Parameters

	0.4	E(0) (MPa)	E  (MPa)	m	  (s)
		UP1	1.226	0.811 	0.080 	(2.0  4
	Si-OH-20-c	DOWN #	0.681			(2  7
		UP2	0.973			(1.4  6
		UP1	2.010	1.364 	0.182 	0.0143 
		DOWN #	0.744	1.000 	0.163 	(4.5  4
	Si-OH-40-c
		UP2	1.480	1.047 	0.111 	(2.7  4
		UP1	0.983	0.70 	0.054 0.003	(4.5  8
	Si-OR-20-c	DOWN #	0.735	0.803 	0.180 	(2.4  6
		UP2	0.922	0.786 		(1.4  7
		UP1	1.246	0.871 	0.086 	(3.1  5
	Si-OR-40-c	DOWN #	0.664	0.853 	0.094 	(7  8
		UP2	1.038	0.856 	0.122 	(2.7  6
		UP1	0.940	0.74 	0.038 	 15
	CB-20-c	DOWN #	0.727		0.084 	 13
		UP2	0.853		0.304 	(1.1 0.210 4
		UP1 #	1.514	E 0 =1.63 0.01 	-0.0740.003 	 15
	CB-40-c	DOWN #	1.209	1.273 0.001	0.117 0.004	(1.3 0.410 7
		UP2	1.411	1.19  0.02	0.053 0.005	 14

deduced from the fit of the data shown in Figure

4

by Equations

[START_REF] Kumar | [END_REF]

. The sign # indicates in Equation (

5

).

ACKNOWLEDGEMENTS

Use of the APS was supported by the DOE, Office of Basis Energy Sciences, under Contract No. W-31-109-Eng-38. The authors thank René Jurk (IFP Dresden) for the preparation of filled elastomers and for the measurement of the tensile stress-elongation curves.

For Table of Contents use only


Philippe Steer
email: philippe.steer@univ-rennes1.fr

Martine Simoes

Rodolphe Cattin

J Bruce

H Shyu

Erosion influences the seismicity of active thrust faults

Assessing seismic hazards remains one of the most challenging scientific issues in Earth sciences. Deep tectonic processes are classically considered as the only persistent mechanism driving the stress loading of active faults over a seismic cycle. Here we show via a mechanical model that erosion also significantly influences the stress loading of thrust faults at the timescale of a seismic cycle. Indeed, erosion rates of about ~0.1 to 20 mm.yr -1 , as documented in Taiwan and in other active compressional orogens, can raise the Coulomb stress by ~0.1 to ~10 bar on the nearby thrust faults over the inter-seismic phase. Mass transfers induced by surface processes in general, during continuous or short-lived and intense events, represent a prominent mechanism for inter-seismic stress loading of faults near the surface. Such stresses are probably sufficient to trigger shallow seismicity or promote the rupture of deep continental earthquakes up to the surface.

INTRODUCTION 2

The evolution of the Earth's topography is dictated by the interactions between tectonics, climate and surface processes (i.e. erosion and sedimentation). Whether this evolution influences tectonic deformation during mountain building has been widely debated. It is now well accepted that surface evolution can drive the localization and intensity of tectonic deformation over geological times [START_REF] Dahlen | Brittle frictional mountain building: 1. Deformation and mechanical energy budget[END_REF][START_REF] Willett | Orogeny and orography: The effects of erosion on the structure of mountain belts[END_REF][START_REF] Beaumont | Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[END_REF] (1-10 Myr). At intermediate time scales (10 kyr -1 Myr), erosion and the subsequent isostatic rebound can favour slip along specific fault planes [START_REF] Cattin | Modeling mountain building and the seismic cycle in the Himalaya of Nepal[END_REF][START_REF] Maniatis | Slip acceleration on normal faults due to erosion and sedimentation-Results from a new three-dimensional numerical model coupling tectonics and landscape evolution[END_REF][START_REF] Calais | Triggering of New Madrid seismicity by late-Pleistocene erosion[END_REF][START_REF] Vernant | Erosion-induced isostatic rebound triggers extension in low convergent mountain ranges[END_REF] . However, the link between surface processes and the stress loading of faults during the seismic cycle (0.1-1 kyr), and in turn the associated deformation mechanisms, remains unsubstantiated.

Faults represent the main mechanical discontinuities of the elastic-brittle Earth's upper crust. They accommodate tectonic deformation by slipping, mostly during earthquakes [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF] . These seismogenic faults are rooted down dip in viscous shear zones [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Doglioni | Role of the brittle-ductile transition on fault activation[END_REF][START_REF] Cowie | Viscous roots of active seismogenic faults revealed by geologic slip rate variations[END_REF] . It is generally accepted that, during the inter-seismic phase (i.e. prior to an earthquake), continuous viscous flow in these deep shear zones leads to the elastic stress loading of active faults closer to failure, and that during the co-seismic phase (i.e. during an earthquake), failure and slip occur along the previously locked fault planes, followed by post-seismic stress relaxation [START_REF] Cattin | Modeling mountain building and the seismic cycle in the Himalaya of Nepal[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF] . Fault failure is commonly defined by the mean of the Coulomb stress change, CFF= η + µ'•ζ n , a function of the fault effective friction µ', the shear η (positive in the direction of slip) and normal ζ n (positive if the fault is unclamped) stress changes [START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF] . Earthquakes can be triggered by tectonic stresses, but also by Coulomb stresses due to episodic and short-lived events such as hydrologic [START_REF] Bettinelli | Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology[END_REF] or snow loading [START_REF] Heki | Seasonal modulation of interseismic strain buildup in northeastern Japan driven by snow loads[END_REF] , nearby earthquakes [START_REF] Reasenberg | Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake[END_REF][START_REF] King | Static stress changes and the triggering of earthquakes[END_REF][START_REF] Stein | The role of stress transfer in earthquake occurrence[END_REF][START_REF] Delescluse | April 2012 intra-oceanic seismicity off Sumatra boosted by the Banda-Aceh megathrust[END_REF] and slow-slip events [START_REF] Segall | Earthquakes triggered by silent slip events on Kilauea volcano, Hawaii[END_REF] .

Here, we show that surface processes significantly contribute to the Coulomb stress loading of thrust faults during the seismic cycle. To illustrate and then demonstrate our point, we consider a mountain range in Taiwan where the rates of erosion [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] and tectonic deformation [START_REF] Yu | Velocity field of GPS stations in the Taiwan area[END_REF][START_REF] Simoes | Investigating the kinematics of mountain building in Taiwan from the spatiotemporal evolution of the foreland basin and western foothills[END_REF] are extremely high and amongst the best-documented in the world. We then investigate how erosion influences the stress loading of thrust faults using a simple model for the seismic cycle.

RESULTS

Coulomb stress changes induced by erosion in Taiwan

Our first model quantifies the Coulomb stress change ΔCFF generated by erosional unloading, as constrained from fluvial suspended sediment load measured over the 30 yr prior to the 1999 M w 7.6 Chi-Chi earthquake in central Taiwan [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] (Fig. 1 and Methods). The 3D velocity field , strain rate ̇ and stress rate ̇ tensors induced by erosion are computed in an elastic halfspace using a Boussinesq approach (Fig. 1C and Methods). We use simplified geometries for active thrust faults located in the foothills of Taiwan 22 (see Methods) and assume a dip angle α of 30°, a 15 km deep brittle-ductile transition [START_REF] Shyu | Neotectonic architecture of Taiwan and its implications for future large earthquakes[END_REF] and an effective friction µ' of 0.5 to compute Coulomb stress changes per unit time (or loading rates) ΔCFF due to erosional unloading on these faults. We find a maximum value of ~4 x 10 -3 bar.yr -1 for the Coulomb stress change ΔCFF induced by erosional unloading on the Liuchia fault system (number 8 on Figure 1) in southwestern Taiwan. Despite a low topographic relief, this area has the highest erosion rates documented in Taiwan (up to 24 mm.yr -1 ), which are proposed to be controlled locally by a low substrate strength, a high storminess and a high seismic moment release rate [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] . However, most of the thrust faults located in the foothills still display a significant ΔCFF of ~0.5 x 10 -3 bar.yr -1 , including the Chelungpu fault (number 3) that ruptured during the Chi-Chi earthquake.

Integrated over a seismic cycle duration [START_REF] Chen | Late Holocene paleoseismicity of the southern part of the Chelungpu fault in central Taiwan: Evidence from the Chushan excavation site[END_REF] of ~500 yr, a ΔCFF of 0.5 to 4 x 10 -3 bar.yr -1 due to erosional unloading gives a net Coulomb stress change of 0.25-2.0 bar. Similar values of Coulomb stress change are documented elsewhere to contribute significantly to the stress loading and dynamics of active faults [START_REF] Bettinelli | Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology[END_REF][START_REF] Heki | Seasonal modulation of interseismic strain buildup in northeastern Japan driven by snow loads[END_REF][START_REF] Reasenberg | Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake[END_REF][START_REF] King | Static stress changes and the triggering of earthquakes[END_REF][START_REF] Stein | The role of stress transfer in earthquake occurrence[END_REF][START_REF] Delescluse | April 2012 intra-oceanic seismicity off Sumatra boosted by the Banda-Aceh megathrust[END_REF][START_REF] Segall | Earthquakes triggered by silent slip events on Kilauea volcano, Hawaii[END_REF] . This suggests that erosional unloading can significantly influence the short-term dynamics of faults.

Erosional unloading modifies the Coulomb stress change on a fault plane in two ways: 1) it decreases the normal stress and unclamps the fault, and 2) it increases the tangential stress (fig. 2a). The increment of stress on a fault plane is proportional to the amount of erosion, but decreases with the square of the distance r between the fault plane and where erosion occurs (fig. 2b). Therefore, the amplitude of erosional Coulomb stress loading on a fault is sensitive 1) to the effective friction µ', which modulates the effect of erosion on the normal stress, and 2) to the fault dip angle α, which decomposes the stresses into fault normal and tangential components and controls the distance of the fault plane to the surface. For instance, a higher effective friction µ' of 0.8 and a lower dip angle α of 15° therefore result in increasing the induced ΔCFF up to ~1

x 10 -2 bar.yr -1 for the Liuchia fault system (Fig. 2c). In addition, the stresses modelled here are invariant with the Young modulus of the material as we are considering a linear elastic material subjected to a surface pressure load (and not to a surface displacement).

Stresses induced by erosion during the seismic cycle

The above computations consider inter-seismic erosion rates calculated from data acquired during the 30 yr preceding the Chi-Chi earthquake [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] . Even though its amplitude relative to co-seismic rock uplift is debated, co-and post-seismic erosional unloading represents a major contribution to erosion in seismic areas [START_REF] Parker | Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth[END_REF][START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF][START_REF] Li | Seismic mountain building: Landslides associated with the 2008 Wenchuan earthquake in the context of a generalized model for earthquake volume balance[END_REF] . In mountain belts with hillslopes close to failure, coseismic ground motion and acceleration can induce a significant amount of landslides [START_REF] Meunier | Topographic site effects and the location of earthquake induced landslides[END_REF] . The sediments produced by these landslides are then transported by rivers mainly during subsequent floods. This post-seismic landscape relaxation phase has a documented potential duration [START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF] of years to decades, one order of magnitude shorter than a complete seismic cycle. Therefore, the contribution of co-and post-seismic erosion to the stress loading of active faults also needs to be evaluated.

To assess the relative contribution of inter-seismic erosion, co-/post-seismic erosion and tectonics to the Coulomb stress loading of faults, we develop a simple model of the seismic cycle that accounts for the effect of both erosion and tectonics (Fig. 3 and Methods). We assume a steady-state landscape over the seismic cycle, i.e. rock uplift rates are balanced by erosion rates Ė over this time scale. Because of the response-time of the geomorphic system to climate or tectonic perturbations and because of their stochastic properties [START_REF] Lague | Discharge, discharge variability, and the bedrock channel profile[END_REF] , this assumption is probably not valid in most settings. However, it offers a simple and self-consistent approach for modelling first-order surface processes during the seismic cycle. A Boussinesq approach is used to compute , ̇ and ̇, while the tectonic stresses and uplift (and therefore erosion) are calculated using dislocations embedded in an elastic half-space [START_REF] Meade | Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space[END_REF] . The effects of tectonic deformation during the inter-and co-seismic phases are accounted for by slip on a deep shear zone and on a shallow brittle fault, respectively [START_REF] Doglioni | Role of the brittle-ductile transition on fault activation[END_REF][START_REF] Vergne | On the use of dislocations to model interseismic strain and stress build-up at intracontinental thrust faults[END_REF] . This seismic cycle model is valid when considering a fault that is fully locked during the inter-seismic phase, as it is proposed for the thrust faults located in the western foothills of Taiwan [START_REF] Loevenbruck | Seismic cycle in Taiwan derived from GPS measurements[END_REF][START_REF] Dominguez | Horizontal coseismic deformation of the 1999 Chi-Chi earthquake measured from SPOT satellite images: Implications for the seismic cycle along the western foothills of central Taiwan[END_REF][START_REF] Hsu | A two-dimensional dislocation model for interseismic deformation of the Taiwan mountain belt[END_REF] .

For comparison with faults in the foothills of Taiwan, we define a reference model with a fault trace length of 80 km and a dip angle of 30°, while keeping all the other mechanical properties identical. We also impose a slip velocity V inter of 40 mm.yr -1 on the shear zone during the inter-seismic phase, and V co of 40 mm.yr -1 on the associated brittle fault during the coseismic phase. This model setup provides only a rough approximation of the seismic cycle at the scale of the whole western foothills of Taiwan. Indeed, deformation is partitioned between several active thrust faults in this area that are probably rooted down at depth into a single decollement with a total slip of ~40 mm.yr -1 (refs. 21, 33). In addition, because our goal is to quantify the co-and post-seismic erosional unloading rates during the landscape relaxation phase following large earthquakes [START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF] , we compute co-seismic slip velocity averaged over the seismic cycle rather than co-seismic instantaneous displacement. Note that these two approaches are strictly equivalent in an elastic model.

In our modelling, inter-and co-seismic rock uplift (and erosion) rates are similar and up to ~20 mm.yr -1 (Fig. 3). We assume no time modulation of inter-and co-seismic erosion and both modeled erosion rates are applied over the entire duration of one seismic cycle. Despite similar erosion rates, co-seismic erosional unloading induces fault Coulomb stress loading ΔCFF E-CO of up to ~8 x 10 -2 bar.yr -1 at very shallow depth (< 1 km), which is about 30 times the maximum stress loading induced by inter-seismic erosion ΔCFF E-INTER . Indeed, the maxima of co-seismic uplift of the surface (and erosion) occurs over the shallow portion of the fault, therefore at a shorter distance to the fault plane than the maxima of inter-seismic erosion (Fig. 3).

Despite a rapid decrease with depth, ΔCFF E-CO is still greater than ~0.5 x 10 -2 bar.yr -1 at a depth of 5 km. Coulomb stress loading ΔCFF E-INTER induced by inter-seismic erosion displays two local maxima of ~0.3 x 10 -2 bar.yr -1 , one located on the deeper part of the fault underneath the maximum of inter-seismic erosion, and the second one close to the fault tip due to its proximity with the surface.

We then compare fault Coulomb stress loading rates induced by erosion to those induced only by tectonics during the inter-seismic phase ΔCFF T-inter . In terms of amplitude, Coulomb stress loading rates due to tectonics are up to two or four orders of magnitude greater than those related to erosion, in particular for the deeper part of the fault. However, at shallower depths (< 5 km), Coulomb stress loading rates due to co-seismic erosion and tectonics are of the same order of magnitude, and the ratio ΔCFF E-co /ΔCFF T-inter even reaches ~20 close to the surface (Fig. 3).

At the contrary, Coulomb stress loading rates associated with inter-seismic erosion, which is maximum on the deeper and shallower part of the fault, do not dominate tectonic stresses. The ratio ΔCFF E-inter /ΔCFF T-inter only reaches ~0.1 at intermediate depths (5-10 km) and ~0.6 close to the surface. Because the upper crust displays a very long stress relaxation time associated to high effective viscosities, these Coulomb stress changes induced by erosion can be accumulated over the time scale of a seismic cycle (~500 yr). On the other hand, increasing the Young modulus from the reference value (10 GPa) by a factor of 2 (20 GPa) or 5 (50 GPa) results in increasing the Coulomb stress change induced by tectonics ΔCFF T-inter by a factor of 2 or 5 and decreasing the ratios ΔCFF E-inter /ΔCFF T-inter and ΔCFF E-co /ΔCFF T-inter by the same factor (Fig. 3h). However, these results still demonstrate 1) that erosion can contribute significantly to thrust fault stress loading during the seismic cycle, and 2) that erosion can even be one of the dominant stress loading mechanisms for the shallower parts of thrust fault planes.

Model sensitivity analysis

To assess the sensitivity of our results to the model parameters, we design a set of models similar to the reference model but with varying values of the effective friction µ' (0.1 to 0.9), the Young modulus E (10, 20 and 50 GPa) and the fault and shear zone dip angle α (15 to 45°). For the sake of simplicity we keep the depth of the brittle-ductile transition at 15 km, which in turn implies that the surface area of the brittle fault increases when decreasing α. Using this modelling approach, erosion rates during the inter-seismic Ė inter and the co-seismic Ė co phases are only sensitive to α (Fig. 4a-b). While Ė inter remains approximately constant around 13 mm.yr -1 , Ė co increases from 11 to 22 mm.yr -1 when α increases from 15 to 45°.

The resulting Coulomb stresses ΔCFF E-inter and ΔCFF E-co on the fault plane are in turn sensitive to 1) the dip angle α, which controls both the distribution and amplitude of Ė inter and Ė co and the distance between the Earth surface and the fault plane, and 2) the effective friction µ' by amplifying the influence of normal stress on the Coulomb stress (Fig. 4c-d). The maximum values of ΔCFF E-inter and ΔCFF E-co obtained on the fault plane show similar distributions in the parameter space. ΔCFF E-inter is minimum (~1 x 10 -3 bar.yr -1 ) for a low effective friction (µ'0.2) combined to a high or low dip angle (15° or 45°), while it increases when increasing µ' and reaches a maximum (~2.7 x 10 -3 bar.yr -1 ) for α around 25°. ΔCFF E-co is minimum (~1.5 x 10 -2 bar.yr -1 ) for a low effective friction (µ'0.2) combined to a high or low dip angle (15° or 45°), while it increases when increasing µ' and reaches a maximum of ~4 x 10 -2 bar.yr -1 for α around Therefore, Coulomb stresses induced by erosion during the seismic cycle represent a significant contribution to fault stress loading, even though their amplitude depends on the properties of the fault (dip angle, effective friction) and of the medium (Young modulus).

DISCUSSION

Based on erosion data from Taiwan [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] , we have demonstrated that the elastic Coulomb stresses induced by erosion are of the order of ~0.5 x 10 -3 bar.yr -1 on the thrust faults located in the western foothills and reach a maximum of ~4 x 10 -3 bar.yr -1 on the Liuchia fault. These results are consistent with the outcomes from a simple model of the seismic cycle of a thrust fault that accounts for the effect of both erosion and tectonics using fault properties and a slip velocity close to the ones inferred for Taiwan [START_REF] Yu | Velocity field of GPS stations in the Taiwan area[END_REF][START_REF] Simoes | Investigating the kinematics of mountain building in Taiwan from the spatiotemporal evolution of the foreland basin and western foothills[END_REF][START_REF] Shyu | Neotectonic architecture of Taiwan and its implications for future large earthquakes[END_REF] . Coulomb stresses induced by inter-seismic ΔCFF E-inter and co-seismic ΔCFF E-co erosion and averaged over the duration of a seismic cycle (~500 yr) reach values of up to ~3 x 10 -3 bar.yr -1 and ~8 x 10 -2 bar.yr -1 , respectively. On the shallower part of thrust faults (< 5 km deep), the ratio of the Coulomb stresses induced by coseismic erosion ΔCFF E-co to the ones induced by tectonic loading ΔCFF T-inter is about equal to 1 for a Young modulus E=10 GPa (~0.2 for E=50 GPa) and even reach a maximum of ~20 closer to the surface (~4 for E=50 GPa). In addition, assuming that co-seismic erosion happens only during a period 10 to 100 times shorter [START_REF] Hovius | Prolonged seismically induced erosion and the mass balance of a large earthquake[END_REF][START_REF] Howarth | Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault[END_REF] (~5-50 yr) than the complete seismic cycle 23 (~500 yr), ΔCFF E-CO increases up to 80 x 10 -1 to 8 bar.yr -1 and the ratio ΔCFF E-co /ΔCFF T-inter to 10 or 100 for E=10 GPa (2 to 20 for E=50 GPa) during the first ~5-50 years following a large earthquake. Large earthquakes with a potential negative mass balance (i.e. erosion greater than uplift), such as the Wenchuan earthquake [START_REF] Parker | Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth[END_REF][START_REF] Li | Seismic mountain building: Landslides associated with the 2008 Wenchuan earthquake in the context of a generalized model for earthquake volume balance[END_REF] , could induce even higher rates of ΔCFF E-co than those predicted by a steady-state model. Our modelling approach imposes that co-seismic erosion is maximum close to the fault trace and above the shallower part of thrust faults (Fig. 3). This result contrasts with most of the observed distribution of earthquake-triggered landslides, with a maximum of landslide density close to the epicentral area and therefore generally above the deeper part of thrust faults [START_REF] Meunier | Topographic site effects and the location of earthquake induced landslides[END_REF] . Therefore, depending on the location of large earthquake hypocenters along the fault plane, our estimates for the contribution of co-seismic erosion to fault stress loading might likely be overestimated. However, our results emphasize that short-lived and intense erosional events associated with efficient sediment transport, such as typhoons, could suddenly increase the Coulomb stress of underlying faults. On longer time-scales, climatic changes or transition from fluvial-to glacial-dominated surface processes could also lead to high transient erosion rates and therefore to transient increases of the fault stress loading due to erosion [START_REF] Cattin | Modeling mountain building and the seismic cycle in the Himalaya of Nepal[END_REF][START_REF] Maniatis | Slip acceleration on normal faults due to erosion and sedimentation-Results from a new three-dimensional numerical model coupling tectonics and landscape evolution[END_REF][START_REF] Calais | Triggering of New Madrid seismicity by late-Pleistocene erosion[END_REF][START_REF] Vernant | Erosion-induced isostatic rebound triggers extension in low convergent mountain ranges[END_REF] . The mechanism proposed in this study is limited neither to convergent settings nor to erosion only, as sediment deposition on the hanging wall of normal faults could also lead to a significant increase in Coulomb stress [START_REF] Maniatis | Slip acceleration on normal faults due to erosion and sedimentation-Results from a new three-dimensional numerical model coupling tectonics and landscape evolution[END_REF] .

Moreover, some less active areas, such as intra-continental faults or old orogens still experience intense erosion and episodic seismic activity [START_REF] Calais | Triggering of New Madrid seismicity by late-Pleistocene erosion[END_REF][START_REF] Vernant | Erosion-induced isostatic rebound triggers extension in low convergent mountain ranges[END_REF] . In the absence of major tectonic deformation, surface processes could significantly contribute to the stress loading of faults in these areas, even when considered independently of the stresses induced by isostatic rebound.

In summary, our results demonstrate that surface processes represent a significant contribution to the Coulomb stress loading of faults during the seismic cycle. In terms of deformation, these additional stresses on the shallower part of fault planes can induce and trigger shallow earthquakes, as illustrated by the seismicity triggered by the large 2013 Bingham Canyon mine landslide [START_REF] Pankow | Massive landslide at Utah copper mine generates wealth of geophysical data[END_REF] , or potentially favour the rupture of large deeper earthquakes up to the surface as, for instance, during the Chi-Chi earthquake [START_REF] Cattin | Why does the co-seismic slip of the 1999 Chi-Chi (Taiwan) earthquake increase progressively northwestward on the plane of rupture?[END_REF] . This offers new perspectives on the mechanisms influencing stress transfers during the seismic cycle, as well as on seismic hazard assessment in areas experiencing rapid erosion. More generally, Coulomb stress loading of faults induced by surface processes over short time scales provides an additional positive feedback between climate, surface processes and tectonics.

METHODS

Seismic cycle model

The deformation model computes the velocity field , strain ̇ and stress ̇ rate tensors induced by surface erosion in a 3D elastic half-space based on the Boussinesq approximation [START_REF] Boussinesq | Équilibre d'élasticité d'un sol isotrope sans pesanteur, supportant différents poids[END_REF] . With respect to this approximation, we assume that the model surface is horizontal (Fig. 1). The effect of such assumption is here limited as the topography of the western foothills of Taiwan is globally lower than 1 km. The model is discretized by cubic cells with a 100 m-resolution, with a Young modulus of E=10 GPa, a Poison ratio of ν=0.25, and a rock density of ρ=2800 kg.m -3 . Velocity, stresses and strain induced by tectonics are simulated by the mean of triangular dislocations [START_REF] Meade | Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space[END_REF][START_REF] Vergne | On the use of dislocations to model interseismic strain and stress build-up at intracontinental thrust faults[END_REF] accounting for the slip 1) of a viscous deep shear zone during the inter-seismic phase and of 2) a frictional fault located in the shallow elastic-brittle crust during the co-seismic phase. The imposed averaged tangential velocities of the shear zone V inter and of the fault V co are both equal to 40 mm.yr -1 . Coulomb stress changes are then computed by projecting the stresses due to surface processes and to tectonics on the fault plane that is discretized with a resolution of 100 m. The extent of fault planes and the dip angle data of the western foothills of Taiwan were simplified from ref. 22. Each fault trace was simplified to a line segment that best reproduces the real fault trace geometry.

Elastic Boussinesq model

( ) (  ) ( ) (  )
Because the model is linear and elastic, the total displacement, stress and strain components for any distribution of surface load are then computed by summation of the displacement, stress and strain components obtained for each individual point load.

Taiwan erosion rates

Erosion rates in Taiwan were calculated from fluvial suspended sediment load measured over the 30 yr prior to the 1999 M w 7.6 Chi-Chi earthquake, considering that the total fluvial sediment load contains 70% of suspended load [START_REF] Doglioni | Role of the brittle-ductile transition on fault activation[END_REF] (Fig. 1b). Erosion rates were calculated assuming that suspended sediment load, measured at 130 gauging stations, represents 70% of the river total sediment load, and that catchment-wide erosion rates correspond to the total sediment load divided by sediment density and by drainage area [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] . The smoothed erosion map of Figure 1 was then obtained using a circular averaging window with a radius of 30 km. The spatial resolution of the erosion map, even though coarse, still allows for resolving potential heterogeneities of ΔCFF induced by erosion between different faults and along each individual fault plane. Topography, simplified main thrust fault systems (red lines and numbers from 1 to 8) and location of the 1999 M w 7.6 Chi-Chi earthquake epicenter (red star) that ruptured the Chelungpu fault (number 3). b) Inter-seismic erosion rates prior to the Chi-Chi earthquake, as calculated from fluvial suspended sediment measurements with a 5-km grid resolution, smoothed at the catchment scale using a circular moving mean with 30-km diameter [START_REF] Dadson | Links between erosion, runoff variability and seismicity in the Taiwan orogen[END_REF] . c) Erosion induced

  stresses induced by tectonics ΔCFF T-inter are also sensitive to the Young modulus E, the ratios ΔCFF E-inter /ΔCFF T-inter and ΔCFF E-co /ΔCFF T-inter are in turn sensitive to α, µ' and E. Increasing the Young modulus from 10 to 20 or 50 GPa increases ΔCFF T-inter and decreases ΔCFF E-inter /ΔCFF T-inter and ΔCFF E-co /ΔCFF T-inter by a factor of 2 or 5,respectively. ΔCFF E-inter /ΔCFF T-inter remains lower than 1 only for all the models tested, independent of the Young modulus E. At the contrary, ΔCFF E-co /ΔCFF T-inter displays a large domain in the parameter space with values greater than 1 (i.e. with at least one element of the fault plane dominated by stresses induced by erosion). For E=10 GPa, only models with low µ' (<0.5) and high α (>35°) are dominated by tectonic stresses (ΔCFF E-co /ΔCFF T-inter <1), while for E=50 GPa, most of the models with α greater than 20-30° are dominated by tectonic stresses.

Figure 1 .

 1 Figure 1. Faults stress loading rates induced by erosion in the foothills of Taiwan. a)

We are grateful to Patience Cowie, Dimitri Lague, Philippe Davy and Frédéric Boudin for discussions. This study was funded by the CNRS-INSU (Centre National de la Recherche Scientifique -Institut National des Sciences de l'Univers, France), the ORCHID program (Partenariats Hubert Curien, France-Taiwan), the EROQUAKE ANR project (Agence Nationale de la Recherche, France) and the Université Rennes 1. This is IPGP contribution number 3569.

We here consider the displacements, stress and strain components generated by a point load F at the surface of a 3D semi-infinite elastic solid of coordinates x, y and z, with z being positive downward. Let For z0, the displacement components are [START_REF] Boussinesq | Équilibre d'élasticité d'un sol isotrope sans pesanteur, supportant différents poids[END_REF][START_REF] Westergaard | Theory of elasticity and plasticity[END_REF] :

)

Assuming infinitesimal deformation, the symmetric Cauchy strain tensor ε is then obtained by differentiating the displacement vector, ( ⁄ ⁄ ) ⁄

. For an isotropic medium, the stress components are then given by the following equation, ,

where δ ij is the Kronecker delta. The 6 stress components are [START_REF] Boussinesq | Équilibre d'élasticité d'un sol isotrope sans pesanteur, supportant différents poids[END_REF][START_REF] Westergaard | Theory of elasticity and plasticity[END_REF] :

) )

AUTHOR CONTRIBUTIONS

P.S. analyzed the data and performed the modelling. All authors contributed equally to the design of the study and to the writing of the paper.

ADDITIONAL INFORMATION

The authors declare no competing financial interests. Supplementary information accompanies this paper. Correspondence and requests for materials should be addressed to P.S.

Coulomb stress loading rates ΔCFF calculated on the fault planes. The horizontal solid black lines represent scale bars of 50 km.


Marc Fouet

Rémi Manczak

Rémi Courson

Marie-Charline Blatché

Karine Reybier

Anne Marie Gué

Multilevel (3D) microfluidic technology for an innovative magnetic cell separation and couting platform

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

1 CNRS, LAAS -2 PHARMADEV, IRD -3 Université de Toulouse -Toulouse, FRANCE. * gue@laas.fr

3D Microfabrication

Magnetic Separation and Trapping

The use of super paramagnetic microbeads that can be functionalized with specific biomolecules (e.g. antibodies) enables micro manipulation, separation and trapping of cells. A magnetic field is created with microcoils, which can be also coupled to magnets to amplify the force (hybrid systems).

Advantages: fast deposition of thick layer and high accuracy alignment. ! Metal deposition and electroplating are used for the coils fabrication, allowing a high resolution (5 x 5 µm wires) and 3D microfluidic channels integration on the same chip.

Laminator

We showed tagged monocytes can be separated from a sample using a set of coils for specific parameters (flow rate, electrical intensity…).

Mixing stage

3D microfluidics -baker's transformation

Separation stage

Vertical magnetic sorting

Detection stage

Electrical impedance measurements

Microfluidic Device

Three dimensional microfluidic structures can be created with dry film lamination techniques. Mixing applications have been investigated, especially with designs dedicated for cells and beads tagging reactions. High integration could be achieved with other functions like magnetophoresis and specific impedance electrical detection: the device is then placed in a chip holder and connected to fluidic controllers and electronics.

Microfluidic chip with 10 coils embedded 18 mm

An interesting application is the diagnosis of infectious diseases based on monocytes subpopulations. An immunotagging is performed: magnetic beads are attached to monocytes.

Tagged Monocyte

Magne:c beads Inlet 1: Sample A serie of microcoils then separate them from the sample, and the different subpopulations are detected with an immunosensor based on electrical impedance measurements.

Inlet 2: Buffer

Introduction

Currently, the technique for the quantitative detection of cells is flow cytometry. This technique has the advantage of being sensitive and reliable but is expensive, time consuming and not suited to both routine screening and point-of-care diagnostics. Miniaturized cell separation devices offer many advantages such as the use of small volumes, portability and low cost.

We propose a new concept of device which, by combining 3D fluid engineering and localized magnetic actuation, enables the full integration of cell tagging, magnetic separation and cell counting in a single device. The labs on chip are manufactured by laminating commercially available photosensitive dry film that fits microfluidic requirements and gives the possibility to build easily 3D microfluidic systems.

We show we can tag efficiently THP1 monocytes and subsequently sort them through magnetic trapping on integrated micro-coils. The separation efficiency is studied at different flow rates. Cell counting capacity, evaluated by using nonfaradic impedance spectroscopy revealed that the cell trapping is selective, depending on the specific antibody grafted and quantitative with the range of detection being 1000 to 30000 infected cells. This range of detection is consistent with the targeted application.

Conclusion

We demonstrated a new concept of devices, which by combining 3D fluid engineering and localized magnetic actuation enables the full integration of a cell tagging, magnetic separation device and cell counting. This approach, whose main drawback is the limitation of flow (and therefore volumes) which can be handled, can be significantly improved by adding external permanent magnets for the generation of high and homogeneous magnetic fields. Such hybrid systems open up the way to more efficient systems but also to the implementation of new kind of functionalities: focusing, shift register, trapping array, manipulation of single cell or chain of clusters.

  13 , R. Manczak 23 , R. Courson 13 , C. Blatché 13 , K. Reybier 23 and A.-M. Gué 13*

  Mask Aligner o Lamination of a dry film photoresist (DFR) DF-1000 series 1 onto a substrate. o Structuration of this DFR by photolithography. o Lamination of a second layer of DFR onto the first structured one and structuration by photolithography. o Repeating these steps allows to create a

  SEM picture of microfabricated spiral coil Impact of flow rate on beads trapping and separa:on efficiency Trapped (magne:cally tagged) THP1 monocytes on a coil Surface modification and Impedance detection o Functionalization of gold microelectrodes with antibody specific to the antigen present on the target cell : after activation of the first layer of SAMs, the electrode surface is modified by the addition of the protein G (saturated by adding BSA) before being coated with specific antibody. o Characterization of electrode modification and cell attachment by impedance spectroscopy. o Determination of the intrinsic limit of detection (LOD) : minimal detectable number of captured cells on the sensor surface. o Optimization of the microelectrode design to maximize the signal to noise ratio and hence the LOD. o Application of a scale cell concentration to determine design sensitivity.