P Fons Van Der

Sophia Ratcliffe

Paloma Ruiz-Benito

Michael Scherer-Lorenzen

Kris Verheyen

Christian Wirth

Miguel A Zavala

Evy Ampoorter

Lander Baeten

Luc Barbaro

Fons Van Der Plas
email: fonsvanderplas@gmail.com

Cristina C Bastias

Jürgen Bauhus

Raquel Benavides

Adam Benneter

Damien Bonal

Olivier Bouriaud

Helge Bruelheide

Filippo Bussotti

Monique Carnol

Bastien Castagneyrol

Yohan Charbonnier

Johannes H C Cornelissen

Jonas Dahlgren

Ewa Checko

Andrea Coppi

Seid Muhie Dawud

Marc Deconchat

Pallieter De Smedt

Hans De Wandeler

Timo Domisch

Leena Finér

Mariangela Fotelli

Arthur Gessler

André Granier

Charlotte Grossiord

Virginie Guyot

Josephine Haase

Stephan Hättenschwiler

Hervé Jactel

Bogdan Jaroszewicz

François-Xavier Joly

Tommaso Jucker

Stephan Kambach

Gerald Kaendler

Jens Kattge

Julia Koricheva

Georges Kunstler

Aleksi Lehtonen

Mario Liebergesell

Peter Manning

Harriet Milligan

Sandra Müller

Bart Muys

Diem Nguyen

Charles Nock

Bettina Ohse

Alain Paquette

Josep Peñuelas

Martina Pollastrini

Kalliopi Radoglou

Karsten Raulund-Rasmussen

Fabian Roger

Rupert Seidl

Federico Selvi

Jan Stenlid

Fernando Valladares

Johan Van Keer

Lars Vesterdal

Markus Fischer

Lars Gamfeldt

Eric Allan

Continental mapping of forest ecosystem functions reveals a high but unrealised potential for forest multifunctionality

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

DATA ACCESSIBILITY STATEMENT

Should the manuscript be accepted, the data supporting the results will be published in Dryad and the data DOI will be included at the end of the article.

ABSTRACT

Humans require multiple services from ecosystems, but it is largely unknown whether trade-offs between ecosystem functions prevent the realization of high ecosystem multifunctionality across spatial scales. Here, we combined a comprehensive dataset (28 ecosystem functions measured on 209 forest plots) with a forest inventory dataset (105,316 plots) to extrapolate and map relationships between various ecosystem multifunctionality measures across Europe. These multifunctionality measures reflected different management objectives, related to timber production, climate regulation and biodiversity conservation/recreation. We found that trade-offs among them were rare across Europe, at both local and continental scales. This suggests a high potential for "win-win" forest management strategies, where overall multifunctionality is maximized. However, across sites, multifunctionality was on average 45.8-49.8% below maximum levels and not necessarily highest in protected areas. Therefore, using one of the most comprehensive assessments so far, our study suggests a high but largely unrealized potential for management to promote multifunctional forests.

INTRODUCTION

One of the greatest challenges in ecology is to understand the effects of global change and nature management on the multiple ecosystem functions on which humans depend (MEA 2005). Such an understanding would help predicting the circumstances under which trade-offs between different ecosystem functions are minimal and therefore when their simultaneous provisioning, i.e. ecosystem multifunctionality [START_REF] Hector | Biodiversity and ecosystem multifunctionality[END_REF][START_REF] Gamfeldt | Multiple functions increase the importance of biodiversity for overall ecosystem functioning[END_REF]), is maximised. Previous studies have identified conditions promoting local-scale ecosystem multifunctionality, e.g. through the maximization of biodiversity [START_REF] Lefcheck | Biodiversity enhances ecosystem multifunctionality across trophic levels and habitats[END_REF].

However, whether such relationships also exist at large spatial scales, and how they vary in space, is less clear [START_REF] Isbell | Linking the influence and dependence of people on biodiversity across scales[END_REF]. Understanding this is essential if ecosystem-functioning studies are to provide policy-relevant advice, because most policy focuses on large scales.

Forests provide a number of functions related to key services such as timber production, climate regulation and recreation [START_REF] Gamfeldt | Higher levels of multiple ecosystem services are found in forests with more tree species[END_REF], and are important for the conservation of many plant and animal species (FAO 2015). Understanding large-scale relationships between different functions is therefore important if we are to find "win-win" management scenarios, which meet different forest management objectives and promote forest multifunctionality.

Quantifying many ecosystem functions at large scales has so far proven challenging. Forests are often managed for a particular subset of functions related to certain ecosystem services (e.g. timber production, climate regulation or nature conservation) that are prioritized by a specific stakeholder group. We aimed to identify areas where functions of all these different sets are high and where trade-offs are weakest. To this end, we combined a multi-site dataset, containing accurate measures of multiple ecosystem functions, with a continental-scale inventory-based dataset with high spatial plot coverage. We extrapolated regional scale relationships between ecosystem functions and their drivers (e.g. forest community composition and climate) to larger spatial scales (Fig. S1) to map both individual ecosystem functions and ecosystem multifunctionality across Europe, in forests without recent intensive management. We then tested for potential trade-offs between sets of functions, at scales relevant for policymakers.

To do this, we developed different measures of multifunctionality corresponding to different management scenarios (Fig. 1). In these, functions related to (sustainable) timber production, climate regulation or biodiversity conservation/recreation were prioritized (Fig. 1).

