INTRODUCTION

Lattice Boltzmann Methods (LBM) are a powerful tool for the simulation of fluid dynamics [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF] . Due to its attractive computational cost 2 , its capacities for massively parallel computing and the ease to deal with complex geometries using multi-level Cartesian grids, these methods have attracted growing interest both in the academic and industrial spheres in the past decade [START_REF] Succi | The lattice boltzmann equation: a new tool for computational fluid-dynamics[END_REF][START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF][START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] .

LBM being initially designed to tackle isothermal flows, extension to thermal flows is today an active topic of investigation in the community. In achieving that goal, the numerical stability of the collision operator, at the heart of the method, used to be a major issue. The single relaxation time Bhatnagar-Gross-Krook model [START_REF] Bhatnagar | A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems[END_REF] , probably the most popular model, lacks stability for shear flows, but more recent models such as multiple relaxation collision [START_REF] Humières | Multiplerelaxation-time lattice Boltzmann models in three dimensions[END_REF] or regularized kernels [START_REF] Malaspinas | Increasing stability and accuracy of the lattice boltzmann scheme: recursivity and regularization[END_REF][START_REF] Jacob | A new hybrid recursive regularised bhatnagargross-krook collision model for lattice boltzmann method-based large eddy simulation[END_REF] significantly improved stability. Another issue is the resolution of energy or temperature equation, which cannot be straightforwardly achieved on low-order lattices [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] . Two main options are available in the literature. The first is the double distribution function (DDF) option that consists in coupling the LBM solver with a second distribution function, whose main order corresponds to either temperature, energy, or enthalpy (see, e.g. [START_REF] Hosseini | Weakly compressible lattice boltzmann simulations of reacting flows with detailed thermo-chemical models[END_REF][START_REF] Saadat | Extended lattice boltzmann model for gas dynamics[END_REF] , for recent studies). A second option is to couple directly LBM with a scalar (temperature, energy, enthalpy), solved in a coupled finite difference solver. This second option was found attractive, as it allows, for a reasonable cost 2 , to include an arbitrary number of additional scalar equations. Following recent successful applications to compressible [START_REF] Feng | A three dimensional lattice model for thermal compressible flow on standard lattices[END_REF][START_REF] Feng | A compressible lattice boltzmann finite volume model for high subsonic and transonic flows on regular lattices[END_REF][START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Guo | An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice[END_REF][START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF] , atmospheric [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Cheylan | Compressible pressure-based latticeboltzmann applied to humid air with phase change[END_REF][START_REF] Jacob | Wind comfort assessment by means of large eddy simulation with lattice boltzmann method in full scale city area[END_REF][START_REF] Merlier | Lattice-boltzmann large-eddy simulation of pollutant dispersion in street canyons including tree planting effects[END_REF][START_REF] Merlier | Lattice-boltzmann large-eddy simulation of pollutant dispersion in complex urban environment with dense gas effect: Model evaluation and flow analysis[END_REF] or reactive flows [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Tayyab | Lattice-boltzmann modelling of a turbulent bluffbody stabilized flame[END_REF] , the second option, often referred to as hybrid LBM is retained for this study. This works aims at investigating their applicability to buoyancy driven flows, such as those encountered in fire-related problems [START_REF] Domino | Predicting large-scale pool fire dynamics using an unsteady flamelet-and large-eddy simulation-based model suite[END_REF] .

The far field of turbulent, axisymmetric, free-plume in a quiescent, unstratified environment, where buoyancy-driven plume exhibits self-similar behavior, play a significant role in various fluid flows of environmental and technological importance, including thermal plumes that arise due to the convective heating on Earth's surface [START_REF] Morton | Buoyant plumes in a moist atmosphere[END_REF][START_REF] Morton | Forced plumes[END_REF][START_REF] Woods | Turbulent plumes in nature[END_REF] , fire protection engineering with problems associated with fire detection, smoke filling rates of indoor spaces, fire venting, fire heating of structural elements of buildings [START_REF] Heskestad | Engineering relations for fire plumes[END_REF] , wildland fires [START_REF] Fromm | Transport of forest fire smoke above the tropopause by supercell convection[END_REF][START_REF] Kahn | Wildfire smoke injection heights: Two perspectives from space[END_REF] , dispersal of volcanic eruptions [START_REF] Pyle | chichester, new york, weinheim, brisbane, singapore, toronto: John wiley & sons. price £85.00 (hard covers). isbn 0 471 93901 3[END_REF][START_REF] Ernst | Sedimentation from turbulent jets and plumes[END_REF][START_REF] Carazzo | On the rise of turbulent plumes: Quantitative effects of variable entrainment for submarine hydrothermal vents, terrestrial and extra terrestrial explosive volcanism[END_REF] , sea ice plumes [START_REF] Skyllingstad | Turbulence beneath sea ice and leads: A coupled sea ice/large-eddy simulation study[END_REF][START_REF] Widell | Salt release from warming sea ice[END_REF] , smoke stacks [START_REF] Briggs | Chimney plumes in neutral and stable surroundings[END_REF] , and cooling tower plume dispersion [START_REF] Fisher | Predicting cooling tower plume dispersion[END_REF] .

Forced plumes or buoyant jets represent a canonical configuration to study such plumes.

They consist by releasing the plumes from a source with some initial momentum. The flow can be decomposed into three regions: (i) a region featuring a jet-like behavior near the source, (ii) a transitional region and (iii) further downstream, the far-field fully-developed buoyancy-driven plume [START_REF] Chen | Vertical buoyant jets: a review of experimental data[END_REF] .

The first far-field plume theories 41-43 assumed a turbulent flow, a point source of buoyancy, the Boussinesq approximation and a dynamic similarity of the mean and turbulent motion at all elevations. Morton at al. [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous sources[END_REF] developed an integral formulation by assuming both "tophat" radial profiles for both velocity and temperature (density) and a point source, and by introducing an entrainment coefficient, α, defined as the ratio of radial velocity at the edge of the plume to the vertical velocity within the plume. Their model predicts correctly the scaling laws for the plume radius, that increases linearly with height, z, as well as for both velocity and temperature rise above the ambient that decay as z -1/3 and z -5/3 , respectively. The weak plume formulation of Morton et al. [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous sources[END_REF] was extended to "strong plumes" by removing both the Boussinesq approximation and the point source assumption through the introduction of a virtual origin and, by considering more realistic Gaussian profiles for both velocity and temperature [START_REF] Heskestad | Engineering relations for fire plumes[END_REF] . This in conjunction with experiments in fire plumes above the flames [START_REF] Heskestad | Dynamics of the fire plume[END_REF] provided expressions for plume radius and centreline velocity and temperature consistent with experimental data [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF][START_REF] Papanicolaou | Investigations of round vertical turbulent buoyant jets[END_REF][START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] . Another important feature of buoyant plumes is that the radial profiles of dimensionless mean velocity and temperature and dimensionless rms turbulent fluctuations of velocity and temperature exhibit a self-similar behaviour with η = r/z as a self-similar variable, where r is the radial coordinate [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF][START_REF] Papanicolaou | Investigations of round vertical turbulent buoyant jets[END_REF][START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF][START_REF] Nakagome | The structure of turbulent diffusion in an axisymmetrical thermal plume[END_REF][START_REF] Dai | Structure of round, fully-developed, buoyant turbulent plumes[END_REF] . Progress was also made in the understanding of the role of buoyancy on the entrainment process and the flow behavior. In particular, Saeed et al. found that buoyancy contributes to enhance the mean kinetic energy budget but also the momentum flux [START_REF] Saeed | Buoyancy-modified entrainment in plumes: Theoretical predictions[END_REF] .

The CFD modeling of forced plumes was also a very active research area. A significant amount of Reynolds-Averaged Navier-Stokes simulations were reported with different formulations of the k-model [START_REF] Nam | Numerical simulation of thermal plumes[END_REF][START_REF] Worthy | Comparison of modified k-epsilon turbulence models for buoyant plumes[END_REF][START_REF] Brescianini | New evaluation of the k-epsilon turbulence model for free buoyant plumes[END_REF][START_REF] Yan | Application of two buoyancy modified k-epsilon turbulence models to different types of buoyant plumes[END_REF] . On the other hand, Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] and Yan [START_REF] Yan | Large eddy simulations of a turbulent thermal plume[END_REF] showed the capability of Large Eddy Simulation (LES) to predict well the plume's puffing, self-preserving and spreading. LES of forced jet were also performed [START_REF] Pant | Evaluation of an energy consistent entrainment model for volumetrically forced jets using large eddy simulations[END_REF] to evaluate the energy-consistent approach for modeling entrainment rate coefficient, α, developed by Kaminski et al. [START_REF] Kaminski | Turbulent entrainment in jets with arbitrary buoyancy[END_REF] and van Reeuwijk and Craske [START_REF] Van Reeuwijk | Energy-consistent entrainment relations for jets and plumes[END_REF] . Direct numerical simulation (DNS) and LES of thermal plumes were also reported [START_REF] Pham | Large-eddy simulation of a pure thermal plume under rotating conditions[END_REF][START_REF] Plourde | Direct numerical simulations of a rapidly expanding thermal plume: structure and entrainment interaction[END_REF][START_REF] Chen | New findings in vorticity dynamics of turbulent buoyant plumes[END_REF][START_REF] Bhaganagar | Numerical investigation of starting turbulent buoyant plumes released in neutral atmosphere[END_REF] . These studies focused mainly on the generation and growth of buoyancy-induced instabilities in the near field that governs the transition from laminar to turbulence. In particular, it was shown that these instabilities have to be fully resolved to capture the dynamics of such purely buoyant thermal plumes [START_REF] Pham | Large-eddy simulation of a pure thermal plume under rotating conditions[END_REF] . In addition, the DNS was found in good agreement with experimental data in the far field [START_REF] Plourde | Direct numerical simulations of a rapidly expanding thermal plume: structure and entrainment interaction[END_REF] . This literature survey reveals that all the numerical investigations were performed by using a low Mach-formulation of the Navier-Stokes equation. To the authors' best knowledge, no attempt to consider lattice-Boltzmann method was reported.

The goal in the current work is to assess the ability of the pressure-based Lattice-Boltzmann method proposed by Farag et al. [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] to correctly predict the behavior of a forced plume in the self-similarity region away from the source. This manuscript is organized as follows. The first Section recalls the target macroscopic equations, as well as the Lattice-Boltzmann algorithm proposed. The second Section presents a number of 2D canonical flow validations, including Rayleigh-Benard and Rayleigh-Taylor instabilities. The third Section presents a large eddy simulation of the plume experimentally described by Shabbir and George [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] , along with comparisons with the large eddy simulation presented by Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] and the theoretical model of Morton [START_REF] Morton | Forced plumes[END_REF][START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous sources[END_REF] . All the simulations were performed using the compressible ProLB code [START_REF]Prolb[END_REF] .

I. LATTICE-BOLTZMANN MODEL FOR COMPRESSIBLE FLOWS

A. Macroscopic governing equations

The flow mass, momentum and energy conservation equations are introduced as:

∂ρ ∂t + ∂ρu i ∂x i = 0 (1)
∂ρu i ∂t + ∂ρu i u j + δ ij p ∂x j = ∂Π ij ∂x j + ρg i , (i = 1, 2, 3) (2)
ρ ∂h ∂t + ρu j ∂h ∂x j = Dp Dt - ∂q j ∂x j + Π ij ∂u i ∂x j , ( 3
)
where ρ is the mass volume, u i is the velocity vector, p is the pressure, g i is the gravitational acceleration, h is the mass enthalpy and δ ij as the Kronecker delta symbol.

Π ij = µ ∂u i ∂x j + ∂u j ∂x i -δ ij 2 3 ∂u k ∂x k , (4)
with µ the dynamic viscosity. Finally, the heat flux q i in the energy equation reads

q i = -λ ∂T ∂x i ( 5
)
where T is the temperature, and λ the heat conductivity, obtained assuming constant Prandtl number:

Pr = c p µ λ = ν α , ( 6
)
where ν is the kinematic viscosity defined as ν = µ/ρ while α is the thermal diffusivity defined as α = λ/ρc p . The system of Eqs (1-3) is fully closed by the choice of an equation of state

p = ρ.r.T, h = c p .T, (7)
with c p the mass heat capacity at constant pressure, and r = 287J kg -1 K -1 is the gas constant for air per unit mass. Note that we assume c p to be constant, since the temperature in the test cases presented does not exceed 600K.

B. Turbulence modeling

For the large eddy simulation presented in Sec. III, a Vreman eddy-viscosity subgrid scale model is used [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF] . The filtered expressions for Eqs. (1-3) are widely reported in the literature (see, e.g. [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] ) and not recalled here. Applying the subgrid-scale model numerically comes down to modifying the viscosity µ through the addition of a turbulent viscosity µ t obtained as:

µ t = ρC B β α ij α ij , (8)
with

α ij = ∂u j ∂x i , (9)
β ij = ∆ 2 m α mi α mj , (10)
B β = β 11 β 22 -β 2 12 + β 11 β 33 -β 2 13 + β 22 β 33 -β 2 23 , (11)
The constant C is related to the Smagorinsky constant C s as C = 2.5C 2 s . ∆ m is the local mesh size. The model is simple to implement and compute as it only requires the local filter width (i.e. mesh size) and the first order derivatives of the velocity field. In our simulation C s = 0.1, following Vreman's recommendation [START_REF] Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF] , which was also adopted in the forced plume large eddy simulation by Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] . The heat flux takes into account the subgrid contribution via a turbulent Prandtl number Pr t .

C. Numerical method: Lattice Boltzmann solver For a complete description of the numerical method, the reader is referred to Farag et al. [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF] . Lattice-Boltzmann methods are derived from a space, time and velocity discretization of the Boltzmann equation [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] . In the present model, the probability density function f i is solved at each point x via the Boltzmann equation discretized as:

f i (x + c i δt, t + δt) = f eq i (x, t) + 1 - δt τ f neq i (x, t) + δt 2 F E i (x, t) (12)
where δt is the time-step, and c i is the i th discrete velocity of the D3Q19 lattice [START_REF] Krüger | The lattice boltzmann method[END_REF] , and F E i is a volume force including gravity and correcting terms, as defined in Appendix A. The equilibrium and off-equilibrium populations (f eq i , f neq i ) are to be defined in Eqs (14, 16).

In Eq. ( 12), the relaxation time τ is related to the dynamic viscosity as:

τ = µ ρc 2 s + δ t 2 , (13)
where c s = δx/( √ 3δt) is the characteristic velocity of the D3Q19 lattice [START_REF] Krüger | The lattice boltzmann method[END_REF] .

The equilibrium function is obtained as:

f eq i = ω i ρθ + H (1) i,α c 2 s ρu α + H (2) i,αβ 2c 4 s ρu α u β + a (3) + a (4) , (14)
where ω i is the D3Q19 weight of discrete velocity c i , H i are the discrete Hermite polynomials, defined in Appendix A, and a [START_REF] Succi | The lattice boltzmann equation: a new tool for computational fluid-dynamics[END_REF] and a [START_REF] Guo | Lattice Boltzmann method and its applications in engineering[END_REF] are third and forth order terms also provided in Appendix A. The reduced temperature θ reads

θ = rT c 2 s . (15)
The off-equilibrium counterpart f neq is obtained as:

f neq i = ω i H (2) i,αβ 2c
= i f i -f eq i + δt 2 F E i H (2) i,αβ . (18)
and the third-order contribution defined in Appendix A.

Finally, the macroscopic variables are reconstructed from f i as:

ρ(t + δ t , x) = i f i (t + δ t , x) -(ρθ)(t, x) + ρ(t, x) (19)
ρu i (t + δ t , x) = i c i f i (t + δ t , x) + δ t 2 F g i . ( 20
)
where F g i is the gravity force term defined in Appendix A. The enthalpy equation ( 3) is solved at the same time using a finite difference discretization under non-conservative form, exactly as presented by Tayyab et al. [START_REF] Tayyab | Experimental and numerical lattice-boltzmann investigation of the darrieus-landau instability[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF] . Second-order consistency to the macroscopic equations (1-3) can be shown via Chapman-Enskog 5 , or Taylor 68,69 expansions.

II. CANONICAL 2D VALIDATIONS

This Section provides canonical validations of the numerical method presented above.

It focuses on two gravity-driven configurations: the Rayleigh-Benard and Rayleigh-Taylor instabilities, since the properties of the flow solver are already validated in a large number of compressible flow configurations in the absence of gravity [START_REF] Boivin | Benchmarking a lattice-boltzmann solver for reactive flows: Is the method worth the effort for combustion?[END_REF][START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: bridging between pressure-based and density-based methods[END_REF][START_REF] Cheylan | Compressible pressure-based latticeboltzmann applied to humid air with phase change[END_REF][START_REF] Tayyab | Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes[END_REF][START_REF] Tayyab | Lattice-boltzmann modelling of a turbulent bluffbody stabilized flame[END_REF][START_REF] Coratger | Large-eddy lattice-boltzmann modelling of transonic flows[END_REF] .

A. Rayleigh-Benard Instability

The Rayleigh-Benard instability is a configuration involving natural convection and heat transfer [START_REF] Bergé | Rayleigh-bénard convection[END_REF] . The Rayleigh number, Ra, describes, on the one hand, the balance between buoyancy versus viscous forces in the momentum equation and, on the other hand, the balance between

Ra = gβ(T H -T C )H 3 αν = Pr gβ(T H -T C )H 3 ν 2 , ( 21
)
where g = 9.81m s -2 is the gravitational acceleration and H = 1 m is the domain size. β is the thermal expansion coefficient that is equal to 1/T ref for an ideal gas with isobaric expansion (i.e. at constant pressure). In the present study, T ref is taken equal to T H .

Setting Pr= 0.71, the viscosity can be deduced from the target Rayleigh numbers (Eq. 21)

of 10 4 , 10 5 and 10 6 , while the thermal conductivity, λ, in the energy equation is obtained from Eq. ( 6).

The domain is discretized with a uniform grid with 256 × 256 cells. The flow is then uniformly initialized as u 0 = 0, T 0 = 300K, ρ 0 = 1.2 kg m -3 , and p = p 0 + ρ 0 gy. The simulation is then carried out until convergence using a time-step, δ t = 6.5 × 10 -6 s. The quantitative agreement is shown in Fig. 3, presenting velocity profiles along the centerlines. The present numerical results are compared with the benchmark solution provided by Ouertatani et al. [START_REF] Ouertatani | Numerical simulation of twodimensional rayleigh-bénard convection in an enclosure[END_REF] . For further validation, the local Nusselt number N u is calculated at the bottom wall as:

Nu = ∂T ∂y y=0 . (22)
It can be observed that both velocity and Nusselt number profiles are in excellent agreement with the reference solutions.

B. Rayleigh-Taylor instability

The Rayleigh-Taylor instability is another classical test case for buoyancy-driven flows due to its practical and fundamental importance. It was investigated extensively in the literature by different numerical methods [START_REF] He | A lattice boltzmann scheme for incompressible multiphase flow and its application in simulation of rayleigh-taylor instability[END_REF][START_REF] Chiappini | Improved lattice boltzmann without parasitic currents for rayleigh-taylor instability[END_REF][START_REF] Hosseini | Lattice boltzmann solver for multiphase flows: Application to high weber and reynolds numbers[END_REF][START_REF] Ding | Diffuse interface model for incompressible two-phase flows with large density ratios[END_REF][START_REF] He | On the three-dimensional rayleigh-taylor instability[END_REF][START_REF] Li | Additional interfacial force in lattice boltzmann models for incompressible multiphase flows[END_REF][START_REF] Tryggvason | Numerical simulations of the rayleigh-taylor instability[END_REF][START_REF] Guermond | A projection fem for variable density incompressible flows[END_REF][START_REF] Lee | A class of high-resolution algorithms for incompressible flows[END_REF][START_REF] Zu | Phase-field-based lattice boltzmann model for incompressible binary fluid systems with density and viscosity contrasts[END_REF] . It consists of two layers of fluids of different densities (ρ H , ρ L ) at rest under gravitational field, as illustrated in Fig. 4. The dynamics of this problem is governed by the Atwood (At) and Reynolds (Re) numbers:

At = ρ H -ρ L ρ H + ρ L , Re = U * L x ν , (23)
where L x is the dimension of the domain in the horizontal direction and U * = √ gL x is a reference velocity.

The investigated configuration was previously studied in the literature [START_REF] He | A lattice boltzmann scheme for incompressible multiphase flow and its application in simulation of rayleigh-taylor instability[END_REF][START_REF] Chiappini | Improved lattice boltzmann without parasitic currents for rayleigh-taylor instability[END_REF][START_REF] Hosseini | Lattice boltzmann solver for multiphase flows: Application to high weber and reynolds numbers[END_REF] with two target Reynolds numbers of 256 and 2048, The domain size is L x × 4L x , discretized with 256 × 1024 (fine mesh) or 128 × 512 grid points (coarse mesh). The heavy (index H) and light (index L) fluids are initially separated by a perturbed interface given by the following equation:

y i (x) = L x 10 cos 2πx L x + 2L x (24)
The fluids initial densities are set to ρ H = 3 kg m -3 and ρ L = 1 kg m -3 , corresponding to At=0.5. The pressure was initialized to account for the gravity field as follows:

p =      p 0 + ρ L gy, 0 ≤ y ≤ y i (x) p 0 + ρ L gy i (x) + ρ H g(y -y i (x)), y i (x) < y (25)
where p 0 is the pressure at y = 0. Finally, the domain size and gravity are set to L x = 0.25m and g = 20m s -2 . The fluid viscosity is obtained from the target Reynolds numbers of 256 and 2048. To investigate the robustness of the method, simulations were carried out on the coarser mesh. An excellent agreement is also obtained with a maximum error less than 2%.

III. LARGE EDDY SIMULATION OF A THERMAL PLUME

This Section presents a large eddy simulation of a buoyant plume, generated by a vertical jet of hot air into a quiescent atmosphere. The source conditions correspond to the experiments of Shabbir and George [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] , summarized in Table I. The plume source diameter, D, the indicates the coarse mesh, (+) for the fine mesh, and (•) for the reference from He et al. [START_REF] He | A lattice boltzmann scheme for incompressible multiphase flow and its application in simulation of rayleigh-taylor instability[END_REF] .

D(m) T a (K) T 0 (K) U 0 (m/s) F 0 (m 4 /s 3 ) M 0 (m 4 /s 3 ) Re = U 0 D ν Fr = U 2 0 gD 0.
based on inflow mean injection velocity, source diameter and kinematic viscosity, is 1273.

The specific momentum, M 0 , buoyancy, F 0 , mass, Q 0 , and the Morton length scale, L M , are defined as:

F 0 = 2πg ∞ 0 U z ∆T T rdr, M 0 = 2π ∞ 0 U 2 z rdr, Q 0 = 2π ∞ 0 U z rdr, L M = M 3/4 0 F 1/2 0 ( 26
)
where r is the radial coordinate. Morton [START_REF] Morton | Forced plumes[END_REF] and Morton and Middleton 83 introduced the source parameter Γ 0 that characterizes the plume as being either lazy (Γ 0 > 1), pure (Γ 0 = 1)

or forced (0 < Γ 0 < 1):

Γ 0 = 5Q 2 0 F 0 8 √ παM 5/2 0 ( 27
)
with α the entrainment coefficient explained in details later on. The value of Γ 0 in our simulation is around 0.9, which indicates a forced plume, having a value of Γ 0 near unity says that the plume is forced but not too much as a result the buoyancy is significant near the source which explains the acceleration zone detailed later.

A. Numerical set-up

The computational domain is a box of size 18D×9D×9D. A uniform mesh, composed of 300×150×150 cells, is considered. The simulation was performed with ProLB on 280 processors. The time-step, based on the sound speed, is δ t = 4.5 × 10 -6 s. In accordance with previous LES of this configuration [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] , the Vreman turbulence model, described in Sec.

I B, is applied with a turbulent Prandtl number of Pr t = 0.3.

The boundary conditions are as follows: at the outlet, a Dirichlet condition is considered for pressure whereas a Neumann condition is applied for other variables with a clip for the axial velocity to prevent any backflow of the plume. Typical inflow/outflow boundary conditions are considered for the vertical sides.

At the inlet, temperature and velocity were imposed to represent a plume source.

For the inlet boundary condition, we followed the strategy of Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF][START_REF] Zhou | Study of density effects in turbulent buoyant jets using large-eddy simulation[END_REF] to ensure a transition from laminar to turbulence at a very short distance of the exit, consistently with the experimental observations of Shabbir and George [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] . This kind of fluctuations work more as perturbations with artificial nature so they are not divergence free. However, this does not represent an important issue because the associated time scale is large compared to the flow turbulent time scales (the fastest time scale of our injection is around 0.2s).

As a consequence, the impact of this synthetic injection vanishes few diameters away from the inlet where we start performing our analysis. It consists in superimposing azimuthal disturbances:

u (r) = AU 0 (r)[(1 - r D ) N n=1 sin(2πf t/n) + r D N n=1 sin(2πf t/n + θ)] (28)
to a mean flow U 0 (r) corresponding to a pipe profile:

U 0 (r) = 1 2 U 0 [1 -tanh(b 2 (2r/D -D/2r))] . (29)
A is the amplitude of the forcing and N = 6 is the number of the modes. f is the frequency of the forcing, that is determined by the jet preferred mode corresponding to a Strouhal number, St = f D/U 0 , of 0.3, leading to f = 4.629Hz. In the mean pipe flow profile, θ is the azimuthal angle and b 2 = 6.25 [START_REF] Michalke | Survey on jet instability theory[END_REF] .

Note that Eq. ( 28) was slightly modified from the original formulation [START_REF] Didden | Unsteady separation in a boundary layer produced by an impinging jet[END_REF][START_REF] Menow | Large-eddy simulations of forced three-dimensional impinging jets[END_REF] , which presented a singularity at the center. The forcing amplitude A = 0.2/ √ 3, corresponds to a RMS fluctuations of 20% for the axial velocity and A = 0.01/ √ 3 corresponds to a RMS fluctuations of 1% for the other two components.

The time-averaged statistics (mean, rms,..) presented hereafter were collected over 15

forcing cycles = 20 s once a statistical steady state was reached. The forcing cycle is defined by the longest period of the sine series in Eq. (28) = 1.3 s.

B. Results and discussion

Qualitative description

Figure 6 illustrates the transition process trough a snapshot of the three dimensional iso-surface for the Q-criterion [START_REF] Jcr | Eddies, streams, and convergence zones in turbulent flows[END_REF] along with temperature and density fields. The Q-criterion is defined as:

Q = 1 2 ( Ω 2 -S 2 ), (30)
where S and Ω are the strain rate and the vorticity tensor, respectively:

Ω = 1 2 (∇u + ∇u T ), (31)
S = 1 2 (∇u -∇u T ), (32)
The Q-criterion defines the areas where the vorticity magnitude is larger than the magnitude of the strain rate, such that Q > 0 indicates the existence of a vortex.The potential core of the plume becomes rapidly turbulent after few diameters from the source which is consistent with the experimental observations of Shabbir and George [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] . The transition occurs due to the growth of azimuthal instabilities that forms large coherent energy containing structures which eventually break down to generate small-scale turbulence.

Figure 7 shows the energy spectrum based on the axial velocity at distance z/D = 4, it

shows the energy cascading reported by the theory of Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF] and that we have the experimentally that the slope change from -5/3 to -3 is due to strong energy feeding as a result of the large plume vortices driven by buoyancy force. This region cannot be clearly identified on the spectrum of temperature fluctuations in Figure 7.

Axial mean quantities

U c = A U z -1/3 F 1/3 0 , T a T c = 1 -A T z -5/3 F 2/3 0 /g, (33)
where A T = 9.4 and A U = 3.4 were fitted from the experimental results (see also in Table II). The centerline velocity in Fig. 8 increases rapidly from its initial value at the inlet to a maximum value of about 1.8 at z/D ≈ 2.5 and then decreases afterwards rapidly to reach values lower than the inflow velocity after about 6 diameters. This behavior was also observed by Lingens et al. [START_REF] Lingens | Instability of buoyant diffusion flames[END_REF] who experimentally investigated buoyant jet diffusion flame. The initial acceleration in the near field is due to the large buoyancy force resulting from large temperature (density) difference between the plume core and the ambient. The rapid deceleration after the peak results from the turbulent mixing of the plume with the surrounding [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous sources[END_REF] who reported that a forced plume will reach a pure plume behavior for z/L M > 5. In addition, our simulation exhibits a good agreement with the experimental profile obtained from the correlation of Shabbir and George 47 (see Eq. 33) in the plume-like region.

Forced plumes becomes plume-like far away from the source in homogeneous environment even if the injected momentum flux is large when the function Γ changes from a value smaller than 1 at the source to a value of 1 in the far field. Four regions were identified in the present simulation: 1) a non-buoyant region where momentum dominates the flow,

2) an acceleration region where the plume is accelerated due to gravity, 3) an intermediate region where influence of initial momentum weakens, and 4) the plume-like region (i.e. selfsimilarity region) where the plume dynamics is dominated solely by the buoyancy forces.

