Yannick Hallez

Jacques Magnaudet

Buoyancy-induced turbulence in a tilted pipe

Keywords: stratified turbulence, turbulent convection, turbulent mixing

Numerical simulation is used to document the statistical structure and better understand energy transfers in a low-Reynolds-number turbulent flow generated by negative axial buoyancy in a long circular tilted pipe under the Boussinesq approximation. The flow is found to exhibit specific features which strikingly contrast with the familiar characteristics of pressure-driven pipe and channel flows. The mean flow, dominated by an axial component exhibiting a uniform shear in the core, also comprises a weak secondary component made of four counter-rotating cells filling the entire cross-section. Within the cross-section, variations of the axial and transverse velocity fluctuations are markedly different, the former reaching its maximum at the edge of the core while the latter two decrease monotonically from the axis to the wall. The negative axial buoyancy component generates long plumes travelling along the pipe, yielding unusually large longitudinal integral length scales. The axial and crosswise mean density variations are shown to be respectively responsible for a quadratic variation of the crosswise shear stress and density flux which both decrease from a maximum on the pipe axis to near-zero values throughout the near-wall region. Although the crosswise buoyancy component is stabilizing everywhere, the crosswise density flux is negative in some peripheral regions, which corresponds to apparent counter-gradient diffusion. Budgets of velocity and density fluctuations variances and of crosswise shear stress and density flux are analysed to explain the above features. A novel two-time algebraic model of the turbulent fluxes is introduced to determine all components of the diffusivity tensor, revealing that they are significantly influenced by axial and crosswise buoyancy effects. The eddy viscosity and eddy diffusivity concepts and the Reynolds analogy are found to work reasonably well within the central part of the section whereas non-local effects cannot be ignored elsewhere.

Introduction

Buoyancy-driven turbulence is encountered in a broad range of fields, such as inertial confinement, supernovae explosions, volcanic or lake eruptions [START_REF] Woods | Turbulent plumes in Nature[END_REF], mixing in unstable atmospheric layers and upper ocean [START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF] or in engineering applications as diverse as combustion, ventilation [START_REF] Linden | The fluid mechanics of natural ventilation[END_REF], fire propagation or oil extraction. Most of the investigations on this generic form of turbulence have focused on Rayleigh-Taylor (RT) turbulence in the unconfined situation where turbulent statistics are inherently time-dependent (except under special circumstances, e.g. [START_REF] Chung | Direct numerical simulation and large-eddy simulation of stationary buoyancy-driven turbulence[END_REF]) and inhomogeneous in the direction of acceleration. Thanks to modern experimental [START_REF] Dalziel | Self-similarity and internal structure of turbulence induced by Rayleigh-Taylor instability[END_REF][START_REF] Ramaprabhu | Experimental investigation of Rayleigh-Taylor mixing at small Atwood numbers[END_REF][START_REF] Banerjee | Detailed measurements of a statistically steady Rayleigh-Taylor mixing layer from small to high Atwood numbers[END_REF]) and computational [START_REF] Ristorcelli | Rayleigh-Taylor turbulence: self-similar analysis and direct numerical simulations[END_REF][START_REF] Livescu | Buoyancy-driven variable-density turbulence[END_REF], 2008) techniques, several fundamental processes at work in RT turbulence have been studied in detail and are now better understood.

In parallel, a research stream aimed at investigating buoyancy-induced mixing in long pipes and channels has developed throughout the last decade. Most of it focused on the vertical case in which the RT turbulence phenomenology combines with confinement effects. While the width of the mixed layer is known to grow quadratically in time during the nonlinear stages of the instability in the unconfined configuration, the confined case exhibits very different evolutions. Since the lateral size of the large-scale structures is constrained, this growth is approximately governed by a diffusive process, so that the mixed layer evolves as the square root of time when the turbulent Reynolds number is low enough [START_REF] Debacq | Self-similar concentration profiles in buoyant mixing of miscible fluids in a vertical tube[END_REF][START_REF] Debacq | Buoyant mixing of miscible fluids of varying viscosities in vertical tubes[END_REF]. At higher Reynolds number, mixing length arguments and experiments in tall channels with a square cross-section rather suggest a growth with a two-fifth power of time [START_REF] Dalziel | Mixing efficiency in high-aspect-ratio Rayleigh-Taylor experiments[END_REF], an evolution confirmed by a hierarchy of one-dimensional models and large-scale computations (Lawrie & Dalziel 2011a,b). Several other investigations inspired by issues related to Rayleigh-Bénard turbulence aimed at detailing the flow structure, using either miscible fluids of different densities in a long cylindrical tube [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF] or differential heating between the top and bottom plates of a long cell with a square cross-section [START_REF] Gibert | Heat convection in a vertical channel: plumes versus turbulent diffusion[END_REF][START_REF] Tisserand | Convection in a vertical channel[END_REF]. These studies, which made extensive use of particle image velocimetry (PIV), provided evidence for the existence of long plumes or 'channels' of heavy (cold) and light (hot) fluid. In particular, thanks to conditional averaging, the latter two revealed that the flow exhibits an S-shaped 'mean' velocity profile during long transients separated by randomly distributed reversals. The plume-like structures were found to be of central importance in the mixing process and flow dynamics, which contrasts with the conclusions of [START_REF] Dalziel | Mixing efficiency in high-aspect-ratio Rayleigh-Taylor experiments[END_REF] who, by means of PIV, observed the large-scale structures involved in the growth of the mixed layer looking rather as 'billows' with a characteristic size in every direction of the order of the channel width. The experiments by [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF] and [START_REF] Gibert | Heat convection in a vertical channel: plumes versus turbulent diffusion[END_REF], as well as the recent computations by [START_REF] Schmidt | Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell[END_REF] in a tall cylindrical Rayleigh-Bénard cell with an imposed axial temperature gradient, all exhibit the inertial (or 'ultimate') scaling where the Nusselt and turbulent Péclet numbers vary as the square root of the Rayleigh number, i.e. the flow statistics do not depend on the molecular diffusivity.

Another series of studies focused on the same problem in long tilted pipes and channels. For a given set of fluid properties and low-turbulent-Reynolds-number conditions, the mean density evolution was again found to follow a diffusive behaviour with an equivalent diffusivity increasing strongly with the tilt angle, Θ, up to a maximum. Beyond that critical angle, the transverse mixing decreases dramatically, since the light and heavy fluids tend to segregate under the effect of the transverse (stabilizing) buoyancy force [START_REF] Séon | Buoyant mixing of miscible fluids in tilted tubes[END_REF]). Similarly, the front velocity of the two currents increases with Θ when the pipe is close to vertical whereas it decreases with the tilt angle at large Θ, exhibiting a clear 'plateau' corresponding to a balance between inertial and gravitational effects at intermediate angles [START_REF] Séon | Buoyancy driven miscible front dynamics in tilted tubes[END_REF][START_REF] Hallez | Effects of channel geometry on buoyancy-driven mixing[END_REF]. At low-to-moderate Reynolds numbers, the effective diffusivity and front velocity were shown to decrease as the inverse of the three-half and three-quarter powers of the Reynolds number, respectively [START_REF] Séon | Front dynamics and macroscopic diffusion in buoyant mixing in a tilted tube[END_REF]. The local characteristics of the turbulent field, especially the distribution of the variance of velocity fluctuations and that of the turbulent shear stress, started to be explored, both experimentally (using PIV) and computationally (with the present code), by [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF]. This study revealed that the near-wall region plays little role in the streamwise mean momentum transfer and showed existence of a four-roll mean secondary flow, the generation mechanism and scaling of which were elucidated by Hallez & Magnaudet (2009). Subsequent measurements by [START_REF] Znaien | Flow structure and momentum transport for buoyancy driven mixing flows in long tubes at different tilt angles[END_REF] revealed that, with large enough density contrasts and Reynolds numbers and small enough Θ, the mean velocity and all second-order moments including the variance of density fluctuations are stationary, although the mean density profile keeps on evolving. Last, Rayleigh-Bénard convection in a long inclined cell was recently considered by [START_REF] Riediger | Heat transport regimes in an inclined channel[END_REF] who, by recording the variations of the mean streamwise temperature gradient with the injected power, were able to identify several distinct flow regimes. In particular, they found 'hard' turbulence (for which the above 'ultimate' scaling holds) to occur for Θ 20 • at large enough Rayleigh numbers, followed by 'soft' turbulence (for which the Nusselt number varies as the square of the Rayleigh number) at larger Θ. For Θ 20 • , the corresponding mean velocity and shear stress profiles are similar to those of [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF], and the smaller Θ the lower the critical Rayleigh number required to reach the 'hard' regime. New data by the same group in the range 5 • Θ 20 • [START_REF] Salort | Turbulent velocity profiles in a tilted heat pipe[END_REF] were interpreted in the light of a mixing length model, but the observed Θ-dependence of the mean velocity profile could only be recovered by altering the usual time scale of Prandtl's original model (i.e. the inverse of the mean shear) with contributions resulting from both the stabilizing and destabilizing buoyancy components. [START_REF] Odier | Fluid mixing in stratified gravity currents: the Prandtl mixing length[END_REF] and Odier, Chen & Ecke (2012) considered a different but related configuration in which an upflowing pressure-driven turbulent light current is confined in between a heavier fluid initially at rest and an upper inclined plate. Using jointly PIV and planar laser-induced fluorescence, they determined the crosswise shear stress and density flux induced by the entrainment process that takes place at the interface. Results indicate that, in this situation where the two buoyancy components are stabilizing, both quantities may be related to the mean shear and mean density gradient via an almost constant mixing length, the value of which was found to scale with the Corrsin length scale to be introduced later.

According to the above review, most available studies in confined geometries focused on the global mixing efficiency and merely characterized it through the temporal evolution of the cross-sectional average density profile. Some of the most recent ones employed PIV to investigate the instantaneous flow field and the large-scale turbulent statistics in the vertical midplane of the channel or pipe, providing mean velocity and shear stress profiles. However none of them explored the complete three-dimensional structure of the velocity and density fields which is of importance in the mixing process, as emphasized by the discovery of the mean secondary flow in the circular pipe configuration [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF]. Also, apart from those by [START_REF] Odier | Fluid mixing in stratified gravity currents: the Prandtl mixing length[END_REF][START_REF] Odier | Understanding and modeling turbulent fluxes and entrainment in a gravity current[END_REF], no joint measurement of the velocity and scalar fields has been reported so far. Finally, only mean streamwise momentum balances and overall scalar balances have been examined. Hence, compared with canonical pressure-driven flows in channels and pipes, for which extensive statistics are available, the present knowledge of confined buoyancy-driven turbulent flows appears to be still in its infancy. The same is true with respect to homogeneous stably stratified sheared turbulence which has been the subject of many detailed experimental and computational investigations (e.g. Chung & Matheou 2012 and references therein), but also with respect to unconfined RT turbulence for which a substantial amount of statistics and energy budgets is now available thanks to recent experiments [START_REF] Banerjee | Detailed measurements of a statistically steady Rayleigh-Taylor mixing layer from small to high Atwood numbers[END_REF]) and extensive computations [START_REF] Cabot | Statistical measurements of scaling and anisotropy of turbulent flows induced by Rayleigh-Taylor instability[END_REF].

Consequently, the purpose of the present study is to better document the specificities of this class of flows by examining various statistical moments and approximate budgets for first-and second-order statistical quantities resulting from a numerical simulation of the low-Reynolds-number buoyancy-driven turbulent flow in a long pipe with Θ = 15 • . This computation may be considered as the counterpart of the experiments reported by [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF][START_REF] Znaien | Flow structure and momentum transport for buoyancy driven mixing flows in long tubes at different tilt angles[END_REF] at this particular angle, since all dynamical and geometrical parameters, including the pipe aspect ratio, are identical; detailed comparisons between experimental first-and second-order large-scale velocity statistics and some of the present numerical predictions may be found in the former paper. The only noticeable difference between the experimental and numerical set-ups lies in the mechanism that drives the small-scale scalar mixing, since we replace the small but non-zero molecular diffusion of the real system by some effect of the numerical scheme to be described below.

The structure of the paper is as follows. Section 2 describes the physical problem, the averaging procedure and the computational approach. A grid convergence study is provided in appendix A to show that the selected grid allows all significant scales of the velocity field to be correctly captured. Section 3 presents the large-scale firstand second-order statistics of the velocity and scalar (i.e. density) fields. The mean momentum and scalar balances and their implications on the crosswise shear stress and scalar flux are discussed in § 4; a rigorous derivation of these balances is provided in appendix B. Section 5 focuses on the characteristic length scales of the flow and the anisotropy they reveal. Various aspects of the energetics of the fluctuating velocity and scalar fields are discussed in § 6 with a particular focus on the production, turbulent transport and confinement mechanisms that control the corresponding variances, and on the influence of the crosswise buoyancy component on the spatial distribution of the crosswise shear stress and scalar flux. An original algebraic model incorporating destabilizing and stabilizing buoyancy effects is introduced in appendix C to relate the various components of the diffusivity tensor to those of the Reynolds stress tensor. A turbulent diffusivity and an eddy viscosity are computed for the crosswise scalar flux and shear stress, respectively, and the validity of the Reynolds analogy is assessed. Section 7 summarizes the main findings and compares some of the most typical statistics of the present flow with those of several well-documented flows of particular relevance.

