Riad Hassani

Patrick Hild

Ioan Ionescu
email: ionescu@univ-savoie.fr

Sufficient conditions of non-uniqueness for the Coulomb friction problem

Keywords: Coulomb friction, elastostatics, non-uniqueness, eigenvalue problem, nite element approximation

We consider the Signorini problem with Coulomb friction in elasticity. Sucient conditions of nonuniqueness are obtained for the continuous model. These conditions are linked to the existence of real eigenvalues of an operator in a Hilbert space. We prove that, under appropriate conditions, real eigenvalues exist for a non-local Coulomb friction model. Finite element approximation of the eigenvalue problem is considered and numerical experiments are performed.

INTRODUCTION

Many applications in solid mechanics involve contact problems between elastic structures. Very often, the Coulomb friction model is chosen in the modelling of the contact phenomena. From a mathematical point of view, the Coulomb frictional contact problem in (continuum) elastostatics causes considerable diculties and is still open. From a mechanical point of view, there is special interest in the investigation of uniqueness of the solutions. The aim of this paper is to shed some light on this question.

The variational formulation of the continuous problem in elastostatics was given by Duvaut and Lions [START_REF] Duvaut | Les in equations en m ecanique et en physique[END_REF]. The rst existence results were obtained by Ne cas et al. in Reference [START_REF] Ne Cas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF] for an innite elastic strip. Thereafter, existence results were obtained for an arbitrary domain [START_REF] Jaru | Contact problems with bounded friction. Coercive case[END_REF][START_REF] Kato | Signorini's problem with friction in linear elasticity[END_REF][START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF]. In all these papers, the existence results hold for small friction coecients and the uniqueness is not discussed. The so-called non-local Coulomb frictional models mollifying the normal stresses were introduced by Duvaut [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF] and developed in References [START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF][START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF]. The smoothing map used in the non-local friction model allows to obtain existence results for any friction coecient. Moreover, uniqueness results can also be established if the friction coecient is small enough [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF][START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF][START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF]. The same type of result (existence for any friction coecient and uniqueness for small friction coecients) was obtained by Klarbring et al. [START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | On friction problems with normal compliance[END_REF] in the case of the normal compliance model, introduced by Oden and Martins [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF][START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF]. Finally, let us remark that the sucient conditions for uniqueness (small friction) given in all the abovementioned papers are not completed by neither sucient conditions for non-uniqueness nor by examples of non-uniqueness.

The discrete (nite element) problem, associated with the continuous static Coulomb friction problem, always admits a solution [START_REF] Haslinger | Approximation of the Signorini problem with friction, obeying the Coulomb law[END_REF][START_REF] Haslinger | Least square method for solving contact problems with friction obeying Coulomb's law[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF] and it is unique if the friction coecient is small enough. Moreover, a convergence result of the nite element model towards the continuous model was established by Haslinger [START_REF] Haslinger | Approximation of the Signorini problem with friction, obeying the Coulomb law[END_REF]. In the nite dimensional context, numerous studies using truss elements have led to examples of non-uniqueness. The early work concerning non-uniqueness was done by Janovsk y [START_REF] Lions | Catastrophic features of Coulomb friction model[END_REF] and was followed by Klarbring who constructed a concrete example of non-uniqueness involving a spring system in Reference [START_REF] Klarbring | Examples of non-uniqueness and non-existence of solutions to quasistatic contact problems with friction[END_REF]. Let us mention that Alart considered the general framework of nite dimensional systems. He obtained in Reference [START_REF] Alart | Crit eres d'injectivit e et de surjectivit e pour certaines applications de R n dans lui même: application a la m ecanique du contact[END_REF] abstract necessary and sucient conditions for uniqueness. In elastostatics, all the uniqueness results in the nite dimensional context are valid for friction coecients lower than a critical value (in the quasi-static case, this does not hold according to the counterexample of Ballard [START_REF] Ballard | A counter-example to uniqueness in quasi-static elastic contact problems with small friction[END_REF]). This critical value always depends on the number of degrees of freedom (on the mesh size when nite elements are used or on the dimension of the system in the case of truss elements). Since this critical value vanishes as the number of degrees of freedom increases, we cannot deduce any result for the continuous problem. Furthermore, the examples of non-uniqueness are specic to the nite dimensional system such that no continuous non-uniqueness example can be constructed from it.

The aim of this paper is to give simple sucient conditions for non-uniqueness of the solution to the continuous Coulomb friction problem which are related to the analysis of an eigenvalue problem. The spectral approach developed here is dierent from the widespread xed-point technique used in the search of solutions to the Coulomb friction problem. To our knowledge, this is the rst preliminary result dealing with non-uniqueness conditions in the continuous context.

After the statement of the problem, we give in Section 3 sucient conditions for nonuniqueness. They deal with a continuous branch of solutions and they do not cover the case of isolated multiple solutions. Only multiple solutions with the same distribution of slip, free and stick zones are considered. These conditions of non-uniqueness require that the friction coecient is a real eigenvalue of a spectral problem. That means that if this spectral problem has a real eigenvalue then the Coulomb friction contact problem is open to non-uniqueness.

In Section 5 we prove the existence of a countable set of complex eigenvalues for the non-local friction model (recalled in Section 4). Moreover, we give there sucient conditions for the existence of at least one real eigenvalue. The eigenvalue problem is approximated in Section 6, and convergence of the nite element method is discussed.

Finally in Section 7, we present some numerical results. First, we implement numerically the eigenvalue problem and we illustrate the convergence of the real eigenvalues. Second we show the non-uniqueness methodology using numerical computations, which unfortunately cannot prove an evidence of non-uniqueness since the convergence results of the nite element model are not established, but which explain quite well the spectral approach proposed in this paper.

PROBLEM STATEMENT

Let an elastic body be given, occupying a domain × R with in R 2 . The generic point in R 3 is denoted x =(x 1 ;x 2 ;x 3 ). We choose plane strain assumptions which means that the displacement eld u =(u 1 ;u 2 ;u 3 ) is vanishing in the Ox 3 direction (u 3 ≡ 0) and u 1 ;u 2 depend only on (x 1 ;x 2 ). The boundary of is assumed to be Lipschitz and is divided as follows: = D ∪ N ∪ C where D ; N and C are three open disjoint parts and meas( D )¿0. We assume that the displacement eld is given on D (i.e. u = U) and that the boundary part N is acted on by a density of surface forces F. The third part is C , which comprises all the points candidate to be in frictional contact with a rigid foundation (see Figure 1). The body is acted upon by a given density of volume forces f . Let n =( n 1 ;n 2 ) represent the unit outward normal vector on and dene the unit tangent vector t =(-n 2 ;n 1 ). We denote by ¿0 the friction coecient on C .

The Coulomb frictional unilateral contact problem consists of nding the displacement eld u :→ R 2 and the stress tensor eld A(u):→ S 2 satisfying (1)-(4):

A(u)=CU(u); div A(u)+f = 0 in (1)
A(u)n = F on N ; u = U on D (2)
where S 2 stands for the space of second order symmetric tensors on R 2 ; U(u)=(∇u +∇ T u)=2 denotes the linearized strain tensor eld, C is a fourth order symmetric and elliptic tensor of linear elasticity and div represents the divergence operator of tensor valued functions.

In order to introduce the equations on C , we adopt the following notation: u = u n n + u t t and A(u)n = n (u)n + t (u)t. The equations modelling contact and friction are as follows on C :

u n 60; n (u)60; n (u)u n = 0 (3)    u t =0⇒| t (u)|6| n (u)| u t =0⇒ t (u)=-| n (u)| u t |u t | (4)
Remark 2.1 Let us mention that the true Coulomb friction law involves the tangential contact velocities and not the tangential displacements. However, a problem analogous to the one discussed here is obtained by time discretization of the quasistatic frictional contact evolution problem. In this case (see Reference [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF]) u; f and F stand for u((i + 1)t); f ((i + 1)t) and F((i + 1)t) respectively and u t has to be replaced by u t ((i + 1)t)u t (it), where t denotes the time step. For simplicity and without any loss of generality only the static case described above will be considered in the following.

The variational formulation of problem ( 1)-( 4) has been obtained by Duvaut and Lions [START_REF] Duvaut | Les in equations en m ecanique et en physique[END_REF]. It consists of nding u verifying:

u ∈ K ad ;a (u; C -u)+j(u; C) -j(u; u)¿L(C -u); ∀C ∈ K ad ( 5
)
where

a(u; C)= (CU(u)): U(C)d;L (C)= f • C d + N F • C d
are dened for any u and C in the standard Sobolev space (H 1 ()) 2 (see Reference [START_REF] Adams | Sobolev Spaces[END_REF]) and the notations • and : stand for the canonical inner products in R 2 and S 2 , respectively. In [START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF], K ad denotes the closed convex set of admissible displacement elds satisfying the non-penetration conditions:

K ad = {C ∈ (H 1 ()) 2 : C = U on D ;v n 60o n C }
The functional j(•; •) given by

j(u; C)=- C n (u)|v t | d (6)
is dened for any C in (H 1 ()) 2 but more regularity is required for u. Two dierent cases when j(u; •) makes sense, are usually considered in the literature. The rst one, which occurs in the continuous problem, involves the space

Ṽ = {C ∈ (H 1 ()) 2 : div A(C) ∈ (L 2 ()) 2 } If u ∈ Ṽ then A(u) belong to H (div; ) and n (u) is an element of H -1=2 ( ) (i.e. the dual of H 1=2 ( )). Since H -1=2 ( )| C is dierent from H -1=2 ( C )
we have to suppose in addition that n (u) ∈ H -1=2 ( C ). With this assumption, (6) makes sense if we replace the integral term by the duality product. For a more precise formulation involving the convenient Sobolev spaces and the set of nonnegative Radon measures, a detailed study can be found in Reference [START_REF]Cimeti ere A. Sur la formulation variationnelle du probl eme de Signorini avec frottement de Coulomb[END_REF]. In the second case, u belongs to a nite element set V h ⊂ (H 1 ()) 2 , which implies that A(u) is at least piecewise continuous so that A(u)n admits a trace on the boundary. In the latter case, the integral notation in [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF] has to be understood in the classical sense.

The rst existence result of ( 1)-( 4) has been proven in Reference [START_REF] Ne Cas | On the solution of the variational inequality to the Signorini problem with small friction[END_REF] when is an innitely long strip and the friction coecient has compact support in C and is suciently small. The extension of these results to domains with smooth boundaries as well as improvements can be found in References [START_REF] Jaru | Contact problems with bounded friction. Coercive case[END_REF] and [START_REF] Kato | Signorini's problem with friction in linear elasticity[END_REF]. More recently in Reference [START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF], existence is stated when the friction coecient is smaller than √ 3 -4=(2-2);denoting Poisson's ratio in (06¡ 1 2 ). To our knowledge there exist neither uniqueness result nor non-uniqueness example of (5) (unless the loads U; f and F are equal to zero).

SUFFICIENT CONDITIONS FOR NON-UNIQUENESS: A SPECTRAL APPROACH

Let us consider an equilibrium position u 0 of the Coulomb frictional contact problem (i.e. a solution of ( 1)-( 4)) supposed to be regular enough. The notation 0 f stands for the points of C which are currently free (separated from the rigid foundation). We denote by 0 s the points of C which are currently in contact but are stuck to the rigid foundation, and by 0 C the points of C which are currently in contact but are candidate to slip. That leads to the following denitions:

0 f = {x ∈ C : u 0 n (x)¡0} (7) 0 s = {x ∈ C : u 0 n (x)=0; | t (u 0 )(x)|¡ -n (u 0 )(x)} (8) 0 C
= {x ∈ C : u 0 n (x)=0; | t (u 0 )(x)| = -n (u 0 )(x)} (9)
Let us adopt the following notation

= ∈ (H 1 ()) 2 such that A()=CU(); div A()=0 in (10) = 0 on 0 D ; A()n = 0 on 0 N ; n = 0 on 0 C ( 11
) t ()= n ()o n 0 C ( 12
)
Remark 3.1 If we choose the commonly used Hooke's law, for homogeneous isotropic materials, given by

ij = E (1 -2)(1 + ) ij kk (u)+ E 1+ ij (u)i n
where E denotes Young's modulus, represents Poisson's ratio and ij is the Kronecker symbol, then the only constitutive constant involved in the eigenvalue problem ( 10)-( 12) is the ratio = =(1 -2). Indeed, the eigenvalues and eigenfunctions are independent of the Young modulus E.

The following theorem states sucient conditions for the non-uniqueness of the equilibrium solution u 0 . Theorem 3.2 Let u 0 be a smooth solution of Coulomb's frictional contact problem ( 1)-( 4) with ¿0 as friction coecient. Let u 1 = u 0 + for some ∈ R and a smooth eigenfunction of ( 10)- [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]. Let us dene the two following cases (i) and (ii):

(i) t (u 0 )(x)60 for all x ∈ 0 C and is the corresponding eigenvalue for . Assume that:

u 1 n (x)¡0; for all x ∈ 0 f ( 13
) n (u 1 )(x)60; for all x ∈ 0 C ( 14
)
| t (u 1 )(x)|¡ -n (u 1 )(x); for all x ∈ 0 s ( 15
)
u 1 t (x)¿0; for all x ∈ 0 C ( 16
)
(ii) t (u 0 )(x)¿0 for all x ∈ 0 C andis the corresponding eigenvalue for . Assume that in addition to ( 13)-( 15), one has:

u 1 t (x)60; for all x ∈ 0 C ( 17
)
If either (i) or (ii) holds then u 1 is another (smooth) solution of ( 1)-( 4).

Proof

Let us rstly remark that u 1 = u 0 + satises Equations ( 1)-( 2) for any ∈ R. Next, we have to check that u 1 veries the frictional contact conditions (3)-( 4). We begin with the unilateral contact conditions (3). If x ∈ 0 f then from [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF] we have u 1 n (x)¡0. Since 0 f ⊂ C , from [START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF] we get A(u 0 )n(x)=0. Having in mind that 0 f ⊂ 0 N and according to [START_REF] Klarbring | On friction problems with normal compliance[END_REF], we deduce n (u 1 )(x)=0. If x ∈ 0 s ∪ 0 C then u 0 n (x)=0 and n (x) = 0, hence u 1 n (x) = 0 and from [START_REF] Haslinger | Approximation of the Signorini problem with friction, obeying the Coulomb law[END_REF] we deduce [START_REF] Jaru | Contact problems with bounded friction. Coercive case[END_REF]. Therefore u 1 satises (3).

If x ∈ 0 s the condition (4) implies u 0 t (x) = 0. From assumption [START_REF] Haslinger | Least square method for solving contact problems with friction obeying Coulomb's law[END_REF] we get

| t (u 1 )(x)| ¡ -n (u 1 )(x) and since 0 s ⊂ 0 D we obtain u 1 t (x)=0. If x ∈ 0 f , then owing to 0 f ⊂ 0 N ,w e have n (u 1 )(x)= t (u 1 )(x) = 0 and (4) is satised. Let x ∈ 0 C .
If case (i) holds, then, for all x ∈ 0 C we have t (u 0 )(x)60 and t (u 0 )(x)= n (u 0 )(x). Using ( 12) with = , we obtain t (u 1 )(x)= n (u 1 )(x) and (4) follows from ( 14) and ( 16). If we consider case (ii) then t (u 0 )(x)=n (u 0 )(x) and from ( 12) with =we nally get t (u 1 )(x)=n (u 1 )(x). Consequently, u 1 satises (4).

In order to apply Theorem 3.2, one has to check the conditions ( 13)-( 17) dealing with the equilibrium u 0 , the eigenfunction and an appropriately chosen value of . The following corollary yields sucient conditions concerning the solution u 0 only. These conditions are more restrictive than those of the previous theorem but easier to handle. Indeed, we shall suppose that all the points of C are in a slipping contact, i.e. 0 s = ∅ and 0 f = ∅.

Corollary 3.3

Let u 0 be a smooth solution of Coulomb's frictional contact problem (1)-(4) with ¿0a s friction coecient. Assume that 0 C = C and that there exist ; ¿0 such that

n (u 0 )(x)6 -for all x ∈ 0 C ( 18
)
Moreover, suppose that one of the following two conditions (i) or (ii) holds: (i) The pair (; ) is a smooth solution of ( 10)-( 12), and

u 0 t (x)¿; for all x ∈ 0 C ( 19
)
(ii) The pair (-; ) is a smooth solution of ( 10)-( 12), and

u 0 t (x)6 -; for all x ∈ 0 C ( 20
)
Then Coulomb's frictional contact problem ( 1)-( 4) admits an innity of solutions. In particular, there exists 0 ¿0 such that u 1 = u 0 + is solution for any satisfying ||6 0 .

Proof Let 0 = sup x∈ C | n ()(x)| and 0 = sup x∈ C | t (x)|. Keeping in mind that n (u 0 + )(x) 6 -+ || 0 we deduce that (14) holds when ||6= 0 .
If condition (i) holds, we can write u 1 t (x)=u 0 t (x)+ t (x)¿ -|| 0 so that the bound ||6= 0 leads to [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF]. Moreover, condition [START_REF] Alart | Crit eres d'injectivit e et de surjectivit e pour certaines applications de R n dans lui même: application a la m ecanique du contact[END_REF] implies that t (u 0 )(x)60 for all x ∈ 0 C .I f we set 0 = min{= 0 ; = 0 } then the rst case in Theorem 3.2 proves the statement of the corollary.