This picture is consistent with the descriptions of Gebhart et al. [START_REF] Gebhart | Buoyancy-induced flows and transport[END_REF] and Chen and Rodi 93 although they did not report the acceleration region (i.e. region 2). Note that the limits of each region in Fig. 9 are defined using the velocity inflection points, consistently with the global behavior of the plume.

Fluctuations quantities

In this section, the axial evolution of the rms values of axial velocity and temperature fluctuations and the cross-correlation between velocity and temperature fluctuations are discussed and compared with experiments [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF][START_REF] Papanicolaou | Investigations of round vertical turbulent buoyant jets[END_REF][START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF][START_REF] Nakagome | The structure of turbulent diffusion in an axisymmetrical thermal plume[END_REF] , for which the fitting parameters are reported in Tab. II.

Figure 10 presents the rms values of axial velocity and temperature fluctuations. As expected, the velocity fluctuations are about 20% at vicinity of the inflow plane and corresponds to the imposed disturbance level. The velocity fluctuations decrease in the potential core region of the plume before starting to increase very abruptly in the laminar to turbu- to turbulence with a realistic energy cascade in the far field where we perform our analysis.

Reference A T A U B T B U (T 2 ) 1/2 /∆T c (u 2 z ) 1/2 /U c u z T /(u 2 z ) 1/2 (T 2 ) 1/2
In the plume-like region, both velocity and temperature fluctuations decrease at a same rate as mean velocity and temperature to ensure constant ratio of u 2 1/2 /U c and T 2 1/2 /(T c -T a ).

