2022)], we numerically model settling clouds produced by localised instantaneous releases of heavy particles in a quiescent fluid. By modelling particles as a field of mass concentration in an equilibrium Eulerian approach, our two-way coupling simulations recover our original experimental observation of a maximum growth rate for clouds laden with particles of finite Rouse number R ≈ 0.22. Consistent with the literature on buoyant vortex rings, our clouds verify the relation α ∝ Γ -2 ∞ between the clouds' growth rate α and their eventually constant circulation Γ∞. As the Rouse number approaches R ≈ 0.22, the baroclinic forcing of the circulation reduces down to a minimum, thus optimising the cloud growth rate α. This analysis highlights the role of the mean flow in the enhanced entrainment of ambient fluid by negatively buoyant clouds. Our results also validate, on the basis of direct comparison with experimental results, the use of a one-fluid two-way coupling numerical model to simulate particle clouds in the limit of weak particle inertia.

I. INTRODUCTION

The present study follows up on a series of 514 systematic laboratory experiments presented in Kriaa et al. [1] which focus on the evolution of instantaneously-released particle-laden clouds settling from rest in initially quiescent water, under the sole action of their buoyancy. By varying the size of particles yet keeping an identical buoyancy for all clouds, the latter proved to grow linearly in depth with a growth rate α that reaches a maximum for a finite particle inertia R ≃ 0.3 ± 0.1, where the Rouse number R = w s /U ref is the ratio of the terminal velocity of an isolated particle w s over the reference fall velocity U ref of the cloud due to the sole action of its buoyancy. This new observation was unexpected because the theory of turbulent thermals [2], commonly used to model such clouds, predicts that all clouds should grow similarly as a salt-water cloud of identical buoyancy. Yet, our measurements revealed that the particle inertia can increase the growth rate by up to 75%. Our experiments did not enable us to answer some remaining questions: Through what mechanism do particle clouds entrain more than salt-water clouds of identical buoyancy? In particular, does the particle inertia alter the mean flow around the cloud or does it modulate the intensity of the turbulence that develops inside the cloud? The aim of the present study is to gain understanding on these questions thanks to complementary Direct Numerical Simulations (DNS). The numerical approach to be adopted should capture the macroscopic physics of particle clouds as observed in our experiments, yet with the minimum ingredients to keep a low numerical cost, with a motivation to later apply the method to a planetary-scale multiphase flow called 'iron snow' [START_REF] Rückriemen | The Fe snow regime in Ganymede's core: A deep-seated dynamo below a stable snow zone[END_REF].

When the size of particles, their volume fraction and the particle-to-fluid density anomaly are low, the fluid governs the motion of particles whose feedback on the flow is negligible, a situation called one-way coupling [START_REF] Balachandar | Turbulent Dispersed Multiphase Flow[END_REF]. Many studies have considered the one-way coupling between particles and fluid motions (e.g. [START_REF] Brandt | Particle-Laden Turbulence: Progress and Perspectives[END_REF][START_REF] Good | Settling regimes of inertial particles in isotropic turbulence[END_REF][START_REF] Yoshimoto | Self-similar clustering of inertial particles in homogeneous turbulence[END_REF][START_REF] Goto | Sweep-Stick Mechanism of Heavy Particle Clustering in Fluid Turbulence[END_REF][START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF]) and evidenced that inertia is a source of non-uniformities in the field of particle concentration [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF]: despite the incompressibility of the fluid phase, inertia allows the disperse phase to be compressible. M. R. Maxey [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF] showed that very small particles behave as passive tracers due to their low inertia; conversely, large particles are insensitive to local modifications of the flow due to their large response time (see also [START_REF] Yoshimoto | Self-similar clustering of inertial particles in homogeneous turbulence[END_REF][START_REF] Ghosh | How turbulence enhances coalescence of settling particles with applications to rain in clouds[END_REF]). Due to this low-pass 'inertial filtering' of timescales, dense particles of intermediate size partly decouple from fluid motions. For some finite inertia they have been observed to optimally couple with the surrounding flow and concentrate in some specific regions. This phenomenon of 'preferential concentration' ( [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF], see also [START_REF] Yoshimoto | Self-similar clustering of inertial particles in homogeneous turbulence[END_REF][START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF][START_REF] Aliseda | Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence[END_REF][START_REF] Salazar | Experimental and numerical investigation of inertial particle clustering in isotropic turbulence[END_REF][START_REF] Toschi | Lagrangian Properties of Particles in Turbulence[END_REF]) leads to the accumulation of small particles in regions of large strain rate or equivalently of low vorticity [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF], possibly due to centrifuging of particles outside of vortices [START_REF] Brandt | Particle-Laden Turbulence: Progress and Perspectives[END_REF][START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF], while larger dense particles rather concentrate in regions of vanishing fluid acceleration [START_REF] Goto | Sweep-Stick Mechanism of Heavy Particle Clustering in Fluid Turbulence[END_REF][START_REF] Bec | Acceleration statistics of heavy particles in turbulence[END_REF].

With the addition of gravity, another heterogeneity originates from the tendency of heavy particles to fall on the side of downward velocity of eddies [START_REF] Wang | Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence[END_REF], which can be grasped by simple advective arguments even in laminar flows (e.g. [START_REF] Nielsen | Turbulence effects on the settling of suspended particles[END_REF]). This phenomenon of 'fast tracking' or 'preferential sweeping' is also optimum for a finite particle inertia [START_REF] Wang | Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence[END_REF]; it has been observed in 2D periodic laminar vortices [START_REF] Maxey | Gravitational Settling of Aerosol Particles in Randomly Oriented Cellular Flow Fields[END_REF] and in turbulent flows (e.g., [START_REF] Good | Settling regimes of inertial particles in isotropic turbulence[END_REF][START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF][START_REF] Hassaini | Scale-to-scale turbulence modification by small settling particles[END_REF][START_REF] Elghobashi | On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification[END_REF][START_REF] Bosse | Small particles in homogeneous turbulence: Settling velocity enhancement by two-way coupling[END_REF]). Additionally, due to their inertia, falling particles lag behind fluid motions so their trajectories are biased in the direction of gravity. As particles cross upwelling and downwelling regions, they spend more time in upwelling regions, increasing their sensitivity to velocity variations in these upwellings. This 'loitering' [START_REF] Nielsen | Turbulence effects on the settling of suspended particles[END_REF] is a third example of preferential sampling of the flow.

All these phenomena of preferential sampling determine the distribution of particles in time and space. This observation is paramount if the flow is altered by the feedback of particles on the fluid, a situation referred to as twoway coupling [START_REF] Balachandar | Turbulent Dispersed Multiphase Flow[END_REF]. Inclusion of this feedback in numerical simulations has proved to be essential as it further modifies the structure of particle-laden flows even in well-controlled idealised isotropic turbulent flows [START_REF] Hassaini | Scale-to-scale turbulence modification by small settling particles[END_REF][START_REF] Elghobashi | On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification[END_REF][START_REF] Bosse | Small particles in homogeneous turbulence: Settling velocity enhancement by two-way coupling[END_REF][START_REF] Monchaux | Settling velocity and preferential concentration of heavy particles under two-way coupling effects in homogeneous turbulence[END_REF], especially in the presence of gravity which breaks the flow isotropy. As an example, experiments [START_REF] Aliseda | Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence[END_REF] and simulations [START_REF] Bosse | Small particles in homogeneous turbulence: Settling velocity enhancement by two-way coupling[END_REF][START_REF] Monchaux | Settling velocity and preferential concentration of heavy particles under two-way coupling effects in homogeneous turbulence[END_REF] have shown that the non-linear modification of the settling velocity due to preferential concentration and preferential sweeping is further favoured by two-way coupling since clusters of particles drag fluid with them as they sweep downward, producing downward acceleration of the fluid which enhances the fluid velocity and therefore the particle settling velocity. Most importantly, even in the presence of low mass loadings, the feedback of particles on the fluid is essential if the flow is driven by the particles themselves, e.g. in downdrafts [START_REF] Kruger | Dynamics of Downdraughts and Cold Pools: An Experimental and Numerical Study[END_REF], turbidity currents [START_REF] Ouillon | Turbidity currents propagating down a slope into a stratified saline ambient fluid[END_REF][START_REF] Necker | High-resolution simulations of particle-driven gravity currents[END_REF], and presumably iron snow [START_REF] Rückriemen | The Fe snow regime in Ganymede's core: A deep-seated dynamo below a stable snow zone[END_REF]. The present work fits in this framework: in our experiments [1] particles were released from rest in quiescent water, hence all fluid motions resulted from the drag exerted by particles on water during their fall.

Thus, in this study we adopt an equilibrium Eulerian approach [START_REF] Balachandar | Turbulent Dispersed Multiphase Flow[END_REF] to model our settling particle clouds at reasonably low numerical cost, while still accounting for the two essential ingredients of (i) the feedback of particles on the fluid through a drag term which forces the flow, and (ii) a differential motion between water and settling particles through a gravitational drift, a formalism already used in the literature to model particle-laden flows [START_REF] Fabregat Tomàs | Numerical simulations of rotating bubble plumes in stratified environments[END_REF][START_REF] Boffetta | The Eulerian description of dilute collisionless suspension[END_REF][START_REF] Reali | Layer formation in sedimentary fingering convection[END_REF][START_REF] Lemus | Modelling Settling-Driven Gravitational Instabilities at the Base of Volcanic Clouds Using the Lattice Boltzmann Method[END_REF][START_REF] Chou | Numerical study of particle-induced Rayleigh-Taylor instability: Effects of particle settling and entrainment[END_REF]. The latter effect is quantified by a Rouse number which is about R = 0.3 for the optimum growth rate in experiments, and which lies below unity for most of our clouds. Our experiments therefore fit in the range of validity of the equilibrium Eulerian model, as pointed out by Boffetta et al. [START_REF] Boffetta | The Eulerian description of dilute collisionless suspension[END_REF]. Note however that we extend our analysis up to R ≃ 3 as in our experiments. Boffetta et al. [START_REF] Boffetta | The Eulerian description of dilute collisionless suspension[END_REF] showed that such particles have so much inertia that they undergo the 'sling effect' i.e. particles tend to converge and thereby form caustics, so that the continuum modelling of particle motions through a unique velocity field locally breaks down. Yet, our results using the equilibrium Eulerian formalism up to R ≃ 3 show good agreement with our experiments, supporting the fact that if the sling effect plays any part in experiments, it still has no strong statistical signature on the macroscopic quantities we measured in experiments and reproduce numerically here.

