E Brun

J Vicente

F Topin

R Occelli

M J Clifton

Microstructure and transport properties of cellular materials: Representative Volume Element

The representative volume element (Rve) plays a central role in efforts to predict the effective thermo-physical and transport properties of heterogeneous materials. A quantitative definition of its size is proposed in this work. It is shown that Rve depends on the morphological or physical property being investigated. The methodology is applied to real samples of open-celled materials (such as metallic foam) whose structure is obtained from X-Ray microtomography. Geometrical Rve measurements have been carried out on a wide range of open-cell structures from 10 ppi up to 40 ppi from different manufacturers. Finally, physical Rve measurements for fluid-flow laws are derived from numerical simulations on representative metallic foam. Finite Volume and Lattice Boltzmann simulations are compared and validated by experimental data.

Introduction

Metallic foams are a class of materials that are attractive for numerous applications, as they present high porosity, low relative density, high thermal conductivity of the solid phase and large accessible surface area per unit volume [START_REF] Ashby | Metal foams: A design guide[END_REF]. Moreover, they also promote mixing and have excellent specific mechanical properties. Metallic foams are thus used in the field of compact heat exchangers, reformers, two-phase cooling systems and spreaders [START_REF] Banhart | [END_REF]3]. Foams have also been used in high-power compact batteries and catalytic-reactor applications such as fuel cell systems [START_REF] Catillon | presented at 19th North American Catalysis Society Meeting[END_REF].

Because of their relative novelty and their complex three-dimensional structure, the transport properties of metal foams have not been exhaustively characterized in contrast with mechanical properties. Classical models widely used for low porosity media are difficult to apply to high porosity materials [START_REF] Vicente | IHTC13, 13TH International Heat Transfer Conference[END_REF].

To model cellular materials, most published work uses idealized and periodic structures (e.g. [START_REF] Boomsma | [END_REF]). Local changes in the structure (e.g. in constriction, strut cross section, surface roughness…) can influence the properties [7,8]. The problem is not only "how to model a cell?" but "how many cells are necessary to represent the macroscopic property studied?"

We propose in this work to measure and compare representative volume elements for geometrical and physical properties. The first part of this article presents the set of open-celled samples. The second part is devoted to the numerical simulation methods used in identifying the size of Rve associated to fluid-flow parameters. The Rve estimation procedure is presented in the third part.

Finally we compare and discuss the geometrical and physical Rve.

Materials