If condition (ii) holds, then [START_REF] Lions | Catastrophic features of Coulomb friction model[END_REF] is satised if ||6= 0 and (20) implies that t (u 0 )(x)¿0 for all x ∈ 0 C . Employing the second case in Theorem 3.2 completes the proof of the corollary.

Remark 3.4

The above results are only sucient conditions for non-uniqueness. They take into consideration only the possibility of existence of multiple solutions having the same distribution of the slip, free and stick zones. Moreover, the above corollary does not cover the case of isolated multiple solutions.

Indeed, as it follows from Corollary 3.3, if the problem is open to non-uniqueness then there exists an innity of solutions located on a continuous branch.

The non-uniqueness conditions considered here imply that the friction coecient (or -) is an eigenvalue of ( 10)- [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]. This eigenvalue problem depends exclusively on the geometry (the domain and the distribution of the dierent types of boundaries) and on the elastic properties incorporated in the operator C (on the Poisson coecient for an isotropic elastic material). As a matter of fact, one can think of a distribution of loads F; f and a displacement eld U such that a solution u 0 of (1)-( 4) for this particular friction coecient satises ( 18)- [START_REF] Alart | Crit eres d'injectivit e et de surjectivit e pour certaines applications de R n dans lui même: application a la m ecanique du contact[END_REF]. We consider, for example, that C is a straight line segment located on the Ox 1axis and that N = ∅. We choose

U(x)= +2 1 - 1 -2 x 2 ; -x 2
for all x =(x 1 ;x 2 ) ∈ D , with ¿0 and f = 0. One can easily check that u 0 (x)=U(x), for all x ∈ is a solution of ( 1)-( 4). Since 0

C = C ; n (u 0 )(x)=-E(1 -)=[(1 -2)(1 + )]¡0 and u 0 t (x)=¿0
, we deduce that the sucient conditions of the corollary hold.

THE NON-LOCAL FRICTION MODEL

There exist several laws 'mollifying' Coulomb's frictional contact model which lead generally to more existence and uniqueness properties. Among these regularization techniques, a special interest is devoted to the non-local procedure introduced in Reference [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF] and developed in References [START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF][START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF]. Moreover, from a physical point of view, this law takes into account some interesting microscopic aspects: the normal pressure n (u) is distributed over a contact area of the deformed asperity (see Reference [START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF] for more arguments). Hence, we consider the non-local normal stress * n (u) given by * n (u

)(x)= C w (|x -y|) n (u)(y)dy C w (|x -y|)dy (21)
where w ; (¿0) stands for a non-negative function with its non-empty support in [-; ] such that x → w (|x|) is an innitely dierentiable function. As for the functional j, the above expression of the non-local normal stress is meaningful in two dierent cases. The rst one concerns the continuous problem when n (u) ∈ H -1=2 ( C ) and the above integral has to be replaced by the duality product between H -1=2 ( C ) and H 1=2 ( C ). The second case happens when using a nite element approximation when (u) is at least piecewise continuous. Another type of smoothing procedure was introduced in Reference [START_REF] Ionescu | Functional and Numerical Problems in Viscoplasticity[END_REF] for friction problems in viscoplasticity. In this case the second order stress tensor eld is averaged in the interior of the domain and its normal trace on the contact boundary provides the non-local normal stress. The denition of the non-local normal stress

* n (u) becomes * n (u)(x)= w (|x -y|)A(u)(y)dy w (|x -y|)dy n(x) • n(x) (22)
Unlike the rst non-local approach in [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF], this second procedure avoids the handling of dual Sobolev spaces such as H -1=2 ( ). Indeed, the latter expression is well dened for any u ∈ (H 1 ()) 2 .

If we replace the above formulas in (4) we get the following 'regularized', non-local friction law on C :

   u t =0 ⇒| t (u)|6| * n (u)| u t =0 ⇒ t (u)=-| * n (u)| u t |u t | (23)
The variational formulation of ( 1)-( 3) and ( 23) is inequality [START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF], the same as in the local friction case in which j is replaced by (see Reference [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF]):

j(u; C)=- C * n (u)|v t |d (24)
From a mathematical point of view, the smoothing map used in the non-local friction model implies compactness properties of the operators involved in the variational approach (5). These properties permit using the Schauder and Tykhonov xed point theorems in order to deduce the existence of at least one solution of the variational inequality [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF][START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF][START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF]. In addition, some uniqueness results can also be obtained for the non-local friction model. As a matter of fact, it was proved in References [START_REF] Duvaut | Equilibre d'un solide elastique avec contact unilat eral et frottement de Coulomb[END_REF][START_REF] Demkovicz | On some existence and uniqueness results in contact problems with nonlocal friction[END_REF][START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF] that there exists a critical friction coecient c such that if ¡ c (i.e. the friction is small) then the solution of ( 5) is unique. As for the local friction case there exist, to our knowledge, no non-uniqueness examples.

The eigenvalue problem corresponding to the non-local friction case is ( 10)-( 11) and

t ()= * n ()o n 0 C ( 25
)
where 0 s and 0 C dened in (8) and ( 9) have to be replaced by

0 s = {x ∈ C : u 0 n (x)=0; | t (u 0 )(x)|¡ - * n (u 0 )(x)} 0 C = {x ∈ C : u 0 n (x)=0; | t (u 0 )(x)| = - * n (u 0 )(x)}
If the normal stress n is replaced by * n then all the sucient conditions for non-uniqueness given in Theorem 3.2 and Corollary 3.3 remain valid.

EXISTENCE OF EIGENVALUES AND EIGENFUNCTIONS

In order to derive the variational formulation of ( 10)-( 11) and ( 25) we consider the subspaces V 0 and Ṽ0 of (H 1 ()) 2 and Ṽ, respectively:

V 0 = {C ∈ (H 1 ()) 2 : C = 0 on 0 D ;v n = 0 on 0 C }; and Ṽ0 = V 0 ∩ Ṽ
Let us introduce the bilinear form b(:; :) given by: b(u; 21) is considered. When adopting the non-local normal stress [START_REF] Adams | Sobolev Spaces[END_REF], this restriction disappears, so that b(u; C) is well dened for all u; C ∈ (H 1 ()) 2 .

The variational formulation of problem ( 10)-( 11) and ( 25) consists of nding ∈ C and 0 = ∈ Ṽ0 such that:

a(; C)=b(; C); ∀C ∈ V 0 ( 26
)
and it can be easily checked that if ∈ C and a non-zero satisfy ( 10)-( 11) and ( 25) then there are a solution of [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF]. Conversely, if ∈ C and a nonzero satisfy [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF], then the pair (; ) is a weak solution of ( 10)-( 11) and ( 25).

Theorem 5.1

The eigenvalues of problem ( 26) consist of a countable set of complex numbers { n } n∈I with n = 0. Each eigenvalue n is of nite algebraic multiplicity. If

I is innite then lim n→∞ | n | = + ∞.

Proof

Let us rst remark that = 0 is not an eigenvalue of [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF]. Otherwise, a(; ) = 0 where is an eigenvector associated with = 0, which contradicts the V 0 -ellipticity of a(:; :). Let us denote by

H = L 2 ( 0 C ) and W 0 = {C ∈ V 0 : div A(C)=0 in ; A(C)n = 0 on 0 N }
and dene P: H → V 0 as follows: for any f ∈ H; P(f) is the unique solution of the variational equality

a(P(f); C)= 0 C fv t d ; ∀C ∈ V 0 :
If we put C ∈ (D()) 2 ⊂ V 0 in the previous equation (the notation D() stands for the space of innitely dierentiable functions with compact support in ), we deduce that div A(P(f)) = 0 in . In the same way, it can be formally checked that A(P(f))n = 0 on 0 N which implies that P(f) ∈ W 0 . Hence, P is a linear continuous operator from H into W 0 . Next, we prove the theorem separately for the two regularization techniques in ( 21) and ( 22).

(i) Case [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF]. The function

* n (C) ∈ H is well dened for any C ∈ W 0 . Set Q:W 0 → H so that Q(C)= * n (C)
, Since Q is a linear and completely continuous operator [9, Theorem 11.2, p. 338] we deduce that T = PQ : W 0 → W 0 is also completely continuous. In order to prove the statement of the theorem, we only have to mention that is a solution in [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF] if and only if 1= is a non-zero eigenvalue for T which is true since T ()= if and only if (; )i s a solution for [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF].

(ii) Case [START_REF] Adams | Sobolev Spaces[END_REF]. The operator Q : V 0 → H , given by Q(C)= * n (C) is well dened. In addition, Q is a linear and completely continuous operator [24, Lemma 1.2, p. 181] and we deduce that T = PQ : V 0 → V 0 is also completely continuous. As a consequence, the proof follows as in case (i).

