Joachim Steenbeek

J Claudet

A Lloret

X Corrales

; Vilas

; Piroddi

; Steenbeek

; Claudet

; Lloret

; Calò

Di Franco

; Font

; Ligas

; Prato

; Sahyoun

; Sartor

; Guidetti

; Coll

To be re-submitted to Ocean and Coastal Management Multi-zone marine protected areas: assessment of ecosystem and fisheries benefits using multiple ecosystem models

Keywords: Northwestern Mediterranean Sea, Ecopath model, Marine Protected Areas, Small-scale fisheries, Recreational fisheries, Ecological indicators

Introduction

Coastal and marine ecosystems provide multiple ecosystem goods and services [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF][START_REF] Martínez | The coasts of our world: Ecological, economic and social importance[END_REF]). However, the provision of these ecosystem goods and services can be degraded as human demand for resources and the associated activities increase in number, frequency, magnitude and spatial coverage (MEA 2005;[START_REF] Worm | Impacts of biodiversity loss on ocean ecosystem services[END_REF][START_REF] Pauly | Catch reconstructions reveal that global marine fisheries catches are higher than reported and declining[END_REF]. Currently, marine ecosystems are under increasing threat from a range of stressors including overfishing, climate change, biological invasions, pollution, aquaculture and habitat degradation, directly or indirectly caused by multiple human activities [START_REF] Costello | A census of marine biodiversity knowledge, resources, and future challenges[END_REF][START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF]. The relationships and interactive effects of multiple stressors are highly complex, and largely unknown or even highly uncertain [START_REF] Crain | Understanding and managing human threats to the coastal marine environment[END_REF][START_REF] Côté | Interactions among ecosystem stressors and their importance in conservation[END_REF]).

In recognition of the complex interrelationship between marine organisms, multiple human activities, and ecosystem functioning and services, more comprehensive frameworks to manage marine ecosystems are required, such as the ecosystem-based management (EBM) [START_REF] Leslie | Confronting the challenges of implementing marine ecosystem-based management[END_REF][START_REF] Long | Key principles of marine ecosystem-based management[END_REF]. The overall objective of this approach is to maintain ecosystems as a whole "in healthy, productive and resilient conditions allowing them to provide the needed ecosystem services to society" [START_REF] Mcleod | Scientific consensus statement on marine ecosystembased management[END_REF]. Within this context, Marine Protected Areas (MPAs hereafter) are considered to be a key management tool for EBM as they are expected to mitigate human impacts, recover exploited resources and support their sustainable use, conserve or restore habitats and biodiversity, maintain and enhance ecosystem services and reduce conflicts between users [START_REF] Lester | Biological effects within no-take marine reserves: a global synthesis[END_REF][START_REF] Halpern | Placing marine protected areas onto the ecosystembased management seascape[END_REF][START_REF] Leenhardt | The role of marine protected areas in providing ecosystem services[END_REF]. A MPA is a discrete area of the sea established to achieve the long-term conservation of natural resources therein [START_REF] Claudet | Marine protected areas: a multidisciplinary approach[END_REF]. The level of protection within MPAs can vary from areas of full protection (FPAs, also called marine reserves or no-take areas), where all extractive activities (e.g., fishing) are prohibited while some nonextractive activities (e.g. diving) can be allowed; to areas of partial protection (PPAs),

where some human activities (e.g. small-scale and recreational fisheries) are allowed but strictly regulated (Horta e Costa, Claudet et al. 2016). Conservation benefits of MPAs largely vary due to their intrinsic features, including their design, the level of enforcement, their conservation goals and organization [START_REF] Guidetti | Italian marine reserve effectiveness: does enforcement matter?[END_REF][START_REF] Franco | Five key attributes can increase marine protected areas performance for small-scale fisheries management[END_REF][START_REF] Giakoumi | Revisiting "success" and "failure" of marine protected areas: a conservation scientist perspective[END_REF][START_REF] Scianna | Organization Science improves management effectiveness of Marine Protected Areas[END_REF]). In addition, networks of MPAs could provide greater benefits than the sum of individual MPAs benefits [START_REF] Gaines | Designing marine reserve networks for both conservation and fisheries management[END_REF][START_REF] Grorud-Colvert | Marine protected area networks: assessing whether the whole is greater than the sum of its parts[END_REF]. A network of MPAs is a system of individual MPAs that operates synergistically to achieve ecological objectives more effectively than individual MPAs could fulfil alone (IUCN-WCPA 2008).

Within the EBM, ecological modelling has been proven to be a suitable tool to merge available ecological information into a coherent form to obtain insights about how ecosystems are structured and impacted by human activities and environment, and delivering ecosystem services (Link 2010;[START_REF] Christensen | Ecosystem approaches to fisheries: a global perspective[END_REF]. Although quantitative models have become a key tool to assess the extent to which MPAs are able to achieve their conservation targets [START_REF] Fulton | Modelling marine protected areas: insights and hurdles[END_REF], the understanding of MPA effects on ecological processes and ecosystem function is still scarce [START_REF] Cheng | Can marine reserves restore lost ecosystem functioning? A global synthesis[END_REF]. For instance, it is mostly unknown if protection benefits inside fully protected areas extend to partially areas and beyond, and if those benefits have the potential to cascade through the ecosystem components. There is little evidence that multi-zone MPAs (i.e., those including both FPAs and PPAs) provide fisheries benefits inside partially protected areas and it is unknown what are the drivers of multi-zone MPAs.

Although considered a hotspot of biodiversity [START_REF] Coll | The biodiversity of the Mediterranean Sea: estimates, patterns, and threats[END_REF]), the Mediterranean Sea has been largely altered by anthropogenic stressors for a long time [START_REF] Coll | The Mediterranean Sea under siege: spatial overlap between marine biodiversity, cumulative threats and marine reserves[END_REF][START_REF] Micheli | Cumulative Human Impacts on Mediterranean and Black Sea Marine Ecosystems: Assessing Current Pressures and Opportunities[END_REF][START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF], that have determined dramatic changes on its ecosystems especially during the last 50 years [START_REF] Azzurro | Tracking signals of change in Mediterranean fish diversity based on local ecological knowledge[END_REF][START_REF] Maynou | Estimating trends of population decline in long-lived marine species in the Mediterranean Sea based on fishers' perceptions[END_REF][START_REF] Piroddi | Historical changes of the Mediterranean Sea ecosystem: modelling the role and impact of primary productivity and fisheries changes over time[END_REF]. In order to address this urgent conservation challenge, a multitude of MPAs (including all forms of nationally and internationally established protected areas) have been implemented in the Mediterranean Sea in the last decades. Within this context, several studies have modelled MPAs in the Mediterranean Sea with the aim to evaluate their effects on the ecosystem and the impacts of small-scale and recreational fisheries on marine resources [START_REF] Albouy | Simulation of the combined effects of artisanal and recreational fisheries on a Mediterranean MPA ecosystem using a trophic model[END_REF][START_REF] Libralato | Food-web traits of protected and exploited areas of the Adriatic Sea[END_REF][START_REF] Valls | Modeling trophic interactions to assess the effects of a marine protected area: case study in the NW Mediterranean Sea[END_REF][START_REF] Prato | Assessing interacting impacts of artisanal and recreational fisheries in a small Marine Protected Area (Portofino, NW Mediterranean Sea)[END_REF]Vilas, Coll et al. Submitted). However, none of the existing studies modelled simultaneously the multizone nature of Mediterranean MPAs nor a series of neighbouring MPAs, potentially constituting a network of protected areas. In addition, the study of the historical effects of MPAs has not been pursued due to difficulties to gather historical data about protection effects. In this study we modelled three relatively well studied MPAs from the Northwestern Mediterranean Sea with a historical perspective: Cerbère-Banyuls, Cap de Creus and Medes Islands. Our aims were to: (1) quantify the main structural and functional traits of the three multi-zone MPAs, and (2) assess how exploitation and protection regimes affected these marine ecosystems over time, using a comparative approach to identify commonalities and differences between MPAs.

Materials & Methods

Study area

We studied three MPAs in the Northwestern Mediterranean Sea: Cerbère-Banyuls (France), Cap de Creus (Spain) and Medes Islands (Spain) (Fig. 1). These MPAs include one Fully Protected Area (FPA) and one (Cerbère-Banyuls and Medes Islands) or more (Cap de Creus) Partially Protected Areas (PPAs) (Fig. 1). While in the FPAs all extractive activities are forbidden, in PPAs both small-scale and recreational fisheries are allowed except in Medes Islands MPA, where recreational fisheries are prohibited.

The three MPAs comprise a similar depth range between ca. 0 and 60 meters, and several vulnerable species such as the Mediterranean seagrass (Posidonia oceanica), dusky grouper (Epinephelus marginatus) and red coral (Corallium rubrum), can be found in these MPAs.