We present results on a wide-ranging set of commercially available metallic foams of different manufacturers and metals. We used cylindrical samples of 40 mm and 16 mm diameter, 13 mm and 6 mm thickness (see Table 1). High-resolution microtomography acquisition was performed on the ID19 beam line of the ESRF. Sample size was maximized while respecting the datavolume constraint and X-ray image resolution. We thus acquired a large number of cells while efficiently detecting hollow struts [8]. In previous studies, we measured useful morphological parameters such as specific surface, porosity, pore diameter distribution and cells orientations [8,[START_REF] Brun | Geometrical measurement of real foams from 3D images[END_REF]. Some of these parameters are reported in Table 1.

It appears clearly that these foams present a rather low specific surface, a large pore diameter and a very high porosity in comparison to classical porous media. We have demonstrated [START_REF] Vicente | IHTC13, 13TH International Heat Transfer Conference[END_REF] that the cells are spatially organized (common mean orientation) but some defects appear in this arrangement. Cells of different orientation and size materialize these defects. They appear to accommodate topological constraints associated to plateau border angles, space filling and constraints produced by manufacturing processes. Thus, Rve should include at least some of these defects and their characteristic size should be comparable to defects distance.

Numerical simulations

We present in this part the numerical methods used to simulate fluid flow in the real open-celled structure obtained from X-ray images. The sample used is the Recemat NC1723 (see Fig. 1 ).

As the Lattice Boltzman method uses periodical velocity conditions (inlet and outlet) the 3D reconstructed foam was double mirrored in the direction of flow. The parallelepipedic sample comprises 198x198x342 cubic elements. To compare results of Lattice Boltzmann simulation and Finite Volume simulation, the same-mirrored sample has been used for the Finite Volume method. To identify the Forchheimer fluid flow parameters, simulations have been carried for 5 different Reynolds numbers ranging from approximately 55 up to 157.

Finite volume method (FV)

Navier Stokes equations are discretized with the finite volume method of Patankar [START_REF] Patankar | Numerical Heat transfer and Fluid Flow[END_REF] with a centered velocity mesh. We use the SIMPLE Algorithm modified by Date [START_REF] Date | [END_REF] with smoothing pressure correction. The advection terms are discretized with a first-order UPWIND scheme.

Boundary conditions are imposed velocity at the inlet and symmetry conditions on the lateral faces. An outlet tranquilization zone of 8 voxels width has been added and classical outflow condition (null velocity gradient) has been used. A null velocity is imposed on the solid matrix voxels. The computation is carried on a cluster of 16 AMD 2.3 GHz with one GB of RAM for each processor. The CPU time for a simulation is about 21hours (steady state computation).

Lattice Boltzmann Simulation

The lattice Boltzmann method (LBM) is another numerical method to simulate viscous fluid flow. The LBM approximates the continuous Boltzmann equation by discretizing physical space with lattice nodes and velocity space by a set of microscopic velocity vectors. Space is discretized into uniform cubic voxels that are either void or solid with a lattice node at the center of each. The time-and space-averaged microscopic molecular movements are modeled using molecular populations called distribution functions, which define the density and velocity at each lattice node. Specific particle interaction rules are set so that the Navier-Stokes equations for fluid flow are satisfied. The flow through the metal foam was calculated using the Lattice Boltzmann method in the Multiple-Relaxation-Time version [12] with a D3Q19 lattice and a numerical viscosity of 10 -3 . As the image was composed of cubic voxels, the conditions imposed at the solid/liquid interfaces were simple bounce-back. The boundary conditions at the inlet and outlet were periodic in velocity with an imposed pressure difference; symmetry conditions were applied on the sides. The calculation was transient with a zero velocity initial condition. Once the flow field had become stable enough (i.e. constant macroscopic averaged velocity), images of the flow field were recorded regularly. The CPU time for a simulation is about 60 hours on a single AMD 2GHz processor.

Results

Cells are organized in successive layers, perpendicular to the main flow direction. We can observe on figure 2 a roughly constant pressure in the pores of each layer. The pressure drops from layer to layer because of cell constrictions and structural organization. The velocity increase at the cell throat and current line separate and merge around the struts. A small anisotropic effect could be inferred as isobars are not fully perpendicular to the main flow direction.

In the range of Reynolds number studied the local flow was always non-stationary. A swirling motion similar to a von Karman vortex street (successive eddies formation and lift off) was observed in the wake of the struts. This behavior explains the importance of inertia effects and the non-darcian behavior of the flow.

Validation with experimental fluid flow laws

In these media, fluid flow is described by the Forchheimer law which relates the pressure drop to the mean superficial velocity in the medium (eq. 1.) and leads to characterization in terms of the permeability K and the inertia coefficient b [7].

eq. 1

To establish these parameters numerically, simulations must be performed for several values of Reynolds number. K and b are then computed by fitting pressure gradient versus superficial velocity data with the Forchheimer model.

A previous study [7] indicates that published experimental fluid-flow parameters for open-celled materials show a wide dispersion. For instance, permeability values for the same medium can differ by up to three orders of magnitude. The final uncertainty over the experimental data presented here is about 50%. The experimental values were determined over a much wider range of Reynolds number (from 10 to 10 000) than the numerical one. Values of calculated inertia coefficient (Table 2) present a small uncertainty and are quite comparable. Calculated permeability values differ from factor two. Globally, LBM results stay closer to experimental data than FVM one (figure 3) but both results are well fitted by second order polynomials. This confirms that the Forchheimer law describes the fluid flow behavior in these media. In the range of Reynolds numbers studied, LBM fits the experimental data better than FVM. The most probable explanation of the discrepancy between the two numerical methods is that FVM results are obtained with a too coarse mesh to capture efficiently velocity gradients end eddies in the vicinity of the struts. Moreover, LBM results are obtained by averaging transient calculations over time while FVM results are calculated at steady state. FVM calculations are done in the stationary regime while some non-stationary effects are present (as shown by transient results) and are not well captured. For Navier Stokes computation in such complex geometry and at relatively low Reynolds number, clearly LBM is more efficient than FVM as it gives more precise results and uses far less CPU time for a given structured mesh. A mesh dependency study is under investigation. We'll use an octree decomposition mesh that will allow a finer resolution near the solid interface in order to efficiently capture local boundary layers and velocity gradients.

Rve determination procedure

A definition of the Rve was recently proposed by Drugan and Willis [13]''It is the smallest material volume element of the composite for which the usual spatially constant (overall modulus) macroscopic constitutive representation is a sufficiently accurate model to represent mean constitutive response''. We propose to qualitatively identify Rve associated to both geometrical and physical macroscopic properties using a statistical approach [14].

The property considered is calculated for a large number of fixed-size boxes (n=100 to 1000) randomly located in the sample. We obtain a distribution of mean values for a given size of box.

We then study the influence of the box size on the standard deviation of these distributions. The size of the Rve is defined as the volume for which the distribution of the mean values presents a standard deviation below a given threshold. We could also define the size of the Rve as the volume for which the maximum difference is below a given threshold; it would be a much robust criteria since Std statistics could hide a MD far greater than Std Dev. In our case both quantities vary similarly and thus STD Dev threshold is sufficient.

For all studied quantities, the standard deviation decrease roughly exponentially with the box size. For a Std Dev of 10 % we observe that maximum differences (MD) between two realizations could exceed 25%. At 5% Std Dev, MD stay around 15%. This value is indeed rather important but the available samples do not permit to be more precise.

Concerning velocity, the method for obtaining Rve is the same as for geometrical properties. But determining the local averaged pressure gradient from the pressure field is more difficult [15,16]. Even if local pressure gradient are not fully perpendicular to main flow direction (Figure 2), this effect is small and we consider that macroscopic pressure is constant on planes perpendicular to the main flow direction (established flow in an isotropic medium). Thus, we calculate mean pressure on planes perpendicular to the main flow direction at several locations, and we estimate the slope of the mean pressure variation by linear regression. We have chosen the coefficient of determination R 2 of the linear regression as a pertinent parameter and we keep our 5% criteria to establish the Rve size. We made these calculations for different sizes of a randomly located box in the medium and we studied the influence of the box volume on the coefficient of determination R 2 of the linear regression.

Rve estimation and discussion

To compare the different foams, the size of the box was normalized by the pore diameter.

Moreover the standard deviation and the maximum difference are normalized with the mean value of the property studied for the considered foam. The figure 4 and5 show the evolution of the standard deviation and the maximum difference in function of the box size.

Geometrical Rve

Porosity and specific surface are estimated from the 3D density images obtained with X-ray tomography. Porosity is trivially obtained by counting the voxels of the different phase. The "Marching Cubes" algorithm uses the 3D density images to obtain the mesh of solid surface [17].

The mesh data triangles are ordered so that each triangle is included into a unique cubic mesh composing the 3D sample. Once the mesh has been generated for the sample, it is then easy to calculate a specific surface area within a small volume taken in the sample: the cubes composing the volume are retrieved and the surface area of the triangles contained in these cubes is summed.

Figure 4 shows the influence of the box size on standard deviation and on the maximum difference of both porosity (a) and specific surface (b) area distribution. For porosity we can define a unique Rve for these samples. In effect, we observe similar standard deviation variations for the several manufacturers. According our 5% criteria the Rve of porosity is a cube of 0.7 pore diameter side.

Similar behavior in standard deviation variation is observed for the Rve of specific surface. A single size of Rve can be defined for all these samples: similar standard deviation variations are observed for the foams from different manufacturers. According our 5% criteria, the Rve for specific surface is a cube roughly 3 dp in side length.

The absence of differences between foams from different manufacturer can be explained by the fact that the manufacturing processes are roughly similar (replication of polymer foam) and lead to the same solid volume fraction (porosity around 90 % is a manufacturer goal). Plateau law gives globally struts shape but local variations in the cells organization induce different shapes of struts inter-connection and strut cross section (roughly curved triangle). These perturbations may be due to several effects (gravity, drainage, deposition, casting…) during manufacturing process. This gives insight on the fact that specific surface Rve size is far greater than porosity one, but further investigation needed to clarify this point.

Physical Rve

Figure 5.a shows the influence of the box size on standard deviation and on the maximum difference of mean velocity. To establish Rve for the mean velocity value we used the velocity field obtained from the computations. No clear influence of the Reynolds number on RVE size can be identified for both methods on the studied range. We observe that for LBM simulation a box of 3 dp side seems to be sufficient to define the mean velocity for these Reynolds numbers.

Velocity magnitudes obtained with FVM present bigger fluctuations in standard deviation.

Therefore Rve for FVM is bigger than for LBM and equal to 5 dp. This indicates that pressure fields given by both numerical methods are similar. On the other hand velocity field are quite different, due to the insufficient mesh condition velocity fluctuations calculated by FVM are apparent on larger scale The estimate of the pressure gradient is quite constant for a box 3 pore diameters in size and the coefficient of determination R 2 is good (>0.99). Similar results have been obtained with the FV method. As a result, 3 pore diameters are sufficient to establish pressure drop Rve for the both methods.

The Rve for the fluid flow properties depends on the biggest Rve (velocity and pressure drop), so Rve for LBM is a box of 3 pore-diameter size, while the FV needs a 5 pore-diameter box.

Conclusion

This work presents the first measurement database of Rve (for both geometrical and fluid-flow properties) in metallic foam. Geometrical data were obtained from 3D X-ray CT and full morphological measurements were carried out. Fluid flow simulations at pore level were carried out in real geometries with two different numerical methods which lead to similar results even if some discrepancies are observed, probably linked to the different treatment of transient terms. A reasonable agreement with experimental data is obtained in terms of permeability and inertia coefficient. Both geometrical and physical Rve characteristic sizes are close to pore size. This work will be extended to conduction and conjugate heat transfer as well as other geometrical parameters in order to gain a better understanding of the impact of the structure of the porous The Rve determination is reported in term of the Pore diameter. This pore diameter is a mean value obtained by the analysis of the whole accessible sample. It would be interesting to study the size of the pore diameter Rve. The Rve size of a studied parameter is clearly depending on the dispersion of these parameters.

Finally, it is clear that using a Std Dev threshold of 1% would lead to a more representative evaluation of Rve size. Unfortunately, the usable field of view of our CT data does not permit to use boxes big enough to satisfy this criterion.

Figure 5 .

 5 Figure 5.b presents the estimation of the pressure gradient by a linear regression and its


Samuel Vercraene
email: samuel.vercraene@insa-lyon.fr

Fabien Lehuédé
email: fabien.lehuede@imt-atlantique.fr

Thibaud Monteiro
email: thibaud.monteiro@insa-lyon.fr

Olivier Péton
email: olivier.peton@imt-atlantique.fr

The Dial-A-Ride Problem with School Bell time Adjustment

Mots-clés : vehicle routing, dial-a-ride problem, paratransit, matheuristic

Medical and Social Institutions (MSI) and specialized schools for students with disabilities use specialized transport services every day. As a result, paratransit represents one of the main costs for specialized schools and is also expensive for public authorities. We investigate a practical case from France, where these establishments typically have around 80 % of their users benefiting from a specialized transportation service. In this study, the considered MSI and specialized schools (called schools in the remaining of this abstract) are interested in studying the benefits of pooling their transport. The organization of such transport systems requires solving a Dial-a-Ride Problem (DARP) to design the best possible routes for vehicles that pick up their passengers at home and drop them off at their destination, provided they can share part of their route [START_REF] Tellez | The fleet size and mix dial-aride problem with reconfigurable vehicle capacity[END_REF]. When several schools have similar start times, it is hard to design feasible routes that can pick up passengers living in the same area and going to separate schools. Hence, modifying the schools' start times, which is known as School Bell time Adjustment (SBA) in school bus routing [START_REF] Park | The school bus routing problem : A review[END_REF], facilitates transportation pooling. We address the optimization problem which consists of simultaneously determining the school start times and vehicle routes and denote it the Dial-A-Ride Problem with School Bell time Adjustment (DARP-SBA).

The Dial-A-Ride Problem with School Bell Time Adjustment

We consider a set of schools denoted by E. The set of users is denoted by U. Each user u ∈ U has a transportation request from a pickup node p u to a delivery node d u . This request is associated with a maximal ride time R u and loads q S u and q W u . In the morning case, each school e ∈ E has a morning interval [a e , b e ] during which all deliveries can be scheduled. Deliveries occur within a so-called dynamic time window of width W e that ends at the school start time.

The DARP-SBA consists in determining a start time H e for each school e ∈ E within its morning interval, selecting a set of vehicles and designing the route of each selected vehicle such that : all transportation requests are served ; all maximum ride times, pickup and depot time windows are satisfied ; all users can be dropped off at their school within its time window ; the capacities of the vehicles are satisfied ; and the overall vehicle, traveling and route duration costs are minimized.

Figure 1 gives a graphical representation of a solution with two routes ω 1 and ω 2 visiting a school e. Note that nodes d 1 , d 2 , d 3 , d 5 and d 6 represent the same physical location. One of the key challenges of this problem is that routes should be synchronized at schools, where a dynamic time window is fixed in time when determining the school bell time.

Matheuristic approach and experiments

To solve the DARP-SBA, we propose a matheuristic which combines the following ingredients : -In each iteration, solutions are generated by a large neighborhood search (LNS) derived from [START_REF] Tellez | The fleet size and mix dial-aride problem with reconfigurable vehicle capacity[END_REF]. It the repair operator solves a DARP in which the schools time windows have been filtered. The complete timing problem is solved by a MILP solver only for complete and promising solutions to validate their feasibility. Feasible routes are collected in a route pool. -On a regular basis, a relaxed, restricted DARP-SBA model is solved for this route pool with a time limit. It selects and combines routes generated at different LNS iterations, possibly adjusting school bell times. -Filtering rules and adaptive strategies are used to ensure a faster exploration of neighborhood and the good use of the LNS and MIP components. Our experimental study considers a set of 34 schools in the area of Lyon and the 575 users transported to these schools. Three types of instances (S,M,L) have been produced by splitting the data set in eight, four and two sub-problems, respectively. To evaluate the potential impact of SBA, we have considered four scenarios. In scenario 0, vehicles cannot visit more than one school. Scenarios 20 and 60 correspond to DARP instances where the delivery time windows are [8 :20, 8 :40] and [8 :00, 9 :00], respectively. Finally, scenario 20/60 corresponds to the DARP-SBA with a dynamic time window of width W e = 20 minutes and a morning interval [a e , b e ] = [8 :00, 9 :00] for each school e ∈ E. Table 1 shows the results of all scenarios aggregated by instance types. A first observation is that the greatest savings can be achieved by pooling transports between schools, with around 10% savings on routing costs. Scenario 20/60 shows that adjusting school bell times enables additional cost reduction by 7% on average. The unrealistic scenario 60 allows for only 5.5% of additional savings with respect to SBA. This work was recently published in [START_REF] Vercraene | The Dial-a-Ride Problem with School Bell Time Adjustment[END_REF].

  Partial time space representation of two routes ω 1 and ω 2 visiting a school e.

					a e	H e -W e	H e	b e
	route ω 1	o +	p 1	p 2	p 3	d 1	d 2	d 3	p 4	o -
	route ω 2	o +		p 5	p 6	d 5	d 6		p 7	o -
	FIG. 1 -


Laëtitia Lemière
email: laetitia.lemiere@inrae.fr

Marie Gosme

Gérard Subsol

Marc Jaeger

FSPM applied to agroforestry system co-design

Keywords: agroforestry, FSPM, co-design, augmented reality, graph theory

Context

Agroforestry systems (AFS) are farming systems that associate trees, crops and animals. They cover a variety of management options and spatial arrangement [START_REF] Dupraz | Agroforesterie, des arbres et des cultures Chieze et al[END_REF]. Therefore, co-design workshops are an efficient way to design AFS. To design the spatial arrangement of species, participants often interact with a physical mockup [Chieze et al, 2021], locating plants to define the system structure. But despite the numerous exchanges between participants, the resulting ecosystem functions are difficult to perceive [START_REF] Potschin | Progress in Physical Geography: Earth and Environment[END_REF].

Building models for agroforestry is a growing concern [START_REF] Burgess | [END_REF]. However, agroforestry models seldom express both structural and functional representations. Some, such as Hi-safe [START_REF] Dupraz | Sustainability[END_REF], focus on functional aspects, extending crop model approaches. Others, such as EcoAF [START_REF] Liagre | 4th World congress on Agroforestry[END_REF], focus on the structural descriptions of the system, applying pattern grammars to describe the spatial arrangement of species. These approaches are difficult to mobilize in co-design workshops, facing in the first case the lack of visualization supports and complex parametrization, and, in the second case, the lack of interactions between elements of the modeled system and the scarcity of ecosystem services that are considered.

In the review of [START_REF] Sereno | [END_REF], Augmented Reality (AR) is shown to be efficient for collaborative work. This technology could allow the co-design workshop participants to share a common vision of the AFS and its functioning. Using traditional functionalstructural plant model (FSPM) with AR during design workshop is promising to link scientific knowledge, visualization, interaction and collaboration.

Our proposal

Our proposal is composed of 3 steps: (1) acquiring and digitizing an AFS mockup such as those that are commonly used in co-design workshops, (2) automatic building of a model of the AFS using graph theory [START_REF] Lemière | 5th World Congress on Agroforestry[END_REF] and, (3) visualizing the system and the associated ecosystem services using AR. We can integrate functionalstructural models at two scales: at the plant scale, i.e. FSPM, and at the agroecosystem scale, leading to a new concept of Functional-Structural Agroecosystem Model (FSAM) as detailed hereunder.

Using FSPM at the plant scale

In step 2, each mockup element becomes a node of the graph and is edge-connected to others depending of their relation (such as adjacency). The element can then be attached to a FSPM. FSPMs offer a realistic and/or volumetric visual representation that can be projected in AR. Depending on their assumptions, models can be integrated into a tool and mobilized for decision making [Surový et al, 2013], or dynamics understanding, as shown for light interception [START_REF] Streit | [END_REF].

Mobilizing a FSPM approach at the system scale

The graph modelling form enables us to expand the FSPM principles at the scale of the whole AFS, building so-called FSAM. In FSAM, functions are based on the AFS pattern and are linked to ecosystem services. With graph structural descriptions, we can search for or identify patterns (such as tree alignments) within the system; we can also design a system by combining existing patterns with desired functions.

Implementation

We are currently implementing a prototype application using these concepts. In particular, steps 2 and 3 (which include function-structure aspects) can be detailed as:

-From the captured mockup, a raster model and a graph model are built -Tree alignment (a simple pattern) is recognized -AFS designers are able to interactively visualize, with AR, the 3D virtual plants (built from FSPM simulations) corresponding to the mockup elements -AFS designers are able to visualize the virtual mockup in situ -AFS designers can visualize the tree growth over the time

Conclusion

With the structural information in graph form, some ecosystem functions can be considered during co-design workshops. Using FSPM at the plant and system scale combined to AR, the current prototype keeps interactivity during co-design workshops while improving visualization and understanding of the system thanks to the system simulation. To do this, we developed a new model (FSAM) using the structures and functions of the AFS. This work has received support from French ANR-16-CONV-0004 #DigitAg program


Fabien Verniest
email: fabienverniest@edu.mnhn.fr

Isabelle Le Viol

Romain Julliard

Laura Dami

Anis Guelmami

Marie Suet

Wed Abdou

Hichem Azafzaf

Nadjiba Bendjedda

Taulant Bino

John J Borg

Luka Božič

Mohamed Dakki

Rhimou El Hamoumi

Vitor Encarnação

Kiraz Erciyas-Yavuz

Khaled Etayeb

Valeri Georgiev

Ayman Hamada

Ohad Hatzofe

Christina Ieronymidou

Tom Langendoen

Tibor Mikuska

Blas Molina

Filipe Moniz

Caroline Moussy

Asmaâ Ouassou

Nicky Petkov

Danae Portolou

Tareq Qaneer

Samir Sayoud

Marko Šćiban

Goran Topić

Danka Uzunova

Gal Vine

Andrej Vizi

Erald Xeka

Marco Zenatello

Elie Gaget

Thomas Galewski

Anticipating the effects of climate warming and natural habitat conversion on waterbird communities to address protection gaps

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

INTRODUCTION

Many studies have documented how recent climate change, and especially temperature increases, has affected biodiversity in the past decades: alteration of phenologies, shifts of distributions and abundances toward poles and higher elevations, disruption of interactions and modifications of structure and composition of communities [START_REF] Kappelle | Effects of climate change on biodiversity: A review and identification of key research issues[END_REF][START_REF] Parmesan | Ecological and Evolutionary Responses to Recent Climate Change[END_REF][START_REF] Root | Fingerprints of global warming on wild animals and plants[END_REF]. Developing future scenarios of biodiversity under climate change has also become a central issue in the last decades. However, in many studies with this aim, the future impact of land-use changes on biodiversity through habitat loss and degradation has been overlooked [START_REF] Santos | The role of land use and land cover change in climate change vulnerability assessments of biodiversity: A systematic review[END_REF][START_REF] Titeux | Biodiversity scenarios neglect future land-use changes[END_REF]; but see [START_REF] Asamoah | Land-use and climate risk assessment for Earth's remaining wilderness[END_REF][START_REF] Newbold | Future effects of climate and land-use change on terrestrial vertebrate community diversity under different scenarios[END_REF] for instance), although it was the main driver of species declines and extinctions during the past decades and is expected to remain so in the coming years [START_REF] Caro | An inconvenient misconception: Climate change is not the principal driver of biodiversity loss[END_REF][START_REF] Ipbes | Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services[END_REF][START_REF] Newbold | Global effects of land use on local terrestrial biodiversity[END_REF]. Therefore, it is crucial to include these two anthropogenic pressures when assessing the future effects of global changes on biodiversity, especially as their effects may interact [START_REF] Mantyka-Pringle | Interactions between climate and habitat loss effects on biodiversity: A systematic review and meta-analysis[END_REF].

Protected areas (PAs) can mitigate the effects of climate change, land-use changes and their interaction on biodiversity. Previous studies demonstrated that PAs can facilitate the colonization of new areas by species shifting their distribution poleward [START_REF] Thomas | Protected areas facilitate species' range expansions[END_REF] and reduce the extinction risk of species at their trailing distribution edge [START_REF] Lehikoinen | Protected areas act as a buffer against detrimental effects of climate change-Evidence from large-scale, long-term abundance data[END_REF][START_REF] Peach | Value of protected areas to avian persistence across 20 years of climate and land-use change[END_REF]. The effectiveness of PAs in mitigating the effects of climate change on biodiversity is partly attributable to their ability to provide high-quality habitats and reduce other anthropogenic pressures [START_REF] Thomas | The performance of protected areas for biodiversity under climate change: Protected Areas Under Climate Change[END_REF], such as habitat loss and degradation [START_REF] Geldmann | Effectiveness of terrestrial protected areas in reducing habitat loss and population declines[END_REF][START_REF] Leberger | Mediterranean wetland conservation in the context of climate and land cover change[END_REF]. Therefore, it is imperativealthough rarely done (but see [START_REF] Verniest | Coupling future climate and land-use projections reveals where to strengthen the protection of Mediterranean Key Biodiversity Areas[END_REF])to explicitly consider these two anthropogenic pressures by anticipating their future intensity and distribution when identifying potential new PAs [START_REF] Heller | Biodiversity management in the face of climate change: A review of 22 years of recommendations[END_REF][START_REF] Pörtner | IPBES-IPCC co-sponsored workshop report on biodiversity and climate change[END_REF], especially since we will need to designate many new PAs to meet ambitious targets of protecting 30% of land by 2030 within the Kunming-Montreal global biodiversity framework [START_REF] Cbd | Kunming-Montreal Global biodiversity framework[END_REF].

The conservation of natural wetlands in the Mediterranean region is of major importance as these ecosystems are key to many species and highly threatened by climate and land-use changes [START_REF] Galewski | Living Mediterranean Report -Monitoring species trends to secure one of the major biodiversity hotspots[END_REF]MWO-2, 2018;[START_REF] Perennou | Evolution of wetlands in Mediterranean region. Book chapter in Water Resources in the Mediterranean Region[END_REF]. Wetlands are among the most threatened ecosystems on the planet, with a loss of 90% of them since 1700 as result of human activities (Convention on Wetlands, 2021; WWF, 2020). In the Mediterranean region, almost half of these ecosystems have been lost since 1970 (MWO-2, 2018) because of anthropogenic pressures such as climate and land-use changes. As these changes are heterogeneously distributed (MWO-2, 2018;[START_REF] Leberger | Mediterranean wetland conservation in the context of climate and land cover change[END_REF] and expected to remain so and to increase in the coming decades in this region [START_REF] Drobinski | How warmer and drier will the Mediterranean region be at the end of the twenty-first century?[END_REF][START_REF] Malek | Global change effects on land management in the Mediterranean region[END_REF][START_REF] Medecc | Climate and Environmental Change in the Mediterranean Basin -Current Situation and Risks for the Future[END_REF][START_REF] Taylor | The future for Mediterranean wetlands: 50 key issues and 50 important conservation research questions[END_REF], anticipating future climate and landuse changes is a major issue for spatial conservation planning of Mediterranean wetlands [START_REF] Taylor | The future for Mediterranean wetlands: 50 key issues and 50 important conservation research questions[END_REF]. This is especially true given that Mediterranean biodiversity is very sensitive to both these pressures [START_REF] Newbold | Tropical and Mediterranean biodiversity is disproportionately sensitive to land-use and climate change[END_REF], is of exceptional value and highly dependent on wetlands (MWO-2, 2018;[START_REF] Geijzendorffer | Mediterranean wetlands: a gradient from natural resilience to a fragile socialecosystem[END_REF][START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF], and Mediterranean PAs network suffers important gaps in wetland conservation [START_REF] Geijzendorffer | A More Effective Ramsar Convention for the Conservation of Mediterranean Wetlands[END_REF][START_REF] Leberger | Mediterranean wetland conservation in the context of climate and land cover change[END_REF][START_REF] Popoff | Gap analysis of the Ramsar site network at 50: Over 150 important Mediterranean sites for wintering waterbirds omitted[END_REF].

Identifying waterbird communities that might experience limitations in thermal adjustment to future climate warming can contribute to effective designation of additional PAs in wetlands. Waterbirds are species that are ecologically dependent on wetlands (Ramsar Convention) for at least part of their annual cycle and of which they are a good indicator of the conservation status. Conservation measures, such as protected areas and management measures, contribute to the conservation of waterbird populations at global scale [START_REF] Amano | Successful conservation of global waterbird populations depends on effective governance[END_REF][START_REF] Wauchope | Protected areas have a mixed impact on waterbirds, but management helps[END_REF] and in the Mediterranean region [START_REF] Gaget | Assessing the effectiveness of the Ramsar Convention in preserving wintering waterbirds in the Mediterranean[END_REF][START_REF] Kleijn | Waterbirds increase more rapidly in Ramsar-designated wetlands than in unprotected wetlands[END_REF] by mitigating the impacts of the anthropogenic pressures they face, such as climate change [START_REF] Gaget | Waterbird communities adjust to climate warming according to conservation policy and species protection status[END_REF][START_REF] Gaget | Assessing the effectiveness of the Ramsar Convention in preserving wintering waterbirds in the Mediterranean[END_REF]. Waterbirds have been extensively monitored for decades through the International Waterbird Census (IWC) [START_REF] Delany | Guidelines for participants in the International waterbird census (IWC)[END_REF] to inform many international conservation policies (e.g., Ramsar Convention, Convention on Migratory Species, Barcelona Convention, and European Union Directives, such as Birds, Water

Framework and Marine Strategy Directives). The IWC data also provided invaluable knowledge to conservation biology (e.g., [START_REF] Amano | Successful conservation of global waterbird populations depends on effective governance[END_REF][START_REF] Johnston | Observed and predicted effects of climate change on species abundance in protected areas[END_REF][START_REF] Wauchope | Protected areas have a mixed impact on waterbirds, but management helps[END_REF]. Using IWC data, Gaget et al. (2020a) found that thermal adjustment of wintering waterbird communities in a climate warming context, i.e., changes in community composition in response to climate warming, is mainly driven by thermal specialists, but is prevented by natural habitat conversion. However, PAs were also found to improve communities adjustment to climate warming [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF][START_REF] Pavón-Jordán | Climate-driven changes in winter abundance of a migratory waterbird in relation to EU protected areas[END_REF].

Despite the potential importance of these results for conservation planning, they have not yet been used to project the potential adjustment of waterbirds to future temperatures.

In this study, we highlight priorities for wetland protection to improve waterbird responses to climate warming in the Mediterranean region. To this end, we assessed the thermal adjustment limitations that non-breeding waterbird communities (151 species) at 2932 IWC sites in 21 Mediterranean countries might experience by the end of the 21 st century because of the interacting effect of climate warming and land-use changes. Assuming that waterbird communities can adjust their community composition to climate warming, but that this adjustment is not fast enough to cope with the velocity of the temperature increase [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF], and that natural habitat conversion is an additional limit to waterbird community adjustment to climate warming (Gaget et al., 2020a), we investigated whether communities that might experience the greatest thermal adjustment limitations to future climate warming are protected, with a particular focus on wetlands of international importance for waterbirds.

To this end, we applied a trait-based Climate Change Vulnerability Assessment framework (CCVA) [START_REF] Pacifici | Assessing species vulnerability to climate change[END_REF], which can be used to prioritize the implementation of spatial conservation measures [START_REF] Foden | Climate change vulnerability assessment of species[END_REF], to combine: i) the thermal specialization of waterbird communities derived from species thermal ranges; ii) the exposure of waterbird communities to climate warming; and iii) their exposure to natural habitat conversion. To ensure the robustness of this prioritization, exposure metrics were computed using late 21 st century (2081 -2100) projections under four recent scenarios [START_REF] O'neill | The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6[END_REF] that have very different spatial patterns of land-use in the Mediterranean region [START_REF] Verniest | Coupling future climate and land-use projections reveals where to strengthen the protection of Mediterranean Key Biodiversity Areas[END_REF].

In a context of adaptation to climate change, this study provides spatial prioritization that highlights wetlands of importance for waterbird communities that might be protected to help mitigate the negative impacts of climate warming on waterbirds and expand the Mediterranean protected areas network.

MATERIAL AND METHODS

Waterbirds monitoring

We used winter abundance data collected from the International Waterbird Census [START_REF] Delany | Guidelines for participants in the International waterbird census (IWC)[END_REF] in Mediterranean countries from 2008 to 2017. IWC surveys were conducted once a year in mid-January by professional and skilled volunteers following a standardized survey protocol and coordinated by Wetlands International (http://www.wetlands.org). Sites with fewer than 3 surveys during the study period or without climate data were deleted, and species that do not overwinter in Mediterranean countries were excluded to avoid vagrant species (see Appendix S1 for more information on species and site selection). The above-mentioned criteria resulted in a dataset of 151 waterbird species monitored at 2932 sites (average number of surveys per site: 6.0) in 21 countries (Appendix S2). Countries were assigned to four subregions: Balkans, Maghreb, Middle-East and Western Europe (Appendix S2).

We used the most up-to-date spatial database of PAs for Mediterranean countries, i.e., the Mediterranean Wetlands Observatory database, that combines national inventories with multiple international data sources, such as the World Database of Protected Areas (https://www.protectedplanet.net), the Common Database on Designated Areas (https://www.eea.europa.eu) and Natura 2000 (https://www.eea.europa.eu). Half of the sites were considered protected (n = 1526, 52.0%), i.e., their geographical features (point, line or polygon) partially or entirely overlapped a PA. One third of protected sites had an IUCN management category from I to IV (n = 534, 18.2%, hereafter referred to as 'strictly protected sites'). Protected sites with other IUCN management category (i.e., V or VI) or no IUCN management category were referred to as 'other protected sites' (n = 992, 33.8%). Sites that did not overlap a PA were referred to as 'non-protected sites' (n = 1406, 48.0%).

We considered 233 sites (7.9%) as of international importance for waterbirds according to the definition of the Ramsar convention (i.e., sites that regularly support one threatened waterbird species, 20,000 waterbirds, or 1% of a population) and based on the results of [START_REF] Popoff | Gap analysis of the Ramsar site network at 50: Over 150 important Mediterranean sites for wintering waterbirds omitted[END_REF].

Thermal specialization of waterbird communities

We evaluated the thermal specialization of waterbirds at the community level (hereafter referred to as 'SpecT°') from species thermal ranges. For each species, we computed the Species Temperature Range (STR) as the difference between the maximum (90% upper) and the minimum (10% lower) temperatures of its wintering thermal range, which was assessed following Gaget et al. (2020aGaget et al. ( ), i.e., using long-term (1950Gaget et al. ( -2000) ) January temperatures across its non-breeding distribution (BirdLife International & HBW, 2017). Therefore, high STR species experience a wide range of winter temperatures while low STR species experience a narrow range of winter temperatures.

For each site survey, we then calculated the Community Temperature Range (CTR), a weighted mean index that informs the average thermal range of species occurring in a community. For a community defined in time and space, the CTR corresponds to the average STR of species within the community weighted by their log(x + 1)-transformed abundance [START_REF] Godet | Waders in winter: Long-term changes of migratory bird assemblages facing climate change[END_REF]. We computed its inverse to describe the thermal specialization of waterbird communities, i.e., SpecT°. Therefore, high SpecT° values indicate that the community has a high thermal specialization and that the average thermal range of a community is low. Finally, we calculated site SpecT° averaged from site surveys SpecT° values.

Environmental data

For both climate and land-use data, we used future projections under four CMIP6 scenarios (SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF] to ensure that the identification of sites for priority protection considers different future pathways. These scenarios combine future socio-economic conditions with radiative forcing levels in 2100 [START_REF] O'neill | The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6[END_REF].

Past (2008 -2017) and future (2081 -2100) mean temperature in January were extracted from the WorldClim 2.1 dataset [START_REF] Fick | WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas[END_REF]; past temperature originally from CRU-TS 4.03, [START_REF] Harris | Updated high-resolution grids of monthly climatic observations -the CRU TS3.10 Dataset[END_REF]. Temperature data was mean-aggregated from a 2.5 arc-minutes resolution to a 15 arc-minutes resolution (approx. 28 x 28 km at the equator) to match land-use data spatial resolution (see below). Following [START_REF] Foden | Climate change vulnerability assessment of species[END_REF], we used for each scenario the future projections from the eight General Circulation Models (GCMs) available to date in WorldClim 2.1 (Appendix S3), and replicated all processing steps described in the following sections eight times (once per GCM) to assess the potential limitations of waterbird communities in adjusting their composition to future climate warming for each scenario and each GCM.

We extracted past (2012 -2013, median years of 2008 -2017) and future (2090, median of 2081 -2100) proportions of land-use categories from the Land Use Harmonization 2 dataset [START_REF] Hurtt | Harmonization of global land use change and management for the period 850-2100 (LUH2) for CMIP6[END_REF]. This dataset provides yearly estimates of the proportion of 12 land-use categories at a 15 arc-minutes spatial resolution.

Exposure of sites to climate warming and natural habitat conversion

To assess the exposure of surveyed sites to climate warming (hereafter 'ExpoT°'), we contrasted past temperatures experienced by communities to future temperatures, weighted by past temperature variations. First, we computed the 2008 -2017 average and standard deviation from the yearly mean temperature in January. Second, we then quantified for each scenario the differences between past and future temperatures using the Standardized Euclidian Distance as in [START_REF] Williams | Projected distributions of novel and disappearing climates by 2100 AD[END_REF]:

𝐸𝑥𝑝𝑜𝑇°𝑖 = √ (𝑎 𝑖 -𝑏 𝑖 ) 2 𝑆 𝑖 2 (1)
where ai and bi are the past and future mean temperatures in January at cell i, and Si is the past standard deviation.

To assess the exposure to natural habitat conversion (hereafter 'ExpoNHC'), we contrasted past and future land-use. We computed the proportion of anthropized habitats by summing the proportion of seven land-use categories: urban land, C3 annual crops, C3 perennial crops, C4 annual crops, C4 perennial crops, C3 nitrogen-fixing crops and managed pasture. Conversely, primary land (forested and non-forested), secondary land (forested and non-forested) and rangeland were considered as natural habitats. The past proportion of anthropized habitats was assessed as the average of years 2012 and 2013, while the future proportion of anthropized habitats was 2090. We then calculated for each scenario ExpoNHC as dissimilarities between the two periods at cell i by subtracting the past to the future proportion of anthropized habitats.

Index of future thermal adjustment limitations

To assess future thermal adjustment limitations of waterbird communities, we combined their thermal specialization (SpecT°) with the two exposure metrics (ExpoT° and ExpoNHC) using a trait-based vulnerability assessment framework [START_REF] Pacifici | Assessing species vulnerability to climate change[END_REF].

Following [START_REF] Huntley | Chapter 5. Using CCVAs and interpreting their results[END_REF], we first converted ExpoNHC, ExpoT°, and SpecT° into three-level ordinal variables (Low, Moderate, High). Because 5% of natural habitat conversion has been shown to prevent community adjustment to climate warming of nonbreeding waterbirds in a Mediterranean context (Gaget et al., 2020a), we used a 5% threshold to convert ExpoNHC into an ordinal variable: Low: inferior to 0%, Moderate: between 0% and 5%, High: superior to 5%. ExpoT° was converted into an ordinal variable using two quantiles as thresholds (provided in Appendix S4). To account for differences in exposure intensity between scenarios, the two quantiles were calculated from the values of all scenarios.

Similarly to ExpoT°, we converted SpecT° into a three-level ordinal variable using the two quantiles method (Appendix S4).

We then computed an index of future thermal adjustment limitations (hereafter 'LimT° index') by combining the three ordinal variables using a rule-based logic framework (Table 1). This approach allowed us to account for interactions between the two exposure metrics and SpecT° [START_REF] Huntley | Chapter 5. Using CCVAs and interpreting their results[END_REF]. Through the combinations used to compute the LimT° index, we assumed that natural habitat conversion was the main factor that could limit community adjustment to climate warming (Gaget et al., 2020a), i.e., that ExpoNHC was the main driver of the LimT° index. As the temperature increase can be too fast for the adjustment of waterbird communities to climate warming [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]), we also considered that higher ExpoT° values would exacerbate limitations caused by natural habitat conversion in community adjustment to climate warming. Finally, community adjustment to climate warming is mainly driven by thermal specialist species under low natural habitat conversion rates, through the colonization of warm-dwelling specialists and the extirpation of cold-dwelling specialists (Gaget et al., 2020a). Therefore, under low ExpoNHC conditions, we assumed that communities with high SpecT° values, i.e., with a higher representation of species with a narrow thermal range than species with a wide thermal range, would adjust to future temperature increase more easily than communities with low SpecT°. In contrast, thermal specialists and generalists both contribute to the non-significant community adjustment to climate warming in high rates of natural habitat conversion conditions, with the extirpation of thermal specialists and the colonization of thermal generalists regardless of their thermal affinities (Gaget et al., 2020a). Therefore, under high ExpoNHC conditions, we considered that communities with high SpecT° would not have a more positive thermal adjustment than communities with low SpecT°.

Following [START_REF] Huntley | Chapter 5. Using CCVAs and interpreting their results[END_REF], we assessed the robustness of our findings to variations in the rule-based logic framework presented below (Table 1) by performing all analyses with slightly different combinations of levels (Low, Moderate, High). This sensitivity analysis resulted in very similar results to those presented in the main body of this article and are provided in Appendix S5.

Data Analysis

For each future scenario, we built a multi-model of the LimT° index by assigning to each site the most frequent LimT° value (i.e., the mode) across the eight GCMs. In case of a tie, the highest level was assigned. This multi-model was used for the statistical analyses detailed below.

We identified protection priorities by comparing LimT° index values between site protection statuses (strictly protected, other protected, non-protected), statuses of international importance for waterbirds (internationally important, or not internationally important) and future scenarios. The effect of each variable was assessed separately and in combination with a subregion effect and their interaction to assess differences between subregions. To this end, we used on each individual variable Proportional Odds Logistic Regression models (POLRs) that can be applied to ordinal response variables such as the LimT° index [START_REF] Guisan | Ordinal response regression models in ecology[END_REF].

Similarly, we compared ExpoT°, ExpoNHC and SpecT° values between protection statuses, statuses of international importance and future scenarios (only for ExpoT° and ExpoNHC) using linear models (LMs) in order to explore the potential drivers of differences of the LimT° index. In addition, to evaluate similarities in spatial patterns of exposure between scenarios, we assessed the correlation between scenarios for ExpoT° and ExpoNHC with a Kendall's rank correlation tests [START_REF] Kendall | Rank correlation methods[END_REF]. Finally, we assessed the correlation between ExpoT°, ExpoNHC and SpecT° for each scenario to evaluate similarities in spatial patterns between the three components of the LimT° index. These analyses were carried out using an average multi-model across the eight GCMs for ExpoT°.

Statistical significance was assessed using 95% confidence intervals for POLR and LMs

and was set at P < 0.05 for correlation tests. POLRs were fitted with the MASS package [START_REF] Venables | Modern Applied Statistics with S. Fourth Edition Springer[END_REF]. Statistical analyses were performed using R 4.1.0 (R Core Team, 2021). Geoprocessing operations were performed using QGIS 3.4.15 (QGIS Development Team, 2020) and R 4.1.0 (R Core Team, 2021), mainly with the sf [START_REF] Pebesma | Simple Features for R: Standardized Support for Spatial Vector Data[END_REF] and raster [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling[END_REF] packages.

RESULTS

Future thermal adjustment limitations of waterbird communities

The LimT° index was 'High' in at least one scenario for more than a third of sites (39.5%, n = 1158) and in all scenarios for only 14 sites (0.5%, Figure 1). Values of the LimT° index increased with the radiative forcing level of scenarios in Western Europe and the Balkans, except for the SSP3-7.0 scenario for which the index took lower values than in SSP2-4.5 and SSP5-8.5 (Appendices S6 and S7). However, we found very different results for other subregions (Appendices S6 and S7), with, for instance, the lowest values of the LimT° index for the Maghreb sites obtained for the SSP5-8.5 scenario. The number of sites with a high LimT° index was lowest for SSP5-8.5 (6.6%, n = 193, Appendices S6 and S8). The variability in LimT° index between GCMs was very low in all scenarios (Appendix S9).

Future thermal adjustment limitations of waterbird communities at sites of international importance

The sites of international importance for waterbirds in all subregions had higher values of LimT° index than other sites in two scenarios (Appendices S10 and S11). More than half of the sites of international importance (n = 131, 56.2%) had a 'High' LimT° index in at least one scenario, and most of them were protected (75.6%, n = 99, Figures 1 and2). Those that were not protected (n = 32) were mostly located in Turkey, the south of the Iberian Peninsula, and the coast of North-West Africa (Figure 2). Only two sites of international importance (0.9%) had a 'High' LimT° index in all scenarios and both were protected (Figure 1): 'Laghi Como:

Garlate: Olginate' (Italy) and 'Lago Maggiore' (Italy).

ExpoT° and ExpoNHC were weakly correlated in SSP1-2.6, SSP2-4.5 and SSP3-7.0 (Mean absolute τ = 0.09, P < 0.001, Appendix S19). In contrast, in SSP5-8.5, ExpoT° and ExpoNHC were strongly and negatively correlated (τ = -0.36, P < 0.001, Appendix S19), suggesting that sites with high ExpoT° have low ExpoNHC.

Thermal specialization of waterbird communities

SpecT° was higher at sites of international importance for waterbirds than at other sites except for Maghreb countries (Appendix S19). SpecT° was higher at strictly protected sites than at other protected sites and non-protected sites in all subregions (Appendix S20). SpecT° was higher for Balkans and Maghreb communities than for Western Europe and Middle-East communities (Appendices S21 and S22).

SpecT° and ExpoT° were very weakly correlated in SSP1-2.6 and SSP2-4.5 (-0.04 < τ < -0.03, P < 0.01, Appendix S19) and were not correlated in SSP3-7.0 and SSP5-8.5. SpecT° and ExpoNHC were very weakly correlated in SSP1-2.6 and SSP5-8.5 (-0.06 < τ < -0.04, P < 0.01, Appendix S19) and were not correlated in SSP2-4.5 and SSP3-7.0.

DISCUSSION

In this study, we identified for the first time the wintering waterbird communities that might experience the greatest limitations to adjust to climate warming by the end of the century, and demonstrated that current Mediterranean PAs, and especially strictly PAs, were overall well located to support their thermal adjustment. However, we identified 490 entirely nonprotected sites whose communities may experience high limitations in thermal adjustment, including 32 sites of international importance (sensu Ramsar Convention). Building on the importance of PAs for future biodiversity conservation in response to both climate warming and natural habitat conversion, our results reinforce the need for PAs designation and acknowledge that decades of wetland conservation efforts are expected to help waterbird responses to climate warming.

Thermal adjustment and vulnerability to climate warming

Our approach presumes that waterbird communities with a better thermal adjustment are less vulnerable to climate warming. Vulnerability to climate change is defined as "the extent to which biodiversity is susceptible to or unable to cope with the adverse effects of climate change" (IPCC, 2007;[START_REF] Foden | Climate change vulnerability assessment of species[END_REF]. As a limited thermal adjustment to warmer conditions could lead to increased climatic debt [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF][START_REF] Soultan | The future distribution of wetland birds breeding in Europe validated against observed changes in distribution[END_REF], it may be considered as a factor of vulnerability to climate change. Wintering waterbird communities adjust to climate warming by species composition changes mainly through thermal specialist species (i.e., high SpecT° species), with the colonization of communities by warm-dwelling specialists and the extirpation of cold-dwelling specialists from waterbird communities (Gaget et al., 2020a). Therefore, communities with a high proportion of thermal specialists should adjust better to climate warming, and thus should be less vulnerable, than communities with a high proportion of thermal generalists (i.e., low SpecT° species).

However, considering future limitations in community thermal adjustment (i.e., LimT° index) as a proxy of vulnerability to climate warming must be done with great caution as its interpretation and estimation raise questions. Thermal generalist species may not adjust to temperature increases because they may be less sensitive and may have a greater adaptive capacity than thermal specialists, and consequently be less vulnerable to climate warming despite their higher climatic debt. Furthermore, changes in SpecT°, a community weightedmean index, only reflect changes in the relative proportion of specialists and generalists, not their presence, and other mechanisms such as interspecific competition and anthropogenic disturbance can also influence thermal adjustment (Liang et al., 2017). Finally, the spatial resolution (approx. 28 x 28 km at the equator) of the land-use data is very coarse and its very low thematic resolution (12 categories) does not include a definition of wetlands. This results in a large uncertainty on the relationship between natural habitat conversion and the monitored wetlands, as well as in a potential overestimation of the exposure to climate warming (ExpoT°) [START_REF] Heikkinen | Fine-grained climate velocities reveal vulnerability of protected areas to climate change[END_REF]. Improved future projections of land-use are required to accurately quantify waterbird exposure to natural habitat conversion.

Nevertheless, these limitations do not diminish the significance of our study because they do not undermine our assumption that communities impacted by high rates of natural habitat conversion are less able to adjust to climate warming (Gaget et al., 2020a). Given the rule-based logic framework used to compute the LimT° index, the limitations based on the thermal specialization of the community also do not invalidate the higher vulnerability of communities exposed to high rates of natural habitat conversion compared to communities not exposed to natural habitat conversion. Finally, the combinations of this framework were somewhat arbitrarily decided, but very similar findings were found with different combinations in the sensitivity analysis (Appendix S5).

Future scenarios of thermal adjustment

Although one might anticipate greater limitations in community thermal adjustment to climate warming at high radiative forcing, this pattern was only partially observed and not even across all subregions. This is because the exposure to land-use changes does not necessarily increase with the radiative forcing level of the scenarios and has very different spatial patterns between CMIP6 scenarios in the Mediterranean region [START_REF] Verniest | Coupling future climate and land-use projections reveals where to strengthen the protection of Mediterranean Key Biodiversity Areas[END_REF]. Therefore, the lower Although it might be appealing to infer that high radiative forcing scenarios are not that problematic, such conclusion must be tempered. Indeed, other CMIP6 scenarios combining the same high radiative forcing levels (i.e., 7.0 W.m -2 and 8.5 W.m -2 ) with other socio-economic conditions (e.g., SSP2), for which data are not yet available, could result in a very different future land-use and thus potentially in higher values of the LimT° index, as

ExpoNHC is its main driver. Furthermore, waterbird community thermal adjustment is far to be perfect and results in a climatic debt even at lower rates of temperature increase [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]. Our research reinforces the importance of considering future projections of habitat loss and degradation, whose results against scenarios are far less straightforward than for climate projections.

Addressing gaps in the protected areas network

Using a trait-based CCVA framework, we highlighted that waterbird communities of wetlands overlapped by strictly PAs (IUCN category from I to IV) would experience greater limitations to adjust to climate warming than non-protected waterbird communities, certainly due to the higher exposure of strictly PAs to natural habitat conversion. We also found that many of the communities that may experience the greatest limitations are protected: among the 14 sites whose waterbird communities had a high LimT° index in all scenarios, 10 located in western Europe (France, Italy, Spain), including the two sites of international importance for waterbirds, are partially or entirely protected. We believe these finding to be encouraging in two respects. First, it implies that PAs are generally located in wetlands where the need to help waterbird communities adjust to new thermal conditions will be most critical. Second, future projections of land-use might overestimate ExpoNHC in PAs, since Integrated Assessment Models (IAMs) [START_REF] Harfoot | Integrated assessment models for ecologists: The present and the future[END_REF] used to produce these projections do not explicitly account for PAs and have a very coarse spatial resolution. This is especially the case for strictly PAs that have more restrictive management than other PAs [START_REF] Dudley | Guidelines for applying protected area management categories[END_REF], and thus, we assume, have a better ability to reduce anthropogenic pressures.

Nevertheless, we assume that many sites that were considered as protected in this study are only partially overlapped by PAs [START_REF] Popoff | Gap analysis of the Ramsar site network at 50: Over 150 important Mediterranean sites for wintering waterbirds omitted[END_REF], although this could not be verified throughout the entire study area because the exact delineation of the IWC monitored sites is not available in many cases [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]. Therefore, we advocate that the delineation of each IWC site be completed and made available. This would enable the expansion of strictly PAs in wetlands hosting waterbird communities that might experience limitations to adjust to climate warming and that are only partially covered by a strict protection status (i.e., strict nature reserve, wilderness area, national park, natural monument or feature, and habitat/species management area, [START_REF] Dudley | Guidelines for applying protected area management categories[END_REF]. In addition, a decline in the effectiveness of PAs to reduce anthropogenic pressures such as natural habitat conversion, which was observed in some regions [START_REF] Geldmann | A global-level assessment of the effectiveness of protected areas at resisting anthropogenic pressures[END_REF], might threaten the ability of PAs to facilitate the thermal adjustment of waterbird communities to climate warming. We thus recommend the implementation of effective management plans in current PAs to ensure their ability to support waterbird responses to climate change [START_REF] Gaget | Protected area characteristics that help waterbirds respond to climate warming[END_REF] and positive waterbird population trends [START_REF] Wauchope | Protected areas have a mixed impact on waterbirds, but management helps[END_REF].

Despite the above-mentioned result, we identified 490 waterbird communities that might experience significant thermal adjustment limitations in at least one scenario and that are not included in the current PAs network (Supplementary Material). Maghreb and Middle-East countries account for 199 of these communities (40.6%), including 22 at sites of international importance for waterbirds, notwithstanding the low number of selected sites in these countries (n = 586, 20.0%). Our results are coherent with [START_REF] Breiner | Setting priorities for climate change adaptation of Critical Sites in the Africa-Eurasian waterbird flyways[END_REF] in which authors found that African and Middle-East Critical Sites for waterbirds are key to support waterbird responses to climate change. Such sites are also important to protect to ensure the connectivity of the wetlands along the Eurasian-African migration flyway [START_REF] Deboelpaep | Highway(s) overhead: Strong differences in wetland connectivity and protected status challenge waterbird migration along the four Palearctic-Afrotropical flyways[END_REF]. Our results urge Mediterranean countries to quickly designate new PAs to protect the 490 communities we identified, with priority given to the 32 wetlands of international importance for waterbirds [START_REF] Popoff | Gap analysis of the Ramsar site network at 50: Over 150 important Mediterranean sites for wintering waterbirds omitted[END_REF]. Protecting these wetlands will help waterbird communities adjust to climate warming [START_REF] Gaget | Benefits of protected areas for nonbreeding waterbirds adjusting their distributions under climate warming[END_REF]

Perspectives

Prioritizing waterbird conservation should not be based only on wintering, i.e., non-breeding, species distribution but also account for the distribution during the breeding and migration periods (e.g., [START_REF] Deboelpaep | Highway(s) overhead: Strong differences in wetland connectivity and protected status challenge waterbird migration along the four Palearctic-Afrotropical flyways[END_REF]. We also recommend assessing the Black diamonds: sites of international importance for waterbirds (protected: n = 99, non-

  values of LimT° index in SSP3-7.0 in Western Europe and the Balkans can be explained by the low values of its main driver, i.e., ExpoNHC, taken in this scenario. The lowest values of the LimT° index in the Maghreb under SSP5-8.5 are due to the strongly negative ExpoNHC values on the coast of Morocco, Algeria, and Tunisia because of the extensive conversion of crops into secondary natural land.

  but might also benefit these ecosystems as a whole, which are among the most threatened on Earth and the least targeted by conservation measures (Abell & Harrison, 2020), as well as human populations (MWO-2, 2018). The designation of new PAs in wetlands whose waterbird communities might experience high thermal adjustment limitations to future climate warming should be facilitated by the new protection targets of the Kunming-Montreal global biodiversity framework (CBD, 2022).

FIGURESFigure 1 :

 1 FIGURES

Figure 2 :

 2 Figure 2: Protected (a) and non-protected (b) sites with high future thermal adjustment

Table 1 :

 1 TABLES Rule-based logic framework used to combine the exposure of sites to climate warming (ExpoT°), the exposure of sites to natural habitat conversion (ExpoNHC) and the thermal specialization of waterbird communities (SpecT°) to compute the LimT° index for each scenario (n = 4) and each GCM (n = 8). x = combination of two metrics, L = Low, M = Moderate, H = High. A community with a high thermal specialization (SpecT° = H) has a low index of future thermal adjustment limitations (LimT° index = L) if exposure to climate warming is high (ExpoT° = H) but exposure to natural habitat conversion is low (ExpoNHC = L). Under high exposure to natural habitat conversion, the index of future thermal adjustment limitations is always high.

	ExpoNHC x ExpoT°

ACKNOWLEDGEMENTS

We thank B. Maas and one anonymous reviewer for comments that improved the quality of the manuscript. We acknowledge the thousands of volunteers and professionals involved in waterbird counts and coordination units of the International Waterbird Census for providing help to national coordinators, centralizing and managing data and making this research possible. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modelling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF.

FUNDING

This work was supported by the Research Institute of Tour du Valat; the "TotalEnergies Fondation"; the Turku Collegium for Science, Medicine and Technology (grant EG); and the "Région Bretagne" [grant number 20SB273U7204_904A3].

KEYWORDS

Community thermal adjustment, Conservation planning, International Waterbird Census, Mediterranean region, Ramsar convention, Wetlands

Future thermal adjustment limitations of waterbird communities in protected areas

The LimT° index was higher for strictly protected sites than non-protected sites in at least three scenarios for all subregions (Appendices S10 and S12). Other protected sites also had also higher values of LimT° index than non-protected sites in at least two scenarios for the Maghreb, the Middle-East and Western Europe (Appendices S10 and S12). Almost half of protected sites (i.e., strictly protected and other protected sites) had a 'High' LimT° index in at least one scenario (n = 668, 43.8%, Figures 1 and2). On the other hand, 490 sites whose communities had a high LimT° index in at least one scenario were not protected (42.3%, Figures 1 and2). A list of these sites is provided as Supplementary Material. More than a third of these sites (40.6%, n = 199) were located in the Maghreb and the Middle-East (Figure 2).

Exposure of sites to climate warming and natural habitat conversion

The values of ExpoT° were higher (lower) at sites of international importance than at other sites in the Balkans and the Middle-East (Maghreb and Western Europe) in all scenarios (Appendix S13). Sites of international importance had higher ExpoNHC than other sites in three scenarios for Western Europe, but the differences for other subregions were less clear (Appendix S14). ExpoT° was overall higher at other protected sites than at non-protected sites except for Maghreb countries (Appendix S15). Strictly protected and other protected sites had overall higher ExpoNHC values than non-protected sites (Appendix S16).

ExpoT° values increased with the radiative forcing level of scenarios in all subregions and were highest for Maghreb sites (Appendix S17). Surveyed sites had overall lower ExpoNHC values in SSP3-7.0 than in other scenarios and had overall highest values in the Middle-East (Appendix S17). Sites had similar ExpoT° spatial patterns between scenarios (0.52 < τ < 0.91, P < 0.01, Appendix S18) but different ExpoNHC spatial patterns between scenarios (0.04 < τ < 0.25, P < 0.01, Appendix S18).

DECLARATION OF INTEREST STATEMENT

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.