The predicted velocity-based turbulence intensity in the plume-like region is lower than those of 0.28 and 0.33 reported by George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] and Shabbir and George [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] , respectively. It is in better agreement with those of 0.25 reported by Papanicolaou and List 46 and Nakagome and Hirata 48 . On the other hand, Fig. 10 shows that the temperature-based turbulence intensity is also consistent with the available data. Figure 10 shows the evolution of the cross correlation between velocity and temperature fluctuations, u T /((u 2 )

1/2 (T 2 ) 1/2 , along the center line. It can be clearly observed that velocity and temperature fields are positively correlated in this type of flows with a predicted nearly constant value in plume-like region of about 0.55. This value is lower that those reported by George and co-workers [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF][START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] in the range 0.6-0.7, averaged to 0.67 (see Table II), and in closer agreement with those of 0.46 and 0.51 reported by Nakagome and Hirata 48 and Papanicolaou and List 46 , respectively.

Self-similarity

An important feature of the mean flow in the fully-developed region of turbulent positively buoyant plumes is the "self-similarity" or "self-preserving" behaviour. The radial mean velocity and temperature profiles follow a Gaussian shape and become wider as the plume rises. The profiles collapse on the same curve when considering appropriate dimensionless variables:

U z U c = exp (-B u r 2 z 2 ) , T -T a T c -T a = exp (-B T r 2 z 2 ) (34)
The coefficients B u and B T are unknown empirical constants that can be obtained by assuming a linear growth of the plume width b 43 :

b z = 6 5 α = const, (35)
The coefficients will then be calculated using B u = z/b u and B T = z/b T , where b u and b T are the plume width defined by the distance from the centerline to the point at which we have 1/e of the centerline values of velocity and temperature, respectively. George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] determined by experiments those coefficients as B u = 55 and B T = 65.