Statement of the problem and computational strategy

The physical problem

We consider an incompressible buoyancy-driven turbulent flow taking place in a long circular pipe of length L, diameter d, inclined at an angle Θ from vertical (figure 1). The Boussinesq approximation is assumed to hold. The pipe is filled with two perfectly miscible, Newtonian viscous fluids having densities ρ 1 and ρ 2 (ρ 2 > ρ 1 ) and the same kinematic viscosity, ν. Their small but non-zero molecular diffusivity is not explicitly taken into account; the way the small-scale scalar mixing is achieved is discussed in the next subsection. In the initial state, the light (respectively heavy) fluid occupies the lower (respectively higher) half of the pipe and the sharp interface that separates the two fluids is perpendicular to its axis. Subsequently, the two fluids mix and the local density ρ of the mixture is related to the local volume fraction (which will generally be referred to as the 'concentration') C of the heavy fluid through

ρ = Cρ 2 + (1 -C)ρ 1 = ρ 0 [1 + 2At(C -1/2)],
(2.1) with ρ 0 = (ρ 2 + ρ 1 )/2 and At = (ρ 2ρ 1 )/(ρ 2 + ρ 1 ). The above situation is characterized by four control parameters, namely the Atwood number, At, the tilt angle, Θ, the flow Reynolds number, Re t = (Atgd 3 ) 1/2 /ν based on the gravitational velocity V t = (Atgd) 1/2 , and the pipe aspect ratio, d/L = 1, g denoting gravity.

Local averages, notation and normalization

We make use of a right-handed Cartesian coordinate system (x, y, z) with x along the pipe axis and z in the crosswise direction defined such that the (x, z) plane is vertical; consequently the y-axis is horizontal. Both x and z are positive upwards (making the two components of gravity negative), with x = 0 (which corresponds to the initial position of the interface) located midway between the end walls and y = z = 0 on the axis.

To define the space-time average of all flow variables, we take into account the fact that the flow field varies over O(d) distances in the y and z directions, whereas it only varies over O(L) distances in the x direction. Therefore, around any axial location x, it is possible to define a spatial window of length x, such that d

x L, over which the velocity and concentration are statistically independent of x. A similar reasoning may be adopted with respect to time, by introducing the characteristic velocity U f at which the fronts of the heavy and light currents separate from each other, and defining the averaging time interval t such that d/U f t L/U f . Hence we define the space-time average of any quantity Ψ as

Ψ (X, y, z, T) = 1 x t t+ t/2 t-t/2 x+ x/2 x-x/2 Ψ (x, y, z, t) dt dx, (2.2)
where the stretched variables X = x and T = t indicate that Ψ still varies over streamwise distances of O(d/ ) and times of O(d/ U f ). We also define the crosssectional average Ψ as

Ψ (X, T) = 1 A S S Ψ (X, y, z, T) dy dz, (2.3)
where S denotes the pipe cross-section with area A S . Fluctuations with respect to the space-time average are denoted with primed lowercase letters, i.e. we have Ψ (x, y, z, t) = Ψ (X, y, z, T) + ψ (x, y, z, t). Note that the above definition of the average does not assume any left/right or top/bottom symmetry or antisymmetry of the statistical quantities, although these properties should emerge over long enough averaging times. Compared with pressure-driven flows in planar channels and pipes or convection between parallel plates, where averaging is also performed in the spanwise or azimuthal direction, the price to pay is that present statistics are frequently much less converged. However, as will become clear later, their convergence is sufficient to explore the physical processes that drive the flow at both large and small scales.

Departures from the expected symmetries may also be used to appreciate 'meandering' phenomena that subsist over long periods of time.

The velocity components along the x, y and z axes are U, V and W and will be referred to as the axial, spanwise and crosswise coordinates and velocities, respectively. Indexes i = 1, 2, 3 will often be used to denote components of vector quantities along the x, y, z directions, respectively. In particular the notation u i will refer to the velocity fluctuation u , v or w for i = 1, 2 or 3, respectively.

In the run whose results are discussed below, we select a tilt angle Θ = 15 • , an Atwood number At = 10 -2 and a flow Reynolds number Re t = 886. The pipe is 176d long, replicating that used in the experiments of [START_REF] Znaien | Experimental and numerical investigations of flow structure and momentum transport in a turbulent buoyancy-driven flow inside a tilted tube[END_REF]. The development of the flow as time proceeds may be observed on the two movies available at http://dx.doi.org/10.1017/jfm.2014.638. The averaging procedure is applied over five 3d-long windows located at x = 0, x = ±4d and x = ±12.5d, except for two-point correlations which are obtained over 8d-long windows to account for long-range correlations. Time averaging is performed over the last 180 000 time steps of the computation which is stopped well before the fronts reach the end walls to avoid any significant finite-length effect. From now on all quantities are presented in dimensionless form. For this purpose, lengths and velocities are normalized by the pipe diameter, d, and the gravitational velocity, V t , respectively. Consequently, times and pressures are normalized by d/V t and ρ 0 V 2 t , respectively.

The results to be discussed below were obtained with the in-house JADIM code by solving the three-dimensional, time-dependent incompressible Navier-Stokes equations in the Boussinesq approximation, together with the equation for the volume fraction of the heavy fluid. The numerical schemes and time-advancement procedures are similar to those reported by [START_REF] Hallez | Effects of channel geometry on buoyancy-driven mixing[END_REF] to which the reader is referred for technical details as well as for examples of computed confined buoyancy-driven flows. The code makes use of a finite volume approach in which the Navier-Stokes equations are spatially discretized with second-order centered schemes and are advanced in time through a projection technique. During the predictor step, a third-order Runge-Kutta algorithm for advective and source terms is combined with a semi-implicit Crank-Nicolson algorithm for viscous terms. Incompressibility is then enforced by solving a Poisson equation for the pressure correction. Provided that the grid is fine enough to capture scales of the order of the Kolmogorov microscale, η, this numerical approach is known to be suitable for performing direct numerical simulation (DNS) of incompressible turbulent flows; a recent example of DNS achieved with the present code is provided in [START_REF] Chouippe | Numerical simulation of bubble dispersion in turbulent Taylor-Couette flow[END_REF].

In the transport equation for C, molecular diffusion is not explicitly taken into account, so that no physical high-wavenumber cutoff exists in the scalar spectrum. Since the computational cells have a finite thickness, small-scale fluctuations of C are not resolved and a subgrid-scale model is in principle required to provide correct estimates of the scalar fluxes exchanged between the resolved and unresolved scales. Instead, we solve the transport equation for C using Zalesak's flux corrected transport (FCT) algorithm [START_REF] Zalesak | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]), a numerical scheme initially devised to capture flow discontinuities, especially shocks. FCT, and more generally flux-limiting schemes, have been shown to provide implicit subgrid-scale models, an idea that led to the so-called implicit large eddy simulation (ILES) approach [START_REF] Boris | New insights into large eddy simulation[END_REF][START_REF] Fureby | Monotonically integrated large eddy simulation of free shear flows[END_REF][START_REF] Hanjalic | One-point closure models for buoyancy-driven turbulent flows[END_REF]Grinstein, Margolin & Rider 2007). This is the strategy that was employed, both for the velocity and density fields, by Lawrie & Dalziel (2011a,b) to study RT mixing in a tall vertical channel. Convincing justifications for the ILES approach applied to scalar transport were recently provided by [START_REF] Verma | On filtering in the viscous-convective subrange for turbulent mixing of high Schmidt number passive scalars[END_REF] who performed fully resolved simulations of the passive transport of a low-diffusivity scalar in homogeneous isotropic turbulence (HIT) and noticed that the spectral flux at scales slightly smaller than η, i.e. near the common limit of the inertial-convective and viscous-convective subranges [START_REF] Batchelor | Small-scale variation of convected quantities like temperature in turbulent fluid. Part 1. General discussion and the case of small conductivity[END_REF], is very small and may easily be replaced by the effect of some high-order upwinded numerical scheme without affecting the large-scale part of the scalar spectrum. A posteriori tests carried out with data obtained through such 'velocity-resolved, scalar-filtered' simulations confirmed their conclusions. Thus, provided that the grid captures O(η) scales on the one hand and that small-scale unresolved scalar fluctuations have a negligibly small effect on the velocity field on the other hand, the present approach is expected to faithfully predict the velocity and scalar fields down to scales of the order of some η.

The results reported below were obtained on a cylindrical grid with 64 × 64 × 2816 grid points in the radial, azimuthal and axial directions, respectively. The grid is uniform in the latter two directions and is refined near the wall in the radial direction to capture the strong normal velocity gradients expected within the boundary layer. The performances of this grid resolution are examined in appendix A, especially by considering a somewhat simpler but closely related problem in which the pipe is vertical and the flow is driven by a prescribed axial density gradient whose magnitude is of the same order as that found near x = 0 in the tilted configuration. Results obtained by doubling the radial and azimuthal resolutions only reveal marginal changes on most quantities of interest, especially on the value of η. This is in line with the results of § 5 which show that η is approximately 0.015 in the core and increases up to 0.018 at the outer edge of the near-wall region: based on these estimates, the above resolution results in a grid spacing, ∆, smaller than 2η, a value known to be a good approximation of the size of the smallest eddies (Pope 2000, p. 242); it also satisfies the theoretical DNS criterion ∆ πη [START_REF] Grötzbach | Spatial resolution requirements for direct numerical simulation of the Rayleigh-Bénard convection[END_REF]. To limit computational costs, we adopted a somewhat coarser resolution in the axial direction with ∆ ≈ 4η everywhere. As shown in appendix A, doubling or even quadrupling this axial resolution leaves the large-scale variances and the large-and small-scale anisotropies nearly unchanged and modifies η by approximately 2 %. The marginal grid dependency observed on η and the fact that the energy density near the high-wavenumber cut-off in the velocity (respectively concentration) spectra (not shown) is typically four (respectively three) orders of magnitude lower than in the low-wavenumber range indicate that all energetically significant components of the turbulent field are adequately resolved on the present grid. These various checks provide evidence that the statistics and balances discussed in the rest of the paper are based on robust data.

Large-scale characteristics

Mean fields

Figure 2(a,b) respectively display the distribution of the mean concentration, C, and streamwise velocity, U, in the cross-section x = 0. Within the core (say |y| 0.3, |z| 0.35), both quantities exhibit a remarkably linear variation with the crosswise coordinate z. Nevertheless, the velocity must return to zero at the wall to satisfy the no-slip condition, which results in the presence of an extremum in the U-profile within each half of any vertical plane, leading to an S-shape z-variation.

The secondary mean flow discussed by Hallez & Magnaudet (2009) is displayed in figure 2(c) and (d). The isovalues indicate that this secondary flow whose maximum intensity is typically 3-5 % that of U has a four-roll structure which fills entirely the pipe section. As the vectors in Figure 2(b) show, it brings fluid from the axis toward the lateral wall in the z = 0 plane and then from the wall toward the axis within the vertical midplane, thus reducing the crosswise concentration gradient ∂ z C. This transfer of heavy/light fluid towards higher/lower regions of the section is clearly reflected in the curved shape of the near-wall isopycnals in the range |z| 0.25 (see figure 2a). The mechanism by which the non-zero (V, W) velocities are generated and sustained (which directly stems from the inhomogeneity of the secondary Reynolds stresses v 2 -w 2 and v w ), as well as their magnitude (which may be shown to vary as Re -1/2 t ) were analysed in detail in the aforementioned paper to which the reader is referred. In figure 2(c) and (d), it may be noticed that the distributions exhibit significant left-right and top-bottom asymmetries, respectively. At places where these two asymmetries do not compensate exactly, the local divergence of the secondary flow takes non-zero values which are balanced by slow asymmetric X-variations of U to satisfy the incompressibility condition. Obviously, these variations cannot be seen in figure 2(b), owing to the local axial averaging; they could only be discerned by comparing several distributions of U taken at distinct positions X. However, they are responsible for the left-right and top-bottom asymmetries that may be observed in (d) W. In all figures, the region corresponding to the small white central disc is excluded because the singularity of the cylindrical coordinate system at y = z = 0 frequently induces inaccuracies in the interpolations (especially the differentiations) involved in the post-processing.

figure 2(b). Although all of these asymmetries would disappear by averaging results at symmetric positions ±y and ±z or by increasing the sampling time, they suggest that the flow exhibits a slow but significant meandering. This will be confirmed in § 5 by examining two-point correlations.