The paper is organised as follows. Section II introduces the equations of motion to model particles as a field of mass concentration verifying an advection-diffusion equation and forcing the flow through a drag term. Section III presents the numerical setup of the three-dimensional numerical simulations. Then, section IV presents the two main regimes of cloud settling; it notably evidences the same maximum of the cloud growth rate α as observed in our experiments [1]. Section V is devoted to the analysis of this effect, showing that particles with a Rouse number close to the optimum impose a weaker baroclinic forcing of the clouds' circulation, resulting in an enhanced growth rate according to the theory of buoyant vortex rings. Further discussion and concluding remarks are presented in section VI. Appendix A provides details on the robustness of the numerical simulations, and appendix B gives numerical results for particle clouds of larger Reynolds number than those presented in the core of this study.

II. EQUATIONS OF MOTION

If a particle has finite inertia and is not neutrally buoyant, it moves with a velocity v p which is different from that of the fluid v in its vicinity, and the particle acceleration verifies its own momentum equation. This momentum equation for a small spherical particle was established in 1983 by [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF] under the assumptions that the particle of radius r p is much smaller than the characteristic macroscopic length scale L of the flow, that the particle Reynolds number based on the slip velocity of the particle is much lower than unity (so the disturbance flow due to the particle can be considered a Stokes flow), and that the diffusive timescale r 2 p /ν is much lower than the advective timescale L/U 0 (ν being the kinematic viscosity of the fluid and U 0 a characteristic velocity scale of the flow). With these assumptions, the leading terms boil down to the particle acceleration, its buoyancy and the Stokes drag exerted by the fluid. Furthermore, in the equilibrium Eulerian formalism, assuming that all particles have very small inertia (limit of vanishingly small Stokes numbers), the particle acceleration can be neglected [START_REF] Balachandar | Turbulent Dispersed Multiphase Flow[END_REF] so that the particle velocity is solely prescribed by the balance between buoyancy and drag [START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF][START_REF] Berk | Dynamics of small heavy particles in homogeneous turbulence: A Lagrangian experimental study[END_REF]. This balance yields

v p = v + w s e z , (1)
as derived in [START_REF] Nasab | Preferential concentration by mechanically driven turbulence in the two-fluid formalism[END_REF], where v is the fluid velocity, v p is the particle velocity, w s is the terminal velocity of the particle in quiescent fluid and e z = g/||g|| is aligned with the gravity field g. In the present study we use this approximation to the momentum equation even for particles with non negligible inertia, and try to assess the accuracy of this approach. To keep a moderate numerical cost, as in previous studies of particle-laden flows [26-30, 34, 35] we model particles as a continuum with a field of concentration C (mass of particles per unit volume) that is advected at the velocity v p , hence, mass conservation reads

∂ t C + v p • ∇C = 0 . (2)
In practice, such modelling neglects particle dispersion at the particle scale. In fact, even in dilute suspensions for which collisions are negligible, and even in the absence of Brownian motions, particles induce long-range perturbations, particle-wake interactions and collective-settling effects in the fluid depending on their particulate Reynolds number (e.g. [START_REF] Guazzelli | A Physical Introduction to Suspension Dynamics[END_REF][START_REF] Subramanian | Evolution of clusters of sedimenting low-Reynolds-number particles with Oseen interactions[END_REF][START_REF] Pignatel | A falling cloud of particles at a small but finite Reynolds number[END_REF][START_REF] Daniel | Clusters of sedimenting high-Reynolds-number particles[END_REF]). These perturbations lead to a dispersion which has been analysed analytically [START_REF] Koch | Hydrodynamic diffusion in a suspension of sedimenting point particles with periodic boundary conditions[END_REF], numerically [START_REF] Ladd | Dynamical simulations of sedimenting spheres[END_REF] and experimentally [START_REF] Ham | Hindered settling and hydrodynamic dispersion in quiescent sedimenting suspensions[END_REF][START_REF] Nicolai | Effect of the vessel size on the hydrodynamic diffusion of sedimenting spheres[END_REF][START_REF] Martin | Hydrodynamic dispersion broadening of a sedimentation front[END_REF][START_REF] Xue | Diffusion, dispersion, and settling of hard spheres[END_REF][START_REF] Lee | Combined effect of sedimentation velocity fluctuation and self-sharpening on interface broadening[END_REF], and which can be approximated by a diffusive process. Together with other ingredients such as concentration gradients and shear [START_REF] Davis | Hydrodynamic diffusion of suspended particles: A symposium[END_REF], these effects are referred to as a hydrodynamic diffusion [START_REF] Guazzelli | Fluctuations and Instability in Sedimentation[END_REF]. They are often approximated at a macroscopic level by a term of diffusion in the mass conservation above [START_REF] Nasab | Preferential concentration by mechanically driven turbulence in the two-fluid formalism[END_REF][START_REF] Martin | Hydrodynamic dispersion broadening of a sedimentation front[END_REF][START_REF] Lee | Combined effect of sedimentation velocity fluctuation and self-sharpening on interface broadening[END_REF][START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF]. In addition, including a diffusive term in the mass conservation equation is necessary to prevent the formation of caustics in the field of particle concentration, which would lead to numerical instabilities in the current Eulerian framework [START_REF] Magnani | Inertial effects in dusty Rayleigh-Taylor turbulence[END_REF]. By introducing the effective particle diffusivity κ p , and using equation (1), the new equation of mass conservation reads

∂ t C + v • ∇C = κ p ∇ 2 C -w s ∂C ∂z , (3)
where the last term accounts for the gravitational drift of particles. Fluid motions are constrained by the condition of incompressibility

∇ • v = 0, (4)
which also prescribes the incompressibility of the field of concentration using equation (1). The fluid velocity v verifies the Navier-Stokes equation for a Newtonian fluid

∂v ∂t + v • ∇v = - 1 ρ f ∇p + ν∇ 2 v + 1 ρ f f drag , (5)
where ρ f is the fluid density, p is the pressure field including the hydrostatic contribution, ν is the constant fluid kinematic viscosity, and f drag is the average drag force exerted by particles on the fluid per unit volume. This drag term in equation ( 5) describes how particles force the flow: the field of concentration and the fluid are now two-way coupled. The present model is derived for supposedly spherical particles of vanishingly small Reynolds number, hence they fall in the Stokes regime with a settling velocity

w Stokes s = 2gr 2 p (ρ p -ρ f ) 9νρ f , (6)
with ρ p the density of a spherical particle of radius r p . Then, using equation (1), the Stokes drag exerted by the fluid on a particle reads -6πρ f νr p w Stokes s e z . Summation on all the particles in the unit volume requires to multiply this individual acceleration by the number of particles per unit volume i.e. 3C/4πr 3 p ρ p . Using Newton's third law, the drag force reads

f drag = C ρ p -ρ f ρ p g . (7)
The previous equations are non-dimensionalised using the diameter of the release cylinder D cyl and the characteristic velocity that clouds can build up when accelerating from rest which is

U ref = g 1 - ρ f ρ 0 D cyl , (8)
where 4) is unmodified while mass conservation now reads

ρ 0 = ρ f + (1 -ρ f /ρ p )m 0 /(4πD
∂ t C + v • ∇C = 1 P e ∇ 2 C -R ∂C ∂z gravitational drift , (9)
where P e = U ref D cyl /κ p is the Péclet number, and R = w s /U ref is the Rouse number which characterises the gravitational drift of particles. When R ≪ 1 particles hardly drift with respect to the fluid and tend to follow fluid motions, whereas particles having R ≫ 1 largely decouple by vertical settling so their trajectories largely differ from those of fluid particles in their vicinity. Using equation ( 7), similar non-dimensionalisation of equation ( 5) leads to

∂v ∂t + v • ∇v = -∇p + 1 Re ∇ 2 v + RiCe z drag term , (10)
where

Re = U ref D cyl /ν is the Reynolds number and Ri = gD cyl (1 -ρ f /ρ p )/U 2 ref is the Richardson number which boils down to a ratio of density contrasts Ri = (1 -ρ f /ρ p )/(1 -ρ f /ρ 0 ) for our specific choice of U ref used to non-dimensionalise the equations.

III. NUMERICAL SETUP

The numerical setup aims at reproducing the experimental conditions for the generation of our particle clouds [1], released with no initial velocity from the top centre of a tank of still water. Every particle cloud is composed of m 0 = 1g of spherical glass beads with a mean radius r p chosen in the range 5µm-1mm, as well as some water, sometimes dyed with rhodamine. All this material is initially contained in a cylinder of diameter D cyl partially immersed in water and sealed by a latex membrane. When an experiment starts, the membrane is ruptured by a needle, it quickly retracts and lets particles settle from rest in the tank with the dyed fluid. Visualisations are performed in a vertical green laser sheet with two cameras: the one with a green filter records the motion of particles only, while the camera with an orange filter records only motions of turbulent eddies dyed with rhodamine, see figure 1a for an illustration. Some reference clouds without particle inertia contained an identical mass excess m 0 = 1g of salt water and were generated in the same way. Their dynamics is characterised by a Rouse number R = 0. Experiments were performed in a tank of depth 1m with horizontal square cross section of 42 × 42cm 2 surface area; the side walls of the tank were considered far enough from particle clouds to have negligible influence on the settling of particles. In the present simulations, the computational domain represents a cube of volume 1m 3 with side length L domain = 31.25D cyl . Figure 1b shows a numerical analog of our experimental clouds in a simulation having R = 0.221 at time t = 34.

To model the instantaneous release of particles from the cylinder, the field of concentration is initialised in a narrow cylindrical region of diameter D cyl , horizontally centred and localised at distance 1.5D cyl below the top of the computational domain. In experiments particles rested at the bottom of the cylinder on the latex membrane as a very thin layer of height H compact . To smooth the initial concentration gradients required to model the initial layer, the field of concentration is initialised in a cylindrical region of height 10H compact , and the lateral and vertical edges of this cylindrical region are smoothed by a hyperbolic tangent with typical length scale H compact . The factor 10 on the initial height of the pile of particles is expected to have negligible influence: indeed, our experiments showed that varying the height of the pile of particles from H compact to 40H compact had negligible impact on the dynamics we aim to model here (see [1] for further details). To compensate for the errors introduced the hyperbolic tangent smoothing, the exact mass of particles is enforced (in dimensional form, m 0 = 1g) by an appropriate rescaling. The uniform fluid is initially motionless, save for a cylindrical enveloppe around particles at time t = 0 in which the velocity field is perturbed by a random uniform infinitesimal noise.

A passive tracer is implemented to mimic rhodamine in simulations. This tracer is initialised exactly like the concentration C and satisfies equation ( 9) with its own Rouse number R = 0 so it does not drift vertically as particles do. Importantly, unlike the field of particle concentration C, the passive tracer does not force any flow (it is absent of equation ( 10)).

On all walls, the velocity field satisfies no slip boundary conditions, while the concentration C and the passive tracer satisfy no-flux Neumann boundary conditions.

Equations ( 4), ( 9) and [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF] are integrated in three dimensions with the solver Basilisk [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] that is second-order accurate in space and time, using a time-splitting pressure-correction discretisation of the Navier-Stokes equation [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF]. The second-order advection scheme of Bell-Colella-Glaz [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF] is used for equation [START_REF] Maxey | The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields[END_REF] and for the two advection for further details).

