Gabriela Zanghelini
email: gabriela.zanghelini@agroparistech.fr

Martine Esteban-Decloux

Stéphane Vitu

Pierre Giampaoli

Violaine Athès

ß-Damascenone highly diluted in hydroalcoholic mixtures: phase equilibrium measurements, thermodynamic modeling and simulation of a batch distillation

Keywords: vapor-liquid equilibrium, batch distillation, simulation, ß-damascenone, NRTL, alcoholic beverages

ß-Damascenone is a carotenoid-derived compound that plays a key role in the characteristic aroma of wines and distilled spirits. Despite the valuable contribution of this compound to the flavor quality of alcoholic beverages, little is known about its volatility and its vaporliquid equilibrium behavior in hydroalcoholic mixtures of different compositions. Such information is essential to be able to choose the right thermodynamic model for ßdamascenone when using a distillation process simulation software to gain a better understanding of the behavior of ß-damascenone during spirit distillation. In this context, this work presents new experimental VLE data for ß-damascenone highly diluted in ethanol + water mixtures. The measurements were carried out in a recirculation ebulliometer operating at 101.3 kPa for ethanol mass fractions between (0.10 and 0.79) g•g -1 . Activity coefficient model NRTL was used to correlate the experimental data and binary interaction parameters. The parameters obtained for the NRTL model were employed to simulate a laboratory-scale batch distillation of a model solution of ß-damascenone, ethanol and water.

The simulation results were well fitted to experimental data. This approach can be further extended to other compounds that contribute to the aroma of distilled spirits and the results can be transferred to simulate larger scale distillations.

INTRODUCTION

Spirit distillation holds the triple role of concentrating volatile aroma compounds from a fermented medium into an ethanol-rich distillate, selecting specific aroma fractions and generating new aroma compounds via chemical reactions [START_REF] Lurton | Production and Aromatic Characteristics[END_REF] . Despite their extremely low concentrations in the hydroalcoholic matrix (with mass fractions between 10 -10 and 10 -3 ), aroma compounds play a primary role in the intricate aroma typicity of distilled beverages.

The contribution of each compound to the overall aroma is dictated by intrinsic factors such as their partition between the hydroalcoholic matrix and the headspace, as well as physical, chemical and structural properties of the different components in the system and the interactions between them [START_REF] Ammari | Flavor Retention and Release from Beverages: A Kinetic and Thermodynamic Perspective[END_REF] . Thence, achieving a desired aroma bouquet in a spirit requires fine-tuning its composition in aroma compounds, which in turn entails controlling the distillation step.

Thermodynamically, alcoholic spirits can be approached as non-ideal multicomponent mixture of aroma compounds highly diluted in an ethanol -water solution. During distillation, compound separation is assured by deviations from a state of vapor-liquid equilibrium (VLE) prompted by a continuous addition of energy through heat. The thermal properties of the system are dictated by the two major components, ethanol and water, whereas the behavior of the minor components in the hydroalcoholic mixture throughout distillation depends on their pure compound properties, i.e., vapor pressure and enthalpy of vaporization, as well as activity coefficients at different ethanol -water compositions, which lead to different volatility values depending on the composition.

Aroma compounds in distilled spirits belong to several chemical families, namely alcohols, esters, carboxylic acids, norisoprenoids, terpenes, aldehydes and ketones. Amongst these families, norisoprenoids are a class of highly odorant carotenoid-derived compounds that greatly contribute to the varietal aroma of wine 3,[START_REF] Ebeler | Wine Chemistry and Flavor: Looking into the Crystal Glass[END_REF] and are often listed among the most odorant components in distilled spirits. Within norisoprenoids, ß-damascenone is a particularly potent odorant, with odor thresholds as low as 0.002 µg•L -1 in water [START_REF] Buttery | Quantitative and Sensory Studies on Tomato Paste Volatiles[END_REF] , 7 µg•L -1 in wine [START_REF] Ickes | Effects of Ethanol on Flavor Perception in Alcoholic Beverages[END_REF] and 0.05 and 10 µg•L -1 in model ethanol-water solutions at, respectively, 12% v/v ethanol [START_REF] Pineau | Which Impact for β-Damascenone on Red Wines Aroma?[END_REF] and 40% v/v ethanol [START_REF] Masuda | Occurence and Formation of Damascenone, Trans-2,6,6-Trimethyl-1-Crotonyl-Cyclohexa-1,3-Diene, in Alcoholic Beverages[END_REF] . It is a natural component of grape which has been reported to impart fruity, apple-like notes in various wines, where it can be found in its free form or bound to glycosylated moieties at concentrations between 0.1 and 10 µg•L -1 9 . Moreover, distillation has been shown to substantially increase the concentrations of ß-damascenone, likely prompted by the high temperatures and acidic conditions involved in the process [START_REF] Zierer | Aroma-Active Compounds in Bartlett Pears and Their Changes during the Manufacturing Process of Bartlett Pear Brandy[END_REF][START_REF] Schreier | Changes in the Composition of Neutral Volatile Components during the Production of Apple Brandy[END_REF][START_REF] Kotseridis | Quantitative Determination of Free and Hydrolytically Liberated β-Damascenone in Red Grapes and Wines Using a Stable Isotope Dilution Assay[END_REF][START_REF] Franitza | Influence of the Production Process on the Key Aroma Compounds of Rum: From Molasses to the Spirit[END_REF][START_REF] Awad | Evolution of Volatile Compounds during the Distillation of Cognac Spirit[END_REF] .

As a result, and given its aforementioned low odor thresholds, ß-damascenone is listed amongst the key odorants responsible for the characteristic aroma of rum [START_REF] Franitza | Characterization of the Key Aroma Compounds in Two Commercial Rums by Means of the Sensomics Approach[END_REF][START_REF] Pino | Characterisation of Odour-Active Compounds in Aged Rum[END_REF] , apricot and apple distillates [START_REF] Genovese | Comparison of the Aroma Compounds in Apricot (Prunus Armeniaca, L. Cv. Pellecchiella) and Apple (Malus Pumila, L. Cv. Annurca) Raw Distillates[END_REF] , spine grape spirit [START_REF] Xiang | Characterization of Odor-Active Compounds in the Head, Heart, and Tail Fractions of Freshly Distilled Spirit from Spine Grape (Vitis Davidii Foex) Wine by Gas Chromatography-Olfactometry and Gas Chromatography-Mass Spectrometry[END_REF] , whisky [START_REF] Poisson | Characterization of the Key Aroma Compounds in an American Bourbon Whisky by Quantitative Measurements, Aroma Recombination, and Omission Studies[END_REF] , tequila [START_REF] González-Robles | The Impact of Maturation on Concentrations of Key Odour Active Compounds Which Determine the Aroma of Tequila[END_REF][START_REF] Lahne | Streamlined Analysis of Potent Odorants in Distilled Alcoholic Beverages: The Case of Tequila[END_REF] , cachaça [START_REF] De Souza | Characterization of Cachaça and Rum Aroma[END_REF] , cognac [START_REF] Uselmann | Decoding the Combinatorial Aroma Code of a Commercial Cognac by Application of the Sensomics Concept and First Insights into Differences from a German Brandy[END_REF][START_REF] Ferrari | Determination of Key Odorant Compounds in Freshly Distilled Cognac Using GC-O, GC-MS, and Sensory Evaluation[END_REF] and calvados [START_REF] Ledauphin | Chemical and Sensorial Aroma Characterization of Freshly Distilled Calvados. 2. Identification of Volatile Compounds and Key Odorants[END_REF] , and is thought to act as an aroma enhancer in wine [START_REF] Pineau | Which Impact for β-Damascenone on Red Wines Aroma?[END_REF][START_REF] Escudero | Analytical Characterization of the Aroma of Five Premium Red Wines. Insights into the Role of Odor Families and the Concept of Fruitiness of Wines[END_REF] and distilled spirits such as cognac [START_REF] Thibaud | Contribution of Volatile Odorous Terpenoid Compounds to Aged Cognac Spirits Aroma in a Context of Multicomponent Odor Mixtures[END_REF] . In light of the above, reaching a better understanding of how the distillation operation affects the behavior of ß-damascenone would be valuable to gain a better control of the final aroma of a distilled spirit.

Process simulation software are promising tools to understand, design and optimize distillation processes for alcoholic spirits. It has been consistently applied over the years to represent industrial scale continuous distillations of cachaça [START_REF] Batista | Computer Simulation Applied to Studying Continuous Spirit Distillation and Product Quality Control[END_REF] , neutral alcohol [START_REF] Decloux | Simulation of a Neutral Spirit Production Plant Using Beer Distillation[END_REF][START_REF] Valderrama | Advances on Modeling and Simulation of Alcoholic Distillation[END_REF][START_REF] Batista | A New Distillation Plant for Neutral Alcohol Production[END_REF][START_REF] Esteban-Decloux | Two Industrial Examples of Coupling Experimentsand Simulations for Increasing Quality and Yieldof Distilled Beverages[END_REF] , armagnac [START_REF] Puentes | Simulation of Spirits G. Zanghelini et[END_REF] and whisky [START_REF] Gaiser | Computer Simulation of a Continuous Whisky Still[END_REF] , as well as batch distillations of cognac [START_REF] Douady | Batch Distillation of Spirits: Experimental Study and Simulation of the Behaviour of Volatile Aroma Compounds[END_REF] , bitter orange distillate [START_REF] Esteban-Decloux | Two Industrial Examples of Coupling Experimentsand Simulations for Increasing Quality and Yieldof Distilled Beverages[END_REF] , cachaça [START_REF] Scanavini | Cachaça Production in a Lab-Scale Alembic: Modeling and Computational Simulation[END_REF][START_REF] Scanavini | Cachaça Distillation Investigated on the Basis of Model Systems[END_REF][START_REF] Soares | A New Pot Still Distillation Model Approach with Parameter Estimation by Multi-Objective Optimization[END_REF] , pisco pear distillate [START_REF] Sacher | Dynamic Modeling and Simulation of an Alembic Pear Wine Distillation[END_REF] and fruit spirits [START_REF] Claus | Defining Still Parameters Using CHEMCAD Batch Distillation Model for Modeling Fruit Spirits Distillations[END_REF] . The accurate simulation of a distillation process requires physical, chemical and thermodynamic pure component properties, as well as vapor-liquid equilibrium data for the different components in the system coupled to a suitable thermodynamic model [START_REF] Valderrama | Advances on Modeling and Simulation of Alcoholic Distillation[END_REF] . While VLE data can often be obtained from extensive databases integrated in the simulation packages or from the literature, they are seldom available for multicomponent systems, particularly in dilute conditions [START_REF] Kontogeorgis | Industrial Requirements for Thermodynamic and Transport Properties[END_REF] . A fast and costefficient alternative could be to employ predictive contribution methods such as UNIFAC, although the reliability and applicability of this approach can be limited for aqueous mixtures and multicomponent systems including highly diluted compounds [START_REF] Fingerhut | Comprehensive Assessment of COSMO-SAC Models for Predictions of Fluid-Phase Equilibria[END_REF][START_REF] Gmehling | A Modified UNIFAC Model. 2. Present Parameter Matrix and Results for Different Thermodynamic Properties[END_REF] . A more reliable, albeit time-consuming alternative is to experimentally measure VLE data using static or dynamic methods.