The radial profiles of mean velocity and temperature from our LES at z/D = 10, 12, 14, 16

and the profiles of George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] form are plotted in Fig. 11. The velocity and the temperature rise above the ambient are normalized by the centreline value. The profiles are plotted versus the non-dimensional radial coordinate r/(z + z 0 ) where z 0 is the virtual origin of the plume. Empirical relationships were reported to estimate the location of the virtual origin [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] . As pointed out by Yang 56 , the location of virtual origin predicted in the simulation can be different from that estimated by empirical formula. Indeed, this location is significantly affected by the transition from laminar to turbulent whose the prediction is a difficult task in LES mainly due to its sensitivity to the plume's source inflow condition.

In the present study, the virtual origin was estimated to collapse the radial profiles in the fully developed region to a single dimensionless Gaussian profile following the methodology proposed by Yan 56 , giving z 0 set equal to 2.3D. It can be observed in Fig. 11 that the self-similarity is well preserved in the simulation and the predicted self-similarity profiles agree well with those reported by Georges et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] .

Entrainment

The mechanism of turbulent mixing which brings air into the buoyant plume is called entrainment. The ideal plume theory is based on both Boussinesq and top-hat radial-profile assumptions and assumes that the mean entrained flux across the edge of the plume E (entrainment rate) is proportional to the local upward velocity W . An air entrainement coefficient is then defined as:

α = E bW (36)
where E, W , b are know as the top-hat variables of entrainment rate, local vertical velocity and plume width defined by Turner [START_REF] Turner | Buoyant plumes and thermals[END_REF] :

b 2 W = ∞ 0 U z rdr, b 2 W 2 = ∞ 0 U 2 z rdr, E = d dz ( ∞ 0 U z rdr) (37)
The Plume width, b, can be calculated as the value at which velocity or temperature reaches a value of 1/e of the centerline value as indicated by Morton et al. [START_REF] Morton | Turbulent gravitational convection from maintained and instantaneous sources[END_REF] . This will be referred to as (method 1) hereafter. It can be also obtained from Eq. 37, as b

2 W/ √ b 2 W 2 .
This second method will referred to as (method 2) hereafter. Figure 13 compares the two methods. The experimental slope obtained by George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] and the numerical prediction obtained by Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] are also plotted in Fig. 13. Both the present predicted velocity and temperature half-widths decrease first due to the "necking" process in the near field, as observed experimentally by Cetegen [START_REF] Cetegen | Behavior of naturally unstable and periodically forced axisymmetric buoyant plumes of helium and helium-air mixtures[END_REF] , before, as expected, increasing almost linearly in plume region. The two methods provide on the whole consistent predictions that agree with both the experimental slope and the numerical results obtained by Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] .

The entrainment coefficient, α, can be calculated using Eq. 35. This method requires the knowledge the width of b, along the plume axis. It can be obtained either from temperature and velocity radial profiles (method 1) or from Eq. 37 (method 2), as discussed previously. Another method was adopted by Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] from Eqs. 36 and 37, leading to

α = E/ √ b 2 W 2
. This method will be referred to as (method 3). Figure 13 shows that the three methods provide consistent results in the far-field. Our results agrees well with the LES of Zhou et al. [START_REF] Zhou | Large-eddy simulation of a turbulent forced plume[END_REF] which settles on a constant value in the far-field, α = 0.09 -0.1. Our predictions of α in the plume region are also close to the value of 0.116 adopted by Morton in his plume model [START_REF] Morton | Forced plumes[END_REF] and the experimental value 0.108 reported by George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] . I):

M = 2π ∞ 0 (U 2 z + u 2 z -v 2 )rdr (38)
The moment flux ratio increases with the height according to the following relationship [START_REF] Denman ; Fischer | mixing in inland and coastal waters. academic press, inc., new york, xiv + 483[END_REF] :

M M 0 = k( x L M ) 4/3 (39)
Different values of 0.35 [START_REF] Denman ; Fischer | mixing in inland and coastal waters. academic press, inc., new york, xiv + 483[END_REF] , 0.34 [START_REF] Shabbir | Experiments on a round turbulent buoyant plume[END_REF] and 0.29 [START_REF] Papanicolaou | Investigations of round vertical turbulent buoyant jets[END_REF] were reported for the coefficient k. The second is the buoyancy flux F , that has to be conserved conserved along the plume height:

F = 2πg ∞ 0 (U z ∆T T + u z T T )rdr (40)
The buoyancy flux is normalized by its injection value, F 0 . Figure 14 shows the evolution of F/F 0 along the plume height. The simulated normalized buoyancy flux evolves around unity, consistently with the theory (solid line). It appears clearly by comparing the solid and dashed lines that the turbulent contribution is essential. When it is disregarded, the buoyancy flux decreases with the axial distance and is no more conserved. The turbulence contribution is predicted around 15% -20% as also noted by Shabbir & George 47 , while George et al. [START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF] and Papanicolaou & List 46 found the contribution to be about 15%.

IV. CONCLUDING REMARKS

A recursive regularized pressure based LBM solver (ProLB) was tested for buoyancy driven flows.

The solver was able to correctly validate the Rayleigh Bénard cavity test case for different Rayleigh numbers Ra = 10 4 , 10 5 , 10 6 , vertical and horizontal velocity profiles as well as Nusslet number profiles at the bottom wall were all in a good agreement with the reference.

As for the Rayleigh Taylor test case, we were able to correctly predict the instantaneous evolution of the positions of bubble and spike for two different Reynolds numbers Re = 256, 2048. In addition, the test case was run on a coarser mesh to test the robustness of the solver and the results were abundantly satisfying.

For the 3d forced plume simulation, which is a critical test case in which the buoyancy is highly coupled with momentum and turbulent mixing, the solver was able to anticipate the correct physics of a thermal plume from numerous aspects listed below:

• The velocity energy spectrum follows the Kolmogorov theoretical slope of -5/3 indicating a proper resolution of the turbulence energy cascading as reported in the literature.

• Axial profiles of mean velocity and temperature were in a good agreement with the experimental data.

• Our forced-plume reaches a plume-like region at around z/L m = 4 ≈ 5 which is consistent with the findings in the literature.

• The axial profiles of rms for velocity and temperature also were in a good agreement with the experiments, we should emphasize that we did not take into account the experimental errors which are significant specially for the second order statistics.

• The cross-correlation between velocity and temperature has a high positive value which compares well with the reported values from experiments and indicates a strong coupling between the velocity and temperature fluctuations due to gravity.

• Self similarity profiles in the far field (i.e. plume-like region) were achieved for both mean and rms of velocity and temperature.

• The growth rate of the plume was examined through the spatial evolution of the plume width. The growth of the plume compared very well with experimental and numerical references.

• The entrainment of fluid form the surrounding was correctly predicted by examining the entrainment coefficient α, and the predictions were in a good agreement with the theoretical, experimental and numerical references.

• Integral quantities, mainly Buoyancy flux and momentum flux, were compared with the experiments and both were in a good agreement, we emphasize about the finding that the turbulent heat flux participates by around 20% in the total buoyancy flux which is consistent with the experiments.

From all the previous points we can conclude that our solver is capable of reproducing the physics of a thermal plume correctly whether the mean values, the second order statistics or even integral quantities through the plume, and that our code can handle any type of flows with variable densities regardless of their complexity. = ω i a (0),eq + H where all the differential operations are performed using first order upwind FD except for the divergence operator for which a second order centered FD scheme was employed. The final expression of the forcing term is then

F E i = ω i 2c 4 s H (2)
i,αβ a F E αβ + F g i , (A30)

where F g i is the gravity force term defined as

F g i = ω i ρg α H (
where g α is the gravity acceleration in the direction α.

Figure 1

 1 Figure 1 depicts the configuration to be simulated. It consists of a square box of dimension 1 m × 1 m initially filled with quiescent air, and surrounded by adiabatic walls on the left and right, and isothermal top and bottom walls, resp. at T C = 299.5K and T H = 300.5K.

Figure 1 :

 1 Figure 1: Schematic of Rayleigh Benard test case

Figure 2

 2 Figure2presents temperature contours as well as streamline patterns, showing a good qualitative agreement with the literature (see[START_REF] Ouertatani | Numerical simulation of twodimensional rayleigh-bénard convection in an enclosure[END_REF] , e.g.): for Ra = 10 4 , the flow is symmetric and dominated by the recirculation in the core region with small eddies near the corner.When increasing the Raleigh number, secondary eddies near the top left and bottom right corners appear and become larger.

Figure 2 :

 2 Figure 2: Rayleigh-Benard instability: Temperature contours (top) and streamlines (bottom), for three Rayleigh numbers (10 4 , 10 5 , 10 6 ), from left to right.

Figure 3 :

 3 Figure 3: Rayleigh-Benard instability. u x along the vertical centerline, u y along the horizontal centerline, and Nusselt number along the bottom wall (from left to right), for Ra=10 4 (solid), Ra=10 5 (dashed) and Ra=10 6 (dot-dashed). Symbols indicate the reference data 72.

Figure 4 :

 4 Figure 4: Schematic of Rayleigh Taylor instability.

Figure 5

 5 Figure 5 represents the density contours obtained for the two Reynolds numbers of 256 and 2048 using the finer mesh. The diagrams of the right of the figure represents the time evolution of the bubble and spike positions. Numerical predictions are compared to reference numerical simulation 73 , showing an excellent agreement.

Figure 5 :

 5 Figure 5: Rayleigh-Taylor instability for Re = 256 (top) and Re = 2048 (bottom). Left: Density contours at different normalized times t.U * /L = 1, 2, 3, 4, 5 obtained for the fine mesh. Right: time evolution of the position of both bubble (solid) and spike dashed). ( )

Figure 6 :

 6 Figure 6: Instantaneous 3D Q-criterion colored by velocity magnitude alongside density and temperature fields

Figure 7 :

 7 Figure 7: Temporal energy spectrum at z/D = 4 for axial velocity (left), temperature (right). Dashed lines indicate the expected characteristic slopes.