terms in equation [START_REF] Falkinhoff | Preferential Concentration of Free-Falling Heavy Particles in Turbulence[END_REF]. The concentration gradient ∂C/∂z is computed with the generalised minmod slope limiter [START_REF]Basilisk -src[END_REF] to reduce spurious oscillations due to sharp concentration gradients at initial times, with negligible impact on our numerical measurements (see appendix A). Due to the large scale separation between the domain size and the initial cloud size (L domain = 31.25D cyl ), an adaptive mesh refinement is adopted. This octree mesh is made of hierarchically organised cubic cells, each refinement of a cell corresponding to a division of this cell in 8 identical cubes. This refinement is based on local values of the concentration C and on the local viscous dissipation, see an illustration in figure 2a. The smallest mesh cell has a size fixed to L domain /1024 while the largest mesh cell has a size L domain /128.

Our main focus is to analyse the influence of gravitational settling on the evolution of particle clouds by varying the Rouse number from R = 0 to R = 3.03. Consequently, with a fixed density ρ p = 2500kg/m 3 for all particles, the Richardson number Ri = 138 does not vary between simulations, and the same conclusion holds for Re and P e. These last two numbers are equal because we assume ν = κ p as a first approximation. The value Re = 454 is chosen by following this conservative estimate: if the particle clouds were highly inertial and if their turbulence was homogeneous and isotropic with a fully-developed energy cascade from the integral cloud scale of order ∼ 3D cyl to the Kolmogorov scale, then the latter would have a size 3D cyl Re -3/4 . By prescribing that this smallest length scale should have a size L domain /1024, this prescribes the value of Re = (L domain /3072D cyl ) -4/3 = 454. This calculation is conservative as no flow structure reaches such a small length scale in our simulations. In fact, the transient cloud formation does not permit the development of an energy cascade down to the Kolmogorov scale. Instead, the clouds we model are close to laminar, with a predominating coherent structure of size ∼ 3D cyl , confirming that all length scales are appropriately resolved. An example of the fields obtained is shown in figure 2b. For completeness, additional simulations at Re = 1183 (which is the Reynolds number of particle clouds in our experiments [1]) are presented in appendix B which show similar results as those presented in the core of the present analysis.

IV. REGIMES OF CLOUD SETTLING

A. Overview

Figure 3 shows snapshots as well as an average image of the field of concentration C(x, y = 0, z, t) in a vertical cross section of the computational domain. They faithfully account for all the qualitative morphological features observed in our experiments (see figure 4 in Kriaa et al. [1]). Particle clouds initially grow as they propagate downward. For the lowest Rouse numbers this regime of growth and dilution lasts over the entire fall of the cloud. For intermediate Rouse numbers, this regime is more prominent (see figure 3c), yet this growth gradually vanishes and stops (see in particular average images in figures 3d, 3e). This transition is due the field of concentration separating from eddies, which is why the field of concentration deforms less and less, and settles more and more vertically.