We also considered a scenario where all functions were valued equally. Our objectives were firstly, to identify "multifunctionality hotspots", i.e. areas with highest multifunctionality.

Secondly, we investigated whether there are synergies (allowing for win-win management) or trade-offs between different multifunctionality measures at both continental and local scales, and how these varied in space. Finally, we investigated whether forest protection status is associated with high multifunctionality, and thus whether potential win-win policies are realized in (protected) forests.

MATERIALS AND METHODS

Our approach to extrapolate ecosystem functioning relationships from regional to continental scales consisted of two main steps (Fig. S1). Firstly, statistical models were fitted to a comprehensive (many ecosystem functions), multi-site dataset ('fitting dataset'). Secondly, these models were extrapolated to a continental-scale dataset containing forest plots distributed across Western Europe ('inventory dataset'). These two datasets share variables related to climate, soils and tree composition, all potential drivers of ecosystem functioning. For three ecosystem functions which were independently measured in the inventory dataset, we cross-validated predicted ecosystem function values. Our approach allowed testing for trade-offs and synergies between individual ecosystem functions and between different multifunctionality measures, at different scales: 1) using all plots (thus including both local and large-scale variation in functions) and 2) within 20×20km localities.

Fitting dataset: design

As part of the EU-FP7 FunDivEUROPE project (www.fundiveurope.eu), which investigates how tree species composition and diversity drive forest ecosystem functioning, 209 30×30 meter plots (Fig. S2) were established. The plots covered six major regions/countries, representing different forest types: 28 boreal forest (Finland), 43 temperate mixed forest (Poland), 38 temperate deciduous forest (Germany), 28 mountainous deciduous forest (Romania), 36 thermophilous deciduous forest (Italy) and 36 Mediterranean mixed forest plots (Spain). These plots covered a broad climatic gradient: mean annual precipitation ranged from 484 to 819mm, mean annual temperature from 1.4 to 14.1°C (WorldClim; [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF] and altitude from 87 to 1404m. Within regions, plots differed in the composition and diversity of regionally common tree species, while site-related factors were similar. Management was either at low intensity or absent [START_REF] Baeten | A novel comparative research platform designed to determine the functional significance of tree species diversity in European forests[END_REF].

Measurement and collation of fitting data

In all plots, we measured 28 different ecosystem characteristics/processes ('ecosystem functions' hereafter) linked to various ecosystem services (see overview in Fig. 1 and methodology in Supplementary Material). For each plot we compiled data on tree species composition (to derive measures of functional and phylogenetic diversity), stand structure, soil pH, altitude and 18 climatic variables. Previous studies demonstrated that climate (Cramer et In each plot, we identified all tree stems ≥7.5cm in diameter at breast height (dbh) to species level. With these data, we calculated total and average tree basal area. In addition, by combining these observations with (1) published trait data [START_REF] Kattge | TRY -a global database of plant traits[END_REF]; Royal Botanic Gardens Kew 2015; see Table S1) representing key life-history strategies [START_REF] Westoby | Plant ecological strategies: some leading dimensions of variation between species[END_REF], and (2) a phylogeny [START_REF] Zanne | Three keys to the radiation of angiosperms into freezing environments[END_REF], we calculated several metrics describing the functional identity, functional diversity and phylogenetic diversity of the tree communities. Firstly, we calculated Community Weighted Means [START_REF] Garnier | Plant functional makers capture ecosystem properties during secondary succession[END_REF]), reflecting functional identities of communities, based on species values for specific leaf area (cm 2 g -1 ), maximum life span (log-transformed; yrs), maximum height (m), wood density (g cm -3 ), seed mass (mg), conifer (proportion) and evergreen (proportion). Secondly, we calculated the functional (trait) diversity within communities as Rao's Quadratic Entropy (Botta-Dukát 2005), for each trait separately and for all traits combined. Finally, we calculated several phylogenetic diversity metrics: Phylogenetic Species Variability, Phylogenetic Species Evenness [START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF]), Faith's Phylogenetic Diversity [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF]) and (abundance-weighted) Mean Phylogenetic Distance [START_REF] Webb | Phylogenies and Community Ecology[END_REF]. As inventory plots differed in size, tree species richness was not investigated, and we selected functional and phylogenetic diversity metrics uncorrelated with species richness.

To represent soil conditions we used pH (methods in Supplementary Materials), as it drives many functions and was the only soil variable available for the inventory dataset. Eighteen variables (see Table S2) related to climate (worldclim data; [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF]) were collated at a 30 seconds spatial resolution. Altitude data were collated from srtm.csi.cgiar.org.

Analysis of the drivers of ecosystem functioning

We used the Random Forest [START_REF] Breiman | Random forests[END_REF]) algorithm to explain ecosystem function variation in the fitting dataset. Random Forest is a machine-learning algorithm, powerful for making predictions (but less suitable in explaining mechanisms) and incorporating both linear and non-linear relationships, as well as interaction effects [START_REF] Strobl | Bias in random forest variable importance measures: illustrations, sources and a solution[END_REF]. It is relatively insensitive to multicollinearity and overfitting [START_REF] Hastie | The elements of statistical learning[END_REF], allowing for the inclusion of many predictors. Initially, we included the 42 predictor variables described above (see also Table S2), describing abiotic conditions, climate, stand structure, functional identity, and functional and phylogenetic diversity. Random Forests were run in R (R Core Team 2013) with the 'randomForest' library [START_REF] Liaw | randomForest: Breiman and Cutler's random forests for classification and regression[END_REF]. Following [START_REF] Seidl | Unraveling the drivers of intensifying forest disturbance regimes in Europe[END_REF], we iteratively removed those variables not reducing the mean square error over random permutations of the same variable. For final Random Forests, we identified, using the 'importance' function, the degree to which the inclusion of each predictor decreases residual model variance.