3.2. Turbulent intensities, Reynolds stresses and scalar fluxes Figure 3 displays the distribution of the root mean square (r.m.s.) values of the concentration fluctuations and those of the three turbulent velocity fluctuations. As figure 3(a) shows, the r.m.s. values of c experience little variation (0.047 c 2 1/2 0.055), except in the immediate vicinity of the wall. They reach a broad maximum around the core, along a ring of radius r = (y 2 + z 2 ) 1/2 ≈ 0.35. The Ellison length scale L E = -c 2 1/2 /∂ z C can be computed from figures 3(a) and 2(a). This scale may be thought of as the typical distance a fluid particle has to travel along the z direction before it either returns towards its equilibrium position or mixes. We find L E ≈ 0.34, which means that it takes nearly one-third of the pipe diameter for the overturning process to be completed. The distribution of the r.m.s. values of u (figure 3b) is very similar to that of c 2 1/2 , except of course close to the wall where it has to return to zero. Again, a broad maximum with u 2 1/2 ≈ 0.31 is observed at the edge of the core,

-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.
u 2 1/2 ; (c) v 2 1/2 ; (d) w 2 1/2 .
the local minimum on the pipe centreline being approximately 0.275. The spanwise and crosswise fluctuations exhibit a markedly different behaviour. They both reach their maximum on the pipe axis and decrease monotonically toward the periphery.

A very similar distribution, with a typical 15 % decrease of the r.m.s. value from the centre to the median radial position r = 0.25 as observed here, was reported in the confined RT configuration (with v 2 = w 2 ) by [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF].

Here, as is customary in homogeneous shear (HS) flows [START_REF] Champagne | Experiments on nearly homogeneous turbulent shear flow[END_REF]Tavoularis & Corrsin 1981a;[START_REF] Tavoularis | Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence[END_REF], the crosswise fluctuation is somewhat smaller than the spanwise one (the maximum of v 2 1/2 is ≈0.19 whereas that of w 2 1/2 is ≈0.15). In addition, some damping of w due to the crosswise stabilizing stratification cannot be ruled out and will indeed be confirmed later. The v 2 1/2 distribution exhibits an almost circular symmetry, except near the wall where it has to fall faster to zero along the y-axis where v is the normal component of the velocity fluctuation than along the z-axis where it is a tangential component; the w 2 1/2 distribution exhibits a similar near-wall anisotropy with the y and z directions interchanged. In the core, the isovalues of w 2 1/2 appear to be slightly squeezed along the y-axis, suggesting some influence of a z-dependent mechanism.

We shall come back to this point in § 6.3. For future purpose it is appropriate to introduce the Brunt-Väisälä frequency N = (-2∂ z C sin Θ) 1/2 and define the buoyancy length scale L B = w 2 1/2 /N which measures the typical distance a fluid particle could rise in the z direction if all of its kinetic energy in that direction were converted into potential energy [START_REF] Smyth | Length scales of turbulence in stably stratified mixing layers[END_REF]. Results from figures 2(a) and 3(d) yield L B ≈ 0.58, so that the ratio L E /L B of the Ellison and buoyancy scales appears to be approximately 0.6. From (b-d) in figure 3, we may also define an 'averaged' r.m.s. turbulent velocity, u, such that

u = q 2 3 1/2 with q 2 = 3 i=1 u 2 i . (3.1)
The turbulent Reynolds number based on the near-axis values of u in the region r 0.1 is Re u = ud/ν ≈ 185.

Figure 4(a-c) show the distribution of the three second-order velocity cross-

correlation coefficients R ij = u i u j /(u 2 i 1/2 u 2 j 1/2
) (i = j). The R 13 coefficient reaches its maximum in the plane z = 0 and then decreases smoothly as |z| increases. It keeps a constant sign essentially everywhere within the section, although small positive values are found very close to the wall around the vertical midplane, typically for |z| 0.43. Since the isovalues are nearly horizontal in the range |y| 0.35, it may be concluded that the lateral wall has a marginal influence on u w in the core. As we saw in figure 2(b), the velocity gradient ∂ z U changes sign around |z| ≈ 0.36 within this central 'slice' of the pipe, whereas no such change is observed in figure 3(a). Hence, u w and ∂ z U have the same sign within the subregion 0.36 |z| 0.43, i.e. the 'eddy viscosity' -u w /∂ z U is negative.

For obvious symmetry reasons, second-order cross-correlations involving the spanwise fluctuation v must vanish in the y = 0 vertical midplane. They also vanish in the z = 0 midplane since U and W are zero there, so that any v fluctuation has the same probability to be associated with positive or negative u or w . Figure 4(a) indicates that R 12 has a four-lobe structure with non-negligible values throughout most of the section. As ∂ y U vanishes for both y = 0 and z = 0, the extrema of u v are located along the four diagonals. Along each of them, the extremum is reached near the edge of the core which corresponds to the region where the influence of the wall on the axial velocity is maximum. Note that the maxima of R 12 in the core are smaller than those of R 13 by at least a factor of four, which suggests that they only play a secondary role in the momentum transfer throughout the core. The correlation coefficient R 23 (figure 4c) also has a four-lobe structure but only exhibits non-negligible values within the peripheral region where the motion of fluid particles is forced to be almost parallel to the wall. This condition yields R 23 = |1| at the wall, as confirmed by figure 4(c). It also dictates the sign of v w : for instance, in the first quadrant, fluid particles moving up faster than the secondary flow have a positive w but their spanwise velocity is negative on average and exceeds V, so that they have a negative v , resulting eventually in the negative sign of v w . The correlation v w is actually directly involved in the generation of the mean secondary flow, as discussed by Hallez & Magnaudet (2009).

Figure 4(d-f ) display the three components of the normalized turbulent fluxes

F i = c u i /(c 2 1/2 u 2 i 1/2
). Their maxima are found to be very dissimilar: F 1 is larger by a factor of 2.5 than F 3 which itself is typically three times larger than F 2 . Figure 4(d) reveals that the axial component experiences little variation throughout the cross-section, except of course in the peripheral region where it goes to zero at the wall. Therefore its behaviour resembles very much that of u 2 , underlining the direct link between density and axial velocity fluctuations. Figure 4(e) indicates that F 2 is negligibly small throughout the core, in line with the negligibly small spanwise

-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.
.0 -0.5 0 0.5 -1.0 y y z z (a) (b) (c) (d) (e) ( f ) FIGURE 4. (Colour online) Cross-correlations coefficients R ij = u i u j /(u 2 i 1/2 u 2 j 1/2 ) (i = j) and normalized fluxes F i = c u i /(c 2 1/2 u 2 i 1/2 ) within the pipe section: (a) R ; (b) R ; (c) R ; (d) F ; (e) F ; (f ) F . gradient ∂ y C.
Although quite 'noisy' because of the low magnitude of the extrema and hence its poorly converged geometry, an eight-lobe structure can be identified at the pipe periphery, with positive and negative values alternating within each quadrant. In the part of each quadrant closest to the vertical midplane, where the largest values are found and the lobes are 'stuck' to the wall, the spanwise gradient ∂ y C is large (typically of the same order as ∂ z C according to figure 2a) and may thus be expected to drive the generation of c v . Close to the midplane z = 0, where the lobes are distinctly detached from the wall, this gradient is much smaller, suggesting that the generation of c v is controlled by other mechanisms such as the spanwise transport of c v by v itself. Figure 4(f ) shows that the distribution of F 3 is quite similar to that of R 13 : it reaches its maximum on the axis an then decreases monotonically towards the wall, keeping a constant sign all along except close to the plane z = 0 where significant negative values are found near the wall. Since ∂ z C is positive everywhere, these negative values correspond to counter-gradient fluxes; their origin will be elucidated in § 6.3.

Disregarding cross-correlations involving v (because their distribution is clearly connected to the wall influence), the observed behaviours may be divided into two families. On the one hand c 2 , u 2 and c u exhibit little variation within the core, with a broad maximum around it. On the other hand, v 2 , w 2 , u w and c w reach their maximum on the axis and decrease monotonically toward the wall. As we shall see later, this radial decrease has several distinct origins, depending on which component is considered.

Mean concentration and momentum balances

The mean concentration (i.e. density) and axial momentum balances are established and discussed in appendix B. This discussion shows that, at leading order in the pipe aspect ratio , all terms involving derivatives with respect to the axial direction are negligible in the latter. The situation is different for the C balance where the forcing term U∂ X C corresponding to the axial transport of the density gradient is found to play a key role.

Streamwise momentum balance

According to appendix B, the mean momentum balance in the axial direction reads at leading order

2 cos Θ C + V∂ y U + W∂ z U + ∂ y u v + ∂ z u w -Re -1 t ∇ 2 ⊥ U ≈ 0, (4.1)
where ∇ ⊥ denotes the two-dimensional gradient in the cross-sectional plane and C(X, y, z, T) = C(X, y, z, T) -C (X, T) its the local transverse stratification. The first contribution, 2 cos Θ C, represents a forcing resulting from the combined action of the axial gravity component and the stabilizing stratification within the cross-sectional plane. The spatial distribution of this contribution is similar to that displayed in figure 2(a) with a uniform shift of 0.5. Figure 5 shows the distribution of the other three groups of terms involved in (4.1) while figure 6(a) shows the profile of all significant terms in the vertical diametrical plane y = 0. As figure 5(c) reveals, the contribution of viscous effects is significant only near the wall, in the top and bottom parts of the section. The secondary flow provides a small, albeit significant, contribution in two distinct regions. First it tends to lower the magnitude of |U| near the wall by bringing low-velocity fluid toward the high-velocity region located around

C; f , W∂ z U; q , ∂ y u v ; s , ∂ z u w ; +, -Re -1 t ∇ 2 ⊥ U). (b) Concentration budget ( p , U∂ X C; u , W∂ z C; q , ∂ y c v ; s , ∂ z c w ).
the plane y = 0. Then it tends to counteract the mean shear ∂ z U by bringing fluid from the high-|U| region located around y = 0, |z| ≈ 0.35 toward the low-momentum region located around the pipe axis. Finally, the momentum flux associated with the Reynolds stresses (figure 5b) provides the contribution that balances most of the forcing term in the core, thus tending to homogenize the flow field. Its distribution is similar to that of the axial velocity, with its two extrema located close to those of U, albeit somewhat closer to the axis (compare figures 5b and 2b). As figure 6(a) indicates, the contribution of the spanwise Reynolds stress u v is far from negligible since it represents approximately 40 % that of u w in the core.

Apart from the role of the secondary flow and spanwise Reynolds stress (which nearly balance each other in the central region, as shown by figure 6a), the present axial momentum balance is qualitatively similar to that encountered in the buoyancy-driven flow in a vertical channel with differentially heated walls (VCDHW). In particular, the linear stratification in the z direction is responsible for the quadratic variation of the main shear stress u w in both cases.

Mean concentration balance

According to the analysis of appendix B, the mean concentration balance reads at leading order in

U∂ X C + V∂ y C + W∂ z C + ∂ y v c + ∂ z w c ≈ 0. (4.2)
Figure 7 shows the distribution of the three groups of terms involved in (4.2); the profile of the various terms along the vertical diametrical plane is shown in figure 6(b). Both figures display the distribution of U∂ X C rather than that of U∂ X C . However, the striking similarity between figure 7(a) and the distribution of U shown in figure 2(b) proves that ∂ X C is essentially constant throughout the cross-section. This confirms that the X-derivative of C is negligibly small, implying ∂ X C ≈ ∂ X C , as concluded in appendix B. In other words, the source that drives the variations of the turbulent fluxes c v and c w across the section (and hence indirectly those of C) is the non-uniform transport of the small axial mean concentration gradient ∂ X C by the axial mean flow U. Although 3-4 times smaller, the contribution of the secondary flow to the transport of C is significant. As shown in figure 7(b), it tends to homogenize the fluid within the section by bringing light (respectively heavy) fluid from the top (respectively bottom) to the centre along the z-axis, and heavy (respectively light) fluid to the top (respectively bottom) along the wall. Finally, the turbulent fluxes bring fluid from the well-mixed near-axis region where C ≈ C toward the high-and low-concentration regions near the top and bottom. This is why the spatial distribution of their two-dimensional divergence (figure 7c) is close to that of the axial velocity, especially in the central 'slice' |y| 0.25. Again, the contribution of the spanwise flux c v is significant, being typically one-quarter that of c w in the core.

Note that the present quadratic variation of c w contrasts with the flat profile of the turbulent heat flux observed in the core of pressure-driven channel or pipe flows or VCDHW flows with isothermal walls. The reason is of course that the source term U∂ X C has no counterpart in those cases, forcing the total crosswise heat flux to be z-independent.

Length scales and dissipation

5.1. Two-point autocorrelations and integral scales Figure 8(a) shows the typical behaviour of the longitudinal two-point velocity correlations, ρ ii (r), in the core within the central vertical strip |y| 0.175. In particular, it reveals that the spanwise movements remain significantly correlated over distances up to r ≈ 5 whereas the axial and crosswise movements are virtually uncorrelated beyond r = 2. Unlike ρ 11 and ρ 33 , ρ 22 exhibits strong oscillations with a marked minimum and a secondary maximum at r ≈ 1.6 and r ≈ 4.2, respectively. These features indicate a slow meandering of the flow in the spanwise direction, in line with the significant left-right asymmetry observed in the V-distribution in figure 2. The above two-point correlations were used to compute the integral length scales, L x ii , through two distinct definitions, namely the integral of the correlation up to its first zero (all correlations exhibit negative values within some range of r), and the specific value of r at which the correlation reaches this first zero. We find L x 22 /L x 33 to be approximately 1.35 with both definitions, while the ratio L x 11 /L x 33 is in the range 2.6-3.9, depending on which definition is considered (the two ratios would respectively be 1.0 and 2.0 in HIT, respectively).