Despite the averred importance of norisoprenoids and particularly ß-damascenone in the sector of alcoholic beverages, only a scant amount of information concerning its thermodynamic properties is available in the literature and no VLE data for the system ßdamascenone + ethanol + water have been reported to date. Furthermore, predictive methods such as UNIFAC are not readily applicable to ß-damascenone since certain elementary groups required to decompose the molecule are missing from the UNIFAC matrix. Therefore, the acquisition of experimental VLE data for ß-damascenone in hydroalcoholic media would be valuable to get a better understanding of its behavior during spirit distillation and, particularly, to be able to represent the compound in computer simulation modules.

In light of the above, this work focuses on the experimental measurement of VLE data for ßdamascenone highly diluted in ethanol + water mixtures at 101.3 kPa. Measurements were conducted in a recirculation ebulliometer for ethanol mass fractions between (0.10 and 0.79) g•g -1 and equilibrium compositions were quantified by gas chromatography and densimetry.

The experimental data were correlated by activity coefficient model NRTL and binary interaction parameters determined for the model. The reliability of the NRTL parameters to accurately simulate the behavior of ß-damascenone during a distillation was tested by simulating laboratory-scale batch distillation of a model solution of ß-damascenone, ethanol and water and comparing simulation results to experimental data.

EXPERIMENTAL PROCEDURE

Employing process simulation software to accurately predict the behavior of ß-damascenone during alcoholic distillations entails several different steps, which are explored in this section: (1) compilation of pure component properties, (2) measurement of vapor-liquid equilibrium data, (3) data correlation by thermodynamic models and (4) design of process simulation modules and (5) experimental validation.

Materials

The different components employed in the experiments are listed in Table 1 with their respective specifications and pure component properties. Deionized water (resistivity 18.2 MΩ-cm at 25 °C) obtained from a Milli-Q purification system was used throughout the experimental work. For simplification, the isomer (E)-ß-damascenone employed in the experiments will be referred to in the text as ß-damascenone. This work, measured at 101.3 kPa (u(T b ) = 0.2 K, u(P) = 0.5 kPa). c Fonseca et al. [START_REF] Fonseca | Experimental Determination of Thermodynamic Properties of Terpene and Aromatic Ketones by Gas Chromatography[END_REF] .

d Arce et al. [START_REF] Arce | Isobaric Vapor-Liquid Equilibria of 1,1-Dimethylethoxy-Butane+methanol or Ethanol+water at 101.32kPa[END_REF] . e Kamihama et al. [START_REF] Kamihama | Isobaric Vapor-Liquid Equilibria for Ethanol + Water + Ethylene Glycol and Its Constituent Three Binary Systems[END_REF] . f Riddick et al. [START_REF] Riddick | Organic Solvents: Physical Properties and Methods of Purification, 4th editio[END_REF] g PhysProp database [START_REF]PhysProp Database[END_REF] . h Hansch et al. 49 .

Pure component properties

Modeling and simulation of chemical processes in commercial software generally requires models for temperature-dependent properties such as vapor pressure, enthalpy of vaporization and heat capacity. For water and ethanol, such data were promptly available in the database of the software Simulis Thermodynamics (ProSim SA, France) [START_REF] Prosim | [END_REF] , which was employed in this work for the thermodynamic modeling and the process simulation. For ßdamascenone, which was missing from the database, pure component properties were gathered from the literature or predicted using group contribution methods, and the compound was subsequently added to the Simulis Thermodynamics database.

The saturated vapor pressure ( ) of a pure compound i at a temperature T/K was calculated from the Antoine equation: The enthalpy of vaporization ( /J•mol -1 ) for ß-damascenone at a temperature T/K was obtained from: in which empirical coefficients b 1 and b 2 were set to -87749.66 and -80.56 (Fonseca et al., 2021 [START_REF] Fonseca | Experimental Determination of Thermodynamic Properties of Terpene and Aromatic Ketones by Gas Chromatography[END_REF] ).

G. Zanghelini
Temperature-dependent heat capacity ( /J•mol -1 •K -1 ) of liquid ß-damascenone was predicted using the group contribution method proposed by Joback and Reid (1987) [START_REF] Joback | Estimation of Pure-Component Properties from Group-Contributions[END_REF] :

With = -49.26, = 1.20, = -8.12•10 -4 and = 2.42•10 -7 .

For ethanol and water, models for and values from the DIPPR database available in the ProSim suite were employed.

Vapor-liquid equilibrium of ß-damascenone in hydroalcoholic mixtures

This subsection describes the apparatus and methodology used for the VLE measurements of ß-damascenone highly diluted in ethanol-water mixtures and the thermodynamic model employed to correlate the experimental data.

G. Zanghelini et al.

Industrial & Engineering Chemistry Research (2022) 7

Experimental VLE measurements

The vapor-liquid equilibrium measurements were conducted in a glass Gillespie-type ebulliometer (Labodest VLE 602, iFisher GmbH, Germany) operating at isobaric and adiabatic conditions. System pressure is assessed by a digital manometer with an accuracy of ±0.5 kPa (P-10 WIKA, iFisher GmbH, Germany) and the temperature of the coexisting liquid and vapor phases is measured in the separation chamber by a Pt-100 platinum probe with an accuracy of ±0.05 K, periodically calibrated against a reference temperature probe (Pt-100 RTD 712, Fluke, France). A schematic diagram of the ebulliometer is provided in experiments by circulating pure ethanol.

Determination of equilibrium compositions

Ethanol concentrations in the liquid and condensed vapor samples were determined from density measurements conducted at 20 °C in an Anton Paar oscillating tube density meter (DMA 4500 M, Anton Paar) with a resolution of (±0.00001) g•cm -3 . A two-point calibration of the device was performed regularly using pure degassed water and air. The calculations of ethanol mass fractions from density at a given temperature were based on the alcoholometric tables from the International Organization of Legal Metrology (OIML) [START_REF]International Alcoholometric Tables, OIML R22[END_REF] . The accuracy of the density measurements and of the ethanol mass concentrations determined using the OIML equation were assessed beforehand by measuring the densities of ethanol-water mixtures of known compositions comprising the entire composition range and confronting the measured values with the literature [START_REF] Zanghelini | Isobaric Vapour-Liquid Equilibrium of α-Terpineol Highly Diluted in Hydroalcoholic Mixtures at 101.3 KPa: Experimental Measurements and Thermodynamic Modeling[END_REF] .

The concentrations of ß-damascenone in the vapor and liquid phases were measured in a gas chromatograph (GC Trace 1300 Series, Thermo Fisher Scientific) equipped with a flame ionization detector and an automatic sampler (TriPlus RSH, Thermo Fisher Scientific). The samples were adjusted to an ethanol mass fraction of w Et = 0.6 prior to analysis to minimize measurement variability between the different samples. All analyses were conducted in triplicate and ß-ionone (CAS number 79-77-6; Sigma) was used as the internal standard at a fixed concentration of 4.8•10 -3 g•g -1 in an ethanol-water mixture at w Et = 0.6 g•g -1 . For each analysis, an aliquot of 0.5 µL of sample was directly injected in splitless mode using a 1 µL precision glass syringe. Both inlet and detector temperatures were set to 250 °C. Helium was used as the carrier gas with a constant flow of 1. Quantification was based on calibration curves of the pure standard. Measurement uncertainties u were computed according to the GUM guidelines [START_REF] Jcgm | Evaluation of Measurement Data -Guide to the Expression of Uncertainty in Measurement[END_REF][START_REF] Taylor | Guidelines for Evaluating and Expressing the Uncertainty of NIST Measurement Results[END_REF] . The detection and quantitation limits for ß-damascenone estimated from the calibration curve were 2.6•10 -5 g•g -1 and 8.0•10 -5 g•g -1 , respectively (or 3.9•10 -6 and 1.2•10 -5 in mole fraction).

Thermodynamic modeling

The equilibrium condition entails an egality of temperatures, pressures and fugacities between the coexisting phases of saturated vapor and liquid. Given the low pressure of the VLE measurements (101.3 kPa), the vapor phase can be adequately described by the ideal gas law and the standard state liquid fugacity can be expressed as a function of pressureindependent activity coefficients [START_REF] Raal | Rigorous Characterization of Static and Dynamic Apparatus for Measuring Limiting Activity Coefficients[END_REF] . The equilibrium data can thus be represented by a modified Raoult's law: in which P/kPa is the equilibrium pressure of the system, y i and x i are, respectively, vapor and liquid equilibrium mole fractions of component i, is the activity coefficient of i in the liquid phase and /kPa is the vapor pressure of pure component i, calculated from eq 1.

