M N Vu

S T Nguyen

M H Vu

A M Tang

V T To

Heat conduction and thermal conductivity of 3D cracked media

Keywords: thermal conductivity, crack, insulator, MR-BEM, self-consistent

This study deals with the heat conduction within a medium containing cracks that are assumed to be perfect insulators. Multi-region boundary element approach is employed to obtain a boundary singular integral equation governing the steady state thermal transfer within this medium. This equation presents the temperature field within the whole cracked body as a function of temperature and rate of heat flow on the domain 's boundary and temperature discontinuity across the cracks. For the particular case of an infinite domain under far-field condition, the temperature field solution is only a function of the cracks temperature's discontinuity. The basic problem of a single crack in an infinite domain is investigated and a closed-form solution is derived for a crack of elliptical plane from this analysis. This solution is the key issue to estimate the effective thermal conductivity of the whole domain by coupling with the classical homogenization schemes. The arbitrary crack form is covered up by using the excluded volume definition. Estimations of effective thermal conductivities stemming from diluted, differential and self-consistent approaches are compared to numerical solution obtained by the finite volume method that is available in literature. This comparison shows that the self-consistent scheme is the most appropriate model to estimate the thermal conductivity of materials containing cracks.

Introduction

The effects of cracks on the thermal conductivity have been investigated on various materials such as optical coatings [START_REF] Redondo | Thermal conductivity of optical coatings[END_REF], graphite materials [START_REF] Pestchanyi | Effective Thermal Conductivity of Graphite Materials with Cracks[END_REF], ceramics [START_REF] Papitha | Effect of micro-cracking on the thermal conductivity and thermal expansion of tialite (Al2TiO5) ceramics[END_REF], thermal barrier coatings [START_REF] Nicholls | Methods to reduce the thermal conductivity of EB-PVD TBCs[END_REF][START_REF] Cernuschi | Modelling of thermal conductivity of porous materials: application to thick thermal barrier coatings[END_REF][START_REF] Sevostianov | Quantitative characterization of microstructures of plasma -sprayed coatings and their conductive and elastic properties[END_REF], fibre-reinforced composites [START_REF] Lu | Thermal conductivity and expansion of cross -ply composites with matrix cracks[END_REF], clay [START_REF] Dao | Anisotropic thermal conductivity of natural Boom Clay[END_REF] and among others. These studies show that the thermal conductivity of material s decreases with the presence of cracks due to the very low thermal conductivity of air void cr eated by cracks. The presence of cracks was evidenced by micrograph [START_REF] Papitha | Effect of micro-cracking on the thermal conductivity and thermal expansion of tialite (Al2TiO5) ceramics[END_REF][START_REF] Cernuschi | Modelling of thermal conductivity of porous materials: application to thick thermal barrier coatings[END_REF], topological sensitivity [START_REF] Bonnet | Sensitivity analysis for shape perturbation of cavity or internal crack using BIE and adjoint variable approach[END_REF][START_REF] Bonnet | Higher-order topological sensitivity for 2-D potential problems. Application to fast identification of inclusions[END_REF] and was quantitatively characterized to interpret the results on thermal conductivity [START_REF] Sevostianov | Quantitative characterization of microstructures of plasma -sprayed coatings and their conductive and elastic properties[END_REF].

X-ray micro-tomography can be considered as a powerful tool to monitor the 3D spatial distribution of cracks in materials [START_REF] Stock | Recent advances in X-ray microtomography applied to materials[END_REF].

To model the heat conduction within materials containing cracks of complex geometry and to take into account the interaction between the cracks, the numerical method imposes naturally. The use of volume discretisation methods, such as finite element method [START_REF] Pouya | Determination of rock mass strength properties by homogenisation[END_REF][START_REF] Saevik | A 3D Computational Study of Effective Medium Methods Applied to Fractured Media[END_REF], finite volume method [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF][START_REF] Mourzenko | Three-dimensional numerical simulation of single-phase transient compressible flows and well -tests in fractured formations[END_REF][START_REF] Bogdanov | Effective permeability of fractured porous media with power-law distribution of fracture sizes[END_REF][START_REF] Mourzenko | Percolation and permeability of fracture networks in excavated damaged zones[END_REF] is confronted to a major difficulty in mesh generation process when the domain contains a numerous cracks of randomly geometries. Contrariwise, boundary element method (BEM) presents an important advantage in mesh reduction [START_REF] Katsikadelis | Boundary elements -Theory and Applications[END_REF][START_REF] Brebbia | Boundary Elements: An Introduction Course[END_REF][START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF]. In practice, the classical BEM becomes also inefficiently for high crack density by the fact of complicated numerical implementation compared to standard BEM for the intact domain.

Atalay et al [START_REF] Atalay | Multi-region heat conduction problems by boundary element method[END_REF] employed the multi-region boundary element method (MR-BEM) to modeling the steady-state heat conduction in a heterogeneous system constituted by several homogeneous piecewise functions. This technique has been also applied for other types of heterogeneous media such as materials in the presence of cracks, in particular in the linear fracture mechanics [START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF][START_REF] Shen | Multi-Region Boundary Element Analysis for Coupled Thermal -Fracturing Processes in Geomaterials[END_REF]. However, only a single or a few cracks are treated and the crack intersection is not considered in this application field. Recently, Pouya [START_REF] Pouya | Three-dimensional flow in fractured porous media: a potential solution ba sed on singular integral equation[END_REF] and Pouya and Vu [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF][START_REF] Pouya | Fluid flow and effective permeability of an infinite matri x containing disk-shaped cracks[END_REF] proposed the closed-form expressions of mass exchange between matrix and cracks at crack intersections for fluid flow modeling. Then, Vu [START_REF] Vu | Modélisation des écoulements dans des milieux poreux fracturés par la méthode des équations intégrales singulières[END_REF] associated this advanced theoretical formulation with MR -BEM to develop the integral boundary equation and its numerical resolution for modeling the fluid flow within a porous domain containing numerous intersecting cracks. This method is now employed to study heat conduction within the cracked media where the crack filled by air whose conductivity is neglected in compared to one of its surrounded material. Therefore, the cracks could be considered as a barrier to the heat flow that is naturally contrary to the case of fluid flow where the crack is more conductive than the porous embedding matrix.

The present work focuses on the steady-state heat flow within a three dimension finite domain containing numerous cracks and its effective thermal conductivity. We make assumption that crack has zero thickness and an infinitely cross thermal resistance, i.e. it acts as a barrier to heat flow. Therefore, the rate of heat flow vanishes on the crack surfaces in its normal direction; and there exists a temperature discontinuity across it. MR-BEM is employed to derive the boundary integral equation that governs heat conduction through the whole domain presenting a high crack density. For this end, the domain is split into several sub-domains by existing cracks and fictitious surfaces, such that a fictitious surface does not intersect cracks and other fictitious surfaces. The unknown temperature on the each face of crack is considered as the boundary condition for the sub-domains. Boundary integral formulation for fluid flow is written for each sub-region and then assembled to obtain only one potential solution for the whole body. In this solution, the temperature field is expressed as function of temperature and heat flow on the domain's boundary and temperature discontinuity across the cracks. This theoretical solution is used subsequently to estimate of effective thermal conductivity of the material with the presence of void cracks. For this purpose, infinite domain is considered. The closed-form solution of temperature discontinuity is derived from the singular integral equation writing for the case of a crack with elliptical shape. The arbitrary crack shape is covered in this study by having recourse to excluded volume definition [START_REF] Adler | Fractures and Fracture Networks[END_REF][START_REF] Charlaix | A criterion for percolation threshold in a random array of plates[END_REF][START_REF] Mourzenko | Trace analysis for fracture networks with anisotropic orientations and heterogeneous distributions[END_REF][START_REF] Mourzenko | Permeability of isotropic and anisotropic fracture networks, from the percolation threshold to very large densities[END_REF]. The irregular shape of cracks could be also taken into account by analysing the strength of the heat flux singularity in vicinity of crack front [START_REF] Sevostianov | Thermal conductivity of a material containing cracks of arbitrary shape[END_REF]. The analytical solution of temperature discontinuity allows estimating the effective thermal conductivity in the framework of various homogenization schemes: dilute, differential and self-consistent. The details of these techniques for thermal properties of porous geomaterial were synthesized in Do et al. [START_REF] Giraud | Effective thermal conductivity of transversely isotropic media with arbitrary oriented ellipsoïdal inhomogeneities[END_REF][START_REF] Gruescu | Effective thermal conductivity of partially saturated porous rocks[END_REF][START_REF] Do | Application des approches d'homogénéisation à l'étude des propriétés thermohydro-mécaniques des roches[END_REF], To et al. [START_REF] To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF][START_REF] To | Closed-form solutions for the effective conductivity of twophase periodic composites with spherical inclusions[END_REF], Nguyen [START_REF] Nguyen | Micromechanical approach for electrical resistivity and conductivity of sandstone[END_REF] and Chen [START_REF] Chen | Estimating effective thermal conductivity of unsatur ated bentonites with consideration of coupled thermo -hydro-mechanical effects[END_REF]. Effective properties of micro-cracked viscoelastic materials are also successfully estimated by coupling the solution of a single crack in infinite domain and the homogenization schemes (Nguyen et al. [START_REF] Nguyen | A Burger model for the effective behavior of a microcracked viscoelastic solid[END_REF][START_REF] Nguyen | Generalized Kelvin model for micro-cracked viscoelastic materials[END_REF]). The numerical solution obtained by Bogdanov et al. [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF], based on the finite volume approach, is finally used to evaluating the accuracy of these three schemes for the material that contains insulating cracks.

Notations

Thermal conduction equations