Figure 8(a) also displays the longitudinal two-point scalar correlation, ρ cc (r). It is immediately seen that it keeps significant positive values over larger distances than any velocity fluctuation. This results in large values of the corresponding axial From top to bottom in the core, the solid lines refer to λ x (u ), λ x (v ), λ x (w ), the dashed lines to λ y (v ), λ y (u ), λ y (w ) and the dash-dotted lines to λ z (w ), λ z (u ), λ z (v ), respectively.

integral scale: using the definition based on the position of the first zero, we find L x cc ≈ 6.75 (to be compared with L x 11 ≈ 2.19), an indication that the large-scale density field contains long plumes travelling along the pipe, in line with the findings of [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF], [START_REF] Gibert | Heat convection in a vertical channel: plumes versus turbulent diffusion[END_REF] and [START_REF] Tisserand | Convection in a vertical channel[END_REF] in vertical pipe and channel configurations.

Taylor microscales and mechanical dissipation

Figure 8(b) displays the profiles along the z-axis of all nine Taylor microscales, λ j (u i ), computed using the classical definition

λ 2 j (u i ) = 2u 2 i (∂u i /∂x j ) 2 , ( 5.1)
(no summation on i). Note that computing the y-and z-derivatives of Cartesian velocity components from data obtained on a cylindrical grid requires geometrical transformations that induce some inaccuracy near the singular axis y = z = 0, leading to the oscillations (respectively colour blobs) visible in the central part of on some curves (respectively maps). Two major features emerge from figure 8(b). One is the large anisotropy of the small-scale flow, with in particular streamwise length scales λ x (u i ) significantly larger than their spanwise and crosswise counterparts for all three velocity components (as shown in appendix A, this is by no means an artifact related to the coarser resolution employed in the axial direction). The second is the weak influence of the wall in the central half of the section: except for λ x (u ) which exhibits a 15 % increase from the axis to |z| = 0.25 and to some extent λ x (w ) which decreases by roughly 10 % over the same distance, all other λ j (u i ) vary only by a few per cent within this range. This indicates that the small-scale structure is nearly unaffected by confinement effects within the core. Examining in more detail the isotropy/anisotropy properties revealed by figure 8(b), it may first be observed that

λ y (u ) ≈ λ y (w ) ≈ λ z (u ) ≈ λ z (v )
and that the latter four length scales are nearly equal i = 1 2 3

i / k 0.57 0.24 0.19 u 2 i /q 2 0.57 0.26 0.17 TABLE 1. Comparison between the normalized diagonal components of the mechanical dissipation and those of the normalized variance of the fluctuating velocity (values are taken in the central region r 0.15).

to ≈ (λ y (v ))/ √ 2 as would be expected for isotropic turbulence. Therefore, some kind of small-scale isotropy is preserved within the section, the only exception being λ z (w ) which is only approximately 80 % of λ y (v ) whereas the two should be equal, were the isotropy in the (y, z) plane strictly satisfied. In contrast all three axial length scales are much larger than those in the y or z directions. In particular λ x (v ) is twice as large as λ z (v ) and λ x (u ) is also twice as large as λ z (w ). Therefore, it may qualitatively be concluded that most of the fluctuating vorticity is concentrated within structures whose axial length is on average twice as large as their transverse size.

Results of figures 8(b) and 3 may be combined to compute the mechanical dissipation rate per unit mass which writes

k = 3 i=1 i = Re -1 t 3 i=1 u 2 i 3 j=1 1 λ 2 j (u i ) , (5.2)
where i denotes the diagonal component of the dissipation tensor along direction i.

As shown in figure 9(g), k exhibits a nearly flat profile in the core with a value close to 0.0275. Then it gently decreases down to a minimum approximately 0.014 reached near the outer edge of the near-wall region. Finally, it sharply increases and peaks at the wall where k ≈ 0.068. Although they all reveal qualitatively similar radial variations, a strong anisotropy exists among the three diagonal components i (not shown), 1 being typically three times larger than the other two in the core. This anisotropy may be compared with that existing among the components of the r.m.s. velocity fluctuations. As shown in table 1, the ratio i / k closely follows the ratio u 2 i /q 2 , giving support to Rotta's model i / k = u 2 i /q 2 (Rotta 1951). Hence, the anisotropy of the large-scale field turns out to have a direct footprint on the smallscale motion. This may be interpreted as a low-Reynolds-number effect, but similar levels of small-scale anisotropy have been reported in uniform shear flows with higher Reynolds numbers (Tavoularis & Corrsin 1981b;[START_REF] Pumir | Persistent small scale anisotropy in homogeneous shear flows[END_REF]. Thus present data may also be seen as suggesting that the small-scale phenomenology is deeply influenced by the mean shear. Extra runs at higher Reynolds number would be necessary to disentangle the two possibilities.

Combining (5.2) with the definition (3.1), we may generalize the classical HIT relationship k = 15Re -1 t u 2 /λ 2 . Then, defining an 'averaged' Taylor microscale, λ, as

λ = Re -1/2 t u √
15/ k , we find λ ≈ 0.165 and

Re λ = Re t λu = Re 1/2 t 5 3 k q 2 ≈ 31. (5.3)
With k at hand, we can compute the three characteristic length scales that allow the range of scales affected by viscosity, stabilizing stratification or shear to be estimated in the framework of Kolmogorov theory. First, the Kolmogorov scale, η = Re -3/4 t -1/4 k , is found to be η ≈ 1.5 × 10 -2 ≈ 0.09λ. The Ozmidov scale, L O = ( k /N 3 ) 1/2 , and the Corrsin scale, L C = { k /(∂ z U) 3 } 1/2 , are indicators of the scale beyond which stabilizing stratification and shear effects affect the local flow structure, respectively. Present data indicate L O ≈ 1.30 and L C ≈ 0.12. These estimates suggest that stratification barely affects the local flow structure at any scale, while the mean shear has a much stronger influence which extends down to structures with a size of a few η. However, it must be kept in mind that the above definitions implicitly assume turbulence to be isotropic, which is far from being the case here. In particular, as shown in table 1, 3 is only approximately 20 % of k , so that L O3 = ( 3 /N 3 ) 1/2 should be a more relevant scale to characterize the possible influence of stabilizing stratification on w . With this definition one finds L O3 ≈ 0.54, which now suggests that stratification has some effect on the large-scale structure of w , although it certainly has a nearly negligible influence on u . Section 6 will confirm this conjecture.

6. Energetics of the turbulent field 6.1. The u 2 i and q 2 balances Results discussed in the previous sections indicate that the turbulent field is strongly anisotropic at all scales and reveal a significant radial inhomogeneity of the diagonal Reynolds stresses u 2 i . To get further insight into the mechanisms that generate and sustain these characteristics, it is necessary to examine the main contributions to the three u 2 i budgets. Numerical data indicate that the time rate-of-change term ∂ T u 2 i and all contributions involving derivatives in the axial direction are negligibly small everywhere. Transport by the secondary flow and viscous diffusion are also found to be very small throughout the core and are significant only in the near-wall region, especially regarding the latter which has to balance i at the wall. Hence, the u 2 i balance in the core essentially reduces to

P i + Φ i -i -∇ ⊥ • J i ≈ 0, (6.1)
where P i is the 'production' term responsible for the conversion of the mean flow kinetic energy or the potential energy into u 2 i , Φ i = p ∂ i u i (no summation on i) is the pressure-strain correlation responsible for energy redistribution between the three components of the fluctuating motion, and the contribution of J i = 1/2u 2 i v ⊥ + p v e y δ i2 + p w e z δ i3 is the lateral flux corresponding to the in-plane transport of u 2 i by turbulent fluctuations, p being the pressure fluctuation and v ⊥ = v e y + w e z the velocity fluctuation within the cross-sectional plane (e y and e z denote the corresponding two unit vectors and δ ij is the Kronecker tensor). The trace of (6.1) yields the turbulent kinetic energy (TKE) balance

P k -k -∇ ⊥ • J k ≈ 0, (6.2)
where

P k = i P i , J k = i J i = {p + (u 2 + v 2 + w 2 )/2}v
⊥ and the detailed contributions to P i write

P 1 = -u v ∂ y U P Uy -u w ∂ z U P Uz -2c u cos Θ P F 1 , (6.3) P 2 = -v 2 ∂ y V -v w ∂ z V, (6.4
)

P 3 = -w 2 ∂ z W -v w ∂ y W -2c w sin Θ P F 3 . (6.5)
The two contributions in P 2 and the first two in P 3 result from secondary flow gradients; their sum, which is involved in P k , is denoted as P sec . Figure 9 displays the spatial distribution of these various terms. The generation mechanisms in the core are clearly dominated by effects of the shear ∂ z U (a) and the axial mass flux c u (d) which have similar magnitudes. Production by the spanwise velocity gradient ∂ y U (b) is negligibly small, except along the two bisectors in the near-wall region. Terms associated with the stabilizing stratification (e) and the secondary flow (c) tend to damp turbulence but have magnitudes typically 15-20 times smaller than the two dominant contributions. In particular, P F3 is less than 3 % of P k , in line with the small value of the shear Richardson number Ri = (N/∂ z U) 2 ≈ 0.046. The negative sign of P sec results from the combination of the already discussed anisotropy v 2 -w 2 > 0 and the stretching/compression around the hyperbolic point y = z = 0 where ∂ y V > 0 and ∂ z W < 0. In contrast, P sec is positive in the top and bottom parts of the near-wall region because P sec ≈ -v 2 ∂ y V there owing to the quadratic near-wall damping of w 2 , and the spanwise velocity V experiences a compression near the hyperbolic points y = 0, z = ±1. Similarly, the vanishing of v 2 and the stretching of W near points y = ±1, z = 0 result in the negative sign of P sec in the corresponding near-wall region.

As shown in figure 9(f ), the total production P k reaches its maximum in the core and decreases gradually toward the wall. This is in stark contrast with the usual case of pressure-driven channel or pipe flows where production is concentrated in the nearwall region and is small in the core. The reason for this difference is of course that the shear stress u w is very weak in the near-wall region in the present case, whereas it is maximum there in pressure-driven flows and decreases linearly to zero in the core. Given the vanishing of the main shear ∂ z U at some crosswise position z = z 0 (y) for any y and the quite flat distribution of P F 1 , the maximums of P k at radial positions r ≈ z 0 are encountered along the y axis, where P U z is still significant. This is why the P k distribution exhibits a pronounced flattening along that axis. Combined with the radial decrease of k at the edge of the core revealed by figure 9(g), the ratio P k / k (figure 9h) is found to exhibit a characteristic 'bat-like' distribution with a maximum approximately 1.75 located near |y| = 0.3 along the y axis and a significantly lower value approximately 1.37 near the axis y = z = 0. Overall, P k / k is larger than unity throughout the entire range of radial locations r 0.43, the integral constraint P k = k being satisfied thanks to the contribution of the thin peripheral layer 0.43 r < 0.5 where P k (which at such locations essentially reduces to P F 1 ) quickly falls to nearzero levels while k reaches its highest values. Again, the present distribution of P k / k contrasts with that of pressure-driven flows where this ratio is larger than unity only within the so-called buffer layer of the wall region.

Given the contrasting orders of magnitude revealed by figure 9, the u 2 i -balances in the core reduce at leading order to

P U z + P F 1 + Φ 1 -1 -∇ ⊥ • J 1 ≈ 0, (6.6) Φ α -α -∇ ⊥ • J α ≈ 0, (6.7)
with α = 2, 3. Hence, turbulent energy is essentially produced in the axial direction by the destablilizing stratification and the mean shear ∂ z U, then redistributed onto the two lateral directions by pressure fluctuations. The distribution of the three Φ i components is displayed in figure 10. They all exhibit a nearly circular symmetry except in the near-wall region. As expected, Φ 1 is negative over most of the section, while Φ 2 and Φ 3 are positive, so that energy is redistributed from u 2 onto v 2 and w 2 . Since the magnitude of Φ 1 decreases from the centre to the periphery, less energy is injected into v and w as r increases. This is one of the two reasons why the variance of both components also decreases monotonically with r, as noticed in figure 3. The redistribution from u to v and w is uneven, with Φ 2 ≈ 1.75 Φ 3 . Hence, it 'feeds' v 2 more than w 2 , a well-known feature of homogenous sheared turbulence [START_REF] Champagne | Experiments on nearly homogeneous turbulent shear flow[END_REF][START_REF] Harris | Further experiments in nearly homogeneous turbulent shear flow[END_REF]. This mechanism goes along with the buoyancy-induced damping of w 2 (the P F 3 term in (6.5)), resulting in a significant anisotropy between v and w , with v 2 ≈ 1.6w 2 (see figure 3). Indeed, although |P F 3 | is less than 3 % of the total TKE production, it is approximately 15 % of Φ 3 near the axis, which significantly contributes to reduce w 2 in the core. After having gradually decreased as r increases, Φ 1 changes sign in the near-wall region r 0.4, where u 2 receives energy from the locally normal fluctuation, i.e. v 2 (respectively w 2 ) near z = 0 (respectively y = 0). The other locally tangential velocity fluctuation, i.e. w 2 (respectively v 2 ) near z = 0 (respectively y = 0), also receives some energy from the locally normal fluctuation. This energy transfer from the wall-normal fluctuation to the tangential fluctuations finds its roots in the kinematic blocking which instantaneously damps any normal velocity fluctuation in the near-wall region [START_REF] Hunt | Free stream turbulence near plane boundaries[END_REF].