Remark 5.2

The technique used in the proof above cannot be used if the non-local assumption is removed. The existence of a countable set of eigenvalues is linked to the compactness of operator T which is assured by the regularized trace operator Q.

The following result ensures, under specic conditions, the existence of at least one real positive or negative eigenvalue for problem [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF] which also minimizes the moduli among all eigenvalues satisfying [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF]. First, we need to dene the convex cone K 0 :

K 0 = {C ∈ W 0 : t (C)¿0o n 0 C }
It is easy to see that each displacement led C of W 0 (and of K 0 ) is determined uniquely by the tangential component t (C) of the stress vector on 0 C . Theorem 5.3 Suppose that one of the two following conditions (i) or (ii) holds:

(i) any C in K 0 satises * n (C)¿0o n 0 C , (ii) any C in K 0 satises * n (C)60o n 0 C .
Then the eigenvalue 0 , minimizing the moduli of all eigenvalues in problem [START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF], is real and its associated eigenvector lies in K 0 . Moreover 0 ¿0 in the case (i) and 0 ¡0i nt h e case (ii).

Let us rst recall a weak form of the Krein-Rutman theorem [START_REF] Krein | Linear operator leaving invariant a cone in a Banach space[END_REF][START_REF] Dautray | Analyse math ematique et calcul num erique pour les sciences et techniques[END_REF][START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF] which we use in the proof of Theorem 5.3.

Theorem (Krein and Rutman [START_REF] Krein | Linear operator leaving invariant a cone in a Banach space[END_REF]) Let X be a Banach space and let K ⊂ X be a convex cone containing 0 (i.e. x+y ∈ K; ∀¿0; ∀¿0;x ∈ K; y ∈ K). Suppose that K is closed, X = K -K and K ∩ (-K)={0}. Let T be a linear operator satisfying T (K) ⊂ K.

If T is compact and its spectral radius r(T ) = 0 then there exists ' ∈ K -{0} such that T (')=r(T )'

Proof

Let us consider operator T = PQ : W 0 → W 0 introduced in the proof of Theorem 5.1. For both non-local frictional approaches, the operator T is compact in the Hilbert space W 0 . Moreover, the closed convex cone

K 0 satises K 0 ∩ (-K 0 )=0 and K 0 -K 0 = W 0 (it suces to write t (C)=( t (C)) + -( t (C)
) -where the notations (:) + and (:) -represent the positive and the negative parts, respectively). We next show that T (K 0 ) ⊂ K 0 . The assumptions of the theorem imply that the operator Q dened by

Q(C)= * n (C) maps K 0 into (L 2 ( 0 C )) + . The operator P dened for all f ∈ L 2 ( 0 C )b y a(P(f); C)= 0 C fv t d ; ∀C ∈ W 0 satises t (P(f)) = f. Hence C ∈ K 0 , which implies that Q(C) ∈ (L 2 ( 0 C )) + and thus T (C) ∈ K 0 .
It follows then from Krein-Rutman's theorem that if T admits a positive spectral radius, then there exists an eigenvalue which is equal to the spectral radius with an associated eigenvector in K 0 .

The case (ii) is handled similarly by using the operator -T .

FINITE ELEMENT APPROXIMATION OF THE EIGENVALUE PROBLEM

The problem we intend to approximate is as follows: nd ∈ C and 0 = ∈ V 0 such that:

a(; C)=b(; C); ∀C ∈ V 0 (27)
which is exactly the eigenvalue problem corresponding to the non-local frictional approach [START_REF] Adams | Sobolev Spaces[END_REF]. Notice that when the regularization procedure ( 21) is adopted, then the convergence analysis is more complicated. A remark at the end of this section explains and gives partial answers to the convergence study in that case. We denote by • 1 the standard norm on (H 1 ()) 2 . Our aim is to approximate the eigenvalues of problem [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF]. Let be given a family of nite dimensional subspaces V 0 h ⊂ V 0 where h denotes the discretization parameter [START_REF] Ciarlet | The nite element method for elliptic problems[END_REF]. The nite dimensional problem consists then of nding h ∈ C and 0 = h ∈ V 0 h such that [START_REF] Fix | Eigenvalue approximation by the nite element method[END_REF][START_REF] Babu Ska | Eigenvalue problems[END_REF]:

a( h ; C h )= h b( h ; C h ); ∀C h ∈ V 0 h (28)
We assume that the following approximation property holds:

lim h→0 inf u h ∈V 0 h u -u h 1 =0; ∀u ∈ V 0
Let -1 be an eigenvalue of T dened in (ii) of the proof in Theorem 5.1. Denoting by I the identity map, there exists a least integer such that Ker(( -1 I -T ) ) = Ker(( -1 I -T ) +1 )=E with dim(E)=m¡∞. The algebraic multiplicity of -1 is m and stands for the ascent of -1 I -T . The set E contains the generalized eigenvectors of T corresponding to -1 . Let T * be the adjoint operator of T dened on the dual space V 0 * . Then -1 is an eigenvalue of T * with algebraic multiplicity m and is also the ascent of -1 I -T * . The notation E * = Ker(( -1 I -T * ) ) stands for the space of generalized eigenvectors of T * associated with -1 . Given two closed subspaces A and B of V 0 , dene the gap between A and B by

(A; B) = max sup u∈A; u 1 =1 inf C∈B u -C 1 ; sup u∈B; u 1=1 inf C∈A u -C 1
Let be an eigenvalue of [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF] and denote by m its algebraic multiplicity. When h tends to zero, there exist exactly m eigenvalues of (28) denoted 1;h ; 2;h ;:::; m; h converging to . Let E h be the direct sum of the generalized eigenspaces associated with 1;h ; 2;h ;:::; m; h and dene

h = sup u∈E; u 1=1 inf C h ∈V 0 h u -C h 1 and * h = sup u∈E * ; u 1=1 inf C h ∈V 0 h u -C h 1
The following theorem, taken from Kolata [START_REF] Kolata | Approximation in variationally posed eigenvalue problems[END_REF], describes the convergence of the nite element approximation.

Theorem 6.1

If h is small enough, the following estimates hold:

- 1 m m i=1 i; h 6 C h * h | -i; h | 6 C( h * h ) 1 ; 16i6m
(E; E h ) 6 C h
where the constant C does not depend on h.

If the rst regularizing approach ( 21) is adopted then the eigenvalue problem becomes: nd ∈ C and 0 = ∈ W 0 such that a(; C)=b(; C), ∀C ∈ V 0 . In such a case there are at least two alternatives for obtaining convergence results. The rst one is to invoke again Kolata's studies in Reference [START_REF] Kolata | Approximation in variationally posed eigenvalue problems[END_REF] which are still valid. It suces then to show two families of nite dimensional subspaces V 0 h and W 0 h of V 0 and W 0 , respectively, where the dimensions of V 0 h and W 0 h are equal. This necessitates to introduce more specic nite element spaces which is out of the scope of this paper. The second possibility is to use a non-conforming nite element approach and approximating V 0 and W 0 with the same nite dimensional space V 0 h although V 0 h ⊂ W 0 . In that case the convergence result requires strong supplementary hypotheses as in Reference [START_REF] Mercier | Eigenvalue approximation by mixed and hybrid methods[END_REF].

NUMERICAL RESULTS

The sucient conditions for non-uniqueness given in Section 3 concern the solution u 0 of the continuous problem (1)-(4). Unfortunately, there are to our knowledge no available analytical examples of non-uniqueness in the continuous framework. That is why we cannot directly verify the sucient conditions given in the Theorem 3.2 or Corollary 3.3. Next we try to illustrate the sucient conditions from a numerical point of view. Since the convenient convergence results for the Coulomb friction model do not exist, the numerical computations cannot stand for a rigorous mathematical proof of the sucient conditions for non-uniqueness. Our aim here is only to illustrate the methodology given in the continuous context.

We consider the eigenvalue problem ( 10)-( 12) where represents a unit square whose partition of the boundary ( 0 C ; 0 D ; 0 N ) is depicted Figure 2. The elastic material is supposed to be isotropic and a Poisson ratio equal to 0.4 is chosen. The nite element discretization is made of uniform quadrilateral meshes of edge size h. The ARPACK library is used to compute the eigenvalues and eigenfunctions of the discrete spectral problem [START_REF] Ciarlet | The nite element method for elliptic problems[END_REF].

Figure 3 depicts the two smallest real eigenvalues denoted 0;h and 1;h as a function of the number of elements on 0 C ( 0;h and 1;h are single eigenvalues). As expected from Theorem 6.1, we note good convergence of these eigenvalues with respect to the mesh size. Note that 0;h , which admits the smallest modulus among all complex eigenvalues converges quite well.

