Keywords: ecological disturbance, coral reef, network analysis, trophic modelling, Ecopath

Ecosystem that appear resilient to a disturbance can hide recovery debts.

Introduction

Despite conservation efforts, ecosystems are increasingly being degraded worldwide [START_REF] Maxwell | Biodiversity: The ravages of guns, nets and bulldozers[END_REF]. Chronic human stressors play an important role in altering ecosystems [START_REF] Vanbergen | Threats to an ecosystem service: pressures on pollinators[END_REF][START_REF] Crowther | Mapping tree density at a global scale[END_REF][START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF], but pulse natural disturbance events can also lead to deep changes in composition of biological communities [START_REF] Dulvy | Extinction vulnerability in marine populations[END_REF][START_REF] Scheffer | Early-warning signals for critical transitions[END_REF]. Such changes can affect ecosystem functioning, thus altering the delivery of ecosystem services [START_REF] Worm | Impacts of biodiversity loss on ocean ecosystem services[END_REF][START_REF] Balvanera | Linking biodiversity and ecosystem services: current uncertainties and the necessary next steps[END_REF]. Managing for sustainability not only requires an understanding of how an impact may propagate through a whole ecosystem but also the identification of recovery pathways at the ecosystem level [START_REF] Lotze | Recovery of marine animal populations and ecosystems[END_REF].

Estimating the recovery time of a specific population or a whole ecosystem following a disturbance is not straightforward [START_REF] Lotze | Recovery of marine animal populations and ecosystems[END_REF]. Recovery is considered as "a return to the initial state" of a specific response variable [START_REF] Lotze | Recovery of marine animal populations and ecosystems[END_REF]. It can be measured through changes in species abundance [START_REF] Alter | Dna evidence for historic population size and past ecosystem impacts of gray whales[END_REF] or cover [START_REF] Osborne | Delayed coral recovery in a warming ocean[END_REF], spawning stock size [START_REF] Hutchings | Marine fish population collapses: consequences for recovery and extinction risk[END_REF], demographic metrics [START_REF] Gerber | The use of demographic sensitivity analysis in marine species conservation planning[END_REF] or more recently changes in trophic structure [START_REF] Hempson | Ecosystem regime shifts disrupt trophic structure[END_REF]).

However, time-lags might exist at ecosystem scale and populations do not necessarily respond in the same pace [START_REF] Tilman | Habitat destruction and the extinction debt[END_REF][START_REF] Dullinger | Europe's other debt crisis caused by the long legacy of future extinctions[END_REF]. Time-delayed response of species to disturbance, also called the "extinction debts" [START_REF] Tilman | Habitat destruction and the extinction debt[END_REF], might seriously underestimate the impact of the disturbance at the whole ecosystem scale while a system might exhibit a "recovery debt" [START_REF] Moreno-Mateos | Anthropogenic ecosystem disturbance and the recovery debt[END_REF].

The recovery debt was recently conceptualized as the resulting deficit in biodiversity and ecosystem functions during the process of recovery [START_REF] Moreno-Mateos | Anthropogenic ecosystem disturbance and the recovery debt[END_REF]. Even if signs of recovery exist in many impacted ecosystems [START_REF] Jones | Rapid recovery of damaged ecosystems[END_REF][START_REF] Graham | Predicting climate-driven regime shifts versus rebound potential in coral reefs[END_REF][START_REF] Cole | Recovery and resilience of tropical forests after disturbance[END_REF], ecosystems might not have fully recovered and will potentially not in the future [START_REF] Moreno-Mateos | Anthropogenic ecosystem disturbance and the recovery debt[END_REF][START_REF] Jones | Restoration and repair of earth's damaged ecosystems[END_REF]). As humans rely on fully functioning ecosystem, recovery debt is expected to affect ecosystem services supply [START_REF] Isbell | The biodiversity-dependent ecosystem service debt[END_REF]) and be of pivotal importance for sustainability of the social-ecological systems [START_REF] Lafuite | Delayed behavioural shifts undermine the sustainability of social-ecological systems[END_REF].

Evaluation of changes in ecosystem structure and functioning following a disturbance can be achieved through the use of a diversity of quantitative ecosystem indicators [START_REF] Heymans | Global patterns in ecological indicators of marine food webs: a modelling approach[END_REF]Steenbeek 2017). Indicators derived from the Ecosystem Network Analysis (ENA) [START_REF] Ulanowicz | An informational synthesis of ecosystem structure and function[END_REF], which use network theory information, have the advantage of quantifying direct and indirect interactions and identifying emergent food-web properties, without reducing the model to only few dominant processes [START_REF] Fath | Ecological network analysis: network construction[END_REF][START_REF] Kones | Are network indices robust indicators of food web functioning? a monte carlo approach[END_REF]). These indicators have been widely used to analysed food-web changes over time (e.g. [START_REF] Heymans | Changes in the northern benguela ecosystem over three decades: 1970s, 1980s, and 1990s[END_REF]2007[START_REF] Tomczak | Ecological network indicators of ecosystem status and change in the baltic sea[END_REF] to compare ecosystems functioning [START_REF] Heymans | Global patterns in ecological indicators of marine food webs: a modelling approach[END_REF] or to test the sensitivity of candidate indicators to a given pressure [START_REF] Bourdaud | New trophic indicators and target values for an ecosystem-based management of fisheries[END_REF][START_REF] Chaalali | From species distributions to ecosystem structure and function: A methodological perspective[END_REF]. Due to the complexity to consider uncertainty in complex food web models, most of older analyses using food web models relied on single indicator values. However, to identify significant changes over time and improve interpretation of model, derived indices uncertainty should be accounted for in model parametrization [START_REF] Fulton | Approaches to end-to-end ecosystem models[END_REF][START_REF] Guesnet | Incorporating food-web parameter uncertainty into ecopath-derived ecological network indicators[END_REF][START_REF] Corrales | Hindcasting the dynamics of an eastern mediterranean marine ecosystem under the impacts of multiple stressors[END_REF], Piroddi et al. 2017).

Coral reefs are a good case studies to trial ecosystem network approaches and candidate indicators aimed at identifying recovery pathways following disturbances. They are biologically and functionally rich ecosystems [START_REF] Bellwood | Limited functional redundancy in high diversity systems: resilience and ecosystem function on coral reefs[END_REF]) that provide vital and valuable services to coastal communities such as coastal protection, food production and recreational activities contributing to social and cultural identity [START_REF] Moberg | Ecological goods and services of coral reef ecosystems[END_REF]. Coral reefs are also among the most vulnerable ecosystems [START_REF] Halpern | Evaluating and ranking the vulnerability of global marine ecosystems to anthropogenic threats[END_REF][START_REF] Holbrook | Reef fishes in biodiversity hotspots are at greatest risk from loss of coral species[END_REF]. They evolve under several natural disturbance regimes, such as extreme climatic events, outbreak of corallivorous crown-of-thorns sea stars (COTS), or mass bleaching events [START_REF] Ban | Evidence for multiple stressor interactions and effects on coral reefs[END_REF][START_REF] Hughes | Global warming and recurrent mass bleaching of corals[END_REF]. In their evolutionary history, coral communities have demonstrated the capacity to recover from these cycles of disturbances [START_REF] Jackson | Pleistocene perspectives on coral reef community structure[END_REF][START_REF] Gilmour | Recovery of an isolated coral reef system following severe disturbance[END_REF][START_REF] Lamy | Three decades of recurrent declines and recoveries in corals belie ongoing change in fish assemblages[END_REF]. Today, global change and increasing demography are rapidly altering the dynamics of theses ecosystems (Bozec andMumby 2015, Williams et al. 2015): Coral reefs are increasingly shifting from a coral-dominated, biologically diverse, productive, states to macroalgae-dominated, poorer, less productive states [START_REF] Gardner | Long-term region-wide declines in caribbean corals[END_REF][START_REF] Bruno | Regional decline of coral cover in the indo-pacific: timing, extent, and subregional comparisons[END_REF][START_REF] Ainsworth | Coral-algal phase shifts alter fish communities and reduce fisheries production[END_REF][START_REF] Hughes | Global warming and recurrent mass bleaching of corals[END_REF]. Seen as major threat, overfishing is one of the most studied stressor on coral reef ecosystem components [START_REF] Arias-González | Trophic models for investigation of fishing effect on coral reef ecosystems[END_REF][START_REF] Weijerman | Comparison of coral reef ecosystems along a fishing pressure gradient[END_REF], Graham et al. 2017). However less is known about how natural disturbance regimes affect exploited reef structure and functioning [START_REF] Ainsworth | Coral-algal phase shifts alter fish communities and reduce fisheries production[END_REF][START_REF] Hempson | Coral reef mesopredators switch prey, shortening food chains, in response to habitat degradation[END_REF]. Furthermore, most coral reef studies have only addressed changes of small subsets of ecosystem interactions [START_REF] Mumby | Phase shifts and the stability of macroalgal communities on caribbean coral reefs[END_REF][START_REF] Adam | How will coral reef fish communities respond to climate-driven disturbances? insight from landscape-scale perturbations[END_REF][START_REF] Pratchett | Changes in biodiversity and functioning of reef fish assemblages following coral bleaching and coral loss[END_REF], neglecting wider off-stage changes on ecosystem structuring and functioning and the potential recovery debt.