In spirit distillation, the behavior of an aroma compound (AC) between a liquid and a vapor phase at equilibrium is described by two properties derived from eq. 1: its absolute volatility or partition coefficient (K AC ), which explains its distribution between the two phases, and its volatility relative to ethanol (α AC/Et ), as expressed by eqs. 5 and 6, respectively.

Pressure-independent activity coefficients were calculated in this work using NRTL (Non-Random Two-Liquid), a semi-empirical model that has been consistently employed for

G. Zanghelini et al.
Industrial & Engineering Chemistry Research (2022) 10 correlating ternary and multicomponent systems of aroma compounds highly diluted in hydroalcoholic mixtures [START_REF] Valderrama | Advances on Modeling and Simulation of Alcoholic Distillation. Part 1: Thermodynamic Modeling[END_REF][START_REF] Puentes | Review and Thermodynamic Modeling with NRTL Model of Vapor-Liquid Equilibria (VLE) of Aroma G. Zanghelini et al. Industrial & Engineering Chemistry Research[END_REF][START_REF] Faúndez | Phase Equilibrium in Binary Aqueous Mixtures of Interest in Alcoholic Distillation Using a Modified PSRK Equation of State[END_REF] . The non-random intermolecular interactions within the mixture are accounted by adjustable binary interaction parameters that are fitted to the experimental data. The model includes two sets of adjustable parameters, and , and the non-randomness parameter , which was set to 0.30 for all binaries as recommended in the literature for ethanol-water mixtures [START_REF] Renon | Local Compositions in Thermodynamic Excess Functions for Liquid Mixtures[END_REF] .

For the binary ethanol (2) + water (3), interaction parameters were set to values from Puentes et al. ( 2018) [START_REF] Puentes | Review and Thermodynamic Modeling with NRTL Model of Vapor-Liquid Equilibria (VLE) of Aroma G. Zanghelini et al. Industrial & Engineering Chemistry Research[END_REF] reported in Table 3. They have been consistently employed by our team for ternary systems involving aroma compounds in hydroalcoholic mixtures with proven reliability. For the binaries ß-damascenone (1) + ethanol (2) and ß-damascenone (1)

+ water (3), temperature-dependent parameters /J•mol -1 •K -1 and /J•mol -1 •K -1 are set to zero. Parameters and were determined from the experimental data using the software Simulis Thermodynamics [START_REF] Prosim | [END_REF] and Excel Solver by minimizing the following objective function:

in which and are, respectively, the absolute volatilities of ß-damascenone determined from the experimental data using eq. 5 and calculated with the NRTL model.

The values for were calculated using the bubble point algorithm of Simulis Thermodynamics, which calculates vapor phase composition (y) and saturated liquid temperature (T) from liquid phase composition (x) and system pressure (P).

Table 3. Binary interaction parameters for the subsystem ethanol (2) + water (3) a .