In addition to the protected zone scheme (FPA and PPA), we included in our study the immediate unprotected areas (UPAs) surrounding each MPA (Fig. 1 and Table 1), in order to add the adjacent non-protected waters where some species found inside the MPAs can theoretically forage and move to (Di [START_REF] Lorenzo | Spillover from marine protected areas to adjacent fisheries has an ecological and a fishery component[END_REF]. The boundaries of the unprotected areas were selected taking into account that they had similar ecological and oceanographic features (i.e., same depth range and habitats) to the MPA, and they constituted part of the Parc (Cap de Creus and Medes Islands) and/or were adjacent to the MPA (Cerbère-Banyuls) (Fig. 1).

Modelling approach -Ecopath and Ecosim

The three MPA models were developed using Ecopath with Ecosim approach (EwE 6.6 version) [START_REF] Christensen | Ecopath with Ecosim: methods, capabilities and limitations[END_REF][START_REF] Christensen | Ecopath with Ecosim version 6 user guide[END_REF]. We used the static Ecopath model to provide a "snapshot" of the ecosystem in terms of flows and biomasses for a defined period of time. To develop the food-web model we used functional groups, which consist of ontogenetic fractions of a species, single species or groups of species sharing common ecological traits. The baseline models of each MPA were developed taking into account the contribution of each zone to the whole MPA based on a previous study (Vilas, Coll et al. Submitted). Afterwards, we used the time-dynamic model Ecosim, which describes the temporal dynamics of species biomass and flows over time by accounting for changes in predation, consumption rate, fishing and the environment [START_REF] Walters | Structuring dynamic models of exploited ecosystems from trophic mass-balance assessments[END_REF][START_REF] Christensen | Ecopath with Ecosim: methods, capabilities and limitations[END_REF]. A short description of the modelling approach is given in the Supplementary Online Material (SOM hereafter) 1; and a detailed information of the EwE approach is documented in [START_REF] Christensen | Ecopath with Ecosim: methods, capabilities and limitations[END_REF] and [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF].

Functional groups and input data

Functional groups

A meta-web structure defined for the Western Mediterranean Sea model (Coll, Steenbeek et al. 2019)

Ecopath input data and balancing procedure

The Ecopath models represented a situation of the Cerbère-Banyuls MPA for 2013, while the Cap de Creus MPA and the Medes Islands MPA represented an average situation for the period 2005-2008 and 2000-2004, respectively. These periods of time to model the MPAs were selected considering the best available biomass data (SOM 2 for details on the parameterization of each functional group). Input parameters were obtained using similar procedures as those documented in [START_REF] Coll | Trophic flows, ecosystem structure and fishing impacts in the South Catalan Sea, Northwestern Mediterranean[END_REF], [START_REF] Corrales | Ecosystem structure and fishing impacts in the northwestern Mediterranean Sea using a food web model within a comparative approach[END_REF], [START_REF] Piroddi | Modelling the Mediterranean marine ecosystem as a whole: addressing the challenge of complexity[END_REF] and Vilas, Coll et al. (Submitted). Biomass and catch data were scaled using data from models of each management unit (FPA, PPA and UPA) developed previously [START_REF] Coll | Report describing the quantitative models[END_REF]Vilas, Coll et al. Submitted) by the percentage of the area of each management unit to the whole MPA (SOM 2 for details on parameterization of each functional group).

Production (P/B, year -1 ) and consumption (Q/B, year -1 ) rates were estimated through the application of empirical equations [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF], taken from literature or from other models developed in the Northwestern Mediterranean Sea [START_REF] Coll | Trophic flows, ecosystem structure and fishing impacts in the South Catalan Sea, Northwestern Mediterranean[END_REF][START_REF] Corrales | Ecosystem structure and fishing impacts in the northwestern Mediterranean Sea using a food web model within a comparative approach[END_REF][START_REF] Coll | Report describing the quantitative models[END_REF].

The trophic information to populate the diet matrix was compiled using published studies on stomach content analyses, giving preference to local or surrounding areas (SOM 2). We used the Diet Matrix Calculator [START_REF] Steenbeek | Diet Calculator -Quick Reference Guide[END_REF], a custom-built extraction tool that automatizes the process of selecting and scaling diet data. Drawing on a large library of published diet studies, the Diet matrix calculator tool selects the most likely suitable diet studies for a specific model area, based on a range of criteria, and generates a diet composition matrix with accompanying pedigree index for each predatory functional group.

Due to the small sizes of the investigated MPA and the capacity of some species to move greatly [START_REF] Gell | Benefits beyond boundaries: the fishery effects of marine reserves[END_REF][START_REF] Grüss | Consequences of adult and juvenile movement for marine protected areas[END_REF], we set a fraction of the diet composition as import for all MPA based on the time that these species spent foraging outside the areas and biological and ecological traits such as size, behavior and ecology of species of each functional group (Table 1 in SOM 3).

To achieve mass-balance of the baseline models, we implemented a manual procedure following a top-down approach (starting from the groups with higher trophic levels) modifying appropriate input parameters (Table 2 in SOM 3) and following the best practice guidelines provided in the literature [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF]). In addition, the PREBAL analysis was used to ensure that the input data complied with general ecologic principles and to guide modifications through the balancing procedure (Link 2010). Detailed explanations of these procedures are given in SOM 3 and the final diet matrixes are provided in SOM 4.

The pedigree routine [START_REF] Christensen | Ecopath with Ecosim: methods, capabilities and limitations[END_REF]) was used to evaluate the quality of the models and the uncertainty associated with the input parameters and to guide modifications in balancing the model. A detailed explanation of this routine is given in SOM 3.

Ecosim historical time series fitting.

The model representing Cerbère-Banyuls MPA ecosystem during the 2013-2017 period was fitted to time series of historical data, while for Cap de Creus MPA and the Medes Islands MPA, the Ecosim model was fitted to time series of 2008-2017 and 2004-2017 periods, respectively (Table 1 in SOM 5). We used historical fishing effort trends to drive the fisheries of the models, while biomass and catch time series were used to calibrate the model and compare predicted to observed results (see SOM 5 for details), respectively, following previous studies in the Mediterranean Sea and best practices of EwE [START_REF] Coll | Ecological role, fishing impact, and management options for the recovery of a Mediterranean endemic skate by means of food web models[END_REF][START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF][START_REF] Corrales | Hindcasting the dynamics of an Eastern Mediterranean marine ecosystem under the impacts of multiple stressors[END_REF].

The models were fitted by using the Stepwise Fitting Procedure [START_REF] Scott | A Stepwise Fitting Procedure for automated fitting of Ecopath with Ecosim models[END_REF].

This procedure tests alternative hypotheses related to the impact of fishing, changes in predator-prey dynamics (vulnerabilities), changes in primary production (production anomalies) or all of the above together (Table 2) [START_REF] Mackinson | Which forcing factors fit? Using ecosystem models to investigate the relative influence of fishing and changes in primary productivity on the dynamics of marine ecosystems[END_REF][START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF]. This procedure estimates different vulnerability parameters and primary productions anomalies to improve model fits by comparing model predictions to observed data using the sum of squares (SS) statistics and to find the statistically "best fit" model based on Akaike's Information Criterion (AIC), which penalizes the overparameterization of the model [START_REF] Mackinson | Which forcing factors fit? Using ecosystem models to investigate the relative influence of fishing and changes in primary productivity on the dynamics of marine ecosystems[END_REF][START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF]. In this study, the maximum number of parameters that could be estimated were 19

(Cerbère-Banyuls), 51 (Cap de Creus) and 48 (Medes Islands), respectively.

Finally, we manually evaluated a range of models fits with the lowest AIC values to identify the best final model, which led to a credible behavior [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF],

following previous studies (e.g., [START_REF] Corrales | Hindcasting the dynamics of an Eastern Mediterranean marine ecosystem under the impacts of multiple stressors[END_REF].

Model analyses and ecological indicators

Ecological indicators of initial conditions

To analyze the food-web structure of the three MPAs, we used the biomasses of selected FGs, trophic flows and trophic levels (TLs) within the flow diagram. The TL was also used to analyze the ecological position of the FG of the three MPA models [START_REF] Lindeman | The trophic-dynamic aspect of ecology[END_REF][START_REF] Stergiou | Feeding habits and trophic levels of Mediterranean fish[END_REF] Several ecological indicators were also computed to describe the structure and functioning of the ecosystems and were divided into four groups using the ECOIND plug-in [START_REF] Coll | Standardized ecological indicators to assess aquatic food webs: The ECOIND software plug-in for Ecopath with Ecosim models[END_REF]):