Figure 8

 8 Figure 8 compares the centerline time-averaged axial velocity U c , and temperature T c , to the experimental data of Shabbir and George 47 , who proposed the following correlation in the plume-like region:

Figure 8 :

 8 Figure 8: Centerline mean axial velocity (left), and Temperature (right) profiles. Solid line for the simulation, symbols for experimental data of Shabbir and George 47

Figure 9 :

 9 Figure 9: Centerline mean non-dimensional axial velocity profile. Solid line for the simulation and symbols for experiment of Shabbir and George 47

Figure 10 :

 10 Figure 10: Center line profiles of r.m.s of axial velocity (left), temperature (center), and the cross-correlation of velocity and temperature fluctuations (right). Solid line for simulation, ( * ) for Shabbir & George 47 , (•) for Geroge et al. 45 , ( ) for Papanicolaou & List 46 , and (+) for Nakagome & Hirata 48 .

Figure 11 :Figure 12 :

 1112 Figure 11: Radial profiles of mean axial velocity (left) ,and mean temperature (right) at four axial positions compared to the experiments of George et al.[START_REF] George | Turbulence measurements in an axisymmetric buoyant plume[END_REF]

Figure 13 :

 13 Figure 13: Evolution with the height of (left) the plume width and (right) the entrainment coefficient α

Fig. 14 compares

 14 our result to these experimental results. Model predictions are in good agreement with the experiments of Shabbir and George 47 and Fisher 96 but overpredict that of Papanicolaou and List 46 .

Figure 14 :

 14 Figure 14: Axial profile of momentum flux ratio (left). The solid line represents the present LES, whereas the symbols represent the experimental data of Fischer et al. 96 (•), (+) for Shabbir & George 47 and ( ) for Papanicolaou & List 46 . Axial profile of buoyancy flux ratio (right), solid line includes the turbulent heat flux while dash line does not.

D

  ≡c iα c iβ c iγ c iδ + c 4 s (δ αβ δ γδ + δ βγ δ δα + δ δα δ βγ ) -c 2 s (c iα c iβ δ γδ + c iβ c iγ δ δα + c iγ c iδ δ αβ + c iδ c iα δ βγ + c iγ c iα δ βδ + c iβ c iδ δ αγ ) . (A13)Any distribution function in the D3Q19r lattice can be written as a (weighted) sum of the contributions from each base polynomial. For instance, the equilibrium distribution in equation (14) reads f eq,19r i

∂

  on the order of the Gauss-Hermite quadrature[START_REF] Krüger | The lattice boltzmann method[END_REF] used in the LB model, an adequate forcing term should be added to achieve a correct viscous stress tensor:a neq αβ ≈ -Π αβ = -µ( ∂u α ∂x β + ∂u β ∂x α -2 D ∂u γ ∂x γ δ αβ ),(A26)with D the spatial dimension. For the D3Q19r basis, the projected forcing term reads asa F E αβ = c 2 s u α αβ + a F D αβ (A27)where a cor αβ is a correction tensor due to the deflection of second order moments of the population introduced by the modification of the mass equation, which can be evaluated asa cor αβ ≡ c 2 s δ αβ ∂(ρ(1 -θ)) ∂t ,(A28)which can be discretized using a backward Euler operator and a F D αβ the correction tensor due to the defect of the lattice at third ordera F D αβ = -x ) ,x (ρu x u y u z ) ,z (ρu x u y u z ) ,y(ρu x u y u z ) ,z (ρu 3 y ) ,y (ρu x u y u z ) ,x(ρu x u y u z ) ,y (ρu x u y u z ) ,x (ρu 3 z )

Table I

 I

: Source parameters of the plume exit mean velocity, U 0 , the hot air temperature, T 0 , and the ambient air temperature, T a , are 6.35 cm, 0.98 m/s, 568 K and 300 K, respectively. The corresponding Reynolds number, Re,

Table II :

 II Summary of mean flow parameters and turbulence intensities for different experiments lence transition region, the initial drop in velocity fluctuations is due to the artificial nature of the fluctuations imposed at the inlet. These artificial fluctuations, without a proper cascade, are dissipated very quickly; however, they constitute the seed for a correct transition

	Shabbir & George 47	9.4	3.4	68	58	0.4	0.33	0.67
	George et al. 45	9.1	3.4	65	55	0.38	0.28	0.67
	Papanicolaou & List 46 14.28 3.85	80	90	0.42	0.25	0.51
	Nakagome & Hirata 48	11.5	3.89 48.1 63	0.36	0.25	0.46

ACKNOWLEDGEMENTS

Centre de Calcul Intensif d'Aix-Marseille and GENCI-TGCC/CINES (Grants 2021-A0092B11951, 2021-A0092A07679) are acknowledged for granting access to their high performance computing resources.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A: Expressions for the LBM solver

In regularized LBM, distribution functions will be constructed using an orthogonal polynomial basis. The basis of the D3Q19r lattice used in the current study consists of 19 polynomials, read 17

i,xxy -H

i,yzz , (A5)

i,zyx

i,xzy

where

where a (0),eq = ρθ , a [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF],eq α = ρu α , a

(2),eq αβ = ρu α u β (A15) a (3r),eq 1 = 3(ρu x u x u y + ρu y u z u z ) , (A16) a (3r),eq 2 = 3(ρu x u z u z + ρu x u y u y ) , (A17) a (3r),eq 3 = 3(ρu y u y u z + ρu x u x u z ) , (A18)

It is worth noting that the the forth-order coefficients (a (4),eq ) are added to improve the isotropicity of the lattice, which could be quite important considering the round jet simulation in the current study.

The third-order off-equilibrium terms are reconstructed recursively from the second-order non-equilibrium tensor as a

(3),neq αβγ = u α a

(2),neq βγ + u β a (2),neq αγ + u γ a

(2),neq αβ


Chaofeng Wang

Haibo Li

Yi Dong

Li Zhao

Gérald Grégori

Yuan Uan Zhao
email: yuanzhao@qdio.ac.cn

Wuchang Zhang
email: wuchangzhang@qdio.ac.cn

Tian Xiao

Gérald Gregori

Planktonic ciliate trait structure variation over Yap, Mariana, and Caroline seamounts in the tropical western Pacific Ocean

Keywords: planktonic ciliate, upwelling, seamount, vertical distribution, Western Pacifi c Ocean

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

INTRODUCTION

Planktonic ciliates (Ciliophora, Spirotrichea, Oligotrichia, and Choreotrichia) are small (size range of 10-200 μm) protists with cilia around the body, including aloricate ciliates and tintinnids [START_REF] Lynn | The Ciliated Protozoa: Characterization, Classifi cation, and Guide to the Literature[END_REF][START_REF] Zhang | An Illustrated Guide to Contemporary Tintinnids in the World[END_REF][START_REF] Zhang | An Illustrated Guide to Marine Planktonic Aloricate Oligotrich Ciliates[END_REF]. Planktonic ciliates are ubiquitous and important consumers of pico-(0.2-2 μm) and nano-(2-20 μm) sized phytoplankton, and are important food sources for metazoan and fi sh larvae [START_REF] Stoecker | Grazing by the jellyfi sh, Aurelia aurita , on microzooplankton[END_REF][START_REF] Dolan | Planktonic ciliates in the Mediterranean Sea: longitudinal trends[END_REF][START_REF] Gómez | Trends on the distribution of ciliates in the open Pacifi c Ocean[END_REF]. As the dominant component of the microzooplankton, the marine planktonic ciliate compartment is a key medium through microbial food web to classical food chain, which plays an important role in material circulation and energy fl ow in the marine ecosystem [START_REF] Azam | The ecological role of water-column microbes in the sea[END_REF][START_REF] Pierce | Ecology of planktonic ciliates in marine food webs[END_REF][START_REF] Calbet | The ciliate-copepod link in marine ecosystems[END_REF].

Trait based approaches was proposed to describe plankton communities in a simple manner [START_REF] Litchman | Trait-based approaches to zooplankton communities[END_REF][START_REF] Kiørboe | A trait-based approach to ocean ecology[END_REF]. The existence of lorica provide some protection as well as additional weight to drag for tintinnids [START_REF] Capriulo | Evolution of the lorica in tintinnids: a possible selective advantage[END_REF]. Tintinnids had a higher proportion in total ciliates in turbulent waters [START_REF] Yu | Diff erences in the vertical distribution and response to freshwater discharge between aloricate ciliates and tintinnids in the East China Sea[END_REF] and waters with higher chlorophyll concentration [START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF]. Size classes of ciliate are important trait as in other plankton [START_REF] Pomerleau | Evaluation of functional trait diversity for marine zooplankton communities in the northeast subarctic Pacifi c Ocean[END_REF][START_REF] Brun | Trait biogeography of marine copepods-an analysis across scales[END_REF]. Arctic waters had a higher proportion of larger cells of aloricate ciliates and of tintinnids with larger lorica opening diameter (LOD) [START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF]. LOD size classes were also used to estimate tintinnid functional redundancy [START_REF] Dolan | Declines in both redundant and trace species characterize the latitudinal diversity gradient in tintinnid ciliates[END_REF]. Plankton in diff erent depth experiences diff erent temperature, food items, and metabolic rate and, therefore, shows diff erent functional traits [START_REF] Prowe | Biogeography of zooplankton feeding strategy[END_REF][START_REF] Teuber | Who is who in the tropical Atlantic? Functional traits, ecophysiological adaptations and life strategies in tropical calanoid copepods[END_REF]. Several types of vertical distribution pattern were reported for tintinnids [START_REF] Kršinić | On vertical distribution of tintinnines (Ciliata, Oligotrichida, Tintinnina) in the open waters of the South Adriatic[END_REF]Wang et al., 2019a).