Binary (i-j) /J•mol -1 /J•mol -1 /J•mol -1 •K -1 /J•mol -1 •K -1 αij (2) - (3

Simulation of batch distillation

This segment presents a computer simulation for a laboratory scale batch distillation that was carried out to assess the capability of the adjusted NRTL binary interaction parameters determined for the system to simulate the behavior of ß-damascenone during a hydroalcoholic distillation.

Experimental setup

The experimental batch distillations were conducted for the ternary mixture in an The initial mixture loaded in the reboiler was prepared by weighing 0.001 g•g -1 of ßdamascenone, 0.250 g•g -1 of ethanol and 0.749 g•g -1 of water, for a total mass of 1508.9 g. As the mixture was brought to a boil, the column was gradually filled with ethanol-rich vapors.

Once the rising vapors reached the condenser, the column was set to operate at total reflux for about 15 minutes, then adjusted to a fixed reflux ratio of 2. The total reflux stage creates a concentration gradient for ethanol within the distillation column, which allows to achieve a better separation of ethanol and ß-damascenone based on their equilibrium behavior.

During the experimental run, distillate samples of approximately 20 g were collected over 13 time and their composition was measured as described in subsection 2.3.2. Heat losses were calculated from preliminary runs with pure water and/or pure ethanol.

Simulation procedure

The simulation module for the laboratory-scale batch distillation was implemented using the software BatchColumn version 1.4.36.0 (ProSim SA, France). The temperature-dependent vapor pressures, enthalpies of vaporization and heat capacities for ß-damascenone were computed using the models detailed in section 2.2. For ethanol and water, pure component properties were calculated using models from the DIPPR database available in the ProSim suite. The vapor-liquid equilibrium approach previously described, and the binary interaction parameters determined in this study were employed for the simulation.

RESULTS AND DISCUSSION

Vapor-liquid equilibrium

Experimental VLE data were measured at 101.3 kPa for ß-damascenone highly diluted in ethanol + water mixtures with ethanol concentrations between (0.10 and 0.79) g•g -1 . The equilibrium composition of each phase for the 16 different mixtures is reported in Table 4 in terms of mole fractions. The corresponding saturated vapor pressures, partition coefficients, relative volatilities and activity coefficients for ß-damascenone are listed in Table 5. The other standard uncertainties are u(T) = 0.2 K, u(P) = 0.5 kPa, u(ρL) = u(ρV) = 0.0002 g•cm -3 , u(Tρ) = 0.03 K, u(x Et ) = u(y Et ) = 0.004. Standard uncertainties are u(T) = 0.2 K, u(P) = 0.5 kPa, u(x Et ) = 0.004.

The binary interaction parameters determined for the two subsystems are listed in Table 6, along with the respective root mean-square deviations (RMSD) for the equilibrium variables , and .

Table 6. Binary interaction parameters regressed by NRTL model and respective root mean-square deviations (RMSD) for the binaries ß-damascenone (1) + ethanol (2) and ß-damascenone (1) + water (3). Figure 3 depicts the experimental partition coefficients and relative volatilities of ßdamascenone in terms of the ethanol mole fraction of the boiling liquid, along with NRTL correlations. Overall, the experimental data are well correlated by the NRTL model for the entire temperature range. For ethanol mole fractions at which , as denoted by the dashed line in Figure 3, the concentrations of ß-damascenone are higher in the vapor phase than in the liquid phase. Analogously, at the region where the compound is more volatile than ethanol. As observed in the figure, ß-damascenone is more concentrated in the vapor phase for x Et < 0.15, becoming more volatile than ethanol when x Et < 0.06. This is reflected in the absolute and relative volatilities of ß-damascenone, which increase substantially with decreasing ethanol mole fractions, reaching maximum predicted values of 65.2 and 5.43 for and , respectively. This high volatility in the water-rich region is accompanied by high activity coefficients (Table 5), which indicate that deviations from ideality in this region are quite pronounced. As ethanol concentration in the liquid phase further increases, values of decrease over 400 times, reaching 0.10 when x Et = 0.59. The behavior observed for ß-damascenone can be explained by two main factors: [1] the higher solubility of most aroma compounds in ethanol than in water, which reduces their concentration in the headspace above solutions that are richer in ethanol; [2] the relatively high hydrophobicity of ß-damascenone, conveyed by a log K ow of 4.20. This high hydrophobicity results in the aroma compound preferably positioning itself with ethanol in the interfacial layer between the aqueous liquid and the headspace, being more prone to evaporation, especially when lower amounts of ethanol are present at the interface [START_REF] Ickes | Effects of Ethanol on Flavor Perception in Alcoholic Beverages[END_REF][START_REF] Tsachaki | Dynamics of Flavour Release from Ethanolic Solutions[END_REF] .

Binary (i-j) Correlation parameters a RMSD b /J•mol -1 /J•mol -1 (1)
Furthermore, the distribution of ß-damascenone between the coexisting vapor and liquid phases at different ethanol compositions is modulated by physicochemical effects involving structural changes associated with the ethanol-water matrix. When small amounts of ethanol are present in aqueous solutions (x Et < 0.05), the ethanol molecules tend to be monodispersed within the hydrogen bonded networks formed by water molecules [START_REF] Conner | Headspace Concentrations of Ethyl Esters at Different Alcoholic Strengths[END_REF][START_REF] D'angelo | Self-Association Behaviour of Alcohols in Diluted Aqueous Solutions (*)[END_REF] . As the solution becomes richer in ethanol, ethanol molecules begin to aggregate to form micelle-like structures that could entrap hydrophobic aroma molecules, thus limiting their release into the headspace [START_REF] Onori | Dynamical and Structural Properties of Water/Alcohol Mixtures[END_REF] . This could explain the trend observed in Figure 3.

Batch distillation and process simulation

The results for the experimental batch distillation are depicted in Figure 4 as hollow shapes. The solid lines in Figure 4 depict the simulated mass concentration profiles of ethanol and ßdamascenone in the distillate (w Et and w ßd , respectively) and the cumulated masses of distillate and of ethanol. In Figure 4(a), a dashed blue line was added to represent ethanol mass fraction in the condenser during the filling and total reflux stages. Given the gradual increase of ethanol mass concentration within the column in the filling stage, the first vapors reach the condenser with a high ethanol concentration, which remains nearly constant during the total reflux stage. As observed in the figure, the experimental data are overall well represented by the simulation module. Minor discrepancies can be observed in Figure 4(a)

near the concentration peak of ß-damascenone, where concentrations for both ßdamascenone and ethanol seem to be overestimated by the simulation. Since ethanol concentrations in the distillate undergo a steep decrease in the region, the peak of ßdamascenone rises rather swiftly and only a limited amount data points were collected in the zone. Furthermore, the experimental concentrations of ß-damascenone are higher than the simulated values for the last 50 min of distillation. It is possible that small amounts of the compound were retained on the walls between the reflux valve and the cooler throughout distillation. If this was case, as the distillate became increasingly aqueous and the liquid flowed at a slower rate, the ß-damascenone retained could have been carried by the distillate into the samples. Nevertheless, the concentration profile of ß-damascenone is sufficiently well represented to consider that the NRTL parameters determined in this work would allow to correctly simulate its behavior during a spirit distillation.

CONCLUSIONS

New VLE data were measured for ß-damascenone highly diluted in ethanol-water mixtures using a recirculating ebulliometer at 101.3 kPa. Pure component properties for ßdamascenone were collected from the literature and computed with predictive methods, and the compound was added to the Simulis Thermodynamics database. The new experimental data were well correlated by the NRTL model and binary interaction parameters were determined for the binaries ß-damascenone-ethanol and ß-damascenone-water. The applicability of the experimental data for simulation purposes was checked by simulating a laboratory-scale batch distillation of a model solution for the ternary system. The binary NRTL parameters allowed to correctly represent the behavior of ß-damascenone during a batch distillation and could be used to predict the behavior of the compound during an industrial-scale distillation of alcoholic beverages. This approach can be extended to other aroma compounds that play a key role in the typicity of distilled spirits.

Figure 1 .

 1 Figure 1. A validation procedure of our apparatus was conducted beforehand by measuring the binary ethanol-water at 101.3 kPa and assessing the thermodynamic consistency of the data, as described by Zanghelini et al. (2022) 54 .

Figure 1 .

 1 Figure 1. Schematic diagram of the Labodest VLE 602 (iFischer) apparatus: (1) mixing chamber, (2) reboiler, (3) "Cottrell pump" or equilibrium chamber, (4) separation chamber, (5) Pt-100 temperature probe, sampling tubes for the (6) vapor and (7) liquid phases. The arrows indicate the sense of circulation: (straight teal arrow), vapor-liquid; (dotted orange arrow), vapor; (dashed blue arrow), liquid.

  2 mL•min -1 . Hydrogen and air flow for the flame jet were respectively 35 and 350 mL•min -1 . The stationary phase was a ZB-Wax Plus capillary column (30 m x 0.32 mm, 0.5 µm). The initial oven temperature was initially set to G. Zanghelini et al. Industrial & Engineering Chemistry Research (2022) 9 held for 5 min. The total running time was 15 min. The chromatographic data were acquired and analyzed using Chromeleon CDS software version 7.2.10 (Thermo Scientific).

Oldershaw-type glass

  distillation column (Vasse) composed of two thermally insulated column sections with five perforated plates each. The columns are connected at the bottom to a reboiler consisting of a triple neck round bottom flask with a volume capacity of 2 L and a 500 W heating mantle with an experimentally determined power of 375 W measured for pure water. A schematic representation of the apparatus is displayed in Figure2. The temperature of the boiling liquid throughout the distillation process was measured by a Pt 100 temperature probe inserted in the flask and pressure drop was assessed by a watercolumn manometer connected to the reboiler. Reflux and sampling were controlled by a magnetic solenoid valve whose opening (sampling) and closing (reflux) were regulated by a control unit.

Figure 2 .

 2 Figure 2. (a) experimental setup and (b) scheme of the laboratory-scale batch distillation column. T B , temperature of the boiling liquid; ΔP, pressure drop; L 1 , reflux; L 11 , liquid to the reboiler; F, feed; D, distillate; V, rising vapor; (red dashed arrow), vapor circulation; (straight blue arrow), liquid circulation.

Figure 3 .

 3 Figure 3. Evolution of (a) the absolute volatility ( ) and (b) the relative volatility (

Figure 4 (

 4 Figure 4(a) shows the mass concentration profiles of ethanol and ß-damascenone throughout

Figure 4 .

 4 Figure 4. Experimental and simulated data for the batch distillation of a model solution consisting of ß-damascenone, ethanol and water. (a) mass concentration profiles of () ethanol and () ß-damascenone in the distillate, (b) cumulated masses of () distillate and of () ethanol as a function of time and (c) () evolution of the temperature in the reboiler over time. ---(solid lines), data simulated with BatchColumn (ProSim SA, France) using the binary NRTL coefficients determined in this work.

Figure 5 22 Figure 5 .

 5225 Figure 5 presents the evolution of the liquid and vapor compositions and the temperature throughout the column at different distillation times of the computer simulation: (a) 24'50'', (b) 40'00'', (c) 90'00'' and (d) 110'00''. Stages 1 and 12 correspond to the condenser and the reboiler, respectively, and stages 2 to 11 correspond to the 10 column plates. Fig. 5(a)represents the beginning of the total reflux period, when a gradient of decreasing ethanol concentration and increasing temperature can be observed from the condenser and upper stages to the reboiler; at this point, ß-damascenone has begun to rise into the bottom of the column, propelled by the low ethanol concentrations in the reboiler and in the lower stages.By the end of the period in which the column operates at total reflux (Fig.5(b)), the column plates are highly concentrated in ethanol, prompting most of the ß-damascenone to be refluxed into the reboiler. A steep increase in concentration for ß-damascenone can be observed again at 90 min (Fig.5(c)), when the ethanol and temperature gradients within the column are more pronounced and ethanol concentration in the reboiler is lower than 0.1 g•g - 1 . After another 20 min of distillation (Fig.5(d)), the column was nearly depleted of ethanol and ß-damascenone, and the temperature of the liquid, which is essentially governed by ethanol composition, is significantly higher.

Table 1 .

 1 Specifications and physical properties of the aroma compounds used in this work: molecular weight (MM/g•mol - 1 ), normal boiling point (T b /K) at 101.3 kPa, log 10 of the octanol/water partition coefficient at T = 298.15 K (log K ow ).

	Compound	(E)-ß-damascenone	ethanol	water
	CAS number	23696-85-7	64-17-5	7732-18-5
	Formula		C 13 H 18 O	C 2 H 6 O	H 2 O
	Structure
	Supplier		Sigma-Aldrich	Carlo Erba
	Mass fraction purity	≥0.992 a	≥0.999 a
	MM/g•mol -1	190.28	46.07	18.01
	T b /K	Exp.		351.45 b	373.15 b
		Lit.	542.22 c	351.44 d,e,f	373.15 d,e,f
	log K ow		4.20 g	-0.31 h	-1.38 h

a determined by gas chromatography by the supplier. b

  et al.

	Industrial & Engineering Chemistry Research (2022)
	in which A i , B i and C i , are empirical coefficients specific to each compound.
	The Antoine coefficients employed for the components in this work are listed in Table 2. For
	ß-damascenone, the coefficients were determined by regression of experimental data
	measured by Nelson (2018) 51 .
	6

Table 2 .

 2 Antoine coefficients for calculating the vapor pressures of each component in the mixture at a given temperature using eq 1. T min /K and T max /K define the temperature interval in which the equation coefficients are valid.

	Compound	A i	B i	C i	T min /K	T max /K	Source
	ß-damascenone	21.13	4261.91	-85.1995	298.15	450	51
	ethanol	23.8046	3803.98	-41.68	270	369	52
	water	23.1963	3816.44	-46.13	284	441	52

Table 4 .

 4 Experimental VLE data for ß-damascenone (ßd) highly diluted in ethanol (Et) + water (w) mixtures. With T/K, equilibrium temperature; ρ/g•cm -3 , density of the mixture in the liquid (L) and condensed vapor (V) phases at 293.15 K; x and y, mole compositions of the liquid and vapor phase, respectively, and u(x ßd ) and u(y ßd ) the standard uncertainties for ßdamascenone.

	T/K	ρ L /g•cm -3	x Et	10 5 xßd	10 5 u(xßd) ρ V /g•cm -3	yEt	10 5 yßd	10 5 u(yßd)
	365.4	0.9861	0.029	7.132	0.027	0.9157	0.276 107.954	0.218
	365.3	0.9847	0.033	7.619	0.027	0.9144	0.282 132.825	0.262
	363.1	0.9808	0.045	8.330	0.027	0.8967	0.349	86.003	0.168
	363.1	0.9792	0.050	8.946	0.028	0.8985	0.342	86.390	0.160
	361.3	0.9759	0.062	9.873	0.028	0.8860	0.393	60.202	0.104
	361.3	0.9751	0.065	9.686	0.027	0.8847	0.399	57.426	0.096
	360.6	0.9734	0.071	7.976	0.023	0.8802	0.419	35.547	0.074
	360.0	0.9708	0.081	8.967	0.025	0.8765	0.435	31.954	0.065
	358.0	0.9600	0.121	11.986	0.028	0.8666	0.483	16.421	0.024
	357.9	0.9587	0.126	11.936	0.027	0.8643	0.494	14.252	0.021
	356.7	0.9480	0.164	10.528	0.026	0.8597	0.518	8.284	0.014
	356.6	0.9456	0.172	10.851	0.026	0.8579	0.527	7.428	0.013
	355.7	0.9287	0.229	14.563	0.028	0.8530	0.555	9.783	0.016
	355.0	0.9171	0.269	11.182	0.024	0.8498	0.573	5.585	0.012
	355.0	0.9149	0.277	11.234	0.024	0.8496	0.574	5.302	0.012
	352.5	0.8468	0.591	18.453	0.028	0.8295	0.698	1.905	0.017

Table 5 .

 5 Experimental VLE variables describing the behavior of ß-damascenone (ßd) highly diluted in hydroalcoholic mixtures at P = 101.3 kPa. With T/K, equilibrium temperature; x Et , ethanol mole fractions in the liquid phase; , saturated vapor pressure of ß-damascenone; Kßd and KEt, partition coefficient of ß-damascenone and ethanol, respectively; αßd /Et , relative volatility of ß-damascenone in relation to ethanol; , activity coefficient of ß-damascenone.

	T/K	x Et	/kPa	Kßd	KEt	αßd /Et
	365.4	0.029	0.37	15.136	9.55	1.59	3.61
	365.3	0.033	0.37	17.433	8.55	2.04	3.68
	363.1	0.045	0.33	10.325	7.77	1.33	3.50
	363.1	0.050	0.33	9.657	6.78	1.42	3.47
	361.3	0.062	0.30	6.098	6.36	0.96	3.31
	361.3	0.065	0.30	5.929	6.16	0.96	3.30
	360.6	0.071	0.29	4.457	5.90	0.75	3.20
	360.0	0.081	0.28	3.564	5.38	0.66	3.11
	358.0	0.121	0.25	1.370	3.98	0.34	2.75
	357.9	0.126	0.25	1.194	3.92	0.30	2.69
	356.7	0.164	0.23	0.787	3.16	0.25	2.54
	356.6	0.172	0.23	0.685	3.07	0.22	2.48
	355.7	0.229	0.22	0.672	2.42	0.28	2.49
	355.0	0.269	0.21	0.499	2.13	0.23	2.38
	355.0	0.277	0.21	0.472	2.07	0.23	2.36
	352.5	0.591	0.18	0.103	1.18	0.09	1.76

Table 7 .

 7 The distillation stage was theoretical distillation stages in the simulation module to integrate the variations in heat losses over time and thus achieve a better adjustment with the experimental data. The condenser was regarded as a total condenser and the liquid in the reboiler was considered to remain a homogeneous mixture throughout distillation.

	G. Zanghelini et al.	Industrial & Engineering Chemistry Research (2022)

Table 7 .

 7 Input parameters of the simulation module.

	Simulation specifications
	Number of trays	12
	Efficiency
	Trays 2-11	0.75
	Reboiler	1.00
	Condenser	1.00
	Liquid holdup
	Trays 2-11	3 mL
	Condenser	5 mL
	Pressure a
	Condenser	103.0 kPa
	Pressure loss	0.8 kPa

a as measured during the experimental distillation.

Acknowledgements

The authors are thankful to Nicolas Descharles for his valuable assistance with the analytical measurements and to N'Guessan Charles Romaric Tano for his aid with distillation experiments.

Funding

This work was supported by the foundation Jean Poupelain.

AUTHOR INFORMATION

Corresponding author

*Gabriela Zanghelini. E-mail: gabriela.zanghelini@agroparistech.fr ORCID Gabriela Zanghelini: 0000-0001-6993-3878

Violaine Athès: 0000-0002-2194-8517 Martine Esteban-Decloux: 0000-0002-8084-2528 Pierre Giampaoli: 0000-0002-6756-0403 Stéphane Vitu: 0000-0002-5442-1279

Notes

The authors declare no competing financial interest.


Mariana Padron
email: mpadron@obs-banyuls.fr.

Guizien

Katell

Modelling the

Keywords: connectivity, metapopulation, seascape genetics, spatially explicit model, sedentary benthic invertebrates, Gulf of Lions

Introduction

Connectivity strongly impacts the dynamics and persistence of marine populations because it influences key processes, from short-term demography to long-term evolution, including resistance to disturbances and environmental changes. Therefore, understanding the development and maintenance of connectivity patterns is essential in marine spatial planning and the proper design and management of Marine Protected Areas (MPAs) [START_REF] Palumbi | Marine reserves and ocean neighborhoods: The spatial scale of marine populations and their management[END_REF][START_REF] Botsford | Connectivity and resilience of coral reef metapopulations in marine protected areas: Matching empirical efforts to predictive needs[END_REF].

However, the quantitative study of connectivity and the identification of connectivity bottlenecks that could compromise the ability for populations to recover after disturbances are current scientific challenges [START_REF] Pineda | Linking larval settlement to larval transport: assumptions, potentials and pitfalls[END_REF][START_REF] Burgess | Beyond connectivity: How empirical methods can quantify population persistence to improve marine protected-area design[END_REF].

For sedentary marine benthic species, understanding connectivity translates into understanding larval dispersal, given that it is during the pelagic larval stage that the exchange of individuals and genes among populations takes place [START_REF] Pineda | Larval transport and dispersal in the coastal ocean and consequences for population connectivity[END_REF][START_REF] Cowen | Larval dispersal and marine population connectivity[END_REF]. Hence, to analyse the patterns of connectivity among populations and study, in more detail, the mechanisms that shape and maintain connectivity patterns through time, the dispersal of organisms should be taken into account [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF][START_REF] Aiken | Environmental fluctuations and asymmetrical dispersal: Generalized stability theory for studying metapopulation persistence and marine protected areas[END_REF].

Different approaches have been undertaken to try to quantify dispersal and connectivity in marine ecosystems. Population genetics has been one of the most frequently used tools to provide direct and indirect measures of connectivity in several marine species [START_REF] Hellberg | Genetic assessment of connectivity among marine populations[END_REF][START_REF] Saenz-Agudelo | Estimating connectivity in marine populations: an empirical evaluation of assignment tests and parentage analysis under different gene flow scenarios[END_REF][START_REF] Lowe | What can genetics tell us about population connectivity?[END_REF]. The dispersal potential of the larvae produced by sedentary benthic species has usually correlated with gene flow and was even used as a predictor of population genetic structure [START_REF] Siegel | Lagrangian descriptions of marine larval dispersion[END_REF].

Using estimates of Pelagic Larval Duration (PLD) as a proxy to measure realised dispersal distance, several studies showed that the correlation between PLD and genetic Isolation by Distance (IBD) was consistent with the expectation that migration was the key factor driving the genetic structuration of marine populations [START_REF] Hellberg | Dependence of gene flow on geographic distance in two solitary corals with different larval dispersal capabilities[END_REF][START_REF] Gilg | The geography of marine larval dispersal-coupling genetics with fine-scale physical oceanography[END_REF]. However, recent data set reviews [START_REF] Riginos | Effects of geography and life history traits on genetic differentiation in benthic marine fishes[END_REF][START_REF] Selkoe | Marine connectivity: A new look at pelagic larval duration and genetic metrics of dispersal[END_REF][START_REF] Faurby | Theoretical limits to the correlation between pelagic larval duration and population genetic structure[END_REF] have now shown that the relationship between PLD and genetic structuration might be weaker than previously thought.

The influence of several processes (e.g. larval behaviour, availability of suitable habitat, and biotic interactions) can cause the patterns of population connectivity to be different from what would be expected from flow integration along PLD alone [START_REF] Pineda | Linking larval settlement to larval transport: assumptions, potentials and pitfalls[END_REF][START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: Turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF]. Among these processes, settlement or recruitment regulation of dispersal alters gene flow [START_REF] Pineda | Larval transport and dispersal in the coastal ocean and consequences for population connectivity[END_REF], in addition to the local demography. This temporal and spatial variability in the strength of gene and demographic flow due to bottlenecks during larval settlement and recruitment implies that there is a set of interacting processes operating across various scales, which introduces complexity into the system. Therefore, inferring connectivity by only one disciplinary approach (either genetic or bio-physic) does not provide the necessary information to properly describe population persistence in a network of MPAs [START_REF] Hastings | Metapopulation dynamics and genetics[END_REF][START_REF] Marko | The complex analytical landscape of gene flow inference[END_REF].

It is therefore necessary, in order to efficiently inform conservation management decisions, to use tools that integrate bio-physical models with information regarding local demographic traits driving genetic linkages that result from larval exchange among populations [START_REF] Werner | Coupled biological and physical models[END_REF][START_REF] Cowen | Larval dispersal and marine population connectivity[END_REF]. However, detecting the genetic signatures caused by the transient dynamics of demographic changes remains a challenge [START_REF] Alcala | Peak and persistent excess of genetic diversity following an abrupt migration increase[END_REF].

Although genetic models incorporating larval dispersal and demographic parameters have recently been developed to quantify connectivity among marine populations [START_REF] Galindo | Seascape genetics: A coupled oceanographicgenetic model predicts population structure of Caribbean corals[END_REF][START_REF] Kool | An object-oriented, individual-based approach for simulating the dynamics of genes in subdivided populations[END_REF][START_REF] Munroe | The role of larval dispersal in metapopulation gene flow: Local population dynamics matter[END_REF], they have some limitations when attempting to understand connectivity patterns among non-equilibrium populations. The model [START_REF] Galindo | Seascape genetics: A coupled oceanographicgenetic model predicts population structure of Caribbean corals[END_REF] developed, for example, does not take into account demographic variability, as population sizes were kept constant over time. The approach [START_REF] Kool | An object-oriented, individual-based approach for simulating the dynamics of genes in subdivided populations[END_REF] developed, although allowing changes in population size, does not account for overlapping generations but, more importantly, being an Individual-Based Model (IBM), population sizes were limited to 10 individuals.

Recently, [START_REF] Munroe | The role of larval dispersal in metapopulation gene flow: Local population dynamics matter[END_REF], presented a modelling case study for Eastern oysters in Delaware Bay (USA) in which they relaxed population size limitations and assigned each population its own demographic dynamics.

In an attempt to generalise Munroe et al.'s (2012) approach, we present a spatially explicit coupled demographic gene flow model that assesses the effect of demography variability on allele frequencies in a marine metapopulation of sedentary benthic species with realistic population densities. The seascape model was used to examine how: 1) demographic traits influence the development of genetic structure in closed populations over time, 2) different patterns of connectivity affect the loss or maintenance of genetic diversity in a metapopulation, 3) fluctuations in local population density help maintain genetic diversity over time and 4) local population stability impacts genetic diversity to provide relevant information for conservation decisions in the Gulf of Lions, France.

Methods

A seascape genetic model for benthic invertebrates with a pelagic larval phase

To simulate the transient demographic effect on local allele frequencies arising from environmental variability modulating connectivity through time, a seascape genetic model was developed. The seascape model relies on the forcing of a spatially explicit dynamical model of allele probability by a spatially explicit metapopulation model (Guizien et al., 2014, adapted from Hastings and[START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF]) (Fig. 1). The combination of these models allows us to project the development of a population genetic structure through time, resulting from demography dynamics driven by local survival and larval transfer on a regional scale.

The spatially explicit metapopulation model simulates the population density Ni(t) at each site i at any time t resulting from the balance of adult survival and input of new recruits resulting from reproduction of adults present in the previous time step. Therefore, Ni represents the adult census size per unit surface, and it is equal to the effective number of breeders per unit surface for the next reproductive event (as defined by [START_REF] Waples | Simple life-history traits explain key effective population size ratios across diverse taxa[END_REF]. It is important to note that such formulation implicitly hypothesized that age at first reproduction and the time lapse between reproductive events are both equal to the model time step, ensuring that effective number of breeders does not vary with population age structure. Larval transfer rate among sites is estimated from a set of realistic connectivity matrices obtained using the larval dispersal model [START_REF] Guizien | Using larval dispersal simulations for marine protected area design: Application to the Gulf of Lions (northwest Mediterranean)[END_REF], and recruitment success (inter-specific competition or predation) is limited by saturating population density (intra-specific competition) (Eq. 1 and 2) [START_REF] Guizien | Vulnerability of marine benthic metapopulations: Implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. with the transfer matrix G(t) defined by:

N i (t+1)=min ( G(t ) N i (t ) , N max ) (1) N i (t+1) =
G ji (t)=a i C ji (t )b j + s jj δ ji
(2) a i = Propagule production rate at site i [larvae per adult] C ji (t )= Propagule transfer rate from site i to site j multiplied by the ratio of surface of site i to surface of site j [no units] b j = Recruitment success at site j [adults per larva] s jj =Adult survivorship rate at site j [no units] δ ji = Kronecker symbol (1 when i=j, 0 otherwise) Maximum population growth rate (Pi) of a closed population in site i is then defined by:

P i =a i b i + s ii (3)
The time it takes for a population at any initial density N o to reach its saturating density N max , from any given initial density was estimated as:

T = ln ( N max N o ) ln ( P i ) (4)
In the seascape genetics model, the frequency of each allele k is calculated at any given time in every site i, as a result of mixing between the surviving individuals with allele k in site i and the input of new recruits with allele k, coming from any other site (Eq. 5). Given the large number of larvae produced by each adult and the high larval survival values assigned (10%) compared to estimated regional retention rates (3%), genetic drift due to larval mortality is neglected and, thus, larval mortality is simulated as a proportional and uniform loss that is applied to the larval pool with allele k produced at each site. The effect of genetic drift is therefore simulated by the random sampling of the alleles lost during two different periods in the life cycle: larval dispersal and recruitment. Simplifying assumptions were that: 1) mating takes place randomly every year among diploid organisms, 2) allele frequencies are identical among sexes, 3) generations are discrete but overlapping, 4) populations are not age-structured, 5) given the relatively short period of time considered in the simulations (100 years), mutation is assumed to be negligible compared to genetic drift and migration, and thus set to zero and, 6) selection can be neglected.