In Figure 4, we show the eigenfunction 0;h , corresponding to 0;h , computed for a mesh of 100 × 100 elements. Notice that all the nodes on 0 C are slipping. Figure 5 represents the normal stress n ( 0;h ) and the tangential slip displacement 0;ht on 0 C for three dierent sizes of meshes. We note the good agreement between these three computations, except maybe, in the neighborhood of the lower right corner of .

The existence of a real eigenvalue for (10)-( 12) is investigated in Section 5. As it was already pointed out in Remark 3.1, the eigenvalues depend only on the ratio = =(1 -2) in the case of a homogeneous and isotropic elastic body. In Figure 6, we observe that for all the rst eigenvalue is real and positive when the computations are performed on a 50 × 50 mesh. Note the sharp variation for small values of the ratio . If for some the positivity of the rst eigenvalue is preserved when h tends to 0 and it converges to an eigenvalue of the continuous problem, then the geometry is open to non-uniqueness according to Remark 3.5.

In order to discuss our non-uniqueness methodology, we choose the same geometry as in the previous test and a Poisson ratio of 0.4. So, we consider an elastic unit square lying on a rigid inclined plane (see Figure 7 numerical results plotted in Figure 8 are obtained using an uniform quadrilateral mesh of 100 × 100 elements.

The friction coecient is chosen to be equal to 0;h =0:84232, the rst eigenvalue computed above for the same mesh. We remark that all the points of C are slipping and therefore we have 0 C = C , 0 D = D , and 0 N = N . Figure 9 represents the normal stress and the tangential displacement on C for three dierent mesh sizes. We remark that the sucient conditions ( 18)-( 19) are satised for the three computations. From the computations we can reasonably expect that these conditions still hold as h tends to 0. If in addition we assume convergence of the nite element solution u 0 h to a solution u 0 of the continuous model ( 1)-( 4) such that L ∞ convergence of the quantities depicted in Figure 9 holds, then the conditions ( 18)-( 19) also hold for the solution u 0 .I f the rst eigenvalue 0;h in Figure 3 converges to a positive eigenvalue 0 of ( 10)-( 12), we choose = 0 and from Corollary 3.3 we could deduce that u 0 is not the unique solution of (1)-(4). Indeed, there exist an appropriate scalar such that u 1 = u 0 + 0 is another solution of (1)-(4). In Figure 10 we try to illustrate this with a very rened mesh (as ne as the computation allows). We have plotted the normal stress and the tangential displacement corresponding to the discrete solutions u 1 h and u 0 h computed with a 200 × 200 mesh. Although the solutions are not quite dierent, we nevertheless observe a local mist (gap) of order of 10%. The Von-Mises stress corresponding to the dierence u 1 hu 0 h is depicted in Figure 11. In this case the dierence (of order of 20%) is concentrated on the lower right corner.

Figure 1 .

 1 Figure 1. Setting of the problem. The domain ; its boundary is divided into three parts: D ; N and C .

  the following eigenvalue problem: Eigenvalue problem. Find ∈ C and 0

Remark 3. 5

 5 If (10)-(12) admits a real eigenvalue then the pair geometry-material is open to the nonuniqueness of the Coulomb frictional contact problem.

Figure 2 .

 2 Figure 2. The geometry, the domain and its boundary with its three open disjoint parts: D , N and C .

Figure 3 .

 3 Figure 3. The two smallest real eigenvalues 0;h and 1;h vs the number of elements on a edge: 1=h.

Figure 4 .

 4 Figure 4. The eigenfunction 0;h computed for a mesh of 100 × 100 elements; all the nodes in 0 C are slipping.

Figure 5 .

 5 Figure 5. (a) The normal stress n( 0;h ) and (b) the tangential slip 0;ht on 0 C for three dierent meshes: 50 × 50, 100 × 100 and 200 × 200 elements.

Figure 6 .

 6 Figure 6. The dependence of the rst eigenvalue on the ratio = =(1 -2).

Figure 7 .

 7 Figure 7. An example of mechanical problem (1) -(4) with non-uniqueness. An elastic square lying on a rigid inclined plane, loaded by the gravity f and by an imposed displacement U.

Figure 8 .

 8 Figure 8. The deformed conguration and the Von-Mises stress eld associated with the solution u 0 h , computed with a 100 × 100 mesh.

  ) as an example of the mechanical problem (1)-(4). The elastic body is loaded by a density of gravity forces f and by an imposed displacement eld U. The following values are used: = Arctan 2 ≃ 63:43 • , = Arctan 2:5 ≃ 68:2 • , |f | = 223:6Nm -3 , |U| =0:005385, and a Young modulus E = 10 GPa. The computations of the numerical solution denoted u 0 h are achieved using the nite element code CASTEM. The

Figure 9 .

 9 Figure 9. The normal stress and the tangential displacement of C , corresponding to the solution u 0 h , computed with three dierent meshes.

Figure 10 .

 10 Figure 10. The normal stress and the tangential displacement corresponding to the solutions u 1 h and u 0 h , computed with a 200 × 200 mesh.

Figure 11 .

 11 Figure 11. The Von-Mises stress corresponding to the dierence u 1 h -u 0 h .

  C in (H 1 ())2 . Concerning the rst variable, b(u; C) makes sense if the non-local normal stress * n (u) can be dened. Hence, b(u; C) is well dened for any u ∈ Ṽ such that n (u) ∈ H -1=2 ( C )o rf o ru ∈ V h (the notation V h represents a nite element type space) if the non-local normal stress (

C)= C * n (u)v t d

for any


Benjamin Adroit
email: badroit@uni-bonn.de

Torsten Wappler

Jean-Frederic Terral

Adam A Ali

Vincent Girard

Keywords: paleoecology, paleoforest structure, Pleistocene, fossil leaves, insect damage

Plants and insects have interacted for millions of years in terrestrial ecosystems and are essential to the function of the food web. These interactions can be classified based on different functional feeding groups and damage types, and provide information about the relationship between plants and insects. We focused this study on the Pleistocene deposit of Bernasso (Hérault, southern France). The Bernasso fossil leaf assemblage has been dated from 2.16-1.96 Ma and corresponds to a deciduous forest comprising both continental and Mediterranean species such as Caprinus orientalis Mill., Parrotia persica C. A. Mey. and Acer monspessulanum L., Sorbus domestica L., respectively. Twenty different plant species are represented by 590 fossil leaves among which more than 36% present insect damage. Thirty-nine different damage types were found with galling being the most prominent. These types of interactions enable discussion about the ecological settings that prevailed at this epoch (forest age, precipitation, temperature). For instance, the mean annual precipitation was certainly higher than today in the region. Moreover, the results of this study offer the possibility to get new insights into the ecological paleoforest community.

Introduction

There is no doubt that human activities in recent decades have impacted the climate [START_REF] Ellis | Long-term effects of a trophic cascade in a large lake ecosystem[END_REF]. Nowadays, some of these anthropogenic impacts have reached a tipping point, generating a "global change" [START_REF] Turner | The Earth as Transformed by Human Action: Global and Regional Changes in the Biosphere over the Past 300 Years[END_REF][START_REF] Heller | Biodiversity management in the face of climate change: A review of 22 years of recommendations[END_REF]. The study of past climate thus can be very helpful to evaluate the impacts of the global change on multitrophic interactions and biodiversity [START_REF] Dowsett | Pliocene three-dimensional global ocean temperature reconstruction[END_REF][START_REF] Salzmann | Climate and environment of a Pliocene warm world[END_REF]. Regionally, for the northern Mediterranean realm and western Europe the highly diverse forest vegetation implies MATs 3-6°C higher than today with a respectively increased in mean annual precipitation (MAP) by 400 mm and 230 mm respectively [START_REF] Utescher | Terrestrial climate evolution in northwest Germany over the last 25 million years[END_REF][START_REF] Contoux | Modelling the mid-Pliocene Warm Period climate with the IPSL coupled model and its atmospheric component LMDZ5A[END_REF]. The current climate change will probably upset ecological processes, especially between vegetation and insects.