In the Western Pacifi c, the trait structures of aloricate size fraction, tintinnid proportion, tintinnid vertical distribution groups, and LOD size-classes were studied by Wang et al. (2019a[START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF]. However, we still know very little about their variations in diff erent hydrographic situations. Previous studies showed that seamount with abrupt structure could alter hydrological conditions at surrounding stations than that off sea areas, such as upwelling [START_REF] Rogers | The biology of seamounts[END_REF][START_REF] Rogers | The biology of seamounts: 25 Years on[END_REF][START_REF] Ma | Environmental characteristics in three seamount areas of the tropical Western Pacifi c Ocean: focusing on nutrients[END_REF]. The response of planktonic ciliate trait structure to these hydrological variations needs to be elaborated. In this paper, we present the variation of planktonic ciliates trait structure in relation to the hydrological characteristics such as the depth of deep Chlorophyll a maximum (DCM) layers, surface temperature, and upwelling during three cruises in the tropical Western Pacifi c Ocean.

MATERIAL AND METHOD

Sampling was conducted during three cruises performed from December 3, 2014 to January 8, 2015 (Yap seamount, transects A and B, winter), March 4 to April 3, 2016 (Mariana seamount, transects C and D, spring) and August 3 to September 7, 2017 (Caroline seamount, transects E and F, summer) aboard R/V Kexue ( Science in Chinese) (Fig. 1) in the Western Pacifi c Ocean where depth was about 5 000 m. The depths at the top of the three seamounts were 300, 30, and 50 m over Yap, Mariana, and Caroline seamounts, respectively. For each seamount, the sampling stations were aligned into two approximately perpendicular transects (Fig. 1).

For each station, vertical profi les of temperature, salinity (conductivity), and chlorophyll a in vivo fl uorescence were obtained from the surface down to 200-m depth (or bottom with depth <200 m) using a conductivity-temperature-pressure sensor (CTD; Sea-Bird Electronics, Bellevue, WA, USA). Water samples were collected at 3 to 7 depths using 12-L Niskin bottles attached to a rosette. The depths sampled were 3 m (surface), 30 m, 50 m (not for Yap seamount), 75 m, 100 m, 150 m, and 200 m (or bottom with depth <200 m). Around the DCM, the sampling depths were adapted to sample the DCM if it was within 10 m of any neighbor sampling depth. All water samples (1 L) were fi xed with 1% acid Lugol's iodine and stored at <4 °C in the dark.

Each water sample was subsequently concentrated to about 100 mL by gently siphoning out supernatant water after at least a 48-h sedimentation in the laboratory. The settling and siphoning processes were repeated to concentrate each sample to a fi nal volume of 20 mL (concentrated sample), then settled in an Utermöhl counting chambers [START_REF] Utermöhl | Zur vervollkommnung der quantitativen phytoplankton-Methodik: mit 1 Tabelle und 15 abbildungen im Text und auf 1 Tafel: methods of collecting plankton for various purposes are discussed[END_REF] for at least 24 h.

For each concentrated sample, the entire volume was examined using an Olympus IX 71 inverted microscope (100× or 400×). The statistical method of counting according to [START_REF] Lund | The inverted microscope method of estimating algal numbers and the statistical basis of estimations by counting[END_REF], and the counting of 10-15 individuals in a sample volume of 1 L implies a counting error higher than 50% [START_REF] Lund | The inverted microscope method of estimating algal numbers and the statistical basis of estimations by counting[END_REF]. Some loricae might have been empty when sampled [START_REF] Kato | Tintinnid ciliates as indicator species of diff erent water masses in the western North Pacifi c Polar Front[END_REF][START_REF] Dolan | Observations of apparent lorica variability in Salpingacantha (Ciliophora: tintinnida) in the Northern Pacifi c and Arctic Oceans[END_REF]. Because mechanic and chemical disturbance associated with collection and fi xation procedures could provoke detachment of the protoplasma from the loricae [START_REF] Paranjape | Cultivation of marine pelagic protozoa[END_REF][START_REF] Alder | Tintinnoinea[END_REF], empty loricae of tintinnids were considered as living cells in our study.

For each species, size (length, width according to the shape) of the cells (for aloricate ciliates) or of the loricae (for tintinnids, especially length and oral diameter) were measured for at least 20 individuals when possible. According to lorica morphology and size, tintinnids were identifi ed to the species level according to references [START_REF] Kofoid | A conspectus of the marine and fresh-water Ciliata belonging to the suborder Tintinnoinea, with descriptions of new species principally from the Agassiz Expedition to the Eastern Tropical Pacifi c, 1904-1905[END_REF][START_REF] Lynn | The Ciliated Protozoa: Characterization, Classifi cation, and Guide to the Literature[END_REF][START_REF] Zhang | An Illustrated Guide to Contemporary Tintinnids in the World[END_REF]Wang et al., 2019a). Ciliate volumes were estimated by using appropriate geometric shapes (cone, ball, and cylinder). Tintinnid carbon biomass was estimated using the equation: C = V l (μm 3 , lorica volume)×0.053+444.5 [START_REF] Verity | Relationships between lorica volume, carbon, nitrogen, and ATP content of tintinnids in Narragansett Bay[END_REF]. The conversion factor of carbon biomass for aloricate ciliates used in this study was 0.19 pg C/μm 3 as defi ned by [START_REF] Putt | An experimentally determined carbon: volume ratio for marine "oligotrichous" ciliates from estuarine and coastal waters[END_REF].

The size of the aloricate ciliates was divided into fractions of 10 μm according to the longest length, according to [START_REF] Lessard | Distribution, abundance and size composition of heterotrophic dinofl agellates and ciliates in the Sargasso Sea near Bermuda[END_REF], [START_REF] Sohrin | Full-depth profi les of prokaryotes, heterotrophic nanofl agellates, and ciliates along a transect from the equatorial to the subarctic central Pacifi c Ocean[END_REF][START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF]. The size-fractions were further clustered into small (10-20 μm), medium (20-30 μm), and large (>30 μm) [START_REF] Sohrin | Full-depth profi les of prokaryotes, heterotrophic nanofl agellates, and ciliates along a transect from the equatorial to the subarctic central Pacifi c Ocean[END_REF]. The vertical distribution of tintinnid species was defi ned according to [START_REF] Kršinić | On vertical distribution of tintinnines (Ciliata, Oligotrichida, Tintinnina) in the open waters of the South Adriatic[END_REF]. According to the sampling depth, four groups were classifi ed: group I (species in waters from 0-to 100-m depth), Ⅱ (species in waters from 50-to 200-m depth), Ⅲ (species in waters deeper than 100 m) and Ⅳ (species throughout the water column, which belong to group V in [START_REF] Kršinić | On vertical distribution of tintinnines (Ciliata, Oligotrichida, Tintinnina) in the open waters of the South Adriatic[END_REF]). The biogeography of tintinnid genera (e.g. neritic, cosmopolitan, and warm water types) were derived from [START_REF] Dolan | The Biology and Ecology of Tintinnid Ciliates: Models for Marine Plankton[END_REF]. The tintinnids LOD were divided into diff erent size classes, which are 4-μm apart (12-16 μm, 16-20 μm, and so on) according to [START_REF] Dolan | Declines in both redundant and trace species characterize the latitudinal diversity gradient in tintinnid ciliates[END_REF]. Redundant tintinnid species and proportion of redundant species calculation also followed [START_REF] Dolan | Declines in both redundant and trace species characterize the latitudinal diversity gradient in tintinnid ciliates[END_REF].

RESULT

Hydrography

Hydrological characteristics over the three seamounts were diff erent (Fig. 2). Horizontally, surface waters over Mariana seamount had a lower temperature (27.7-28.6 °C, average 28.1±0.2 °C) and a higher salinity (34.5-34.7, average 34.7±0.1) than surface waters over Yap (temperature 28.9-29.3 °C, average 29.1±0.1 °C, salinity 33.2-34.0, average 33.8±0.2) and Caroline (temperature 29.8-30.1 °C, average 30.5±0.2 °C, salinity 33.0-33.8, average 33.6±0.2) seamounts. The depths of 28 °C isotherm appeared at upper 100 m, 30 m, and 50 m over the Caroline, Mariana, and Yap, respectively. The chlorophyll a (Chl a ) in vivo fl uorescence values were less than 0.1 at surface layers over the three seamounts (Fig. 2). Water column was stratifi ed according to temperature, salinity, and Chl a vertical distribution (Fig. 2). There was a weak upwelling in Mariana seamount at the top 50 m with temperature higher than 27.7 °C, salinity lower than 34.7. The halocline and DCM layers were in about the same depths over three seamounts (Fig. 2). Compared to the DCM depths, Caroline (ranged from 120-150 m, average 126.7±8.3 m) and Mariana (80-110 m, average 107.5±5.6 m) seamounts were deeper than Yap (60-90 m, average 76.1±3.6 m) seamount (Fig. 2). The depth of weak upwelling around the top of Mariana seamount was shallower than 50 m, which had no eff ect on the DCM depths (average 107.5±5.6 m).

Ciliate abundance and biomass

The vertical profi les of ciliate average abundance showed bimodal (in the surface and at the DCM layers) patterns over Yap, Mariana, and Caroline seamounts . However, there were some 1 8 9 9 8 7 6 5 0 1 2 3 4 1 6 1 5 1 4 1 3 1 0 1 1 1 2 1 6 1 5 1 4 1 3 1 2 1 1 0 1 7 1 8 1 9 2 0 2 1 1 7 1 8 1 9 2 0 2 1 1 diff erences in some details.

) S t. 1 2 2 2 1 1 1 8 9 9 8 7 6 5 0 1 2 3 4 1 6 1 5 1 4 1 3 1 0 1 1 1 2 1 6 1 5 1 4 1 3 1 2 1 1 0 1 7 1 8 1 9 2 0 2 1 1 7 1 8 1 9 2 0 2 1 1 0 3 3 3 0 0 4
Ranges of ciliate abundance and biomass at surface (130-323 inds./L, 0.24-0.57 μg C/L) and DCM layers (130-331 inds./L, 0.17-0.64 μg C/L) over Yap seamount were lower than Mariana (surface 103-384 inds./L, 0.11-0.87 μg C/L, DCM 161-405 inds./L, 0.17-0.72 μg C/L) and Caroline (surface 174-380 inds./L, 0.23-1.06 μg C/L, DCM 141-405 inds./L, 0.14-0.77 μg C/L) seamounts. The integrated abundance (the upper 200 m or bottom when depth <200 m) over Yap ((71.7-196.3)×10 3 inds./m 2 , average (147.1±32.7)×10 3 inds./m 2 ) was lower than Mariana ((109.6-254.6)×10 3 inds./m 2 , average (163.3±36.4)×10 3 inds./m 2 ) and Caroline ((150.5- 238.7)×10 3 inds./m 2 , average (195.3±25.7)×10 3 inds./ m 2 ) seamounts (Table 1), but the integrated biomass over Yap (0.18-0.33 mg C/m 2 , average 0.26±0.05 mg C/m 2 ) seamount was similar to Mariana (0.17-0.41 mg C/m 2 , average 0.26±0.06 mg C/m 2 ) seamount, which both lower than Caroline (0.23-0.41 mg C/m 2 , average 0.33±0.05 mg C/m 2 ) seamount (Table 1).

Vertically, the highest tintinnid abundance appeared at DCM layers over Yap and Mariana seamounts. While over Caroline seamount with the deepest DCM (average 126.7±8.3 m), the highest tintinnid abundance occurred at 100-m depth instead of the DCM layer. Yap seamount had the lowest value of tintinnid average abundance proportion to total ciliates (average 4.8%±2.9%) compared to Mariana (average 6.6%±3.9%) and Caroline (average 6.4%±3.1%) seamounts. Although the highest tintinnid abundance proportion to total ciliates appeared at around 100-m depth over all three seamounts, this value over Mariana (12.0%) and Caroline (11.5%) seamounts were about twice higher than Yap (6.4%) seamount (Fig. 4).