(1) Biomass-based: calculated from the biomass of components included in the foodweb model, they could provide valuable information to evaluate MPAs effectiveness [START_REF] Micheli | Trajectories and correlates of community change in no-take marine reserves[END_REF][START_REF] Claudet | Marine reserves: size and age do matter[END_REF]. We included Total Biomass (TB, t•km -2 •year -1 ), Biomass of Fish species (FB, t•km -2 •year -1 ), and Kempton's Q diversity index (QI).

(2) Trophic-based: reflect the TLs for different groups of the food web, provide information on the structure of the ecosystem and are used to quantify the impact of fishing [START_REF] Rochet | Which community indicators can measure the impact of fishing? A review and proposals[END_REF]. We selected TL of the community (TLc), TL of the community, TL of the community including organisms with TL ≥ 3.25 (TL3.25) and TL of the community including organisms with TL ≥ 4 (TL4).

(3) Species and size-based: based on species traits and conservation status, they could offer insights of the effects of MPAs [START_REF] Claudet | Marine reserves: fish life history and ecological traits matter[END_REF]. We selected biomass of IUCN-endangered species biomass in the community (t•km -2 •year -1 ), mean length of fish in the community (ML, cm) and mean life span of fish community (MLS, year).

(4) Catch-based: based on catch, they reflect the fishing strategy of the fisheries and are used to quantify the impact of fishing [START_REF] Hilborn | Quantitative fisheries stock assessment: choice, dynamics and uncertainty[END_REF][START_REF] Pauly | Fishing down marine food webs[END_REF]). We included total catch (TC t•km -2 •year -1 ), trophic level of the catch (TLc), intrinsic vulnerability index of catch (VI).

The Mixed Trophic Impact (MTI) analysis was performed to assess the direct and indirect impact in the food web that a hypothetical increase in the biomass of one functional group would have on the biomasses of all the other functional groups in the food web, including the fishing fleets [START_REF] Ulanowicz | Mixed trophic impacts in ecosystems[END_REF][START_REF] Christensen | Ecopath with Ecosim version 6 user guide[END_REF]). To evaluate the impact of small-scale fisheries on the MPAs, the MTI was used to quantify the direct and indirect impact of each fishery on the functional groups for the studied MPAs, and their potential competition and trade-offs between them.

Finally, to identify the key species within the ecosystem (both keystone and structuring species), we calculated the keystoneness index developed by [START_REF] Valls | Keystone species: toward an operational concept for marine biodiversity conservation[END_REF] of the most important reef functional groups. A keystone group is defined as a predator species with a large and broad impact on the food web despite its low biomass [START_REF] Paine | Food web complexity and species diversity[END_REF][START_REF] Paine | A note on trophic complexity and community stability[END_REF][START_REF] Valls | Keystone species: toward an operational concept for marine biodiversity conservation[END_REF].

Time-dynamic analyses

Once the fitting procedure was completed, we used the best fit model to examine biomass and catch time series predicted by the model to explore the dynamics of selected functional groups. We selected four target species due to their role in terms of biomasses and commercial interest: the common two-banded seabream (Diplodus vulgaris), white seabream (D. sargus sargus), common dentex (Dentex dentex) and groupers (this group was mainly represented by Epinephelus marginatus). These species play an important role in the ecosystem as high trophic level predators or intermediate trophic species and are of great importance in small-scale and recreational fisheries. In addition, the previous selection of ecological indicators from the ECOIND plug-in was used to describe ecological changes in the ecosystem over time [START_REF] Coll | Standardized ecological indicators to assess aquatic food webs: The ECOIND software plug-in for Ecopath with Ecosim models[END_REF].

The Monte Carlo (MC) uncertainty routine and the ECOSAMPLER plug-in were employed to assess the impact of uncertainty in Ecopath input parameters on Ecosim simulations [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF][START_REF] Coll | Standardized ecological indicators to assess aquatic food webs: The ECOIND software plug-in for Ecopath with Ecosim models[END_REF][START_REF] Steenbeek | Ecosampler: A new approach to assessing parameter uncertainty in Ecopath with Ecosim[END_REF]). We ran 500 MC simulations based on the sensibility of Ecopath input parameters obtained from the pedigree routine (SOM 3). Results from the MC simulations were used to plot the confidence intervals of the selected ecological indicators in Ecopath and to plot the 5 th and 95 th percentile confidence intervals for the Ecosim outputs. Finally, we used the Spearman's rank correlation to evaluate the presence of significant trends (either increasing or decreasing) in model outputs (biomass, catch and ecological indicators), following previous studies (e.g., [START_REF] Corrales | Hindcasting the dynamics of an Eastern Mediterranean marine ecosystem under the impacts of multiple stressors[END_REF]).

Results

Structure and functional traits of protected areas

The pedigree index values of the three MPA models showed similar values, ranging from 0.41 to 0.51. The highest pedigree values were obtained for Cerbère-Banyuls (0.51), followed by Medes Islands (0.45) and Cap de Creus (0.41).

The visualization of trophic links and flows between functional groups highlighted the complexity of these coastal ecosystems due to the large number of trophic links between functional groups and the important role of detritus (FG 65) and other macro-benthos (FG 51) in transferring energy up to the food web (Fig. 2).

The functional groups of the models ranged from trophic level (TL) of 1 for primary producers (FG. 60-67) and detritus groups (FG 66-67) to TL = 4.2 for both groups of dolphins (FG 1-2) (Fig. 2 and Table 3 in SOM 3). Invertebrates groups were classified with a TL between 2 and 3.5, with benthopelagic cephalopods (FG 41) showing higher TLs. Fish had TLs between 3 and 4, with the exception of salema (FG 33) and mugilidae (FG 34), which showed lower TL due to their herbivorous and detritivore behaviors. Overall, similar

TLs were found for the three MPAs.

Total and fish biomass displayed similar patterns, with the highest biomass values found for Cerbère-Banyuls and Medes Islands and lower values in Cap de Creus (Fig. 3).

Conversely, Kempton's Q Index presented a higher value in Cap de Creus, followed by Medes Islands and Banyuls-Cerbère (Fig. 3). With the exception of TL of the community 4, trophic based indicators also presented higher values for Cerbère-Banyuls, followed by Medes Islands and Cap de Creus (Fig. 3). Species and size-based indicators showed that Cerbère-Banyuls presented the highest values for ML of fish community and IUCN species B, followed by Cap de Creus and Medes Islands (Fig. 3), while MLS of fish community was higher in Cap de Creus, followed by Cerbère-Banyuls and Medes Islands (Fig. 3).

Total catch showed similar values between Cap de Creus and Medes Islands, while Cerbère-Banyuls had the lowest value (Fig. 3). TL of the catch and the Intrinsic Vulnerability Index presented similar values in Cerbère-Banyuls and Cap de Creus, while

Medes Islands had the lowest value (Fig. 3).

The keystoneness index identified groupers (FG 27), common dentex (FG 24), other commercial medium demersal fishes (FG 31) and non-commercial medium demersal fishes (FG 32) as keystones species in the three MPAs, followed by red scorpionfish (Scorpaena scrofa) (FG 25) and Scorpaenidae (FG 26) (Fig. 4).

The MTI analysis highlighted that the small-scale fisheries had the most widespread negative impact on many FG of all MPAs, especially in Cerbère-Banyuls, while the impact of recreational fisheries was more prominent in Cap de Creus and Medes Islands (Fig. 5).

This analysis showed strong negative impact of fisheries on target species such as other large pelagic species (FG 8), common dentex (FG 24), and groupers (FG 27) while competitors or preys of those species may be positively impacted such as white seabream (FG 22), common two-banded seabream (FG 23) and brown meagre (Sciaena umbra) (FG 28) (Fig. 5). Results also highlighted that each fishery had a negative impact on itself due to self-competition for resources according to the MTI results and that competition between fisheries was complex (Fig 5). While small-scale fisheries had positive impacts on recreational fisheries in Cerbère-Banyuls MPA, recreational fisheries had slightly negative impacts on small-scale fisheries (Fig. 5a). On the contrary, small-scale fisheries had negative impacts on recreational fisheries in Cap de Creus and Medes Islands MPAs, while recreational fisheries had positive impacts on the small-scale fisheries (Fig. 5b andc).

Ecological impacts of protection over time

The best-fitted food wed temporal models were obtained when trophic interactions, fishing, and primary production anomaly were included in the model configuration for all MPAs (Step 8 in Table 3). However, for all MPAs the best model was not able to reproduce the trends of white seabream, common two-banded seabream, common dentex and groupers satisfactorily, which we selected as target groups of this study. Therefore, we moved through the fitting procedure analysis to find the model that was able to reproduce the trends of most of the groups (and specifically the target groups) and was highly significant.