During larval dispersal, the pool of larvae with allele k produced by each site, remaining in the region (accounting for regional retention rate), is randomly sampled before being re-distributed among the destination sites according to the connectivity matrix (Eq. 6 and 7). The number of larvae with allele k that will survive dispersal is computed by sorting out the rank of the allele carried by the surviving larvae. During recruitment, saturating density limits the number of arriving larvae that can be recruited at each site, therefore, a random subset of alleles carried by the pool of larvae arriving at each site is actually recruited.

p ik (t +1) = ( s ii N i( t ) p ik (t ) ) +b j ∑ j C ij ∑ l C lj S l S j n ' jk i N i( t +1) (5)
where n ' jk is a random number between 0 and n jk , with the constraint that ∑

k n' jk i =n ' j i n ' j i = A j a j ( ∑ l C lj S l S j ) (6) n jk =p jk A j a j (7) p ik (t +1) = Frequency of allele k at site i at time t+1 [no units] p ik (t ) = Frequency of allele k at site i at time t [no units] p jk (t)= Frequency of allele k at site j at time t [no units] S i = Surface area of site i [m 2 ] A j =

Test cases

The simulations were designed to evaluate the changes in allele frequencies due to variations in species demographic parameters and connectivity structure. Although this modelling approach applies to any sedentary species with a reproductive dispersive stage, to test the sensitivity of the model, we focused on marine benthic species with four different life expectancies