Forest inventory data

We combined data from 163,451 plots of the National Forest Inventories (NFIs) of Spain (59,048 plots), France [START_REF] Naidoo | Global mapping of ecosystem services and conservation priorities[END_REF]844), Wallonia (Belgium, 1,238), Germany [START_REF] Ratcliffe | Modes of functional diversity control on tree productivity across the European continent[END_REF]832), Sweden [START_REF] Elith | Species distribution models: ecological explanation and prediction across space and time[END_REF]212) and Finland [START_REF] Bateman | Bringing economic services into economic decision-making: land-use in the United Kingdom[END_REF]456). NFIs contained data on individual trees in each plot, including species identity, dbh and basal area. Furthermore, estimates of timber production (increase in tree basal area per hectare per year), tree biomass and tree recruitment (tree saplings per hectare) were available for many plots. To ensure that data from different NFIs were comparable to the fitting dataset plots, we only included trees with dbh ≥7.5cm. Furthermore, we only included the 105,316 plots that were at low to mid-altitudes (<1500m), without indication of recent logging, and dominated by one of the 'target' species of the fitting dataset [START_REF] Baeten | A novel comparative research platform designed to determine the functional significance of tree species diversity in European forests[END_REF].

We calculated the same climate, functional identity and functional and phylogenetic diversity variables for the NFI dataset as for the fitting dataset. Soil pH, calculated for the top 10cm of the soil at 1km 2 resolution, was obtained from the ESDAC database [START_REF] Panagos | European Soil Data Centre: Response to European policy support and public data requirements[END_REF]). These variables had similar ranges as in the fitting dataset (Table S3).

Extrapolating and mapping ecosystem functions across Europe

We used the 'predict' function in R to predict values of each ecosystem function in inventory plots, based on the Random Forests (built using the fitting dataset with independently collected FunDivEUROPE data; [START_REF] Baeten | A novel comparative research platform designed to determine the functional significance of tree species diversity in European forests[END_REF]) and the climate, functional identity, diversity (of the most recent survey) and abiotic conditions in the inventory plots. To determine the accuracy of our predictions, we correlated the three ecosystem functions (timber production, tree biomass and tree recruitment) that were measured in inventory plots with the values predicted by the Random Forests. We did the validations across all plots at continental scale (local and large scale variation) and within (only local variation) and among (only large-scale variation) 20×20km grid cells ('localities') containing ≥20 plots. In addition, we compared observed correlations between ecosystem functions with extrapolated ones. We also compared the average values for tree biomass and recruitment between fitting and inventory datasets (productivity was not comparable as it was measured in different units). To investigate how mapped functions changed across latitude, we fitted linear models with linear and quadratic effects of latitude as predictors.

Calculating multifunctionality and quantifying trade-offs

We used the 'threshold-approach' [START_REF] Gamfeldt | Multiple functions increase the importance of biodiversity for overall ecosystem functioning[END_REF]) to calculate ecosystem multifunctionality for each inventory plot, based on the predicted values of individual ecosystem functions. Ecosystem multifunctionality was measured at both local and continental scales and defined as the number of functions exceeding a threshold. The threshold was defined as the proportion (25%, 50% (default threshold reported in main results), 75% or 90%) of the 'maximum' value observed for that function, either within a 20×20km locality (local scale) or across Europe (continental scale). The maximum was defined as the 97.5 th percentile of observed functioning across plots, thus removing extreme outliers. For a concrete example on quantifying multifunctionality, we refer to Fig. S3. We excluded ecosystem functions that (a) had poor Ranfom Forest fit, with R 2 (correlation between observed and predicted) values <0.20 (default analysis; Fig. 1C) and (b), as a sensitivity analysis, also those which had a low validation R 2 (see results: tree recruitment and the related function of seedling growth). As a further sensitivity analysis, we calculated ecosystem multifunctionality using Random Forest R 2 values as weights.

We also calculated multifunctionality according to various management objectives, following Allan et al. (2015). In these measures, we gave different weightings to the various ecosystem functions, according to their presumed importance (based on a consensus of expert opinions of all authors) for delivering the ecosystem services required for the given objective (Fig. 1). The equal weights measure described above corresponds with most previous studies (e.g. [START_REF] Lefcheck | Biodiversity enhances ecosystem multifunctionality across trophic levels and habitats[END_REF]). In the measures representing management objectives, functions were weighted with loadings ranging from 0 (unimportant) to 1 (high importance). Functions directly related to the objective received a weight of 1, i.e. timber production and quality for 'timber production multifunctionality', carbon sequestration-related functions for 'climate regulation' and functions directly measuring biodiversity (e.g. bird/understory diversity) for 'biodiversity conservation/recreation'. Other functions were weighted 0.25; 0.50 or 0.75, depending on their relevance (Fig. 1). We also quantified a 'narrow-sense' biodiversity conservation measure, where only functions directly measuring biodiversity were included, with weights of 1 (Fig. 1).