Owing to the axial buoyancy contribution P F 1 , P 1 is still large at the position where Φ 1 vanishes. Being 'fed' by potential energy and redistributing only a marginal amount of energy onto the other velocity components or even gaining some energy from one of them around that position allows u 2 to reach a maximum at the edge of the core, typically at r ≈ 0.3.

It may be shown that the sum P i + Φ ii is slightly negative in the case of u 2 while it is positive for the other two components, which also makes its trace positive. Thus, J 1 is directed towards the centre of the pipe while J 2 and J 3 (and J k ) are directed towards the wall. A crude Fickian behaviour of the type J i = -K∇ ⊥ u 2 i where K is a positive constant would correctly predict the sign of the three fluxes since u 2 (respectively v 2 and w 2 ) increases (respectively decrease) from the axis towards the periphery. The positive sign of J k throughout the core is a key feature of the large-scale phenomenology of the present flow: owing to its statistical steadiness and invariance along the axial direction, it is the only mechanism by which it can accommodate the positive imbalance between P k and k by transferring the corresponding energy excess from the near-axis region toward the periphery. Since this transfer is entirely secured by v and w over most of the core (from the axis until approximately the inflection point of the u 2 profile), the magnitude of J 2 and J 3 increases with r in order for the net transfer to remain oriented toward the wall everywhere. This is the second reason that makes v and w decrease monotonically with r, and the larger r the stronger this decrease.

The distribution of the Φ i components reveals no footprint of the flattened shape of P k along the y-axis. Their nearly circular symmetry rather suggests that, beyond the central part of the core, say r 0.15, they are to a large extent controlled by the distance to the wall. This is not unlikely since the aforementioned kinematic blocking is known to act over a distance of the order of the integral scale [START_REF] Hunt | Free stream turbulence near plane boundaries[END_REF], and the transverse integral length scales L x 22 and L x 33 are about half a pipe radius. This implies in turn that, besides radial turbulent transport, wall confinement controls the radial distribution of v 2 and w 2 over a good part of the section. This view is reinforced by the comparison of present results with those obtained in DNS of VCDHW flows [START_REF] Boudjemadi | Budgets of turbulent stresses and fluxes in a vertical slot natural convection flow at Rayleigh Ra = 10 5 and 5[END_REF][START_REF] Versteegh | Turbulent budgets of natural convection in an infinite, differentially heated, vertical channel[END_REF], where it was found that the v 2 and w 2 profiles (corresponding in these flows to the spanwise direction of infinite extension and the wall-normal direction, respectively) have similar shapes for r 0.15 but depart from each other beyond that point, with a much flatter decrease of v 2 . These observations make it clear that the distribution of v 2 and w 2 in buoyancy-driven confined flows is directly influenced by the channel geometry. This contrasts with pressure-driven flows in which the distributions of v 2 and w 2 in planar channel and circular pipe configurations are very similar [START_REF] Eggels | Fully developed turbulent pipe flow: a comparison between direct numerical simulation and experiment[END_REF]. The reason for this striking difference is that turbulence production concentrates near the wall in the latter (so that curvature effects in a circular pipe have a weak influence because the characteristic length scales involved in the turbulence generation process are much smaller than the pipe radius), whereas it concentrates in the core in buoyancy-driven flows. Hence, the large-scale structures produced in this region are influenced by the geometry of the entire section, making the v and w distributions very similar in the present cylindrically symmetric case and significantly different in large-aspect-ratio VCDHW flows.

6.2. The c 2 balance and the scalar-to-mechanical time ratio Although molecular diffusivity is not explicitly taken into account in the computations, the spatial discretization and the FCT algorithm employed to solve the concentration equation result in some truncation error, say f , which may be seen as a local source/sink of C. The fluctuations of this error, f , then yield a term 2f c in the c 2 balance. This term, say -χ, provides the scalar dissipation rate required to balance all other contributions, especially the sources of c resulting from the mean concentration gradients ∂ X C and ∂ z C. Again, the balance may be simplified because numerical data confirm that the time rate-of-change term ∂ T c 2 and the axial transport terms {U∂ X c 2 + ∂ X c 2 u } are negligibly small everywhere, whereas transport by the secondary flow is significant only in the near-wall region. Hence, the c 2 balance in the core reduces to

-2 c u ∂ X C P c 2 X -2c w ∂ z C P c 2 z -∇ ⊥ • c 2 v ⊥ -χ ≈ 0. (6.8)
The two contributions P c 2 X and P c 2 z to the production term P c 2 are displayed in figure 11(a,b). Although the axial gradient ∂ X C is small, its transport by the axial velocity fluctuation u results in a contribution of the same order as that due to the crosswise gradient ∂ z C because u is significantly larger than w and is better correlated with c . Owing to the weak variations of the magnitude of c in the core, the turbulent transport term (not shown) only provides a small negative (≈ -0.7 × 10 -4 ) contribution to (6.8). This flux tends to increase the variance of c at the edge of the core, which is consistent with the broad maximum of c 2 found in the range 0.35 r 0.4 (see figure 3). Although we cannot compute the scalar Taylor microscales since we do not resolve the small-scale fluctuations of C, we can obtain χ as the sum of the other three terms in (6.8). The result is displayed in figure 11(c).

As expected, χ is positive everywhere and exhibits values approximately 8.9 × 10 -4 in the central region. Owing to its definition, it falls to zero at the wall. Indeed, if a finite molecular diffusivity D were considered, the Schmidt number Sc = ν/D would also be finite and one would have χ = (Re t Sc) -1 {2(∂c /∂x j ) 2 -∂ 2 c 2 /∂x 2 j }, i.e. χ would be zero at the wall if the latter were an isopycnal surface. In other words, the quantity χ computed here may be identified with the scalar dissipation that would be observed with a real pair of miscible fluids only sufficiently far from the wall.

The scalar-to-mechanical time ratio T χ = c 2 k /(χq 2 ) is the local ratio of the turnover times of the energy-containing velocity and scalar eddies. Using the above values of χ, one finds T χ ≈ 0.51 in the central region. In HS flows with negligible buoyancy effects, T χ is typically in the range 0.25-0.3 [START_REF] Rogers | An algebraic model for the turbulent flux of a passive scalar[END_REF][START_REF] Kaltenbach | Large-eddy simulation of homogeneous turbulence and diffusion in stably stratified shear flow[END_REF], whereas it is in the range 0.2-0.25 in the core of VCDHW flows with turbulent Reynolds numbers comparable to that of the present flow [START_REF] Boudjemadi | Budgets of turbulent stresses and fluxes in a vertical slot natural convection flow at Rayleigh Ra = 10 5 and 5[END_REF][START_REF] Versteegh | Turbulent budgets of natural convection in an infinite, differentially heated, vertical channel[END_REF][START_REF] Dol | A DNS-based thermal second-moment closure for buoyant convection at vertical walls[END_REF]. The reason why T χ is typically twice as large here can be qualitatively understood by noting that the scalar time scale, t c , may be written as t c = l c /u, where l c is the characteristic length of large-scale scalar fluctuations. In the present flow these fluctuations are in good part due to long axial plumes generated by the ∂ X C gradient. In contrast, in the aforementioned cases, the typical size of large-scale scalar blobs is of the order of that of energetic eddies, l u = u 3 / k , because velocity and scalar fluctuations are both controlled by the mean shear (no vertical temperature gradient exists in VCDHW flows between isothermal plates and turbulence is mostly generated by shear in that case). These arguments suggest that the ratio T χ = l c /l u should indeed be significantly larger in the present flow. They are supported by the comparison of the scalar and axial velocity integral length scales which are such that L x cc /L x 11 ≈ 3.0 here, while the results of [START_REF] Versteegh | Turbulent budgets of natural convection in an infinite, differentially heated, vertical channel[END_REF] indicate L x cc /L x 11 ≈ 0.4 (considering in both cases the definition based on the position of the first zero). DNS of Rayleigh-Bénard convection with Prandtl numbers approximately 0.7 and Reynolds numbers comparable to the present one also yield T χ ≈ 0.5 on the midplane of the cell [START_REF] Hanjalic | One-point closure models for buoyancy-driven turbulent flows[END_REF][START_REF] Otic | Analysis and modelling of the temperature variance equation in turbulent natural convection for low-Prandtl-number fluids[END_REF]. Thus, buoyancy-driven turbulence appears to be consistently characterized by scalar-to-mechanical time ratios typically two times larger than those found in classical shear-driven situations.

Shear stress and concentration fluxes

The shear stress u w and the concentration fluxes c w and c u play a key role in the flow. Here we try to better understand how they interplay with the primary velocity and concentration gradients and the variances of turbulent fluctuations. The balance equation for Ψ = u w , c u or c w may be written symbolically in the form P Ψ + Γ Ψ = 0, where P Ψ stands for the generation/destruction contribution involving only secondorder turbulent moments and mean concentration or velocity gradients, and Γ Ψ gathers all other contributions, especially those involving gradients of pressure, velocity or concentration fluctuations and transport terms. Figure 12 displays the three P Ψ maps. Figure 12(b) indicates that P c u keeps a constant negative sign throughout the flow and experiences only modest variations in the central region. This is in line with the negative sign and flat profile of F 1 observed in figure 4 (see also the behaviour of P F1 in figure 9). The other two maps reveal a very different behaviour of P u w and P c w . , (6.9)

P c u -10.2/-8.2 = -u 2 ∂ X C -1.4/-1.7 -u w ∂ z C -2.0/-0.6 -c w ∂ z U -2.5/-1.1 -2c 2 cos Θ -4.3/-4.8
, (6.10)

P c w +2.2/+0.3 = -u w ∂ X C +0.3/+0.1 -w 2 ∂ z C +3.0/+1.5 -2c 2 sin Θ -1.1/-1.3 , (6.11)
where the two numbers below each term refer to its value near the axis and at |z| = 0.30, respectively, all values having been multiplied by 10 3 . All contributions in P c u have the same sign and the mild decrease of |P c u | from the axis to the periphery results from the quadratic decrease of |u w | and c w , partially compensated by the small increase of the axial velocity variance u 2 and the concentration variance c 2 which is found to provide the largest contribution. In (6.9) and (6.11) we note that the main contribution is proportional to w 2 . As far as production is large, which is the case in the core, the relation between u w (respectively c w ) and P u w (respectively P c w ) is mostly local, i.e. turbulent transport plays only a secondary role. Hence, the z 2 -variation of u w and c w implies that w 2 also varies quadratically with z. This is the reason why the w 2 distribution (figure 3d) is somewhat squeezed along the y-axis.

Still in (6.9) and (6.11), a buoyancy-induced contribution resulting from the pipe inclination (the sin Θ term) is involved, the sign of which is opposite to that of all other terms. Combined with the decrease of u w and c w with |z| resulting from (4.1) and (4.2) and that of w 2 resulting from the wall influence through Φ 3 , this buoyancy effect is strong enough to change the sign of P u w and P c w at the edge of the core. In (6.11), the ratio of the last two terms, i.e. the relative contribution of buoyancy and crosswise mean concentration gradient, may be written in the form (L E /L B ) 2 , where L E and L B are the Ellison and buoyancy length scales introduced in § 3.2, respectively. Similarly, the ratio of the buoyancy and shear-induced contributions in (6.9) may be shown to be of O(Ri 1/2 L E /L B ). With Ri ≈ 0.046 and L E /L B ≈ 0.6, this ratio is found to be approximately 0.12 in the near-axis region while the former is typically three times larger. Hence, stabilizing buoyancy has a non-negligible influence on both budgets, with a stronger magnitude on the crosswise mass flux c w . That u w and c w are able to keep their sign unchanged over large regions where that of P u w and P c w has reversed owing to these buoyancy effects indicates that the contribution of lateral transport terms involved in Γ Ψ is significant and mitigates the stabilizing stratification effect. Qualitatively similar features were noticed by [START_REF] Garg | Stably stratified turbulent channel flows. I. Stratification regimes and turbulence suppression mechanism[END_REF] in a stably stratified horizontal channel flow with isothermal walls (see their figure 16). However, the modest variations of w 2 with |z| in the core of pressure-driven flows and the fact that c w keeps nearly constant values owing to the isothermal wall condition (in contrast with the quadratic variation imposed here by (4.2)) strongly reduce the influence of stratification effects and prevent P u w and P c w from changing sign in that case.