Interactions between plants and insects have played an important role in the structure of the food web and the study of their evolution could help to reconstruct climate dynamics [START_REF] May | How many species are there on Earth?[END_REF][START_REF] Wilf | Insect-damaged fossil leaves record food web response to ancient climate change and extinction[END_REF]. Plant-insect interactions have also been shown to virulently impact entire modern ecosystems [START_REF] Cornelissen | Climate change and its effects on terrestrial insects and herbivory patterns[END_REF]. Distinctive signs found on fossil leaves are thus important to understand the evolution of interactions [START_REF] Labandeira | Insect fluid-feeding on Upper Pennsylvanian tree ferns (Palaeodictyoptera, Marattiales) and the early history of the piercing-and-sucking functional feeding group[END_REF][START_REF] Wilf | Insect herbivory, plant defense, and early Cenozoic climate change[END_REF][START_REF] Wilf | Decoupled Plant and Insect Diversity After the End-Cretaceous Extinction[END_REF][START_REF] Currano | Sharply increased insect herbivory during the Paleocene-Eocene Thermal Maximum[END_REF][START_REF] Wappler | No post-Cretaceous ecosystem depression in European forests? Rich insect-feeding damage on diverse middle Palaeocene plants, Menat, France[END_REF][START_REF] Wappler | Testing for the Effects and Consequences of Mid Paleogene Climate Change on Insect Herbivory[END_REF][START_REF] Labandeira | A paleobiologic perspective on plant-insect interactions[END_REF]. Therefore, the fossil record offers a unique possibility for conducting meta-analysis and for analyzing the dependence between climate parameters and plant-insect interactions through geologic time [START_REF] Wilf | Insect herbivory, plant defense, and early Cenozoic climate change[END_REF][START_REF] Currano | Fossil insect folivory tracks paleotemperature for six million years[END_REF][START_REF] Wappler | Insect herbivory close to the Oligocene-Miocene transition -A quantitative analysis[END_REF][START_REF] Knor | Plant-arthropod associations from the Early Miocene of the Most Basin in North Bohemia-Palaeoecological and palaeoclimatological implications[END_REF][START_REF] Wappler | Testing for the Effects and Consequences of Mid Paleogene Climate Change on Insect Herbivory[END_REF][START_REF] Labandeira | The Fossil Record of Plant-Insect Dynamics[END_REF]. Bernasso, a locality near Lunas (Hérault, southern France), provides a new snapshot of early Pleistocene multitrophic interactions between insects and plants.

Previous studies on pollen records [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF]) and fossil leaf assemblages [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF] suggest that Bernasso's climate was warmer and wetter than today. However, a recent study based on CLAMP analysis, still not published yet, could suggest a climate certainly wetter but cooler than today (unpublished data).

Moreover, Bernasso provides a rare glimpse into multitrophic interactions between insects and plants in a time period slightly wetter and under different temperatures, as well as a similar geographic configuration to today [START_REF] Rohling | Antarctic temperature and global sea level closely coupled over the past five glacial cycles[END_REF][START_REF] Raymo | Departures from eustasy in Pliocene sea-level records[END_REF]. [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF] showed that Bernasso fossil leaves fossilized in a lacustrine deposit, represented a deciduous forest mainly dominated by Parrotia persica C. A. Mey., Carpinus orientalis Mill and other important species as e.g.

Zelkova ungerii

