E Guilmineau
email: emmanuel.guilmineau@ec-nantes.fr

M Visonneau
email: michel.visonneau@ec-nantes.fr

FP15-AERO2022-guilmineau EVALUATION OF TURBULENCE MODELLING FOR SUBSONIC VORTEX INTERACTION ON MULTI DELTA WING CONFIGURATION

In the present work, numerical simulations of the flow around an aircraft at several angles of attack and with a moderate drift angle are performed with the ISIS-CFD flow solver which is an incompressible unsteady Reynolds-averaged Navier-Stokes method. The turbulence modelling is performed using the linear isotropic k-ω SST model and the hybrid RANS/LES DDES-SST model. For every turbulence closure, the anisotropic adaptive grid refinement functionality is used. Aerodynamic coefficients as well as flow field data are compared with experimental measurements. The main conclusion is that the numerical simulations with the DDES-SST model are in better agreement with the experimental data than the simulations with the k-ω SST.

INTRODUCTION

High-agile high performance aircraft configurations have to cover a wide scope of performance and maneuverability requirements from an supersonic cruise to a high maneuverability at subsonic conditions [START_REF] Hitzel | Combat Aircraft Vortex Interaction -Design, Pysics & CFD-Tools[END_REF]. Hybrid-delta wings present many assets for the aerodynamic design of high agility military aircraft. The flow over such wings is dominated by multiple vortices. When the angle of attack increases, the vortices start interacting with each other enabling a more stable vortex system. For high angles of attack, the vortex breaks down resulting in an increasing pitching moment and a decrease in the gradient of lift coefficient [START_REF] Breitsamter | Unsteady flow phenomena associated with leading-edge vortices[END_REF]. An extensive synthesis of CFD simulations for the F16-XL configurations, which has a triple-deltawing configuration, has been conducted [START_REF] Luckring | a Synthesis of Hybrid RANS/LES CFD Results for F-16XL Aircraft Aerodynamics[END_REF] using hybrid RANS/LES methods which predicted correctly the vortical characteristics.

A low-aspect ratio multiple-swept wing fuselage is studied in the NATO STO AVT-316 task group on "Vortex Interaction Effects Relevant to Military Air Vehicle Performance". This collaborative research has been initiated to assess Computational Fluid Dynamics (CFD) methods to predict vortex-interaction effects and extend the understating of vortex-interaction flow physics for these problems through numerical and physical experimentation.

In this paper, we conduct a numerical simulation of the flow past the wing-fuselage combination at Mach number 0.15 for several angles of attack with a side slip angle, by means of RANS and hybrid RANS/LES simulations. A comparison will be carried out between the numerical simulations and experimental data.

TEST CASE

Figure 1 illustrates the test configurations [START_REF] Hitzel | Vortex Flow Aerodynamic Challenges in the Design Space for Future Fighter Aircraft[END_REF][START_REF] Pfnür | Leading-Edge Vortex Interactions at a Generic Multiple Swept-Wing Aircraft[END_REF]. The first one is a triple delta wing configuration, the socalled NA1-W1, see Figure 1(a), with three different consecutive wing sections featuring a varying leading-edge sweep. The inboard and outboard wing sections exhibit a medium leading-edge sweep of ϕ 1 = ϕ 3 = 52.5°and the mid-board wing section is highly swept at ϕ 2 = 75.0°. The second one is a double delta configuration, the socalled NA1-W2, see Figure 1(b), with two different consecutive wing sections of varying leading-edge sweeps characterized by a high leading-edge sweep of ϕ 2 = 75.0°i n the inboard wing section and a medium leading-edge sweep of ϕ 3 = 52.5°in the outboard wing section. Both configurations have the same total length of l tot = 1.160 m and the same root chord of c r = 0.802 m.

Experiments are conducted in the wind tunnel of Technical University of Munich (TUM) for a speed of V ∞ = 51.97 m/s, which leads to a Reynolds number, based on c r , of 2.36 × 10 6 and a Mach number of 0.15. Two side

FLOW SOLVER