We finally choose a model fit with 12 vulnerabilities and 3 spline points, 15 vulnerabilities and 3 spline points, and 11 vulnerabilities and 5 spline points for Cerbère-Banyuls, Cap de Creus and Medes Islands models respectively, as the best options (Step 8 in Table 3).

Observed biomass and catch time series were satisfactorily reproduced by model predictions for most of the target groups (Fig. 6 and SOM 6, 7 and 8) when using the chosen fitted model. The temporal models showed a non-significant biomass pattern of white seabream in Cerbère-Banyuls, firstly decreasing and later increasing, and in Cap de Creus, where firstly increased and later decreased, while in Medes Islands it presented a significant decreasing biomass trend (Fig 6). The models showed similar biomass patterns for common two-banded seabream, highlighting a non-significant trend in Cerbère-Banyuls and Cap de Creus and a significant decreasing trend in Medes Islands (Fig. 6). We observed a significant decreasing biomass trend of common dentex in Cerbère-Banyuls, while in Cap de Creus and Medes this group showed non-significant biomass trends (Fig. 6). The results highlighted a non-significant biomass pattern of groupers in Cerbère-Banyuls, while they significantly increased in Cap de Creus and Medes Islands (Fig 6). However, in Medes Islands the model did not capture well the overall declining trend of observations for groupers biomass.

Regarding the temporal changes of ecological indicators, the Kempton's Index showed a non-significant trend in Cerbère-Banyuls and Cap de Creus, while in Medes Islands, the Kempton's Index significantly increased during the simulated period (Fig. 7). The TL of the community presented a non-significant pattern in Cerbère-Banyuls and Cap de Creus, while in Medes Islands significantly decreased (Fig. 7). IUCN species biomass significantly declined in Cerbère-Banyuls, while in Cap de Creus and Medes Islands presented nonsignificant patterns (Fig. 7). Total catches highlighted non-significant trends in Cerbère-Banyuls and Cap de Creus (in the last one firstly increasing and then decreasing), while in Medes Islands, total catches significantly decreased (Fig. 7).

Discussion

Structure and functional traits of protected areas

Overall, ecological indicators showed similar patterns, with the highest values in Banyuls-Cerbère MPA, followed by Medes Islands and Cap de Creus MPAs. This may be related to differences in the ecological effectiveness of the three MPAs, which is partly explained by MPA design, management and implementation features (e.g. extent of area protected, fully protected area enforcement, time since protection, MPA design) [START_REF] Claudet | Marine reserves: size and age do matter[END_REF][START_REF] Guidetti | Italian marine reserve effectiveness: does enforcement matter?[END_REF][START_REF] Edgar | Global conservation outcomes depend on marine protected areas with five key features[END_REF][START_REF] Giakoumi | Ecological effects of full and partial protection in the crowded Mediterranean Sea: a regional meta-analysis[END_REF][START_REF] Franco | Linking home ranges to protected area size: The case study of the Mediterranean Sea[END_REF]. The lack of enforcement is one of the most relevant issues concerning MPAs in the Mediterranean context [START_REF] Fenberg | The science of European marine reserves: Status, efficacy, and future needs[END_REF]. Within this context, while Cerbère-Banyuls and Medes MPAs are considered to have a high level of enforcement [START_REF] Sala | The structure of Mediterranean rocky reef ecosystems across environmental and human gradients, and conservation implications[END_REF][START_REF] Giakoumi | Ecological effects of full and partial protection in the crowded Mediterranean Sea: a regional meta-analysis[END_REF]), Cap de Creus has been considered a MPA only on paper due to a lack of sufficient enforcement (Lloret, Zaragoza et al. 2008;Lloret, Zaragoza et al. 2008). In fact, total catches were higher in Cap the Creus than in the other two MPA; with Cerbère-Banyuls presenting the lowest total catches. Total catches are expected to be higher in well designed and enforced MPAs [START_REF] Halpern | Spillover from marine reserves and the replenishment of fished stocks[END_REF][START_REF] Lorenzo | Spillover from marine protected areas to adjacent fisheries has an ecological and a fishery component[END_REF], so our results could be explained by the allowed and/or real fishing effort and the small effects of these MPAs in PPAs and adjacent areas. Horta e Costa, Claudet et al. (2016) presented a novel classification system for MPAs by scoring allowed uses in each management unit based on their impacts on biodiversity. In this scale, Cerbère-Banyuls obtained a rate of 4.7, being a highly protected area and Medes 6.4 being less well protected (Horta e Costa, Claudet et al. 2016). Although Cap de Creus was not included in this study, we could expect a higher MPA index because of its smaller FPA and allowed uses. These features (smallest FPA and weak enforcement) and the fact that Cap de Creus is the newest (Table 1) MPA in the study, could suggest that the year of establishment, enforcement and the size of protected areas have a strong effect on MPA effects, in line with other studies [START_REF] Claudet | Marine reserves: size and age do matter[END_REF][START_REF] Guidetti | Italian marine reserve effectiveness: does enforcement matter?[END_REF].

Our results also indicated that both the small-scale and recreational fisheries can have a notable impact affecting organisms from lower to higher trophic levels. According to our study, small-scale fisheries had the largest negative ecological impacts. In fact, this fishery tends to develop its main activity inside the MPA or in surrounding areas [START_REF] Goñi | Spillover from six western Mediterranean marine protected areas: evidence from artisanal fisheries[END_REF][START_REF] Stelzenmüller | Spatial assessment of fishing effort around European marine reserves: implications for successful fisheries management[END_REF]. Recreational fisheries seemed to have larger impacts in the Spanish MPAs (Cap de Creus and Medes Islands MPA), while they had a negligible impact in Cerbère-Banyuls MPA, although the impact of recreational fisheries could be similar than the small-scale fisheries, as highlighted in coastal waters of Cap de Creus by Lloret, Zaragoza et al. (2008). This could be related to cultural reasons, differences in enforcement of each MPA and/or the reliability of catch data between study sites. In addition, our results highlighted that despite the recreational fisheries seemed to have an overall lower negative impact, it had strong negative impacts on vulnerable species such as common dentex and groupers (this group was mainly represented by Epinephelus marginatus) in Cap de Creus and Medes Islands MPAs, in line of the observed impact of small-scale and recreational fisheries on vulnerable species [START_REF] Lloret | Recreational and small-scale fisheries may pose a threat to vulnerable species in coastal and offshore waters of the western Mediterranean[END_REF], which is an important fact to consider when establishing management plans for these species. Our results also highlighted the competition for resources between both (smallscale and recreational) fisheries. Competition between both fisheries has become an important issue in coastal areas (especially in surrounding areas of MPAs), as they target similar species and fishing grounds [START_REF] Chuenpagdee | Transdisciplinarity for small-scale fisheries governance[END_REF][START_REF] Lloret | Small-scale coastal fisheries in European Seas are not what they were: Ecological, social and economic changes[END_REF]. In small and touristic areas such as the ones in the present study, a strong competition between small-scale and recreational fishers could be expected due to the current economic, societal and environmental challenges of the small scale fisheries and the increasing importance of recreational fisheries as a leisure activity in these areas [START_REF] Gómez | The Small-Scale Fisheries Guidelines as a Tool for Marine Stewardship: The Case of Cap de Creus Marine Protected Area, Spain. The Small-Scale Fisheries Guidelines[END_REF].

Ecological impacts of protection over time

The temporal dynamic model of Cerbère-Banyuls predicted an increasing biomass pattern of white seabream and common two-banded seabream while common dentex and groupers largely decreased. The decrease of common dentex and groupers could be related to the impact of fishing, although fishing effort did not increase during the study period and initial fishing mortalities were not high. The increase of both seabream species could be explained by the reduction of their predators and/or by some other processes (e.g., increase in resources). In fact, common dentex and groupers are the most important predators of these species in terms of total consumption.

In contrast, the temporal dynamic model of Cap de Creus predicted a decreasing biomass pattern of white seabream, common two-banded seabream and common dentex while groupers increased. These results suggest a recovery of groupers population in Cap de Creus [START_REF] Hereu Fina | Seguiment del medi marí al Parc Natural del Cap de Creus i al Parc Natural del Montgrí, les Illes Medes i el Baix Ter[END_REF], despite the low level of enforcement, that could cause biomass reductions of their preys such as both seabream species. Similarly, the temporal dynamic model of Medes Islands highlighted decreasing trends for white seabream, common two-banded seabream and common dentex, while groupers slightly increased. These results evidence that under protection the food-web effects can play an important role and there are winners and losers as a result of the new ecological state, where results of protection can include cascading effects of predators on preys species [START_REF] Edgar | Global conservation outcomes depend on marine protected areas with five key features[END_REF][START_REF] Cheng | Can marine reserves restore lost ecosystem functioning? A global synthesis[END_REF] as well as an increase in competitive interactions [START_REF] Micheli | Trajectories and correlates of community change in no-take marine reserves[END_REF].