Here, we investigate how a pulse disturbance affecting coral communities can cascade through a whole coral reef ecosystem and assess ecosystem transformations along its recovery pathways. We use as case study, the coral reef ecosystem of Moorea Island (French Polynesia) that undergone a serious decline in coral cover and shift in fish species composition following a crown-of-thorns sea stars outbreak and a cyclone [START_REF] Lamy | Understanding the spatio-temporal response of coral reef fish communities to natural disturbances: insights from beta-diversity decomposition[END_REF]2016)

Methods

Studied ecosystem and environmental disturbances

Moorea is a volcanic high island located in the central south Pacific, 20 km north-west of Tahiti. This island is enclosed by an offshore barrier reef, the fore reef, that defines a shallow lagoon. In the past several decades, the offshore coral communities were severely impacted by repeated disturbances, including cyclones, bleaching events and crown-of-thorns sea stars (COTS) outbreaks [START_REF] Adjeroud | Recurrent disturbances, recovery trajectories, and resilience of coral assemblages on a south central pacific reef[END_REF][START_REF] Lamy | Three decades of recurrent declines and recoveries in corals belie ongoing change in fish assemblages[END_REF].

Each event have resulted in a drastic decline of living coral cover, as well as a turnover in fish communities [START_REF] Adam | Herbivory, connectivity, and ecosystem resilience: response of a coral reef to a large-scale perturbation[END_REF][START_REF] Lamy | Understanding the spatio-temporal response of coral reef fish communities to natural disturbances: insights from beta-diversity decomposition[END_REF]. However, until now, Moorea's offshore reef has displayed a high recovery potential by returning to a coral dominated ecosystems within a decade [START_REF] Adjeroud | Recurrent disturbances, recovery trajectories, and resilience of coral assemblages on a south central pacific reef[END_REF][START_REF] Adam | Herbivory, connectivity, and ecosystem resilience: response of a coral reef to a large-scale perturbation[END_REF][START_REF] Lamy | Three decades of recurrent declines and recoveries in corals belie ongoing change in fish assemblages[END_REF].

Here, we focus on the last COTS outbreak (2007)(2008)(2009)(2010) and cyclone (2010) that occurred in Moorea. Immediately after those disturbances, coral cover in the fore reef was close to 0% but recovered after to levels around 40% similar to those before the environmental disturbances within about a decade (Fig. 1). Filamentous turf algae cover was inversely related to coral cover, and herbivorous fishes (such as parrot-fishes) doubled in biomass after the disturbance, thus controlling macro-algae development on the fore reef.

Coral communities inside the lagoon were not directly affected by the cyclone or COTS outbreak, so limited changes in coral cover and fish communities were observed in the lagoon (Fig. 1). Unlike the fore reef, the lagoon is subjected to land-based and marine anthropogenic drivers. However, over the time period considered here, no clear interactions were identified between these pressures and the ecosystem state [START_REF] Leenhardt | Exploring social-ecological dynamics of a coral reef resource system using participatory modeling and empirical data[END_REF].

Ecosystem model

Three food web models were developed to represent ecosystem states before (2006), during (2010) and after (2016) the environmental disturbance (COTS outbreak and cyclone), respectively. In order to properly assess changes in ecosystem structure and network functioning, models were built with the same structure and area modelled using the Ecopath with Ecosim (EwE) modelling approach (Christensen andPauly 1992, Christensen et al. 2005). Ecopath with Ecosim are powerful frameworks to capture complex ecological processes like those encountered in coral reefs. The Ecopath module of EwE enables to develop a static representation of the food-web in a given ecosystem for a specific period of time. The three Moorea Ecopath models developed in this study were structured using forty-six functional groups (Table . 1). Since the Ecopath model, the static module of the EwE approach, does not explicitly include any spatial dimension, we Ecoystems accounted for the spatial structures in the ecosystem by assigning, where relevant, functional groups either to the fore reef or to the lagoon (see similar approaches in [START_REF] Guénette | Assessing the contribution of marine protected areas to the trophic functioning of ecosystems: a model for the banc d'arguin and the mauritanian shelf[END_REF][START_REF] Piroddi | Modelling the mediterranean marine ecosystem as a whole: addressing the challenge of complexity[END_REF]2017). All benthic species and small territorial fishes were considered separately for the fore reef and lagoon, while larger and often mobile species were assumed to interact between the two habitats.

Each functional group was represented by a production and consumption equation that ensure the mass balance within the system. The production for the group was modelled as function of its total mortality so that:

(1)

where is the production ( ) for the group , the catch ( and detailed information is located in the Supplementary Material Part I.

Incorporating uncertainty into the ecosystem model

Model uncertainty was incorporated into initial parametrization options. For each input parameter (i.e. B, P/B, Q/B and DC ), a log-normal probability distribution was applied following the pre-defined table of the Pedigree Index [START_REF] Pauly | Ecopath, ecosim, and ecospace as tools for evaluating ecosystem impact of fisheries[END_REF][START_REF] Funtowicz | Uncertainty and quality in science for policy[END_REF][START_REF] Christensen | Ecopath with ecosim: a user's guide[END_REF]). The Pedigree module included in EwE allows for a score to be assigned to each individual input variable based on quality of data origin.

According to [START_REF] Christensen | Ecopath with ecosim: a user's guide[END_REF], this score is then converted to a default value of a confidence interval expressed as . For instance, reef fish biomass data, which are locally collected by annual surveys, were assigned a Pedigree score of 5 conventionally assuming an uncertainty of 30% (Supplementary Information Table.6). We then generated, for each period, an ensemble of 1000 Ecopath models, meeting the ecotrophic efficiency balance (EE 1), with input parameters varying within the uncertainty limits, using a Monte Carlo resampling routine. These Ecopath models were used to calculate the output result indicators and their respective confidence limits. This procedure used a suite of Matlab functions provided by [START_REF] Guesnet | Incorporating food-web parameter uncertainty into ecopath-derived ecological network indicators[END_REF].

Analysing models outputs

Trophic spectra analysis

Once food web models are balanced and results are checked for consistency, changes in ecosystem structure and functioning were firstly assessed through trophic spectra i.e. the biomass distribution across trophic levels [START_REF] Gascuel | The trophic spectrum: theory and application as an ecosystem indicator[END_REF][START_REF] Gascuel | Ecotroph: modelling marine ecosystem functioning and impact of fishing[END_REF].

This allowed us to analyse how the overall structure of the fore reef and lagoon sub-systems changed over time. The trophic spectra were calculated , and disturbance models for both habitats (i.e. fore reef and lagoon) using the Ecotroph routine [START_REF] Colléter | An introduction to the ecotroph r package: Analyzing aquatic ecosystem trophic networks[END_REF]) in R and the Ecopath results. For a given model and a given habitat, the trophic spectrum is built including each functional group of that habitat (i.e.

sedentary and mobile groups) and represented as the sum of all groups' spectra. Since no data were available on the relative contribution of mobile groups on each habitat, we accounted for differential contribution according to the total surface area of each habitat.

As fore reef accounts for 10% of the total area of the system, biomass of mobile species were assumed equal to 10% of their total biomass in the fore reef and 90% was assigned inside the lagoon which is inline with observed fish biomass Table .1.

Ecosystem network analysis

Changes in ecosystem structure and functioning between the , and periods were also assessed by comparing trajectories of 12 ecosystem network analysis indicators (Table .2). Here we considered three indicator categories: flow related indicators, path related indicators and structure related indicators.

Flow related indicators:

In this category we computed main flows of the ecosystem related to the Total System Throughput. The Total System Throughput ( ) is measured as the total flow magnitude of a network and has been defined as the size the system [START_REF] Finn | Measures of ecosystem structure and function derived from analysis of flows[END_REF]). Related to we computed the sum of all flows coming to the ( ), lost through the ( ) and coming from the ( ).

The Ecotrophic Efficiency is defined as the fraction of the production of group that is used in a system [START_REF] Christensen | Ecopath with ecosim: a user's guide[END_REF]. The mean Ecotrophic Efficiency ( ) was thus computed as the average of the ecotrophic efficiency of all functional groups with a trophic level (TL) higher than 2.

Path related indicators:

These indicators have been computed to analyse all direct and indirect pathways of the system. The Average Path Length ( ) is the average number of groups through which each inflow passes weighted by the size of the inflows [START_REF] Kones | Are network indices robust indicators of food web functioning? a monte carlo approach[END_REF]. was proposed by [START_REF] Finn | Flow analysis of models of the hubbard brook ecosystem[END_REF] to quantify the cycling capacity of an ecosystem and is computed as the total throughflow that is recycled to the total system throughflow. The to ratio of the communities with TL 2 ( ), averaged at the network scale, informs about the kinetic of trophic flows [START_REF] Gascuel | Trophic flow kinetics in marine ecosystems: toward a theoretical approach to ecosystem functioning[END_REF]). The Efficiency Cumulated Indicator ( ) was recently proposed by [START_REF] Maureaud | Global change in the trophic functioning of marine food webs[END_REF]

Ecoystems

The state of indicators among the three periods were compared using permutational ANOVAs [START_REF] Anderson | A new method for non-parametric multivariate analysis of variance[END_REF]). Significant differences and confidence limits were further considered using the multiple comparison procedure of Tukey-Kramer method based on means comparisons [START_REF] Tukey | Comparing individual means in the analysis of variance[END_REF].

Quantifying the recovery debt

The recovery debt, defined as the interim loss that accumulates during the recovery process was assessed, for each of our ecosystem indicators, at the whole ecosystem scale. (5)

3 Results

Overall changes in trophic structure and functioning

Results of flow diagram allows for two sub-webs to be considered: the lagoon and the fore reef sub-webs, characterized by specific bentho-pelagic pathways, and connected through pelagic pathways (Fig. 2). The highest degree of connection are coming from the pools of detritus and turf which are the largest source of energy for the benthos and herbivorous fishes. Even if corals represent high biomass in the study area, they only marginally contribute to the transfer of energy in the foodweb as few connections stem from them.

The benthic fauna biomass is substantially more represented in the lagoon than in the fore reef. The highest trophic levels are represented by mobile predators (such as Groupers, Jackfish, Large piscivores and Sharks) with trophic values up to 4. Changes between models representing different time periods are mainly observed in the fore reef, with a biomass decrease of coral communities and an increase of detritus and algae during the disturbance, resulting in a significantly higher bottom activity.

Biomass is unequally distributed across trophic levels (Fig. 3). Herbivore and detritivore species (2 TLs 2.5) always hold the largest biomasses in the fore reef and lagoon, while high trophic level communities (TLs 3.5) are less represented. In both habitats the biomass of predators decreases during the disturbance period. In contrast intermediate trophic levels seems to benefit from the disturbance, with increasing biomass.

Nevertheless, this increase does not affect the same trophic levels in the two habitat sub-webs. In the fore reef it is significant between trophic levels 2.5 and 3.0 (mainly omnivores), while in the lagoon the increase during disturbance affects trophic levels ranging from 3.0 to 3.5 (mainly generalist carnivores). After the disturbance the biomass of intermediate TLs seems to have recovered in the fore reef, while biomass of trophic levels between 2.5 and 3 continues to decrease in the lagoon. In both habitats the decrease in biomass is also observed for predators (TLs 3.5).

Changes in network flows

The disturbance of the coral reef ecosystem also affected the ecosystem flows differently in the fore reef and lagoon habitats (Fig. 4). In the lagoon the Flow to Detritus (

), total Dissipation ( ) and Consumption ( ) were larger (+6%, +3% and +2.5% respectively) during the disturbance period compared to before and after the disturbance. Such changes were not observed in the fore reef, with a constant low flow to detritus and constant high rate for dissipation and consumption. The Mean Ecotrophic Efficiency ( ) was oppositely affected in the lagoon (5.4% increasing) and the fore reef (2.9% decreasing) during the disturbance, with a return to its initial value in both systems.

Changes in network path

In the lagoon, the Average Path Length ( ) and the recycling capacity of the system (Finn's Cycling Indice, ) decreased by 21% and 36%, respectively, during the disturbance, with no (or slight) return to previous values after the disturbance (Fig. 5). The total Production over Biomass ( ), remained unchanged in the fore reef but decreased by 20% in the lagoon during the disturbance before returning to its initial value after. The Efficiency Cumulating indice ( ), increased by 10.8% in the fore reef while it decreased inside the lagoon along the disturbance pathway.

Changes in network structure

In the lagoon, the average Mutual Information ( ), measuring the degree of a system organization, increased by 6.5% during and after the disturbance (Fig. 6) while no changes were observed in the fore reef. The pattern was opposite for the Conditional Uncertainty (

), measuring the diversity of the system, with a decrease in the fore reef (-9.7 %) and no change in the lagoon. The Realized Uncertainty ( ), that is the ratio of AMI over H, increased (+4.7%) in the lagoon during and after the disturbance while it increased (+4.4%) only after the disturbance in the fore reef. The System Omnivory Index ( ) responded in opposite ways in the fore reef and in the lagoon with a decrease (-26%) and an increase (+10.7%), respectively, during the disturbance.

Recovery debt

After 6 years of recovery, highest debt is found in the cycling processes (28.85% FCI). The Average Path Length (APL) is 17.43% shorter compared to the pre-disturbance value and the turnover ( ) is 10% slower. The whole ecosystem shows high transfer efficiency debt (9.17%) as well as constraint structure debt (8.61% , 5.58% and 5.54%

).

Discussion

Fish biomass and coral cover indicators, as used in previous ecosystem assessment (e.g [START_REF] Adam | Herbivory, connectivity, and ecosystem resilience: response of a coral reef to a large-scale perturbation[END_REF]2014[START_REF] Lamy | Understanding the spatio-temporal response of coral reef fish communities to natural disturbances: insights from beta-diversity decomposition[END_REF] suggest that the ecosystem has fully recovered from the last disturbances. However, using an ecosystem network approach, we show for the first time, that the studied ecosystem highlights an important recovery debt as evidenced by its changes in ecosystem structural and functional properties.

Combined top-down and bottom-up controls buffer the disturbance in the fore reef functioning

In many ecosystems, stress is expected to alter energy fluxes in food webs [START_REF] Gascuel | Trophic flow kinetics in marine ecosystems: toward a theoretical approach to ecosystem functioning[END_REF][START_REF] Coll | Decadal changes in a nw mediterranean sea food web in relation to fishing exploitation[END_REF][START_REF] Trebilco | Ecosystem ecology: size-based constraints on the pyramids of life[END_REF][START_REF] Maureaud | Global change in the trophic functioning of marine food webs[END_REF]. In coral reefs, this has been evidenced for instance through a shift from long-lived species to high turnover species during a coral-algae transition state [START_REF] Ainsworth | Coral-algal phase shifts alter fish communities and reduce fisheries production[END_REF]). Here we show that, even if biomass and flow changed during the depletion of coral cover, relative flows ( , and ) and cycling indicators ( , and ) in the fore reef were not significantly changed over time suggesting that its functioning remained almost unaffected by the disturbance. Indeed, following the pulse natural disturbance, Moorea's fore reef experienced trophic cascading effects due to combined effects of top-down and bottom-up controls that buffered the effect of the disturbance on the fore reef general functioning. [START_REF] Hempson | Coral reef mesopredators switch prey, shortening food chains, in response to habitat degradation[END_REF] recently highlighted that during an algae transition-state reefs displayed a concave trophic structure characterized by a decline in mid trophic levels (from 2.5 to 3), due to the loss of coral habitat making these mid trophic level organisms more vulnerable to the predation by upper trophic levels. Our trophic pyramid (Fig. 3) shows contrasting results. The decline in coral-dependent fishes were compensated at the same time by an increase in other small territorial fishes (such as planktivorous damselfishes and omnivores fishes) that benefited from the release of predation from the mobile predators and consequently buffer the effect of the habitat loss on that trophic class. As a consequence of this biomass re-organization, fore reef turn over rate remains almost constant. The increase of mobile herbivores with high turn over was compensated at the same time by an increase in planktivores reef fishes with low turn over, and the cycling capacity was also unaffected since detritivorous species (such as macrobenthos and omnivorous fishes) were not reduced.

Trophic cascades following the declines of predators have been observed in many marine environments [START_REF] Estes | Killer whale predation on sea otters linking oceanic and nearshore ecosystems[END_REF][START_REF] Dulvy | Coral reef cascades and the indirect effects of predator removal by exploitation[END_REF][START_REF] Baum | Cascading top-down effects of changing oceanic predator abundances[END_REF]) but a debate exits on the top-down control in coral reefs [START_REF] Ruppert | Caught in the middle: combined impacts of shark removal and coral loss on the fish communities of coral reefs[END_REF], that can be explained by the high diversity of trophic interactions in these systems [START_REF] Thompson | Trophic levels and trophic tangles: the prevalence of omnivory in real food webs[END_REF] weakening possible cascade effects [START_REF] Casey | A test of trophic cascade theory: fish and benthic assemblages across a predator density gradient on coral reefs[END_REF]). Here we add evidence that the removal of predators can lead to an increase in planktivorous dameselfishes, but we also evidence that this top-down effect can also counterpart the bottom-up control due to coral decline.