Near the vertical diametrical plane, the radial turbulent flux that maintains the positive sign of c w within regions where P c w is negative corresponds to the contribution -∂ z c w 2 . In contrast, near the plane z = 0, this radial flux is directed along the y-axis and therefore corresponds to the term -∂ y c v w . Given the weak degree of correlation between v and w near z = 0 (see figure 4c), the latter flux is much smaller than the former. Hence, unlike the behaviour observed around the z-axis, c w is forced to follow directly the spanwise variations of its production rate around the y-axis. This is why it is negative in the near-wall region where the negative, buoyancy-induced, contribution -2c 2 sin Θ dominates P c w . In the regions where c w is negative, so is the corresponding contribution P c 2 z in (6.8), thus decreasing the local production of the concentration variance. This is why, within the ring 0.35 r 0.4 where c 2 reaches its maximum whatever the azimuthal position, this maximum is slightly lower in these regions, as may be observed in figure 3(a). Owing to the strong correlation between c and u , this feature is also reflected in the distribution of u (figure 3b).

Eddy viscosity and diffusivities

It has been known for long [START_REF] Batchelor | The Theory of Homogeneous Turbulence[END_REF]) that, in homogeneous anisotropic turbulence, the simplest rational relation between the flux of a passive scalar, c u i , and the mean scalar gradient, ∂ j C, involves a second-order non-symmetric diffusivity tensor D ij , i.e. c u i = -D ij ∂ j C (summation on j implied). Approximate models aimed at expressing the various components of D ij as a function of the Reynolds stress tensor components u i u j and mean velocity gradients ∂ j U i have been derived using either a Lagrangian formalism [START_REF] Riley | The relation of turbulent diffusivities to Lagrangian velocity statistics for the simplest shear flow[END_REF][START_REF] Tavoularis | Effects of shear on the turbulent diffusivity transfer[END_REF] or a direct modelling of the unknown terms in the c u i balance [START_REF] Freeman | Tensor diffusivity of a trace constituent in a stratified boundary layer[END_REF][START_REF] Rogers | An algebraic model for the turbulent flux of a passive scalar[END_REF]. Models generally indicate that, in a simple shear flow with U(z) = z∂ z U, the non-zero off-diagonal components of D ij are negative and such that D 13 < D 31 < 0, an inequality confirmed by simulations (e.g. [START_REF] Kaltenbach | Large-eddy simulation of homogeneous turbulence and diffusion in stably stratified shear flow[END_REF]). In the present case, the various components of the diffusivity tensor cannot be extracted directly from the numerical results because each of the fluxes c u and c w depends on the two concentration gradients ∂ X C and ∂ z C. Hence, a model is required to assess the magnitude of the individual components of D ij and thus better understand the main characteristics of the coupling between velocity and density fluctuations, especially the role of stabilizing or destabilizing buoyancy effects. To this end, a simple model based on a two-time approach is developed in appendix C. This model, which is a direct extension of that of [START_REF] Rogers | An algebraic model for the turbulent flux of a passive scalar[END_REF] to the buoyant case, involves one time scale to relate the scalar variance to its dissipation rate and another time scale to characterize the relaxation of scalar fluxes. Its main outcome is the explicit expression of the diffusivity tensor (C 7) which involves the 'buoyancy' tensor B ij defined in (C 3), whose qualitative effects in vertical and horizontal shear flows are made apparent in (C 8)-(C 9) and (C 10)-(C 11), respectively.

A crude estimate of the two characteristic times required to apply this model to present data may be obtained from (6.8), (6.10) and (6.11) by neglecting non-local effects and defining these time scales as the ratio of the corresponding second-order moment (e.g. c 2 ) over its production rate (e.g. P c 2 ). Considering in each case the value of the production rate in the near-axis region, (6.8) yields T c ≈ 4.95 while (6.10) and (6.11) provide τ ≈ 0.95 and τ ≈ 0.93, respectively. Based on these results, we select the approximate values T c = 5.0 and τ = 1.0 to assess qualitatively the model predictions. This yields a ≈ 1.34, b ≈ 2.48, c ≈ 0.05, d ≈ 0.82 and D ≈ 1.22 in (C 6). Making use of the near-axis values of the second-order moments involved in (C 7), we then obtain

D ij ≈   0.113 0 -0.062 0 0.036 0 -0.007 0 0.015   .
(6.12)

These estimates yield c u ≈ -0.010, c w ≈ 0.002 and c 2 1/2 ≈ 0.049 (obtained through the assumed relation c 2 = (1/2)T c P c 2 ), i.e. F 1 ≈ -0.76 and F 3 ≈ 0.29. This compares well with the numerical results F 1 ≈ -0.75, F 3 ≈ 0.29 and c 2 1/2 ≈ 0.047 provided by figures 4(d,f ) and 3(a), respectively. The estimate (6.12) reveals that D ij is strongly anisotropic: D 11 is about eight times larger than D 33 and |D 13 | is more than half of D 11 . Moreover, since ∂ z C is much larger than ∂ X C, the contribution D 13 ∂ z C provides approximately 80 % of the axial flux c u . This makes it clear that the two components D 11 and D 13 are required to correctly predict c u . In contrast, owing to the smallness of D 31 and that of ∂ X C, the corresponding contribution only provides approximately 5 % of c w , so that the crosswise scalar flux is essentially related to the crosswise mean concentration gradient. If c were assumed to be a passive scalar, other things being equal, the model would predict D 11 ≈ 0.090, D 13 ≈ -0.037, D 31 ≈ -0.015 and D 33 ≈ 0.022, which would yield F 1 ≈ -0.46 and F 3 ≈ 0.40, with now c 2 1/2 ≈ 0.051. This is a clear indication that buoyancy effects play a central role in the generation of the scalar fluxes (through the c 2 terms in (6.10) and (6.11)), be they destabilizing (as they enhance c u as shown by (C 8)), or stabilizing (as they reduce c w and to some extent c u as shown by (C 11) and (C 10), respectively). Since D 31 is small compared with the other components of the diffusivity tensor, there is evidence to suggest that c w is almost independent of the axial concentration gradient. On this basis, we define an eddy diffusivity for the crosswise flux as

D T = -c w ∂ z C . (6.13)
Similarly, (6.9) shows that the scalar fluxes only marginally contribute to the generation of the shear stress u w in the near-axis region because they almost counterbalance each other. This is less true at the edge of the core, owing to the decrease of c w and that of the main production term w 2 ∂ z U. Nevertheless, the change of sign of P u w and that of the mean shear ∂ z U take place close to each other.

Based on these remarks, it is reasonable to consider that u w may be directly related to the local value of ∂ z U whatever z, although this relation will obviously fail in the near-wall region where ∂ z U and u w share the same sign. By so doing we define the eddy viscosity ν T as