Ecological indicators play an important role within the EBM framework as they are able to describe the condition of the ecosystem, its components and its functioning, to evaluate the impact of human activities in marine ecosystems and to inform management decisions [START_REF] Shin | Using indicators for evaluating, comparing, and communicating the ecological status of exploited marine ecosystems. 1. The IndiSeas project[END_REF][START_REF] Tam | Towards ecosystem-based management: identifying operational food-web indicators for marine ecosystems[END_REF]). In our study, ecological indicators showed overall contrasting results, as some indicators indicated recovery (e.g. Kempton's Q in the three MPAs) while others showed degradation patterns (e.g., TL of the community in the three MPAs). These contrasting results could respond to the limitation in historical data to represent ecosystem dynamics in the study area and could indicate limited recovery of species and ecosystems in the three MPAs.

The most notable effects of MPAs are an increase in abundance, average size and biomass inside protected areas [START_REF] Lester | Biological effects within no-take marine reserves: a global synthesis[END_REF][START_REF] Giakoumi | Ecological effects of full and partial protection in the crowded Mediterranean Sea: a regional meta-analysis[END_REF]. Due to the increased density inside the MPA, adults and juveniles may move to adjacent areas (spillover effect) and therefore, higher fisheries yields in adjacent areas are expected [START_REF] Goñi | Spillover from six western Mediterranean marine protected areas: evidence from artisanal fisheries[END_REF][START_REF] Lorenzo | Spillover from marine protected areas to adjacent fisheries has an ecological and a fishery component[END_REF]. Larger abundance and sizes could entail higher reproductive potential of species, and eggs and larvae from recovered populations could then be exported to external unprotected locations, including adjacent ones [START_REF] Gell | Benefits beyond boundaries: the fishery effects of marine reserves[END_REF][START_REF] Harrison | Larval export from marine reserves and the recruitment benefit for fish and fisheries[END_REF]). In addition, MPAs could restore ecosystem functioning [START_REF] Cheng | Can marine reserves restore lost ecosystem functioning? A global synthesis[END_REF]. Although some of our temporal

Data gaps, limitations and uncertainties

Overall the input data used was of acceptable quality compared to the distribution of pedigree values in other existing models [START_REF] Morissette | Complexity, cost and quality of ecosystem models and their impact on resilience: a comparative analysis, with emphasis on marine mammals and the Gulf of St[END_REF][START_REF] Lassalle | A toolbox to evaluate data reliability for whole-ecosystem models: Application on the Bay of Biscay continental shelf food-web model[END_REF], although they are among the lowest values in the Mediterranean Sea [START_REF] Corrales | Ecosystem structure and fishing impacts in the northwestern Mediterranean Sea using a food web model within a comparative approach[END_REF]. This could be due to the challenges of modelling small and local coastal areas, where specific data is a strong requirement. In fact, limitations of available data are a common concern for most Mediterranean MPAs [START_REF] Prato | Assessing interacting impacts of artisanal and recreational fisheries in a small Marine Protected Area (Portofino, NW Mediterranean Sea)[END_REF]Vilas, Coll et al. Submitted) despite there are monitoring programs inside these MPAs.

In general, there is a lack of biomass estimates for many functional groups, especially regarding benthic invertebrates. One of the main hurdles to evaluate the impact of fishing activities is to obtain realistic estimates of total catch (official and Illegal, Unregulated and Unreported) [START_REF] Pauly | Reported'versus' likely'fisheries catches of four Mediterranean countries[END_REF][START_REF] Pauly | Catch reconstructions reveal that global marine fisheries catches are higher than reported and declining[END_REF]. In our study this challenge is even higher due to the importance of recreational fishers and illegal fishing in coastal areas, specifically in MPA and surrounding areas (e.g., Lloret, Zaragoza et al. 2008;[START_REF] Ben Lamine | Comparing commercial, recreational and illegal coastal fishery catches and their economic values: A survey from the southern Mediterranean Sea[END_REF].

The capacity of ecosystem models to replicate trends increases with data availability and Despite these limitations, the models presented in this study were constructed using the best available data, including several sources of information including unpublished and ad hoc field data, expert knowledge and published data, and following the best practices in ecosystem modelling development [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF]. The present models represent the first attempt to develop temporal MPA models that are nearly located and confirm the capability of EwE to evaluate protection effects within an ecological perspective in small MPAs even in cases of data limitations.

Concluding remarks

Human activities have been concentrated in coastal ecosystems, resulting in a major modification of these marine areas [START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF]. Since coastal areas account for a large amount of the global ecosystem services [START_REF] Costanza | The value of the world's ecosystem services and natural capital[END_REF],

there is a need for the recovery of coastal marine ecosystems in order to ensure their ecological role. In addition, in recent decades, marine ecosystems have been increasingly

impacted by other stressors, directly or indirectly induced by multiple anthropogenic activities [START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF]. For example, sea warming and the spread of nonindigenous species represent a major challenge for conservation in Mediterranean ecosystems, especially in the Eastern Mediterranean Sea [START_REF] Lejeusne | Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea[END_REF][START_REF] Katsanevakis | Invading the Mediterranean Sea: biodiversity patterns shaped by human activities[END_REF]. Organisms and ecosystems that have been already

impacted by fishing are more vulnerable to these additional stressors (Occhipinti-Ambrogi and Savini 2003; [START_REF] Poloczanska | Responses of marine organisms to climate change across oceans[END_REF]. MPAs, when properly managed, have demonstrated to be an effective tool to protect target species and habitats to local stressors such as fishing [START_REF] Sala | No-take marine reserves are the most effective protected areas in the ocean[END_REF], while their role in promoting resilience to regional and global stressors such biological invasions and climate change is debated [START_REF] Giakoumi | Assessing the effects of marine protected areas on biological invasions: a global review[END_REF][START_REF] Roberts | Marine reserves can mitigate and promote adaptation to climate change[END_REF][START_REF] Giakoumi | Exploring the relationships between marine protected areas and invasive fish in the world's most invaded sea[END_REF]. In highly impacted and crowded areas like the Northwestern Mediterranean Sea, the establishment of well-designed networks of well-enforced MPA is necessary to achieve "clean, healthy and productive" oceans (Good Environmental Status) according to the Marine Strategy Framework Directive [START_REF] Fenberg | The science of European marine reserves: Status, efficacy, and future needs[END_REF]. Therefore, there is a need to better understand how current management options can contribute to the recovery and conservation of marine ecosystems and resources.

Within the EBM, ecosystem models and ecological forecasts have become an essential analytical and decision-making tool despite their limitations due to high uncertainties and complex ecosystem characteristics [START_REF] Link | Dealing with uncertainty in ecosystem models: the paradox of use for living marine resource management[END_REF][START_REF] Collie | Ecosystem models for fisheries management: finding the sweet spot[END_REF][START_REF] Maris | Prediction in ecology: promises, obstacles and clarifications[END_REF]). They have the potential to provide insights of possible future impacts on marine ecosystems and can offer guidance to decision-makers by evaluating the trade-off between different management units and identify those measures that have the potential to meet conservation objectives [START_REF] Fulton | A multi-model approach to engaging stakeholder and modellers in complex environmental problems[END_REF][START_REF] Acosta | Using scenarios and models to inform decision making in policy design and implementation. Methodological assessment of scenarios and models of biodiversity and ecosystem services[END_REF]. However, in order to better assess MPA effectiveness and design, there is a need to take into account economic and social attributes of the municipalities of the MPAs [START_REF] Bennett | Why local people do not support conservation: Community perceptions of marine protected area livelihood impacts, governance and management in Thailand[END_REF][START_REF] Clarke | An investigation into the socio-cultural dimension of determining MPA effectiveness[END_REF]. Therefore, the integration of the social and economic components into ecosystem models needs to be promoted. Also future research should try to assess how the benefits potentially delivered by MPAs in terms of species recovery and increase of fisheries catches actually affect human well-being of local communities inhabiting within and around MPAs.

This study is the baseline to develop future scenarios of alternative management options in order to maximize the impacts of MPAs to their surrounding areas and fisheries sustainability by alternative MPAs configurations. In addition, due to the vicinity of the three MPAs along a latitudinal gradient, this study is a part of a nested modelling approach with different geographic scales with the aim to assess the current effects of the actual MPA network and to perform simulations of alternative MPA network configuration (Coll, Steenbeek et al. 2019).

FPA= Fully protected Area, PPA= partially protected area, UPA= unprotected area flanking the MPA.

Table 2. Overall fitting procedure applied to the three MPA models following the methodology suggested by [START_REF] Mackinson | Which forcing factors fit? Using ecosystem models to investigate the relative influence of fishing and changes in primary productivity on the dynamics of marine ecosystems[END_REF] and [START_REF] Heymans | Best practice in Ecopath with Ecosim food-web models for ecosystem-based management[END_REF].