The off-stage impacts on the lagoon functioning

While the lagoon has no direct impacts of the studied environmental disturbances, our network analysis reveals deep indirect changes in its functioning through time. During the disturbance, biomass of mobile species such as herbivores and invertebrate mobile feeders was considerably increased disrupting its functioning. Contrary to the fore reef where there were changes in bottom composition, the lagoon sub-web did not experience a bottom up re-organization that likely buffered the increase of flows entering the system.

Following the disturbances, the cycling and transfer of biomass in the lagoon were strongly affected. The loss of biomass of omnivores decreased the recycling processes and the food web chain was shortened compared to the before period (as indicated by the Average Path Length). This decrease is likely the result of the increase of the generalist carnivores that feed on both invertebrates and fishes. This generates a stress in the lagoon sub-system that is highly dependent on matter recycling for its overall functioning.

The overlooked impact on the entire ecosystem structure

Although "channelling" networks are increasingly observed in ecosystems enduring anthropogenic stressors [START_REF] Tomczak | Ecological network indicators of ecosystem status and change in the baltic sea[END_REF][START_REF] Ainsworth | Coral-algal phase shifts alter fish communities and reduce fisheries production[END_REF], we show that environmental disturbances can also lead to a simplification of food webs. Following the disturbance, mobile herbivorous species substantially increased in response to the food availability in the fore reef. As most energy are concentrated at low trophic levels, this high biomass production resulted in more organized flow pathways in the fore reef (as indicating by the Conditional Uncertainty and the Realized Uncertainty) and less Ecoystems diversified (as indicating by the System Omnivory Index). The decrease in flow diversity is also likely the result of the lost of habitat structural complexity, with associated loss of coral-dependent species.

Beyond the structural change of the fore reef network, we show that the flows into the lagoon displayed a more constrained organisation during the disturbance (as indicated by Average Mutual Information). Contrary to the fore reef, this degree of organisation was not related to a decreased diversity of flows, since the System Omnivory Index tended to increase during the disturbance (and was reduced in the fore reef). This constrained organisation of flows in the lagoon reefs was likely the result of an increase in mobile fish biomass entering into the system that created new possible pathways especially through the increase of more generalist carnivores that take advantage during the fore reef disturbance.

While previous studies always showed a decrease of ecosystem diversity in regime shifted reefs [START_REF] Ainsworth | Coral-algal phase shifts alter fish communities and reduce fisheries production[END_REF][START_REF] Hempson | Coral reef mesopredators switch prey, shortening food chains, in response to habitat degradation[END_REF], our results strongly suggest than even in case of recovered reefs, diversity of flows can be altered, leading to a more food chain-like system with evident changes at the structural and functioning levels.

The ecosystem recovery debt

We show that the system highlights a recovery debt. The structural and functional traits of the lagoon and fore reef sub-webs did not fully recover, although abundance did [START_REF] Adam | Herbivory, connectivity, and ecosystem resilience: response of a coral reef to a large-scale perturbation[END_REF]2014[START_REF] Lamy | Understanding the spatio-temporal response of coral reef fish communities to natural disturbances: insights from beta-diversity decomposition[END_REF].

We evidence here that disturbance of Moorea fore reef was propagated into the lagoon functioning, with overall recycling capacity reduced about 29% compared to levels prior the disturbance, after 6 years of recovery (Fig. 5 and7). With a reduction in detritivorous species (like most of the benthos), the ecosystem is therefore switching from a productive state, mainly based on detrital pathways to a less productive state with slower nutrient cycling. Detrital pathways are an important source of organic matter entering the coral reef food web that greatly contribute to the transfer efficiency of primary production across the food web [START_REF] Cebrian | Patterns in the fate of production in plant communities[END_REF][START_REF] Cebrian | Variability and control of carbon consumption, export, and accumulation in marine communities[END_REF]. As cycling capacity is still low transfer efficiency was also reduced and remained 9% less efficient in the whole system than before disturbance (Fig. 5). As Efficiency Cumulating Indice only focus on transfers coming from growth and reproductive processes [START_REF] Maureaud | Global change in the trophic functioning of marine food webs[END_REF], our results suggest that communities still suffer higher metabolic losses than before the disturbance. Following the disturbance, higher predation mortality was exercised on communities (especially inside the lagoon), which were stressed and responded with higher respiration. As a result of this predation, small bodied size fishes were replaced during the disturbance by slow growing larger Ecoystems species (like generalist carnivores), resulting in a high debt in ecosystem turnover (about 10% slower).

Structure related indicators revealed that both sub-webs changed from a web-like structure to a more chain-like system, with no apparent recovery. Eleven years after the disturbance, the whole ecosystem has lost [5-8]% of its internal structure flexibility, exhibiting more constrained flows than pre-disturbance state. These loses of ecosystem structural and functional traits may jeopardize future capacity of the system to recover from future natural or anthropogenic disturbance. The disturbance therefore created a disequilibrium between the organized (AMI) and non-organized (DR) part of the system, which can alter its capacity to outlast future disturbances [START_REF] Saint-Béat | Trophic networks: How do theories link ecosystem structure and functioning to stability properties? a review[END_REF].

Conclusion

Here we show that even when ecosystems appear resilient an ecosystem network approach can reveal a recovery debt. As such a debt can pertain to the intrinsic functioning of the ecosystem this can have serious implications for ecosystem services delivery. Those implications can be even stronger in coral reefs where societies highly depend on ecosystem services for their livelihoods. Monitoring programs, impact assessments and resilience based management should not only account for species-based indicators but also track ecosystem network indicators that better reveal the functioning of ecosystems.

Acknowledgements

Service .1 for more details on functional groups.

Table 1: Basic estimates of the three models ( , and the disturbance) for each functional group. Parameters in bold are parameters estimated by the model.

  but that currently shows large signs of recovery. Based on time series, we developed three food-web models , and the disturbance, respectively and assessed trajectories of ecosystem structure and functioning using ecosystem network indicators. The variability of those indicators was accounted for by incorporating uncertainty distributions of model parameters in 1000 simulations for each model.

  to detritus related to the other mortality ( ), with the ecotrophic efficiency of the group i.e. the proportion of the production that is used (consumed,harvested or accumulated) in the system. Consumption for the group ( an Ecopath model requires at least five key inputs parameters: the biomass ( ), production/biomass ratio ( ), consumption/biomass ratio ( ), trophic interactions in a form of a diet matrix ( ) linking functional groups, and fisheries catches. Most of biomass data were obtained from observations (i.e. the CRIOBE annual transect surveys) and production and consumption over biomass ratios were calculated using empirical relationships or using values from the scientific literature on similar ecosystems (Supplementary Information Part I). Diet composition of sedentary fish groups was estimated from a database summarizing 282 dietary information of reef fishes available from the literature (Kulbicki, unpublished data) and adapted to biomass in the habitat of the study. Since mobile predators are using resource on both the fore reef and the lagoon habitats, their diet was divided according to the relative prey density, so that: (3) where is the predation exerting by species on species in the habitat and the biomass of the prey on habitat . Diet composition of benthic invertebrates was taken from Wabnitz et al. (2010). Annual catches were indirectly estimated from local consumption of food fishes assessed from social surveys (Fabre et al. 2018), and were assumed to increase from 2006 to 2016 as demography increased (Insee-ISPF 2007; 2012). Basic input and data sources used to parametrized the models are listed in Table.

  1

  to quantify the fraction of production passing from TL = 2.0 to TL = 4.0. is computed as the product of all partial Production to Consumption ratio within each trophic class of the entire foodweb.Structure related indicators:This category refers to the whole network interactions. The Average Mutual Information () comes from the communication theory and is used to measure the amount of constraint in the network organization[START_REF] Ulanowicz | Quantitative methods for ecological network analysis[END_REF]). 'Constraint' of a system refers to any flow patterns that are well organized, non random. The stability of a network can be measured as the difference between the amount of its uncertainty ( , which can refers to any indirect flow patterns), and of its constraint ( ) and is called the Conditional Uncertainty ( )[START_REF] Latham | Network flow analysis algorithms[END_REF]. The Realized Uncertainty ( ) is the to and refers to the proportion of the total uncertainty[START_REF] Latham | Network flow analysis algorithms[END_REF]. The System Omnivory Index ( ) measures the degree of feeding on several trophic levels and reflects the complexity of a system linkage (Ulanowicz and Abarca-Arenas 1997). These indicators were computed on the basis of the Kones et al. (2009)'s Matlab function. They were calculated separately for the lagoon and the fore reef habitats by considering mobile species as common species in both habitat.