ν T = -u w ∂ z U . (6.14)
The distribution of D T and ν T is shown in figure 13 together with that of the turbulent Schmidt number Sc T = ν T /D T . Near the vertical diametrical plane, ν T is seen to decrease gradually from values approximately 0.012 on the pipe axis down to values three times smaller for |z| ≈ 0.34 (this decrease is actually quadratic since it follows that of u w ). A mixing length, l T , may be defined through the relation ν T = l 2 T ∂ z U. Since the mean shear is almost constant in the range |y| 0.2, |z| 0.3, the resulting distribution of l T (not shown) is qualitatively similar to that of √ ν T in this region, i.e. it is linear, with values decreasing from 0.095 on the axis to 0.065 at |z| = 0.3. Keeping in mind that the Corrsin scale near the axis is L C ≈ 0.12, it turns out that l T is slightly smaller than L C in the central part of the core. This is in qualitative agreement with the findings of [START_REF] Odier | Fluid mixing in stratified gravity currents: the Prandtl mixing length[END_REF][START_REF] Odier | Understanding and modeling turbulent fluxes and entrainment in a gravity current[END_REF] who, in their 'gravity current' experiment, found l T /L C ≈ 0.7. However, this proportionality does not subsist at the edge of the core (r ≈ 0. countergradient fluxes c w < 0 found in the region |y| 0.35 near the plane z = 0.

In contrast with the above observations, comparison of figures 13(a) and 4(b) reveals two qualitatively different distributions. In particular, within the strip |z| 0.2, values of ν T in the range 0.2 |y| 0.35 are larger than those encountered in the near-axis region, a direct consequence of the high degree of correlation between u and w in that region and of the decrease of the mean shear as |y| increases. Since local values of w 2 are significantly smaller than near the axis, whereas c u is slightly larger, this spanwise variation cannot be explained on the basis of that of P u w /∂ z U (see (6.9)), implying that the distribution of u w is affected by non-local processes. Consequently, the Reynolds analogy is no longer valid at such spanwise locations, as shown by the sharp increase of the turbulent Schmidt number.

7. Final discussion and concluding remarks 7.1. Salient features of the flow dynamics The results discussed throughout this paper helped identify and quantify specific features of buoyancy-induced turbulence in a tilted pipe which result from the combined effects of shear, destabilizing and stabilizing buoyancy forces and wall Reference The generating mechanism of the entire flow field resides in the initial destabilizing stratification. Therefore, turbulence in the present flow shares important features with RT turbulence, with due consideration for the additional influence of wall confinement as exemplified in the experiments of [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF]. At large scale, this manifests itself in a strong anisotropy of the Reynolds stress tensor: turbulence being entirely generated along the axial direction through the action of the axial flux c u , r.m.s. values of u are significantly larger than those of v and w which are only produced through redistribution by pressure fluctuations. This generation process also implies a strong degree of correlation between c and u which is reflected in the closely resembling spatial distributions of their r.m.s. values, or equivalently in the high value of the correlation coefficient F 1 which is approximately 0.75 near the pipe axis (see table 2) and is larger than 0.5 everywhere except within a thin wall layer. A characteristic feature of the buoyancy-generating mechanism is the presence of long intermittent plumes travelling along the pipe, which result in large values of the scalar integral length scale in the axial direction. These plumes also yield unusually large scalar-to-mechanical time ratios (T χ ) in the core, with values nearly twice as large as those observed in HS flows of in VCDHW flows but similar to those found in Rayleigh-Bénard convection at comparable Reynolds numbers. At small scale, the footprint of RT turbulence is found in the unusually large anisotropy exhibited by the axial Taylor microscales λ 1 (u i ) which are typically twice as large as their crosswise and spanwise counterparts.

w 2 /u 2 v 2 /u 2 R 13 F 1 F 3 P k / k T χ Present (Θ = 15 • ) S ≈ 6.0,
Inclination is responsible for both the mean shear ∂ z U and the crosswise concentration gradient ∂ z C. Combined with the axial gravity component, the latter induces a quadratic z-dependence of the shear stress u w . Similarly, combined with the small axial concentration gradient ∂ X C, the mean shear (more precisely the resulting linearly varying axial velocity) is responsible for the quadratic z-dependence of the crosswise mass flux c w . These two features make the large-scale structure in the present flow very different from that of pressure-driven flows (where u w varies linearly with the distance to the wall) and horizontal stably stratified shear flows, Rayleigh-Bénard or VCDHW convection (where c w is essentially constant outside the diffusive boundary layers). In particular, the quadratic variation of u w imposes that turbulence production is maximum in the core rather than in the near-wall region. As a substantial part of the transfer across scales is local, the mechanical dissipation rate mirrors this distribution, reaching its maximum on the axis and decreasing to a minimum at the edge of the core. Owing to the shear-induced anisotropy, energy redistribution toward the spanwise and crosswise velocity fluctuations is uneven, yielding v 2 -w 2 > 0 in line with known results in HS flows (see table 2). Finally, as the dominant contribution to the production rate of u w and c w is proportional to w 2 and transfers are mostly local in the core, w 2 is also forced to vary quadratically with z.

Up to this point, the wall has virtually no role in the flow structure apart from forcing the velocity components, hence the turbulent shear stresses and scalar fluxes, to vanish on it. However, owing to the kinematic blocking it imposes, it controls the redistribution process within a large part of the cross-section, say for r 0.2, reducing drastically the rate at which energy can be transferred from u to v and w as r increases. This influence results in a nearly circular distribution of the pressurestrain correlations Φ i in the outer part of the core, which in turn tends to impose a nearly isotropic decrease to v 2 and w 2 as r increases. The small values of the axial component Φ 1 at the edge of the core combine with the buoyancy-induced generation which directly injects energy into the axial fluctuation u to make the r.m.s. value of the latter reach a maximum there. In contrast to pressure-driven flows, the near-wall region, say 0.4 r 0.5, plays a minor role in the overall dynamics: it is essentially the layer where most turbulent energy is dissipated, and the place where the secondary flow is generated through the radial variations of the secondary Reynolds stresses inhomogeneity (Hallez & Magnaudet 2009).

Owing to the modest value of the shear Richardson number, the stabilizing stratification associated with the combined presence of the crosswise concentration gradient ∂ z C and the crosswise gravity component leaves the TKE balance nearly unaffected. However its damping effect on w is significant and table 2 suggests that w 2 /u 2 is lower by 15-25 % than in non-buoyant shear flows with a comparable shear number S = q 2 ∂ z U/ k . Although the correlation coefficients R 13 and F 3 are normalized quantities, i.e. the w damping is taken into account in their definition, they exhibit typical 10-25 % and 35-40 % reductions compared with values in non-buoyant flows, respectively. The reason for these reductions and for their difference in magnitude stems from the direct influence of stabilizing buoyancy effects in the u w and c w budgets, which was shown to be of O(Ri 1/2 L E /L B ) and O(L E /L B ) 2 , respectively, implying a stronger effect on the crosswise mass flux. Helped by the radial decrease of w , buoyancy is actually strong enough to reverse the sign of the production term at the edge of the core (r ≈ 0.35) in both budgets. This has little influence on the shear stress which is forced to nearly vanish there anyhow, owing to the vanishing of the mean shear. In contrast it has a quite spectacular effect on c w whose sign at larger r depends on the efficiency of the radial transport by turbulent fluctuations. From the z-axis up to the two bisectors of the (y, z) plane, the net transport of c w by the crosswise fluctuation w is strong enough to keep this sign unchanged. This is not the case around the y-axis, owing to the poor correlation between c w and the spanwise fluctuation v . This is why significant negative values of c w , i.e. countergradient fluxes, take place in the near-wall region comprised between the horizontal axis and the two bisectors.

As figure 9(h) revealed, turbulence production significantly exceeds dissipation throughout the core, with values ranging approximately from 1.25 to 1.75. Since axial and temporal variations of statistical quantities are negligibly small, this imbalance implies existence of radial turbulent fluxes directed towards the pipe periphery. A parallel can be made between this phenomenology and that of HS flows. Indeed, it is now established that under certain conditions the latter may reach an asymptotic self-preserving state where the small-scale structure (especially the Taylor microscales) and normalized large-scale quantities such as the Reynolds stress ratios presented in table 2 are independent of the streamwise position [START_REF] Tavoularis | Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence[END_REF][START_REF] George | The self-preservation of homogeneous shear flow turbulence[END_REF]. In contrast, the un-normalized Reynolds stresses, dissipation rate and integral length scales continuously grow downstream (exponentially) as a result of the imbalance P kk > 0. In the present case, such a growth is prevented by confinement which limits the lateral size of the energy-containing eddies, forcing all statistical moments to become independent of the axial position. Then the various radial fluxes resulting from turbulent transport are the direct counterparts of the streamwise transport contributions in HS flows, e.g. ∇ ⊥ • J k in (6.2) is the counterpart of U(z)∂(q 2 /2)/∂x. Hence, provided a Fickian relation between J k and q 2 approximately holds, a downstream growth of the TKE in the latter situation corresponds to a negative curvature of the q 2 -profile, i.e. a radial decrease of the TKE, in the present case. Here the Taylor microscales are almost uniform in the central region r 0.2 so that the present flow within this region may be seen as the confined counterpart of a HS flow, and the radial variations of the statistical moments observed therein are just the only way the flow can deal with the positive local imbalance between turbulence production and dissipation.

Final remarks

Data analysed in this paper revealed a rich phenomenology because the flow involves a remarkable combination of fundamental mechanisms, due especially to the interplay of negative and positive buoyancy and shear. Although we now have a pretty clear view of the resulting turbulent structure and how the various mechanisms fashion it, several important challenges remain.

Obviously one can consider increasing the Reynolds so as to reduce low-Reynoldsnumber effects and obtain a broader core. Increasing the spatial resolution by an order of magnitude is easily achievable and would make the analysis of the turbulent structure at Re λ = O(10 2 ) feasible. Increasing the Atwood number would allow non-Boussinesq effects, which are expected to induce top-down asymmetries within the pipe section, to be explored. Nevertheless we think the most urgent step should rather consider varying the pipe inclination, following the experimental investigations of [START_REF] Séon | Buoyant mixing of miscible fluids in tilted tubes[END_REF] and [START_REF] Znaien | Flow structure and momentum transport for buoyancy driven mixing flows in long tubes at different tilt angles[END_REF][START_REF] Riediger | Heat transport regimes in an inclined channel[END_REF] in the context of heat transport. It may be inferred from (4.1) and (4.2), (6.9) and (6.11) that ∂ z U, ∂ z C, u w and c w vary linearly with sin Θ for small inclination angles. Therefore, the behaviour of flows in pipes with 0 < Θ < 15 • may be expected to be intermediate between the one described here and the confined RT configuration considered experimentally by [START_REF] Cholemari | Axially homogeneous, zero mean flow buoyancy-driven turbulence in a vertical pipe[END_REF] and numerically by [START_REF] Schmidt | Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell[END_REF]. Even this purely RT case deserves more exploration since the available two studies provide only few detailed statistics. The situation is much less clear at larger inclinations, although the structure of the turbulent field is expected to be increasingly controlled by shear as Θ increases. However, the driving force at any angle is provided by the small axial density gradient which vanishes when the pipe becomes horizontal, except in the early stages which correspond to a lock-exchange configuration. Hence, one expects relaminarization to take place at some finite tilt angle whose value depends on the Atwood and Reynolds numbers. Investigating this relaminarization process should be the subject of an upcoming work.

Last, developing an 'engineering-oriented' model capable of predicting the mean parameters ∂ X C, ∂ z C and ∂ z U together with the main components of the diffusivity tensor D ij and the eddy viscosity ν T at any angle is worth the effort. An attempt based on the mixing length approximation was recently reported by [START_REF] Salort | Turbulent velocity profiles in a tilted heat pipe[END_REF] in the case of a prescribed axial density gradient. This model emphasizes the crucial role of the stabilizing stratification as Θ increases and introduces empirical corrections to the usual mixing length definition to cope with this influence. We believe that the approach involving the complete diffusivity and 'buoyancy' tensors introduced in § 6.4, supplemented with experimental and numerical results to specify the characteristic time and length scales of the dominant processes as a function of Θ, At and Re t , should help extend and rationalize this type of model. When ∂ X C is unknown as in the present case, the components of the diffusivity tensor may be used to close the transport equation (B 4) and determine the average density ρ (X, T). This approach will be the subject of a future paper. TABLE 3. Influence of grid resolution on some characteristic time-and volume-averaged statistics. As all horizontal directions are statistically equivalent, index h refers to any of them and one has q 2 = u 2 + 2 u 2 h and k = x + 2 h . Grid G1 is the one used throughout the paper. concentration at each time step. We selected a diameter-to-height aspect ratio of the pipe = 1/2 and kept the Reynolds number Re t equal to 886 as in the rest of the paper. We set the driving concentration gradient C x to 2.02 × 10 -2 , a value close to that found near x = 0 in the case the pipe is tilted with an angle of 15 • . Table 3 gathers the most significant volume-averaged statistical flow characteristics obtained by integrating the corresponding radial profiles. Grid G1 is that used in the rest of the paper. Grids G2 and G3 were obtained by doubling the resolution in the radial and azimuthal directions, respectively, whereas the axial resolution was doubled and quadrupled in grids G4 and G5, respectively. The Reynolds number Re u based on the fluctuating velocity defined in (3.1) and the Kolmogorov microscale η = Re -3/4 t k

-1/4 are used to characterize the strength of the large-scale motions and the size of the small-scale eddies, respectively. Table 3 indicates that, among grids G1-G5, these quantities vary by less than 4 %. The relative magnitudes of the horizontal component of the velocity variance, u 2 h , and dissipation rate, h , are used to qualify the flow anisotropy at large and small scales, respectively. These quantities experience negligible variations from one grid to another, except that the latter decreases by nearly 15 % on grid G2 which has the highest radial resolution and allows the large values of the radial velocity gradients between neighbouring small-scale upflowing and downflowing plumes to be slightly better captured in the core. We note in passing that these results display a strong qualitative resemblance with those of § 5 for the tilted case. Indeed, they show that the present flow is highly anisotropic at all scales, as 70 % of the energy and 53-60 % of the dissipation (depending on which grid is considered) are contained in the axial component. The next column considers specifically the ratio of the axial and horizontal Taylor microscales λ x (u ) and λ h (u h ). While this ratio should be unity in HIT, table 3 shows that it consistently reaches much higher values approximately 2.2 (±5 %) in the present flow, which again demonstrates a marked small-scale anisotropy between the axial and cross-sectional directions. The volume-averaged energy balance resulting from (6.2)-(6.5) results in the constraint -2 c u = k , provided that the flow is statistically stationary. Hence the Kolmogorov microscale can also be inferred from the axial buoyancy flux, yielding η c u = Re -3/4 t (-2 c u ) -1/4 . The last column shows that the ratio η c u / η only departs from unity by approximately 5 % on grid G1 and gets closer to one as the grid is refined by a factor of two in each direction. The results in this last column are important because they compare a large-scale estimate of the Kolmogorov microscale, η c u based on c u which is dominated by the contributions of large-scale density fluctuations, and η which is directly related to k and hence is determined by the small-scale velocity gradients. That the two definitions yield consistent estimates indicates that the numerical dissipation, which (together with the residual time dependence) is responsible for their difference, is only a small fraction of k .

We also performed an additional run with an = 1/4 geometry using the same resolution as in grid G1, i.e. a 64 3 grid in that case, other things being equal. We found the Reynolds number to be Re u = 219 and the Kolmogorov microscale to be 0.013, i.e. the former (respectively latter) is somewhat lower (respectively larger) than with = 1/2. The corresponding ratio η c u / η is 1.018, i.e. it is closer to unity than with grid G1. Not surprisingly this test shows that, for a given resolution, the above constraint is better satisfied when the turbulent intensity is reduced. Overall, this convergence study indicates that, although some variations are observed on small-scale quantities when the grid is refined, the results obtained on grid G1 are robust, both at large and small scales, a robustness that includes small-scale anisotropy properties. Moreover, it must be kept in mind that in the case the pipe is tilted by 15 • , the turbulent Reynolds number is Re u = 185, i.e. it is 20 % lower than in the above test case, and the Kolmogorov microscale ranges from 0.015 on the axis to 0.018 for r = 0.4, i.e. it is at the very least 25 % larger than in the convergence study. Therefore, data obtained with grid G1 which we discuss throughout the paper are expected to have an even better accuracy than those obtained in this test case.

Finally, to make sure that numerical dissipation does not subtly mitigate or even hide the effects of the axial under-resolution, we first computed the one-dimensional energy spectra for the five grids G1-G5 and averaged them over all (y, z) positions within the cross-section to improve statistical convergence. All spectra (not shown) display an excellent collapse for axial wavenumbers k x k M = 20. Beyond that wavenumber, spectra obtained with grids G1-G3 (N x = 32) depart from those obtained with the refined grids G4 and G5. On the standard grid G1, wavelengths 4 x and 5 x ( x being the axial grid spacing) correspond to k x ≈ 25.1 and k x ≈ 20.1, respectively, which implies that eddies of size l M = π/k M = 2.5 x are well resolved while those with l = 2 x are already affected by numerical errors. Since the ratio l M /∆ x only depends on the numerical schemes and not on the physical configuration, we conclude that in the tilted case numerical errors start to manifest themselves for eddies whose size is approximately 9-10η. Then we computed the one-dimensional spectra E i (k x ) corresponding to each velocity component u i in the tilted configuration at the same position as the two-point correlations displayed in figure 8(a). For each of these spectra (not shown), we evaluated the integrals

I i = k max 0 k 2 x E i (k x ) dk x and J i = k M 0 k 2 x E i (k x ) dk x ,
where k max = π/ x. The difference ∆ i = I i -J i represents the contribution to (∂u i /∂x) 2 of the axial wavenumbers belonging to the range k M k x k max . Only a fraction of this contribution results from numerical errors but, since no convergence study was carried out on the tilted configuration, this fraction is unknown. Hence, ∆ i provides a crude overestimate of the numerical error due to the axial under-resolution. Using the definition (5.1) of the Taylor microscales and considering that u 2 i itself is almost unaffected by the axial resolution (as confirmed by table 3 in the test problem), the relative contribution of the range k M k x k max to λ x (u i )/λ x (u i ) is half ∆ i /I i , from which we could conclude that λ x (u )/λ x (u ) ≈ 8 %, λ x (v )/λ x (v ) ≈ 11.5 % and λ x (w )/λ x (w ) ≈ 13 %. Standard uncertainty evaluation then indicates that the relative contribution of these axial wavenumbers to 1 (through the term (∂u /∂x) 2 ) is ≈1.2 %. We performed similar estimates with 2 and 3 and concluded that the corresponding contributions are ≈7.6 % and ≈9.2 %, respectively. The latter two estimates were obtained by noting that, owing to continuity, axial derivatives are not only involved in (∂v /∂x) 2 and (∂w /∂x) 2 but also in (∂v /∂y) 2 and (∂w /∂z) 2 in such a way that the relative contribution of the range k M k x k max to λ y (v ) and λ z (w ) is at worst the same as that to λ x (u ). Finally, gathering previous estimates, we found that the relative contribution of that wavenumber range to k is approximately 4.3 %. Although the above percentages are not totally negligible, they are small enough to conclude that the spectral range affected by the axial under-resolution, hence the associated numerical dissipation, only marginally contributes to the small-scale statistics discussed in § 5.

  FIGURE 1. (Colour online) Flow in the tilted pipe: sketch of the global geometry (left) and snapshot of isopycnals surfaces in a short median part of the vertical diametrical plane (y = 0, |x| 5.5d) in the final stage of the simulation (right).

  FIGURE 2. (Colour online) Isocontours of mean concentration and velocity profiles within the pipe cross-section: (a) C; (b) U with vectors showing the secondary flow; (c) V;(d) W. In all figures, the region corresponding to the small white central disc is excluded because the singularity of the cylindrical coordinate system at y = z = 0 frequently induces inaccuracies in the interpolations (especially the differentiations) involved in the post-processing.

  FIGURE 3. (Colour online) Root mean square values of the velocity and concentration fluctuations within the pipe section: (a) c 2 1/2 ; (b) u 2 1/2 ; (c) v 2 1/2 ; (d) w 2 1/2 .

FIGURE 6 .

 6 FIGURE 5. (Colour online) Contributions to the axial mean momentum balance within the cross-sectional plane: (a) transport by the secondary flow; (b) contribution of Reynolds stresses; (c) viscous diffusion.

  FIGURE 7. (Colour online) Contributions to the mean concentration balance in the cross-sectional plane: (a) axial transport; (b) transport by the secondary flow; (c) turbulent transport.

  FIGURE 8. (Colour online) Two-point correlations and Taylor microscales. (a)Longitudinal two-point correlations resulting from an average taken over two symmetric arcs of a circle whose middles are located on the z-axis at z = ±0.25 and whose length covers a 90 • angle in the azimuthal direction (corresponding to 16 grid points): solid line, ρ 11 ; dash-dotted line, ρ 22 ; dashed line, ρ 33 ; solid line with crosses, ρ cc . (b) Profiles of the Taylor microscales along the z-axis. From top to bottom in the core, the solid lines refer to λ x (u ), λ x (v ), λ x (w ), the dashed lines to λ y (v ), λ y (u ), λ y (w ) and the dash-dotted lines to λ z (w ), λ z (u ), λ z (v ), respectively.

  FIGURE 9. (Colour online) Contributions to the production/destruction of the turbulent energy within the pipe section: (a) P Uz ; (b) P Uy ; (c) P sec ; (d) P F 1 ; (e) P F 3 ; (f ) TKE production, P k = P Uz + P Uy + P sec + P F 1 + P F 3 ; (g) dissipation rate k (for r 0.47, the near-wall values of k are out of scale, the peak value at the wall being 0.068); (h) ratio P k / k of TKE production over dissipation.

  FIGURE 10. (Colour online) Distribution of the three diagonal components of the pressurestrain correlation φ i in the cross-sectional plane. (a) Φ 1 ; (b) Φ 2 ; (c) Φ 3 .

  FIGURE 11. (Colour online) The three leading terms in the concentration variance budget. (a) the axial production term P c 2 X ; (b) the crosswise production term P c 2 z ; (c) the scalar dissipation rate χ.

  FIGURE 12. (Colour online) The total production/destruction contribution to the budget of (a) the shear stress u w , (b) and (c) the concentration fluxes c u and c w , respectively.

  FIGURE 13. (Colour online) Distribution of (a) the eddy viscosity, (b) the eddy diffusivity and (c) the turbulent Schmidt number in the pipe section; white zones correspond to subregions where the quantity under consideration is out of scale or undefined owing to the smallness or even vanishing of the corresponding mean gradient (in a and b), or where the eddy diffusivity is near-zero or negative (in c).

  The time scale ratio S, sometimes referred to as the shear number, is defined as S = q 2 ∂ z U/ k . confinement. The most salient of these features are summarized below. Some large-scale statistical data available in the literature for non-buoyant or stably stratified HS flows, VCDHW and confined RT flows, are gathered in table 2 and help discuss the influence of some of the above mechanisms in the present case.

	Ri ≈ 0.046 0.30 0.42 0.36 0.35 HS4 (Θ = 90 • ) S ≈ 5.0, Ri = 0.13 HS1 (averaged), Ri = 0.0 HS2 S ≈ 12.5, Ri ≈ 0 HS3 S ≈ 5.0, Ri = 0.0 0.27 HS5 S ≈ 11.0, Ri = 0.0 0.32 HS6 (Θ = 90 • ) S ≈ 11.0, Ri = 0.058 0.22 VCDHW (Θ = 0 • ) S ≈ 7.1 0.42 0.50 RT (confined, Θ = 0 • ) TABLE 2. Comparison of some statistics in the core with some reference data from 0.47 -0.36 -0.75 0.29 1.37 0.51 0.53 -0.48 --1.29 -0.53 -0.45 -0.59 0.45 1.75 0.34 0.58 -0.49 -0.60 0.47 1.55 0.28 0.54 -0.37 -0.60 0.32 1.20 0.20 ? -0.52 -0.70 0.47 1.87 -? -0.47 -0.70 0.40 1.38 ? 0.58 -0.40 -0.57 0.45 1.25 0.20 0.50 0 ? 0 ? ? HS flow, VCDHW or confined RT flows, in the limit of the Boussinesq approximation (At 1). Data are from Tavoularis & Karnik (1989) (HS1), Tavoularis & Corrsin (1981a) (HS2), Kaltenbach et al. (1994) (HS3 and HS4), Holt, Koseff & Ferziger (1992) (HS5 and HS6), Versteegh & Nieuwstadt (1998) (VCDHW) and Cholemari & Arakeri (2009) (RT), respectively.

† Email address for correspondence: magnau@imft.fr

Acknowledgements

Part of this research was funded by the Agence Nationale de la Recherche in the framework of the GIMIC project under grant No. ANR-07-BLAN-0181. The HPC resources were provided by CALMIP and GENCI-CINES under grant 2009-025035. We acknowledge the invaluable help of Annaïg Pedrono with the optimization of the JADIM code.

Supplementary movies

Supplementary movies are available at http://dx.doi.org/10.1017/jfm.2014.638.

Appendix A. Grid convergence

To assess the quality of the numerical data, we performed a grid convergence study not directly on the flow discussed in the paper but on a closely related one in which converged statistics may be obtained at a much more reasonable cost. In this problem, which was recently considered by [START_REF] Schmidt | Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell[END_REF] in the context of Rayleigh-Bénard convection at unit Prandtl/Schmidt number, the pipe is vertical and the flow is driven by a prescribed axial density gradient, i.e. we have C(x, y, z, t) = x C x + c (x, y, z, t). Periodic boundary conditions are imposed on velocity and concentration fluctuations on the top and bottom sections, while the flow still obeys a no-slip condition on the lateral wall. In this situation, statistical homogeneity in horizontal planes allows spatial averaging to be performed over both axial and azimuthal positions, making it possible to obtain well-converged statistics in much less time than in the tilted configuration. Since the pipe is not closed, the mean flowrate U is not necessarily zero, which may cause a drift in the cross-sectional averaged density owing to the source term U C x in the C-balance. To prevent the occurrence of such features, we followed the strategy employed by [START_REF] Schmidt | Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell[END_REF], i.e. we removed the cross-sectional averages of the axial velocity and Grid number Given the definitions of the space-time averages introduced in § 2.2, the incompressibility condition implies

where V ⊥ = (V, W) denotes the secondary mean flow and

is the two-dimensional gradient within the cross-sectional plane. Note that although the secondary velocities are of O(Re -1/2 t ) (Hallez & Magnaudet 2009), their gradients in the y and z directions may locally be of O(1) because they are confined within the near-wall region whose thickness is of O(Re -1/2 t ). Hence, only the combination of the two gradients involved in the two-dimensional divergence is of order . Integrating (B 1) over any cross-section implies ∂ X U = 0, which indicates that the averaged flow rate does not vary along the pipe. Given the presence of the end walls this flow rate is zero, so that U = 0. (B 2)

Using (B 1) and invoking the incompressibility of the instantaneous velocity field, the mean concentration balance may be written in the form

where v ⊥ = (v , w ) denotes the projection of the velocity fluctuation within the crosssectional (y, z) plane. At this point it is convenient to write C(X, y, z, T) = C (X, T) + C(X, y, z, T), where the deviation C represents the mean stratification in the (y, z) plane. Noting that U C = 0 and integrating within that plane yields

which is the balance governing the changes of the average mixture density. The difference between (B 3) and (B 4) may then be written in the form

This balance governs the variations of the turbulent flux c v ⊥ within the pipe cross-section. If we provisionally set = 0 in (B 5), only the transport term associated with the secondary flow and the last term in the right-hand side subsist. Furthermore, if we consider a two-dimensional geometry in which the mean flow is purely axial, (B 5) reduces to ∂ z c w = 0 which, owing to the vanishing of the velocity fluctuation at the wall, yields c w = 0 everywhere. Hence at least one O( ) term is required to avoid this unphysical degeneracy. In the situation considered in this paper, Re 1/2 t ≈ 0.15 so that the whole set of terms within the first two parentheses (which is of O( )) may be neglected compared with the third term of the left-hand side (which is of O(Re -1/2 t )). Moreover, assuming the generalized Fickian relationship 2) and O( C), both of which are also negligibly small compared with the third term. Therefore, it may be concluded that (B 5) reduces at leading order to