Table 3. Results of the fitting procedure of the three MPA of Cerbère-Banyuls, Cap de

Creus and Medes Islands. The table shows the statistically "best" model for each step. Vs = number of vulnerabilities estimated, PPsp = number of primary production spline points, k = number of parameters (Vs + PPsp), %IF = improved fit compared to the baseline AICc.

The "best" models are highlighted in bold.

Figure legends

  results showed potential recoveries (some target species, increases of predators and declines of prey, and some ecological indicators), our study evidences that they are still far from what we would expect from the temporal protection effects of MPAs, highlighting an overall modest historical positive effect of protection on these MPAs. These results are in line with Hereu Fina, Aspillaga Cuevas et al. (2017), which have pointed out illegal fishing as one of the main reasons for a lack of recovery and have called for further enhance the enforcement in Mediterranean MPAs.

Figure. 1 .

 1 Figure. 1. The study area encompassing the three multi-zone MPAs in the Northwestern Mediterranean Sea with the fully protected areas (FPAs), partially protected areas (PPAs) and unprotected neighboring areas (UPAs).

Figure. 2 .

 2 Figure. 2. Flow diagram of Cerbère-Banyuls (a), Cap de Creus (b) and Medes Island (c) MPA models. The numbers identify the functional groups of the model (listed in SOM 2).The size of each circle is proportional to the biomass of the functional group. The thickness of the connecting lines is proportional to the magnitude of their trophic flows.

Figure. 3 .

 3 Figure. 3. Ecological indicators estimated for the three multi-zone MPA models of Cerbère-Banyuls, Cap de Creus and Medes Islands. Boxplot shows the distribution of values for an ecological indicator derived from the Monte Carlo routine while the dot represents the value of the indicator in the baseline Ecopath balanced model.

Figure. 4 .

 4 Figure. 4.Functional groups plotted against Keystone Index[START_REF] Valls | Keystone species: toward an operational concept for marine biodiversity conservation[END_REF] and trophic level for Cerbère-Banyuls (a), Cap de Creus (b) and Medes (c) multi-zone MPA models.The numbers identify the functional groups of the model (listed in SOM 2). The size of the circles is proportional to the biomass of the functional group.

Figure. 5 .

 5 Figure. 5. Mixed Trophic Impact (MTI) analysis of the three MPA applied to the fisheries in a) Cerbère-Banyuls, b) Cap de Creus, and c) Medes Islands multi-zone MPA models.

Figure. 6 .

 6 Figure. 6. Predicted (solid lines) versus observed (dots) biomass (t•km -2 ) values for target species of Cerbère-Banyuls, Cap de Creus and Medes Islands multi-zone MPAs models.Blue shadows represent the 5% and 95% percentiles obtained using the Monte Carlo routine. Rho and p-values result from Spearman's rank correlation test.

Figure. 7 .

 7 Figure. 7. Temporal trends of ecological indicators of Cerbère-Banyuls, Cap de Creus and Medes Islands multi-zone MPAs models. Blue shadows represent the 5% and 95% percentiles obtained using the Monte Carlo routine. Rho and p-values result from Spearman's rank correlation test.

  quality[START_REF] Giron-Nava | Quantitative argument for long-term ecological monitoring[END_REF]. Within this context, in Cerbère-Banyuls a maximum of four data points of biomass per functional group were available for the fitting procedure. The Medes Islands MPA model had the longest available time series, and thus presented a better and more informative fit, while the Cap de Creus MPA model showed an intermediate situation. This limited the capability to calibrate and validate the models and to track ecosystem dynamic status. Therefore, our results highlight the need to further monitor MPAs within the Mediterranean Sea and should inform future scientific research objectives in the area.

Acknowledgements