(2, 5, 10 and 20 years) that reproduce each year. The values of minimum recruitment success were varied in order to test the effect of maximum population growth in the maintenance of genetic diversity over time. Minimum recruitment success was also the same in all populations (10%), although it was limited by saturating spatial density, which was set to Nmax=6250 individuals/m 2 per population.

Six groups of simulations were performed (Table 1). All simulations were run for 100 iterations, one iteration being the time lapse between two reproductive events. The initial genetic structure was described by ten different alleles present simultaneously in every site. In order to remove the potential bias due a specific initial condition for the genetic structure, simulations were repeated at least 100 times with different and random initial conditions of allele frequencies.

The first group consisted of simulations of a closed population (no larval inputs) in which life expectancy was raised from two to 20 years. Groups 2, 3A, 3B, 4A and 4B, consisted of simulations of a metapopulation of a species with a life expectancy of 10 years, distributed among 32 sites along the Gulf of Lions (northwestern Mediterranean) (Fig. 2), in which spatial structure was driven by connectivity only. These groups of simulations correspond to the test case of dominating species of the soft-bottom benthic communities in the Gulf of Lions (polychaetes). The continuous sandy-bottom habitat was described by sites of approximately the same width (~ 20 km) contiguously distributed along the coast spanning the 10 to 30 meters isobaths. Based on these criteria, the areas of the 32 sites ranged from 6.2 to 121 km 2 . Survival and fecundity rates were the same among populations within the metapopulation and were kept constant over time in all of the simulations.

Larval transfer between the 32 sites was computed from Lagrangian larval dispersal simulations based on regional hydrodynamic simulations performed with the coastal circulation model SYMPHONIE [START_REF] Marsaleix | Energy conservation issues in sigma-coordinate free-surface ocean models[END_REF]. Larval dispersal was simulated by releasing larvae at the average site depth (20 m). Connectivity matrices containing larval transfer probabilities were produced for 20 sets of 10-day long spawning periods in 2004 and 2006 with an average PLD of four weeks [START_REF] Guizien | Vulnerability of marine benthic metapopulations: Implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. A total of 20 different connectivity matrices were used for this analysis.

Simulations of Group 2 were performed to evaluate the effect of the distinct connectivity structure on heterozygosity. In this case, connectivity was set deterministically using the same connectivity matrix every year (either one of the 20 variants) and, thus, simulations were only run 100 times.

Groups 3A, 3B, 4A and 4B were performed in order to test the effect of recruitment (3A and 3B) and habitat loss (4A and 4B) on metapopulation size and the maintenance of heterozygosity within a metapopulation over time. Metapopulation size was represented as average regional coverage, defined as the ratio between regional average population density and saturating density. Connectivity was set stochastically using a different connectivity matrix among the 20 variants every year and, thus, simulations were run 500 times to reach the convergence of the mean when combining stochasticity in connectivity and initial condition. Group 3A simulated gene flow among all populations, and recruitment success was constant every year. Group 3B explored a scenario where recruitment failure (no recruitment) occurred every three years. Group 4A simulated the loss of four populations (sites 17, 18, 19 and 20). These populations were identified as being essential for the regional persistence of species at minimum recruitment success in metapopulation modelling, comparing the vulnerability to habitat loss of four ports in the Gulf of Lions [START_REF] Guizien | Vulnerability of marine benthic metapopulations: Implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. Nonetheless, it is important to mention that for the present simulations, the values of minimum recruitment success and fecundity were set to ensure longterm persistence of the metapopulation after removing those four sites, unlike in [START_REF] Guizien | Vulnerability of marine benthic metapopulations: Implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. Group 4B evaluated a scenario where only one population (site 21) was lost from the metapopulation. This site was characterised by consistently showing the highest values of betweenness in the metapopulation when applying graph theory analysis to connectivity matrices (Costa, pers. comm.).

Metrics and statistics

Heterozygosity and allelic richness are two measures of genetic diversity commonly used in the population genetics literature, and in the present study both were used. Heterozygosity (Hei), was estimated taking into account all individuals at each site i over time (Eq. 7).

❑ ❑ () ∑ ❑ ❑ (❑ ❑ ( )) ❑ (7)
where p ik (t ) = frequency of the kth allele of n alleles in site i at time t [no units].

However, heterozygosity and allelic richness describe the state of genetic diversity of a population at a particular time and not its dynamics. In order to measure the speed at which genetic diversity varied as a function of demographic parameters (regulating genetic drift) and connectivity structure, new metrics describing the initial and average slope of the evolution of genetic diversity over time are introduced: 1) Allele Drifting Time (Td), which is the moment at which allelic diversity starts decreasing, meaning the time when the first allele is lost from the population or the metapopulation. 2) Allele Fixation Time (Tf), which is the time lapse between the first loss of an allele (Td) and the fixation of any allele (Fig. 3).

These metrics were defined on allelic richness rather than heterozygosity to avoid ambiguity in defining genetic diversity reference points and given the known sensitivity of allelic richness to the loss of rare alleles (Allendorf 1986).

❑ ❑ () ∑ ❑ ❑ (❑ ❑ ()) ❑ p ik (t )

Results

Do differences in demographic parameters affect genetic diversity in closed populations?