Spach., Acer Monspessulanum L., Fraxinus ornus L., Sorbus domestica L. and Carya minor Nutt. Today, there is no exact modern analog for the Bernasso paleoforest anywhere in the world. Bernasso was a combination between typical Mediterranean plants species (A. monspessulanum, F. ornus…) and some Caucasian tree species such as P. persica (its last natural population is located around the Caspian Sea in Iran; [START_REF] Akhani | Plant biodiversity of Hyrcanian relict forests, N Iran: an overview of the flora, vegetation, palaeoecology and conservation[END_REF].

The Bernasso deposits were formed 2.16-1.96 Ma ago, just after the beginning of glacial/interglacial cycles (-2.6Ma). Pollen assemblages show that during the time lapse of deposition, two glaciations occurred separated by an interglacial event [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF][START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF]. Presently, studies show the possible impact of the glacial/interglacial cycles on precise species [START_REF] Dechaine | Using coalescent simulations to test the impact of quaternary climate cycles on divergence in an alpine plant-insect association[END_REF][START_REF] Delucia | Climate Change: Resetting Plant-Insect Interactions[END_REF]. For example, the genetic diversification of Myrtus communis in the Mediterranean region was possible because of a glaciation event around 2 Ma ago [START_REF] Médail | La biodiversité végétale méditerranéenne, organization et evolution[END_REF]. In parallel, other studies suggest that the variability in abundance and diversity of plant-insect interactions could be due to temperature fluctuations [START_REF] Currano | Fossil insect folivory tracks paleotemperature for six million years[END_REF], plant defense strategies [START_REF] Coley | Herbivory and plant defenses in tropical forests[END_REF], and age of the forest [START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF]. Furthermore, [START_REF] Nicole | Les relations des insectes phytophages avec leurs plantes hôtes[END_REF] has shown that there is also a strong correlation between some plant species and interactions amongst them.

In the present study, interaction structure of the paleoforest of Bernasso was investigated. This study corresponds to one of the first that focus on the herbivory structure of a forest ecosystem at the beginning of glaciation cycles. We tried to highlight, with the example of Bernasso, if interaction structures were already in agreement with climatic conditions or if forest ecosystems were still slowly acclimatized to the recurrent climatic fluctuation of that period.

Material & Methods

Study area

The lacustrine Bernasso deposits (43°43'34"N, 3°15' 25"E) (Figure 1) are located 5km to the east of Lunas (Hérault, southern France) at an altitude of 512m asl [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF][START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF]. Bernasso is characteristic of the supramediterranean bioclimatic zone. The MAT is around 12.05°C and the MAP is 1070mm [START_REF] Kessler | Météo de la France. Tous les climats localité par localité[END_REF]. This region is now mostly dominated by deciduous taxa with maples (Acer monspessulanum, A. campestre L.), oaks (Quercus pubescens Willd) and ashes (Fraxinus oxyphylla L., F. excelsior L.). Bernasso lacustrine deposits, composed of diatomites, were formed 2 Ma ago.

They are rich in plant remains (pollen and fossil leaves) [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF][START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF]. Basaltic flows (due to the activity of the Escandorgue volcano) shut a canyon valley and allowed the development of a lake in which diatom floras developed. Lacustrine deposits correspond to the early Pleistocene [START_REF] Hilgen | Astronomical calibration of Gauss to Matuyama sapropels in the Mediterranean and implication for the geomagnetic polarity timescale[END_REF].

Dating was done by the K/Ar analysis [START_REF] Ildefonse | Mise en évidence de la transition paléomagnétique Gauss-Matuyama dans les formations volcaniques de l'Escandorgue, Hérault, France[END_REF] on the diatomite and completed with cyclostratigraphy (Suc and Popescu, 2005) and paleomagnetism work [START_REF] Ambert | Révision chronostratigraphique de la séquence paléobotanique de Bernasso (Escandorgue, Midi de la France)[END_REF]. For more details, see [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF].

The Bernasso palaeoflora was mainly dominated by Carpinus orientalis, Parrotia persica, Zelkova ungeri and Carya minor and indicate that modern Caucasian flora could serve as models for such vegetation [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF].

Mediterranean species such as Acer monspessulanum, Sorbus domestica, Tilia sp. and Vitis sp. are examples of the macrofossils [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF], Buxus L., Olea L., Ericaceae and Poaceae [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF] are also preserved in the microfossils.

Pollen analyses have shown a succession of three different climatic periods [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF]. The base and top of the sequence correspond to cold periods (glaciations) and the regional vegetation was steppic. The intermediate period was warmer (interglacial) and allows the development of woody vegetation. The leaf fossils came from this last climatic period.

Beside plant micro-and macrofossils, Bernasso diatomites provided few others fossils. For example, only two fossil insects have been found at Bernasso (a coleoptera and an orthoptera) (pers. obs. B.Adroit, 2013).

Fossil leaves collection

The collection includes more than 800 fossil leaves, collected between 1995 and 2012. Only 590 were used for this study exclusively because the complete or almost completely preserved specimens were chosen. All the specimens are conserved in the leaf paleontological collection of Institut des Sciences de l'Evolution de Montpellier, France (ISEM). The collection is split into two sets according to the kind of preservation of the fossil. The first group contains 321 specimens (54.4% of all) preserved on the extracted rocks (hereafter Ro). The second group contains 269 specimens (45.6% of all) totally removed from the diatomite with the peeling method (hereafter BDC). With this method, the dried fossil leaf is fixated on a plastic film and placed between two glass slides with Canada balsam (fir resin). Only terrestrial angiosperms were taken into account in this study. The aquatic species Ceratophyllum L., ferns and conifers were excluded. For the complete flora list and illustrations, see [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF].

Fossil analysis and classification of damages

The fossils were observed under a Leica MZ7.5 equipped with a Leica DFC300Fx camera. The BDC fossils were investigated under an Olympus MX51 microscope. All of the plant-insect interactions were identified using the Guide to Insect Damage Types (DT) by [START_REF] Labandeira | Guide to insect (and other) damage types on compressed plant fossils[END_REF]. Functional feeding group (FFG) and host specificity (HS) defined by [START_REF] Labandeira | Guide to insect (and other) damage types on compressed plant fossils[END_REF] were determined for each damage type indicated, host specificity is an index assigned for every DT between 1 and 3 (the higher the HS is, the more specialized the damage is). Two kinds of interactions are recognized: specialized damages that include galling, mining, piercing & sucking, and oviposition (damages with mainly HS=2/3) and external damages that comprise holes feeding, margin feeding, surface feeding, and skeletonization (damages with mainly HS=1). All data were compiled in a spreadsheet including leaf species, leaf quantity, DT, FFG (Table S1).

Statistic analysis

The R software (R Core Team, 2012) and its vegan package [START_REF] Oksanen | Package "vegan[END_REF] were used to perform rarefaction test in order to analyze diversity of FFG according to leaf quantity. A rarefaction test was also performed on the proportion of specialized damage types in regard to all damages. As the data does not follow a normal distribution, a Wilcoxon test was used when necessary to check if there was a significant difference between two species. At last, a Chi-squared test (α < 0.1%) was used to compare the proportion of DTs to the kind of fossils (Ro or BDC). The goal was to determine if the type of fossil preservation influenced conservation and/or detection of DTs. The statistical analyses were performed on leaf species represented by at least 20 specimens. Seven species were analyzed by statistics.

Results

3.1 Damage diversity and distribution on the leaf species 37% of the Bernasso leaves (220/590) displayed interactions with insects. 39 different DTs were identified. These DTs are organized into external foliage feeding (with seven types of holes feeding, eight types of margin feeding, eight types of skeletonization, three types of surface feeding and one undetermined DT) and specialized foliage feeding (with four types of mining, seven types of galling and a single type of oviposition) (Figure 2). Oviposition is represented by a unique occurrence detected on a leaf of P. persica. Figure 3 illustrates the distribution of the seven FFGs within the different leaf species. P. persica and C. orientalis (the most abundant plant species) have the most important quantities of damages. P. persica has a noticeable dominance of galling while C. orientalis has a dominance of holes feeding and margin feeding. Galling represents half of the specialized damages (p < 2.2e-16) (Figure 4) and its diversity is seven DTs. Moreover, specialized damages (21 DTs), in which galling is included, are significantly more diversified than external damage (18DTs) (p = 2.534e-12).

At Bernasso four FFGs represent 90% of the total DTs: galling being 30% of the total interaction observed, hole feeding 24%, margin feeding 18% and skeletonization 18%. Galling is statistically more present than hole feeding (Wilcoxtest: p < 2.2e-16), margin feeding (Wilcox-test: p = 1.56e-12) and skeletonization (Wilcox-test: p < 2.2e-16). Hole feeding and margin feeding are statistically equally distributed (Wilcox-test: p = 0.1881) and skeletonization is less representative than hole feeding among the Bernasso fossil leaf collection (Wilcox-test: p = 3.116e-05).

No statistically significant difference between margin feeding and skeletonization is observable (Wilcox-test: p.value= 0.078). Concerning the distribution of these four FFGs over the palaeoflora (Figure 3), galling is represented in 19 different species but is abundant only in four leaf species: P. persica (34% of the galling), A. monsepessulanum (10%), S. domestica (10%) and F. ornus (9%). Galling is more important on P. persica than A. monsepessulanum (p = 3.472e -05 ), S. domestica (p = 3.472e -05 ) and F. ornus (p = 0.01231). The most important galling DT is DT32 (60% of all galling). It is also the main galling found on P. persica and A. monspessulanum (these 2 species contain 64% of DT32). Another important galling type is DT80 that represents 23% of all galling. Among this damage type DT80, 53% of its occurrences are present on S. domestica and F. ornus plant species. Skeletonization is present on 15 different species and 45% of this FFG is represented by DT16. This FFG is widely distributed but leaves of Acer have the largest proportion (36%) of total skeletonization. Lastly, hole feeding and margin feeding are respectively found on 11 and 14 species. These two FFGs are homogeneously present on all leaf species. Hole feeding is in correspondence to two DTs DT01 (25%) and DT02 (29%), whereas 63% of margin feeding is represented by DT12 only.

Conservation bias

First of all, there are statistically more fossil leaves within the Ro (54%) than the BDC (46%) dataset (Chi-squared test: p = 0.03; α = 1%; ddl = 1).

BDC fossil leaves contain 70% of the total DTs detected while Ro fossils contain only 30% of the DTs (Figure 5). This difference is significant (Chi-squared test: p < 0.01; α = 0.1%; ddl = 1) and is mostly due to two FFGs: margin feeding and galling (Figure 5). These two FFGs are at least twice more abundant within BDC fossils than Ro fossils (Chi-squared test: p < 0.01; α= 0.1%; ddl = 1). Concerning the two other important FFGs (hole feeding and skeletonization), the difference between BDC and Ro fossil is not significant (Chi-squared test: p = 0.27 and p = 0.26; α = 0.1%; ddl = 1).

Discussion

Representative and taphonomic bias

Rarefaction data shows that all the species exhibit damage. The sampling effort during the collecting field trips was sufficient for representing plant diversity from 300 fossil leaves (19.21 ± 0.74) Also the significant difference between the kinds of fossil preparation (BDC and Ro) could result from biases generated by the BDC preparation. This particular method could lead to an overestimation of margin feeding. During the transfer, parts of leaves could be torn. Moreover, the Canada balsam tends to blacken the margin of the leaves, mimicking the reaction after an insect feeding attack. On the contrary, an under-estimation of galling could occur for the Ro fossils. Most galls were more readily observed on BDC fossils, possibly due to the transparency of these fossils.

Damage abundance on the paleoforest

The two most abundant species P. persica and C. orientalis have had 40% of their leaves attacked. However, FFGs on these two dominant species are completely different. C. orientalis mostly contains hole feeding and margin feeding that are general types of damage caused by polyphagous insects [START_REF] Currano | Sharply increased insect herbivory during the Paleocene-Eocene Thermal Maximum[END_REF][START_REF] Wilf | Insect-damaged fossil leaves record food web response to ancient climate change and extinction[END_REF][START_REF] Wappler | Insect herbivory close to the Oligocene-Miocene transition -A quantitative analysis[END_REF]. P. persica mainly have galling, a specialized interaction caused by oligophagous insects [START_REF] Currano | Fossil insect folivory tracks paleotemperature for six million years[END_REF]. P. persica also has a lot of generalist interactions (e.g., DT01, DT02, DT12), probably produced by polyphagous insects [START_REF] Currano | Sharply increased insect herbivory during the Paleocene-Eocene Thermal Maximum[END_REF]. These results are not exactly congruent with the hypothesis of [START_REF] Currano | Fossil insect folivory tracks paleotemperature for six million years[END_REF] that in modern ecosystems insects feed more on the dominant plant species [START_REF] Novotny | Why Are There So Many Species of Herbivorous Insects in Tropical Rainforests?[END_REF]. Proportionally P. persica and C. orientalis are not the two most damaged species. The size of these two species can provide some explanation. Modern P. persica and C. orientalis in Iran can grow to the height of 15-20m (Abedi M., pers. comm. 2015). If it is estimated that sizes of these two species were the same at Bernasso 2 Ma ago, P. persica and C. orientalis should have been the dominate trees of the canopy in the Bernasso paleoforest and their leaves were subject to more intense sunlight than other species. This probably resulted in a more efficient photosynthesis for P. persica and C. orientalis and therefore an increase of C/N ratio in their leaves [START_REF] Bazzaz | The Response of natural ecosytems to the rising global CO2 levels[END_REF][START_REF] Barradas | Responses of CO2 assimilation to changes in irradiance: laboratory and field data and a model for beans (Phaseolus vulgaris L.)[END_REF]. These two phenomena are both unfavorable for herbivorous organisms [START_REF] Yamasaki | Temporal and spatial variations in leaf herbivory within a canopy of Fagus crenata[END_REF]. Furthermore, P. persica and C. orientalis have both specialist and generalist interactions. According to [START_REF] Currano | Sharply increased insect herbivory during the Paleocene-Eocene Thermal Maximum[END_REF], the physiology and/or anatomy of these two plants should be directly involved in their diversity of interactions.

Proportionally A. monspessulanum and S. domestica are more attacked than P. persica and C. orientalis. 59% of A. monspessulanum leaves and 54% of S. domestica leaves possess damages. These proportions have to be compared to the number of leaves of each species. C. orientalis and P. persica are represented by more than 100 specimens, A. monspessulanum and S. domestica are only represented by ca. 20 specimens. It is necessary to have more A. monspessulanum and S. domestica specimens to statistically test these differences. However, it seems that the majority of insects adapted to feed on the two dominant species (C. orientalis and P. persica) and more specialized insects adapted for feeding on the other less represented plant species (A. monspessulanum and S. domestica).

Community structure of the paleoforest: ecological parameters

Specialized damages, especially galling, are an interesting tool to discuss a paleoforest community structure [START_REF] Fernandes | The adaptive significance of insect gall distribution: survivorship of species in xeric and mesic habitats[END_REF][START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF]. Figure 6 compares the quantity and the diversity of galling at Bernasso according to the theories of [START_REF] Fernandes | The adaptive significance of insect gall distribution: survivorship of species in xeric and mesic habitats[END_REF][START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF]. Three ecological parameters were analyzed: the age of Bernasso forest, the temperature and the precipitation at Bernasso ca. 2 Ma ago. The aim of figure 6 is to synthetize the following interpretation.

At Bernasso, specialized damages, mainly due to galling, are important but undiversified (Figure 6). DTs 32 and 80 represent nearly 90% of the total galling diversity and they were found on almost all species (Table S1). Moreover their HS is low [START_REF] Labandeira | Guide to insect (and other) damage types on compressed plant fossils[END_REF] even if galling, in general, represents a specialized FFG (e.g., [START_REF] Knor | Galls and gall makers on plant leaves from the lower Miocene (Burdigalian) of the Czech Republic: Systematic and palaeoecological implications[END_REF]. This low diversity connected to the large quantity of galling is in contrast with [START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF] (Figure 6). Indeed, when compared to the modern forests in the Hyrcanian region (Iran) [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF], Bernasso could be considered as a late-successional forest. P. persica, Z. carpinifolia and C. orientalis are known in Iran to be the main species of late-successional forests [START_REF] Djamali | Relations pluie polliniquevégétation sur un transect forêt-steppe dans le parc national du Golestan, nord-est de l'Iran[END_REF]. According to [START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF], in old forest such as Bernasso, there should be an expected important quantity and diversity of galling. In the case of Bernasso, the large quantity of galling may be explained by the age of the forest (Figure 6). Also, [START_REF] Fernandes | The adaptive significance of insect gall distribution: survivorship of species in xeric and mesic habitats[END_REF] suggest that a cold temperature implies a low diversity and quantity of galling. Thus, the low diversity of galling at Bernasso could suggest low mean annual temperature (Figure 6). Lastly, precipitation may also be correlated with diversity and quantity of galling. A high ratio of precipitation seems to imply a low quantity and diversity of galling [START_REF] Fernandes | The adaptive significance of insect gall distribution: survivorship of species in xeric and mesic habitats[END_REF]. Thus, the high mean annual precipitation [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF][START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF] could also explain the low diversity of galling at Bernasso. All these information and data may support the results of the unpublished CLAMP data.

More generally, Bernasso fossils (590 specimens/20 species) have a particularly high diversity of different DTs (39 DTs), which is comparable to other fossil localities. The Oligocene locality of Quegstein (404 specimens/36 species) has 4.5 ±1.6 DT for 50 leaves when Bernasso have 12.23 ± 2.38 DT for 50 leaves [START_REF] Wappler | Insect herbivory close to the Oligocene-Miocene transition -A quantitative analysis[END_REF]. The locality of Rott, another Oligocene site (2476 specimens/112 species) only has 7.5 ± 2.1 DT for 50 leaves [START_REF] Wappler | Insect herbivory close to the Oligocene-Miocene transition -A quantitative analysis[END_REF]. Even if it is difficult to quantify insect diversity (because it is not always correlated to the diversity of damage; [START_REF] Currano | Sharply increased insect herbivory during the Paleocene-Eocene Thermal Maximum[END_REF], the large quantity of DTs at Bernasso could suggest an increased diversity of insects. However, 85% (34 DTs) of damage types are rare and only 15% of them (6DTs: DT01,02,12,16,32,80) are dominant. [START_REF] Carvalho | Insect Leaf-Chewing Damage Tracks Herbivore Richness in Modern and Ancient Forests[END_REF] argued that the abundance of rare DTs could mainly be the result of the activity of multi-damaging insects. This may indicate that the Bernasso DT diversity could be biased by the presence and diversity of multi-damaging insects [START_REF] Novotny | Low host specificity of herbivorous insects in a tropical forest[END_REF][START_REF] Currano | Fossil insect folivory tracks paleotemperature for six million years[END_REF].

Conclusion

This study set up a snapshot of plant-insect interactions during the onset of the glacial/interglacial cycles. The Bernasso fossil leaf assemblage contained a large diversity of DTs, even if only 6 DTs represent more than 56% of the total of interactions. The most important damage evidenced is galling, but it is represented mostly by a unique DT (DT32). Based on previous studies [START_REF] Fernandes | The adaptive significance of insect gall distribution: survivorship of species in xeric and mesic habitats[END_REF][START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF]), it appears that galling is undiversified despite its abundance.

The overwhelming presence of galling at Bernasso can be explained by the community structure of a late-successional forest. Furthermore, the low diversity of galling could suggest that the mean annual precipitation was relatively high as indicated in previous studies [START_REF] Suc | Analyse pollinique de dépots Plio-Pléistocènes du sud du massif basaltique de l'Escandorgue (Site de Bernasso, Lunas, Hérault, France)[END_REF][START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF] and the temperatures was cooler than modern day and could confirm the recent unpublished results (unpublished data) Finally, a Pliocene locality from southeast France (such as Armissan [Aude, France]; Saporta, 1893) would make for an interesting comparison to Bernasso and to bring out evolutionary trends or perturbation of plant-insect interactions due to installation of glaciation cycles. Also, it could be interesting to study galling patterns in the current Mediterranean and Caucasian forests to compare with Bernasso paleoforest. This level of comparison should highlight links between galling, temperature and precipitation and therefore, by extension, relationships between paleoclimatic fluctuations and changes in plant-insect interactions.

Acknowledgments

We would like to give special thanks to Samantha Moody (Bonn University, Germany) for her considerable support in improving the English grammar of this manuscript. Thank you to Anthony Jijina (Laramie, WY) and Gentry Gatlett (Oxford, OH) for their help in the understanding of the article. The present study was financially supported by grants of the Deutsche Forschungsgemeinschaft (DFG, grant no WA 1492/8-1;11-1). This is a publication of the Institut des Sciences de l'Evolution de Montpellier (Unité Mixte de Recherche 5554 du Centre National de la Recherche Scientifique) no. 200x-00x. [START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF]. of leaves with galling for the bigger size, and less than 20% for the smaller.

Information about galling was obtained comparing Bernasso Galling and data from Fernandes and Price (1996) and [START_REF] Fernandes | Gall-inducting insect species richness as indicators of forest age and health[END_REF]. The symbols + and -respectively means high and low.

3 ecological parameters were used: {Left to the right} Age of Bernasso's forest, Mean annual temperature at Bernasso (MAT), Mean annual precipitation at Bernasso (MAP).

Figure 1 Figure

 1 Figure 1

Figure 1 :

 1 Figure 1: Map of the Hérault region showing the fossil locality of Bernasso (simplified from[START_REF] Leroy | Latest Pliocene pollen and leaf floras from Bernasso palaeolake (Escandorgue Massif, Hérault, France)[END_REF].

Figure 2 :

 2 Figure 2: Representative insect damages on Bernasso leaves. A-Galling (DT32) on Parrotia persica preserved by Canada balsam method (BDC fossil); enlarged in B, DT32 are the small black dot, plenty on this leaf. C-Galling on Parrotia persica preserved by (Ro fossil); enlarged in D. E-Skeletonization (DT20) on Acer opalus; F-Hole (DT78) on Carpinus orientalis; G-Margin feeding (DT12) on Parrotia persica; H-Surface feeding (DT130) on Carpinus orientalis; enlarged in I; J-Mining (DT90) on Zelkova ungeri; enlarged in K. Scale bars: black= 0.5cm / stripes= 0.25cm.

Figure 3 :

 3 Figure 3: Interactions distribution (FFG) according Bernasso plant species. Box on the right is a highlight of the big one behind. Species in bold are the main 7 species with more than 20 specimens. Due to the diversity of maples, the genus Acer was individualized.

Figure 4 :

 4 Figure 4: Mean diversity of all damage types, specialized damage types and galling damages on Bernasso floras. The box-plots are based on the rarefaction data for each group of interactions considered for the database.

Figure 5 :

 5 Figure 5: Distribution of FFGs according to the fossil preparation. The grey bars represent the fossils preserved with Canada balsam method; the black bars represent the fossils preserved on the rock. The numbers are the exact quantity of damages. Chi-squared test enables to see the significance of the difference between Ro and BDC per FFG.***: alpha =0.1%. Lack of "***" means there is no significant difference.

Figure 6 :

 6 Figure 6: Overview of Bernasso forest. The different trees represent the 7 main species with more than 20 specimens. The purple circles near each tree represent the galling. The 2 different sizes of the circles indicate that there are more than 20%