This work was funded by EU Research Project SAFENET ("Sustainable Fisheries in EU Mediterranean waters through Network of MPAs). Call for proposal MARE/2014/41, Grant Agreement n. 721708.

Author Contributions

XC, DV, CP, JS and MC designed the study and analyses. CP, JS, JC, JL, AC, ADF, TF, AL, GP, RS, PS, PG and MC contributed with material and analysis tools. XC, DV and MC developed the models. XC, DV, CP, JS, JC, JL, AC, ADF; AL, GP, PG and MC interpreted the results. XC wrote the manuscript with assistance from DV, JC and MC and input from CP, JS, JL, AC, ADF, TF, AL, GP, RS, PS and PG. Table 1. Surface area (km²) covered by management units (MU) and year of creation of each Marine Protected Area (MPA) of Cerbère-Banyuls, Cap de Creus and Medes Islands.

Tables legends


Stéphane Jaffard

Stéphane Seuret

Herwig Wendt
email: herwig.wendt@irit.fr

Patrice Abry

Roberto Leonarduzzi

A multivariate multifractal analysis of lacunary wavelet series

Multivariate multifractal analysis proposes to estimate the multivariate multifractal spectrum of several signals as a way to reveal the correlations between their singularity sets. Recently discovered failures of the "natural" multifractal formalism challenge its foundations, and raise questions on the information it provides. We investigate these questions by supplying some general results on the multivariate multifractal analysis of processes satisfying the large intersection property, and illustrate it on random lacunary wavelet series.

I. INTRODUCTION

Context: Multifractal analysis. The purpose of multifractal analysis is to analyze functions or measures (random or deterministic) by estimating the size of the sets of points where they have a given pointwise regularity. Mathematically, this information is encapsulated through the multifractal spectrum D(H), which yields the Hausdorff dimensions of the sets E(H) where a pointwise regularity exponent takes the value H [START_REF] Parisi | Fully developed turbulence and intermittency[END_REF]. In practice, for a signal or an image, a direct estimation of D(H) is unfeasible; the purpose of multifractal formalisms is to estimate it through the computation of a Legendre spectrum which, for univariate data, yields an upper bound of D(H) [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF], and coincides with it for important classes of models (cf. e.g., [START_REF] Abry | Irregularities and scaling in signal and image processing: Multifractal analysis[END_REF] and ref. therein). Thus, the derivation of practically computable Legendre spectra as close as possible to D(H) is a central issue. Related works. Recent applications involve simultaneous recordings of large collections of signals or images, thus calling for a multivariate multifractal analysis. When several functions X 1 (x), • • • , X M (x) are involved, the natural mathematical quantity considered is the multivariate multifractal spectrum D(H 1 , • • • H M ), which yields the Hausdorff dimensions of the sets of points where the regularity exponent of each signal X i takes a given value H i . A natural multivariate multifractal formalism was proposed in [START_REF] Meneveau | Joint multifractal measures -theory and applications to turbulence[END_REF], initially in the context of turbulence, and has been used in a wide range of applications. However, its validity was recently questioned, and a dramatic turning point was the discovery that, in general, this formalism does not yield an upper bound for the multivariate spectrum [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF]; counterexamples are supplied by classical models, such as independent Lévy processes, or couples of anti-correlated binomial cascades. This unexpected failure raises the problems of understanding under which conditions the formalism is valid (i.e. yields an upper bound) and which information it supplies when such conditions are not fulfilled (or cannot be verified), in particular what it reveals about the correlations between the locations of the singularities of the components X i . Goals, contributions and outline. We investigate these questions following the guideline supplied by a simple stochastic model: Lacunary wavelet series (LWS). Their sample paths were shown to be multifractal [START_REF] Jaffard | On lacunary wavelet series[END_REF], and we study couples of such processes, either independent or with tuned correlations. In particular, we show that, in this case too, the natural multivariate formalism fails. This confirms a major source of failure when the singularity sets of the processes satisfy the large intersection property; we back this intuition by general mathematical results on multivariate processes which share this property. In Section II, we recall the basic notions in multivariate multifractal analysis. In Section III we derive mathematical results on the bivariate multifractal analysis of processes satisfying the large intersection property. In Section IV we study, theoretically and numerically, how the bivariate spectrum of LWS changes when their correlation is tuned.

II. MULTIVARIATE MULTIFRACTAL ANALYSIS

For simplicity, in the following, we only consider the bivariate case (M = 2) where functions are defined on R. Higher dimensional extensions easily follow.

A. Multifractal spectrum

Multifractal analysis quantifies the fluctuations of the regularity of a signal X(t). Pointwise regularity of functions is usually measured using the Hölder exponent, h(t), defined as follows: X ∈ C α (t), α ≥ 0, if there exist a polynomial P t with deg(P t ) < α and a constant C > 0 such that

|X(t + a) -P t (t + a)| ≤ C|a| α , |a| → 0, (1)
h(t) is the largest α such that (1) is satisfied,

h(t) sup{α : X ∈ C α (t)} ≥ 0, (2)
see e.g. [START_REF] Abry | New exponents for pointwise singularity classification[END_REF] for alternative notions of pointwise regularity. Let h(t) (h 1 (t), h 2 (t)) denote the Hölder exponents of the components of the bivariate signal X = (X 1 , X 2 ). Let H = (H 1 , H 2 ); the bivariate multifractal spectrum D(H 1 , H 2 ) of X is defined as (see e. g. [START_REF] Meneveau | Joint multifractal measures -theory and applications to turbulence[END_REF], [START_REF] Jaffard | Multivariate multifractal analysis[END_REF])

D(H) dim H (E H ), where E H = {t : h(t) = H} (3)
and dim H denotes the Hausdorff dimension (by convention dim H (∅) = -∞). It provides a global, geometric description of the pointwise regularity properties of the components X i ; here and in the following the subscript i refers to the i-th component (i = 1, 2). Its precise shape and width quantify key information regarding correlations between the exponents h i (t) [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF], [START_REF] Wendt | Assessing cross-dependencies using bivariate multifractal analysis[END_REF], [START_REF] Leonarduzzi | Multifractal characterization for bivariate data[END_REF].

B. Wavelet leader bivariate multifractal formalism

The estimation of the multifractal spectrum requires numerically robust formulas: multifractal formalisms. A univariate formalism constructed from the multiscale statistics of wavelet leaders [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF], [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF], has recently been extended to the multivariate setting, see e.g. [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF], [START_REF] Jaffard | Multivariate multifractal analysis[END_REF]- [START_REF] Leonarduzzi | Multifractal characterization for bivariate data[END_REF], [START_REF] Slimane | Baire typical results for mixed hölder spectra on product of continuous besov or oscillation spaces[END_REF], and is recalled now.

Let N > 0, and let ψ denote a "mother wavelet", i.e. an oscillating well-localized function in

C N , such that {ψ j,k (t) = 2 -j/2 ψ(2 -j t -k)} (j,k)∈Z 2 forms an orthonormal basis of L 2 (R). The wavelet coefficients of a univariate function f are d f (j, k) = 2 -j/2 ψ j,k |f . The wavelet leaders of f are f (j, k) sup λ ⊂3λ j,k |d f (λ )|, where λ j,k = [k2 j (k + 1)2 j )
denotes the dyadic interval of size 2 j and 3λ j,k stands for the union of λ j,k with its two closest neighbors [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF], [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. In the bivariate case, the structure functions S(q 1 , q 2 , j) and scaling function ζ(q 1 , q 2 ) of X = (X 1 , X 2 ) are defined by

S(q 1 , q 2 , j) 2 j k L X1 (j, k) q1 L X2 (j, k) q2 ∼ c q 2 jζ(q1,q2) (4)
(for j → -∞), where the order q i of empirical moments is taken in some interval including zero. In contradistinction with the univariate case, the bivariate Legendre spectrum L(H) = inf q (1 + q, H -ζ(q)), q = (q 1 , q 2 ), [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF] does not always yield an upper bound for D(H), see [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF] for a discussion of conditions under which it is the case.

The case where X 1 and X 2 are independent was studied in [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF]: Under a stationarity hypothesis, Legendre spectra satisfy

L(H 1 , H 2 ) = L 1 (H 1 ) + L 2 (H 2 ) -1. ( 6
)
This formula is backed by the remark that, applied to multifractal spectra, the corresponding formula

D(H 1 , H 2 ) = D 1 (H 1 ) + D 2 (H 2 ) -1 (7)
follows the generic rule that, for smooth sets, codimensions of intersections of sets add up [START_REF] Mattila | Geometry of Sets and Measures in Euclidean Spaces: Fractals and Rectifiability[END_REF] (note that, if (7) yields a negative dimension it is understood that the corresponding set is empty). However, a different formula applies to models for which the lower singularity sets E ≤ (H) = {x ∈ R : h(x) ≤ H} are sets with large intersections, i.e. any couple among them satisfies dim H (A ∩ B) = min(dim H A, dim H B), see [START_REF] Falconer | Sets with large intersection properties[END_REF], [START_REF] Durand | Sets with large intersection and ubiquity[END_REF]. Though the E H are not sets with large intersection, this leads to the alternative conjectured formula for the bivariate spectrum of independent processes of this type:

D(H 1 , H 2 ) = min(D 1 (H 1 ), D 2 (H 2 )). (8)

III. LARGE INTERSECTION PROCESSES

For several standard models such as Lévy processes, lacunary and random wavelet series, the sets E ≤ (H) are sets of large intersection, see [START_REF] Aubry | Random wavelet series[END_REF]- [START_REF] Durand | Singularity sets of Lévy processes[END_REF] and ref. therein. The reason is that these sets are obtained as limsup of intervals centered at equidistributed points; as a consequence the following property holds: For s = dim H (E ≤ (H)), they belong to the G s class of Falconer, which is composed of sets A which are G δ sets (i.e. countable intersections of open sets) and such that, for any countable sequence of similarity transformations f i , dim H (∩f i (A)) ≥ s. An important consequence is that there exists a s-dimensional Hausdorff measure of A which is positive, see [START_REF] Falconer | Sets with large intersection properties[END_REF], [START_REF] Durand | Sets with large intersection and ubiquity[END_REF]. We now investigate which implications on bivariate spectra hold in this setting.

We consider a couple of functions (X 1 , X 2 ), deterministic or random, without assumptions on their correlations in the random case. We assume that there exists an interval

[H -, H + ] such that: ∀H ∈ [H -, H + ], the set E ≤ (H) belongs to G s with s = D ≤ (H) (in
≤ (H 1 , H 2 ) dim H ({t : h 1 (t) ≤ H 1 and h 2 (t) ≤ H 2 }) satisfies D ≤ (H 1 , H 2 ) = dim H (E ≤ 1 (H 1 ) ∩ E ≤ 2 (H 2 )) = min D ≤ 1 (H 1 ), D ≤ 2 (H 2 ) .
In particular, it can not yield information on correlations between the two processes. On the opposite, determining the bivariate spectrum D(H 1 , H 2 ) requires information on these correlations. However, we will see that, under an additional assumption, a part of this spectrum does not depend on these correlations. We now assume that, for each process X i (i = 1, 2) the function H → D ≤ i (H) is continuous and strictly increasing on [H - i , H + i ] (these assumptions are verified by the examples we mentioned). The corresponding process will be called a large intersection process. It follows that, for i = 1, 2,

D i (H) = D ≤ i (H).
We define the ridge of the bivariate spectrum as the set:

{(H 1 , H 2 ) : D ≤ 1 (H 1 ) = D ≤ 2 (H 2 )}.
Note that, if the two processes share the same law, then the ridge simply is a segment of the diagonal

H 1 = H 2 .
Theorem 1: Under the previous assumptions, the bivariate spectrum of (X 1 , X 2 ) satisfies that, for any (H 1 , H 2 ) which belongs to the ridge,

D(H 1 , H 2 ) = D 1 (H 1 )(= D 2 (H 2 )). Proof: The upper bound is straightforward since E(H 1 , H 2 ) is included in E 1 (H 1 ). Conversely, the G s class assump- tion implies that E ≤ (H 1 , H 2 ) has a positive D ≤ (H 1 , H 2 )- dimensional Hausdorff measure; E(H 1 , H 2 ) contains E ≤ (H 1 , H 2 ) - H<H1 E ≤ (H, H 2 ) - H<H2 E ≤ (H 1 , H).
By the increasing spectrum assumption, for H < H 1 (resp. H < H 2 ), the sets E ≤ (H, H 2 ) and E ≤ (H 1 , H) have a vanishing D ≤ (H 1 , H 2 )-dimensional Hausdorff measure; since the unions are countable, E(H 1 , H 2 ) has a positive D ≤ (H 1 , H 2 )dimensional Hausdorff measure; hence the lower bound holds.

Remark: The support of the bivariate spectrum can boil down to the ridge, for instance when the functions X 1 and X 2 coincide. Thus the determination of the bivariate spectrum outside of the ridge necessarily requires information on the correlation between the two processes.

IV. LACUNARY WAVELET SERIES

We start by recalling the model studied in [START_REF] Jaffard | On lacunary wavelet series[END_REF]. For the sake of simplicity, the random process X is constructed on [0, 1]. Let η ∈ (0, 1) and α > 0. At scale j ≤ 0, [2 -ηj ] dyadic intervals (λ ⊂ [0, 1]) are drawn at random, and for these intervals, c 1 λ = 2 αj ; the other wavelet coefficients are set to 0. Additionally, the coefficients are independent at each scale. The parameters of the series are its lacunarity η and its uniform regularity H min = α. The multifractal spectrum of X is

D(H) = ηH/α if H ∈ [α, α/η] = -∞ else.
We now consider two LWS X 1 , X 2 of parameters (η 1 , α 1 ) and (η 2 , α 2 ). The results of Section III apply: Without any correlation assumption,

∀(H 1 , H 2 ) ∈ [α 1 , α 1 /η 1 ] × [α 2 , α 2 /η 2 ], D ≤ (H 1 , H 2 ) = min η 1 H 1 α 1 , η 2 H 2 α 2 . ( 9
)
The ridge is included inside the line η1H1 α1 = η2H2 α2 and, on the ridge, D(H 1 , H 2 ) is given by [START_REF] Wendt | Assessing cross-dependencies using bivariate multifractal analysis[END_REF].

Additionally, since the sets

E 1 (H 1 ) and E 2 (H 2 ) are empty if H 1 < α 1 or if H 2 < α 2 , it follows that E ≤ (α 1 , α 2 ) = E(α 1 , α 2 ), so that D(α 1 , α 2 ) = D ≤ (α 1 , α 2 ) (note that (α 1 , α 2 ) is not in the ridge except if η 1 = η 2 ).
1) Independent LWS: If X 1 and X 2 are independent, then the assumptions under which the derivation of ( 6) is performed in [START_REF] Jaffard | Multifractal formalisms for multivariate analysis[END_REF] are verified, so that, if

(H 1 , H 2 ) ∈ [α 1 , α 1 /η 1 ] × [α 2 , α 2 /η 2 ], then L(H 1 , H 2 ) = η 1 H 1 α 1 + η 2 H 2 α 2 -1 (10)
(and L(H 1 , H 2 ) = -∞ else). Therefore L(H 1 , H 2 ) does not supply an upper bound for D(H 1 , H 2 ), see Fig. 1 (c).

2) Correlated lacunary wavelet series: We consider a first LWS X 1 of parameters (α, η), and construct a second random process X 2 , correlated with X 1 : We pick ω > η and β > 0. For each interval λ picked in the construction of X 1 , we consider the unique dyadic interval λ of scale j = [ωj] ( where j is the scale of λ) which contains λ; then we set c 2 λ = 2 β j , and all other wavelet coefficients of X 2 are set to 0. Considered independently of X 1 , X 2 is another lacunary wavelet series of lacunarity parameter η = η/ω and uniform regularity β (note that, strictly speaking, this is not true because the number of nonvanishing wavelet coefficients at a given scale is random and may slightly differ from [2 ηj ] but, clearly, this does not alter the conclusions concerning the spectra).

Remark: The case ω = 1 and α = β corresponds to the case where X 1 = X 2 . The case ω = 1 and α = β corresponds to the case where X 1 is analyzed jointly with its fractional primitive of order β -α.

We now recall how the pointwise regularity of X 1 is derived in [START_REF] Jaffard | On lacunary wavelet series[END_REF]. Let Λ 1 denote the random collection of dyadic intervals indexing nonvanishing wavelet coefficients of X 1 , and Λ 2 denote those of X 2 . Let θ ∈ [η, 1]; we consider "dilated" intervals constructed as follows: λ θ has the same center as λ ∈ Λ 1 , and width 2 θj , and we denote by F θ the limsup of the λ θ (i.e. the set of points that belong to an infinite number of the λ θ for λ ∈ Λ 1 ). The key result that allows to derive the multifractal spectrum is that, if

x ∈ F θ , then h 1 (x) ≤ α/θ and if x / ∈ F θ , then h 1 (x) ≥ α/θ. It follows that x ∈ G θ := θ <θ F θ - θ >θ F θ ⇐⇒ h 1 (x) = α/θ,
and the large intersection process assumption implies that, if

θ ∈ [η, 1], then dim(G θ ) = η/θ.
In order to derive the bivariate spectrum of X 1 and X 2 , we separate two cases:

• Let θ ∈ [ω, 1]. If x ∈ F θ ,
then we already know that h 1 (x) ≤ α/θ, and x belongs to an infinite number of the λ (or one of its two closest neighbors), so that h 2 (x) = β.

• Let θ ∈ [η, ω]. If x ∈ F θ ,
then as before, h 1 (x) ≤ α/θ, and x belongs to an infinite number of the λ θ , which can be seen as an interval of same center as λ and of width | λ| θ/ω so that h 2 (x) ≤ βω/θ. Lower bounds follow in the same way, and we obtain that the regularity at x of both series is determined by the G θ to which x belongs. For a given θ ∈ [ω, 1], the dimension of the corresponding set is η/θ. The following result follows, where the spectrum is parametrized as a function of H 1 .

Theorem 2: The bivariate multifractal spectrum of two correlated lacunary wavelet series of parameters (α, η, ω) is supported by the two following (not colinear) straight lines:

if H 1 ∈ α, α ω , H 2 = β and D(H 1 , H 2 ) = H 1 η α if H 1 ∈ α ω , α η , H 2 = H 1 βω α and D(H 1 , H 2 ) = H 1 η α .
V. NUMERICAL ILLUSTRATIONS Monte Carlo simulation. We illustrate the above theoretical results using 100 realizations of independent LWS of sample size 2 16 , for two different sets of parameters: (α 1 , η 1 ) = (0.6, 0.85) and (α 2 , η 2 ) = (0.6, 0.85) for X 1 and X 2 , respectively, corresponding to a multifractal spectrum with relatively narrow support, and (α 1 , η 1 ) = (0.5, 0.4) and (α 2 , η 2 ) = (0.4, 0.4) for X 1 and X 2 , respectively, which corresponds to a multifractal spectrum with larger support and marginal spectra that do not coincide. Sample paths and increments for the two parameter settings are plotted in Fig1. 1 (a). Estimation is performed by linear regression of log 2 S(q 1 , q 2 , j) for scales j from -10 to -6 and -7 to -4 for the two parameter sets, respectively, for a sufficiently wide range of positive and negative moments (q 1 , q 2 ) so as to cover the theoretical support of the multifractal spectra. Numerical process synthesis and estimation procedures have been implemented by ourselves and are available upon request. Structure functions. Log-log plots of average structure functions S(q 1 , q 2 , j) for several values of (q 1 , q 2 ) are plotted as a function of scale j in Fig. 1 (b). They show that the power laws postulated in Eq. ( 4) hold for a large range of scales, and in particular for the range over which linear regressions are performed for the estimation of ζ(q 1 , q 2 ) and L(H).

Theoretical spectra and Legendre spectra. Fig. 1 (c) plots averages of Legendre spectra estimates (5) (red surfaces and lines) together with theoretical spectra, once corresponding to the theoretical expression [START_REF] Leonarduzzi | Multifractal characterization for bivariate data[END_REF] for the Legendre spectrum (left columns) and once to that for the multifractal spectrum (9) (right columns), respectively. Clearly, the predictions for the Legendre spectrum and theoretical bivariate multifractal spectrum do not coincide, while the marginal spectra (solid lines) do. Moreover, the estimated Legendre spectra provide a better upper bound for expression [START_REF] Leonarduzzi | Multifractal characterization for bivariate data[END_REF] than for expression [START_REF] Wendt | Assessing cross-dependencies using bivariate multifractal analysis[END_REF], which further corroborates the above theoretical results. The fact that the estimated and theoretical Legendre spectra do not match perfectly is not surprising since it is well documented that numerical estimation for processes with lacunarity is notoriously difficult.

VI. CONCLUSION

The present study confirms the failure of the multifractal formalism for couples of processes which satisfy the large intersection property. In this case, the Legendre spectrum does not bound the multifractal spectrum, and fails to provide information describing the cross-correlations between components, which was its original goal. Indeed, it is shown that a part of the bivariate spectrum is independent of the correlations between the components. These theoretical results are confirmed by numerical results obtained on sample paths of LWS. This raises two important questions: How could a non-trivial upper bound of the multifractal spectrum be recovered in practice? and which information is supplied by the bivariate Legendre spectrum? The example of correlated LWS yields new insights on the key problem of understanding how the correlations of several processes are reflected in their multivariate spectrum.

ACKNOWLEDGMENT

Work supported by ANR-16-CE33-0020 MultiFracs, France.

(α 1 (α 1 Fig. 1 .

 111 Fig. 1. Increments and sample paths (a), and structure functions (b) for bivariate LWS for two different parameter settings, and estimated (red) vs. theoretical (blue) multifractal spectra (c); the theoretical spectra (in blue) plotted in (c.1) and (c.3) correspond to expression (10) for the Legendre spectrum, those in (c.2) and (c.4) correspond to the expression (9) for the multifractal spectrum.

  the examples we mentioned, the range [H -, H + ] coincides with [H min , H max ] of all Hölder exponents taken by the process). It follows that (in this range of exponents) the bivariate lower spectrum D