Aloricate ciliate size-fraction

Abundance proportion of diff erent size-fractions of aloricate ciliates from surface to 200-m depth revealed that individuals of 10-20-μm size-fraction occupied more percentage (larger than 30% at each depth) over the three seamounts (Fig. 5). Vertically, the abundance proportion of aloricate ciliates in the small sizefraction (10-20 μm) increased from surface down to 200-m depth, and the medium size-fraction (20-30 μm) abundance proportion was relatively constant throughout the water column over the three seamounts (Fig. 5). Despite the above-mentioned similarities, there were also some dissimilarities. First, at the surface layer of Yap seamount, the large size-fraction (>30 μm) aloricate ciliate abundance proportion (44.4%) were higher than that of Mariana (32.8%) and Caroline (36.1%) seamounts. Second, at surface and 30-m depth of Yap seamount area, the large sizefraction (>30 μm) abundance proportion was larger than the small size-fraction (10-20 μm). While over both Mariana and Caroline seamounts, the large sizefraction (>30 μm) abundance proportion was lower than the small size-fraction (10-20 μm) (Fig. 5).

Tintinnid assemblage

Tintinnid species composition

Range of tintinnid abundance over Yap (0-26 inds./L) seamount was lower than Mariana (0-33 inds./L) and Caroline (1-76 inds./L) seamounts. Totally, 89 tintinnid species from 32 genera were found in all stations over three seamounts (Fig. 6; Supplementary Fig. S1; Table 2). Among them, 64 tintinnid species from 28 genera, 73 tintinnid species

Ⅲ Ⅱ C E. tenuis I C A. lata Ⅲ Ⅳ Ⅳ C E. tubulosus Ⅳ Ⅳ Ⅳ C A. minutissima Ⅳ Ⅲ Ⅱ C Metacylis conica Ⅱ Ⅱ C Albatrossiella agassizi Ⅲ Ⅲ Ⅲ - M. sanyahensis I I I C Amphorellopsis acantharus Ⅱ Ⅲ C Ormosella apsteini Ⅲ Ⅲ Ⅲ - A. acuta I I C O. bresslaui Ⅲ Ⅲ Ⅲ - A. amphora I C O. schmidti Ⅲ - A. gaarderae I C O. schweyeri Ⅲ - A. laackmanni Ⅲ Ⅲ C Parundella aculeata Ⅱ Ⅱ Ⅱ C A. minor I Ⅳ C P. diffi cilis Ⅲ C A. quadrilineata I I I C P. infl ata Ⅲ Ⅱ C Ascampbelliella armilla Ⅳ Ⅳ Ⅳ W Petalotricha major Ⅱ Ⅱ W A. retusa I I I W Poroecus apiculatus I Ⅱ Ⅱ W Brandtiella palliata Ⅱ Ⅱ Ⅱ W P. curtus I Ⅳ W Canthafi ella brevis Ⅳ Ⅳ W Proplectella claparedei Ⅳ Ⅳ Ⅳ W C. pyramidata Ⅲ Ⅲ Ⅱ W P. parva Ⅱ Ⅱ Ⅱ W C. truncata Ⅲ Ⅲ W P. perpusilla I Ⅳ Ⅳ W Climacocylis scalaria I I W P. urna Ⅳ Ⅱ W Codonaria oceanica Ⅱ Ⅳ Ⅱ W Protorhabdonella curta Ⅳ Ⅳ Ⅳ C Codonella amphorella Ⅲ Ⅱ C P. simplex I I Ⅳ C C. cuspidata Ⅱ Ⅱ C Rhabdonella amor I Ⅳ Ⅳ W C. galea Ⅱ C R. conica I I W C. tropica I C R. cornucopia I W Codonellopsis meridionalis Ⅱ Ⅲ Ⅲ C R. elegans I I I W C. robusta I Ⅳ C R. exilis I I I W Coxliella laciniosa Ⅱ Ⅱ - R. sanyahensis I I I W Cyttarocylis acutiformis Ⅱ Ⅱ W Rhabdonellopsis apophysata I I W C. eucecryphalus Ⅳ Ⅱ Ⅱ W Salpingella acuminata Ⅳ Ⅳ Ⅳ C Dadayiella cuspis Ⅱ C S. curta Ⅱ Ⅱ Ⅳ C D. ganymedes Ⅳ Ⅳ Ⅳ C S. decurtata Ⅱ Ⅱ Ⅱ C D. pachytoecus Ⅲ C S. faurei Ⅳ Ⅳ Ⅳ C D. striata Ⅲ C S. lineata Ⅱ Ⅱ Ⅱ C Dictyocysta pacifi ca Ⅱ C S. rotundata Ⅱ Ⅱ C D. reticulata Ⅳ I I C Steenstrupiella gracilis I I I C D. spinosa I Ⅲ C S. intumescens Ⅱ Ⅱ Ⅳ C
Epiplocylis acuminata

I I I C S. robusta I I I C E. constricta I Ⅳ Ⅱ W S. steenstrupi I I I C E. undella I Ⅱ W Undella hemispherica Ⅱ W Epiplocyloides reticulata I Ⅳ Ⅳ W U. turgida Ⅳ I Ⅳ W

Eutintinnus apertus

I I I C Xystonella treforti Ⅳ Ⅱ W E. fraknoi I I I C Xystonellopsis brandti Ⅱ Ⅱ Ⅱ W E. hasleae Ⅳ Ⅳ Ⅳ C X. cyclas Ⅲ W E. lusus-undae I I Ⅳ C X. favata Ⅲ W E. pacifi cus I I I C X. heros Ⅱ Ⅲ W E. stramentus I I I C X. paradoxa Ⅱ Ⅱ Ⅱ W
BDP: biogeographic distribution patterns according to [START_REF] Dolan | The Biology and Ecology of Tintinnid Ciliates: Models for Marine Plankton[END_REF]; I: Group I, only appeared in 0-100 m; Ⅱ: Group Ⅱ, only appeared in 50-200 m;

Ⅲ: Group Ⅲ, only appeared in >100 m; Ⅳ: Group Ⅳ, appeared in every depth; -: unclassifi ed in BDP; C: cosmopolitan; W: warm water.

from 31 genera and 76 tintinnid species from 33 genera were recorded over Yap, Mariana, and Caroline seamounts, respectively (Table 2). Fifty tintinnid species appeared in all three seamounts areas. Salpingella faurei was the top dominant species over Mariana and Caroline seamounts. Abundance over Mariana (ranged from 0-11 inds./L) and Caroline (ranged from 0-5 inds./L) seamounts were higher than Yap (ranged from 0-3 inds./L) seamount.

Tintinnid vertical distribution

Based on the vertical distribution of abundance, tintinnid species were divided into four groups over every seamount (Table 2). Groups I and Ⅳ were dominant component from surface to 100 m contributing to more than 60% of the total tintinnid abundance proportion over all three seamounts. Group I species abundance was high in surface, and then decreased to 100-m depth. Group Ⅳ species abundance increased from surface to DCM (Yap and Mariana seamounts) or 100 m (Caroline seamount), then decreased to 200-m depth (Fig. 7). The group Ⅳ abundance and abundance proportion at surface layer over Mariana (3.9 inds./L; 36.9%) seamount were higher than Yap (1.6 inds./L; 14.9%) and Caroline (3.2 inds./L; 28.7%) seamounts (Fig. 7). At 30-m depth, group Ⅳ species abundance was higher than group I over Mariana seamount (Fig. 7), while over Yap and Caroline seamounts, group Ⅳ species abundance was lower than group I.

Lorica oral diameter (LOD) size-classes and redundancy of tintinnid groups

The number of size-classes for the tintinnid species LOD over Yap (13), Mariana (15), and Caroline (16) seamounts were similar (Fig. 8). The LOD size-class with a maximum species richness (28-32 μm) was also same over the three seamounts. Although the number of redundant species over Caroline seamount (61 species) was higher than Yap (50 species) and Mariana (59 species) seamounts, the proportion of redundant species was similar over Yap (79.4%), Mariana (79.7%), and Caroline (79.2%) seamounts.

The LOD size-class with the highest abundance proportion was 24-28 μm over the three seamounts (Fig. 8). Steenstrupiella gracilis , Canthariella brevis , and Ascampbelliella armilla were the main contributors to the higher abundance proportion of the 24-28-μm LOD size-class for the three seamounts (Fig. 6).

The abundance proportion of 12-16-μm LOD sizeclass was much higher over Mariana and Caroline seamounts compared to Yap seamount. Salpingella faurei was the main contributor of the 12-16-μm LOD size-class (Figs.6 & 8). Average maximum abundance of Salpingella faurei in DCM layer over Caroline (2.1 inds./L) seamount or at 150-m depth for Mariana (3.3 inds./L) seamount were higher than at Yap (0.9 inds./L, 100-m depth) seamount.

Correlation between ciliate abundance and environmental factors

Diff erent ciliate group (aloricate ciliates, tintinnids, and total ciliates) abundances had diff erent correlation with environment factors (depth, temperature, salinity, and Chl a ) (Table 3). The total ciliate abundance had a positive correlation with temperature and Chl a in vivo fl uorescence, but negative correlation with depth over the three seamounts (Table 3). Ciliate abundance had negative correlation with salinity over Yap and Caroline seamounts, but a positive correlation over Mariana seamount. Aloricate ciliate and tintinnid abundance had signifi cant positive correlation with temperature over the three seamounts (Table 3). Over Mariana seamount, aloricate ciliate and tintinnid had signifi cant positive correlation with salinity, while it changed into a signifi cant negative correlation over Caroline seamount (Table 3). Tintinnid abundance over Yap, Mariana, and Caroline seamounts had signifi cant positive correlation with Chl a (Table 3). The ciliate vertical distribution patterns were the same over three seamounts with abundance peaks in surface and DCM layers. The food items of the ciliate might be the reason for the two groups: Prochlorococcus and picoeukaryotes had maximum abundance in DCM, Synechococcus and heterotrophic bacteria had maximum abundance in upper 50-m depth [START_REF] Zhao | Comparison in the distribution of microbial food web components in the Y3 and M2 seamounts in the tropical western Pacifi c[END_REF][START_REF] Zhao | Virioplankton distribution in the tropical western Pacifi c Ocean in the vicinity of a seamount[END_REF].