Relationships between multifunctionality measures can either be caused by large-scale climatic/biogeographical factors (e.g. temperature gradients) or local-scale factors (e.g. management, soil conditions). Therefore, using Pearson correlations, we tested for trade-offs and synergies, at both continental (all plots) and local scales (within localities with >10 plots). With t-tests we investigated whether local-scale correlations, differed from zero.

Several functions had high weights in multiple multifunctionality measures, reflecting their relevance for different ecosystem services (Fig. 1B). Raw correlation coefficients between multifunctionality measures are therefore inflated by this overlap. To remove this effect, we calculated a null expectation for the correlation-coefficients by reshuffling ecosystem function values, without replacement, across plots 100 times. This eliminated any correlations among functions, while maintaining the original distribution of values. With these resampled ecosystem functions, we again calculated the different multifunctionality measures, and the average and 95% confidence intervals of the correlations between them. We calculated correlationcoefficients corrected for overlap in functions by subtracting expected values (in the absence of correlations among functions) from observed ones. As a sensitivity analysis, we repeated these analyses only including plots located within those 150 localities in which validations of both timber production and tree biomass were adequate (both r>0.1).

Comparing multifunctionality between protected versus non-protected forests

In total, 11.8% of the inventory plots were within protected areas which, depending on the NFI, indicated either that forestry activities were restricted (Germany, Sweden) or that the plot was in a National Park or nature reserve (Finland, France, Spain, Wallonia), see Supplementary Material for more detailed information. Within each country, we investigated, for each measure, whether local-scale multifunctionality was higher inside versus outside protected areas, using Welch's t-tests.

RESULTS

Explaining variation in ecosystem functioning

On average, across the different ecosystem functions in our fitting dataset, Random Forests explained 40.7% of the total variation. The explained variation in ecosystem functions ranged from high (timber production: 72.5%; resistance to insect herbivory: 67.6%) to low (browsing resistance: 2.4%, Fig. 1C). The single most important explanatory factor (i.e. with lowest residual variance) varied between the functions. For sixteen functions it was a climate variable, for six a functional identity variable, for two altitude, for two a functional diversity variable and soil pH and average stem diameter for one each (Fig. 1C; Table S4).

Three ecosystem functions allowed for validation of predicted values in inventory plots.

For timber production and tree biomass, across all plots, predicted values correlated reasonably well with observed values, with 'extrapolation' R 2 values (correlation between predicted and observed values in inventory plots) of 0.219 and 0.280, respectively. For tree recruitment the R 2 was only 0.040; Fig. S4. Validations generally worked best at large spatial scales and less well at local scales. Correlations between predicted and observed values of timber production, tree biomass and tree recruitment were, respectively, 0.390; 0.472 and 0.027 across 20×20km localities, and on average 0.127 (range: 0-0.976); 0.124 (range: 0-0.971) and 0.091 (range: 0-0.967) within localities. Absolute values of tree biomass were similar between NFI observations and Random Forest predictions, but for tree recruitment the values differed (Fig. S5). For more information on model validations, see Supplementary Material (S3).

Levels of ecosystem functioning and multifunctionality throughout Western Europe

After removing ecosystem functions poorly explained by the Random Forests (R 2 <0.2; see Fig. 1C), we predicted levels of 22 ecosystem functions for the inventory plots (Fig. S6).

Many of the mapped functions showed clear continental trends. For example, some (e.g. timber production) had highest levels in central Western Europe, while others had highest values in boreal (e.g. timber quality) or Mediterranean (e.g. bat diversity) regions (Fig. S6; Table S5).

Most functions tended to be highest at mid-latitudes. Consequently, most continental-scale multifunctionality measures were highest in central Western Europe (multifunctionality hotspots) and lowest in southern Europe (Fig. 2). When only diversity measures were considered (narrowsense biodiversity conservation), multifunctionality was also high in southern/central Spain and parts of Scandinavia. These patterns were broadly similar when functions with a high proportion of explained variance were weighted more heavily (Fig. S7). As expected, local-scale multifunctionality values did not show any large-scale spatial patterns (Fig. S8). Local multifunctionality scores were on average 45.8%, 47.1%, 49.2%, 49.8% and 47.8% below their maximum possible score (i.e. all functions above the 50% threshold) in the timber production, climate regulation, broad-sense and narrow-sense biodiversity conservation and overall multifunctionality scenario, respectively, and higher than 90% of the maximum possible score in 97, 49, 49 and 11,625 plots (out of 105,316 plots) in the timber production, climate regulation, broad-sense and narrow-sense biodiversity conservation scenario, respectively, whereas it exceeded 90% and 80% of maximum overall multifunctionality in only 3 and 446 plots respectively (Fig. 2B). Importantly, while ecosystem functions varied strongly at the continental scale (with 97.5 percentile values being on average 42.8% higher than mean values), there was also substantial variation within localities, with 97.5 percentile values being on average 12.6% higher than mean values (Table S6).

Trade-offs and synergies

Pairwise correlations between individual functions were positive on average at both scales, although correlations were weaker at local (𝑟 = 0.012) than at continental scales (𝑟 = 0.021), probably due to lower variation in functioning within localities (Table S6). Moderately to strongly positive correlations (r>0.3; n = 57 (continental-scale) and 22 (local scale)) outnumbered negative (r< -0.3; n = 45 (continental-scale) and 14 (local-scale)) correlations (Table S7,8). At the continental scale, correlations between timber production and tree biomass were similar for observed (r = 0.55) and extrapolated (r = 0.65) values. However, within localities this match was weaker (𝑟 = 0.63 observed and 0.24 predicted), with fits generally best in France and central/southern Spain, and weaker in Germany and northeast Spain (Fig. S9).