Figures 4A and 4B display the loss of genetic diversity in closed populations at initial saturating capacity as a function of the ratio between mortality (being complementary to unity of survival) and recruitment success. This ratio, having the dimensions of larvae/adult/year, depicts the origin of recruitment regulation (sensu by [START_REF] Hixon | Population regulation: Historical context and contemporary challenges of open vs. closed systems[END_REF] and informs on the relative contribution of intrinsic vs. extrinsic drivers: high values indicate recruitment regulation due to low recruitment success driven by inter-specific competition or predation, while low ratio values indicate recruitment regulation due to space limitation at a high population density (intra-specific competition). The observed relationship between the loss of genetic diversity (measured as Td) and recruitment success in these simulations, as it is defined in Equation 2, also reflects the relationship between Td and local retention in the case of an isolated population.

Fig. 4A shows that the number of reproductive events after which the population begins to loose its genetic diversity (Td/Re) increases with the ratio between mortality and recruitment success. This demonstrates a reduction of genetic drift when recruitment success decreases for any given life expectancy. However, genetic drift is inhibited at high levels of recruitment success, as depicted by low and almost constant Td/Re for the lowest values of the ratio between mortality and recruitment success. Such minimum Td/Re is of the order of magnitude of the number of reproductive events in a life expectancy. For a recruitment frequency of 0.18, short-lived species (two years) could start loosing their genetic diversity after three years, while it would take around 14 years for a long-lived species (20 years). This reflects the effect of genetic drift inhibition by the reduction of reproductive events due to the lack of space caused by the presence of adults as life expectancy increases.

Similarly, the maintenance of genetic diversity, measured by the number of reproductive events before allele fixation (Tf/Re), in closed populations at an initial saturating density increases with life expectancy, regardless of whether the population growth rate is low (<2.5) or high (>2.5) (Fig. 5). This confirms that genetic drift is limited by population inertia due to recruitment inhibition as a result of space limitation, as neither fecundity nor recruitment success affects the relationship.

However, when dealing with overlapping generations, it must be taken into account that the reduction of genetic drift by space limitation during the first population turnover decreases with species longevity. This is exhibited by the different hierarchies between life expectancy curves when representing the time at which the population starts to loose its genetic diversity in terms of population turnover (adimensionalising Td by the life expectancy of the species) (Fig. 4B). Species with short life expectancies begin to experience a reduction of their genetic diversity after a higher number of population turnovers than do long-lived species, as fewer reproductive events occur during a short life expectancy compared to during a long one. This result confirms that genetic drift can be controlled by effective reproduction (which may be limited by long life expectancy or low recruitment success) and its frequency along the species lifespan in a context of finite population size.

Different initial population densities also affect the time at which a population starts loosing genetic diversity (Fig. 6A) -for any value of the ratio between mortality and recruitment success, the lower the initial density of the population, the slower the loss of genetic diversity. This result points out that genetic drift is slower when considerable space is available for recruitment compared to at saturating density where space availability is regulated by mortality. Quantifying space availability by the time it will take a population to reach saturating density (6250 individuals/m 2 ) shows that it would take twice as long for a population with an initial density of 50 individuals/m 2 compared to a population with an initial density of 1000 individuals/m 2 (Fig. 6B).

Allele drifting time increases linearly with the time necessary to reach saturating density.

The population growth rate triggers the effect of initial population density when it comes to the maintenance of genetic diversity, measured as Tf (Fig. 5). While in fast growing populations, the maintenance of genetic diversity is not sensitive to initial density, in slow growing populations, the maintenance of genetic diversity lasts longer and is more variable when the initial population size is small (50 individuals/m 2 ) compared to when it is at saturating density. Such sensitivity to initial population size is amplified at small life expectancies.

Overall, in closed populations, genetic drift is the quickest when the population is at saturating density with high recruitment success, slowing down at low recruitment success or when population size is lower than saturating density. This result demonstrates the importance of transient demography in regulating the speed of genetic drift.

Do different patterns of connectivity affect the maintenance of genetic diversity in a metapopulation?

Simulated heterozygosity patterns among populations forming a metapopulation varied according to each connectivity matrix. The two most contrasting spatial patterns with low heterozygosity values either in the west (A) or in the east (B) of the Gulf of Lions are shown in Fig. 7. A plausible explanation is the different spatial patterns of local retention (defined as the proportion of larvae produced locally that remain in the same spatial unit) vs. import (estimated as the proportion of all larvae produced among all sites that settled within the focal site). For connectivity matrix 1, the ratio between local retention and import was higher in the west than in the east (0.074 and 0.027, respectively). On the contrary, for connectivity matrix 2 the values of the same ratio were higher in the east compared to the west (0.069 and 0.0085, respectively).

However, simulations demonstrated that the relationship between Td and import do not exhibit any trend with no linear correlation (R 2 = 0.011; P < 0.01) (Figure S1). Similarly, Td was not linearly correlated to local retention rate values obtained from the analysis of the 20 connectivity matrices (R 2 = 0.031; P > 0.05). However, Td was delimited by a piecewise function enabling the populations (sites 17-20) responsible for bringing demographic persistence to the metapopulation at minimum recruitment success divides it into an eastern and western group, causing a slightly faster decrease of genetic diversity in the latter group than in the former. However, Scenario D

shows how the loss of only one site, identified as essential for maintaining the integrity of the connectivity graph, indeed disrupts the genetic connectivity of the metapopulation by augmenting the differentiation between the eastern and western groups in terms of the loss of heterozygosity, with drifting starting earlier in the eastern part. In Scenarios C and B, recruitment success was set in order to ensure that the loss of sites within the metapopulation would still allow for species persistence. It should be noticed that in both scenarios, a similar constant metapopulation size was maintained, even if not at saturating density. The effect of migration in this case will also lead to the loss of genetic diversity in the metapopulation over time, although at a slower rate.

Discussion

The present study highlights that genetic drift increases considerably with settlement and post-settlement limitations due to space availability when the effective number of breeders is approaches saturating capacity, and causes a strong reduction of effective population size.

Early authors had already shown how changes in genetic diversity are related to effective population size [START_REF] Wright | Evolution in Mendelian populations[END_REF][START_REF] Kimura | Random genetic drift in multi-allelic locus[END_REF]). Nonetheless, it has been demonstrated theoretically that life cycle processes, like overlapping generations [START_REF] Jorde | Temporal allele frequency change and estimation of effective size in populations with overlapping generations[END_REF], and demographic parameters can also have an effect on the amount of temporal allele frequency change and should be taken into account [START_REF] Ryman | Minimizing adverse effects of fish culture: Understanding the genetics of populations with overlapping generations[END_REF]. Our results highlight the role of life cycle processes, other than age at first reproduction [START_REF] Lee | The influence of persistent individual differences and age at maturity on effective population size[END_REF], in reducing effective population size and increasing the effect of genetic drift. In closed populations, recruitment intensity regulates the speed of genetic drift on genetic diversity, demonstrating the relevance of transient demography in shaping the genetic structuration of marine species. In open populations, genetic drift is slowed down as a result of larval exchange. The effect of recruitment failure and the loss of key source populations on heterozygosity confirm that transient demographic fluctuations, remaining above the risk of regional extirpation, actually help maintain genetic diversity over time.

Recently, several meta-analyses and empirical studies evaluating the relationship between PLD and population genetic differentiation across different marine species indicated a weak correlation between the two [START_REF] Selkoe | Marine connectivity: A new look at pelagic larval duration and genetic metrics of dispersal[END_REF][START_REF] Faurby | Theoretical limits to the correlation between pelagic larval duration and population genetic structure[END_REF]. In some cases, the high uncertainty in the estimates of larval dispersal and measures of population differentiation were suggested as the main causes for the unexpected results, while other studies suggested that differential mortalities and levels of abundance among populations were more relevant [START_REF] Munroe | The role of larval dispersal in metapopulation gene flow: Local population dynamics matter[END_REF]. The results from our simulations support these findings, while further demonstrating that considering other biological parameters and the inclusion of transient dynamics caused by other demographic parameters, such as life expectancy, fecundity, local retention and recruitment success, would help explain better the variability of genetic connectivity in a metapopulation. [START_REF] Dawson | Dispersal potential and population genetic structure in the marine intertidal of the eastern North Pacific[END_REF], found that differences in fecundity, population size and PLD could account for the differences in population genetic structure between eight co-distributed rocky intertidal invertebrate species in the eastern North Pacific. They found that population size was positively related to population genetic structure and they attributed the variability among species to differences in recruitment and demography. Furthermore, they highlighted that the differences in migration potential found between theoretical and empirical data could be attributed to the role of genetic drift and natural selection. Here, we corroborate their hypotheses by demonstrating how the strength of genetic drift is influenced by recruitment success, local retention and fluctuations in population density due to habitat fragmentation and recruitment.

Our results confirm that genetic drift is controlled by effective reproduction, as it is observed when genetic diversity decreases due to a reduction of recruitment success caused by space limitation, rather than changes in the number of reproductive adults in the population. This is particularly the case when the filtering during reproduction applies to the fate of the larvae and not on the adults contributing to reproduction. We demonstrate that a long life expectancy and low recruitment success regulate the temporal variation of allele frequencies and reduce the loss of genetic diversity over time. This provides an alternative explanation to empirical evidence of strong genetic structure in long-lived species, when recruitment success is limited by space availability. In fact, it has already been shown how a long-lived species, such as the red coral

Corallium rubrum, can exhibit patterns of strong genetic differentiation at small spatial scales (tens of metres) in the Mediterranean Sea, suggesting that this could result from recruitment regulation by low population turnover [START_REF] Costantini | Genetic structuring of the temperate gorgonian coral (Corallium rubrum) across the western Mediterranean Sea revealed by microsatellites and nuclear sequences[END_REF].

Population regulation occurs when at least one demographic rate is density-dependent [START_REF] Hixon | Population regulation: Historical context and contemporary challenges of open vs. closed systems[END_REF]. Clear evidence of genetic drift accelerating due to recruitment regulation is observed in our simulations for closed populations at the initial saturating density. Furthermore, the evaluation of the ratio between mortality and recruitment success allowed for disentangling the relative contribution of intrinsic (intra-specific competition) vs. extrinsic drivers (inter-specific competition or predation).