Although ciliates trait structure was proposed by previous studies [START_REF] Kršinić | On vertical distribution of tintinnines (Ciliata, Oligotrichida, Tintinnina) in the open waters of the South Adriatic[END_REF][START_REF] Dolan | Declines in both redundant and trace species characterize the latitudinal diversity gradient in tintinnid ciliates[END_REF][START_REF] Yu | Diff erences in the vertical distribution and response to freshwater discharge between aloricate ciliates and tintinnids in the East China Sea[END_REF][START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF], there were few studies on their variations. The variations of size classes and tintinnid proportion were compared through a large space scale in the Arctic, Subarctic, and tropical waters [START_REF] Wang | Diff erence of planktonic ciliate communities of the tropical West Pacifi c, the Bering Sea and the Arctic Ocean[END_REF]. The present study contributed to the study on ciliate trait structure variations in the tropical waters with diff erent hydrological conditions.

Large ciliates in the surface layer over waters with shallow DCM depth

The DCM depth over Yap seamount was shallower than other two seamounts, and the integrated abundance over Yap seamount was lower than Mariana seamount, while the integrated biomass over Yap seamount was similar to Mariana seamount (Table 2). The large size-fraction (>30 μm) abundance proportion was also higher than Mariana and Caroline seamounts. For tintinnids, the 12-16-μm LOD sizeclass had low abundance proportion in the shallow DCM depth water over Yap seamount. These results indicate that the largest aloricate ciliates and tintinnids were dominant component over waters with a shallow DCM depth. This phenomenon also occurred over Cobb seamount where the ciliate trait structure showed a low abundance but a large size [START_REF] Sime-Ngando | Ciliated protozoan communities over Cobb Seamount: increase in biomass and spatial patchiness[END_REF].

Relationships between tintinnids and upwelling

A weak upwelling was put in evidence over Mariana seamount. Previous studies showed that nutrients (e.g., NO 3 -N, dissolved inorganic nitrogen (DIN), PO 4 -P, and SiO 3 -Si)) had obvious uplifts around this seamount [START_REF] Ma | Environmental characteristics in three seamount areas of the tropical Western Pacifi c Ocean: focusing on nutrients[END_REF]. However, no seamount eff ect was evidenced on phytoplankton [START_REF] Zhao | Comparison in the distribution of microbial food web components in the Y3 and M2 seamounts in the tropical western Pacifi c[END_REF][START_REF] Dai | The seamount eff ect on phytoplankton in the tropical western Pacifi c[END_REF]. The reason might be the shallowness of the weak upwelling. The nutrient concentrations were high in waters deeper than 100 m over three seamounts [START_REF] Ma | Environmental characteristics in three seamount areas of the tropical Western Pacifi c Ocean: focusing on nutrients[END_REF]. Therefore, the weak upwelling in the upper 50 m could not cause a strong increase of nutrient concentration in surface water, which could not eventually enhance phytoplankton growth signifi cantly. Aloricate ciliate group also did not have any obvious increase around the seamount compared to off seamount. For tintinnids, group Ⅳ (deep-water dweller) had a higher abundance proportion than group I (surface water dweller) at 30-m depth. This observation was diff erent for Yap and Caroline seamounts, and along a transect in the tropical Western Pacifi c Ocean (Wang et al., 2019a). This might refl ect a seamount eff ect over Mariana seamount.

Tintinnids and DCM depth

Our results showed that DCM depth infl uenced tintinnid trait structure. When the DCM depths were shallower than 100 m, the tintinnid highest abundances appeared at these DCM depths. While when DCM depths were deeper than 100 m, they were in about 100-m depth (Fig. 4). This result was consistent with a transect in the tropical Western Pacifi c Ocean described in Wang et al. (2019a). The around 100-m depth might be a deepest boundary for tintinnid maximum abundance. In the tropical seas, DCM is a common feature, with depths varied from 50 to 150 m [START_REF] Martin | Latitudinal variation in plankton size spectra in the Atlantic Ocean[END_REF]. Whether these relationships reoccur in diff erent parts of the oceans need more investigations.

Salpingella faurei is a cosmopolitan deep water species and abundant when DCM was deeper than 100-m depth over Mariana and Caroline seamounts. This species distributed from tropical seas to Arctic Oceans [START_REF] Dolan | The Biology and Ecology of Tintinnid Ciliates: Models for Marine Plankton[END_REF][START_REF] Li | Tintinnid diversity in the tropical west pacifi c ocean[END_REF]Wang et al., 2019a, b). In the Bering Sea and the Arctic Ocean, S . faurei mainly distributed in the DCM layers. While at two stations with temperature and salinity well mixed in the Bering Sea, this species could also be brought to surface waters (Wang et al., 2019b).

CONCLUSION

In this study, we compared the planktonic ciliate trait structure over Yap, Caroline, and Mariana seamounts. We evidenced that weak upwelling and DCM depths infl uence the planktonic ciliate trait structure. At the surface layer with shallower DCM, the large size-fraction (>30 μm) abundance proportion to aloricate ciliates increased. The weak upwelling evidenced on Mariana seamount could take group Ⅳ tintinnid species into surface waters, which lead to a higher abundance of this group compared to group I at 30-m depth. The depths of maximum tintinnid abundance were around 100 m even when DCM deeper than 100 m. Despite of the depth variation of DCM, tintinnid abundance proportion to total ciliate was highest at about 100-m depth. Our results represent a snapshot infl uence of the upwelling and the deeper DCM to ciliate trait structure. It constitutes a baseline for further comparative studies on the temporal and spatial environmental infl uence on ciliate in other seas.

DATA AVAILABILITY STATEMENT

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

ACKNOWLEDGMENT

Special thanks to the great eff orts of the captain and crew of R/V Kexue during the cruises over Yap, Mariana, and Caroline seamounts.

Fig. 1

 1 Fig.1 Survey stations and transects over Yap (a), Mariana (b), and Caroline (c) seamounts in the tropical Western Pacifi c Ocean Red dashed box means transects A, C and E.

Fig. 2

 2 Fig.2 Vertical distribution of temperature ( T ), salinity ( S ), and chlorophyll a (Chl a ) in vivo fl uorescence from surface down to 200-m depth Black dots are the sampling points. St. means station.

Fig. 3

 3 Fig.3 Vertical distribution of ciliate abundance and biomass from the surface down to 200-m depthBlack dots are the sampling points. St. means station.

Fig. 4

 4 Fig.4 Vertical distribution of the average planktonic ciliate (aloricate ciliates and tintinnids) abundance in every sampling layers (a), with a zoom for tintinnids (b), and contribution of the two groups to the total ciliate abundance (c) DCM: deep chlorophyll a maximum layer (black dashed box).

Fig. 5

 5 Fig.5 Average abundance proportion (AP) of the diff erent sized aloricate ciliates in every sampling layersAloricate ciliate size fraction were divided into 10 μm increments (a), which further classifi ed into small (10-20 μm), medium (20-30 μm) and large (>30 μm) (b).

Fig. 6

 6 Fig.6 Photomicrographs of several tintinnid species found in this study a. Steenstrupiella steenstrupi ; b. Steenstrupiella intumescens ; c. Cyttarocylis eucecryphalus ; d. Salpingella faurei ; e. Epiplocylis constricta ; f. Canthariella brevis ; g. Steenstrupiella gracilis , h. Ascampbelliella armilla ; i. Eutintinnus hasleae .

Fig. 7

 7 Fig.7 Vertical distribution of diff erent types of tintinnid abundances and its abundance proportion (AP) in every sampling layersGroup I: only appeared in 0-100 m; Group Ⅱ: only appeared in 50-200 m; Group Ⅲ: only appeared in >100 m; Group Ⅳ: appeared in every depth.

Fig. 8

 8 Fig.8 Number of species ( N ) and its abundance proportion (AP) for each lorica oral diameter (LOD) size-class of tintinnid species

Table 1 Ciliate integrated abundance and biomass from surface down to 200-m depth over Yap, Mariana, and Caroline seamounts

 1

	Seamount	AR (inds./L)	IAR (×10 3 inds./m 2 ) AIA (×10 3 inds./m 2 )	BR (μg C/L)	IBR (mg C/m 2 )	AIB (mg C/m 2 )
	Yap	24-331	71.7-196.3	147.1±32.7	0.03-0.64	0.18-0.33	0.26±0.05
	Mariana	33-405	109.6-254.6	163.3±36.4	0.02-0.87	0.17-0.41	0.26±0.06
	Caroline	41-405	150.5-238.7	195.3±25.7	0.02-1.06	0.23-0.41	0.33±0.05
	AR: abundance range; IAR: integrated abundance range; AIA: average integrated abundance; BR: biomass range; IBR: integrated biomass range; AIB:
	average integrated biomass.

Table 2 Summary of all tintinnid species over Yap, Mariana, and Caroline seamounts

 2

	Species	Yap	Mariana Caroline	BDP	Species	Yap	Mariana Caroline	BDP
	Acanthostomella conicoides	Ⅱ

Table 3 Spearman's rank correlation between the planktonic ciliate abundance and depth, temperature, salinity, and chlorophyll a (Chl a ) in vivo fl uorescence

 3

	Seamount Abundance (inds./L)	Depth (m)	Temperature (°C)	Salinity	Chl a (μg/L)
		Aloricate ciliate	-0.732 **	0.739 **	-0.059 0.205
	Yap	Tintinnid	-0.560 **	0.535 **	0.206 0.348 **
		Total ciliate	-0.727 **	0.735 **	-0.042 0.220 *
		Aloricate ciliate	-0.501 **	0.514 **	0.543 ** 0.008
	Mariana	Tintinnid	-0.164 *	0.177 *	0.545 ** 0.351 **
		Total ciliate	-0.492 **	0.505 **	0.560 ** 0.033
		Aloricate ciliate	-0.448 **	0.411 **	-0.352 ** 0.235 **
	Caroline	Tintinnid	-0.126	0.183 *	-0.261 ** 0.180 *
		Total ciliate	-0.420 **	0.390 **	-0.351 ** 0.235 **
	** : P <0.01,

* : P <0.05, t -test.

* Supported by the Science & Technology Basic Resources Investigation Program of China (No. 2017FY100803), the Strategic Priority Research Program of the Chinese Academy of Sciences (No. XDB42000000), the China Postdoctoral Science Foundation (No. 2020M672149), the National Key Research and Development Program of China (No. 2017YFA0603204), the CNRS-NSFC Joint Research Projects Program (No. NSFC 41711530149), the 2017-2019 Sino-French Cai Yuanpei Program, and the International Cooperation Project-Dynamics and Function of Marine Microorganisms (No. ICP-DYF2M): Insight from Physics and Remote Sensing

Electronic supplementary material

Supplementary material (Supplementary Fig. S1) is available in the online version of this article at https://doi. org/10.1007/s00343-021-0476-4.