As different multifunctionality variables had similar continental-scale patterns (Fig. 2), continental-scale correlations between most measures were positive (Table 1). Only correlations between narrow-sense biodiversity conservation and both timber production (r = -0.13) and climate regulation multifunctionality (r = 0.01) were not. These correlations became more positive at more extreme (25 and 90%) multifunctionality thresholds (Table S9-S11).

Within localities, similar patterns were found. Relationships between timber production, climate regulation and broad-sense biodiversity conservation/recreation were positive, whereas relationships between narrow-sense biodiversity conservation and other multifunctionalty variables were close to zero, or negative, on average (Fig. 3, Table 1). Negative relationships largely disappeared when multifunctionality was based on 25% or 90% thresholds (Table S9-S11). Importantly, positive relationships between timber production and climate regulation multifunctionality, and to a lesser extent between timber production/climate regulation multifunctionality and broad-sense biodiversity conservation/recreation multifunctionality, were very widespread across Europe (Fig. 3).

We used null models to investigate whether observed correlations between multifunctionality variables were larger than expected. Relationships between multifunctionality variables were to a large extent driven by functions contributing to multiple multifunctionality variables, as observed minus expected correlation-coefficients were often close to zero (Fig. 3, Table 1). Nevertheless, at a continental scale, relationships between timber production, climate regulation and broad-sense biodiversity conservation multifunctionality remained significantly positive (all P<0.05). At the local scale, relationships between timber production and climate regulation multifunctionality also remained significantly (although weakly) positive, whereas relationships between timber production and the biodiversity conservation measures became significantly, weakly, negative. In sensitivity analyses these patterns hardly changed when (i) recruitment-related functions were omitted from multifunctionality measures, (ii) ecosystem functions with a high Random Forest fit had proportionally higher loadings in multifunctionality measures, or (iii) only plots from localities with high validation R 2 values of Random Forests explaining timber production and tree biomass were included (Table 1). Negative relationships largely disappeared when multifunctionality was quantified based on 25% or 90% thresholds (Table S9-S11). Importantly, functional overlap-corrected correlation-coefficients between different ecosystem multifunctionality scenarios varied greatly, from positive to negative, throughout localities (Fig. 3).

Multifunctionality inside versus outside protected areas

Local-scale associations between values of multifunctionality and protection status differed widely between countries and scenarios (Fig. 4). In Spain and Germany, timber production and climate regulation multifunctionality were lower inside protected areas, whereas the opposite was observed in France. In Germany, biodiversity conservation-related multifunctionality was highest inside protected areas, whereas in France the opposite was found.

These results were largely insensitive to the way in which multifunctionality was quantified (Table S12).

DISCUSSION

In our study trade-offs between groups of functions were rare in European forests, at both continental and local scales. We found synergies between individual ecosystem functions and few trade-offs between multifunctionality measures focused on timber production, climate regulation and biodiversity conservation/recreation. When corrected for overlap in functions among scenarios, some relationships were weakly positive throughout most of Europe (timber production versus climate regulation), some were weakly negative (timber production versus biodiversity conservation/recreation) and some were close to zero (climate regulation versus biodiversity conservation/recreation).The lack of strong trade-offs indicates that functions related to (sustainable) timber production can go hand in hand with functions related to services such as biodiversity conservation. Mapping local trade-offs and synergies across Europe revealed substantial variation in these relationships, showing that strong synergies are realized in a few environments. While biodiversity and timber production are currently maximised in some forests, suggesting a "win-win" for conservation and commercial forestry, across plots, average multifunctionality values were almost 50% below maximum possible levels, and the proportion of forest plots providing high levels of 'overall multifunctionality' (where timber production, climate regulation and biodiversity conservation are all maximized) was very small. Hence, while forest management has the potential to realize high multifunctionality, this is currently not common. Most multifunctionality measures had many ecosystem functions in common, as some ecosystem functions are valued under a range of different management objectives (e.g. [START_REF] Chan | Conservation planning for ecosystem services[END_REF]Allan et al. 2015). Relationships between different multifunctionality measures were generally much more strongly positive if not corrected for this functional overlap. While these raw correlations are statistically spurious (as the different measures partly contain the same data), they can be highly relevant for management. For instance, tree growth is important for both timber production and climate regulation, which suggests that forest management promoting tree growth will maximize both services. Our results therefore suggest many possibilities for win-win forest management strategies.