In open systems, it was already suggested that the amount of local retention could have a strong influence on the relative importance of demographic vs. connectivity parameters in the dynamics of metapopulations [START_REF] Figueira | Connectivity or demography: Defining sources and sinks in coral reef fish metapopulations[END_REF]. Our results confirm those findings and indicate that local retention rate values larger than 2% diminish the influence of the connectivity structure over genetic drift, as genetic drift accelerates in those populations close to saturating capacity. This is particularly relevant for the spatial planning of MPAs. Whether an MPA is created for the conservation of biodiversity or for fishery management, its main purpose is to maintain the persistence of its populations. In the case of marine metapopulations, local retention is considered the currency of persistence: the higher the local retention, the closer the population is to its saturating capacity [START_REF] Burgess | Beyond connectivity: How empirical methods can quantify population persistence to improve marine protected-area design[END_REF]. Nonetheless, it is important to note that, as Figueira The maintenance of genetic diversity in the metapopulation over at least 100 years in the case of regular recruitment failure, creating regional demographic fluctuations, was not unexpected. These results, incorporating realistic spatially explicit connectivity matrices, extend previous theoretical studies that found that demographic instability and fluctuations in migration produce major changes in the evolution of genetic variation and population differentiation [START_REF] Whitlock | Temporal fluctuations in demographic parameters and the genetic variance among populations[END_REF]. Demographic disturbances due to population bottlenecks have been shown to induce strong genetic drift. In the present study, demographic disturbances around saturating capacity can reversely reduce the effect of genetic drift. Furthermore, our results suggests that there is a potential for opposing effects to population disturbances depending on the condition of the population prior to the disturbance: 1) accelerating genetic drift when population census size is high with recruitment limitation, or 2) decelerating genetic drift when population census size is low with no recruitment limitation.

Given that disturbances can cause variations in key demographic and biological processes, disturbance history may be one of the major drivers shaping the patterns of genetic diversity in many natural populations [START_REF] Banks | How does ecological disturbance influence genetic diversity?[END_REF]. The pattern of genetic diversity maintenance over time observed in our simulations has been reported in an empirical study on natural populations of frogs, where no losses of genetic diversity after a disturbance were found when survival was high or when recovering populations recruited many individuals from multiple sources [START_REF] Spear | Genetic structure among coastal tailed frog populations at Mount St. Helens is moderated by post-disturbance management[END_REF]. In addition, [START_REF] Larson | Measuring gene flow among populations having high levels of genetic fragmentation[END_REF] suggested, based on a survey of 22 species of salamanders, that historical influences are much more important in describing genetic variation patterns than is recent migration.

Previous studies in the region identified the spatial structuration of population vulnerability to habitat loss driven by hydrodynamical connectivity. Four sites essential for species regional persistence at minimum recruitment success were identified from metapopulation modelling around the port of Sète (sites 17-20, [START_REF] Guizien | Vulnerability of marine benthic metapopulations: Implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. However, analysis of the same connectivity graph used in the metapopulation model depicted those four sites as high bridging centrality, spanning across the two subclusters of the region. Another site (site 21) was instead identified for maintaining network integrity, as depicted by the highest betweeness (Costa, pers. comm). Here, we show that both the loss of four essential populations and the loss of the highest betweenness site divide the metapopulation into two genetically distinct subpopulations (eastern and western). Nevertheless, the loss of only site 21 leads to a higher level of genetic differentiation in the metapopulation. This result suggests that the highest betweeness of hydrodynamical connectivity graph could be interpreted for relaying gene flow and maintaining genetic connectivity, similar to how [START_REF] Rozenfeld | Network analysis identifies weak and strong links in a metapopulation system[END_REF] interpreted the highest betweenness of genetic distance. This example encourages further investigation of the potential of graph analysis of hydrodynamical connectivity matrices to identify sites important for conservation of genetic diversity.

On a local scale, high local retention, which helps maintain a species in isolated protected marine areas, accelerates genetic drift and as a result, accelerates the reduction of genetic diversity.

At the network level, the site having the most important role in heterozygosity maintenance (highest betweeness) did not correspond to those sites identified as essential for regional persistence in the metapopulation. The combination of these findings poses a question for management: in order to protect species persistence and genetic diversity, should we be seeking a compromise between preserving populations that are low in terms of diversity but are also demographically stable vs. highly diverse but unstable populations? The present model, by combining demography and gene flow in realistic conditions allows to 1) disentangle the relative contribution of local demography and environmental connectivity in shaping seascape genetics, and 2) perform in silico evaluations of different scenarios for marine spatial planning.

  Vector of spatial density of individuals at each site at time t+1 [individual per m 2 ] N i (t) = Vector of spatial density of individuals at each site at time t [individual per m 2 ] N max = Saturating population density based on space limitation [individual per m 2 ]

  Density of individuals in site j [individual per m 2 ] n ' jk i = Number of larvae produced per unit surface in site j having allele k that remain in the region and will survive until sexual maturity in site i n ' j i = Number of larvae produced per unit surface in site j that remain in the region and survive until sexual maturity in site i n jk = Number of larvae with allele k produced in site j per unit surface

(

  2009) suggested, under the same flow structure, differences in PLD among species could cause distinct patterns of local retention for different species in the same place, making biodiversity conservation efforts even more difficult. In any case, the loss of genetic variability through genetic drift at a high local retention rate can diminish future adaptability to a changing environment, which might be detrimental to local conservation efforts, strongly advocating for a network of MPAs formed by sites with lower local retention rate.

Figure 1 .

 1 Figure 1. Schematic representation of seascape genetic model. The frequency of each allele k is calculated at any given time in every site i, as a result of mixing between the surviving individuals with allele k in site I and to the input of new recruits with allele k, coming from any other site.

Figure 2 .

 2 Figure 2. Effect of different life expectancies on the loss of genetic diversity, measure as allele drifting time (T d ), for an isolated population at initial saturating capacity (6250 ind. m 22 ). Adimentionalizing allele drifting time by the number of reproductive events (T d /Re) (a) represents the differences in the number of generations it takes for a species with different life expectancies to start loosing its genetic diversity as a function of the ratio between mortality and recruitment success. In contrast, adimentionalizing allele drifting time by the life expectancy (T d /Le) of a species (b) provides a different vision corresponding to population turnover.

Figure 3 .

 3 Figure 3. Effect of population growth rate over allele fixation time (T f ) for three different initial conditions of population abundance (50, 1000, and 6250 ind. m -2 ) as a function of species life expectancy (Le) in a closed population.

Figure 4 .

 4 Figure 4. Effect of different initial population abundances (50, 1000, and 6250 ind. m -2 ) on the loss of genetic diversity, measured as allele drifting time (T d ), as a function of (a) the ratio between mortality and recruitment success [(1-s)/b], and (b) the time to reach saturating capacity (T). T d is adimensionalized by the life expectancy (T d /Le) of a species that lives 10 years, thus representing the number of population turnovers to start loosing genetic diversity in an isolated population.

Figure 5 .

 5 Figure 5. Spatial pattern of genetic diversity caused by distinct connectivity structures. Variability of heterozygosity for 32 sites forming a metapopulation for a species with a life expectancy of 5 years.

Figure 6 .

 6 Figure 6. Relationship between local retention (LR; defined as the proportion of larvae produced in one population that remain in the same population) and the number of population turnovers to start loosing genetic diversity (T d /Le), in a metapopulation over time. Values correspond to 20 connectivity matrices, and species with three different life expectancies (2, 5, and 20 years). Dashed arrow indicates the threshold at which the effect of LR on T d /Le begins to decrease.

Figure 7 .

 7 Figure 7. Effect of fluctuations in local population stability on heterozygosity over time. Scenarios: (a) Simulates gene flow among all of the populations forming a metapopulation with constant recruitment success every year. (b) Simulates recruitment failure (no recruitment) every 3 years. (c) Simulates the loss of four source sites (sites 17-20) in the metapopulation, while keeping recruitment success constant every year. (d) Simulates the loss of one key population (site 21) for maintaining genetic connectivity in the metapopulation.

Table 1 .

 1 Settings of six simulation groups. Connectivity, when present, was set either deterministically (D) using the same matrix every year (either one of the 20 variants) or stochastically (S) picking a different matrix among the 20 variants every year. Disturbance scenarios were only applied to the groups with stochastic connectivity.

		Group 1	Group 2 Group 3A Group 3B	Group 4A	Group 4B
	Connectivity	NO	D	S	S	S	S
	Initial population density	50, 1000, Nmax	50	50	50	50	50
	(Individuals/m 2 )

Acknowledgments

This work was (co-)funded through a MARES Grant. MARES is a Joint Doctorate programme selected under Erasmus Mundus coordinated by Ghent University (FPA 2011-0016). Check www.mares-eu.org for extra information. This work was also partly funded by the French National Program LITEAU IV of the Ministère de l'Ecologie et de l'Environnement Durable under project RocConnect-Connectivité des habitats rocheux fragmentés du Golfe du Lion (PI, K. Guizien, Project Number 12-MUTS-LITEAU-1-CDS-013). We would like to thank Proof-readingservice.com for English edition and corrections of the manuscript, and Dr. Lorenzo Bramanti for an early critical review of the manuscript.

discrimination between two groups (Fig. 8) -populations with low levels of local retention (< 2%) that show high variability in the time they begin to loose their genetic diversity (Td/Le ranging from 0.2 to 50) and populations with high local retention (> 2%) where genetic drift is much quicker and less variable (Td/Le ranging from 0.2 to 12). This indicates that the speed of genetic drift is regulated locally only at a high local retention rate (larger than 2%) while, for local retention rates lower than 2%, the strength of genetic drift depends on the regional structure of connectivity. According to the direction of the transfer and the genetic characterisation of the contributing populations, genetic drift will then be accelerated or slowed down in the focal population.

How does local population stability affect heterozygosity in a metapopulation?

When looking into the variation of the average heterozygosity in the metapopulation through time for a species with a life expectancy of 10 years, it becomes clear that changes in population density have a major effect on the loss of genetic diversity in a metapopulation (Fig. 9).

The rate at which a metapopulation reaches a stable regional coverage influences its possibility to maintain genetic diversity over time, even in the presence of migration.

Overall, we observe that as soon as regional coverage reaches 80%, heterozygosity values decrease rapidly and genetic diversity is lost. Scenario A shows how, with constant recruitment success, the effect of migration among populations homogenises allele frequencies over time and leads to the fixation of the same allele in all populations over a relatively short period of time. Most interestingly, in Scenario B, we find that high variability in metapopulation regional coverage due to recurrent failures in recruitment actually helps maintain genetic diversity over time. Stochastic connectivity leading to demographical fluctuations in populations enables the maintenance of genetic diversity.

On the other hand, Scenario C exemplifies how, in the Gulf of Lions, removing