The first term in the left-hand side of (B 6) may be thought of as a forcing resulting from the transport of the density variations of the mixture along the pipe. Owing to the existence of a constant mean shear, the magnitude of this forcing varies linearly with z, which in turn forces the turbulent fluxes to vary quadratically.

The average density ρ at the axial position x = X/ is ρ (X, T) = 1 + 2At( C (X, T) -1/2). Thus, we define the hydrostatic contribution P(x, X, z, T) to the pressure field P as

and the average pressure deviation as P(X, y, z, T) = P(x, X, y, z, T) -P(x, X, z, T).

With this definition, the axial variations of P take place over O( -1 ) distances (i.e. P depends on X but not on x), in contrast with those of P and P, because the mean effect of the axial gravity component has been removed. The mean momentum balance in the streamwise direction then reads

Terms within parentheses are of O( ) or even smaller. In contrast with (B 5), these terms can all be neglected because the first term in the right-hand side provides an O(1) forcing. Transport terms associated with the secondary flow cannot be neglected because the axial velocity U falls from its O(1) maximum to zero in the O(Re -1/2 t )-thick wall region, so that the corresponding velocity gradients are of O(Re 1/2 t ); this remark also implies that the viscous term is significant within this wall region. Therefore, at leading order

Again, the linear z-variation of the stabilizing stratification yields a quadratic variation of the Reynolds stress u w .

Equations (B 3) and (B 8), or their approximate forms (B 6) and (B 9), are the key equations of the mean flow. Comparison of the latter two emphasizes the fact that the mean momentum balance only involves processes taking place within the crosssectional (y, z) plane, while the mean concentration balance crucially depends on the axial variations of the mean concentration, i.e. of the mixture density.

Appendix C. A simple algebraic model for the turbulent fluxes in the presence of buoyancy effects Considering the form of the production term in (6.8), we assume that the scalar variance c 2 obeys

Neglecting turbulent transport, T c may be thought of as the time scale such that the scalar dissipation χ equals 2c 2 /T c . Injecting (C 1) into (6.10) and (6.11) yields the generic form of the production term in the c u i balance as

where the 'buoyancy' tensor B is defined as

In the present case, the two positive components B i1 may be thought of as the square of the growth rate of concentration disturbances resulting from the negative buoyancy due to ∂ x C, while the two negative components B i3 are, in absolute value, the square of the characteristic frequencies of oscillations associated with the positive buoyancy due to ∂ z C. In particular, -B 33 is the square of the Brunt-Väisälä frequency N. Now, let us write the c u i -budget in the form P ci + Γ ci = 0 (in HS flows, Γ ci is dominated by the contribution -c ∂ i p ). Following [START_REF] Rogers | An algebraic model for the turbulent flux of a passive scalar[END_REF], we assume that Γ ci may be written in the form Γ ci = -c u i /τ , where τ is a second time scale characterizing the relationship between the scalar fluxes and their production rates. Then the c u i balance formally reduces to the algebraic equation

with K ik = ∂ z Uδ i1 δ k3 + δ ik /τ -T c B ik . Buoyancy effects render (C 4) nonlinear with respect to concentration gradients. However, we are interested in predictions of the diffusivity tensor for specified values of the shear stress and of all components B ik , so that (C 4) is actually a 3 × 3 linear system which can be inverted to obtain the coefficients of the diffusivity tensor in the form D ij = A ik u k u j . This yields

). The non-buoyant case is recovered by setting T c = 0 in (C 5), which yields equations (3.12)-(3.16) of [START_REF] Rogers | An algebraic model for the turbulent flux of a passive scalar[END_REF]. Defining

the explicit form of the diffusivity tensor in the present flow writes