The in-house solver ISIS-CFD developed by the authors, also available as a part of the FINE™/Marine computing suite worldwide distributed by Cadence Design Systems, is an incompressible multiphase unsteady Reynoldsaveraged Navier-Stokes (URANS) solver mainly devoted to marine hydrodynamics. It is based on a fully-unstructured (face-based) finite volume discretization with specific functionalities needed for multiphase flows and industrial applications [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF][START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF]. The method features several sophisticated turbulence models: apart from the classical two-equation k-ε and k-ω models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available [START_REF] Cecora | Differential Reynolds-Stress Modeling for Aeronautics[END_REF][START_REF] Deng | Comparison of explicit algebraic stress models and second-order turbulence closures for steady flow around ships[END_REF][START_REF] Duvigneau | On the role played by turbulence closures in hull ship optimization at model and full scale[END_REF]. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid RANS/LES turbulence models based on Detached Eddy Simulation (DES-SST, DDES-SST, IDDES-SST) are also implemented and have been thoroughly validated on automotive flows characterized by large separations [START_REF] Guilmineau | Assessment of hybrid RANS-LES formulations for flow simulation around the Ahmed body[END_REF] and ships at steady drift [START_REF] Visonneau | Local flow around a surface combatant at various static drift conditions: The role played by turbulence closures[END_REF]. Moreover, the solver accepts sliding and overset grids and features an anisotropic adaptive grid refinement functionality [START_REF] Mozaffari | Average-based mesh adaptation for hybrid RANS/LES simulation of complexe flows[END_REF][START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in freesurface flow simulation[END_REF] applied to unstructured hexahedral meshes.

NUMERICAL SIMULATION SET-UP

For this study, only the side-slip condition β = 5°is investigated. The turbulence models used are k-ω SST 2003 for the RANS model and DDES-SST for the hybrid RANS/LES modelisation. The time step used with the k-ω SST model is ∆t = 1.5×10 -4 s and ∆t = 5×10 -5 s with the DDES-SST model. The averaging time, tU ∞ /c r , varies between 30 and 60 depending on the case.

The meshes are generated using Hexpress™, an automatic unstructured mesh generator. This software generates meshes containing only hexahedrals. For the aircraft and the sting, no-slip boundary conditions are imposed and the wall normal resolution has y + below 1.

During the simulation, the average-based mesh adaptation procedure, included in ISIS-CFD, is also used, choosing as refinement criterion the flux component Hessian. The grid is adapted every 10 time steps. As the flow depends on the turbulence model, the number of cells varies depending on the turbulence model and the angle of attack α. Thus the number of cells varies between 65M and 130M. Fig. 2 shows an example of the mesh at cross-section x/c r = 0.592 for α = 24°and β = 5°for the geometry NA1-W2.

RESULTS

Forces and moments

Fig. 3 presents the evolution of the longitudinal aerodynamic coefficients (the drag coefficient C D , the lift coefficient C L and the pitching moment coefficient C my ) versus the angle of attack, α, for both configurations. For the NA1-W1 geometry, see Fig. 3(a), the numerical results obtained with the k-ω SST and DDES-SST turbulence models are in good agreement with the experimental data. In the experiments, the vortex breakdown occurs at α = 15°. For the NA1-W2 geometry, see Fig. 3(b), these aerodynamic coefficients are well predicted by both The lateral aerodynamic coefficients (the side-force coefficient C S , the pitching moment coefficient C mx and the yawing moment coefficient C mz ) are presented in Fig. 4. For the NA1-W1 geometry, see Fig. 4(a), the agreement between the experimental data and the numerical results is less satisfactory than that obtained for the longitudinal aerodynamic coefficients. The side-force coefficient C S and the yawing moment coefficient C mz are relatively well predicted by both turbulence models while the prediction of the rolling moment coefficient C mx is inaccurate. The k-ω SST model is unable to predict the decay and the change of the slope observed experimentally around 16°. The DDES-SST model predicts the correct negative value at α = 16°. For α = 18°, C mx increases as in experiments but for the higher angle of attack, the value decreases instead of keeping at constant value close to zero as in the experiments. For the very high value of angle of attack, α = 32°, the two turbulence models predict the positive value of C mx , while C S is hidden under C mx for DDES-SST. For the NA1-W2 geometry, see Fig. 4(b), both turbulence models predict coefficients in agreement with the experimental data for the angles of attack α below 20°. For the higher angles, α ≥ 24°, the k-ω SST model is unable to predict the monotonous decay of C mz and for the C S prediction, the value is in good agreement with the experimental data for α = 24°but it is strongly overestimated at 32°. As with the NA1-W1 wing profile, the behaviour of the rolling moment, with its change of sign around 20°, is not predicted with the k-ω SST model which predicts a decay. With the DDES-SST model, all lateral aerodynamic coefficients are well predicted and in good agreement with the experimental data, notably the change of sign of the rolling moment.

Angle of attack α = 16°F

ig. 5 shows the vortical structures for NA1-W1 and NA1-W2 geometries. The vortex developing in the front 0.592 is presented in Fig. 6 for the NA1-W1 geometry.

In the experiments, for the NA1-W1 geometry on the leeward side, the IBV and MBV vortices have the characteristics of a jet-type vortex with a high axial velocity in the core of these vortices. On the windward side, the MBV vortex is also a jet-type vortex while the IBV has the characteristics of a wake-type vortex. For the numerical simulations, at windward side, the DDES-SST model predicts the low level of axial velocity and thus the wake-type vortex. It is not the case with the k-ω SST model. At leeward side, no turbulence model is able to predict the correct level of axial velocity in the IBV core. Both turbulence closures show two vortices very close to each other with opposite axial velocity. This feature is not observed in the experiments. Fig. 7 presents the dimensionless axial velocity distribution at x/c r = 0.592 for the NA1-W2 geometry. For this geometry, the experiments show more symmetric vortices on each wings contrary to the NA1-W1 geometry. All vortices on the leeward side and the windward side have the jet-type characteristics. The k-ω SST model predicts a symmetric solution as in the experiments. However, the velocity in the IBV and MBV vortex core are lower than that obtained in the experiments. With the DDES-SST model, the level of axial velocity in the IBV vortex core at the leeward side is in agreement with the experimental data but at the windward side, the level is lower than the experiments. Thus the DDES-SST model predicts a asymmetric behaviour.

Fig. 8 shows the dimensionless axial velocity in the IBV and MBV vortices on the leeward side. In the experiments for the NA1-W1 geometry, the IBV vortex has the characteristics of a jet-type vortex with a plateau of high velocity in the vortex core over the total length of the wing. The k-ω SST model predicts this plateau but the level is twice as small as the experimental value. With the DDES-SST model, the velocity starts with a value around 2.0 and then decreases until x/c r ∼ 0.6 to reach a value around 0.5. For the geometry NA1-W2, the IBV vortex is again characterised by a plateau with an unchanged level of core velocity. Both numerical simulations predict a plateau but for this geometry, the DDES-SST model is able to predict the correct level of velocity, contrary to the k-ω SST model. Fig. 9 shows the dimensionless axial velocity in the IBV and MBV vortices on the windward side. For the NA1-W1 geometry, until x/c r ∼ 0.6, the experiments depict a decay of the IBV core velocity and then a slow increase until the trailing edge of the wing. The DDES-SST model predicts very well the behaviour of the IBV core velocity. The k-ω SST does not correctly predict this decrease. For the NA1-W2 geometry, the experiments reveal again a plateau with a high core velocity. Both numerical simulations are not able to predict this plateau. Instead, a decrease of the core velocity is predicted.

Angle of attack α = 24°F

ig. 10 shows the vortical structures for both geometries, NA1-W1 and NA1-W2, at 24°. For the geometry NA1-W1, the DDES-SST model predicts more intense vortical structures on both leeward and windward wings. For this angle of attack, a flow reversal is observed in the inboard vortex, all along its trajectory on the windward wing while flow reversal occurs only close to the trailing edge of the leeward wing for both turbulence closures.

For the geometry NA1-W2, the DDES-SST model predicts a more intense and stable IBV vortex on the leeward wing, leading to a more asymmetric vortical structures. Fig. 11 depicts the dimensionless axial velocity distribution at x/c r = 0.592 for the NA1-W1 wing configuration. For the NA1-W1 configuration, the TUM experiments reveal an asymmetry of the vortices developing on each wing. On the leeward side, both IBV and MBV are clearly jet-type vortices with a strong axial velocity in the core of these vortices. The IBV vortex is similar between both numerical simulations but the jet-type vortex is not captured. Both turbulence closures predict the existence of two zones very close to each other corresponding to two co-rotating vortices, the upper one with a velocity overshoot and the lower one with a velocity undershoot. Only the DDES-SST model is able to reproduce the behaviour of the MBV vortex. On the windward side, the IBV vortex observed experimentally is a wake-type vortex while the MBV vortex is a jet-type vortex. Both turbulence closures are able to capture the strong decrease of the axial velocity in the core of the IBV vortex. Moreover, both numerical simulations slightly overpredict the core velocity in the MBV vortex.

Fig. 12 presents the dimensionless axial velocity distribution at x/c r = 0.592 for the NA1-W2 wing configuration. As for the NA1 W1 geometry, the TUM experiments show an asymmetric situation characterised for the IBV vortex by a jet-type vortex on the leeward side and a wake-type vortex on the windward side while the MBV vortices keep their jet-type configuration. For the numerical simulations, at windward side, the level of the axial velocity in the IBV, obtained the k-ω SST model, is not in agreement with the experiments. Moreover, this model predicts an IBV vortex surrounded by a high axial velocity. With the DDES-SST model, the area of the low velocity in the IBV is lower than that obtained in the experiments. At the leeward side, the IBV vortex is under-estimated with the k-ω SST model. The DDES-SST model captures the right asymmetry of the flow even if the area of the high velocity is lower than that observed experimentally.

Fig. 13 shows the dimensionless axial velocity in the cores of the IBV vortex and MBV vortex on the leeward side. For the NA1-W1 geometry, both numerical simulations are not able to predict the high velocity during the IBV progression until x/c r ∼ 0.6. At the end of the wing, the DDES-SST model provides a lower core axial velocity close to zero. For the NA1-W2 geometry, the IBV vortex is a jet-type vortex in the experiments and the numerical simulations reproduce this feature. However the k-ω SST model underestimates the axial core velocity and this plateau is maintained up to x/c r = 0.6, and downstream, the velocity decreases. With the DDES-SST model, the velocity in the core of the IBV vortex is overestimated close to the onset of the wing. However, the high velocity plateau is maintained up to the trailing edge of wing and it is in very good agreement with the experimental data.

The dimensionless axial velocity in the cores of the IBV vortex and MBV vortex on the windward side is presented in Fig. 14. For the NA1-W1 geometry, the experiments show a flow reversal occurring in the core of the IBV vortex. This behaviour is well captured by the two numerical simulations. The MBV core velocity is overestimated by both turbulence closures. For the NA1-W2 geometry, the vortex breakdown of the IBV in the exper- iments occurs for x/c r around 0.5 and it is followed by a region of negative velocity up to the trailing edge of the wing. The k-ω SST model underestimates the axial velocity until the vortex breakdown which is somewhat smeared, compared to the experiments. The DDES-SST model provides a globally better prediction of the axial velocity in the IBV vortex even if the velocity is not negative after the vortex breakdown.

Angle of attack α = 32°F

ig. 15 shows the vortical structures for both geometries and both turbulence closures. For the NA1-W1 geometry, a massive separation at the windward side is predicted by both turbulence models. However, this separation appears more intense with the DDES-SST model. At the leeward side, the result depends on the turbulence closure used. With the DDES-SST model, the flow reversal in the IBV vortex starts at mid-wing while with the k-ω SST, it starts close to the onset of the wing. For the NA1-W2 geometry, the two numerical simulations provide a different solution. At the windward side, the DDES-SST model predicts a massive separation as for the NA1-W1 geometry while the k-ω SST model predicts a longitudinal vortex with a zone of flow reversal. On the leeward wing, both numerical simulations predict longitudinal vortices as for the previous angle of attack and the vortices ob- The dimensionless axial velocity in the cores of the IBV vortex and MBV vortex on the leeward side is presented in Fig. 18. For the NA1-W1 geometry, the experiments show a high velocity at the onset of the wing which decreases drastically until x/c r around 0.3 reaching zero. After this position the velocity keeps a negative value. Both numerical simulations are not able to predict this high velocity at the onset. However, the DDES-SST model predicts very well the plateau with a zero or slightly negative velocity in the vortex core while the model k-ω SST predicts a value that oscillates between 0.5 and 1.0. For the NA1-W2 geometry and the IBV vortex, the experiments show a plateau with a high core velocity until x/c r ∼ 0.7 before a slow decay until the trailing edge of the wing. The DDES-SST model predicts the same behaviour but the core velocity is overestimated. With the k-ω SST model, a plateau is predicted but the core velocity is underestimated and the decay starts at x/c r ∼ = 0.55 to reach a velocity 0.5 while in the experiments, the measured value is 2.0.

CONCLUSIONS

In this paper, a study of the flow around two different wing configurations, the so-called NA1-W1 and NA1-W2 geometries, for several angles of attack with a drift angle, has been presented. Two turbulence closures were used: the k-ω SST 2003 and the hybrid RANS/LES DDES-SST models. Numerical simulations were run using flow-adapted grids to capture automatically the vortical structures during their evolution. This functionality is a very precious specificity of the flow solver ISIS-CFD for this category of compoutations. Comparisons of aerodynamic coefficients, local flow details and vortex core analysis were conducted with the TUM experiments. The main conclusion is that the NA1-W2 geometry appears easier to simulate than the NA1-W1 geometry with the available turbulence closures, numerical simulations being in better agreement with the experiments and DDES-SST model being the most reliable turbulence closure. The change of sign of the rolling moment is indeed a challenge for turbulence models. For the NA1-W2 geometry, the DDES-SST model is the only one able to predict this change of sign around the angle of attack α = 20°. This change of sign of the rolling moment is due to an asymmetric behaviour of longitudinal vortices on the leeward and windward sides. For the NA1-W1 geometry, the rolling moment change of sign is not satisfactorily predicted by both turbulence closures. This incorrect prediction is probably due to the early interaction between the fuselage vortex and the IBV vortex, not far from the onset of the wing. Finally, while the use of hybrid RANS/LES models is the only viable option to accurately predict the change of rolling stability, the computations should obligatorily be performed with anisotropic adaptive mesh refinement, since this is the only possibility to maintain the fine grid density in the core of the various interacting vortices and keep the computational effort for industrial studies within high but reasonable limits.

ACKNOWLEDGMENT

The computations were performed using HPC resources from GENCI (Grand Equipement National de Calcul Intensif) which is gratefully acknowledged.

  (a) Model NA1-W1. (b) Model NA1-W2.

Figure 1 :

 1 Figure 1: Wing planforms.

  (a) k-ω SST. (b) DDES-SST.

Figure 2 :

 2 Figure 2: NA1-W2: α = 24°, β = 5°-View of the flow adapted grids at the cross-section x/c r = 0.592.

  (a) NA1-W1 geometry. (b) NA1-W2 geometry.

Figure 3 :

 3 Figure 3: Longitudinal aerodynamics coefficients versus the angle of attack α for both geometries.

  (a) NA1-W1 geometry. (b) NA1-W2 geometry.

Figure 4 :

 4 Figure 4: Lateral aerodynamics coefficients versus the angle of attack α for both geometries.

  (a) NA1-W1: k-ω SST. (b) NA1-W1: DDES-SST. (c) NA1-W2: k-ω SST. (d) NA1-W2: DDES-SST.

Figure 5 :

 5 Figure 5: α = 16°-Vortical structures for the NA-W1 and NA1-W2 configurations.

Figure 6 :

 6 Figure 6: NA1-W1: α = 16°, β = 5°-Axial component of the velocity at x/c r = 0.592.

Figure 7 :

 7 Figure 7: NA1-W2: α = 16°, β = 5°-Axial component of the velocity at x/c r = 0.592.

  (a) NA1-W1 geometry. (b) NA1 W2 geometry.

Figure 8 :

 8 Figure 8: Leeward wing: α = 16°, β = 5°-Dimensionless axial velocity for the IBV and MBV of the NA1-W1 and NA1-W2 configurations.

  (a) NA1-W1 geometry. (b) NA1 W2 geometry.

Figure 9 :

 9 Figure 9: Windward wing: α = 16°, β = 5°-Dimensionless axial velocity for the IBV and MBV of the NA1-W1 and NA1-W2 configurations.

  (a) NA1-W1: k-ω SST. (b) NA1-W1: DDES-SST. (c) NA1-W2: k-ω SST. (d) NA1-W2: DDES-SST.

Figure 10 :

 10 Figure 10: α = 24°-Vortical structures for the NA-W1 and NA1-W2 configurations.

  (a) k-ω SST. (b) DDES-SST. (c) TUM Experiments.

Figure 11 :

 11 Figure 11: NA1-W1: α = 24°, β = 5°-Axial component of the velocity at x/c r = 0.592.

  (a) k-ω SST. (b) DDES-SST. (c) TUM Experiments.

Figure 12 :

 12 Figure 12: NA1-W2: α = 24°, β = 5°-Axial component of the velocity at x/c r = 0.592.

  (a) NA1-W1 geometry. (b) NA1 W2 geometry.

Figure 13 :

 13 Figure 13: Leeward wing: α = 24°, β = 5°-Dimensionless axial velocity for the IBV and MBV of the NA1-W1 and NA1-W2 configurations.

  (a) NA1-W1 geometry. (b) NA1 W2 geometry.

Figure 14 :

 14 Figure 14: Windward wing: α = 24°, β = 5°-Dimensionless axial velocity for the IBV and MBV of the NA1-W1 and NA1-W2 configurations.

  (a) NA1-W1: k-ω SST. (b) NA1-W1: DDES-SST. (c) NA1-W2: k-ω SST. (d) NA1-W2: DDES-SST.

Figure 15 :

 15 Figure 15: α = 32°-Vortical structures for the NA-W1 and NA1-W2 configurations.

Figure 16 :

 16 Figure 16: NA1-W1: α = 32°, β = 5°-Axial component of the velocity at x/c r = 0.592.

Figure 17 :

 17 Figure 17: NA1-W2: α = 32°, β = 5°-Axial component of the velocity at x/c r = 0.592.

Figure 18 :

 18 Figure 18: Leeward wing: α = 32°, β = 5°-Dimensionless axial velocity for the IBV and MBV of the NA1-W1 and NA1-W2 configurations.


Sylvain Blouet
email: sylvain.blouet@ville-agde.fr

Lorenzo Bramanti

Katell Guizien

Artificial reefs geographical location matters more than shape, age and depth for sessile invertebrate colonization in the Gulf of Lion (NorthWestern Mediterranean Sea)

come   L'archive ouverte pluridisciplinaire

Introduction

The decline of fish stocks and natural marine habitat degradation resulting from human exploitation have been documented worldwide for decades [START_REF] Claudet | Human-driven impacts on marine habitats: A regional meta-analysis in the Mediterranean Sea[END_REF]Jackson, 2001;[START_REF] Pauly | Towards sustainability in world fisheries[END_REF].

Artificial reefs (ARs) have been primarily implemented to reduce the pressure of fisheries in coastal areas, complementing other management tools such as marine protected areas or regulatory measures such as fishing licenses [START_REF] Claudet | Marine protected areas and artificial reefs: A review of the interactions between management and scientific studies[END_REF][START_REF] Seaman | Artificial habitats and the restoration of degraded marine ecosystems and fisheries[END_REF][START_REF] Wilson | Restoration of Hong Kong fisheries through deployment of artificial reefs in marine protected areas[END_REF]. Moreover, ARs could provide economic benefits such as recreational and traditional fishing and scuba diving [START_REF] Chen | Recreational Benefits of Ecosystem Services on and around Artificial Reefs: A Case Study in Penghu, Taiwan[END_REF]. Beneficial effects such as increase in fish biomass and capture efficiency near ARs have been reported (reviewed by [START_REF] Bohnsack | Artificial reef research: A review with recommendations for future priorities[END_REF][START_REF] Tessier | A case study of artificial reefs as a potential tool for maintaining artisanal fisheries in the French Mediterranean Sea[END_REF] but led to a debate on the effects of ARs fishery, opposing attraction vs production [START_REF] Grossman | Do Artificial Reefs Increase Regional Fish Production? A Review of Existing Data[END_REF]. The fish attraction argument is based on the quick colonization by fishes and mobile invertebrates [START_REF] Powers | Estimating enhancement of fish production by offshore artificial reefs: uncertainty exhibited by divergent scenarios[END_REF][START_REF] Relini | Trophic relationships between fishes and an artificial reef[END_REF][START_REF] Santos | A fourteen-year overview of the fish assemblages and yield of the two oldest Algarve artificial reefs (southern Portugal)[END_REF]. The fish production argument is based on the hypotheses of a better protection against predators and an increase in available substrate area for larval establishment thanks to habitat complexification and an increase of available trophic resource [START_REF] Bohnsack | Are High Densities of Fishes at Artificial Reefs the Result of Habitat Limitation or Behavioral Preference?[END_REF]. In natural rocky habitats, benthic invertebrates play an essential role in fish trophic networks [START_REF] Ardizzone | Colonization and disappearance of Mytilus galloprovincialis Lam. on an artificial habitat in the Mediterranean sea[END_REF][START_REF] Martens | Marine Biodiversity: Patterns and Processes, Assessment, Threats, Management and Conservation[END_REF]. ARs trophic network showed similarity with natural rocky habitat, with dominance of filter-feeders using phytoplanktonic primary production and fish predation on crustacean colonizing the ARs (Cresson, 2013). Moreover, AR deployed in sandy areas are expected to enhance fish productivity given that epifauna secondary production per ARs unit surface has been estimated to be 30 times greater than that of natural sandy infauna [START_REF] Steimle | Benthic macrofauna productivity enhancement by an artificial reef in Delaware Bay, USA[END_REF]. However, those estimates were made shortly after immersion and do not prove the long-term fish production in ARs and supporting this argument would require extending data on colonization by benthic invertebrates in the long-term [START_REF] Svane | On the Problems of Epibioses, Fouling and Artificial Reefs, a Review[END_REF]. Indeed, age since deployment has been described as a key factor to explain ARs coverage by benthic invertebrates [START_REF] Svane | On the Problems of Epibioses, Fouling and Artificial Reefs, a Review[END_REF]. The assemblages of benthic communities are expected to change over time in a succession between pioneer and specialist species [START_REF] Connell | Mechanisms of Succession in Natural Communities and Their Role in Community Stability and Organization[END_REF]. In contrast with pioneer species, specialist ones have slower colonization dynamics, because of their lower fecundity [START_REF] Fava | Role of Recruitment Processes in Structuring Coralligenous Benthic Assemblages in the Northern Adriatic Continental Shelf[END_REF]. However, after colonization, specialists are expected to outcompete pioneer species due to their more efficient use of environmental resources [START_REF] Connell | Mechanisms of Succession in Natural Communities and Their Role in Community Stability and Organization[END_REF]. Among these, light availability is an essential factor shaping marine benthic communities across the water depth gradient (Odum, 1971). Several studies have shown a decrease in the density of benthic invertebrates with depth on ARs (Lewbel et al., 1986;[START_REF] Moura | Effect of depth and reef structure on early macrobenthic communities of the Algarve artificial reefs (southern Portugal)[END_REF]Shinn and Wicklund, 1989;[START_REF] Van Der Stap | Marine Fouling Assemblages on Offshore Gas Platforms in the Southern North Sea: Effects of Depth and Distance from Shore on Biodiversity[END_REF] explained by the decrease in light intensity (Relini et al., 1994). The structural complexity has also been put forward as important characteristics linked to ARs efficiency in ecological restoration [START_REF] Strain | Eco-engineering urban infrastructure for marine and coastal biodiversity: Which interventions have the greatest ecological benefit?[END_REF]. Structural complexity increases available surface for colonization and niches diversity with various shelter and light exposure conditions, the latter being related to different benthic assemblage compositions [START_REF] Glasby | Surface composition and orientation interact to affect subtidal epibiota[END_REF]Glasby, 1999) and higher productivity [START_REF] Vivier | Marine artificial reefs, a meta-analysis of their design, objectives and effectiveness[END_REF]. The recent 3D printing techniques using concrete allow the design of ARs mimicking natural habitats [START_REF] Ly | Optimisation of 3D printed concrete for artificial reefs: Biofouling and mechanical analysis[END_REF]. However, those studies concerned short-term colonization (<3.5 years) (Wendt et al., 1989) hence based on pioneer species with high dispersal capacities which long term colonization is likely regulated by post-settlement processes such as competition, predation and physical disturbance [START_REF] Todd | Larval supply and recruitment of benthic invertebrates: do larvae always disperse as much as we believe?[END_REF]. In contrast to mobile species, sessile benthic invertebrates can only colonize reefs after larval dispersal which is limited by reproduction frequency [START_REF] Thorson | Reproductive and larval ecology of marine bottom invertebrates[END_REF]. Colonization implies thus an effective dispersal between natural areas and ARs, which depends on source population spatial distribution, species fecundity, dispersive larval traits and ocean circulation. Nevertheless, until now dispersal drivers have been disregarded while colonization disparities among ARs may result from differences in both larval connectivity (which in turn depends on fecundity, dispersal capacities and adult distribution in the natural habitat) and post-recruitment processes.

The objective of the present study was to test the hypothesis that the geographical location of ARs deployment with respect to the natural habitat can condition ARs colonization after more than 10 years. To this aim we investigated the effects of local (shape, depth and immersion duration) and regional (geographic area) factors on the presence of five species of sessile invertebrates with different life history traits, endemics to the Gulf of Lion (GoL) (Northwestern Mediterranean Sea) and frequently found on natural hard substrates.

Methods

Study area and spatially stratified sampling design

The study area extended along 160 km of the GoL coastline (Figure 1). The GoL is a wide micro-tidal continental shelf dominated by soft-bottom habitat with few small rocky habitat patches of less than 20 km 2 . The GoL is a homogeneous and isolated hydrodynamic unit [START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF], delimited by the northern current [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF]. Between 1985 and2009, 763 ARs with different shapes or material for a total volume of 37 575 m 3 [START_REF] Tessier | Assessment of French artificial reefs: due to limitations of research, trends may be misleading[END_REF] have been deployed along the GoL coastline over 66 km² of state concessions (Cepralmar, 2015) between 10 and 35 m depth (Figure 1, [START_REF] Blouet | Database of the location and typology of artificial reefs in the Gulf of Lion (NW Mediterranean Sea)[END_REF]. The ARs deployed in GoL represent 40% of the total AR volume in France [START_REF] Tessier | Assessment of French artificial reefs: due to limitations of research, trends may be misleading[END_REF].

In the present study we examined the oldest (in the 80s) and the youngest (in the 2000s) ARs with the most common shapes (namely pipe, steel cage, Bonna, Comin, pyramid, heap of telegraph poles, and concrete box, Figure 1). For the analysis of data, we followed a stratified sampling design.

To this aim, the GoL coastline was regularly divided into 6 geographical sectors separated by a distance ranging from 12 to 117 km, with a median value of 49 km (AGM referring to Aigues-Mortes, AGD to Agde, VLR to Valras, GRU to Gruissan, LEU to Leucate, CST to Canet/Saint-Cyprien sectors, Figure 1). Each sector included ARs deployed during either the first (1985, CST geographical sector), the last (2002-2009, AGM, VLR, GRU, LEU geographical sectors) or both deployment periods (AGD), and of different shapes. In each geographical sector, two sites were defined except in AGD (5 sites, AGD1, AGD2, AGD3, AGD4, AGD5) and GRU (3 sites, GRU1, GRU2, GRU3). To test for location effect, sites were separated by a distance ranging from 2.1 to 11.7 km (median value of 7.5 km) except in GRU, where two of the three sites (GRU1 and GRU2) were defined in the same location to account for different shapes. Hence, depending on the geographical sector, sites may differ either by AR shape, immersion depth or immersion duration (Figure 1). Due to this set up, the effect of these three factors was only tested locally (see Statistical analysis). In each site (except in the geographical sector GRU), three sampling units separated by a distance between 4 m and 3.6 km (median distance of 251 m) were set out by pooling neighboring ARs reaching a minimum surface of 89 m 2 per sampling unit and totalling a minimum developed surface of 306 m 2 per site (Table 1: Supplementary Material). In the geographical sector GRU, GRU2 and GRU3 included only one sampling unit because the surface of a single AR in these sites already yielded 459 m 2 (heap of telegraph poles; Table 1: Supplementary Material). Such large continuous sampling units in each site aimed at limiting the effect of recruitment spatial variability over distances from 100s meters to kilometers arising from the nonuniformity of the flow of larvae [START_REF] Daigle | Fine-Scale Distribution and Spatial Variability of Benthic Invertebrate Larvae in an Open Coastal Embayment in Nova Scotia, Canada[END_REF][START_REF] Glasby | Surface composition and orientation interact to affect subtidal epibiota[END_REF][START_REF] Simpson | Large scale variability in the structure of sessile invertebrate assemblages in artificial habitats reveals the importance of local-scale processes[END_REF][START_REF] Smale | Extreme spatial variability in sessile assemblage development in subtidal habitats off southwest Australia (southeast Indian Ocean)[END_REF]. Such a spatial scale is consistent with the spatial scale of flow homogeneity obtained in simulations over GoL soft-bottom habitat [START_REF] Briton | High-resolution modelling of ocean circulation can reveal retention spots important for biodiversity conservation[END_REF]. This inventory methodology enabled us to test for the existence of structuring factors at the local and regional scales. Developed reef surface was calculated on the basis of technical specifications data present in the state concession documents taking into account only the colonizable surface (surfaces in contact with the sediment were excluded).

Colonization assessment and species selection

Assessment of ARs colonization was carried out by autonomous scuba-diving in 2020 by direct visual census counting the number of individuals of the five target species on the entire surface of ARs (7295 m²; Table1: Supplementary Material). On the ARs, all individuals older than one year (size > 2 cm) have been recorded. Among the species listed in previous ARs inventories in the GoL, we selected five species that were present in most inventories, easy to identify by direct visual census and spanning different phyla with contrasting life-history traits (Créocean, 2003(Créocean, & 2004; Table 1). We selected two gorgonians Eunicella singularis (Esper, 1791) and Leptogorgia sarmentosa (Esper, 1789), one bryozoan Pentapora fascialis (Pallas, 1766), one annelida Sabella spallanzanii (Gmelin, 1791) and one ascidian Halocynthia papillosa (Linnaeus, 1797) (Figure 2). The five species have a similar wide natural repartition area along European coasts ranging from 1 m to 250 m depth [START_REF] Giangrande | Utilization of the filter feeder polychaete Sabella[END_REF][START_REF] Gori | Spatial distribution patterns of the gorgonians Eunicella singularis, Paramuricea clavata, and Leptogorgia sarmentosa (Cape of Creus, Northwestern Mediterranean Sea)[END_REF][START_REF] Ponti | Mediterranean gorgonian forests: Distribution patterns and ecological roles. 3rd Mediterranean symposium on the conservation of coralligenous & other calcareous bioconcretions (SPA/RAC)At[END_REF][START_REF] Turner | The early development of epifaunal assemblages on artificial substrata at two intertidal sites on an exposed rocky shore in St. Andrews Bay, N.E. Scotland[END_REF][START_REF] Weinberg | The Life Cycle of a Gorgonian: Eunicella Singularis (Esper, 1794)[END_REF]. In addition, S. spallanzanii has been recorded along the coasts of Brazil, Australia and New Zealand where it is classified as an invasive non-indigenous species [START_REF] Currie | Reproduction and distribution of the invasive European fanworm Sabella spallanzanii (Polychaeta: Sabellidae) in Port Phillip Bay, Victoria, Australia[END_REF]. The five species are present in the rocky habitat of the NW Mediterranean Sea (Laubier, 1966;True, 1970;Hong, 1980). However, in the GoL, where natural rocky habitat covers uneven surfaces within the 6 geographical sectors (from 3,123 10 7 m 2 for the AGM sector to 5 10 5 m 2 for the LEU sector), the five species display different spatial distributions (Dutrieux et al., 2005;Dalias et al., 2011;[START_REF] Guizien | Database of sea fans spatial distribution in the northwestern Mediterranean sea[END_REF]S. Blouet personal observation). E. singularis is frequently observed and abundant throughout the GoL (from the AGM sector to the CST sector), while L. sarmentosa, less abundant, is present mainly in the center of the GoL (AGD, VLR, LEU, CST). P. fascialis is abundant in the west of the GoL (AGD, LEU, CST, and south of CST). The distribution of H. papillosa is not well known, however the species has been observed in all the rocky areas of the GoL. S. spallanzanii is present but rare in natural rocky habitat. Nevertheless, S. spallanzanii is very abundant in lagoons, ports and marinas of the GoL (S. Blouet personal observation) which has been indicated as preferred habitat of the species [START_REF] Currie | Reproduction and distribution of the invasive European fanworm Sabella spallanzanii (Polychaeta: Sabellidae) in Port Phillip Bay, Victoria, Australia[END_REF].

The five species display different life-history traits. All five species reproduce once a year in different seasons and with different strategies (Table 1).

S. spallanzani reproduces in January-February, when water temperature is the coldest. The species displays multiple reproductive strategies: internal fertilization, with larvae brooded either inside or outside the mineral tube secreted around the body, and external fertilization broadcast spawning [START_REF] Giangrande | Gametogenesis and larval development in Sabella spallanzanii (Polychaeta: Sabellidae) from the Mediterranean Sea[END_REF]. In addition, asexual reproduction by fission has been reported [START_REF] Read | Arrival of the alien fanworm Sabella spallanzanii (Gmelin, 1791) (Polychaeta: Sabellidae) in two New Zealand harbours[END_REF]. S. spallanzani releases lecithotrophic larvae with a planktonic larval duration (PLD) of about 4 weeks [START_REF] Giangrande | Gametogenesis and larval development in Sabella spallanzanii (Polychaeta: Sabellidae) from the Mediterranean Sea[END_REF]. Its life span can exceed 5 years, with sexual maturity after one year (Giangrande and Petraroli 1994;[START_REF] Giangrande | Gametogenesis and larval development in Sabella spallanzanii (Polychaeta: Sabellidae) from the Mediterranean Sea[END_REF].

Like most gorgonians, E. singularis releases lecithotrophic larvae in early summer (June to August). Even though larval competency period can reach up to 2 months [START_REF] Guizien | Survive or swim: different relationships between migration potential and larval size in three sympatric Mediterranean octocorals[END_REF][START_REF] Zelli | Settlement dynamics and recruitment responses of Mediterranean gorgonians larvae to different crustose coralline algae species[END_REF], PLDs ranging from 7 to 14 days best explained gene flow among E. singularis natural populations dwelling in the fragmented rocky habitat of the GoL [START_REF] Padron | Passive larval transport explains recent gene flow in a Mediterranean gorgonian[END_REF]. E. singularis life span can reach 25-30 years with sexual maturity before 6 years [START_REF] Gori | Spatial variability in reproductive cycle of the gorgonians Paramuricea clavata and Eunicella singularis (Anthozoa, Octocorallia) in the Western Mediterranean Sea[END_REF][START_REF] Weinberg | The Life Cycle of a Gorgonian: Eunicella Singularis (Esper, 1794)[END_REF].

The other gorgonian, L. sarmentosa also releases lecithotrophic larvae but in the late summer (September to October) and the PLD is unknown [START_REF] Rossi | Reproductive features and gonad development cycle of the soft bottom-gravel gorgonian Leptogorgia sarmentosa (Esper, 1791) in the NW Mediterranean Sea[END_REF]. L. sarmentosa life span can reach 20 years with female sexual maturity within 2-3 years after settlement [START_REF] Rossi | Reproductive features and gonad development cycle of the soft bottom-gravel gorgonian Leptogorgia sarmentosa (Esper, 1791) in the NW Mediterranean Sea[END_REF].

H. papillosa is a simultaneous hermaphrodite which releases larvae in late summer (September to October; [START_REF] Becerro | Reproductive cycles of the ascidians Microcosmus sabatieri and Halocynthia papillosa in the Northwestern Mediterranean[END_REF], presumably lecithotrophic. The PLD of H. papillosa larvae is unknown but PLD shorter than 12 hours has been consistently reported for other solitary ascidian species [START_REF] Ayre | Genetic evidence for contrasting patterns of dispersal in solitary and colonial ascidians[END_REF]. We did not find any data about the age at sexual maturity and the life span of H. papillosa. However, the ascidians are considered as highly invasive, particularly because of their rapid growth and early sexual maturity [START_REF] Zhan | Ascidians as models for studying invasion success[END_REF], with some species such as Ciona intestionalis complex and Ciona savigniy, reaching sexual maturity at the age of 2 months [START_REF] Zhan | Ascidians as models for studying invasion success[END_REF] and continuous spawning (Carver et al., 2003).

P. fascialis displays both sexual and asexual reproduction. During sexual reproduction, most bryozoans release lecithotrophic larvae which settle after a few minutes or a few hours, rarely beyond several days [START_REF] Keough | Patterns of recruitment of sessile invertebrates in two subtidal habitats[END_REF]. P. fascialis larval release has been inferred to happen in June based on recruitment observations (Cocito et al., 1998a). Asexual reproduction happens by colony fragmentation or budding extension (Cocito et al., 1998b). Individual life span is estimated to be about 10 years with early sexual maturity after 2 years (Cocito et al., 1998b).

Statistical analysis

We examined to which extent ARs colonization measured by the presence/absence of the five species and on the dissimilarity between their co-occurrence assemblages is affected by the AR shape, location and timing of deployment (age). Due to the GoL AR deployment set up, the effect of some factors was tested only locally. We first verified that shape was not affecting colonization, by testing the effect of shape in one geographical sector where different shapes were deployed at same depth and time (VLR1/VLR2, Table 1: Supplementary Material). Second, the effect of factors age (2 levels: 1985 and the 2000s) and depth (2 levels, >20 m, <20 m depth) was tested within the AGD geographical sector only. Factor age was tested in 3 sites at <20 m (AGD1/ AGD2 & AGD3) and factor depth was tested in 4 sites deployed in 2009 (AGD2&AGD3 / AGD4&AGD5). Third, factor location was tested both at regional and local levels (geographical sector being the regional factor, site being the local factor). To avoid any confounding effect due to age or depth, the regional factor vs the local factor were tested on the 5 geographical sectors where AR were deployed during the 2000s immersion phase and at <20 m depth only (AGM, AGD, VLR, GRU, LEU: 5x2 levels).

A Jaccard similarity matrix was built on presence/absence data across all pairwise sampling units used in each test (Table 2: supplementary material). Two multivariate analyses were performed. Non-parametric multivariate analysis of variance with permutation was applied to test for the effects of shape, age, depth and geographical sector on species assemblages (NP-manova: [START_REF] Anderson | A new method for non-parametric multivariate analysis of variance[END_REF][START_REF] Zar | Biostatistical analysis[END_REF]. Site was considered as a random factor, nested either in depth, age or the geographical sector. Another multivariate analysis was performed to cluster most similar species assemblages in the sector of AGD (SIMPROF: [START_REF] Clarke | Testing of null hypotheses in exploratory community analyses: similarity profiles and biota-environment linkage[END_REF] (Table 2: supplementary material).

When significant differences between the five species co-occurrence assemblages were detected for a factor, a non-parametric univariate analysis (Kruskall-Wallis) was performed for each species independently to detect the species driving the difference (Table 2: supplementary material). Fisher posthoc test was used to identify the site where the difference arose. A same p-value of 0.05 was taken for detecting significant differences. Analyses were performed with Matlab software using the Fathom package for multivariate analyses (Jones, 2014) and the Matlab statistics toolbox for univariate analyses.

Results

Artificial reefs colonization by the five target species at regional scale Among the five target species, S. spallanzanii was the only one whose presence was recorded in all the sampling units and geographical sectors (Figure 3). In only one out of the 16 sites, it was the only species detected. H. papillosa and L. sarmentosa were detected in five of the six geographical sectors (not present in GRU and AGM respectively). P. fascialis and E. singularis were detected in 3 of the 6 geographical sectors (CST, LEU, and AGD for P. fascialis and AGM, AGD and LEU for E. singularis). Finally, E. singularis was the least frequently observed species, being detected in only 5 out of the 16 sites. In all geographical sectors, at least three of the five target species were detected, except GRU where only two of the five species were detected. Assemblages of two species were found in two sites out of the 16 (sectors AGD and GRU), assemblages of three species were found in 7 sites, assemblages of four species were found in 2 sites and assemblages of five species were found in 3 sites.

Shape effect on five target species at local scale

No difference in species assemblages was found in VLR between pipe and steel cage (NP-MANOVA, F=9.08 e-19, P>0.05; Table 2).

Similarly, in the geographical sector GRU, the same assemblages were detected on sites differing by ARs shape (pipes and a heap of telegraphical poles) in the same location (GRU1 and GRU2). Conversely, different assemblages were detected between two sites (GRU2 and GRU3) separated by 9 km although ARs shape was the same (heaps of telegraphical poles). Due to the absence of replication in unit samplings, it was not possible to perform a statistical test on the effect of shape in this geographical sector. Age and depth effects on five target species at local scale Despite all five target species being detected on ARs in the AGD sector, assemblage composition among sites differed (Figure 3). Clustering of sampling units within the 5 sites (AGD1, AGD2, AGD3, AGD4, AGD5) in AGD identified 2 clusters (SIMPROF: P<0.05; Figure 4). The two sites (AGD2 and AGD3) with same age (2009), depth range (less than 20 m) and reef shape (pipe) were attributed to different clusters.

In fact, one cluster grouped ARs of different ages at a same depth (1985 in site AGD1 and 2009 in site AGD2) while the other cluster grouped ARs of the same age but at different depths (less than 20 m in site AGD3 and more than 20 m in sites AGD4 and AGD5). In both clusters, different AR shapes were found (steel cage and pipes in one cluster, pipes, Comin and Bonna in the other cluster) (Figure 4). The geographic distance between the two clusters (AGD4-AGD5) and (AGD3-AGD4-AGD5) was 7.5 km and the median value of the intra-cluster geographic distance was 3 km.

Multivariate analysis of variance confirmed that neither age (NP-MANOVA, F=1.43 P<0.05; Table 3) nor depth (NP-MANOVA, F=1.37 P<0.05; Table 4) explained site differences in the five species assemblages found on ARs in AGD (P<0.05; Tables 3 and4). Differences among the five sites in AGD were due to different ARs colonization by three species (Kruskall-Wallis: E. singularis, P. fascialis, H. papillosa; all P<0.05; Table 5). Site AGD3 differed from sites AGD5, AGD1 and AGD2 by the presence of E. singularis, sites AGD4 and AGD5 differed from sites AGD1 and AGD2 due to the presence of P. fascialis and the site AGD1 differed from sites AGD5, AGD3 and AGD4 by the presence of H. papillosa (Post-hoc tests, Table 3: supplementary material).

Geographical effect on five target species at local and regional scales

The five species co-occurrence assemblages on ARs deployed in the same period and at the same depth were significantly different at both regional and local scales (NP-MANOVA): geographical sector F=5.09 P<0.05; site (geographical sector) F=2.78 P<0.05; Table 6). These differences were due to different colonization of ARs by three of the five species, E. singularis, L. sarmentosa and P. fascialis. For the latter two species, regional differences (Kruskall-Wallis: geographical sector P=0.0002) were more significant than local differences (Kruskall-Wallis: site P=0.001 for L. sarmentosa and P=0.005 for P. fascialis; Table 7). Both species were not detected in the north of the GoL (AGM). In contrast, for E. singularis, local differences (Kruskall-Wallis: site P=0.017) were more significant than regional ones (geographical sector P=0.036; Table 7), the species being detected in geographical sectors in the north, center and south of GoL.

Discussion

The resilience of marine natural biodiversity, including benthic invertebrates, is essential to ensure the sustainability of ecosystem functions. Increasing habitat surface is key to support the natural biodiversity resilience, which in the case of the rocky habitat can be achieved by effectively integrated ARs. The effective integration of ARs into the rocky habitat network starts with their colonization by species from the natural habitat. In the present study the five species selected from the natural rocky habitat of the GoL colonized differently the ARs spread along the GoL coastline ten years after their deployment. Locally, neither age, immersion, depth or shape of the ARs significantly affected colonization patterns.

Colonization of ARs are expected to evolve toward a stable state comparable to that of the natural environment, through the succession of opportunistic species (wide dispersal, high fertility, low tolerance of reduced resource levels, short life-spans, minimal dietary specialisation) followed by specialized species (limited dispersal, slow growth to a large size at maturity, delayed and limited reproduction, optimization to reduced resources and long life spans; [START_REF] Platt | Natural Disturbances and Directional Replacement of Species[END_REF]Faurie et al., 2003). Monitoring of ARs short-term colonization (<3 years) have indeed shown a dominance of pioneer species (hydroids, serpulids, barnacles and bivalves), most of them having life history traits typical of opportunistic species [START_REF] Fariñas-Franco | Community convergence and recruitment of keystone species as performance indicators of artificial reefs[END_REF][START_REF] Glasby | Interactive effects of shading and proximity to the seafloor on the development of subtidal epibiotic assemblages[END_REF][START_REF] Moura | Effect of depth and reef structure on early macrobenthic communities of the Algarve artificial reefs (southern Portugal)[END_REF]Pamintuan and Ali, 1994;[START_REF] Ponti | Benthic assemblages on artificial reefs in the northwestern Adriatic Sea: Does structure type and age matter?[END_REF]Relini et al., 1994;[START_REF] Spagnolo | Patterns of colonization and succession of benthic assemblages in two artificial substrates[END_REF][START_REF] Toledo | Ecological succession of benthic organisms on niche-type artificial reefs[END_REF].

Long-term studies confirmed successions in ARs colonization [START_REF] Burt | Benthic development on large-scale engineered reefs: A comparison of communities among breakwaters of different age and natural reefs[END_REF][START_REF] Nicoletti | Long-term changes in a benthic assemblage associated with artificial reefs[END_REF][START_REF] Perkol-Finkel | Recruitment of benthic organisms onto a planned artificial reef: shifts in community structure one decade post-deployment[END_REF][START_REF] Whomersley | Long-term dynamics of fouling communities found on offshore installations in the North Sea[END_REF], but none have described saturation [START_REF] Svane | On the Problems of Epibioses, Fouling and Artificial Reefs, a Review[END_REF]. In the Tyrrhenian Sea, Nicolleti et al. ( 2007) described colonization in 5 distinct phases: (i) a first recruitment by pioneer species (hydroids, serpulids, barnacles and bivalves) during the firsts months after immersion, followed by phases of (ii) cover dominance, (iii) regression and (iv) absence of Mytilus galloprovincialis for more than 10 years. The installation of diverse bio-builders bryozoans characteristic of the natural environment was recorded after 20 years only (v). Our study shows that biobuilder engineering species such as bryozoans (P. fascialis) and gorgonians (E. singularis, L. sarmentosa) colonized ARs as early as 10 years after their deployment, without significant difference between 10 years and 35 years old ARs. However, the presence of S. spallanzanii described as an opportunistic species (sexual precocity, various reproduction modes, rapid growth, short lived; [START_REF] Giangrande | Gametogenesis and larval development in Sabella spallanzanii (Polychaeta: Sabellidae) from the Mediterranean Sea[END_REF] on all ARs independently of their age of deployment suggests that ARs did not yet reach a stable state comparable to the natural environment. Thus, the presence of bio-builders is not a sufficient indicator of the ARs naturalization to the local biodiversity. The GoL's ARs being located in the sandy coastal zone are likely regularly disturbed by sediment deposits due to river delivery or/and their resuspension by either trawling activities or the mechanical action of the swell [START_REF] Dufois | Impact of winter storms on sediment erosion in the Rhone River prodelta and fate of sediment in the Gulf of Lions (North Western Mediterranean Sea)[END_REF][START_REF] Durrieu De Madron | Trawling-induced resuspension and dispersal of muddy sediments and dissolved elements in the Gulf of Lion (NW Mediterranean)[END_REF][START_REF] Ulses | Suspended sediment transport in the Gulf of Lions (NW Mediterranean): Impact of extreme storms and floods[END_REF]. Testing the impact of swell and sediment deposit on ARs requires exploring the colonization of ARs along a gradient of depth and distance from the coast ( [START_REF] Van Der Stap | Marine Fouling Assemblages on Offshore Gas Platforms in the Southern North Sea: Effects of Depth and Distance from Shore on Biodiversity[END_REF]. However, current ARs deployment in the GoL ranged from 15 to 30 m depth and within 3 miles from the coast and did not allow testing for differential effect of sediment disturbances as swell impact occurs every year in this area [START_REF] Guizien | Spatial variability of wave conditions in the Gulf of Lions (NW Mediterranean Sea)[END_REF]. Testing the impact of sediment disturbance on ARs colonization would require exploring reefs deployed deeper than 50 m, such as the anchorages of the floating wind farm that will be placed in the GoL in the next future (https://info-efgl.fr/le-projet/le-parc/#). Light is also expected to be an important factor structuring benthic assemblages, along a depth gradient in natural and artificial environments (Glasby, 1999a;Glasby, 1999b;[START_REF] Svane | On the Problems of Epibioses, Fouling and Artificial Reefs, a Review[END_REF] . Absence of depth effect in the present study, although in the GoL light intensity strongly attenuates within the upper 30 m of the water column (Durrieu de Madron et al., 2011), was potentially a bias due to the five species selected in the present study whose distributions are not strongly structured by light intensity.

Another factor which has been shown to drive the intensity of ARs colonization is structural complexity (see [START_REF] Bohnsack | Artificial reef research: A review with recommendations for future priorities[END_REF] for a review). Nevertheless, there is no consensus about the relationship between complexity and subtidal benthic invertebrates abundance, due to potential bias in controlling the surface and scale in ARs of different complexity (Rouse et al., 2019;[START_REF] Strain | Eco-engineering urban infrastructure for marine and coastal biodiversity: Which interventions have the greatest ecological benefit?[END_REF]. The similarity in the 5 species co-occurrence between different reef shapes at the same depth and of the same age found in the present study suggests that shape is less important than the geographical location in AR colonization by benthic invertebrates when controlling the colonized surface. However, shape is an imprecise measure of structural complexity. The latter is rarely assessed quantitatively and can be described by different metrics which may be similar for apparently different shapes (such as steel cages and pipes; [START_REF] Riera | Vers une construction raisonnée d'une nouvelle génération de récifs artificiels : Analyses comparatives des facteurs intrinsèques favorisant leur colonisation de la micro à la macro-échelle 192[END_REF].

Benthic invertebrate assemblages result from complex processes that operate at multiple spatial and temporal scales [START_REF] Smale | Extreme spatial variability in sessile assemblage development in subtidal habitats off southwest Australia (southeast Indian Ocean)[END_REF]. At the regional scale, larval availability can become a major factor explaining colonization success [START_REF] Padron | Passive larval transport explains recent gene flow in a Mediterranean gorgonian[END_REF]. Change in the composition of assemblages during the early colonization of artificial substrates by benthic invertebrates has been attributed to the availability and abundance of larvae during the seasons rather than a sequence of distinct succession [START_REF] Bramanti | Settlement and recruitment: The first stages in the life cycle of two epibenthic suspension feeders (Corallium rubrum and Anomia ephippium)[END_REF][START_REF] Turner | The early development of epifaunal assemblages on artificial substrata at two intertidal sites on an exposed rocky shore in St. Andrews Bay, N.E. Scotland[END_REF]. The larval behaviour (buoyancy and motility) and the characteristics of the biological cycle of the species (spawning timing and PLD) can play a key role in determining the dispersal distance [START_REF] Todd | Larval supply and recruitment of benthic invertebrates: do larvae always disperse as much as we believe?[END_REF], and consequently the possibility to reach habitat suitable for settlement. Dispersion distance is generally correlated with PLD, thus a species with a long PLD is supposed to colonize habitats further away than species with a shorter PLD [START_REF] Shanks | Barnacle settlement versus recruitment as indicators of larval delivery. II. Time-series analysis and hypothesized delivery mechanisms[END_REF]. In this study, the five species were chosen among different phyla known for their contrasting planktonic durations, swimming abilities and larval dispersal periods, although these larval traits are only known accurately for E. singularis [START_REF] Guizien | Survive or swim: different relationships between migration potential and larval size in three sympatric Mediterranean octocorals[END_REF][START_REF] Zelli | Settlement dynamics and recruitment responses of Mediterranean gorgonians larvae to different crustose coralline algae species[END_REF]. P. fascialis and H. papillosa, the two species with short PLD (<24h and <12h, respectively) colonized ARs located close to the natural habitat where they are present (<4.8 km and <10 km, respectively). The coastal circulation of the GoL allows such dispersal distance over periods of a few days [START_REF] Guizien | Using larval dispersal simulations for marine protected area design: Application to the Gulf of Lions (northwest Mediterranean)[END_REF]. S. spallanzanii, which has a PLD larger than 3 weeks, colonized all the inventoried ARs, in line with a dispersal distance of 40 km after 3 weeks [START_REF] Guizien | Using larval dispersal simulations for marine protected area design: Application to the Gulf of Lions (northwest Mediterranean)[END_REF]. In contrast, E. singularis did not colonize all ARs within geographical sectors of 30 km width where the species is present in the natural habitat, although a 2 week PLD was expected to enable such dispersal [START_REF] Padron | Passive larval transport explains recent gene flow in a Mediterranean gorgonian[END_REF]. The other gorgonian species, L. sarmentosa colonized more ARs located within distances of less than 30 km from its natural habitat than E. singularis while the PLD of the two species are presumably the same. This suggests that other factors influence the colonization failure of E. singularis.

Reproductive traits are another key to the success in colonising new settings [START_REF] Stearns | Life history evolution: successes, limitations, and prospects[END_REF]. In this regard, E. singularis colonization potential could be limited by its low fecundity (~25-40 larvae.cm -1 of colony branch; [START_REF] Ribes | Cycle of gonadal development in Eunicella singularis (Cnidaria: Octocorallia): trends in sexual reproduction in gorgonians: Gametogenic cycle in a temperate gorgonian[END_REF]Theodor, 1967) compared to that of L. sarmentosa (~75 larvae.cm -1 of colony branch; [START_REF] Rossi | Net negative growth detected in a population of Leptogorgia sarmentosa: quantifying the biomass loss in a benthic soft bottom-gravel gorgonian[END_REF][START_REF] Rossi | Reproductive features and gonad development cycle of the soft bottom-gravel gorgonian Leptogorgia sarmentosa (Esper, 1791) in the NW Mediterranean Sea[END_REF]. The wide colonization of ARs by S. spallanzanii is in line with its reproductive traits typical of opportunistic species (early sexual maturity, high fecundity with more than 50 000 eggs per female, [START_REF] Currie | Reproduction and distribution of the invasive European fanworm Sabella spallanzanii (Polychaeta: Sabellidae) in Port Phillip Bay, Victoria, Australia[END_REF], a fertilization close to 100%, [START_REF] Giangrande | Gametogenesis and larval development in Sabella spallanzanii (Polychaeta: Sabellidae) from the Mediterranean Sea[END_REF]. Since arriving in the Pacific Ocean, S. spallanzanii has been declared one of the ten priority pest species in the marine environment by the Australian authorities and classified as an invasive species [START_REF] Hayes | Historical Overfishing and the Recent Collapse of Coastal Ecosystems[END_REF]. Similarly to S. spallanzanii, H. papillosa colonized nearly all ARs located within its 10 km dispersal distance from the natural habitat. Within the ascidian class, a wide disparity in species fecundity has been reported [START_REF] Pandian | Reproduction and Development in Echinodermata and Prochordata[END_REF]. This suggests H. papillosa reproductive traits would be close to those of invasive ascidians [START_REF] Zhan | Ascidians as models for studying invasion success[END_REF].

Ultimately, understanding ARs colonization requires a precise mapping of source populations in the natural environment. To this respect, the abundance of S. spallanzani on ARs is surprising, as the species is not abundant in the natural rocky habitat of the GoL. For this species, other sources of larval supply than natural settings should be considered, such as the numerous ports and marinas along the coast of the GoL, as S. spallanzanii is very tolerant to environmental conditions [START_REF] Currie | Reproduction and distribution of the invasive European fanworm Sabella spallanzanii (Polychaeta: Sabellidae) in Port Phillip Bay, Victoria, Australia[END_REF]. In this case of intense colonization by an endemic benthic invertebrate species, ARs apparently extended its metapopulation, acting as stepping stones for further larval dispersal beyond its natural current colonization limits [START_REF] Bishop | Effects of ocean sprawl on ecological connectivity: impacts and solutions[END_REF][START_REF] Wang | Global warming and artificial shorelines reshape seashore biogeography[END_REF].

In the GoL, the 14 500 m 3 of ARs deployed 30 years ago are now decommissioned and the relevance of their removal is currently debated. Connectivity between natural populations has been shown to support species resilience after disturbances in fragmented habitat, and could be extended to ARs [START_REF] Fahrig | Effects of Habitat Fragmentation on Biodiversity[END_REF]. However, ARs may also facilitate the spread of non-indigenous species introduced with maritime traffic in ports [START_REF] Glasby | Nonindigenous biota on artificial structures: Could habitat creation facilitate biological invasions?[END_REF].

The present study advocates accounting for the geographical arrangement in planning ARs deployment to extend the habitat of hard bottom benthic invertebrate natural populations.

Figure 1 :

 1 Figure 1: Map showing the layout of the 6 geographical sectors (AGM, AGD, VLR, GRU, LEU, CST) and the 16 sites prospected (indicated in dark blue) for artificial reef sampling. Areas indicated in light blue corresponded to artificial reef deployment areas. Areas indicated in black corresponded to natural hardsubstrate areas. Some of the AR shape illustrations were taken from Tessier et al., 2014.

Figure 2 :

 2 Figure 2: Photographs of the five species inventoried on AR a) Sabella spallanzanii, b) Eunicella singularis, c) Leptogorgia sarmentosa, d) Halocynthia papillosa and e) Pentapora fascialis; all © Blouet sylvain

Figure 3 :

 3 Figure 3: Map showing the five species co-occurence assemblages inventoried on ARs in the 16 sites (AGM1, AGM2, AGD1, AGD2, AGD3, AGD4, AGD5, VLR1, VLR2, GRU1, GRU2, GRU3, LEU1, LEU2, CST1, CST2) in the Gulf of Lion.

Figure 4 :

 4 Figure 4: Diversity and dissimilarity. Upper part: diversity of species assemblages in the sampling units of the five sites (AGD1, AGD2, AGD3, AGD4, AGD5) in AGD sector together with the type of AR, depth and years of deployment. Lower part: dendrogram obtained by group average clustering based on the Jaccard dissimilarity index using the presence /absence of species (P=0.04 at 62% of dissimilarity). The red dotted line separates the two clusters identified by the analysis.

Table 1 :

 1 Life history traits (spawning period, relative fecundity, life expectancy, age at sexual maturity, larval type, planktonic larval duration (PLD), relative abundance in the GoL) for the five species inventoried on ARs.

	Species	Spawning	Fecundity	Life	Age at	Larval type	Planktonic	Abundance in
		period		expectancy	sexual		larval duration	the GoL
					maturity		(PLD)
	Leptogorgia	June -August	**	20 years	2-3 years	Lecithotrophic	1-4 weeks	**
	sarmentosa						(estimated
							from other
							gorgonians)
	Eunicella	July -August	*	25-30 years	6 years	Lecithotrophic	1-4 weeks (in	***
	singularis						aquarium)
	Pentapora	June	**	10 years	2 years	Lecithotrophic	<1 days	**
	fascialis
	Sabella	January -	***	5 years	1 years	Lecithotrophic	21 days (in	*
	spallanzani	February					aquarium)
	Halocynthia	September -	***	unknow	2 months	lecithotrophic	<12 hours	**
	papillosa	October	(estimated		(estimated	(estimated	(estimated
			from other		from other	from other	from other
			ascidians)		ascidians)	ascidians)	ascidians)
			Symbol legend: * = low, **=medium ***=strong

Table 2 :

 2 Results of multivariate NP-Manova testing the effect of the shape between steel cage and pipe on the presence/absence assemblage. Sites tested: VLR1 vs VLR2. Significant (P<0.05) values in bold.

	Multivariate measures

Table 3 :

 3 Results of multivariate NP-Manova testing the interactive effects of the age (of deployment) and site (nested in Age) on the presence/absence assemblage. Sites tested: AGD1 vs (AGD2 vs AGD3).Significant (P<0.05) values in bold.

	Multivariate measures

Table 4 :

 4 Results of multivariate NP-Manova testing the interactive effects of the depth (of deployment) and site (nested in depth) on the presence/absence assemblage. Sites tested: (AGD2 vs AGD3) vs (AGD4 vs AGD5). Significant (P<0.05) values in bold.

	Multivariate measures

Table5: Results of univariate Kruskall-Wallis testing the effect of location of deployment on the presence/absence by species. Sites tested: all sites in AGD. Significant (P<0.05) values in bold.

Table 6 :

 6 Results of multivariate NP-Manova testing the interactive effects of geographical sector and site (nested in geographical sector) on the presence/absence assemblage.

	Sites tested: (AGM1 vs AGM2)

Table 7 :

 7 Results of univariate Kruskall-Wallis testing the effects of geographical sector and site on the presence/absence by species. Sites tested: AGM1 vs AGM2 vs AGD2 vs AGD3 vs VLR1 vs VLR2 vs GRU1 vs GRU2 vs GRU3 vs LEU1 vs LEU2. Significant (P<0.05) values in bold.

	Univariate measures

Peer Community Journal, Vol. 2 (2022), article e24 https://doi.org/10.24072/pcjournal.107

Acknowledgements

This work was funded by the Agence de l'Eau Rhône-Méditerranée-Corse under project ICONE -Impacts actuels et potentiels de la CONnectivité Ecologique ajoutée par les récifs artificiels sur la biodiversité fixée des substrats durs du Golfe du Lion (PI, K. Guizien, AAP 2016). The authors gratefully acknowledge the helpful assistance during sampling of the staff of the Aire Marine Protégée Agathoise. A preprint version of this article has been reviewed and recommended by Peer Community In Ecology (https://doi.org/10.24072/pci.ecology.100093)

Data accessibility

Data are available online: https://doi.org/10.5281/zenodo.5567023

Supplementary material

Script and codes are available online: https://doi.org/10.1101/2021.10.08.463669

Supplementary

Authors contribution

SB and KG conceived the study, SB carried out sampling and statistical analysis. All contributed to manuscript writing.

Conflict of interest disclosure

The authors of this article declare that they have no financial conflict of interest with the content of this article.