Following

  Moreno-Mateos et al. (2017), we quantified the recovery debt as: (4) with the value of the variable before the disturbance (i.e. the reference value), the recovery period (i.e. from 2010 to 2016) and AUC the Area Under the Curve. AUC were assessed using an exponential approximation. If the value in 2010 exceeded the pre-disturbance value we proceeded to an inverse-transformation as proposed by Moreno-Mateos et al. (2017). The recovery debt was then expressed as percentage:

  dObservation CORAIL from CRIOBE kindly provided the ecological monitoring data. This work was made possible through financial support from ANR (ANR-14-CE03-0001-01) and Fondation de France (INTHENSE). MC was partially funded by the European Unions Horizon research program grant agreement No 689518 for the MERCES project.

Figure 1 :

 1 Figure 1: Biomass and cover of key reef species groups, inside and outside the lagoon for sedentary species and across those two subsystems for mobile species, along the disturbance (COTS outbreak that started in late 2006 and cyclone in 2010). Before, During and After disturbance periods used in subsequent analyses are highlighted. SeeTable.1 for more details on functional groups.

Figure 3 :

 3 Figure 3: Trophic spectra of biomass (in ) for the fore reef and lagoon subsystems , and disturbance.

Figure 4 :

 4 Figure 4: Flow related indicators trajectories along the disturbance. Coloured indicators show significant variations (pvalue 0.05) among the 3 periods. Blue arrows indicate indicators showing signs of recovery.

Figure 5 :

 5 Figure 5: Path related indicators trajectories along the disturbance. Coloured indicators show significant variations (pvalue 0.05) among the 3 periods. Blue arrows indicate indicators showing signs of recovery.

Figure 6 :

 6 Figure 6: Structure related indicators trajectories along the disturbance. Coloured indicators show significant variations (pvalue 0.05) among the 3 periods. Blue arrows indicate indicators showing signs of recovery

Ecoystems

, and the disturbance). Each node represent functional groups within the system and its size is logarithmically proportional to its biomass. Each line represents the trophic interaction among two given groups; the width of the line is scaled to flow proportion and colours are proportional to the magnitude of the flux. The light grey horizontal lines indicate trophic levels of the groups.

Ecoystems

Table 3: Recovery debt per annum estimated at the whole ecosystem scale.


В Е Мизонов

Иван Александрович Балагуров

Cendrine Henri Berthiaux
email: berthiau@enstimac.fr

C Gatumel

Вадим Евгеньевич

Cendrine Gatume

Сыпучий Материал

Vadim E Mizonov
email: mizonov46@mail.ru

Ivan A Balagurov

INFLUENCE OF MULTILAYER LOADING OF DISSIMILAR PARTICULATE SOLIDS INTO A BATCH MIXER ON MIXING KINETICS AND CAPACITY

Keywords: смешивание, сегрегация, многослойная загрузка, цепь Маркова, оптимизация, качество смеси, время смешивания particulate solids, mixing, segregation, multilayer loading, Markov chain, optimization, mixing quality, mixing time

ВЛИЯНИЕ МНОГОСЛОЙНОЙ ЗАГРУЗКИ РАЗНОРОДНЫХ ЗЕРНИСТЫХ МАТЕРИАЛОВ В СМЕСИТЕЛЬ ПЕРИОДИЧЕСКОГО ДЕЙСТВИЯ НА КИНЕТИКУ И ПРОИЗВОДИТЕЛЬНОСТЬ СМЕШИВАНИЯ

Целью настоящей работы является демонстрация нескольких оптимизационных задач периодического смешивания склонных к сегрегации сыпучих материалов, которые могут быть поставлены и решены на основе моделей, основанных на теории цепей Маркова. Для того, чтобы повысить адекватность моделирования и исключить некоторые физические противоречия, которые возникают при использовании линейных моделей, предложена нелинейная Марковская модель перемешивания сегрегирующих компонентов. Исследовано влияние начального

INTRODUCTION

Mixing of powders and granular materials is of central importance for the quality and performance of a wide range of products. It is emphasized in [START_REF] Bridgwater | Mixing of particles and powders: Where next? Particuology[END_REF][START_REF] Bridgwater | Mixing of powders and granular materials by mechanical means A perspective[END_REF] that the design and operation of the mixing process are very difficult, being largely based on judgments rather than science. The next stage of development is to build on the emerging knowledge and methods so that the basics for such design can be laid down. Then this design can become predictable with operation giving effective control of performance. One of the key problems in mixing of dissimilar granular materials is their segregation into each other. The segregation occurs due to differences in physical properties of the components, such as particle size, density, shape, etc. Action of gravity that is always present in mixing is different on different sorts of particles and also leads to their segregation. At no segregation the achievement of homogeneous mixture is the problem of mixing time alone. Very often it is virtually impos-sible to achieve the state of homogeneous mixture if segregation occurs. First, the homogeneity of a mixture increases, then passes its maximum and decreases again. There are a lot of studies, mostly experimental, on influence of the segregation effect on mixture quality (for instance, the papers [START_REF] Jha | Percolation segregation in binary size mixtures of spherical and angular-shaped particles of different densities[END_REF][START_REF] Tang | Segregation quantification of twocomponent particulate mixtures: effect of particle size, density, shape, and surface texture[END_REF][START_REF] Puri | Percolation segregation of multi-size and multi-component particulate materials[END_REF] and others). However, the effect of segregation on mixing kinetics is studied much less.

Modeling the process mathematically helps to understand it better and to search for the ways how it can be improved. However, it is important to choose an appropriate mathematical tool that allows to set up the optimization problems clear and to get optimal solutions that are not too far from possible practical realization.

According to the authors' viewpoint such tool is the theory of Markov chains, which is rather native to the process of mixing because both are related to the evolution of the state of a stochastic system. The basic idea of the Markov chain approach consists in separating the operating volume of the mixer into small but finite zones (cells) and then observing the evolution of the key component concentration in the zones at discrete moments of time after a small but finite time step between them. Application of the approach to describe the state of mixture after passing a static mixer was published in [START_REF] Wang | Axial mixing of grains in a motionless Sulzer (Koch) mixer[END_REF]. However, in this work, the approach missed the evolution of process parameters and did not present the physical features of the mixing zone. In further works by the authors [START_REF] Wang | Stochastic modeling of segregation in a motionless mixer[END_REF][START_REF] Fan | Numerical and experimental simulation studies on the mixing of particulate solids and the synthesis of a mixing system[END_REF], the model, in which transitions only to the neighboring cells were allowed, was developed.

The general strategy of application of the theory of Markov chains to model different processes in powder technology was described in our previous papers [START_REF] Berthiaux | Applications of Markov Chains in Particulate Process Engineering: A Review[END_REF][START_REF] Berthiaux | Application of the theory of Markov chains to model different processes in particle technology[END_REF][START_REF] Berthiaux | Modelling Continuous Powder Mixing by Means of the Theory of Markov Chains[END_REF]. However, these and all other works use the theory namely to describe the process and practically never use the theory to optimize it.

Presented below, the results of the theory application to modeling the kinetic of the binary mixture formation are mainly based on the approach described in [START_REF] Mizonov | Theoretical search for solutions to minimize negative influence of segregation in mixing of particulate solids[END_REF] where it was used to model and optimize the mixing of two dissimilar components of particulate solids. It was mentioned in this work that the structured loading of dissimilar components into a batch mixer can give a considerable gain in comparison to the common loading of them as two separated layers: one component to the bottom of a mixer, and another (segregating downward) one to its top. The detailed study of the approach to improve mixing quality and capacity is the objective of the paper. THEORY Suppose that a batch mixing zone is presented as a one-dimensional array of m perfectly mixed cells. The key segregating component distribution over the cells can be described by the column state vector S = {S j } of the size m×1. As far as the objective of the study is a binary mixture of particulate solids, the downward segregation of a key component always goes together with the upward segregation of a basic component in accordance with the continuity equation. In each cell, the relative content of the basic component is equal to one minus the relative content of the key component. Thus, in order to estimate the mixture quality, it is enough to know how the key component is distributed over the operating volume of a mixer. Namely this distribution and its characteristics will be the objective function of optimization.

Also suppose that we observe the state of the process at discrete moments of time t k = (k -1)Δt where Δt is the transition duration, and k is the transition number that can be interpreted as the discrete analogue of time. In this case, the evolution of the key component state (i.e., mixing kinetics) can be described by the recurrent matrix equation:

k+1 k S = PS , (1)
where P is the matrix of transition probabilities that distributes S j over the cells at each time step, or transition. The j-th column of it belongs to the j-th cell and contains the transition probabilities to go to the neighboring cells and to stay within this cell. These probabilities have the symmetrical part d that is related to the pure quasi-diffusion, which always leads to flattening of the component distribution, and the nonsymmetrical part v that is related to segregation, which leads to non-homogeneity. The values of d and v can be calculated as: d = DΔt/Δx 2 , v = VΔt/Δx where D is the dispersion coefficient, V is the rate of segregation, Δx is the cell height. Thus, the matrix P is a tri-diagonal matrix. In the linear models, the matrix is kept constant from one transition to another. However, it is a very rough assumption to keep the value of v constsnt. It can lead not only to quantitative errors but also to qualitative contradictions. Indeed, a segregating component cannot segregate inside itself, and, if the (j+1)-th cell is filled with it fully (contains S max of it), the value of v should be equal to zero. On the contrary, if the (j+1)-th cell does not contain it at all, the value of v should be maximal. It is natural to introduce the assumption that the dependence between v and S j+1 k is linear:

k j+1 k j 0 max S v = v (1 - ) S , (2)
where v 0 is the rate of segregation into the pure basic component and S max can be supposed to be equal to 1. Thus, the model becomes non-linear because the matrix P depends on the current state vector S k and varies at each time transition. Under this assumption the entries of the matrix P can be calculated as:

k k j, j+1 0 j+1 P = v (1 -S ) + d , j=1,…,m-1 (3) k j+1, j P = d , j=1,…,m-1 (4) m k j, j i, j i=1,i j P = 1 - P , i=1,…,m. (5)
The technological objective of mixing of solids is to obtain the mixture, in which a key component is distributed over the mixture volume as homogeneously as possible, or meets technological requirements on the mixture non-homogeneity. If the nonhomogeneity is characterized by the standard deviation of a key component distribution σ, this condition can be formulated as σ→σ min , or σ ≤ σ c . If the parameters of the mixture agitation in a mixer are given (i.e., the values of d and v 0 are defined), the only way to optimize the process is to find the optimal mixing time, at which σ = σ min . Under these conditions, there is no ways to influence on the matrix P but it is possible to influence on the initial state vector S 0 that is being formed during the loading of components into a mixer.

The common way to begin the mixing process is to place a segregating key component to the top of a mixing zone above a basic component and to run the mixer, i.e., to introduce agitation of the mixture by mechanical action (vibration, blades, etc.). However, such one-time loading does not allow reaching the homogeneous mixture, and the minimum standard deviation very often is higher than it is necessary due to technological requirements. Formally, there exists the theoretical way to reach the homogeneous mixture of segregating components. It is necessary to take from the components particle by particle and to put them into the mixing volume in the proper order. However, the capacity of such mixing process is close to zero. A compromise settlement is to load the com-ponents into the mixing volume layer by layer like sandwiches, in which the number layers of the components keep the same proportion as in the total mixture. In this case, the process of loading becomes a sort of pre-mixing, and only short time post-mixing is required after that in order to reach the necessary homogeneity.

In order to estimate the efficiency of multilayer loading the numerical experiments were carried out with the model described above. The results of them are shown in Fig. 1. The number of cells was taken equal to 48, i.e., one cell was the scale of the probe to measure the content of the key component, which total content in the mixture was 1:5. The calculations were done for v 0 = 0. The case 1 is related to the one time loading. The best mixing quality is being reached after 162 time transitions when σ min = 0.0824. The cases 2, 3 and 4 are related to the multilayer loading as 2, 4 and 8 sandwiches respectively. It can be seen that the growth of the sandwiches number leads to the considerable decrease of σ min (the increase of the mixture quality) and optimum mixing time k opt (the increase of mixing capacity). If we compare the key component distribution at the optimum mixing time for the one time loading (Case 1) and the same for the loading as 8 sandwiches (Case 4) we can see a big difference. At one time loading the best mixture that can be obtained remains rather non-homogeneous, while it is close to the perfect one at loading as 8 sandwiches. It is also important that this mixture can be reached more than 10 times faster. Thus, the more layers we have in the mixer load, the better. Despite of the fact that this conclusion is more or less obvious, it gets the quantitative estimation on the basis of the proposed model of mixing.

However, the process of multilayer loading take more and more time with the growing number of layers, or sandwiches that obviously influences on the total capacity of mixing. Let us suppose that the loading time (measured in the number of time transition of duration Δt) is directly proportional to the number of sandwiches J. In this case, the duration of the mixing cycle "loading plus mixing" K can be calculated as: K=k mix +βJ, [START_REF] Wang | Axial mixing of grains in a motionless Sulzer (Koch) mixer[END_REF] where k mix is the mixing time proper, β is the proportional coefficient. In order to compare the different cases of loading shown in Fig. 1, let us also suppose that the mixing homogeneity σ min = 0.0824, which can be reached at one-time loading, meets the technological requirement, and the only objective of multilayer loading is to increase the mixing capacity. The required mixing times for different multilayer loading can be easily found for each case from the kinetic curves in Fig. 1.

Fig. 2 illustrates how the total mixing time K varies with the number of sandwiched at different value of β.

It can be seen from the graphs that at very fast loading (β=1) the best case is to have 8 sandwiches while at slow loading (β = 20) 2 sandwiches are preferable. The optimum number of sandwiches exists depending on β.

EXPERIMENTAL

The objective of this section is not to present the results of detailed experimental study of a mixer with multilayer loading but just to demonstrate an example that the multilayer loading can be realized in practice. One of possible technical solutions of multilayer loading is presented schematically in Fig. 3.

A movable cylindrical basket separated by vertical partitions into two sections proportionally to the volume of components to be mixed is inserted into the cylindrical mixing chamber of a vibration mixer. The basket has a perforated bottom but also can be without the bottom at all. The basket has the possibility to rotate around the vertical axis and to move up in the vertical direction. Actually, the partitions play the key role, and they can be adjusted to the vertical shaft even without the basket. The components are loaded into the sections formed by the partitions up to the identical level. After that the basket is brought in motion. It rotates and, at the same time, moves up leaving the components layer by layer on the bottom of the mixing chamber. After all the material lefts the basket, we get as many sandwiches as many turns the basket does during its motion to the upper point. Then the basket should be moved out, and the formed multilayer mixture should get a short time vibration agitation to finish mixing.

In order to check such loading process in practice, the experimental set-up was made. Its photograph is shown in Fig. 4a. The mixing chamber was imitated by the transparent glass, and the partitions were installed at the vertical shaft without basket. The shaft was driven by two independent electrical motors with regulated frequency of rotation placed on the frame. The fine ground salt and buckwheat in equal proportion were used as the test materials. Fig. 4b gives an example of the layers formation. It can be seen that the proposed device for multilayer loading runs the process rather well. After placing the mixing glass on a vibration table and applying a short time vibration action the mixture becomes very close to the homogeneous one (like for the case 3 or 4 in Fig. 1).

CONCLUSIONS AND PERSPECTIVES

It is shown that transition from one-time loading of dissimilar particulate solids in a batch vibration mixer to multilayer loading of them is a considerable reserve of improving the mixing quality and productivity. It is also shown that if the loading time and mixing time are taken into account together, the optimum number of layers exists that provide the maximum productivity under mixing quality being equal. The technical possibility to realize the multilayer loading is demonstrated experimentally that opens the ways to improve the design of vibration mixers. According to our viewpoint, the further development of the approach lies in transition to ternary and more components mixtures in modeling and to complete experimental validation of such process at a pilot scale mixer.

  2 and d = 0.3. The initial distributions of the key component are shown at the top of the graph, the mixing kineticsin the middle, and the distributions of the key component after mixing during the optimal mixing time are shown at the bottom of it.

Fig. 1 .

 1 Fig.1. Influence of initial distribution on mixing kinetics and component distribution at maximum homogeneity of the mixture Рис. 1. Влияние начального рапределения на кинетику смешения и распределение компонентов при максимальной однородности смеси

Fig. 2 .

 2 Fig. 2. Influence of the number of sandwiches on the total mixing time at different value of β Рис. 2. Влияние количества сандвичей на общее время смешения при различном значении β

Fig. 3 .Fig. 4 .

 34 Fig.3. Schematic presentation of loading device for multilayer loading Рис. 3. Схема загрузочного устройства при многослойной загрузке

ACKNOWLEDGEMENT

This work is supported by the Russian Foundation for Basic Research (project 15-08-01684).


Lucie Lefort

Romain Bonabe De Rougé

Pascal Schetelat

Thomas Berthou

Peter Riederer

Bruno Duplessis

Eric Peirano

Development of a methodology of validation for urban building energy models and application to French residential consumption

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

Introduction