The cloud evolution can be quantified by tracking the vertical position z f (t) of its front, defined as the lowermost position of the iso-contour of concentration C(x, y, z, t) = 10 -8 (this low value ensures that measurements are independent of the specific threshold), and by tracking the vertical position z(t) of the centre of mass, defined as the weighted average of the vertical coordinate z in the whole computational domain as z(t) = zC/ C. Figure 4 shows the evolution of these two positions in time; for all figures we use dashed lines if R > 0.221 and solid lines otherwise, with grey shades diverging from the thick darkest line of R = 0.221. For low Rouse numbers (R < 1), after a short phase of acceleration, the positions z f (t) and z(t) evidence a concave evolution revealing the cloud deceleration. For intermediate Rouse numbers, a smooth transition occurs that leads to a regime of constant settling velocity when particles separate from fluid motions and settle in quiescent liquid. For the largest Rouse numbers (R > 1), clouds almost immediately transition from the regime of acceleration to the regime of constant settling velocity. The literature [1, [START_REF] Rahimipour | Dynamic behaviour of particle clouds[END_REF][START_REF] Bush | Particle clouds in homogeneous and stratified environments[END_REF][START_REF] Lai | Modeling and experiments of polydisperse particle clouds[END_REF] distinguishes between the first regime of turbulent thermal, or equivalently of buoyant vortex ring, and the second regime of swarm. All these observations are consistent with our experimental observations, which were obtained by tracking the front of the light intensity reflected by the glass spheres in time, which is analogous to z f here. In the next sections, the two regimes of buoyant vortex ring and swarm are described separately.

B. Buoyant vortex ring regime (aka thermal regime)

Due to the lower cloud Reynolds number in numerical simulations than in experiments (see appendix B for results at higher Reynolds number), numerical results evidence clear buoyant vortex rings whose toroidal structure is less apparent in our more turbulent thermals in the laboratory. Yet, several studies have shown that turbulent thermals (i.e. instantaneous releases of a finite volume of buoyant fluid with no initial momentum) are formed like buoyant vortex rings [START_REF] Pottebaum | The pinch-off process in a starting buoyant plume[END_REF][START_REF] Gharib | A universal time scale for vortex ring formation[END_REF], and have a very similar structure [START_REF] Bond | Effects of initial geometry on the development of thermals[END_REF][START_REF] Lherm | Thermal and chemical partitioning dynamics during the differentiation of terrestrial planets[END_REF][START_REF] Lecoanet | Entrainment in Resolved, Dry Thermals[END_REF][START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF] and dynamics [START_REF] Turner | Buoyant vortex rings[END_REF][START_REF] Landeau | Experiments on the fragmentation of a buoyant liquid volume in another liquid[END_REF], to such an extent that thermals have been considered as buoyant vortex rings whose circulation is fully determined by their buoyancy [START_REF] Bond | Impact of buoyancy on vortex ring development in the near field[END_REF][START_REF] Zhao | On the internal vorticity and density structures of miscible thermals[END_REF]. Consistently, turbulent thermals at the cap of starting plumes [START_REF] Turner | The 'starting plume' in neutral surroundings[END_REF], immiscible thermals [START_REF] Landeau | Experiments on the fragmentation of a buoyant liquid volume in another liquid[END_REF] and settling particle clouds [START_REF] Ruggaber | Dynamics of sediment clouds related to open-water sediment disposal[END_REF][START_REF] Lai | Two-phase modeling of sediment clouds[END_REF][START_REF] Moghadaripour | Experimental Study of Particle Clouds in Stagnant Water[END_REF] have all been successfully modelled as vortex rings. In fact, different models of turbulent thermals (for example, references [2, 70]) and models of buoyant vortex rings [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF][START_REF] Turner | Buoyant vortex rings[END_REF][START_REF] Nikulin | Analytical model of motion of turbulent vortex rings in an incompressible fluid[END_REF] lead to the same scalings after the transient formation of these structures. Given the morphology of clouds in figure 3, we present the essential equations governing the evolution of buoyant vortex rings, which will prove enlightening to analyse the clouds growth.

When buoyant particles start settling, the cloud initially rolls up as a buoyant vortex ring. After a transient, one observes that vorticity concentrates inside a toroidal core. It is often assumed, and appropriate in our simulations as we shall see, that the vortex ring is sufficiently thin-cored to guarantee the absence of any buoyant material along the vortex centreline (axis of symmetry of the vortex ring, parallel to e z and passing through the centre of symmetry of the torus), and that no vorticity diffuses through this line, leading to the conclusion that the circulation Γ of the vortex ring remains constant after the short transient of spin-up (see equation 4 in reference [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF]), as verified in numerical simulations for a cloud Reynolds number of 630 and 6300 [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF]. Consequently, the initial spin-up of the vortex is paramount because it sets the ultimately constant value of the circulation which plays a key role in determining the cloud growth rate.

The impulse of a thin-cored vortex ring of radius R under the Boussinesq approximation reads πρ f ΓR 2 [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF][START_REF] Turner | Buoyant vortex rings[END_REF]. This impulse varies in time due to buoyancy which is the sole external force, hence we have

d dt πρ f ΓR 2 = m 0 g , (11)
from which it is clear that in the absence of buoyancy, the impulse would be constant hence the vortex ring would keep a constant radius, as verified in the literature [START_REF] Gharib | A universal time scale for vortex ring formation[END_REF][START_REF] Didden | On the formation of vortex rings: Rolling-up and production of circulation[END_REF] for vortex rings of Reynolds number up to 10 4 [START_REF] Didden | On the formation of vortex rings: Rolling-up and production of circulation[END_REF]. Assuming a constant circulation, the previous equation simplifies to

R(t) 2 -R 2 (t = 0) = m 0 g πρ f Γ t . (12)
To verify this scaling, the cloud radius is measured in numerical simulations with the quantity σ h , which is the horizontal standard deviation of the particles' spatial distribution with respect to the cloud centre of mass

σ h (t) = C(x ′ , y ′ , z ′ , t)(x ′2 + y ′2 ) C(x ′ , y ′ , z ′ , t) , (13)
where the origin of coordinates (x ′ , y ′ , z ′ ) corresponds to the cloud centre of mass, and integrals are computed in the whole computational domain. From this definition, the scaling ( 12) is verified in figure 5: aside from slight oscillations, clouds of low Rouse number (R ≤ 0.221) eventually grow as σ h ∼ t 1/2 as evidenced by the plateau of the compensated quantity [σ 2 hσ 2 h (t = 0)]/t. Conversely, clouds with large Rouse numbers (R > 0.221) grow slower than σ h ∼ t 1/2 so the quantity [σ 2 hσ 2 h (t = 0)]/t eventually decreases in time, meaning that these clouds do not follow the scaling [START_REF] Aliseda | Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence[END_REF]. Note that according to equation ( 12), the plateaus in figure 5 are inversely proportional to the circulation Γ, suggesting the existence of a minimum of circulation for intermediate Rouse numbers close to R = 0.221; this observation will receive considerable attention in section V. The clouds' growth is generally analysed along depth z rather than in time. Figure 6a shows the evolution of the radius σ h (z) for all clouds. For R ≪ 1 clouds tend to grow linearly in depth after a short transient. For R = 0.221 a noticeable decrease of the slope dσ h /dz is visible at depth z = 10, as already observable in figure 3c, which is due to the transition to the swarm regime. For large Rouse numbers, typically above unity, it is manifest that the slope dσ h /dz is never constant, it keeps decreasing as clouds fall deeper.

More importantly, figure 6a shows that the growth rate dσ h /dz(z) tends to be maximum for an intermediate Rouse number close to R = 0.221. To compare this growth rate with the entrainment rate measured in our experiments over a 45cm-deep field of view, the value of dσ h /dz(z) is computed in simulations, then averaged in the range z < 0.45L domain . The resulting average value is denoted α and results are shown in figure 6b. Consistently with our experiments (see figure 6 in [1]), we observe an asymptote towards a constant growth rate as R → 0, a local maximum for a finite Rouse number which lies in the range R ≃ 0.3 ± 0.1 determined from experiments, and a clear decrease of the growth rate as R increases beyond R = 0.3. The maximum amplitude of the enhancement α(R = 0.221)/α salt = 3.20 ± 0.83 is larger than the value found in experiments α(R ≃ 0.3)/α salt = 1.75 ± 0.30. This is likely due to the difference of cloud Reynolds number between experiments (Re = 1183) and simulations (Re = 454). In appendix B, numerical results for clouds of Reynolds number Re = 1183 yield α(R = 0.221)/α salt = 2.20 ± 0.42 as a new maximum, which is in agreement with our experiments. The important point here is that the present Eulerian approach is capable of reproducing the physical effect we observed in laboratory experiments.

We now turn to the swarm regime before analysing the role of the gravitational drift in enhancing the growth rate α in Section V.

C. Swarm regime

During the gradual transition from the buoyant vortex ring to the swarm regime, particles increasingly decouple from fluid motions. This gradual decoupling is analysed thanks to the concentration C tracer of the passive tracer, which is implemented in four new simulations of Rouse numbers R = 0.100, 0.221, 0.498, 0.885. As a wish to perform similar processings as in our experiments where the flow was visualised in a vertical laser sheet, the decoupling is illustrated with space-time diagrams in figures 7a-7d after averaging the fields C and C tracer along x in the plane y = 0 (this averaging process is made explicit with the brackets ⟨•⟩ x ). The space-time diagrams reveal that increasing R leads to a faster separation between the field ⟨C⟩ x (in blue-red colours) and the passive tracer ⟨C tracer ⟩ x (in yellow-red contours) which is left behind particles. By construction, the sole ingredient responsible for this decoupling is the gravitational drift in equation ( 9) which only concerns the field C, not the passive tracer. A quantification of the decoupling between C and C tracer is possible by defining a correlation coefficient which computes the overlap between C and C tracer as a percentage of the total region occupied by these two fields. In practice, minimum thresholds are applied on the fields, whose values ensure the convergence of the correlation coefficient. As in figure 7, the correlation coefficient is computed in the plane y = 0 (denoted Oxz below) as

C C,Ctracer (t) = Oxz ξ 1 (t)dxdz Oxz ξ 2 (t)dxdz , ( 14
)
where the booleans (ξ 1 (t), ξ 2 (t)) are defined as

ξ 1 (t) = 1, if C(t) > C 1 (t) and C tracer (t) > C 2 (t) 0, otherwise ; ξ 2 (t) = 1, if C(t) > C 1 (t) or C tracer (t) > C 2 (t) 0, otherwise , (15)
with C 1 (t) = max{C(t)} × 10 -4 , C 2 (t) = max{C tracer (t)} × 10 -4 . The time evolution of C C,Ctracer (t) is shown in figure 8a for the same four simulations as above. The correlation coefficient is initially equal to unity because both fields of concentration are identical. Then, the correlation coefficient decreases as particles gradually shift away from the passive tracer due to the gravitational drift, all the faster as the Rouse number is larger. Finally figure 8b confirms that the settling velocity of swarms approaches a constant value after separation, equal to the individual settling velocity which is equal to the Rouse number R in our dimensionless units (see inlet in figure 8b). Figure 6a already showed that the growth rate of swarms starts reducing after separation (the concave deflection of σ h (z) is most visible for R = 0.221), so swarms ultimately fall with constant velocity, retaining a bowl shape without deforming, as shown by the snapshots in figure 3e for R = 3.03.

V. ROLE OF THE ROUSE NUMBER ON THE ENHANCED GROWTH RATE

The linear growth in depth of buoyant clouds is usually described as resulting from entrainment of ambient fluid into the cloud, leading to its growth and dilution, and consequently its deceleration through mixing drag. One reference to describe the process of entrainment for turbulent thermals is the model of Morton et al. [2]. It is based on the entrainment hypothesis which states that the inflow velocity at the interface of the turbulent thermal is proportional to the vertical velocity of the cloud's centre of mass; this inflow velocity is considered to be produced by turbulent motions.

Yet, our observation of a maximum growth rate α even for moderate Reynolds numbers (Re = 454 here) suggests that the enhanced entrainment due to the finite gravitational drift finds an origin in the large-scale buoyancy-induced mean flow (which, in the case of turbulent particle clouds, is obtained by an average over realisations) rather than in turbulent fluctuations. Actually, this question of the origin of entrainment in turbulent flows has long been debated: does it originate from the large-scale mean flow incorporating ambient fluid into the turbulent region through 'engulfment' ? Or from small-scale fluctuations and diffusive processes which mix the ambient material in the turbulent structure close to its interface through 'nibbling' ? While Mathew and Basu [START_REF] Mathew | Some characteristics of entrainment at a cylindrical turbulence boundary[END_REF] observed that mixing in a cylindrical turbulent jet seemed to be driven by nibbling close to the jet interface, Townsend [START_REF] Townsend | LXXXI. The eddy viscosity in turbulent shear flow[END_REF] showed that large-scale eddies of increasing intensity produce more energy in the turbulent wake past a cylinder and consequently favour a larger growth rate of this wake, suggesting that the mean flow drives entrainment through engulfment. Discriminating between engulfment and nibbling can be complex because large and small scales may be insufficiently separated at moderate Reynolds numbers, and because fluxes at both scales can be connected through some relationships [START_REF] Mathew | Some characteristics of entrainment at a cylindrical turbulence boundary[END_REF]. For example, in [START_REF] Odier | Understanding and modeling turbulent fluxes and entrainment in a gravity current[END_REF], Reynolds stresses at the mixing interface of a gravity current are modelled based on Prandtl mixing length theory, thus enabling, through the use of a large-scale quantity based on the mean flow, the description of entrainment in the mixing layer by fundamentally local fluxes.

Importantly, the reason for this ambiguity is that both processes contribute to entrainment, as evidenced by Fox [START_REF] Fox | Forced plume in a stratified fluid[END_REF] who derived the equations of evolution of a self-similar plume while considering the equation of conservation of energy, hence lifting the constraint of modelling entrainment to guarantee a closure of the equations. In doing so, he showed that entrainment depends both on the Reynolds stress and on a contribution from the mean flow due to buoyancy. This was confirmed by Reeuwijk and Craske [START_REF] Van Reeuwijk | Energy-consistent entrainment relations for jets and plumes[END_REF] who carried this analysis further and showed that a third contribution comes from possible deviations from self-similarity in the streamwise direction. By analysing data from the literature, Reeuwijk and Craske [START_REF] Van Reeuwijk | Energy-consistent entrainment relations for jets and plumes[END_REF] showed that the term of turbulence production due to shear hardly varies between a pure jet and a pure plume, even though plumes have a larger growth rate than jets. Hence this last difference between jets and plumes is attributed to the contribution of the mean flow due to buoyancy, as later confirmed by Reeuwijk et al. [START_REF] Van Reeuwijk | Turbulent transport and entrainment in jets and plumes: A DNS study[END_REF] in DNS. These conclusions about entrainment apply similarly for plumes and thermals, as pointed out by Landeau et al. [START_REF] Landeau | Experiments on the fragmentation of a buoyant liquid volume in another liquid[END_REF]. From the model of Morton et al. [2], Landeau et al. [START_REF] Landeau | Experiments on the fragmentation of a buoyant liquid volume in another liquid[END_REF] proved experimentally that the growth rate of an immiscible thermal verifies a linear relationship with respect to the thermal's Richardson number, which is exactly analogous to the entrainment model of Priestley and Ball [START_REF] Priestley | Continuous convection from an isolated source of heat[END_REF] for plumes, see [START_REF] Van Reeuwijk | Energy-consistent entrainment relations for jets and plumes[END_REF] for details. Consistently with these conclusions, Lecoanet et al. [START_REF] Lecoanet | Entrainment in Resolved, Dry Thermals[END_REF] showed that turbulence only enhances entrainment by 20% between turbulent thermals of Reynolds number 630 and 6300, whereas an artificial sudden shut off of buoyancy during the same numerical simulations drastically reduces the entrainment rate of turbulent thermals [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF], highlighting the driving role of buoyancy in entrainment.

Consequently, in the present section we focus on mean flow azimuthally-averaged quantities to try and understand the optimum growth rate of particle clouds for a Rouse number around R ≃ 0.22. Past studies have shown that in the absence of buoyancy, the circulation of a vortex ring generated from a nozzle increases from zero to a constant value during the transient rolling-up of the viscous boundary layer in the nozzle [START_REF] Didden | On the formation of vortex rings: Rolling-up and production of circulation[END_REF], with this constant circulation increasing as the ratio of the nozzle length over its diameter increases. When the vortex ring is buoyant, buoyancy provides a new contribution to the total circulation, which was fully derived by Mc Kim et al. [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF] for a thin-cored Boussinesq vortex ring, yielding the following scaling for the growth rate α (see [START_REF] Bond | Effects of initial geometry on the development of thermals[END_REF][START_REF] Nikulin | Analytical model of motion of turbulent vortex rings in an incompressible fluid[END_REF] but also [START_REF] Lecoanet | Entrainment in Resolved, Dry Thermals[END_REF][START_REF] Turner | Buoyant vortex rings[END_REF][START_REF] Scorer | Experiments on convection of isolated masses of buoyant fluid[END_REF] for turbulent thermals)

α ∝ m 0 g ρ f Γ 2 ∞ , (16)
with a proportionality constant accounting for the cloud added mass and its morphology [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF], and with Γ ∞ the asymptotically constant value of the circulation. The key result here is that since all clouds undergo the same buoyancy force m 0 g in our simulations, their growth rate α is dictated by Γ ∞ only. Consequently, vortex rings of lower circulation should have a larger growth rate -crudely speaking, they should entrain more, as consistently observed in [START_REF] Bond | Effects of initial geometry on the development of thermals[END_REF][START_REF] Landeau | Experiments on the fragmentation of a buoyant liquid volume in another liquid[END_REF].

To verify this scaling, the vortex ring circulation is computed a posteriori for several snapshots after (i) interpolating the mesh on a regular cartesian grid, (ii) averaging the azimuthal vorticity ω θ along the azimuth e θ in radial and vertical bins (when ambiguity is possible, we denote azimuthally-averaged quantities with an overline such as ω θ ), (iii) integrating the resulting average azimuthal vorticity over the whole radial extent and in the range z ∈ [0.3, L domain -0.3] (top and bottom walls are removed from this range to avoid integration of vorticity near those boundaries). Results are shown in figure 9a. We verify that the circulation produced by the particle clouds ultimately reaches a plateau, except when R ≫ 1 i.e. for clouds which are not expected to behave as buoyant vortex rings due to separation. From previous sections, figure 9a already suggests that the circulation Γ ∞ is all the lower as the cloud grows faster; this is confirmed in figure 9b where the scaling α ∝ Γ -2 of equation ( 16) is in excellent agreement with the best fit (see the solid dark line in figure 9b) of measurements of the circulation when the cloud centre of mass is at mid-depth in the computational domain i.e. with the definition Γ ∞ ≡ Γ(z = L domain /2).

Another remarkable observation is the modification of the structure of the vortical core with increasing the Rouse number, as illustrated in figure 10. When R = 4.28 × 10 -2 (figure 10a) the vortex core is neatly defined, centered around a maximum of azimuthal vorticity whose structure is at first order isotropic in a plane (e r , e z ) of fixed azimuth. Conversely when R = 0.221 (figure 10c), the vortex core is made up of sheets of vorticity of alternate sign and varying intensity which occupy a much larger region than observed for R = 4.28 × 10 -2 at the same depth. Let us show that as R gets closer to 0.221 the vortical core induces circulation farther and farther away, thus expanding the region of entrainment through the toroidal mean flow. After defining the core centroid as the barycentre of the azimuthal vorticity, the following normalised circulation is computed

Γ * (r) = 1 πr 2 2π 0 r 0 < ω θ > θ r ′ dr ′ dθ, (17)
where (r ′ , r, θ) here correspond to polar coordinates centered on the circulation centroid. The size of the vortex core is defined as argmax{Γ * (r)} i.e. as the radial distance from the circulation centroid where the circulation Γ * (r) is maximum. Consistently, figure 11a shows that this circulation Γ * spreads further away from the circulation centroid when R = 0.221, meaning that this vortex ring induces velocity farther away, therefore incorporating more ambient fluid within the particle cloud, hence the latter grows faster. This trend is even clearer in figure 11b where the core size is computed at the same depth L domain /2 for all clouds: one verifies that the core extension is maximum when the Rouse number is closest to R = 0.221.

The key question is then: How come vortex rings have a wider core and a lower circulation when the particle Rouse number is closer to 0.221? The vortex circulation is only produced during a short initial transient, mainly by the baroclinic torque as long as its contribution along the vortex centreline is non-negligible. As soon as the buoyant material has spun up and widened sufficiently, one can define a closed contour encircling the core where no vorticity diffuses and no buoyant material is present, so that circulation is conserved (e.g. [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF]). The vorticity equation along the azimuth e θ reads

Dω θ Dt = ω θ v r r stretching + baroclinic torque g ′ ρ f ∂C ∂r + ν 1 r ∂ ∂r r ∂ω θ ∂r + ∂ 2 ω θ ∂z 2 diffusion , (18)
which shows that azimuthal vorticity is produced by vortex stretching, the baroclinic torque and diffusion of vorticity. The quantity ω θ v r is at first order symmetrical around the vortex core, so for a thin-cored vortex ring having a radius much larger than the vortex core, the stretching term should be vanishingly small, as previously argued by other authors [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF]. Then, if diffusion is assumed negligible, most forcing is expected to originate from the baroclinic torque. This is especially true at initial times when the stretching and diffusion terms vanish while the baroclinic torque remains finite. The baroclinic torque is therefore the leading source of circulation at initial times [START_REF] Mckim | Buoyancy-driven entrainment in dry thermals[END_REF]. To verify this, the dimensionless baroclinic torque Ri∂ r C along e θ is first averaged along the azimuth (the resulting axisymmetric torque is denoted Ri∂ r C). Then the axisymmetric torque is integrated in the plane (e r , e z ) and integrated in time until t = 10 when we observe that the torque has vanished for all simulations. Results are shown in figure 12. We verify that Rouse numbers close to R = 0.221, which correspond to the largest growth rate α (figure 6b) and lowest circulation (figure 9a), also correspond to the lowest baroclinic forcing. This observation is robust: integrating the baroclinic torque in time even just up to t = 1 modifies the value of the integrated torque, but leaves the curve in figure 12 unchanged. As an indication, the errorbars in figure 12 show the little influence of integrating the baroclinic torque up to t = 8 or t = 10 which respectively correspond to the lower and upper bound of errorbars. The picture that emerges from these results is the following: under the assumption that baroclinicity is the leading forcing of the vortex rings' circulation, the maximum of entrainment capacity of particle clouds with a Rouse number close to R = 0.221 seems to be due to the gravitational drift and two-way coupling of particles with the fluid which reduces the baroclinic torque, thus reducing the cloud circulation. Since all clouds undergo an identical buoyancy force m 0 g, these same clouds have a larger growth rate α as predicted by equation ( 16). Similar results have been obtained at a larger Reynolds number Re = 1183 (see Appendix B), in good agreement with experiments. Note however that the role played by fluctuations in the limit of very large Reynolds numbers remains to be explored and might have an influence on our conclusions derived from moderate Reynolds numbers only.

VI. CONCLUDING DISCUSSION

The previous section showed that the gravitational drift modifies the distribution in space and time of the field of particle concentration C(x, y, z, t) compared to that of a passive tracer. This modification alters the forcing by the drag force, minimises the baroclinic torque and concurs to a maximum growth rate α for a Rouse number around R = 0.221. All these modifications are observed at moderate Reynolds numbers and notably quantified by the azimuthally-averaged circulation and baroclinic torque; this is consistent with the literature pointing towards the leading role of the mean flow and buoyancy in controlling entrainment and the growth of thermals (see section V).

A key conclusion is that the present Eulerian two-way coupling numerical simulations successfully reproduce our experimental observation of a maximum growth rate α/α salt of particle clouds for a Rouse number R ≈ 0.22 lying within the experimental range R ≃ 0.3 ± 0.1. While results at Re = 454 yield a maximum growth rate slightly above the experimental value α(R ≃ 0.3)/α salt = 1.75 ± 0.30, results at Re = 1183 lie in the experimental range within uncertainty margins. A systematic study with varying Reynolds numbers might clarify whether the mechanism identified in this paper persists in the presence of intense turbulent fluctuations.

Our results raise a new question: How does the gravitational drift of particles contribute to reducing the baroclinic forcing ? Some light could be shed on this matter by analysing the properties of the flow induced by a canonical laminar vortex ring while the field of concentration drifts radially outward until separation, but the evolution of the structure of the vortex core in figure 10 suggests that the feedback of particles on the vortex ring itself probably plays a non-negligible part. Furthermore, even though the robust agreement between our experiments and the present results supports the responsibility of the mean flow in the maximisation of α for a finite Rouse number, it remains to be investigated whether other physical ingredients could be at play in experiments, in particular velocity fluctuations due to turbulence at much higher Reynolds numbers than considered here.

In the present numerical simulations, fluctuations are very low compared to the mean flow due to the low Reynolds number Re = 454 at the scale of the particle cloud. Our experiments, on the opposite, were characterised by a Reynolds number Re = 1183. Even though this value is too low to have a well-developed turbulent flow with a clear separation of scales between the integral cloud scale and the scale of the smallest dissipative eddies, our experimental clouds evidenced some more fluctuations than in numerical simulations. On one hand, these may modify the cloud circulation during its transient increase, hence during a limited amount of time. On the other hand, after this transient, entrainment can be increased by the term of production of turbulent kinetic energy (TKE), as shown by Reeuwijk and Craske [START_REF] Van Reeuwijk | Energy-consistent entrainment relations for jets and plumes[END_REF] in their entrainment relations. This production term may differ for one-phase turbulent thermals vs. particle-laden turbulent thermals due to turbulence modulation by particles [START_REF] Brandt | Particle-Laden Turbulence: Progress and Perspectives[END_REF], as observed in simulations [START_REF] Elghobashi | On the two-way interaction between homogeneous turbulence and dispersed solid particles. I: Turbulence modification[END_REF][START_REF] Bosse | Small particles in homogeneous turbulence: Settling velocity enhancement by two-way coupling[END_REF] and experiments [START_REF] Hassaini | Scale-to-scale turbulence modification by small settling particles[END_REF][START_REF] Berk | Dynamics of small heavy particles in homogeneous turbulence: A Lagrangian experimental study[END_REF]. These studies notably showed a redistribution of energy from small to large wave numbers known as 'pivoting', which may favour nibbling-like entrainment at small scales rather than engulfment by the mean flow. Consequently, the possible enhancement of α by fluctuations in a more vigorous turbulent flow cannot be ruled out, calling for further investigation with a dedicated larger experimental setup and numerical simulations with a clear separation of scales. Answering these questions also probably requires a more advanced model such as a two-fluids approach where the particles have their own velocity field (e.g. [START_REF] Nasab | Preferential concentration by mechanically driven turbulence in the two-fluid formalism[END_REF][START_REF] Magnani | Inertial effects in dusty Rayleigh-Taylor turbulence[END_REF][START_REF] Nakamura | Linear stability analysis of bubble-induced convection in a horizontal liquid layer[END_REF]), or a point-force Lagrangian model (e.g. [START_REF] Chou | Numerical study of particle-induced Rayleigh-Taylor instability: Effects of particle settling and entrainment[END_REF][START_REF] Yamamoto | Numerical simulation of concentration interface in stratified suspension: Continuum-particle transition[END_REF][START_REF] Climent | Large-Scale Simulations of Bubble-Induced Convection in a Liquid Layer[END_REF][START_REF] Mazzitelli | Evolution of energy in flow driven by rising bubbles[END_REF]). I: Variability of the macroscopic quantities α, ż and the work of the drag force Cv z in the range z < 0.45L domain when varying the Schmidt number, the size h min of the finest mesh cell, or the numerical scheme used to compute the concentration gradient.

transfer of energy from particles to the fluid during the cloud fall. These three quantities are averaged when the cloud position verifies z(t) < 0.45L domain , which is the range we analysed with our experiments [1]. Reported values evidence little to negligible impact of the three numerical parameters on average measurements of α, ż and Cv z .

  FIG. 1: (a) Visualisation of an experimental particle cloud in a vertical laser sheet with particles in grey (set of polydisperse particles of Rouse number R = 0.308 ± 0.080) and dye (rhodamine) in orange. (b) Numerical analog in the plane y = 0 of the previous photograph with a grey field of concentration C modelling particles and an orange tracer concentration C tracer modelling dye (cloud of Reynolds number Re = 1183 as in the experiment (a), with a fixed Rouse number R = 0.221). The larger vertical spread of particles in (a) compared to (b) is due to the polydispersity of the former, whereas simulations are performed for monodisperse particle clouds (see Kriaa et al. [1]for further details).

  FIG. 2: (a) Adaptive mesh refinement on the field of concentration C for R = 0.221 in the plane y = 0 at time t = 35. (b) Bird's eye view of the 3D structure of C(x, y, z, t = 33.75) for R = 0.221 in blue-red colours (colorbar and opacity in the bottom left-hand corner). Streamlines show the toroidal velocity field, with blue-red colours for the velocity magnitude ||v|| (colorbar in the top left-hand corner).

FIG. 3 :

 3 FIG. 3: Each row shows snapshots of the field of concentration C(t, x, z) in the plane y = 0 for a different Rouse number, as well as the average of 40 snapshots taken with a constant timestep ∆t = 1 over the fall.

FIG. 4 :

 4 FIG. 4: For several Rouse numbers, evolution in time of the vertical position of (a) the centre of mass of particles, and of (b) the cloud front. Lines are dashed if R > 0.221 and solid otherwise.

FIG. 5 :

 5 FIG. 5: Time evolution of the increment of square radius compensated with time [σ 2 h (t)σ 2 h (t = 0)]/t, measured in simulations for several Rouse numbers. Lines are dashed if R > 0.221 and solid otherwise.

FIG. 6 :

 6 FIG. 6: (a) Evolution of the clouds' radius σ h along depth. Lines are dashed if R > 0.221 and solid otherwise. (b) Growth rate α computed in the range z < 0.45L domain , divided by the reference value α salt of a salt-water cloud i.e.of a cloud with no particle settling (R = 0).

FIG. 7 :

 7 FIG. 7: Hovmöller diagrams for four increasing Rouse numbers showing the gradual decoupling in the plane y = 0 between the field of concentration ⟨C⟩ x (in blue-red colours) and the passive tracer ⟨C tracer ⟩ x (in yellow-red contours).The values shown for both concentrations are horizontal averages along x in the plane y = 0 at each time step.

FIG. 8 :

 8 FIG. 8: (a) Time evolution of the correlation coefficient C C,Ctracer (t) between particle and passive tracer concentrations for R ∈ {0.100, 0.221, 0.498, 0.885}, from top to bottom. (b) Evolution in depth of the vertical velocity of the clouds' centre of mass for all simulations. Lines are dashed if R > 0.221 and solid otherwise. The inset shows the normalised velocity ż/R as a function of time t only for clouds verifying R ≥ 0.221.

FIG. 9 :

 9 FIG. 9: (a) Time evolution of the axisymmetric circulation in the whole computational domain. Lines are dashed if R > 0.221 and solid otherwise. (b) Correlation between the growth rate α averaged over the range z < 0.45, and measurements of Γ ∞ . The solid dark line is the linear least square fit of ln(α) vs. ln(Γ ∞ ).

FIG. 10 :

 10 FIG. 10: Each row shows snapshots of the azimuthally-averaged vorticity ω θ in blue-red colours, and the azimuthally-averaged concentration C with dashed contours in yellow-red shades (see colorbars on the left-hand side). Each snapshot is visualised in the plane (r, z) in the range 3 < z < 22 and 0 < r = x 2 + y 2 < 10.

FIG. 11 :

 11 FIG. 11: (a) Time evolution of the size of the vortical core for all clouds. Lines are dashed if R > 0.221 and solid otherwise. (b) Size of the core at depth z = L domain /2 for all clouds. The dark arrow corresponds to the value for R = 0.

FIG. 12 :

 12 FIG. 12: Results of time and volume integration up until t = 10 of the axisymmetric baroclinic torque as a function of the Rouse number.

Appendix B :

 : Key results for clouds of Reynolds number Re = 1183 For completeness, we briefly present numerical results for clouds of Reynolds number Re = 1183 as in our experiments [1]. Figure 13 shows time averages in the plane y = 0 of the field of concentration C for a set of clouds with varying Rouse numbers. More fluctuations do appear at low Rouse numbers (figures 13a-13c), whereas clouds of large Rouse numbers still evidence a thin bowl shape due to rapid particle separation from fluid motions (figures 13d,13e). Importantly, these images show the existence of a maximum growth rate as R → 0.221, as quantitatively confirmed by figure 14, similarly as in section IV B and consistently with our experiments [1]. The maximum amplitude of the enhancement α(R = 0.221)/α salt = 2.20 ± 0.42 is lower than the one measured for Re = 454 in section IV B, and in good agreement with the maximum measured in experiments α(R ≃ 0.3)/α salt = 1.75 ± 0.30. A systematic study as a function of the Reynolds number might clarify what determines the amplitude of the optimum α/α salt , but this is beyond the scope of the present work which focuses on the origin of this amplification, and on the capacity of the numerical model to reproduce it.

FIG. 13 :

 13 FIG. 13: Field of concentration C in the plane y = 0 averaged over up to 20 snapshots taken with a constant timestep ∆t = 2.5 over the cloud fall for clouds of Reynolds number Re = 1183.

FIG. 14 :

 14 FIG.14:Growth rate α computed in the range z < 0.45L domain , divided by the reference value α salt of a salt-water cloud i.e. of a cloud with no particle settling (R = 0). The Reynolds number is Re = 1183.

  Two-way coupling Eulerian numerical simulations of particle clouds settling in a quiescent fluid Quentin Kriaa, * Benjamin Favier, † and Michael Le Bars ‡ Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France (Dated: October 25, 2023) To get a deeper understanding of our laboratory experiments [Kriaa et al., Phys. Rev. Fluids 7(12) 124302 (

  3 cyl /3) is the typical initial effective density of particle clouds (see Ref. [1] for details). Time t, pressure p and concentration C are respectively non-dimensionalised by the advective timescale D cyl /U ref , the characteristic dynamic pressure ρ f U 2 ref , and the fluid density ρ f . With the dimensionless variables, equation (

TABLE

VII. ACKNOWLEDGEMENTS

This work was supported by the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by CNES. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. This work was granted access to the HPC resources of IDRIS under the allocation 2022-A0120407543 and 2023-A0140407543 made by GENCI.

Appendix A: Robustness of numerical measurements

We verified that numerical measurements of key quantities are invariant with respect to three numerical parameters: (i) the Schmidt number Sc = ν/κ p , (ii) the size h min of the finest mesh cell and (iii) the numerical scheme implemented to compute the concentration gradient. For all simulations presented in the core of this study, the concentration gradient was computed with the generalised minmod slope limiter [START_REF]Basilisk -src[END_REF], which reads along the direction e z

with θ a scalar ranging between 1 (the most dissipative scheme) and 2 (the least dissipative scheme), and C j∈N is the evaluation of C in a mesh cell j along the direction e z . The default value θ = 1.3 of Basilisk was adopted. We verified that adopting a second-order centered scheme

does not alter our measurements. Results are presented in table I, providing the average and standard deviation of (a) the growth rate α, (b) the cloud vertical velocity ż and (c) the work of the drag term Cv z responsible for the


Sylvie Masse

Guillaume Laurent

Thibaud Coradin

Andrzej Pisera

•

  Sponge spicules are biocomposite materials composed of a siliceous skeleton embedded in an organic matrix • Inner skeleton is of various shapes depending on the taxon, but surprisingly 29 Si HPDec ssNMR signature is quite ever the same and a condensation degree of ca. 0.93 is usually observed • 13 C CP-MAS ssNMR should be more promising to discriminate the samples but due to the low abundance of13 C and low content of carbon into the cleaned samples, the signal is too poor to get detailed assignment and to explore the organicmineral interface through 2D heterocorrelation NMR mapping • Dynamic Nuclear Polarization (DNP-enhanced ssNMR) should be a helpful technique to go further in the comprehension of the biomineralization process of sponges Specific limitations of our system: • Only << 1 wt% C into acid-cleaned samples • No isotopic enrichment possible • Only methodological and instrumental developments could help to get the signal faster A. Pisera on board of the IRD ship for an expedition in New Caledonia


Nicolas Courbin
email: ncourbin@gmail.com

Emmanuelle Dortel
email: emmanuelle.dortel@gmail.com

David Grémillet
email: david.gremillet@cebc.cnrs.fr

Jean-Dominique Lebreton
email: jean-dominique.lebreton@cefe.cnrs.fr

Aurélien Besnard
email: aurelien.besnard@cefe.cnrs.fr

à la gestion et à la conservation des populations d'oiseaux marins nicheurs du littoral français métropolitain Note on demography to assist in the management and conservation of breeding seabird populations along the French mainland coast

200-250 mots

synthétisé les paramètres démographiques connus pour 27 des 28 espèces d'oiseaux marins nicheurs en France. Toutes ces espèces ont une importante longévité (10 à plus de 50 ans) due à une survie adulte élevée (attendue > 0.90), une maturité sexuelle différée (2 à 9 ans), et une fécondité faible à moyenne (~1 à 3 oeufs). Le taux de croissance de ces populations est ~10 fois plus influencé par le taux de survie adulte que par des changements dans les paramètres reproducteurs, bien que les espèces les plus fécondes puissent compenser en partie une faible survie adulte. C'est ce que nous démontrons en modélisant la dynamique de population de deux espèces aux trait d'histoire de vie contrastés. Modifier la trajectoire des populations par des actions de gestion améliorant la fécondité localement n'est possible que si le succès reproducteur est déjà faible et si l'amélioration est forte. Finalement, la conservation des oiseaux marins nicheurs français devrait reposer essentiellement sur des actions favorisant la survie adulte. Cependant, les risques de mortalité identifiés sont globalement liés aux captures accidentelles, aux collision avec les éoliennes et à l'appauvrissement de la ressource, largement hors du cadre d'actions des gestionnaires locaux. Des mesures coordonnées au niveau régional, national et international améliorant la survie des espèces s'imposent.

MOTS CLÉS

Espèce longévive, menaces anthropiques, mesures de gestion, stratégie démographique, survie adulte.

ABSTRACT (+détaillé que le français)

The intent of this note is to present a general framework based on the demography of seabird populations in order to identify relevant management leverages for their conservation. First, we reviewed the influence of life history traits on demography and presented the demographic drivers expected along the slow-fast continuum hypothesis. Then, we summarized the known demographic parameters for 27 of the 28 seabird species breeding in France. All these species are on the slow extremity of the slow-fast continuum and have a long life span (10 to over 50 years) linked to a high adult survival (expected > 0.90), delayed sexual maturity (2 to 9 years), and low to medium fecundity (~1 to 3 eggs). The growth rate of these populations is ~10 times more influenced by adult survival than by changes in reproductive parameters, although the most fecund species (cormorants, gulls and terns) may partly compensate for low adult survival. We demonstrate this by modelling the population dynamics of two species with contrasting life history traits, a very long-lived and few productive species, the Scopoli's shearwater (Calonectris diomedea Scopoli, 1769), and a less long-lived and very fecund species, the European shag (Gulosus aristotelis Linnaeus, 1761). Modifying the trajectory of populations through management actions that improve fecundity locally is only possible if reproductive success is already low and if the improvement is strong. Finally, the conservation of French breeding seabirds should be based essentially on actions that favor adult survival. However, the risks of mortality identified are globally linked to accidental captures, oil spills, collisions with wind turbines and the depletion of the resource, largely outside the scope of actions by local managers.

Migratory species also face human disturbances at-sea off foreign countries. Coordinated measures at regional, national and international level to improve the survival of the species are needed.

KEYWORDS

Long-lived species, anthropic threats, management measures, demographic strategy, adult survival.

La dynamique de population est dépendante des traits d'histoire de vie

Les espèces d'oiseaux et de mammifères présentent une grande diversité de traits d'histoire de vie (croissance, reproduction et survie), qui s'échelonnent le long d'un gradient lent-rapide d'histoires de vie [START_REF] Stearns | The influence of size and phylogeny on patterns of covariation among life-history traits in mammals[END_REF][START_REF] Pontier | An analysis of demographic tactics in birds and mammals[END_REF][START_REF] Promislow | Living fast and dying young: a comparative analysis of life-history variation among mammals[END_REF]. Les espèces dites rapides sont caractérisées par une maturité sexuelle précoce, une forte fécondité et la production de nombreux jeunes à chaque évènement de reproduction, de faibles taux de survie et par conséquence une courte durée de vie. À l'extrémité lente du gradient, on retrouve les espèces longévives, à maturité sexuelle tardive, à faible fécondité, produisant généralement un seul jeune par évènement de reproduction, et caractérisées par des probabilités de survie juvénile, immatures et adulte élevées. La position des espèces le long du gradient d'histoires de vie est corrélée avec leur temps de génération [START_REF] Yoccoz | Generation time: a reliable metric to measure life-history variation among mammalian populations[END_REF] et le ratio entre le taux de fécondité et l'âge à la première reproduction [START_REF] Oli | The relative importance of life-history variables to population growth rate: Cole's prediction revisited[END_REF]. Elle est aussi fortement corrélée à la masse corporelle et à la phylogénie mais l'existence de ce continuum persiste après prise en compte de ces deux éléments [START_REF] Berger | Life histories, axes of variation in[END_REF].

Pour les populations d'espèces dites rapides, le taux de croissance de la population, généralement noté λ en dynamique des populations, est très fortement influencé par des variations des paramètres reproducteurs (âge à la première reproduction et taux de fécondité) et il est peu impacté par les variations des taux de survie [START_REF] Clobert | Bird population dynamics, management, and conservation: the role of mathematical modeling. Dans: Bird population studies: their relevance to conservation and management[END_REF]. Au contraire, dans les populations d'espèces dites lentes, λ est largement influencé par les taux de survie juvénile et adulte, et peu par des changements dans les paramètres reproducteurs [START_REF] Clobert | Bird population dynamics, management, and conservation: the role of mathematical modeling. Dans: Bird population studies: their relevance to conservation and management[END_REF]. En se basant sur la comparaison de 138 populations de mammifères, [START_REF] Oli | The fast-slow continuum and mammalian life-history patterns: an empirical evaluation[END_REF] a empiriquement confirmé ces prédictions. Les populations d'espèces lentes sont aussi souvent caractérisées par un taux de survie adulte relativement stable entre les années (hypothèse de canalisation des traits ayant le plus d'impact sur la dynamique des populations, [START_REF] Festa-Bianchet | Population dynamics of large herbivores: variable recruitment with constant adult survival[END_REF], mais par une forte variabilité inter-annuelle du taux de survie juvénile, due à des processus stochastiques (ex. des fluctuations environnementales) ou de la densité-dépendance [START_REF] Festa-Bianchet | Population dynamics of large herbivores: variable recruitment with constant adult survival[END_REF].

Les stratégies démographiques des oiseaux marins nicheurs du littoral français

D'après le 5 ème recensement national des oiseaux marins nicheurs en France métropolitaine conduit de 2009 à 2012, il y a 28 espèces d'oiseaux marins nicheurs, dont 24 espèces littorales (Cadiou et al. 2015). Ces espèces se reproduisent, la plupart du temps, en colonie et sont très fidèles à leur site de reproduction (elles reviennent nicher sur la même colonie d'année en année, voire sur leur colonie de naissance). D'un point de vue démographique, les oiseaux marins nicheurs du littoral français sont généralement des espèces présentant une stratégie lente caractérisée par (Tableau 1) :

-une importante longévité, avec un gradient allant de 10 ans à plus de 50 ans ; -une maturité sexuelle différée : allant de 2-3 à 5-6 ans et même jusqu'à 9 ans, avec un âge de recrutement augmentant généralement avec la longévité de l'espèce ; -une fécondité faible à moyenne: 1 seul oeuf pour les espèces plus longévives, mais certaines espèces en produisent plusieurs, 2-3 chez les mouettes et goélands, jusqu'à 6 chez les cormorans ; -des probabilités de survie adultes interannuelles élevées (attendu > 0,90 à > 0,95 selon les espèces).

Bien que toutes ces espèces d'oiseaux marins nicheurs du littoral français soient longévives, on peut les hiérarchiser selon un gradient de stratégies de reproduction associé à un investissement plus ou moins important dans la productivité ou la survie adulte. Un premier groupe d'espèce composé par le fulmar boréal (Fulmarus glacialis Linnaeus, 1761), les puffins, l'océanite tempête (Hydrobates pelagicus Linnaeus, 1758), le fou de Bassan (Morus bassanus Linnaeus, 1758), le guillemot de Troïl (Uria aalge Pontoppidan, 1763), le pingouin torda (Alca torda Linnaeus, 1758) et le macareux moine (Fratercula arctica Linnaeus, 1758), est constitué d'espèces très longévives (20-30 ans), qui se reproduisent pour la première fois tardivement (4 à 9 ans), et pondent un seul oeuf par évènement de reproduction (Tableau 1). Ces espèces investissent donc peu dans la reproduction mais beaucoup plus dans la croissance et la survie. Au contraire les cormorans, les mouettes, les goélands, les sternes et les guifettes se reproduisent plus tôt (2-3 ans en général), vivent généralement moins longtemps (10-20 ans) et produisent entre 2 à 6 oeufs selon les espèces (Tableau 1). Ces espèces investissent donc plus fortement dans la productivité au détriment de la survie, mais restent malgré tous des espèces dans la partie lente du gradient lent-rapide d'histoire de vie à l'échelle de l'ensemble des oiseaux.

Étant donné la différence d'investissement dans la reproduction en fonction des espèces d'oiseaux marins, on peut s'attendre à des différences dans l'influence relative des paramètres démographiques sur la dynamique de population. Nous proposons ici d'illustrer ces différences en étudiant l'impact de la modulation des paramètres démographiques (fécondité, survie juvénile et survie adulte) sur la dynamique des populations de deux espèces nichant sur le littoral français et aux histoires de vie contrastées. Nous avons modélisé la dynamique des populations à l'aide de modèles matriciels populationnels [START_REF] Caswell | Matrix population models[END_REF].

Cas d'une espèce très longévive : le puffin de Scopoli (Calonectris diomedea Scopoli, 1769)

Le puffin de Scopoli est une espèce longévive (20-25 ans) ne produisant qu'un seul oeuf par an, et sa population semble stable en Méditerranée (Cadiou et al. 2015). Le cycle de vie de l'espèce (10 stades) et les paramètres démographiques des populations des archipels de Marseille ont été caractérisés dans différentes études de Capture-Marquage-Recapture [START_REF] Genovart | Varying demographic impacts of different fisheries on three Mediterranean seabird species[END_REF], Courbin et al. 2018a, Fig. 1). Sur le Frioul, le temps de génération estimé est supérieur à 17 ans, le succès reproducteur (nombre de jeune produit par couple reproducteur par an) est de 0,55, le taux de survie annuelle est de 0,45 pour les juvéniles, et de 0,88 pour les immatures et les adultes (Courbin et al. 2018a). Étant donné la faible survie adulte observée (0,88, soit une valeur très inférieure au 0,95 attendue pour une espèce très longévive, e.g. l'albatros à sourcils noirs Thalassarche melanophris Temminck, 1828 ; [START_REF] Jenouvrier | Climate change and functional traits affect population dynamics of a long-lived seabird[END_REF], associée à une faible productivité, on s'attendrait à un déclin très rapide de la population de puffins. Une forte immigration doit avoir lieu pour expliquer le maintien des populations de Marseille (Courbin et al. 2018a), qui pourrait se combiner à une émigration non-négligeable et des captures accidentelles dans les engins de pêche (pouvant expliquer en partie la faible survie apparente).

L'utilisation d'un modèle matriciel populationnel avec différentes valeurs pour les paramètres démographiques indique, comme attendu pour une espèce très longévive, que la dynamique de population du puffin de Scopoli de l'archipel du Frioul est largement pilotée par la survie adulte et très peu influencée par le succès reproducteur (Fig. 2, voir aussi Courbin et al. 2018a pour les analyses de sensitivité et d'élasticité du λ aux variations des paramètres démographiques). Les variations de la survie juvénile n'affectent aussi que peu la trajectoire de la population. L'augmentation de 1 % du taux de croissance de la population ne nécessite qu'une augmentation de 2 % de la survie adulte, mais une augmentation de 23 % de la survie juvénile ou de la fécondité. Le faible potentiel reproducteur du puffin de Scopoli ne lui permet donc pas de compenser une faible survie adulte (Fig. 2B).

Cas d'une espèce moins longévive : le cormoran huppé (Gulosus aristotelis Linnaeus, 1761)

Le cormoran huppé a une espérance de vie moindre (15 ans) et une plus forte fécondité (1 à 6 oeufs par couvée) en comparaison avec le puffin. Il présente une tendance à l'augmentation des effectifs sur les façades Manche-Atlantique et Méditerranée (Cadiou et al. 2015). Le cycle de vie de l'espèce (5 stades), relativement court par rapport à celui du puffin, et les paramètres démographiques de la population de cormorans huppés du Mor Braz en Bretagne ont été caractérisés dans la littérature (Fortin et al. 2003, Genovart et al. 2017, Fig. 3). Le succès reproducteur moyen (fécondité) est de 1,00 (voir 1,50) et présente une forte variabilité inter-annuelle. La survie juvénile est de 0,44, la survie pour la deuxième année est de 0,76, et la survie adulte est de 0,81 [START_REF] Fortin | The demography and ecology of the European shag Phalacrocorax aristotelis in Mor Braz, France[END_REF].

Le modèle matriciel populationnel indique un temps de génération de 7,6 ans pour cette population. L'utilisation de différentes valeurs pour les paramètres démographiques indique, comme attendu pour une espèce relativement moins longévive et avec un plus fort potentiel de productivité, que la dynamique de la population du cormoran huppé du Mor Braz reste de loin pilotée par la survie adulte mais que la survie juvénile et le succès reproducteur ont une contribution plus importante chez cette espèce que chez les espèces plus longévives comme le puffin (Fig. 4). L'augmentation de 1 % du taux de croissance de la population nécessite une augmentation de 1,3 % de la survie adulte, ou une augmentation de 8 % de la survie juvénile ou de 9 % de la fécondité. La Figure 4B 1 données de longévité maximum issues de https://euring.org/data-and-codes/longevity-list.

30

Figure 1 .

 1 Figure 1. Cycle de vie du puffin de Scopoli (recensement en pré-reproduction) utilisé pour

Figure 2 .

 2 Figure 2. (A) Projections déterministes de la dynamique de population de puffins de Scopoli

Figure 3 .

 3 Figure 3. Cycle de vie du cormoran huppé (recensement en pré-reproduction) utilisé pour

Figure 4 .

 4 Figure 4. (A) Projections déterministes de la dynamique de population de cormoran huppé

Figure 1 .

 1 Figure 1.

Figure 2

 2 Figure 2.

  33

Figure 4 .

 4 Figure 4.

  Tableau 1. Traits d'histoires de vie et caractéristiques démographiques des populations d'oiseaux marins nicheurs de France. Des résultats provenant de sites d'études situés hors de France sont également mentionnés.

				-	0,5-1	-	-	0,850 Espagne	Carboneras 2015
	Goéland railleur (Chroicocephalus genei	23, 31	2-3	2	0,660	-	0,700	0,810 France (Côte	Doxa et al. 2013
	Brème, 1839)								méditerranéennes)
	Goéland d'Audouin (Ichthyaetus audouinii	21, 25	2-3	3	0,380	0,622	0,822	0,892 Espagne	Genovart et al. 2017
	Payraudeau, 1826)
	Espèce Goéland cendré (Larus canus Linnaeus, 1758) 10, 34 Longévité (an) moyenne, max 1 Goéland brun (Larus fuscus Linnaeus, 1758) 26, 35 Fulmar boréal (Fulmarus glacialis Linnaeus, 1761) 20-30, 44 Goéland argenté (Larus argentatus 10-20, 35 Pontoppidan, 1763) Goéland leucophée (Larus michahellis 10, 20 Puffin de Scopoli (Calonectris diomedea 20-25 Naumann, JF, 1840) Goéland marin (Larus marinus Linnaeus, 10-20, 33 Scopoli, 1769) Puffin des Anglais (Puffinus puffinus 20-25, 51 1758) Brünnich, 1764)	Couvée (oeuf) 3 3 3 2-3 1 3 1	Recrutement (an) 3 5 9 4 9 5 5 -5 5	Fécondité Survie juvénile (1ère année) 0,543 0,410 0,530 0,820 0,190-0,510 -0,900--1,500 0,419 -0,920 0,798 0,660 0,450 1,600 0,460 0,740 -0,697 -1,139	Survie immature 0,710 0,885 (avant repro) 0,880 -0,260 0,834 0,880 ---	Survie 0,828 Royaume-Uni Région 0,885 Royaume-Uni adulte 0,957 Bretagne/Manche/Mer 0,881 Bretagne du Nord 0,936 Royaume-Uni 0,834 Royaume-Uni 0,880 France (Ile Frioul) 0,900 Pologne 0,953 Bretagne (Molène) 0,870 Royaume-Uni 0,930 Royaume-Uni	Horswill & Robinson 2015, JNCC Source Horswill & Robinson 2015, Nager & O'Hanlon 2016 Pons & Migot 1995, Cadiou et al. 2017 Kerbiriou et al. 2012, Gallien 2019 Horswill & Robinson 2015 Horswill & Robinson 2015 Skorka et al. 2005 Courbin et al. 2018a Linard & Monnat 1990, Horswill & Robinson 2015 Cadiou & Yesou 2006 Horswill & Robinson 2015, Nager & O'Hanlon 2016
	Puffin yelkouan (Puffinus yelkouan Acerbi, Sterne caugek (Thalasseus sandvicensis	20, 23 24, 31	1 1-2	4 3-4	0,500 0,490-	0,180 0,250	0,825 0,830	0,920 France (Ile Riou) 0,940 Atlantique/Bretagne/Ma Courbin et al. 2018a Robinson 2010, Jacob 2021
	1827) Latham, 1787)				0,700				nche/Grande-Bretagne
	Océanite tempête (Hydrobates pelagicus	20, 33	1	4 3	0,530 0,702	-0,358	-0,741	0,820 Bretagne (Molène) 0,898 Royaume-Uni	Cadiou 2001, Matović et al. 2017 Horswill & Robinson 2015
	Linnaeus, 1758) 1813) Sterne de Dougall (Sterna dougallii Montagu,	10, 24	1-2	3	0,510 0,770	--	-0,735	0,735 Espagne (Côte méditerranéenne) Royaume-Uni 0,855 Bretagne/Manche	Tavecchia et al. 2006, Sanz-Aguilar et al. 2009 Ratcliffe et al. 2008, Jacob 2021
	Fou de Bassan (Morus bassanus Linnaeus, 1758) Sterne pierregarin (Sterna hirundo Linnaeus, 1758)	16-20, 37 10, 33	1 2-3	5 5 3 3-4	0,969-0,540 0,700 0,430-0,570 0,764	0,350 0,424 -0,441	0,880 0,829-0,895 0,280 0,850	0,950 0,919 0,900 0,883 Royaume-Uni Bretagne (Sept-Iles) Royaume-Uni Atlantique/Manche Allemagne	Grémillet et al. 2006, Cadiou et al. 2017, Jacob 2021 Mission Migration 2019 Szostek & Becker 2012 Horswill & Robinson 2015 Horswill & Robinson 2015
	Guillemot de Troïl (Uria aalge Pontoppidan, 1763) Sterne naine (Sternula albifrons Pallas, 1764) 21, 24 32, 43	1 1-3	6 3	0,672 0,590 0,255 0,180-	0,560 0,601	0,854 0,578	0,939 Royaume-Uni Bretagne (Cap Fréhel) 0,899 Bretagne/Italie	Horswill & Robinson 2015, JNCC Cadiou et al. 2017 Tavecchia et al. 2006, Jacob 2021
	Pingouin torda (Alca torda Linnaeus, 1758)	20, 42	1	5 2	0,630 0,518	--	--	0,890 Canada 0,800 Royaume-Uni	Lavers 2007, Lavers & Jones 2007 Horswill & Robinson 2015
	Sterne hansel (Gelochelidon nilotica Gmelin,	-, 16	2-3	5 5	0,570 0.230-	0,630 -	0,630 -	0,895 Royaume-Uni -Bretagne	Horswill & Robinson 2015 Schreiber & Burger 2002,Quaintenne 2020
	Macareux moine (Fratercula arctica Linnaeus, 1789)	22, 45	1	6	0,620 0.790	0,780	-	0,890 Iles écossaises	Saether & Bakke 2000, Miles et al. 2015
	1758) Guifette moustac (Chlidonias hybrida Pallas,	10	3-5	5 2	0,617 1,900	0,709 -	0,783 0,540	0,906 Royaume-Uni 0,800 Pologne	Horswill & Robinson 2015 Ledwon et al. 2014
	Grand cormoran (Phalacrocorax carbo 1811)	20, 32	4 (3-6) 2	2,500	0,500	0,870	0,900 Danemark	Henaux et al. 2007, Frederiksen et al. 2011
	Linnaeus, 1758) Guifette noire (Chlidonias niger Linnaeus,	17, 21	3	3 2	1,985 -	0,540 0,440	0,540 0,570	0,868 Royaume-Uni 0,870 USA	Horswill & Robinson 2015 Naugle 2004
	Cormoran huppé (Gulosus aristotelis 1758)	15, 31	4 (1-6) 3	1,000-	0,440	0,760	0,810 Bretagne/Manche	Fortin et al. 2013, Cadiou et Leicher 2020
	Linnaeus, 1761)				1,500
				2	1,303	0,513	0,737	0,858 Royaume-Uni	Horswill & Robinson 2015
	Mouette tridactyle (Rissa tridactyla Linnaeus,	21, 28	1-3	3	0,651-	-	-	0,804 Bretagne	Frederiksen et al. 2005, Cadiou et al. 2017
	1758)				0,750
				4	0,690	0,790	0,854	0,854 Royaume-Uni	Horswill & Robinson 2015
	Mouette rieuse (Chroicocephalus ridibundus	10-15, 33	2-3	2-3	0,605-	0,213	0,184	0,860 France (Forez)	Niel & Lebreton 2005, Péron et al. 2010
	Linnaeus, 1766)				1,732
				2	0,625	-	-	0,825 Royaume-Uni	Horswill & Robinson 2015
	Mouette melanocéphale (Ichthyaetus	15, 22	2-3	2	1,000-	-	-	-	France	Isenmann et al. 2004
	melanocephalus Temminck, 1820)				1,600
										29