Let us consider a three-dimensional domain Ω embedding m cracks numbered by j and denoted by Γ j (j=1,m) (Fig. 1). In the mathematical model, cracks are represented by a smooth surface function z j (s) from 3 → 2 , where s = (s 1 ,s 2 ) is the curvilinear parameters. Γ designates the set of cracks Γ = Γ j . The matrix without cracks corresponds to Ω-Γ that is assumed to be homogeneous with a conductivity tensor . The heat conduction in the matrix is governed by Fourier's law:

xΩ-Γ ( ) . ( ) q x T x     (1)
where q(x) and T(x) are local rate of heat flow and temperature fields, respectively.

The conservation of energy in the matrix reads:

xΩ-Γ . ( ) 0 qx  (2)
Substituting (1) into (2) yields Laplace's equation for the case of an isotropic matrix, i.e = :

xΩ-Γ ( ) 0 Tx  (3)
Solving this equation gives the solution of steady state heat diffusion in the matrix -. This solution must verify the condition on the cracks and the boundary conditions. As the cracks are non-conducting and act as barriers to thermal flow, the conditions on the crack surfaces are given as:

xΓ ( ). ( ) 0 q x n x  or ( ). ( ) 0 T x n x  (4)
where n(x) is the unit normal vector at point x on cracks.

Two boundary conditions, namely Dirichlet and Neumann conditions are prescribed respectively on  1 ,  2 of the domain boundary  (= 1 + 2 ):

1 x   0 ( ) ( ) T x T x  (5) 2 x   0 - ( ). ( ) ( ) T x n x q x   (6)
3 Potential solution

Multi-region boundary elements approach

This sub-section is devoted to establish a singular integral equation for the temperature field in the whole domain  that satisfies the Laplace equation ( 3), the condition on crack surfaces (4) and the boundary conditions ( 5), [START_REF] Sevostianov | Quantitative characterization of microstructures of plasma -sprayed coatings and their conductive and elastic properties[END_REF]. For this purpose, the MR-BEM is had recourse. As mentioned in the introduction section, this approach exhibits a great advantage in mesh generation term for three-dimensional domain in the presence of numerous surfaces of discontinuities.

The domain  is first partitioned into n sub-domains  i by the fictitious boundaries that pass through the crack parts, as seen in Fig. 2. A crack  of normal n is constituted by two parallel surfaces  + and  -with two normal vectors n + , n -(n + = -n -= n) (Fig. 3). On the crack surface, the temperature is discontinuous between its two faces and the normal heat rate vanishes, i.e. q.n = 0. As sub-domain  i does not contain the crack, therefore, equation (3) governs the heat conduction within it. Applying the BEM procedure [START_REF] Katsikadelis | Boundary elements -Theory and Applications[END_REF][START_REF] Brebbia | Boundary Elements: An Introduction Course[END_REF][START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF] for this Laplace's equation yields:

Ω i ( ) ( ) ( ) ( , ). ( , ) ( ). i i i c T T x G x nds G x T x nds             (7)
where n is the normal to the subdomain boundary  i ; c i () a coefficient that depends on the position of  relative to the subdomain

 i : c i ()=1 if ii      ; c i ()=0 if i    and 0<c i ()<1 if i    .
For the latter case when the field point on the boundary of subdomain, the coefficient c i depends on the local boundary geometry at this point and it takes the value of 1/2 for smooth boundary. In general, this coefficient can be determined by the following equation:

Ω i 0 ( ) lim ( , ) i S c G x dx       (8)
where S  is an infinitesimal spherical surface of centre  and radius  enclosed in the solid  i .

The fundamental function G(x,) and its derivations for three-dimensional steady-state heat transfer are expressed such as [START_REF] Katsikadelis | Boundary elements -Theory and Applications[END_REF][START_REF] Brebbia | Boundary Elements: An Introduction Course[END_REF][START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF]: ( , ) ( ).

1 ( , ) 4
s s i i i c T T x G x nds G x T x nds             1 1 ( ) ( , ). ( , ) ( ). i i k k n n T x G x n ds G x T x n ds                   8 2 2 ( ) ( , ). ( , ) ( ). i i k k n n T x G x n ds G x T x n ds                   (10) ' ' 3 3 ( ) ( , ). ( , ) ( ). i i k k n n T x G x n ds G x T x n ds                 ' ' 4 4
( ) ( , ). ( , ) ( ).

i i k k n n T x G x n ds G x T x n ds                
As a matter of fact, the field point  can be located inside, outside or on the boundary of domain  i . The relative position of this point relative to the subdomain is represented by the factor c i . Therefore, we can repeat the equation ( 10) for all subdomains but for a field point  to obtain n boundary integral equations. The sum of these n equations results in:

1 ( ) ( ) ( ) ( , ). ( , ) ( ). n i i c T T x G x nds G x T x nds               (11) 1 ( ) ( ) ( , )
.

k m k T x T x G x nds           
To show how we obtain the equation ( 11), we consider a fictitious surface  ij that are shared by two subdomains  i and  j . Across the fictitious surface, the temperature and its gradient or thermal flux rate are continuous, therefore, the seventh, eighth, ninth, tenth terms in the boundary integral equation writing at the field point  for the subdomain  i will be eliminated by the same terms in the one writing for its adjacent subdomain  j . Moreover, the rate of heat flow is none on the crack surface, thus the fourth and sixth terms in the boundary integral equation for each subdomain are vanished.

By noting

1 n i i cc    , it
is easy to see that the coefficient c is always equal to 1 when the field point is inside the whole cracked domain. Therefore, the factor c of the cracked domain  can be represented by one of the un-cracked domain bounded by .

Introducing the temperature discontinuity across the cracks into equation [START_REF] Stock | Recent advances in X-ray microtomography applied to materials[END_REF] yields the boundary integral equation for temperature field within the cracked media:

( ) ( ) ( ) ( , ). ( , ) ( ). c T T x G x nds G x T x nds             (12) 1 [[ ( ) ( , ). k m k T x G x nds       
Let us consider now the particular case of an infinite medium surrounding cracks under the far-field condition T  (x). As the field point  always inside the domain, thus the boundary factor c is equal to 1. The equation ( 12) reduces to:

1 ( ) ( ) [[ ( ) ( , ). k m x k T T T x G x n ds            (13)
In heat conduction modelling point of view, numerical resolving equation ( 12) covers equation ( 13) by making a large domain's boundary, where the far-field condition T  (x) is prescribed, compared to cracked zone. The heat conduction solution is used then to calculate the effective thermal conductivity of cracked domain. However, the direct derivation of equation ( 13) that reduce the complication of numerical implementation in comparison with the one of equation( 12), allows also the determination of thermal conductivity for a cracked zone [START_REF] Pouya | Numerical modelling of steady-state flow in 2D cracked anisotropic porous media by singular integral equations method[END_REF].

The heat conduction solution [START_REF] Pouya | Determination of rock mass strength properties by homogenisation[END_REF] obtained under a singular integral equation form is suitable for both numerical simulation and analytical development. This provides an interesting material that can be useful for those who performs numerical simulation based on boundary element method by beneficing the mesh reduction from 3D to 2D. Besides, this solution allows deriving the closed-form solution for heat conduction in some particular cases. As showed subsequently, the analytical resolution of equation [START_REF] Saevik | A 3D Computational Study of Effective Medium Methods Applied to Fractured Media[END_REF] in the case of a single fracture within an infinite medium, coupling with homogenization methods affords us another ways to estimate the thermal conductivity of material with the presence of cracks.

Closed-form solution of temperature discontinuity

This section is devoted to derive the closed-form solution of temperature discontinuity across the a single crack with in an infinite matrix under far -field condition T ∞ ()=A.. Deriving equation ( 13) with respect to  variable reads:

2 1 ( ) ( ) ( , ) [[ ( ) k m k x T T G x T x ds n n n n                      (14)
Considering the field point  located on the fracture surfaces at curvilinear parameter s, the left hand side of equation ( 14) vanishes according to the condition (4). Thus, the equation ( 14) becomes:

( ). ( ) ( ( ). ( , ))( ( ). ( , ))

( ) ( ) 3 4 ( , ) ( , ) j j n s n x n s R s x n x R s x x s ds r s x r s x             (15) with () () x Tx x n     ; ( ) [[ ( ) x T x    ( 16
)
This equation is a hyper--singular integral equation in which the temperature discontinuity on the fractures (x) is unknown function and the hypersingular kernel behaves as r -3 at the singular point (s=x). The hyper-singular kernels in equation ( 15) are also arisen in symmetry Galerkin boundary element approximation for solving tw o or three-dimensional Laplace's equations. The algorithms based on limit process are presented in Sutradhar et al. [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF] for evaluating these hyper-singular integrals. In such process, the integral involving derivatives of the fundamental solution are defined as limits from the interior.

In order to upscale the effective permeability, a basic problem is considered. This is to derive the temperature field around a single insulating crack within an infinite matrix under far-field condition. The crack is modelled by an elliptical form D with halfdiameters designated by

d 1 , d 2 (d 1  d 2 )
and two corresponding principal directions e 1 , e 2 denote elliptical fracture principal directions and its unit normal e 3 . A local coordinate system is defined such as the origin located at the centre of the crack and its axes parallel to e 1 , e 2 and e 3 . The equation of ellipse D in this coordinate system reads (Fig. 4):

x.B.x = 1 [START_REF] Mourzenko | Percolation and permeability of fracture networks in excavated damaged zones[END_REF] with x.e 3 = 0 or x=(x 1 ,x 2 ,0) and

 1 = 1/d 1 2 ,  2 = 1/d 2 2 , B =  1 e 1 e 1 +  2 e 2 e 2 (18)
Fig. 4 A elliptical penny-shaped crack within an infinite matrix For a single planar crack, ( ). ( , ) ( ). ( , ) 0 n s R s x n x R s x  and n(s).n(x) = 1, thus, equation ( 15) is rewritten as follows:

xD 12 1 2 1 2 3/2 22 1 1 2 2 ( , ) 1 ( , ) 4 ( ) ( ) D ss x x ds ds s x s x            ( 19
)
This equation is equivalent to:

xD 1 2 1 2 1 2 1/2 22 1 1 2 2 1 4 ( , ) ( , ) ( ) ( ) s D x x s s ds ds s x s x            (20)
Since the temperature is continuous at the crack border, i.e. (x 1 ,x 2 ) subject to the following condition:

x∂D 12 ( , ) 0 xx  (21)
Integrating by part and the condition (21) lead to:

xD 1 2 1 2 1 2 1/2 22 1 1 2 2 1 4 ( , ) ( , ). ( ) ( ) D x x s s ds ds s x s x           (22)
In the case that the far-field condition is a constant pressure gradient A, i.e T  (x)=A.x, hence,

( ) . x Ae 

, equation ( 22) is rewritten as:

x D 3 3 1 2 1 2 1/2 22 1 1 2 2 1 4 . ( , ). ( ) ( ) D A e s s ds ds s x s x           ( 23
)
This equation is well-known in fracture mechanic (Bui [42]; Guidera and Lardner [START_REF] Guidera | Penny-shaped crack[END_REF]) and its solution is expressed as follows:

x D 2 3 3 2 ( ) . 1 . . () d x A e x x E     B ( 24
)
where

2 2 1 1       d d
and E() is the complete elliptic integral of the second kind:

  /2 22 0 1 s i n          Ed (0   1) (25)
The semi-analytical or iterative numerical expression of the complete elliptic integrals can be found in mathematical handbooks of Abramowitz and Stegun [START_REF] Abramowitz | The Process of the Arithmetic-Geometric Mean[END_REF].

Considering the particular case of a circle of radius d 1 = d 2 = R, equation ( 24) reads:

x D 2 3 4 ( ) . . x A e R x x     (26)
By noting  the integral of discontinuity temperature function over the crack, we obtain:

3 8 [[ ( ) . ( ). 3 D T x ndx R n n A       (27)
This closed-form solution is an important element to evaluate the effective thermal conductivity of cracked media by using the homogenization approaches. Three schemes, namely dilute, differential and self-consistent, will be considered in the section 4.

Thermal conductivity

This section presents the homogenisation schemes to estimate the effective thermal conductivity of a medium embedding cracks that are non-conducting and acts as barriers to the heat flow. For this purpose, a domain , containing several cracks  j prescribed on its boundary by a linear temperature T = A.x, i.e. a constant temperature gradient T=A, is considered.

The average rate of heat flow and average temperature gradient over  are formulated as:

1 () Q q x dx     ; 1 ( ). ( ) G T x n x dx A     (28)
The effective thermal conductivity  eff is given globally by Fourrier's law:

.

eff QG  (29)
The divergence theorems for the temperature function in a domain containing discontinuities Γ j reads:

. [[ .

j j T dx T n dx T n dx            (30)
Substituting equations ( 1) and ( 30) into the average rate of heat flow Q yields

1 [[ . j j Q A T n dx              (31)
Equations ( 29) and ( 31) implies the reduction of thermal flow rate, by the temperature

discontinuities across cracks 1 [[ . j j T n dx     
, is a linear function of average gradient of temperature G, i.e. there exists a tensor  such as:

1 [[ . . j j T n dx G        (32)
Replacing Q in equation ( 29) by its expression in equation ( 31) and [START_REF] Giraud | Effective thermal conductivity of transversely isotropic media with arbitrary oriented ellipsoïdal inhomogeneities[END_REF] gives

() eff      (33)
Therefore, the effective thermal conductivity can be determined when the temperature discontinuity on cracks are known. For the general case, solving equation ( 12) by numerical methods, such collocation method [START_REF] Bonnet | Equations Intégrales et Éléments de Frontière[END_REF], Symmetry Garlerkin Boundary Element method [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF], etc, allows obtaining the solution of temperature discontinuity field on cracks. Besides, the effective thermal conductivity  eff can be estimated by using the theoretical up-scaling scheme based on the analytical solution [START_REF] Adler | Fractures and Fracture Networks[END_REF]. Three schemes are considered in the following and its results are then analyzed and compared with numerical solution in the literature.

Dilute scheme

Let us consider m fracture families are embedded in a porous domain . Each family Γ i has a density  i = N i / (i.e N i fractures within ), orientation n i (unit normal vector) and radius R i (Fig. 5). In the dilute framework, each crack is supposed to be subjected to the same far-field condition, i.e. undisturbed by the presence of other cracks. This approach presents a simple and accurate estimation of equivalent properties for material containing inhomogeneities when the crack density is small. By neglecting fracture interaction, the term  in equation ( 32) is expressed such as:

3 8 () 3 m ii ii i R n n      (34)
The dimensionless crack density is introduced to cover other fracture shapes:

0 i i i V    (35)
In this relation, V 0i is a reference volume around a fracture of family i. Since the effective conductivity resulted in the current work will be compared to one obtained by Bogdanov et al. [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF], we take thus 23 0ii

VR 

to unify the notation of dimensionless crack density between these two work. According to Adler and Thovert [START_REF] Adler | Fractures and Fracture Networks[END_REF] and Charlaix et al. [START_REF] Charlaix | A criterion for percolation threshold in a random array of plates[END_REF], for an isotropic distribution of fractures having the same circular shape of radius R, the volume 23 ex VR  is the excluded volume that is defined as one around an object of the system where other objects do not appear in order not to intersect it. As the exclusion volume characterizes all objects in the system as an ensemble , i.e. it is not given for each family independently, the volume V 0i does not have sense of excluded volume with the presence of several fracture families with different sizes and orientations. A general expression of V ex for an anisotropic system could be referred in Mourzenko et al. [START_REF] Mourzenko | Trace analysis for fracture networks with anisotropic orientations and heterogeneous distributions[END_REF][START_REF] Mourzenko | Permeability of isotropic and anisotropic fracture networks, from the percolation threshold to very large densities[END_REF], where its calculation includes orientations, shapes and sizes of all fractions within the fracture network.

Substituting [START_REF] To | A numerical-analytical coupling computational method for homogenization of effective thermal conductivity of periodic composites[END_REF] into (34) yields:

2 8 () 3 m ii i i nn       (36)
For the case of an isotropic distribution of fractures with density  having identical radius R, the average value 1 3 nn   , the effective conductivity tensor is isotropic  eff = eff , and: eff      [START_REF] Nguyen | Micromechanical approach for electrical resistivity and conductivity of sandstone[END_REF] This formulation is identical to one given by Shafiro and Kachanov [START_REF] Shafiro | Anisotropic effective conductivity of mater ial with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF].

Differential scheme

It is worth noting that the crack density  is a discrete function, however, the dimensionless crack density  is a continuous function since the crack size, i.e R varies from zero. Therefore, differential scheme can be applied to upscale the effective thermal conductivity of cracked materials. This method consists in estimate the thermal conductivity explicitly from an initial material through a series of incremental additions [START_REF] Mclaughlin | A study of the differential scheme for composite materials[END_REF][START_REF] Norris | A differential scheme for the effective moduli of composites[END_REF]. Supposing d the adding crack density, the effective thermal conductivity at a considered step is determined from one at the previous step by equation [START_REF] Nguyen | Micromechanical approach for electrical resistivity and conductivity of sandstone[END_REF]:

1 2 8 1 9 ii d           (38)
Furthermore,

1 2 8 9 ii i d        (39)
and then,

Self-consistent scheme

As seen in equation [START_REF] Sevostianov | Thermal conductivity of a material containing cracks of arbitrary shape[END_REF] the average rate of heat flow Q is reduced by Q * =A vis-a-vis one of un-cracked medium. Regardless the two previous approaches, the self-consistent scheme considers the cracks surrounding by matrix having the effective conductivity. In other words, Q * = eff A. Replacing Q * by this expression into [START_REF] Sevostianov | Thermal conductivity of a material containing cracks of arbitrary shape[END_REF] and then (29) yields:

2 1 8 1 9 eff      (42)
For the case  << 1, the first order expansion of equation ( 42) results the formulation [START_REF] Nguyen | Micromechanical approach for electrical resistivity and conductivity of sandstone[END_REF].The results stemming from differential scheme (equation ( 41)) and self-consistent one (equation ( 41)) are compared with numerical result obtained by Bogdanov et al. [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF]. This comparison allows choosing which theoretical approach is better in the case of media containing insulating crack.

Bogdanov et al. [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF] investigated the effective property of finite porous media embedding an isotropic distribution of hexagonal fracture network based on a threedimensional finite volume method. Consider the insulating fracture have a conductivity << 1 and a cross section resistance >>1. A least squares fit of all numerical results allows them to introduce an expression of effective property versus  as follows: 2 1.0010 0.0604 0.0014

eff        (  8) (43)
According to equations ( 37), ( 41), ( 42), [START_REF] Guidera | Penny-shaped crack[END_REF], variations of effective thermal conductivities as function of crack density  are presented in Fig. 6. Obviously, two approximated curves, namely differential model [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF] and self-consistent one [START_REF] Bui | Mécanique de la rupture fragile[END_REF], are tangent to the curve of Shafiro and Kachanov [START_REF] Shafiro | Anisotropic effective conductivity of mater ial with nonrandomly oriented inclusions of divers ellipsoidal shapes[END_REF] at =0. It can be note that the selfconsistent model [START_REF] Bui | Mécanique de la rupture fragile[END_REF] is better than differential one [START_REF] Sutradhar | Symmetric Galerkin Boundary Element Method[END_REF] and diluted one [START_REF] Nguyen | Micromechanical approach for electrical resistivity and conductivity of sandstone[END_REF] in comparison with numerical model of Bogdanov et al. [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF]. In practical point of view, the characterization of crack system (size, distribution, orientation, connectivity) within a material is not an easy task, even impossible. Besides, the measured thermal conductivity on sample in laboratory and particularly in -situ is usually supposed to be isotropic, i.e. a randomly crack distribution. Therefore, the assumption of simple crack shape, as well as an isotropic cracks distribution is habitually made to simulate the experiment measurement and the real engineering processes. As a conclusion, although the current work considers a simple case of a system of equally sized disc shapes with isotropic distribution for cracks, but the selfconsistent model [START_REF] Bui | Mécanique de la rupture fragile[END_REF] with the validation against to numerical model [START_REF] Bogdanov | Effective permeability of fractured porous media in steady state flow[END_REF] seems to be very useful in engineering fields. This model could be applied to estimate rapidly the thermal conductivity of cracked media or to determine the crack parameters from the measurement by inverse analysis.

Conclusion

Steady state heat transfer within media containing cracks, which act as barrier to heat flow, is investigated in this paper. The cracks are mathematically modelled by a two dimensional smooth surface of zero thickness. Multi-region boundary element method is used to obtain the potential solution that presents the temperature solution in the whole domain as function of thermal variables on the boundary, namely temperature and rate of heat flow, and temperature discontinuity on the cracks. This boundary integral equation allows reducing the dimension problem from 3D to 2D in numerical analysis. Fast multipoles method [START_REF] Liu | Fast Multipole Boundary Element Method: Theory and Applications in Engineering[END_REF] dealing with this equation is on-going to model accurately the heat transfer within this material.

Applying the boundary integral equation to an infinite medium that contains a crack of elliptical shape, allows deriving the closed-form solution of temperature discontinuity across this crack surface. Based on this solution, three homogenization schemes (dilute, differential and self-consistent) are employed to estimate the effective thermal conductivity. The excluded volume is had recourse to extend the arbitrary form of crack. These three models are compared to numerical solution proposed in the literature to evaluating the validity of homogenization estimations. This comparison shows that the self-consistent model is the most appropriate to estimate the thermal condu ctivity of material containing insulating cracks.

Fig. 1 :

 1 Fig. 1: Heat conduction within a medium embedding cracks .

Fig. 2 :Fig. 3 :

 23 Fig.2: MR-BEM approach to study heat conduction within cracked media

  with R=x; r = R =  -x, nx and n  the normal to crack surface at x and  The integrals on the right hand side of equation (7) are the sum of this integral on the different parts of boundary of this subdomain, i.e. on the s i  , n 1i positive fracture surface segments k   , n 2i negative fracture surface segments k   , n 3i fictitious positive lines ' k   and n 4i fictitious negative lines '

35

 35

Fig. 5 :

 5 Fig. 5: Medium containing a weak crack density.

Fig. 6 .

 6 Fig. 6. Comparison homogenisation models of effective thermal conductivity of cracked medium (developed within this paper) with results obtained by using the finite volume approach proposed by Bogdanov et al. (2003).

Acknowledgements

The authors wish to thank the anonymous reviewers for their constructive comments on an earlier version of this article.


Yoann Poher
email: yoann.poher@imbe.fr

Philippe Ponel

Frédéric Guiter

Valérie Andrieu-Ponel

Frédéric Médail

Erosion of insect diversity in response to 7000 years of relative sea-level rise on a small Mediterranean island

Keywords:

We have investigated the potential effects of global sea-level rise on Mediterranean coastal wetlands by studying the Coleoptera and pollen fossil remains in a 7000-year sedimentary record, which we obtained from a coastal marshy area on a small Mediterranean island (Cavallo, southern Corsica). Using beetle structural diversity and plant composition as recorded prior to marine and human influences as a 'past analogue', we reconstructed the impact of the Holocene relative sea-level rise on the coastal ecosystem. Our results show that beetle species richness and diversity were highest when freshwater was predominant, which was the case until about 6200 years ago. We also found that a major increase in salinity had occurred over the last 5300 years, experiencing a peak rate of increase at about 3700 years ago. These changes are clearly reflected in the fossil records of the following key taxa: halophilous beetles (Ochthebius sp., Pterostichus cursor), halophilous plants (Chenopodiaceae, Tamarix) and non-pollen palynomorphs (microforaminiferal linings). In particular, we note that the majority (60%) of wetland beetle fauna became locally extinct in response to the salinity changes, and these changes were exacerbated by the recent aggravation of human pressures on the island. The major part of this diversity loss occurred 3700 years ago, when the relative Mediterranean sea-level rose above -1.5 ± 0.3 meters. These findings demonstrate the value of fossil beetle assemblage analysis as a diagnostic for the response of coastal wetland biodiversity to past salinity increases, and serve as a means of forecasting the effects of sea-level rise in the future. The conservation of inland freshwater bodies could ultimately prove essential to preserving freshwater insect diversity in threatened coastal environments.

Introduction

Coastal ecosystems such as wetlands are major biodiversity hotspots. Despite being the subject of intense research, the full extent to which these fragile environments are vulnerable to future global sea-level rise remains poorly understood. Mediterranean coastal wetlands are particularly challenging to study because they are subject to threats from human activities and land-use changes [START_REF] Beltrame | Mediterranean wetlands: outlook. First Mediterranean wetlands observatory report-Technical report-2012[END_REF]) in addition to relative sea-level rise (RSLR) induced by current climate warming. RSLR is particularly damaging because it can lead to the intensification of coastal erosion and to the submersion of low elevation islands and coastal wetlands [START_REF] Nicholls | Coastal degradation as an indicator of global change[END_REF][START_REF] Wong | Coastal systems and low-lying areas[END_REF].

Fortunately, there are multiple methods by which these myriad effects can be teased apart and analysed. One approach is based on predictive global sea-level change models, and has been employed in recent studies to suggest that several oceanic island hotspots will likely suffer a great loss of terrestrial biodiversity (Bellard et al. 2013a, b;Courchamp et al. 2014). Although quite informative, such prospective studies do not consider local parameters (which may differ between islands) together with the 'historical background'. For this reason they should be thought of as providing more of a general idea rather than a highly detailed view of the situation.

A more nuanced perspective can be obtained using the alternate approach, which is to begin by reconstructing the ecosystem dynamics of Mediterranean wetlands in the distant past [START_REF] Anderson | Linking palaeoenvironmental data and models to understand the past and to predict the future[END_REF][START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF][START_REF] Willis | Biodiversity baselines, thresholds and resilience: testing predictions and assumptions using palaeoecological data[END_REF]. However, this tends to be complicated by the fact that ecological monitoring records do not typically go back more than a few decades. Therefore, palaeoecological investigations on sedimentary archives, through multi-disciplinary approach, are the key to understanding how coastal ecosystems have responded to past RSLR. This knowledge can then be used to shed light on future RSLR scenarios [START_REF] Woodroffe | Sea-level rise and coastal change: the past as a guide to the future[END_REF]. This approach, though potentially quite effective, appears relatively infrequently in the body of literature dealing with the impact of past RSLR on the biodiversity of coastal wetlands [START_REF] Virah-Sawmy | Threshold response of Madagascar's littoral forest to sealevel rise[END_REF].

In this paper we present a multi-proxy study, using beetles, pollen, sedimentology and radiometric dating to analyse a sedimentary profile from the small Corsican island of Cavallo (Lavezzi archipelago). We chose these particular proxies and this specific record because (1), fossil beetles and pollen are excellent markers of past environmental conditions and landscape dynamics, (2) Corsica has historically been the source of several paleoecological records and has a landscape whose Holocene history is relatively well known (e.g., [START_REF] Reille | Origine de la végétation actuelle de la Corse sud-orientale; analyse pollinique de cinq marais côtiers[END_REF][START_REF] Reille | New pollen-analytical researches in Corsica: the problem of Quercus ilex L. and Erica arborea L., the origin of Pinus halepensis Miller forests[END_REF], and (3) past RSLR is well documented in this part of the Mediterranean [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF][START_REF] Lambeck | Sea-level change in the Mediterranean sea since the LGM: model predictions for tectonically stable areas[END_REF][START_REF] Morhange | Changes of relative sea level during the past 5000 years in the ancient harbor of Marseilles, Southern France[END_REF][START_REF] Vacchi | Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: sea-level variability and improvements in the definition of the isostatic signal[END_REF][START_REF] Vella | Relative sea-level rise and neotectonic events during the last 6500 yr on the southern eastern Rhône delta, France[END_REF].

The overall aim of the present study is to (1) reconstruct past ecosystem dynamics and past Coleoptera diversity changes in relation to Holocene sea-level rise, (2) provide a 'past analogue' suitable for predicting the condition of the area by 2050 (Galassi and Spada 2014), and (3) propose strategies for insect conservation within the framework of future sea-level rise in the Mediterranean basin.

Materials and methods

Study site, sampling and chronology

Cavallo Island [N 41°22′06′′/E 9°15′45′′; area: 120 ha; maximal elevation: 32 meters above sea level (m.a.s.l.)] is a low-elevation island of the Lavezzi archipelago (Southern Corsica, France; Fig. 1a). The waters surrounding Cavallo are shallow, with depths rarely exceeding 10 m. A coastal retrodunal pond (Greco pond; 0.48 ha, 1 m.a.s.l.), located on the northern part of the island, was sampled using a Russian corer. A 340 cm length sediment profile (CAVA-1) was obtained. After description of its lithology, CAVA-1 was subsampled into seventy 5 cm-thick slices and 38 circa 2 cm 3 subsamples respectively for coleopteran and pollen analyses.

Six AMS radiocarbon measurements were performed on bulk peat sediment samples (Poznan Radiocarbon Dating Laboratory, Poland). Ages were calibrated to calendar years before present (cal. yr B.P.; present fixed at 1950) and an interpolated age-depth model was created using the 'CLAM' R package [START_REF] Blaauw | Methods and code for 'classical' age-modelling of radiocarbon sequences[END_REF]) from the R 2.13.2 software (R Development Core Team 2011) (Online Resources 1, 2). The fine sand layers intercalated in the peat profile (likely corresponding to instant aeolian deposits) were not taken into account in the age-depth modelling.

Coleopteran fossil analysis

Fossil insects were extracted using the standard method (Coope 1986) and Coleoptera remains were identified by direct comparison with modern specimens from a reference collection (see Online Resource 3). The term 'taxa' refers to beetle fossils identified at any taxonomic rank (e.g., species, genus, family). The palaeoecological interpretation of the fossil assemblages is based on the species' modern biological and geographical distribution derived from the literature available for the region [START_REF] Balachowsky | Faune de France 50[END_REF][START_REF] Caillol | Catalogue des Coléoptères de Provence d'après des documents recueillis et groupés, 1er partie[END_REF][START_REF] Coiffait | Coléoptères Staphylinides de la région paléartique occidentale, V, Sous famille Paederinae, Tribu Paederini 2, Sous famille Euaesthetinae[END_REF][START_REF] Constantin | Coléoptères Dasytidae de France. Société linnéenne de Lyon, Lyon Coope GR[END_REF][START_REF] Delobel | Les plantes hôtes des bruches (Coleoptera Bruchidae) de la faune de France, une analyse critique[END_REF][START_REF] Delobel | Les Coléoptères des denrées alimentaires entreposées dans les régions chaudes[END_REF][START_REF] Guignot | Effect of macrophytes on Phytophilous macroinvertebrate community: a review[END_REF][START_REF] Hoffman | Coléoptères Curculionides (première partie)[END_REF][START_REF] Hoffman | Coléoptères Curculionides (deuxième partie)[END_REF][START_REF] Hoffman | Coléoptères Curculionides (troisième partie)[END_REF][START_REF] Jeannel | Coléoptères Carabiques première partie[END_REF][START_REF] Lechevalier | Coléoptères Carabiques deuxième partie, Faune de France 40[END_REF]Koch 1989aKoch , b, 1992;;[START_REF] Paulian | Faune des Coléoptères de France[END_REF]Saint-Claire Deville 1914;[START_REF] Soldati | Faunistique des Coléoptères Tenebrionidae de Corse, résultats d'une deuxième champagne de prospections[END_REF]Thérond 1975;[START_REF] Tronquet | Climate change, sea-level rise, and conservation: keeping island biodiversity afloat[END_REF]. These data were compared with the BugsCEP database [START_REF] Buckland | BugsCEP Coleopteran ecology package[END_REF] which is a rich source of information on European beetle fauna. The nomenclature and taxonomic order are based on [START_REF] Tronquet | Climate change, sea-level rise, and conservation: keeping island biodiversity afloat[END_REF].

In order to provide a robust final dataset for supporting interpretations, a method derived from [START_REF] Poirier | Influence of hydro-sedimentary factors on mollusc death assemblages in a temperate mixed tide-and-wave dominated coastal environment: implications for the fossil record[END_REF] was applied to the initial dataset (a 70 samples 9 133 beetle taxa matrix, see Online Resource 4). The original method, which was initially used on mollusc death assemblages, was designed to obtain a statistically meaningful dataset of n stations 9 p species using selection criteria. Then the stations were classified into homogeneous groups according to the taxonomic composition of the death assemblages using hierarchical agglomerative cluster analysis. For the present study we employ the same approach and use it to investigate the taxonomic composition and temporal distribution In this category are lumped together taxa whose ecology is not sufficiently known and which do not provide any precise ecological information (not presented in Fig. 3a)

along the profile of Coleoptera fossil assemblages. Thus, samples and taxa that successively match the four following requirements were excluded: (1) taxa identified at family rank, (2), taxa that occurred in less than five samples, (3) taxa for which the sum of relative abundance did not exceed 1% and (4) samples devoid of fossil beetles. A final dataset of 68 samples × 23 taxa consisting of the most abundant and widely represented beetles and the richest assemblages was obtained. Dissimilarity matrices, as parameterised by Bray-Curtis distances, were then compiled for both samples (Samples-matrix) and taxa (Taxamatrix) using the 'rioja' R package [START_REF] Juggins | The rioja package: Analysis of Quaternary Science Data, R package, version 0[END_REF]. A stratigraphically-constrained cluster analysis was carried out on the Samples-matrix via the CONISS agglomeration method [START_REF] Grimm | CONISS: a fortran 77 program for stratigraphically constrained cluster analysis by the method of incremental sum of squares[END_REF]) to obtain significant insect fauna units (IFUs) and an unconstrained cluster analysis was performed on Taxa-matrix using the UPGMA (Unweighted Pair Group Method Average) method to obtain homogeneous ecological groups. The crosscomparison between the two clusters is given in Fig. 2. Entomological data are expressed in minimal number of individuals (MNI). Fossil beetles have been plotted in a simplified diagram (Fig. 3a, drawn with C2 software; Juggins 2007) according to ecological habitats (Table 1). The specific richness, i.e., the total number of identified beetle taxa and the Shannon diversity index [START_REF] Shannon | The mathematical theory of communication[END_REF] were calculated for each sample. We also calculated the Pielou's evenness index which is sensitive to the relative evenness of taxa in the samples [START_REF] Pielou | The measurement of diversity in different types of biological collections[END_REF]. The value is minimal when one taxon dominates the fossil assemblage from one sample and maximal when all the taxa are present in equal abundance within the assemblage.

Pollen and non-pollen palynomorph analyses

Extraction of pollen fossils and non-pollen palynomorphs (NPP including spores, microforaminiferal linings and algae such as Botryococcus and Spirogyra) from the sediment follows the standard method [START_REF] Faegri | Textbook of pollen analysis[END_REF]. Pollen and NPP counts were performed using a Leica microscope at × 500 magnification oil immersion. Identifications were made with the standard atlases and keys (Beug 2004;[START_REF] Reille | Pollen et spores d'Europe et d'Afrique du Nord, 2nd edn. Laboratoire de Botanique Historique et Palynologie[END_REF]; see Online Resource 5). A minimum of 300 pollen grains were counted, beyond this number, only new pollen types were counted. Aquatic and helophytic plants (Alisma, Cyperaceae, Lotus, Myriophyllum type spicatum, Myriophyllum type alterniflorum, Nymphaea, Sparganium-Typha type) and NPP were excluded of the total pollen sum. To facilitate the crosscomparison with insect data, some pollen and NPP have been merged into ecological categories (see Online Resource 6) and plotted in Fig. 3a.

Past Relative sea-level framework

A recent review based on 917 relative sea-level data-points provides a detailed Holocene relative sea-level framework for the western Mediterranean [START_REF] Vacchi | Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: sea-level variability and improvements in the definition of the isostatic signal[END_REF], which made it possible to re-evaluate relative sea-level data published for southeastern France [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF][START_REF] Morhange | Changes of relative sea level during the past 5000 years in the ancient harbor of Marseilles, Southern France[END_REF][START_REF] Vella | Relative sea-level rise and neotectonic events during the last 6500 yr on the southern eastern Rhône delta, France[END_REF] and for Corsica [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF]. These corrected data were plotted on Fig. 3b and compared to the predictive model of sea-level on the French Mediterranean coast from [START_REF] Lambeck | Sea-level change in the Mediterranean sea since the LGM: model predictions for tectonically stable areas[END_REF].

Results

Lithology

Sediment profile CAVA-1 is composed of two sedimentary units both interspersed with marine sandy layers (see Online Resource 2). Between 490 and 266 cm, the sediment consists of homogenous fine peat. Fine sandy peat is found deposited between 262 and 160 cm. It must be noted that the curves of total MNI and of Coleoptera specific richness show a similar trend throughout the record (Fig. 3a), which raises the question of preservation and taphonomy. This point will be addressed again later.

Beetle assemblages and diversity

From a total of 2039 arthropod fossils belonging to 144 taxa, 133 Coleoptera taxa were identified (see Online Resource 4). Thirty-nine percent of them were identified down to the species level, 48% to the genus and 13% to a higher taxonomic level. The fossil beetles were dominated by aquatics (35 taxa), by hygrophilous (32 taxa) and then by open envi-ronment and woodland/scrubland-dependent beetles (both 13 taxa). The richness varies from 0 to 29 taxa in all samples (11 taxa in average) and in terms of abundance, each taxa is represented by 1-173 individuals along the sediment core (8 individuals on average). Curves of specific richness and of Shannon diversity index have a similar trend since number of taxa and heterogeneity of the diversity are positively correlated (r = 0.8929; P < 0.00001; n = 70). The Pielou's evenness index suggests that the less diversified samples are also those whose taxa have an equal abundance within the assemblage (Fig. 3a).

The constrained cluster analyses sorted samples into four insect fauna units (Fig. 2). The three basal insect fauna units (IFU_1 to IFU_3) correspond to the lower sedimentary Fig. 2 Constrained hierarchical cluster analysis and unconstrained cluster analysis performed respectively on samples-matrix with the Bray-Curtis distance method and on the Taxa-matrix with the UPGMA aggregation method: grey-shaded areas indicate clusters of beetle groups unit (homogenous peat) and IFU 4 corresponds to the upper sedimentary unit (sandy peat). The unconstrained cluster analysis performed on Taxa-matrix allowed us to define eight main beetle assemblages (Fig. 2) detailed in Table 2. The comparison between the two cluster approaches makes it possible to characterise the temporal distribution of each taxa within the eight main beetle assemblages for each IFU. The three aquatic or semi-aquatic taxa of the 'Ochthebius assemblage' are present in each IFU with high relative abundances in IFU_1 and IFU_3. The 'Bembidion octomaculatum assemblage' is the most diversified, since it encompasses three species and species-groups that are found at the water margin (Dryops, Stenus, Bembidion octomaculatum) and five aquatic beetles. All of these taxa are recorded together in IFU_1, IFU_2 and IFU_3 but the relative abundances are generally higher in the first. The two groups of taxa in the 'Onthophagus assemblage' (Onthophagus and Bruchinae) are recorded together only in IFU_1 and IFU_3. The specific composition of the 'Pterostichus cursor assemblage' is the most ecologically incoherent since it combines Apioninae taxa which are commonly found on herbaceous vegetation and P. cursor, a halotolerant hygrophilous beetle living on the margins of brackish water. These taxa were sporadically recorded but they were most frequent in IFU_3. The 'Hypoborus ficus, Bagous and Berosus assemblages' comprise a single taxon. Except for H. ficus (which is almost-exclusively recorded in IFU_2), Bagous and Berosus are more or less continuously recorded in every IFU.

History of the coastal wetland inferred from the pollen and beetle fossil assemblages

A total of 100 plant types belonging to 64 families and 14 NPP were identified via pollen analysis (see Online Resource 5). The most abundant pollen types of aquatic and hygrophilous plants were Sparganium/Typha, Myriophyllum type spicatum, M. type alterniflorum, Radiola linearis type, Isoetes and Cyperaceae. The description and interpretation of the coleopteran, pollen and NPP assemblages, presented in Fig. 3a, are described for each IFU. During this period, the relative abundances and MNI of aquatic and hygrophilous beetles decrease (respectively 24 and 12 taxa) while diversity remains high. The persistence of wetland taxa of the 'B. octomaculatum assemblage' indicates the persistence of freshwater conditions (Fig. 3a). Aquatic beetle diversity briefly decreases (disappearance and reduction in number of taxa belonging to the 'Hydrochus and Ochthebius assemblages') possibly caused by eutrophication of the pond, as this is suggested by an up to 13% increase in the salt-intolerant macrophyte Myriophyllum spicatum.

IFU_3 (5311-3090 cal. yr B.P)

This unit is marked by an increase in microforaminiferal linings concentrations and regular occurrences of sand dune beetles, suggesting that the mire is subject to an increasing marine influence (Fig. 3a). Beginning around 4400 cal. yr B.P., the increase in groups of halotolerant beetles (e.g. Pterostichus cursor and specimens of the genus Ochthebius) and Chenopodiaceae (which probably are halophilous species) also indicates that salinity undoubtedly increased locally (Fig. 3a). Indeed, among the 8 Chenopods currently present on Cavallo Island [START_REF] Médail | Contributions à l'inventaire de la biodiversité terrestre de l'île de Cavallo[END_REF], the majority of taxa are halophytes (i.e. Atriplex halimus, Atriplex prostrata, Beta vulgaris subsp. maritime, Halimione portulacoides, Salicornia patula, Salsola kali subsp. tragus), the others being halo-resistant (Chenopodium album) or non-native species recently introduced to the island (Chenopodium ambrosioides).

The diversification of wetland beetle fauna (34 taxa) which occurs between 4600 and 4200 cal. yr B.P., simultaneously to the development of Sparganium-Typha (up to 66%) is followed by a collapse of these wetland beetle fauna (e.g. 'B. octomaculatum assemblage') from 3700 cal. yr B.P. coinciding with the maximal concentration of microforaminiferal linings. During this period, the pond is under a strong marine influence and salt water intrusions are probably occurring.

IFU_4 (3090-345 cal. yr B.P.)

This last unit is marked by the impoverishment of all ecological groups (38 taxa identified in the entire insect fauna unit). From 3100 cal. yr B.P., the increase of Tamarix frequencies (up to 23%) is concomitant with the decrease of microforaminiferal linings (Fig. 3a), suggesting that the formation of a sand-barrier might have disconnected the pond from marine input. This sand-barrier system (favourable for Tamarix development) probably causes a modification of water circulation and creates mixed water sources beginning at 2000 cal. yr B.P., which allows for the persistence of ubiquist (e.g., Berosus, Bagous) or halophilous aquatic beetles of the genus Ochthebius and the development of hygrophilous plants (Cyperaceae). At that time, the morphology of the Greco pond was probably comparable to its current state.

Past changes to the relative sea-level

Figure 3b is a compilation of data on RSLR available in the literature for the north-west Mediterranean basin and modified according to the last study by [START_REF] Vacchi | Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: sea-level variability and improvements in the definition of the isostatic signal[END_REF]. First, during the period encompassed by IFU_1, the relative sea-level is estimated as rising from -6.4 ± 0.5 to -4.7 ± 0.5 m.a.s.l. between ~7200 and 6200 cal. yr B.P. This level is reconstructed at -2.8 ± 0.5 m.a.s.l. at the end of IFU_2 (Vella and Provansal 2000) (Fig. 3b). Whereas the sea-level data for the Rhone Delta show a plateau between ~5000 and 4000 cal. yr B.P. [START_REF] Vella | Relative sea-level rise and neotectonic events during the last 6500 yr on the southern eastern Rhône delta, France[END_REF], data from northern Corsica and La Ciotat indicate respectively a continuous rise until -1.6 ± 0.4 and -1.3 ± 0.4 m.a.s.l. at ~4000 cal. yr B.P. This suggests a rapid rise of about 5.1 m during IFU_1, IFU_2 and the first part of IFU_3. Finally, during IFU_4, the sea-level estimated from northern Corsica and southern France coasts indicates a rise of -1.0 ± 0.4 over the last ~3100 cal. yr B.P. [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF][START_REF] Morhange | Changes of relative sea level during the past 5000 years in the ancient harbor of Marseilles, Southern France[END_REF]. The pattern of RSLR indicates a significant decrease in rates of rise for the second part of IFU_3 and all of IFU_4. The last slowing of the marine transgression during IFU_4 is concomitant with the formation of a sand-barrier and the expansion of Tamarix (Figs. 3a,b).

Discussion

Beetle assemblages as bioindicators of salinity changes

In reconstructing sea-level variations along coastal marine environments, previous studies have focused on some bioindicator assemblages (e.g., foraminifera, ostracods, diatoms, molluscs). Our results suggest that wetland Coleoptera could represent excellent bioindicators of salinity changes because of their narrow ecological requirements. For example, in the North Sea coastal plain, the increase in number of halotolerant beetles is interpreted as being the result of increasing salinity from marine incursions [START_REF] Ponel | Insect evidence for environmental and climate changes from Younger Dryas to Sub-Boreal in a river floodplain at St-Momelin (St-Omer basin, northern France), Coleoptera and Trichoptera[END_REF]). On Cavallo Island, the record of three dulcicolous aquatic beetles from the 'Hydrochus assemblage' and the presence of salt-intolerant plants (Myriophyllum type spicatum) indicate that freshwater conditions prevailed until ~5500 cal. yr B.P. Following this, the composition of wetland beetle fauna suggests that until 3700 cal. yr B.P. the salinity rise caused by marine transgression may have resulted in the disappearance of several aquatic beetles (out of 35 fossil aquatic taxa, 94% of them are present in IFU_1 and IFU_2 and 54% are in IFU_3), whereas hygrophilous beetles were less impacted (out of a total of 32 hygrophilous beetles, 63% of them are present in IFU_1 and IFU_2 and 59% are in IFU_3) (Fig. 3a). Saltmarsh and salt land beetle diversity remains stable but the relative abundance of some halophilous taxa (Pterostichus cursor, Ochthebius) increases until 3700 cal. yr B.P., suggesting that this well-adapted fauna was able to find suitable salinity conditions (Fig. 2).

◄ Fig. 3 a Selected ecological groups of beetles (underlined letters, in minimal number of individuals: MNI), pollen (percentage of the pollen calculation sum) and NPP (concentration of elements by grams) from Greco pond. Taxa included into ecological groups are given in Online Resources 4 and 6. b Relative sealevel curve in the NW Mediterranean Sea from observations in SE France [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF][START_REF] Morhange | Changes of relative sea level during the past 5000 years in the ancient harbor of Marseilles, Southern France[END_REF][START_REF] Vella | Relative sea-level rise and neotectonic events during the last 6500 yr on the southern eastern Rhône delta, France[END_REF], Corsica (Laborel et al. 1994) and from eustatic and glacio-hydroisostatic model (Lambeck and Prucell 2005) The Holocene relative sea-level rise as a driver of coastal wetland biodiversity?

The succession observed in the beetle assemblages fits well with the reconstruction of past RSLR based on several studies (Fig. 3b): during the period covered by IFU_1 and IFU_2, freshwater indicators are present while the sea-level is remains below -2.8 ± 0.5 m.a.s.l. [START_REF] Vella | Relative sea-level rise and neotectonic events during the last 6500 yr on the southern eastern Rhône delta, France[END_REF]. When the marine influence is at its maximum in IFU_3, relative sea-level is estimated at -1.5 ± 0.3 m.a.s.l. in southern France [START_REF] Morhange | Changes of relative sea level during the past 5000 years in the ancient harbor of Marseilles, Southern France[END_REF], corrected by Vacchi et al. 2016) and estimates vary between -1.3 ± 0.4 and -1.1 ± 0.4 m.a.s.l., at respectively ~3800 and 3500 cal. yr B.P. in northern Corsica [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF], corrected by Vacchi et al. 2016) (see Fig. 3b). After 3700 cal. yr B.P., only 31% of all aquatic taxa, 38% of all wetland taxa and 58% of all halotolerant taxa remain. This indicates that the wetland insect richness and abundance decline greatly when freshwater habitats are impacted by a sharp increase in salinity due to a shift of the coastline. However, the potential impact of taphonomy on the past biodiversity record at Cavallo remains questionable. It must be noted that the drastic lowering of total MNI and Coleoptera specific richness curves is concomitant with sedimentological and fossil evidence of marine influences (Fig. 3a). Therefore, changes from freshwater to brackish water could have impacted the deposition and preservation processes. Moreover, the slight decrease in sedimentation rate from 262 cm onward might have led to a degradation of beetle macroremains prior to their burial in the sediment (see Online Resource 2). If this hypothesis were valid, all beetle assemblages would have been poorly preserved, which is not the case. As examples, occurrences of Bagous, Berosus or Carpelimus are more-or-less ubiquitous throughout (Fig. 2). There is also a nearly continuous presence of Ochthebius along IFU_4. Moreover, the fossils of this last insect unit do not show the visible marks of degradation commonly found in cases where there are taphonomical problems (i.e., holes, corrosion marks and thinning down of macroremains). Pollen analysis clearly revealed local aquatic vegetation changes and salinity increases in IFU_3 and IFU_4 (Fig. 3a). This rather suggests that the diversity decrease of the Greco pond is a consequence of adverse conditions for fresh water aquatic and hygrophilous beetles.

Schallenberg et al. ( 2003) also found reduced diversity and altered structure in invertebrate communities in coastal freshwater wetlands following marine intrusions. Other studies at shorter temporal scales generally show an inverse relationship between the increase of salinity and taxa richness of macroinvertebrates in lakes [START_REF] Hammer | Saline Lake Ecosystems of the World[END_REF]), wetlands [START_REF] Pinder | Occurrence of aquatic invertebrates of the wheat belt region of western Australia in relation to salinity[END_REF] or rivers [START_REF] Velasco | Response of biotic communities to salinity changes in a Mediterranean hypersaline stream[END_REF]. [START_REF] Pinder | Occurrence of aquatic invertebrates of the wheat belt region of western Australia in relation to salinity[END_REF] registered a decline in total wetland species richness as salinity increased (above 4.1 g l -1 ) but when salt-tolerant species were excluded, a lower level of salinity (from 2.6 g l -1 ) led to a decline in species richness. This pattern is comparable to our results, since fresh water beetles of the 'Hydrochus assemblage' first decline in IFU_2 even though the taxa richness remains relatively high until 3700 cal. yr B.P. (Figs. 2 and3a).

It must be also noted that another water body was previously located at the center of Cavallo: this swamp has been described in the IGN map dating from 1866 as covering around 1.6% of the island surface (Fig. 1b). There is no doubt that such a large wetland has played a major role as 'refuge' for aquatic and wetland fauna when environmental conditions became unfavourable in the coastal pond.

Did climate or humans play a role in the observed biodiversity changes?

The important question remains as to whether other drivers, such as human alteration or climatic changes, in addition to Holocene RSLR have influenced beetle diversity. In the central Mediterranean, high lake levels from Italy [START_REF] Magny | Possible complexity of the climatic event around 4300-3800 cal. BP in the central and western Mediterranean[END_REF][START_REF] Magny | North-south palaeohydrological contrasts in the central Mediterranean during the Holocene: tentative synthesis and working hypotheses[END_REF], low fire activity in the western Mediterranean [START_REF] Vannière | Circum-Mediterranean fire activity and climate changes during the mid-Holocene environmental transition (8500-2500 cal. BP)[END_REF]) and pollen-inferred climate reconstructions across the Italian Peninsula [START_REF] Peyron | Contrasting patterns of climatic changes during the Holocene across the Italian Peninsula reconstructed from pollen data[END_REF] suggest prolonged wet phase from ~4300 to 3800 cal. yr B.P. However, even though this episode could perhaps explain the transitory diversification of wetland beetles and the scarce presence of aquatic freshwater plants around ca. 4400 cal. B.P., climatic factors apparently played a minor role. Notwithstanding the salinisation process from ~4100 cal. yr B.P. which could have possibly been buffered by wetter conditions, the impoverishment of the wetland beetle fauna after 3700 cal. yr B.P. seems mainly due to the RSLR because respectively 69 and 63% of aquatic and hygrophilous fossil beetles are only recorded at times prior to marine intrusions.

On the other hand, it is well known that human activities exacerbated by sea-level changes have triggered coastal instability [START_REF] Anthony | Human influence and the changing geomorphology of Mediterranean deltas and coasts over the last 6000 years: from progradation to destruction phase?[END_REF]). In the Mediterranean region, human impacts led to major environmental changes that are reflected in the coleopteran fauna dating back to Neolithic times (Andrieu-Ponel and Ponel 1999) and during the Antiquity (Andrieu- [START_REF] Andrieu-Ponel | Palaeoenvironments and cultural landscapes of the last 2000 years reconstructed from pollen and Coleopteran records in the Lower Rhône Valley, southern France[END_REF]. On Cavallo Island, significant records of open environment and coprophilous beetles, ruderal-anthropogenic plants and coprophilous fungi suggest successive phases of human activities during the Neolithic period (from 6993 to 6200 cal. yr B.P.), the Chalcolithic period (from 5000 to 4000 cal. yr B.P.) and since the Iron Age (from 2700 cal. yr B.P. onward). This is confirmed by several traces of human occupation in the Lavezzi archipelago since the Chalcolithic time according to [START_REF] Agostini | Recherches archéologiques dans l'île Cavallu (Bonifacio, Corse) 1972-1977[END_REF]. The opening of the vegetation (increase of herbaceous pollen) and the decline of saproxylic beetles which are clearly recorded during the second anthropogenic phase from 5000 to 4000 cal. yr B.P. suggest that this long-time human occupation has impacted the landscape. Nevertheless, our results suggest that there is no negative impact on wetland and aquatic Coleoptera, since a marked abundance increase in these communities is recorded around ca. 4400 cal. B.P. (Fig. 3a). It is likely that the decline in trees and shrubs induced by humans should have caused an increase in soil erosion and the progressive filling in of the pond, leading to the reed bed (Sparganium-Typha) expansion on the pond margin along with the associated hygrophilous beetle fauna. Finally, the most 'deleterious' human impact has certainly been the recent urbanisation of Cavallo since the 1970s (see Fig. 1c). The filling in of the main swamp during the construction of an airstrip is probably partially involved in the differences between current and past biodiversity, but deconstructing its specific impact is beyond the scope of this paper. Moreover, this urbanisation predated all modern entomological inventories [START_REF] Cocquempot | Liste des Arthropodes terrestres recensés dans les réserves naturelles des îles Cerbicale et Lavezzi (France, Corse-du-Sud)[END_REF][START_REF] Médail | Contributions à l'inventaire de la biodiversité terrestre de l'île de Cavallo[END_REF][START_REF] Médail | Contributions à l'inventaire des arthropodes terrestres et aquatiques sur l'île de Cavallo (Archipel Lavezzi, Bonifacio, Corse du Sud). Note naturaliste, Initiative pour les Petites Îles de Méditerranée (PIM)[END_REF].

Past local extinction of some remarkable beetles

Currently 283 beetles species are known to inhabit the Lavezzi archipelago [START_REF] Cocquempot | Liste des Arthropodes terrestres recensés dans les réserves naturelles des îles Cerbicale et Lavezzi (France, Corse-du-Sud)[END_REF][START_REF] Médail | Contributions à l'inventaire de la biodiversité terrestre de l'île de Cavallo[END_REF][START_REF] Médail | Contributions à l'inventaire des arthropodes terrestres et aquatiques sur l'île de Cavallo (Archipel Lavezzi, Bonifacio, Corse du Sud). Note naturaliste, Initiative pour les Petites Îles de Méditerranée (PIM)[END_REF]. The comparison of palaeoecological data and modern inventories shows that 63 of them were found in the fossil assemblages whereas 1 family, 1 subfamily, 31 genera and 37 species of fossil beetles were not reported. Although this could be tempered by the fact that the modern Coleoptera fauna of Cavallo is probably incompletely inventoried, 60% of the hitherto unknown taxa are aquatic and hygrophilous beetles (respectively 24 and 18 taxa on a total of 70). Both are relatively easy to detect, since the Greco pond is the only extant body of water on the island. This suggests that many beetles became locally extinct during a short time span (less than 7000 years).

Even though RSLR has directly affected the aquatic and hygrophilous coleopteran diversity through a major change in salinity, it has also induced an indirect impoverishment of beetle fauna through the wetland vegetation composition changes. Indeed, many authors have shown a close relationship between macrophytes communities and the diversity, abundance and distribution patterns of invertebrates [START_REF] Ali | Importance of aquatic macrophyte for invertebrate diversity in large subtropical reservoir[END_REF][START_REF] Walker | Macrophyte presence and growth form influence macro invertebrate community structure[END_REF]) since they play a major role in maintaining the freshwater ecosystems by determining the physical structure of these systems (Chambers et al. 2008), by increasing habitat complexity (Thomaz and Ribiero da Cunha 2010), and by constituting a food resource for invertebrate populations [START_REF] Guignot | Effect of macrophytes on Phytophilous macroinvertebrate community: a review[END_REF]. Thus, the disappearance of salt-intolerant plants (i.e. Alisma-type, Lythrum type junceum, L. type salicaria, Myriophyllum type spicatum, Radiola linearis-type) has been tantamount to a loss of habitat for stenotopic species (Fig. 3a). Such is the case of Hylobius transversovittatus and Phytobius leucogaster which respectively live on Lythrum salicaria and Myriophyllum sp. [START_REF] Koch | Biological evidence of sea-level rise during the last 4500 years on the rocky coasts of continental southwestern France and Corsica[END_REF] that were only observed between 6500 and 4500 cal. yr B.P.

The Holocene RSLR cannot by itself obviously explain the apparent extinction of hygrophilous taxa of the reed bed such as Protapalochrus flavolimbatus and Oodes gracilis (Koch 1989a;Saint Claire Deville 1914) since their habitat is still present on the Greco pond. The recent human impact on the Island, which has caused particular damage to Cavallo's central swamp, has undoubtedly played a major role in the local extinction of these stenotopic species. As already evidenced in southern France, drainage and expansion of cultivation might have contributed to regional extinction of Odacantha melanura and Cerapheles lateplagiatus, two species that occur in very similar biotopes [START_REF] Andrieu-Ponel | Human impact on Mediterranean wetland Coleoptera: an historical perspective at Tourves (Var, France)[END_REF]. As a whole, the fossil occurrences of all these now-absent wetland species suggest that their past regional distributions were greater than they are today, as is true of their biotopes.

Conclusion

Our multidisciplinary investigation has revealed that the past ecosystem dynamics and biodiversity changes in the coastal wetland of Cavallo Island were significantly driven by 7000 years of RSLR. Coleoptera fossil assemblages responded rapidly to increasing salinity, culminating in the local extinction of freshwater aquatic beetles and the regression of hygrophilous and riparian freshwater beetles. The regression/disappearance of fresh-water helophytes and the subsequent increase in salty ground plants have played also a major role in these insect diversity changes. Furthermore, the recent destruction of the central swamp on Cavallo, which may have played a role as a reservoir for the freshwater insect community, has undoubtedly worsened the resilience of the island ecosystems that have faced past RSLR.

It can therefore be concluded that the effects of ongoing RLSR on coastal diversity could be more severe in areas destabilised by humans: because the Mediterranean coastline is experiencing increasing human pressure on natural habitats (due to urbanisation, tourism, concretisation of the littoral; see [START_REF] Beltrame | Mediterranean wetlands: outlook. First Mediterranean wetlands observatory report-Technical report-2012[END_REF], the protection of coastal freshwater bodies and inland marshes (especially on small continental islands; see Médail 2013; PIM initiative: Initiative pour les petites îles de Méditerranée) represents a high priority. Such a habitat loss has already resulted in a drastic retraction of the geographical range of many aquatic beetles over the past hundred years (Ribera 2000). Thus, a Mediterranean sea-level rise of 9.8-25.6 cm by 2050, as estimated by [START_REF] Galassi | Sea-level rise in the Mediterranean Sea by 2050: roles of terrestrial ice melt, steric effects and glacial isostatic adjustment[END_REF], represents an added threat on top of generalised direct human pressure: coastal freshwater wetlands will probably suffer salinity increases while inland water bodies, which for now are distant from the coastline, will no longer play the role of refugia.

Facing future RLSR, comparison between prospective studies such as [START_REF] Irmler | Zonation of ground beetles (Coleoptera: Carabidae) and spiders (Araneida) in salt marshes at the north and the Baltic sea and the impact of the predicted sea level increase[END_REF] and retrospective studies such as the present work will be valuable in improving estimates of how coastal wetland insect diversity will be likely to change their composition and adapt. This is going to be an important field of research in the Mediterranean basin where some 12,500 islands of all sizes have been identified, most of them biodiversity hotspots [START_REF] Médail | Hotspots revisited: earth's biologically richest and most endangered terrestrial eco regions[END_REF][START_REF] Médail | The unique nature of Mediterranean island floras and the future of plant conservation. In: Cardona Pons E, Estaún Clarisό I[END_REF]. The future trajectories of these island ecosystems remain questionable, but there is no doubt that their capability to be resilient to RSLR depends on the persistence of natural habitats beyond the coastline. The major challenge for decision-makers will be to preserve and protect these future refugia (for instance, by reducing damming, draining, fragmentation, pollution and isolation of source habitats).

Fig. 1 a

 1 Fig. 1 a Location map of the Lavezzi archipelago (Southern Corsica, France), showing Cavallo Island; b IGN map from 1866 (http://www.geoportail.gouv.fr, SARTENE); c photo of Cavallo Island, taken in 2002 (http://www.geoportail.gouv.fr, aerial photography mission CP02000032_2002_fd2A2B250_0223, picture n°223)

Table 1

 1 Groups of Coleoptera according to their ecological affinities

	Ecological groups	Ecological features and non-exhaustive beetle examples
	Aquatic Coleoptera	Beetles living in standing and/or running water at some stage of their
		development (e.g., taxa belonging to the Dytiscidae or Hydrophilidae
		families) that can be halotolerant or not
	Salt marsh/Salt land Coleoptera	Halotolerant aquatic beetles (e.g., Hydrovatus cuspidatus) or non-
		aquatic beetles living on salty land (e.g., Bledius unicornis)
	Woodland and scrubland	Beetles living in forest or shrubby environments (e.g., Scobicia
	Coleoptera	chevrieri, Triodonta cf. cribellata) that can be xylophagous or not
	Sand dune Coleoptera	Psammophilous beetles adapted to sand biotope (e.g., Catomus
		consentaneus, Platytomus laevistriatus)
	Driftwood-dependent Coleoptera Beetles strictly associated with driftwood beached along the coastline
		(e.g., Mesites pallidipennis, Pselactus spadix)
	Wetland and waterside Coleopte-	Non-aquatic beetles associated with wetlands and marshes (e.g., Dryops,
	ra	Phalacrus). A few of these beetles were also halotolerant
	Open environment Coleoptera	Beetles living in grasslands (e.g., Bruchidius, Longitarsus)
	Decayed vegetable matter-	Non-coprophilous beetles feeding or living in decayed vegetal matter
	dependent Coleoptera	(not presented on the Fig. 3a)
	Coprophilous/coprophagous	Beetles feeding directly upon dung (e.g., Caccobius schreberi,
	Coleoptera	Onthophagus taurus) or predatory, feeding upon small coprophagous
		organisms such as Diptera larvae (e.g., Oxytelus piceus). A few of
		these beetles were classed both in this group and that of open
		environment or woodland and scrubland Coleoptera
	No information

Table 2

 2 Coleopteran assemblages obtained by the unconstrained cluster analysis performed on Taxa-matrix (see Fig.2) Beetle diversity is maximal (up to 32 taxa per sample). Aquatic and hygrophilous beetles prevail (respectively 27 and 15 taxa) and some taxa of the 'Ochthebius, B. octomaculatum and Hydrochus assemblages' reach their highest relative abundances (Fig.2). Wetland environment is indicated by freshwater plants (e.g. Myriophyllum type spicatum, Radiola type linoides) and dulcicolous beetles of the 'Hydrochus assemblage'(Figs. 2, 3a) suggesting that no marine influence has significantly disturbed the pond ecosystem during this period although salt marsh and salt land coleopteran are also well-represented in this IFU. These beetles probably indicate the presence of salty environments on Cavallo Island not far from the freshwater pond.

	Coleoptera assemblage	Taxa included and ecological features
	Hypoborus ficus	Hypoborus ficus (saproxylic beetle digging galleries in the branches of Ficus
	assemblage	carica)
	Ochthebius assemblage	Ochthebius (aquatic genus where the majority of species are halophilous);
		Contacyphon (subaquatic); Carpelimus (wetland margin beetle)
	Bembidion octomaculatum Dryops; Stenus; Bembidion octomaculatum (3 dump ground beetles);
	assemblage	Limnoxenus niger; Hygrotus; Enochrus (3 aquatic halotolerants);
		Colymbetes fuscus; Noterus (2 ubiquitous aquatic beetles)
	Bagous assemblage	Bagous (aquatic beetle of standing or slow-moving waters)
	Berosus assemblage	Berosus (aquatic beetle of standing waters)
		hygrophilous
	assemblage	beetle)

Onthophagus assemblage Onthophagus (coprophilous beetle); Bruchinae (open ground beetle) Hydrochus assemblage Hydrochus; Helochares; Limnebius (3 dulcicolous aquatic beetles); Scopaeus (hygrophilous beetle); Mesites (drift-wood dependent beetle) Pterostichus cursor Apioninae (open-land beetle); Pterostichus cursor (halotolerant

Acknowledgements

This study was incorporated within the framework of the CoP 2 A (Corsican Palaeoclimate, Palaeoenvironments & Anthropization) project, which is supported by the Laboratoire d'Excellence Objectif-Terre Bassin méditerranéen (Labex OT-Med) of Aix Marseille University (ANR-11-LABEX-0061), and by the DyPaCo (Dynamique des paléoenvironnements de la Corse, convention no.15/005) project which is supported by the Office de l'environnement de la Corse (OEC)/Conservatoire botanique national de Corse. The fieldwork was funded by the French government (projet Investissements d'Avenir) within the Initiative d'excellence A*MIDEX/MEDNET of Aix Marseille University (ANR-11-IDEX-0001-02), and performed during a field school of the Master SET-SBEM of Aix Marseille University.

The authors wish to thank the Conservatoire Botanique National de Corse and its director Laetitia Hugot for their constant support of this research. We also thank the Association pour la protection de l'environnement de l'île de Cavallo (APEIC), notably Mrs. Matthieu Bidali and Michel Orlanducci for their permission to work on the private island of Cavallo. Finally, thanks are due to the two anonymous reviewers for their constructive remarks.