With an urgent need to act against climate change, more and more cities set ambitious environmental targets. In order to achieve such goals, reducing energy consumption in buildings seems to be a keystone. In France, for example, tertiary and residential sectors represented 45% of the national final energy consumption in 2015 (ADEME, 2015). To assist in lowering these consumptions at the city scale, urban building energy models (UBEM) have been developed over the last fifteen years [START_REF] Reinhart | Urban building energy modeling -A review of a nascent field[END_REF]. UBEM compute the energy consumption of every building in a district and take into account interactions among them, like, for example, solar shading. The aim of UBEM is to provide cities with guidance to reduce final energy consumption not only for some specific buildings but also at the scale of a whole district. However, the benefits brought by UBEM are balanced with certain challenges that may jeopardize their reliability. A first challenge is the data enrichment process. It is much more difficult to gather data required for the simulation of one hundred buildings than for a single one. Information is often missing and assumptions must then be made. This practice of data enrichment creates an error that is not yet well-quantified. A second challenge is to take into account new phenomena that appear at the district scale, like the phenomenon of urban heat islands. Finally, although a building energy model (BEM) has to be precise for each building it simulates, this is not necessarily the case for a UBEM. It may be acceptable for a UBEM to be statistically correct at the scale of a district but to be rather inaccurate at the scale of a single building. The challenge is then to identify the spatial and time scales on which the UBEM is reliable. A thorough validation process would enable to verify the reliability of UBEM and would give them more credits as decision-support tools for cities. However, such validation process is not yet implemented within the UBEM research community and is therefore needed. The aim of this paper is to present a validation methodology designed specifically for UBEM. This methodology is described in the next section before being applied to two UBEM tools in a second section. Results obtained by the two UBEM tools throughout the validation process are then being discussed.

Development of a validation methodology for urban building energy models

Most UBEM stem from building energy models (BEM). However, at the difference with the BEM, uSIM2020 -Building to Buildings: Urban and Community Energy Modelling, November 12th, 2020 there is no consensus yet on their validation. In this section, we review the existing validation attempts that can be found in the literature and we draw a list of requirements for a validation methodology. Then, we describe a first proposal of validation methodology based on UBEM use cases.

Existing attempts of validation

BEM validation

Building energy model can be seen as a cornerstone of UBEM. Therefore it seems natural to look at the existing validation process implemented for BEM. A lot of research has been carried on BEM validation among which were the BESTEST and the PASSYS project.

In the BESTEST [START_REF] Judkoff | Twenty years on!: Updating the iea bestest building thermal fabric test cases for ashrae standard 140[END_REF], several sets of perfectly described buildings are provided and inter-comparison of models are computed on these well-known buildings. In addition, some empirical validation has been proposed with measurements on test cells. [START_REF] Jensen | Validation of building energy simulation programs: a methodology[END_REF] gives a precise description of the methodology developed within the PASSYS project. The proposal is to split up the validation process into two main parts: the validation of single processes and the validation of the whole model. Jensen encourages namely to use analytical verification for single process validation and to use empirical validation for whole model validation. The validation of thermal models used within UBEM through one of these two validation methods would not suffice to consider a UBEM tool as properly validated for several reasons:

Additional physics and additional energy systems are modelled at the urban scale; An enrichment process is often used as a prior to UBEM and needs to be validated as well; It may be acceptable for UBEM to be accurate only at certain spatial scales, larger than the building scale.

However, some principles from these methodologies should be used as guidance for the development of a UBEM validation methodology. For example, separating single process validation and the whole model validation as advised by [START_REF] Jensen | Validation of building energy simulation programs: a methodology[END_REF] is of great interest for UBEM as well. Furthermore, in the same idea as the PASSYS project, the aim of the present work is to develop a validation methodology rather than validating a single tool.

Regarding validation attempts on UBEM, two approaches have been observed so far: comparative validation and empirical validation.

Comparative validation

Inter-model comparison has started within two international research projects, the IEA Annex 60 (Wetter, 2017) and the IBPSA Project 1 [START_REF] Saelens | Towards a DESTEST: a District Energy Simulation Test Developed in IBPSA Project 1[END_REF] which is still ongoing. In both projects, the taken approach is quite similar to the BESTEST [START_REF] Judkoff | Twenty years on!: Updating the iea bestest building thermal fabric test cases for ashrae standard 140[END_REF]. A synthetic district has been designed to serve as a District Energy Simulation Test (DESTEST). Several UBEM tools are being tested on this DESTEST and their results are compared. This setup enables to control all the simulation parameters and thus proceed with the single process validation as advocated by [START_REF] Jensen | Validation of building energy simulation programs: a methodology[END_REF]. However, since the district used is not real, it is not possible to compare the obtained results with measured data nor to test the enrichment process used by each UBEM tool.

Empirical validation

The aim of empirical validation is to ensure that UBEM tools represent the real-life urban projects accurately. However, such process is made difficult due to a lack of data at district and urban scales. Some examples can be found in the literature like in [START_REF] Sokol | Validation of a Bayesian-based method for defining residential archetypes in urban building energy models[END_REF] or in [START_REF] Fonseca | Integrated model for characterization of spatiotemporal build-ing energy consumption patterns in neighborhoods and city districts[END_REF]. These attempts give a first idea of UBEM performances on real projects. However, datasets used in these validation processes are rarely publicly available thus preventing other research teams from going through the same validation process. Furthermore, there is no agreement so far on the metrics to be used to measure the error on UBEM outputs. This prevent the reader from comparing results obtained in different papers. Nonetheless, some requirements can be drawn from these existing attempts of UBEM validation.

Requirements for a UBEM validation methodology

Some lessons can be learned from this brief literature review and it can be summarized as follows. A validation methodology dedicated to UBEM should comply with the following requirements: Specifying the spatial and time scales on which a given UBEM tool is reliable; Defining metrics to properly measure the UBEMs accuracy throughout the validation process;

Identifying the error generated by the model setting;

Ensuring applicability and replicability of the methodology in different contexts (country, building sector etc.); Proposing test cases to apply the developed validation method.

A first proposal of methodology is presented hereafter based on these requisites.

A validation method based on UBEM's use cases

The methodology developed here offers to inventory UBEMs use cases and to identify, for each use case, critical UBEM outputs. The aim is to verify that the outputs of interest for each use case are accu-uSIM2020 -Building to Buildings: Urban and Community Energy Modelling, November 12th, 2020 rate enough to provide a real aid to the decisionmaker. Such outputs of interest are called key performance indicators (KPI) in the rest of the paper. Once use cases and KPI have been identified, validation datasets are used in order to study, for each use case, the error distribution on the KPI. This methodology is described in Figure 1.

In an ideal case, several batches of measured data would be available, enabling to obtain an error distribution for each KPI of each use case and UBEM tool. The main idea behind this search for an error distribution is to provide the UBEM user with a confidence interval on the [START_REF] Kontokosta | A data-driven predictive model of city-scale energy use in buildings[END_REF].

Application of the proposed methodology to two UBEM tools

Validation framework

A first use case was chosen to apply and test the proposed validation methodology. This use case can be described as territorial energy diagnosis. It consists in being able to give an accurate outlook of a district energy consumption. In this study, the annual electricity and gas consumptions of residential buildings were studied at the district scale. This use case may be seen as a preliminary study for other use cases since a first step for UBEM is to provide accurate representation of existing consumptions before evaluating different scenarios like the retrofitting of buildings or the development of renewable energies. The choice of this use case was also driven by the publicly available datasets. This use case was applied on the French department of Seine-et-Marne at the scale of the IRIS. The IRIS is an administrative division of the French territory.

Each IRIS gathers around 2000 inhabitants. All residential buildings located in about a hundred of IRIS were simulated over a year and their annual electricity and gas consumptions aggregated at the IRIS scale were compared to measured data. Simulations were completed with two different UBEM tools, Dimosim [START_REF] Riederer | Development of a simulation platform for the evaluation of district energy system performances[END_REF] and Smart-E [START_REF] Berthou | SMART-E: A TOOL FOR ENERGY DEMAND SIMULATION AND OPTIMIZATION AT THE CITY SCALE[END_REF], both relying on automated data enrichment processes, and results were compared. This case study deals only with the residential sector since less databases are available for the simulation of the tertiary sector in France. Datasets used for the simulation and for the validation are presented hereafter.

Simulation datasets

Three databases were used for simulation: BD TOPO® (IGN, 2019) This database provides a map of the whole French territory where each building is described by a footprint and a height. Additional information is also given per building such as:

-Building function (residential, commercial, etc.) -Construction year -Number of floors -Number of dwellings Though the footprint and the height are given for every building, the additional information is not always available. This additional information is obtained through a joint between tax record files and the map. Information is missing when the joint has not been successful. The BD TOPO® is not publicly available yet but it can be granted for free for any research project. It should become completely public by 2022. PHEBUS (Conseil National de l'Information Statistique, 2013): This database is built upon a survey completed in 2013 for 5000 households and for each household, their appliances and their energy bills were reported. Following this first step, energy performance audits were completed in 2000 out of these 5000 households. Fichiers Détails logements (FDL) (INSEE, 2016): This database originates from the population census. Each row corresponds to 1 to -Dwelling's number of occupants; -Socio-professional status of the household; -Surface category of the dwelling; -Type of dwelling (individual housing or collective housing); -Construction year of the building; -Main heating energy; -Main heating system type (centralized or decentralized).

Validation dataset

The validation dataset employed for this use case is the Données Locales de l'Energie database (SDES, 2020). This database originates from a French law for the ecological transition that makes compulsory for energy suppliers to disclose annual energy consumptions data aggregated at the IRIS scale and detailed by sector. Annual electricity and gas consumption of the residential sector are thus available at the IRIS scale. These data are publicly-available and are updated every year. The related disclosure law was recently amended compelling energy suppliers to publish energy consumption data at the building scale, given some privacy threshold for the residential sector. These datasets requires some preprocessing but are of great interest for the validation of UBEM.

Two different UBEMs, one validation process

The proposed validation method aims to be applicable to a large range of UBEM. Two first UBEM tools were used here to test this method: Dimosim [START_REF] Riederer | Development of a simulation platform for the evaluation of district energy system performances[END_REF] and Smart-E [START_REF] Berthou | SMART-E: A TOOL FOR ENERGY DEMAND SIMULATION AND OPTIMIZATION AT THE CITY SCALE[END_REF]. These tools are very different from each other.

A first difference comes from the data enrichment process. Smart-E uses internal hypotheses to complete missing data while Dimosim is coupled with a tool designed for the enrichment process: Qiriel (Ansanay- [START_REF] Ansanay-Alex | Statistical and Stochastic Modelling of French Households and their Energy Consuming Activities[END_REF].

A second difference comes from the datasets used as inputs data and as parameters. Dimosim takes as input data the geometry of the buildings retrieved from BD TOPO® and then tries to assess thermal properties and household characteristics thanks to PHEBUS and FDL. Conversely, Smart-E takes as input data the description of households provided by FDL that includes the surface of the dwellings but not the geometry of the buildings and then tries to allocate façade areas and roof areas thanks to BD TOPO®. BD TOPO® and FDL both give information on dwellings but since FDL comes from the census, it is anonymized to prevent identifying the households and the buildings. Therefore, the link between these two databases is not straightforward and assumptions must be made. Furthermore, information provided by these databases is not always consistent. For example, the number of dwellings recorded in BD TOPO® and FDL is different. A third difference is the number of thermal zones defined for the simulation. Several thermal zones per building can be defined in Dimosim but for the sake of simplicity only one thermal zone per building was used in this case. Smart-E uses two thermal zones per dwelling by default, one for the heated space and one for the non heated space. [START_REF] Booth | Handling uncertainty in housing stock models[END_REF] identified the most influential parameters on housing stock modelling uncertainties.

The mean value of these parameters distributions used for each UBEM tool in this study are summarized in Table 1. These values are the average on all the simulated dwellings weighted by the dwellings surfaces. These mean values have been displayed to illustrate the difference of parameterization between the two tools. For example, the share of dwellings heated with gas and the share of dwellings heated with electricity vary from an IRIS to another and are often slightly different between Dimosim and Smart-E due to different enrichment processes.

Results & Discussion

KPI and metrics

The proposed methodology was tested on 93 IRIS (around 90000 dwellings), all located in the French department of Seine-et-Marne. The studied KPI in this case were the annual electricity and gas consumptions at the IRIS scale but also the annual thermosensitivity for gas and electricity at the IRIS scale.

The thermosensitivity is defined as the increase of electricity or gas consumption due to the decrease of the outdoor temperature. Its estimation for each tool was done following the methodology provided with the DLE [START_REF] Enedis | Estimation des données de thermosensibilité et de part thermosensible[END_REF]. Measured data were available in the DLE for each of these KPI and for each IRIS. The error on annual consumption and thermosensitivity was measured through a percentage of error (PE) defined in Equation 1.

P E X = X Simulated -X measured X measured • 100 (1)
The results are presented in Figures 2 and4.

Analysing the error distribution with respect to input data

The error distribution obtained with Dimosim on annual consumptions is shifted toward the positive values, meaning that Dimosim mainly overestimates the electricity and gas consumptions. Conversely, Smart-E tends to underestimate these consumptions.

In absolute values, Dimosim seems to generate a greater error than Smart-E on this particular use case. This performance gap could come from the fact that Smart-E is calibrated in order to fit national electricity consumptions while there is no similar process in Dimosim.

Regarding the error on the thermosensitivity, it seems that Dimosim also overestimates this phenomenon for both gas and electricity. An hypothesis is then that Dimosim overestimates the heating needs in general.

When comparing inputs parameters between Smart-E and Dimosim, the fraction of space heated looks like a good candidate to help reducing this error. In order to confirm this hypothesis, simulations should be run again with a different space-heated ratio and results should then be compared. When taking a closer look at the error generated by Smart-E, the error on the electrical thermosensitivity seems surprisingly high compare to the error on the annual electricity consumption. In order to get a better understanding of all these error distributions and to try to explain phenomenon like the one observed on the thermosensitivity for Smart-E, an analysis of the correlation between the error and a large set of inputs parameters have been completed. Correlations were estimated through the Pearson correlation coefficient. This coefficient is defined in equation 2 for two variable X and Y, where cov(X,Y) is the covariance of X and Y, σ X is the standard deviation of X and σ Y is the standard deviation of Y. A coefficient of 1 indicates a linear relationship between variables while a coefficient close to 0 indicates that there is no linear relationship (but a non-linear relationship is still possible between variables).

ρ X,Y = cov(X, Y ) σ X • σ Y (2)
Coefficient correlations between the percentage of error on annual gas and electricity consumptions, gas thermosensitivity and electrical thermosensitivity and inputs parameters are presented in Figure 3.

In general, the obtained coefficients are quite low, most of them being below 0.5 in absolute value. It is then difficult to draw influential parameters from this study. However, some trends can be identified. Thus, the two UBEM tools seem to behave quite differently in terms of error. Smart-E presents a coefficient of 0.64 between the percentage of error on the gas consumption and the one on the electricity consumption but such correlation does not appear for Dimosim. Furthermore, the higher coefficients can be found for Smart-E between the percentages of error (gas and electricity) and the window and wall U-values. In Dimosim, a relatively high correlation is observed be- tween the percentage of error on electricity and the share of dwellings using electricity as heating energy.

Finally, a mild correlation between Dimosim and Smart-E errors was identified as well, questioning the quality of databases used for simulation and validation and the error they might bring in. Thus, this analysis enabled to identify levers to reduce the overall error of both tools.

uSIM2020 -Building to Buildings: Urban and Community Energy Modelling, November 12th, 2020

Conclusion

A methodology of validation designed specifically for UBEM was presented in this paper. This method aims to verify that UBEM are accurate enough to provide a real aid in the decision-making process. To do so, the validation procedure is broken down based on UBEM use cases and error distributions are studied for KPI of each use cases. This method was tested with two first UBEM tools on a French test case. The analysis of the obtained error distributions enabled to identify weaknesses of both UBEM tools and hence improve their performances. The development of this methodology is still ongoing and future work should be carried on UBEM outputs at smaller spatial and time scales. The error analysis and the identification of critical inputs parameters in the output error generation should also be strengthen.

Lastly, an important step would be to verify the applicability of the methods to a large range of UBEM tools. Since the available datasets for simulation influence greatly the enrichment process and UBEM behaviour in general, it would be of great interest to apply this methodology in other countries than France and with other UBEM tools as well. The authors would then be very interested in collaborating with other research teams on this topic.

Figure 1 :

 1 Figure 1: Proposed validation methodology

  Figure 2: Annual measured and simulated consumptions of gas and electricity

Figure 4 :

 4 Figure 4: Percentage of error on annual energy consumptions and on thermosensitivity

Table 1 :

 1 Influential parameters on results uncertainty

	Parameters			Dimosim Smart-E
	Fraction of space heated	100 %	77 %
	Heating temperature set-point (C)	19.6	19.6
	Share of dwellings with gas as heating energy	48.3 %	50.6%
	Share of dwellings with
	electricity as heating en-	38.2 %	35.7 %
	ergy
	Efficiency heaters	of	electric	1.0	1.0
	Efficiency of gas boilers	0.78	0.69
	Window (W/m 2 K)	U-value	3.21	2.69
	Wall U-value (W/m 2 K)	0.96	1.38
	Window-to-wall ratio	0.25	0.15
	Air-change rate (including
	infiltrations and ventila-	0.6	0.52
	tion) (vol/h)
	Number of simulated dwellings	89876	93717
	Average dwelling (m 2 ) surface	per	105	92