Our multifunctionality variables were intended to represent the bundle of functions needed to meet certain forest management objectives (following Allan et al. 2015). They should therefore be more useful to managers than traditional multifunctionality metrics that assume equal importance of each ecosystem function. However, they could be further improved to consider how multiple functions are related to final ecosystem services, using production functions, and then services can be valued in monetary or other units to calculate the overall benefits supplied by different management scenarios (e.g. Nelson et al. 2009;[START_REF] Bateman | Bringing economic services into economic decision-making: land-use in the United Kingdom[END_REF]. Ultimately, sustainable ecosystem management needs to minimize trade-offs between ecosystem benefits for different stakeholders [START_REF] Díaz | The IPBES conceptual framework -connecting nature and people[END_REF] and our targeted multifunctionality metrics represent a step towards quantifying and mapping these trade-offs at large scales.

Other studies, performed in grasslands (e.g. [START_REF] Lavorel | Using plant functional traits to understand the landscape distribution of multiple ecosystem services[END_REF] or across different ecosystems or land-use types [START_REF] Chan | Conservation planning for ecosystem services[END_REF]) have documented strong trade-offs between ecosystem functions and services, especially between productivity-related functions and those associated with biodiversity conservation or recreation. However, in forests, relationships between tree biomass and the biodiversity of associated taxa often show more mixed patterns [START_REF] Jukes | The influence of stand structure and composition on diversity of canopy Coleoptera in coniferous plantations in Britain[END_REF]. For example, the positive relationship between tree productivity and bird diversity in our data could be due to the strong dependence of specialist species on forests with many old trees (Gil-Tena et al. 2007), while the trade-off between productivity and understorey biomass may be driven by light competition between trees and understorey plants. When biodiversity conservation multifunctionality was quantified using only the four direct measures of biodiversity, weakly negative relationships with timber production and climate regulation multifunctionality were found. Their approximately equal strength at continental and local scales (Table 1) suggests that the relationship was primarily driven by local-scale factors, such as stand composition. The negative response of understorey plants to tree growth is likely responsible for this trade-off, as it is difficult to maximize timber production whilst maintaining an open canopy.

We also found that protected forests were not necessarily associated with high local-scale ecosystem multifunctionality. In Spain, several multifunctionality measures were in fact lower inside protected areas. In other countries, patterns were more mixed, but overall multifunctionality was never highest inside protected areas. Importantly, associations between forest protection status and multifunctionality were unlikely to be driven by climate, as localscale climatic variation is low within our 20×20km regions. Associations between local-scale multifunctionality and protection status seem therefore to be driven by local factors, such as tree diversity or composition. However, it is uncertain whether these observed relationships are causal, as forests were likely not designated to be protected at random. For example, they may have had low productivity and particular tree compositions before they were protected.

Furthermore, services such as the conservation of forest specialist species were not quantified, but these could be high inside protected areas. Many protected areas were only established relatively recently [START_REF] Paillet | Quantifying the recovery of old-growth attributes in forest reserves: a first reference for France[END_REF], so protected forests may still be recovering from past management. Finally, we only investigated forests without evidence of recent logging activity, which may have reduced the contrast between protected and non-protected areas. Regardless, although our results suggest a high potential for win-win forest management scenarios, the simultaneous maximization of timber production, climate regulation and biodiversity has not yet been realized within protected areas.

Our results also provide evidence that climate drives large-scale variation in many ecosystem functions and the synergies between them. Many functions, such as tree biomass or litter production, had highest levels in central Western Europe (Fig. S6) and some synergies between multifunctionality scenarios were stronger at continental than at local scales. A strong continental-scale synergy between earthworm biomass and litter decomposition (Table S7) may have arisen because they were both strongly associated with climate (Table S4). The correlation was also present at the local scale (Table S8), suggesting additional direct links between them.

While earlier studies have already shown the importance of climate for functions such as primary production and carbon sequestration (e.g. [START_REF] Cramer | Global response of the terrestrial ecosystem structure and function to CO2 and climate change: results from six dynamic global vegetation models[END_REF], our more comprehensive study shows that climate may be a driver of many more ecosystem functions, such as earthworm or microbial biomass. The fact that so many functions appear related to climate, especially to wet season precipitation (Table S4), may have important implications. For example, timber production multifunctionality was lower in dry climates, suggesting detrimental effects of projected future decreases in precipitation (IPCC 2014). However, while our approach is powerful in describing patterns, it is not suited to identify underlying processes. Therefore, more research on the causality of climate-ecosystem functioning relationships (e.g. [START_REF] Elith | Species distribution models: ecological explanation and prediction across space and time[END_REF], have a much stronger tradition in this respect. Three ecosystem functions could be validated with independent observations, which showed that: (1) validations were generally adequate for timber production and tree biomass, but not for tree recruitment, (2) validations worked best at large spatial scales, whereas at local scales there was large variation in their accuracy but (3) relationships between different multifunctionality variables were insensitive to the inclusion of localities where the validation was less well supported. Our approach is therefore promising, but we emphasize that validations could only be carried out for those three ecosystem functions for which independent inventory data was available, so future validations of other functions are needed. Local-scale data related to soil fertility or management could thus further improve the accuracy of ecosystem function predictions.

Our study presents a new approach to quantify ecosystem functioning at scales relevant for policy makers. The increasing availability of large datasets on ecosystem functioning from integrated projects means our approach may become increasingly feasible for other systems and regions. A further possibility would be to combine local-scale ecosystem functioning datasets with remote sensing data to map services at large scales. Remote sensing approaches have successfully predicted some ecosystem functions, but have difficulties with other functions, such as soil processes (de [START_REF] De Araujo | Remote sensing of ecosystem services: a systematic review[END_REF]. By combining data on forest and climate attributes with remotely sensed parameters, we could map ecosystem functions even more accurately in the future. Our study is a first step in reaching the ultimate goal of predicting how future ecosystem functioning and service provision will be altered by ongoing global trends, such as climate change (IPCC 2014), eutrophication and acidification [START_REF] Galloway | Transformation of the nitrogen cycle: Recent trends, questions, and potential solutions[END_REF] or land-use change [START_REF] Newbold | Global effects of land use on local terrestrial biodiversity[END_REF]. Future studies could combine our approach with models on climate change (e.g. IPCC scenarios), biodiversity change (e.g. [START_REF] Isbell | The biodiversity-dependent ecosystem service debt[END_REF] or management scenarios to investigate the impacts of these global trends for the future functioning and service provisioning of forests and other ecosystems.

In conclusion, our study, among most comprehensive overviews of forest ecosystem functioning to date, showed that different measures of forest multifunctionality tend not to tradeoff with each other, at both local and continental scales. Within some areas there were strong synergies between different multifunctionality measures, indicating that even though they are currently uncommon, "win-win" forest management strategies are possible and could be promoted in the future. However, we also found that multifunctionality is often not higher inside than outside protected areas. Our study therefore suggests a high but unrealized potential for multifunctionality in European forests.

RUNNING TITLE:

  continental mapping of multifunctionality KEYWORDS: biodiversity, climate, ecosystem multifunctionality, ecosystem services, forest, FunDivEUROPE, large-scale, phylogenetic diversity, tree communities, upscaling TYPE OF ARTICLE: letter

  Studies have used exhaustive remote sensing or ground-based measurements (e.g. Prince & Goward 1995; Ratcliffe et al. 2016), mechanistic models (e.g. McGuire et al. 2001), indirect measures (e.g. where certain habitat types are assumed to promote certain functions; Maskell et al. 2013) or a combination of these (Maes et al. 2012; Mouchet et al. 2017) to quantify single or multiple functions at large spatial extents. However, for some important functions, such as biological pest control or timber quality, large scale maps have not yet been developed, limiting our understanding of ecosystem functioning synergies and trade-offs. In contrast, many localscale studies, such as biodiversity experiments (e.g. Hector & Bagchi 2007; Zavaleta et al. 2009) or comparative studies (Lavorel et al. 2011), have accurately quantified a large number of functions. Extrapolating these small-scale observations to larger scales could increase our understanding of the drivers of ecosystem functioning trade-offs and the resulting provision of ecosystem multifunctionality.

FIGURESFigure 1 .

 1 FIGURES

  importance for the management objective. Note that in the climate regulation scenario, loadings of the decomposition variables are negative. C: Proportion of variance of ecosystem functions explained by Random Forests. Letters above the bars indicate which type of predictor was most important in explaining variation: C = climate-related; I = functional identity-related; P = pH; A = altitude; D = biodiversity-related; S = stand structure related. In further analyses, only those functions with R 2 values above 0.2 (dashed horizontal line) were included.

Figure 2 . 8 .FrequencyFigure 3 . 12 *

 28312 Figure 2. While high values of continental-scale multifunctionality (A, C-F) in central Europe

Figure 4 .

 4 Figure 4. Local-scale ecosystem multifunctionality is generally not higher inside protected areas,

  al.

	2001), soil pH (Foy 1992), functional community composition (Diaz et al. 2004) and tree
	diversity (van der Plas et al. 2016; Liang et al. 2016) can all drive (forest) functioning.

  De Boeck et al. 2008; Šímová & Storch 2017) is needed to predict ecosystem responses to climate change. Extrapolations are still relatively rare in ecosystem functioning studies (but see Lee et al.

	2000; Isbell et al. 2014; Manning et al. 2015), although other subtopics of ecology, such as
	species distribution research

ACKNOWLEDGEMENTS

This paper is a joint effort of the working group 'Scaling biodiversity-ecosystem functioning relations: a synthesis based on the FunDivEUROPE research platforms' on the 24 th -26 th November 2014 in Leipzig, Germany, kindly supported by sDiv, the Synthesis Centre of the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, funded by the German Research Foundation (FZT 118). The FunDivEUROPE project received funding from the European Union's Seventh Programme (FP7/2007-2013) under grant agreement No. 265171.

We thank the MAGRAMA for access to the Spanish Forest Inventory, the Johann Heinrich von Thünen-Institut for access to the German National Forest Inventories, the Natural Resources Institute Finland (LUKE) for making the Finnish NFI data available, the Swedish University of Agricultural Sciences for making the Swedish NFI data available, and Hugues Lecomte, from the Walloon Forest Inventory, for access to the Walloon NFI data. The study was supported by the TRY initiative on plant traits (http://www.trydb.org). The TRY initiative and database is hosted, developed and maintained at the Max Planck Institute for Biogeochemistry, Jena, Germany.

TRY is/has been supported by DIVERSITAS, IGBP, the Global Land Project, the UK Natural Environment Research Council (NERC) through its program QUEST (Quantifying and Understanding the Earth System), the French Foundation for Biodiversity Research (FRB), and GIS "Climat, Environnement et Société" France.

AUTHOR CONTRIBUTIONS

FvdP, EA, LG, MF, SR, PRB, MSL, CW, LB, RB, HB, TJ, SK, GK, CN, BO, AP and FR developed the ideas of this study at a workshop organized by EA and LG. FvdP, SR and PRB analysed the data. All authors, except FvdP, EA, MF, SKa, PM, BO, AP and FR contributed to the data collection. FvdP wrote the manuscript. All authors contributed in editing the manuscript. Table 1. Correlations between values of different multifunctionality measures at both continental and local scales and both across all plots and within countries. Here, multifunctionality was based on a 50% threshold level. Correlations were also quantified after correcting for the overlap in ecosystem functions between multifunctionality measures. This is indicated as 'no functional overlap' or 'no FO' in the table. As sensitivity analyses, correlations were also calculated based on (a) multifunctionality measures in which recruitment-related functions were excluded, (b) multifunctionality measures in which loadings of ecosystem functions was proportional to Random Forest R 2 values and (c) only those plots within 20x20 km grid cells with a high validation R 2 (>0.10) for timber production and tree biomass. Significant correlations are shown in bold. TP = timber production, CR = climate regulation, BCB = broad-sense biodiversity conservation and BCN = narrow-sense biodiversity conservation.


E Menier

M A Bucci

M Yagoubi

L Mathelin

M Schoenauer

CD-ROM Complementary Deep -Reduced Order Model

Proper Orthogonal Decomposition can be used to extract the principal modes (V) from a physical simulation :

The flow approximates as a linear combination of a reduced number (r) of modes :

An equation for the coefficients of this simplified form is obtained through Galerkin Projection of the original system :

∂u(x, t) ∂t = g (u(x, t)), V T u = α =⇒ dα dt = V T g (u(x, t))
The dynamics g () of the reduced model are computed from the approximate solution, which introduces error :

dα dt ≈ V T g ( ū) = V T g (u(x, t))

ROM Correction

To address classical POD-Galerkin model's shortcomings, we propose to add a correction term to their dynamics :

dα dt = V T g (u) = V T g ( ū) + f
This correction term depends on information outside of the POD basis:

Corrected ROM:

g (u) = g (VV T u + (I -VV T )u) ≈ g (V α) + g V α (I -VV T )u
Following Takens's theorem, we can retrieve the information lost during the projection on the POD basis by considering the past states of the system :

DDE

Delay Differential Equations can be used to aggregate information from the past in a time continuous manner :

dx dt = f (t, x, y ), y (t) = t -∞ x(τ )e λ(τ -t) dτ, λ ∈ R +
These equations are solved as an augmented ODE system :

dx dt dy dt = f (t, x, y ) x -λy
We use Neural Networks to learn the ROM's correction from the memory y :

dα dt dy dt = V T g (V α) + N N (y ) α - λy
The information available in memory y is limited by its dimension :

y (t) = t -∞ α(τ )e λ(τ -t) dτ =⇒ dim(y ) = dim(α)
To address the bottleneck in the memory, we propose to add an encoder network to expand the dimension of the memory :

dα dt dy dt = V T g (V α) + N N (y ) E(α) - λy
Leading to the final CD-ROM formulation :

Training data for the true ROM trajectories can be obtained from the snapshot data :

αt i = V T u DNS,t i
The model can then be trained through adjoint backpropagation2 :

L(α) = tn t 0 α t i -αt i 2
The test trajectory presents statistics similar to the original attractor :

Probability density function of the amplitude of each mode over a trajectory of ≈ 1000 seconds. The CD-ROM results are presented in blue while the DNS is in black.

The test trajectory presents statistics similar to the original attractor :

Time Horizons

We can learn a different λ for each dimension of the memory :

y (t) = t -∞ E (x(s))e Λ(s-t) ds, Λ =    λ 1 . . . λ N    , N = dim(y )
Each λ i defines the time horizon of each memory dimension, which can be computed as :

τ i = 1 λ i
Sine waves are a good test-case to assess the ability of the CD-ROM to learn the time-horizons : With a wave of period T , using a delay of τ = T 4 yields a perfect circle in phase space.

B. R. Noack et al., "A hierarchy of low-dimensional models for the transient and post-transient cylinder wake", Journal of Fluid Mechanics 497, 335-363 (2003).

R. T. Q. Chen et al., "Neural ordinary differential equations", (2019).

N. Deng et al., "Low-order model for successive bifurcations of the fluidic pinball", Journal of Fluid Mechanics 884, 10.1017/jfm.2019.959 (2019).

Pinball Results

R. T. Q. Chen et al., "Neural ordinary differential equations", (2019).

Cylinder Results

/ 21

Chaotic Pinball

We tested the approach on a more challenging case, the chaotic pinball 3 : Reminder : 8 modes are enough to capture 99 % of the cylinder flow's energy.

Unstable ROM

Taking a very reduced number (r = 10) of modes for the pinball flow yields a very unstable model.

Trajectory Fit

We use a superposition of two waves of periods T 1 = 50, T 2 = 168 :

We need a memory of dimension 3 to represent the dynamics of the wave.

Learned Horizons

Starting from varying initial conditions, the model learns similar horizons.

Summary & Future Work

Summary :

• We propose a novel, time-continuous, neural architecture for modeling partially observable systems.

• It provides a computationally inexpensive and data-efficient correction to POD-Galerkin models, illustrated here in different flow configurations.

• Relevant and physically consistent time-constants are extracted from observations of the system.

Computational Costs

Although the CD-ROM increases the computational costs, it remains less expensive than bigger models :

Pinball Training

This model can then be progressively corrected and stabilised through our approach :

Comparison with Reservoir Computing Echo State Networks are based on random matrices W in and W :

Only the output operator is optimised to map the output sequence (O) to the memories :

Comparison with RC

Ignoring the non-linearity, the ESN becomes a linear dynamical system :

The operator W = W -I ∆t eigenvalues have negative real part. It plays the same memory dissipation role as the Λ matrix in our architecture :

Neural ODE

With Neural ODEs 4 , we can back-propagate gradients through the simulation steps :

ROM NODE

In our case, the Neural ODE corresponds to the complete ROM :

Because the model is fully expressed with differentiable objects, the required jacobians can be easily evaluated through automatic differentiation.

Cylinder Hyperparameters