Mathilde Tavares
email: mathilde.tavares@u-pem.fr

Stéphane Vincent

Meryem Ould-Rouiss

Jean-Luc Estivalezes

Stephane Vincent

Stephane Vincent

A priori

Introduction

Two-phase flows are involved in many industrial and environmental applications such as spray formation, wave breaking or oil transportation. The modeling and simulation of complex interactions between turbulence and interface in multiphase flows remain challenging. As in single phase flow, large Eddy Simulation (LES) solve the problem of large scale while smaller subgrid scales are computed through physical or mathematical models in two-phase flows, involving this way several unclosed subgrid terms characteristic of turbulent and interfacial motion subgrid correlations. Several works [START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase les filtering[END_REF], [START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], have established a hierarchy of unclosed subgrid terms through a priori studies to improve LES simulation of two-phase flows. The main objective of the present work is to characterize the hierarchy of subgrid terms for two-phase flows in the case of a phase inversion flow. In the following sections the 1-fluid model is presented and filtered to highlight the unclosed LES terms for two-phase flows. Then, we present the reference simulation, a phase separation flow from a numerical benchmark [START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF]. Finally the results of a priori tests are discussed to evaluate the relative importance of LES subgrid terms in the reference case.

Numerical modeling

Navier Stokes equations/1-fluid model for two-phase flows

The two-phase flow numerical modeling is restricted to immiscible, incompressible and isothermal fluids. For two non miscible fluid where an index k refers to one phase if k = 1 and the other phase if k = 2, the dynamic of each phase k is governed by the Navier-Stokes equations. Thanks to jump relations, mass and momentum conservation is ensured at the interface. By adding these boundary conditions to NavierStokes equations, it is possible to build in the framework of the 1-fluid formalism of [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF] a unique set of equation valid in both phase by introducing the volume fraction function C which is equal to 1 in one phase an 0 elsewhere and describes the interface evolution through a material advection equation (1.1).

∇.u = 0 (1.1a) ∂ ρu ∂t + ∇.(ρu ⊗ u) = -∇P + ρg + ∇.(2µS) + σ κn i δ i (1.1b) ∂C ∂t + u.∇C = 0 (1.1c)
u is the velocity, P the pressure, t the time, g the gravity vector, ρ the density, µ the dynamic viscosity. S k = 1 2 (∇u k + ∇ T u k ) is the viscous stress tensor. σ κn i is the surface tension force where σ denotes the surface tension coefficient, κ the local curvature and n i the normal to the interface. Local quantities such as density ρ or viscosity µ can be defined by using the volume fraction function C.

ρ = ρ 2 + (1 -C)ρ 1 (1.2a) µ = µ 2 + (1 -C)µ 1 (1.2b)
The curvature κ depends on the interface topology. It can modeled using the Continuum Surface formulation of [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]:

κn i δ i = ∇. ∇C ||∇C|| ∇C, n = ∇C ||∇C|| , n i δ i = ∇C (1.3)

Filtered Navier-Stokes equations for two-phase flow

LES formalism required the application of a low-pass frequency filtering to the field to resolve large turbulent scales while the effects of smaller scales are modeled [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF].

The filtering operation of a hydrodynamic field is defined on a domain Ω by:

Φ(x,t) = G * Φ → Φ(x,t) = Ω t -∞ G( ∆ (x,t), x -x ,t -t )Φ(x ,t )dx dt (1.4)
∆ is the cutoff lengthscale of the filter, G a formal operator. It can be assumed that the filtering operator commutes with time and spatial derivatives. Then, the volume filtering of incompressible Navier-Stokes equations depend on { ū, P, C, ρ, μ, κ} and can be written as:

∇. ū = 0 (1.5a) ∂ ρ ū ∂t + ∇.( ρ ū ⊗ ū) -ρg + ∇ P -∇.(2 μ S) -σ κ∇ C = -∇.(τ conv -τ di f f ) + τ super f - ∂ τ temps ∂t (1.5b) ∂ C ∂t + ū.∇ C = τ inter f (1.5c)
This approach is divergence free but required modeling of 5 subgrid terms. It can be more convenient to use the Favre average or mass weighted filtering for the velocity field such as ũ = ρu ρ (

This approach is based on { ũ, P, C, ρ, μ, κ}. With this formulation, the velocity field is no more divergence free with the additional source term τ inter f (ρ 2 -ρ 1 ) ρ but only 4 subgrid terms require modeling. A summary of the subgrid term for the classic volume filtering and the phase weighted filtering are given in the table (1.1).

Subgrid terms

Volume filter Mass weighted filter (Favre)

∂ τ temps /∂t ∂ ρu ∂t -∂ ρ ū ∂t (ρ 2 -ρ 1 ) ũτ inter f ∇.τ conv ∇.(ρu ⊗ u -ρ ū ⊗ ū) ∇.( ρ( u ⊗ u -ũ ⊗ ũ)) ∇.τ di f f ∇.(2µS -2µS) ∇.(2µS -2 μ S) τ super f σ (κ∇C -κ∇ C) σ (κ∇C -κ∇ C) τ inter f u.∇C -ū.∇ C u.∇C -ũ∇ C
Table 1.1: Subgrid terms for two-phase flow 1.3 A priori filtering of phase separation flow

Description of the problem

Our study is based on the numerical benchmark of [START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF] where a parametric study on macroscopic quantities such as kinetic or potential energy have been performed with different codes and meshes to characterize vorticity generation on a phase separation flow between water and oil (1.1a). We have chosen to deal only with the turbulent Table 1.2: Fluids characteristic for the phase inversion problem results of DyJeat (Dynamic of Jet ATomization), a computational in-house code developed at ONERA/DMAE. DyJeat is implemented with the finite volume method for the discretization of Navier-Stokes equations on a staggered mesh.The Ghost-Fluid method is used to deal with surface tension forces, density and viscosity jumps at the interface. The Level-Set approach is also used for tracking interfaces and a Level-Set/VOF coupling is ensured to improve mass conservation. Details of implementation and validation can be found in [START_REF] Vincent | A phase inversion benchmark for multiscale multiphase flows[END_REF]. All the fields have been taken for a single time at the enstrophy peak (1.1b) corresponding to the time where the turbulence is maximum. Three meshes (128 3 , 256 3 , 512 3 ) have been used in order to study the effect of the mesh on the subgrid contributions. Moreover, we have used two different filter size, G2 and G4, and two filter's type, the volume filtering and the mass weighted filtering to evaluate the influence of the type and the size of the filter. The size of the filter G2 is for instance: In agreement with [START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], [START_REF] Vincent | Numerical simulation of phase separation and a priori two-phase les filtering[END_REF], our a priori analysis shows that the subgrid viscous term τ di f f is always negligible. Besides, τ super f , the subgrid contribution linked to the surface tension, is the dominant term. Indeed, capillarity drives the phase inversion mechanism at the enstrophy peak with important rupture and coalescence phenomena, increasing the contribution of τ super f . The convective τ conv and interfacial τ inter f subgrid terms have the same order of magnitude but remain significant. We observe a strong influence of the mesh on unclosed LES terms. Indeed, all the subgrid contributions decrease when the mesh size becomes smaller except for τ super f . On the other hand, the figure (1.3b) shows that subgrid contributions are bigger when we use volume filtering than phase weighted average as shown in literature. According to [START_REF] Labourasse | Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests[END_REF], using volume filtering process overestimates the convective transfer across the interface whereas phase weighted underestimates it. The figure (1.3a) shows that the filter size (G2 and G4) has almost the same effect than the mesh size on the subgrid contributions and for the smaller filter size, unclosed LES terms are weaker even for τ super f in our case.

G (2) (Φ i jk ) = x i+1 x i-1 y j+1 y j-1 z k+1 z k-1 Φ(x,

Conclusion

In this study, we have shown that the contribution linked to the surface tension has to be taken into consideration. In the case of the phase inversion flow, at the enstrophy peak, τ super f and τ inter f has been the predominant term, whatever the mesh considered. It's worth noting that the mesh size here is four times smaller than those studied in literature whereas the Weber is six times smaller. The capillarity effects are then more important and better resolved with our finer grid, increasing the part of the capillary subgrid term. To conclude, this term cannot be neglected and it is then necessary to model it properly to improve LES multiphase flow simulation, in addition to inertial terms. It has to be noted that the τ inter f is also not negligible and can lead to non divergence free filtered LES fields.

Fig. 1 . 1 :

 11 Fig. 1.1: Initial conditions of phase inversion flow problem (1.1a) and DNS solution at enstrophy peak (1.1b)

3 Fig. 1 . 2 :

 312 Fig. 1.2: Subgrid term hierarchy for meshes 128 3 (1.2a) and 512 3 (1.2b) with volume filtering filter along x, y, z directions

Fig. 1 . 3 :

 13 Fig. 1.3: Subgrid term hierarchy for meshes 512 3 depending on the size of filter (1.3a) and the type of the filter (1.3b) along z direction


A Revell

I Afgan

A E A Ali

M Camps Santasmasas

T Craft

A De Rosis

J Holgate

D Laurence

B E O Iyamabo

A Mole

B Owen

M Savoie

A Skillen

J Wang

X Zhang

Coupled Hybrid RANS-LES Research at The University of Manchester

This short communication summarises recent developments and ongoing work on Hybrid RANS-LES Methods (HRLM) at The University of Manchester.

Our work is focussed primarily on the development of coupled zonal methods, investigating improvements and novel applications of embedded turbulence simulation via the Synthetic Eddy Method. These activities are centred around the premise that no single modelling methodology will be generally applicable. Not only does this apply to different test cases, but for particular regions of the flow within a single test case. Similarly, when factoring in the computational resources available to the modernday CFD practitioner, a framework for multifidelity approaches seems pragmatic. The work presented herein is intended primarily as a window through which we outline our current activities, mostly driven by PhD projects, many of which with industrial projects in mind. We also introduce recent developments focused on using the lattice Boltzmann method as part of a coupled Hybrid RANS LES approach.

Introduction

Hybrid RANS-LES Methods offer compromise to the accuracy-cost-speed trilemma for simulations at moderate-to-high Reynolds number, and can generally be separated into two distinct classes, 1) Zonal and 2) Non-Zonal (Global). The latter has become more widely used than the former over the past decades, in particular due to the prevalence of Detached Eddy Simulation proposed by Strelets et al. [START_REF] Strelets | Detached eddy simulation of massively separated flows[END_REF], which has received a lot of attention for industrial applications [START_REF] Haase | DESider-A European Effort on Hybrid RANS-LES Modelling: Results of the European-Union Funded Project[END_REF]. On the other hand, Zonal approaches such the framework proposed by Deck [START_REF] Deck | Zonal-detached-eddy simulation of the flow around a high-lift configuration[END_REF] can offer different advantages. These methods generally employ sub-domains which are then solved individually with a separate set of transport equations. Unlike non-zonal methods, there is a clear and conscious demarcation between mean RANS and transient LES regions. In the transfer to LES regions, instantaneous fluctuations must be generated using information taken from the upstream flow; whereas the mean/modelled turbulent quantities must be recovered from the LES field upon re-entering a RANS region.

At Reynolds numbers of practical engineering interest, the incoming flow will often be turbulent. The generation of turbulent inflow boundary conditions for scale-resolving simulations poses a formidable challenge. A popular solution to this problem is the Synthetic Eddy Method (SEM) of Jarrin et al. [START_REF] Jarrin | A synthetic-eddy-method for generating inflow conditions for large-eddy simulations[END_REF] whereby Lagrangian spots (or synthetic eddies) are generated to represent fluctuations, with the intensity and sign set to satisfy prescribed first and second order statistics. In practise, these low-order statistics may be obtained by experiment, theoretical approximation, or as is perhaps the most common, by a separate or precursor RANS study. Embedded LES (ELES) is a natural evolution of this concept, whereby the RANS data generated in one mesh is used directly for the generation of synthetic turbulence and subsequent turbulent simulation in another. An early implementation of ELES is presented by Cokljat et al. [START_REF] Cokljat | Embedded LES Methodology for General-Purpose CFD Solvers[END_REF] who considered a number of internal flow applications and a number of more recent implementations have followed. Notably Li et al. al. [START_REF] Li | Validation of a Window-Embedded RANS/LES Method Based on Synthetic Turbulence[END_REF] validated a two-region ELES solver consisting of a single RANS and implicit LES region and Anupindi & Sandberg [START_REF] Anupindi | Implementation and Evaluation of an Embedded LES-RANS Solver[END_REF] implemented an ELES model within OpenFOAM. In a slightly different approach, Vonlanthen et al. [START_REF] Vonlanthen | Assessment of a one-way nesting procedure for obstacle resolved large eddy simulation of the ABL[END_REF] introduced a oneway nesting procedure which embeds a highlyresolved LES within a low-resolution LES. For an extensive review of Embedded LES approaches the reader is referred to Holgate et al. [START_REF] Holgate | A review of embedded large eddy simulation for internal flows[END_REF].

Through work at the University of Manchester, the SEM has evolved in several directions; a divergence-free version was proposed by Poletto et al. [START_REF] Poletto | A New Divergence Free Synthetic Eddy Method for the Reproduction of Inlet Flow Conditions for LES[END_REF] and a more efficient formulation was presented by Skillen et al. [START_REF] Skillen | Accuracy and efficiency improvements in synthetic eddy methods[END_REF]. More recent work has focused on the codecode coupling aspects, wherein mixed fidelity solvers can be used interchangeably. In particular we have developed ideas based on the Two-Velocity method of Uribe et al. [START_REF] Uribe | Computation of flow in a 3d diffuser using a two-velocity field hybrid rans/les[END_REF], and the consistent Dual-Mesh approach of Xiao and Jenny [START_REF] Xiao | A consistent dualmesh framework for hybrid LES/RANS modeling[END_REF]. These methods solve RANS and LES models simultaneously on separate grids, with added 'drift terms' in the momentum, pressure and turbulent transport equa- tions used to relax the flow statistics towards a consistent solution.

As a final point, it is worth noting that in all the above, zonal methods are enhanced when the upstream prediction of turbulent Reynolds stresses is good. As such there is a rational case for the use and further refinement of advanced turbulence modelling tools in this framework. We have been using the Elliptic Blending model of Manceau and Hanjalić [START_REF] Manceau | Elliptic blending model: A new near-wall reynolds-stress turbulence closure[END_REF] with its demonstrated advantages complex flows [START_REF] Billard | Application of recently developed elliptic blending based models to separated flows[END_REF], as well as models based on the stressstrain lag approach of Revell et al. [START_REF] Revell | Turbulence modelling of unsteady turbulent flows using the stress strain lag model[END_REF], recently reformulated for elliptic blending [START_REF] Lardeau | Development of an elliptic-blending lag model for industrial applications[END_REF]. This is highlighted in recent work demonstrating potential for multiple embedded regions, as shown in 1, where faster recovery is observed with the more advanced turbulence model.

Recent work with the SEM

As recently demonstrated by Skillen et al. [START_REF] Skillen | Accuracy and efficiency improvements in synthetic eddy methods[END_REF], with improved normalisation of the fluctuations on the inlet plane, the Synthetic Eddy Method (SEM) can accomplish a short recovery of the turbulent statistics. However, it remains sensitive to the prescribed variation of the length scale σ to produce realistic turbulence. The length scale of the eddies is generally assigned the following form where C σ is a scaling factor, discussed further below.

σ i = min(C σ R 3/2
ii /ε , σ max ). Defining the term Eddy Simulation (ES) as applying to both DES and LES methods, and following extensive testing for internal flows, we now consider application to a spatially developing boundary layer case in the fully turbulent regime. A schematic showing how the domain is split between mean and resolved (ES) domain is shown in 2. In this case we use Improved Delayed Detached Eddy Simulations (IDDES) in place of the LES in order to further reduce the computational cost according to Shur et al. [START_REF] Shur | A hybrid rans-les approach with delayed-des and wall-modelled les capabilities[END_REF]. The Reynolds number based on the momentum thickness at the start of the LES region is Re θ = 3040. In order to assess the development length, the evolution of the skin friction coefficient, C f , along the length of the ES domain is shown in 2. Some sensitivity to imposed length scale variation is observed as also identified in ear-lier work for internal flows [START_REF] Holgate | A review of embedded large eddy simulation for internal flows[END_REF], in an attempt to minimise the development length. In this instance increasing the length scale by setting C σ = 2 improved the development length. The figure also demonstrates that using an anisotropic length scale definition is a significant improvement over an isotropic definition.

First and second order statistics of the flow at a location downstream of the SEM were also assessed. These are given for a Reynolds number based on the momentum thickness of Re θ = 5200 and compared to experimental data of DeGraaff & Eaton [START_REF] Degraaf | Reynolds number scaling of the flat-plate turbulent boundary layer[END_REF]. Although a degree of log-layer mismatch is present, the close agreement of the statistics suggest that the relationship between the IDDES and the SEM is working. As an example of a more complex test for SEM, we consider the case defined by Liu et al. [START_REF] Liu | Interaction between a spatially growing turbulent boundary layer and embedded streamwise vortices[END_REF] for the superposition of a pair of counterrotating vortices in a turbulent boundary layer, directly relevant to flow downstream of vortex generators and other flow control devices. Initial work has been performed on a 'Common Flow Down' vortex pair, generated by superimposing Batchelor vortices onto the boundary layer inlet before applying the SEM. 4 shows a visualisation of the flow for this case using Q-iso surfaces.

Dual Mesh with Heat transfer

In this work, the dual-mesh approach [START_REF] Xiao | A consistent dualmesh framework for hybrid LES/RANS modeling[END_REF] was extended and applied for the first time to natural convection flows, following work by Tunstall et al. [START_REF] Tunstall | CFD and turbulence modelling for Nuclear Plant Thermal hydraulics systems[END_REF]. The two cases reported here are a high Rayleigh number differentially heated square cavity flow and the flow in a cylindrical annuli. Detailed descriptions of these test cases can be found in [START_REF] Sebilleau | Direct Numerical Simulation of an air-filled differentially heated square cavity with Rayleigh numbers up to 1011[END_REF] and [START_REF] Addad | Quasi-DNS of natural convection flow in a cylindrical annuli with an optimal polyhedral mesh refinement[END_REF].

In the dual-mesh approach, two separate meshes are produced, tailored to the specific needs of each solver; the RANS grid is refined in the near-wall regions and the LES grid is designed to be under-resolved near the walls of the domain. The two simulations are 'drifted' towards each other using source terms added to the momentum, temperature and turbulence equations of the RANS and the LES. A novel criterion is introduced to predict the locations at which the LES is under-resolved. At these locations, the mean velocity and temperature, total turbulent kinetic energy and total temperature variance (the latter two including resolved and modelled contributions) of the LES solution are corrected towards the corresponding RANS quantities. At remaining locations, the RANS solution is corrected towards the solution of the LES. More details about the formulation of the method can be found in [START_REF] Xiao | A consistent dualmesh framework for hybrid LES/RANS modeling[END_REF] and [START_REF] Tunstall | CFD and turbulence modelling for Nuclear Plant Thermal hydraulics systems[END_REF].

A new criterion has been introduced based on the ratio of the turbulence length scales to the grid size; designed to account for the presence of both turbulent and laminar regions within the flow domain. Reasonable results for the square cavity and cylindrical annuli flows were obtained when the new criterion was used. Results for the square cavity are shown in 5 which displays the LES zone weight (σ L ) as well as the instantaneous velocity. This quantity is equal to 0 at locations where RANS drives the LES, and equals unity at all other locations. It can be observed that the new criterion returns a rational variation, with σ L =0 restricted to the near-wall regions of the boundary layers. This is in contrast with the original formulation of the model and demostrates the criteria's ability to discern nearwall flow from laminarization that occurs in the core of the cavity due to the stable stratification.

With regards to the cylindrical annuli flow, a schematic of the flow geometry is shown in 6. In this figure, the tangential velocity component at a line located at an angle of 60 • is also shown and compared to quasi-direct numerical simulation data from [START_REF] Addad | Quasi-DNS of natural convection flow in a cylindrical annuli with an optimal polyhedral mesh refinement[END_REF]. It can be seen that RANS and LES results from the dual-mesh framework (solid lines) offer improved prediction when compared to the RANS and LES results computed alone (dashed lines).

Subdomain Wall Function

In addition to the work presented in the previous section, we are developing a modified approach wherein the RANS mesh covers only a narrow near wall region, termed a subdomain wall function (SWF). This method shares much with the previous section and is similar to the wall-modelled LES (WMLES) of [START_REF] Balaras | Twolayer approximate boundary conditions for largeeddy simulations[END_REF] [START_REF] Balaras | Twolayer approximate boundary conditions for largeeddy simulations[END_REF]. As portrayed in 3, the primary domain covers the entire region of the flow field and computes the transport equations of LES. The second domain, termed here the subdomain, overlaps the near-wall area of the coarse LES grid and solves the RANS equations. The two domains are coupled to exchange information between grids. The SWF makes use of ideas of Xiao and Jenny [START_REF] Xiao | A consistent dualmesh framework for hybrid LES/RANS modeling[END_REF] by applying the wall function as a weak volumetric source term in the near-wall LES region beyond the first cell at the wall. This source term corrects the under-resolved near-wall LES fields. In return, the LES grid supplies information to the interface of the RANS grid to complete the boundary conditions of the RANS subdomain. A consistent coupling at the interface is established by computing the partial time average of the instantaneous LES fields. The use of the partial mean differentiates the method from previous WMLES approaches which instead couple instantaneous LES directly to the RANS values.

The partial time-average of LES fields of interest is computed using the exponentially weighed average (EWA), where the averaging time period is sufficiently long to ensure the assumption that the partial timeaveraged filtered LES velocity is equivalent to the Reynolds-averaged velocity ( U i EW A ≈ U i RAN S ). This assumption enables the calculation of the resolved fluctuations about the EWA of the filtered velocity as u i = U i -U i EW A . The LES stress tensor is determined as the summation of the resolved and modelled stresses about the EWA of the filtered velocity τ ij = u i u j + τ r ij . Hence, the LES stress tensor is assumed to be equivalent to the RANS stress tensor as

τ ij EW A ≈ u i u j RAN S .
The momentum transport equation for the LES domain is defined as:

∂U i ∂t + ∂U i U j ∂x j = - ∂P ∂x i + 2ν ∂ ∂x j S ij - ∂τ r ij ∂x j + Q i , ( 1
)
where Q i is the drift term which is defined as:

Q i = (1 -σ L ) U i RAN S -U i EW A γ l1 + G i γ l2 . (2)
The variable σ L is the zoning parameter that determines the areas of the LES grid that are under-resolved. The drift term is only active in the area of the LES grid that overlaps the RANS subdomain. The left term in the square bracket of the drift term relaxes the partial mean filtered LES velocity towards the velocity field of the RANS subdomain. Also, the right term G i rescales the trace of the LES stress tensor towards the RANS turbulent kinetic energy. The RANS subdomain sends information to the LES grid to enable the computation of the drift term. On the other hand, the subdomain receives partial time-averaged LES fields of velocity and pressure at the interface to complete the boundary conditions of the RANS. This work makes use of eddy viscosity models in the subdomain, and so the interface receives the trace of the LES stress tensor (k RAN S = 0.5 τ ij EW A ) and the LES dissipation rate

( LES = 2νS ij S ij -τ SGS ij S ij ).
The performance of the new method is applied to the Re τ = 10, 595 flow through periodic hills as shown in 7, which displays contour plots for both LES and RANS regions. Present results are then compared on 8 with wall-resolved LES of [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] which employed a mesh of 13M cells; noting that the present mesh as less than 0.6M cells, a factor of 20 less. For reference results for LES using the same coarse grid without the SWF are also included. The benefit of the wall function is clear in the near-wall velocity field of the LES towards the RANS, and overall SWF results are in excellent agreement with the reference data. Nowadays, the lattice Boltzmann method (LBM) is a consolidated approach to simulate fluid flows. Its popularity at least in part stems from the intrinsic simplicity of the Bhatnagar-Gross-Krook (BGK) collision operator [START_REF] Bhatnagar | A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems[END_REF], where all the populations are forced to relax towards a discrete equilibrium state derived by applying a Gauss-Hermite quadrature to the continuous Maxwellian distribution. However, the BGK model is well known to be prone to numerical instability in the lowviscosity regime, thus becoming unsuitable for the prediction of turbulent flows. To cope with this problem, a family of collision models based on the relaxation of central moments (CMs) was introduced in 2006 [START_REF] Geier | Cascaded digital lattice Boltzmann automata for high Reynolds number flow[END_REF]. Here, the collision process shows pyramidal hierarchical topology, where the post-collision state of CMs is constructed starting from the lowest order, and then proceeding in ascending sequence, hence the name cascaded lattice Boltzmann method. It has been demonstrated to outperform both BGK and multiple-relaxation-time [29] LBM in terms of stability (the latter based on absolute moments).

More recently, we approached central moments from a different viewpoint. Given a certain lattice, our methodology consists of building a transformation matrix allowing us to move from the space of populations to the one of central moments and vice-versa. The resultant algorithm loses the above-mentioned pyramidal cascaded structure and, as a conse-quence, it can be interpreted as a non-cascaded way to apply the collision step in CMs space [START_REF] Rosis | Nonorthogonal central-momentsbased lattice boltzmann scheme in three dimensions[END_REF][START_REF] De Rosis | Role of higherorder hermite polynomials in the central-momentsbased lattice boltzmann framework[END_REF]. Interestingly, it is also shown that the present method recovers the behaviour of the cascaded LBM while allowing the derivation of forcing terms in an a systematic and straightforward manner. This was thoroughly demonstrated by successfully recovering different sets of governing equations with this approach, hence allowing the simulation of a rich variety of physics problems such as magnetohydrodynamics [START_REF] De Rosis | Advanced lattice boltzmann scheme for highreynolds-number magneto-hydrodynamic flows[END_REF] and multicomponent flows [START_REF] Saito | Color-gradient lattice boltzmann model with nonorthogonal central moments: Hydrodynamic melt-jet breakup simulations[END_REF] among others. As an example, 9 shows an instantaneous snapshot from a turbulent magnetohydrodynamics case.

The group has recently developed a lattice Boltzmann (LB) code primarily for the purpose of testing new approaches. LUMA (LBM at The University of Manchester) [START_REF] Harwood | LUMA: A many-core, FluidStructure Interaction solver based on the Lattice-Boltzmann Method[END_REF] offers grid refinement, turbulence and fluid-structure interaction, and a GPU version is now in development based on our earlier work [START_REF] Mawson | Memory transfer optimization for a lattice boltzmann solver on kepler architecture nvidia gpus[END_REF]. Validation results for turbulent flow are displayed in 9 for the periodic cube array case of [START_REF] Meinders | Vortex structure and heat transfer in turbulent flow over a wallmounted matrix of cubes[END_REF]. Ongoing work is in progress to further develop our LBM capabilities in the context of hybrid RANS LES methods.

A Dual-Solver using LBM

A natural evolution of the dual mesh approach is to use different codes for each region. In this work, focused on developing fast prediction for turbulent flow around urban geometries, the solver in the region of interest is a lattice Boltzmann (LB) solver based on [START_REF] Harwood | LUMA: A many-core, FluidStructure Interaction solver based on the Lattice-Boltzmann Method[END_REF] running on GPUs. Since the LB region is a small volume of the domain the hybrid model overcomes the mesh and memory related drawbacks of using GPU accelerated LB, while maintaining its inherent speed and accuracy. The solver for the remainder of the domain is a finite volume RANS solver running on CPU, which takes ad-vantage of the mesh flexibility and lower memory consumption of the NS algorithm. Moreover, the RANS-LBLES model exploits the availability of heterogeneous CPU/GPU hardware on consumer devices. The method uses two fully-coupled subdomains as shown in 10, wherein the LBLES and the RANS regions exchange required information at overlapping interfaces, marked with numbers in 10. Odd numbers denote the LBLES boundaries that receive data from the RANS solver, while the even numbers denote RANS boundaries that receive data from the LBLES solver. The volume between them is the overlap region, which is part of both subdomains. The external boundaries (boundaries 1, 3 and 5) incorporate the Synthetic Eddy Method (SEM) to generate the instantaneous velocity as needed by the LES. The two subdomains are coupled using the parallel explicit algorithm in the preCICE libraries [START_REF] Bungartz | pre-CICE A fully parallel library for multiphysics surface coupling[END_REF]. In this framework it is noted that eventually unsteady mean flow i.e. URANS could be coupled with LES in the same way.

The SEM inlet boundary has been implemented within a lattice Boltzmann solver for the first time, and tested fo a channel flow at Re τ = 180, where the domain is discretised using a uniform Cartesian grid with y + = 4. 10 displays the streamwise components of the instantaneous velocity u x /u τ .

Conclusions and Ongoing work

This communication has provided a brief overview of recent and ongoing activities in our group at The University of Manchester in the context of Hybrid RANS LES methods; the interested reader is invited to contact us where more details are desired. These activities underpin our primary objective to develop a set of computational tools which enable fast and efficient coupling between different levels of modelling fidelity. The resulting series of coupled Hybrid RANS LES methods are applicable across a range of applications. Furthermore, we have demonstrated the potential for combining lattice Boltzmann method in this framework, particularly given its inherent scalability on GPU processors. Work in our group is also in progress on the more applied side, in particular testing their potential for application to complex configurations which cover a range of lengthscales. Examples include the assessment of wind turbines on environmental landscapes [START_REF] Hewitt | Building a digital wind farm[END_REF], assessment of individual buildings in context of cityscale flows, application to ground vehicles aerodynamics [START_REF] Ashton | Key factors in the use of ddes for the flow around a simplified car[END_REF][START_REF] Ashton | Assessment of rans and des methods for realistic automotive models[END_REF] and cardiovascular flow in and around the heart [START_REF] Al-Azawy | Assessment of turbulence models for pulsatile flow inside a heart pump[END_REF][START_REF] Owen | Computational hemodynamics of abdominal aortic aneurysms: three-dimensional ultrasound versus computed tomography[END_REF]. One ongoing study is motivated by noise levels in road vehicles, wherein we will apply the dual solver approach in section 6 to a generic side mirror (GSM) model [START_REF] Siegert | Numerical Simulation of aeroacustic sound generated by generic bodies placed on a plate: Part II -Prediction of Radiated Sound Pressure[END_REF], consisting of a quarter sphere mounted on a half cylinder. Thus far, work has focused on using Delayed Detached Eddy Simulation (DDES) but will soon incorporate LBM also; initial results are displayed in 11. Other ongoing work will apply these methods to truck 0 in tandem, with an exper-

Figure 1 :

 1 Figure 1: Proof of concept for ELES, Top) Two-region nested LES computation coloured by instantaneous velocity. Btm) Skin friction coefficient normalised by periodic LES value through the entire channel.

1 Figure 2 :

 12 Figure 2: Top) Schematic of TBL domain. Btm) Downstream recovery of skin friction coefficient C f .

Figure 3 :

 3 Figure 3: Velocity and Reynolds stress profiles at Re θ = 5200. Solid : Total stress, Dashed : Resolved stress, Dotted : Modelled stress

Figure 4 :

 4 Figure 4: Q-iso surfaces. Top: View from inlet, Bottom: Side view. Green: Instantaneous, Red and blue: Mean.

Figure 5 :

 5 Figure 5: Square Cavity Case. Distribution of the LES zone weight and Instantaneous Velocity

Figure 6 :

 6 Figure 6: Top) Schematic of the cylindrical annuli and Instantaneous Velocity from simulation. Btm) Tangential velocity along radius at 60 • .

Figure 7 :

 7 Figure 7: Top) Schematic of the subdomain wall function for LES. Btm) Contours of the LES and RANS domains predicted by the subdomain wall function

Figure 8 :Figure 9 :

 89 Figure 8: Profiles of the mean streamwise velocity in the LES domain. Red line: SWF; Green line: Pure LES; Blue line: WMLES; Black markers: reference data

Figure 10 :

 10 Figure 10: Top) Dual solver concept Btm) Streamwise instantaneous velocity on vertical section of a channel flow using SEM with LBM.

Figure 11 :

 11 Figure 11: Ongoing work: Top) Isosurface of Q-criterion for flow around GSM, Btm) Experimental validation case for truck aerodynamics, Embedded simulation for cardiovascular flow.

Acknowledgment

We acknowledge support from the UK Turbulence Consortium (EPSRC Grant EP/G069581/1) and The UK Consortium on Mesoscale Engineering Sciences (Grant No. EP/L00030X/1), both of which provided computational access to ARCHER, UK. The authors also acknowledge assistance from IT Services and the use of the Computational Shared Facility at The University of Manchester.


The author proposes an expression for the equilibrium nearbed concentration C a based on scales of turbulent bursts. This is a significant improvement over existing empirical formulations (VanRijn 1984;Zyserman and Fredsoe 1994). He gives a relation for the bed sediment-entrainment function E in which the outer-scale law is used for the determination of the normalized turbulent bursting period. However, as pointed out by the author, a direct validation of the proposed near-bed sediment entrainment function was not possible, because quantitative data for sediment entrainment under different particle size and hydraulic conditions is lacking. Instead, the author has undertaken a calibration to determine the bursting parameter based on existing data.

ϩ T /A B C
Recently, the discussers have investigated a method to estimate the dynamics of suspended sediment transport capacity from instantaneous mass flux profile measurements by applying a conditional sampling technique (Hurther and Lemmin, unpublished manuscript, 2000). Based on the measured bursting characteristics, the discussers will herein verify and discuss the author's expression of the near-bed equilibrium concentration.

The data set used to discuss these points was obtained by [START_REF] Cellino | Experimental study of suspension flow in open channels[END_REF] in open-channel laboratory experiments using the Acoustic Particle Flux Profiler (Shen and Lemmin 1999a). Detailed information concerning the experimental setup is given in [START_REF] Cellino | Experimental study of suspension flow in open channels[END_REF]. Uniform suspension flows under capacity charge conditions were investigated as shown in Table 5.

In the entrainment flux equation [( 12)], the author uses the bursting period estimated by conditionally sampling the instantaneous velocity field with a half-value shear stress threshold, as defined by Nezu and Nakagawa (1993). The discussers will refer to this sampling condition as ''NN50'' herein.

The author assumes an equilibrium condition [(17)] between the entrainment flux E [( 12)] and the deposition flux D [( 14)]. It is important to note that the entrainment flux E is conditionally sampled while the deposition flux D is not. As will be seen, this difference between the two terms has a number of important implications.

In writing ( 12), ( 14), and ( 17), the author assumes that the particle flux other than the vertical is negligible. Coherent structures are observed to be three-dimensional patterns (Grass et al. 1993) in which nonuniformities can exist. This implies a December 1999, Vol. 125, No. 12, by Z. Cao (Paper 14113). that turbulent mass fluxes other than the vertical one may have to be considered in the diffusion equation.

In order to quantify the suspended sediment transport capacity of the flow structures related to the bursting phenomenon, the discussers investigated experimentally the order of magnitude of the different terms of the conditionally sampled sediment advection-diffusion equation, written as

Ѩ͗c͘

Ѩ͗cu͘ Ѩ͗cw͘

H H H 1 1 1 1 = ϩ (25) ͫ ͬ Ѩz v Ѩx Ѩz 0
where denotes a temporal average over instantaneous val-͗ ͘ H 1 ues, during which the instantaneous shear stress was greater than H 1 times the local mean shear stress value. The ''hole size'' H 1 is defined in Nezu and Nakagawa (1993). Since low sediment concentration situations were investigated, the effect of the suspension on the particle settling velocity was not considered.

The profiles of the first and second term of the right-hand side of (25) as functions of the selection criterion H 1 are presented in Fig. 10. Even for H 1 = 3, which corresponds to the NN50 shear stress level, the longitudinal gradient of the horizontal mass flux is obviously negligible compared with the vertical gradient of the vertical mass flux. The discussers' measurements show that a relation between the conditionally sampled vertical turbulent mass flux and the deposition flux can be considered as a good approximation, thus confirming the author's assumption.

However, in ( 17), the author implies that the total deposition flux at the equilibrium near-bed level is caused by the conditionally sampled entrainment due to bursts that are stronger than NN50. A priori, that hypothesis is not obvious, since there is no physical reason that weaker bursts should not contribute to the deposit. To justify this assumption it would be necessary to prove that the contribution of the weaker bursts is insignificant.

Therefore, the entrainment capacity of those flow structures which are delimited by the threshold level NN50 has to be quantified with respect to the total entrainment flux. In order to investigate this point, we start by calculating the dependence of different relative parameters on the threshold parameter H 1 at depth z/h = 0.08, which is close to the bottom but above the author's equilibrium depth. The parameters R1, R2, and T5 represent the mean shear stress contribution, the mean vertical mass flux, and the time fraction of the unselected flow part, respectively.

Concerning the mean contributions of shear stress and vertical mass flux, the relative vertical mass flux R2 is found to be lower than the relative shear stress contribution R1 (Fig. 11). As discussed in Hurther and Lemmin (unpublished manuscript, 2000), this difference originates from the decorrelation between downward mass fluxes and sweeps moving towards the wall, while ejection events are always highly correlated with ascendant mass flux events over the entire flow depth. This observation indicates that, in the vicinity of the bottom of the flow, the sediment entrainment may arise from pressure fluctuations at the bed, whereas the sediment resuspension process is highly correlated with shear stress ejection events. As indicated in Fig. 11, the NN50 sampling condition corresponds to a value of 3 for H 1 . This is in good agreement with the results given by Nezu and Nakagawa (1993). The data in Fig. 11 demonstrate further that the entrainment at the near-bed level, related to those bursts that obey the NN50 condition as the burst delimiting value, corresponds to only 40% of total the vertical mass flux. It is interesting to note that this 40% entrainment occurs during only 10% of the total time. This is further discussed in Hurther and Lemmin (unpublished manuscript, 2000).

In order to determine the percentage of the total entrainment at the author's near-bed equilibrium depth, the discussers calculated the concentration profiles relative to the near-bed equilibrium concentration for H 1 equal to 3 corresponding to the NN50 sampling condition and for all runs (Fig. 12). The value of C a was obtained from independent suction samples (Shen and Lemmin 1999b). Extrapolating our results to the near-bed equilibrium depth, only about 60% of the total entrainment is found at that level independent of the hydraulic conditions. Individual experiments scatter around this mean value, due to the very strong gradient of the curves in this depth range. This is in contradiction to the author's assumption in (17).

Under the NN50 sampling condition that the author imposes, our measurements discussed above show that the corresponding deposition flux accounts for only 60% of the total deposition flux. As a consequence, based on the results of our measurements, the outer-scale formulation of the volumetric near-bed concentration for low sediment concentration [the author's (19)] should be corrected with a factor of 0.6 as follows:

A C d ( Ϫ )U C 0 c ϱ 0.6C Х (26) a ϩ T w h B 0 c
Obviously, (18) should be corrected in the same way.

VALIDATION OF PROPOSED CORRECTION

As mentioned above, in the discussed paper, the parameter has been determined through calibration. The correction (26) for our data set. First, the normalized bursting period is directly esti-ϩ T B mated from our measurements. For the results summarized in Table 5, a scaling on the outer flow variables was used, applying the author's method based on the NN50 condition. For all runs, we find a bursting period that varies between 3.4 and 3.8. This range is in good agreement with the experimental results previously given by Laufer and Narayanan (1971) and Jackson (1976).

As mentioned by Cao (1997), ( 26) can be written as follows:

A C d U i 0 ϱ 0.6C Х (27) a ϩ T w h B 0 with Ϫ c A = A (28) 1 C c or A = A (29) 2 C c
where A i = mean surface portion of the bursts per unit bed area of those bursts containing enough energy to lift up sediments from the bed (expressed by the NN50 condition); and A c = mean surface portion of all bursts per unit bed area. Cao (1997) proposed ( 28) and ( 29) to relate the last two parameters and found that (28) gives better results [i.e., the combination of ( 16) with (28) results in ( 26)]. If ( 29) is used instead of (28), the following expression is found for the mean near-bed equilibrium concentration:

A C d U C 0 ϱ 0.6C Х (30) a ϩ T w h B 0 c
The validation of the proposed corrections can be achieved by comparing the surface portion of the bursts per unit bed area, A C , estimated with ( 19), ( 26), and (30). Based on the results of Kline et al. (1967) and Kim et al. (1971), Cao (1997) calculated a value of A C roughly equal to 0.02. Fig. 13 shows the values of A C calculated with measured ϩ T B and using ( 19), ( 26), and ( 30) for all investigated runs. It can be seen that results using ( 30), which includes the correction, are in close agreement with the expected value of 0.02, while those obtained using ( 19) give a value of A C that is significantly higher.

SUMMARY AND CONCLUSION

The present discussion is concerned with the formulation of the near-bed equilibrium concentration C a developed by the author considering turbulent bursts scales in the near-wall flow region.

The discussers have used direct instantaneous mass flux measurements in several suspension flow conditions to verify different aspects of the author's analysis. We have shown that the conditionally sampled vertical entrainment flux is by far the dominant contributor to entrainment fluxes, in agreement with the author's assumption. However, when the sampling condition NN50 is applied to the vertical entrainment fluxes, the resulting entrainment flux is not in equilibrium with the total deposition flux as assumed in the author's (17).

Therefore, the discussers propose to introduce a correction factor in his near-bed equilibrium concentration formulation. Finally, the validity of our proposed correction has been checked by calculating the burst surface portion values for one particle size but for different hydraulic conditions. The corrected predictions are found to be in very good agreement with results given in the literature.

An important point in this context is whether the correction the discussers propose is universal. Here we have only worked with one particle size and with different hydraulic conditions. The author's analysis indicates that the parameter is a

ϩ T /A B C
strong function of particle size, particularly for small-sized particles. However, his Fig. 6(a) shows low scatter between the experimental results and the fitted curve based on (22a).

From this observation, one may expect that the value of the correction factor will not vary too much. Obviously our data do not fit onto the curve in Fig. 6(a). Therefore, (22a) will also have to be modified in order to take into account the correction once more data for other particle sizes are available.

Grass, A. J., Stuart, R. J., and Mansour-Tehrani, M. (1993). ''Common vortical structure of turbulent flows over smooth and rough boundaries.'' AIAA J., 31(5), 837-847. Kim, H. T., Kline, S. J., and Reynolds, W. C. (1971). ''The production of turbulence near a smooth wall in a turbulent boundary layer.'' J. Fluid Mech., 50, 133-160. Kline, S. J., Reynolds, W. C., Schraub, F. A., and Runstadler, P. W. (1967). ''The structure of turbulent boundary layers.'' J. Fluid Mech., 30(4), 741-773. Laufer, J., and Narayanan, M. A. B. (1971). ''Mean period of the turbulent production mechanism in a boundary layer. '' Phys. of Fluids, 14, 82-183. Shen, C., and Lemmin, U. (1999a). ''Application of an acoustic particle flux profiler in particle-laden open-channel flow.'' J. Hydr. Res., Delft, The Netherlands, 37(3), 407-419. Shen, C., and Lemmin, U. (1999b). ''Improvements in acoustic sediment concentration profiling using a LMS compensation algorithm.'' IEEE Oceanic Engrg., 23(2), 96-104.

Closure by Z. Cao 4

The writer is grateful to the discussers for their interest in and discussion of the research reported in the paper.

The discussers first raised a question that the turbulent mass fluxes in addition to the vertical component may have to be considered due to the three-dimensional structure of coherent motions. Using their experimental datasets, they then confirm the vertical mass flux is by far the dominant. This conclusion has essentially been well known, as in most previous formulations for sediment transport the streamwise and spanwise fluxes are neglected (for steady uniform flows).

The major concern of the discussers appears to be related to the equilibrium condition, (17), used by this writer. In this connection, it is necessary to emphasize that in equilibrium sediment-laden flows [(17)] must be satisfied, characterizing the balance between the time-averaged entrainment and deposition fluxes. Otherwise, there will be aggradation or degradation of the channel bed, which is not the problem being considered. Therefore, the discussers' concern should be interpreted more appropriately, as whether the dominated role of bursting in entraining sediment can be justified. The discussers show that, at the near-bed level z/h = 0.08 (far above the equilibrium reference level used by this writer), turbulent bursting contributes only 40% of the total entrainment flux and suggest that pressure fluctuations near the bed may contribute to entrainment. Further, the discussers extrapolate the sediment concentration profiles to the writer's reference level a = 10d (equivalent to z/h = 0.0192) and mention that the contribution to entrainment is only 60%. Based on this extrapolation, a correction factor of 0.6 is introduced in (26). Unfortunately, the discusser's extrapolation to the equilibrium reference level z/h = 0.0192 is clearly unreliable because of the highly variable gradient of sediment concentration closer to the bed as well as the apparent scatter of their concentration data. Thus, their correction factor is questionable. At best it is a conjecture, though this writer would confess that weaker bursts should lead to entrainment to a much less extent than energetic bursts. The evidence of fractional contribution of bursting to entrainment at a higher level z/h = 0.08 can only be understood to be the situation locally. Moreover, it is noted that the percentage of contribution of bursting to the total entrainment tends to increase closer to the bed, as observed by the discussers. This trend at least qualitatively supports the premise that turbulent bursting contributes most to entrainment closer to the bed.

4

Prof. and Dir., Inst. of Hydr. Res., Tianjin Univ., Tianjin 300072, China; presently, Visiting Prof., Facu. of Sci. and Engrg., Ritsumeikan Univ., Noji-Higashi 1-1-1, Kusatsu, Shiga 525-8577, Japan. E-mail: zcv98138@se.ritsumei.ac.jp Starting from their extrapolation-based correction factor, the discussers further show that the value of unit bursting area A C can be rendered to be closer to 0.02 than was estimated by this writer based on turbulent bursting of single-phase flows (Cao 1997a). However, it must not be forgotten that the discussers' experiments are for capacity flows. The near-bed concentrations are higher than 0.008 in volume (i.e., 21.31 kg/ m 3 ). In these cases, the interaction between the fluid and the dispersed particle phases is influential [START_REF] Crowe | Numerical models for two-phase turbulent flows[END_REF]. It is recognized that the bursting-particle interaction is as yet poorly known, which is one of the major reasons that this writer has chosen to calibrate the bursting variables ϩ T /A B C using existing measured datasets rather than directly using a value of 0.02 for A C and a value around 5.0 for . The dis-ϩ T B cussers use the value of A C = 0.02 for clear-water flows to back up their extrapolation for capacity flows, which to this writer is unbelievable.

Quantifying bed sediment entrainment continues to be one of the fundamental impediments to refined modeling of sediment transport in turbulent flows. The last four decades have seen enhanced understanding of the physics of sediment entrainment. Turbulent bursting has been experimentally found to play a central role in picking up sediment. This background forms the physical basis of the research, in which bed sediment entrainment has been linked to the time-averaged characteristics of bursting. To the best of this writer's knowledge, this is the first formulation of its kind in the context of fluvial sediment transport. It must, however, be pointed out that the development along this line has had to be based on current understanding of the process and is still in its infancy. High-quality measurement, although extremely difficult close to the bed, is evidently needed. The writer very much appreciates the discussers' viewpoints and would be happy to see their experimental findings relevant to this research.

FIG. 10 .

 10 FIG. 10. Conditionally Sampled Terms of Sediment Diffusion Equation for Run Q50: (a) Vertical Mass Flux Gradient; (b) Horizontal Mass Flux Gradient

2

  Res. Asst., Lab. de Recherches Hydrauliques, Ecole Polytechnique Fe ´de ´rale de Lausanne, CH-1015 Lausanne, Switzerland.

3

  Res. Assoc., Lab. de Recherches Hydrauliques, Ecole Polytechnique Fe ´de ´rale de Lausanne, CH-1015 Lausanne, Switzerland.

FIG. 12 .

 12 FIG. 12. Transport Capacity of Bursts Sampled at Shear Stress Half-Value Threshold Level H 1 = 3

  FIG. 13. Comparison of Surface Portion per Unit Area of Nursts versus Shields Parameter at Equilibrium Depth, Calculated with Cao's [Eq. (19)] Model and Herein Corrected Models [Eqs. (26) and (30)] nificant discrepancy between the results from his model based on burst scales and the results from other models. However, in order to determine whether the proposed correction factor provides a better representation of the physical processes, we can analyze the results of the author's calibration. Since the values of obtained by (22a) are based ϩ T /A B C on the calculations using (19), we can compare them directly to the parameter determined with the herein proposed ϩ T /A B C

TABLE 5 .

 5 Hydraulic Parameters for Experiments

	Run	Q (m 3 /s)	h (cm)	U (cm/s)	u * (cm/s)	S (ϫ10 Ϫ3 )	R h (ϫ10 3 )	F h	ϩ k S	d (mm)

S (kg/m 3 ) 0 (mm/s) c (ϫ10 Ϫ2 ) (ϫ10 Ϫ1 )

C a (kg/m 3 ) ϩ

ACKNOWLEDGMENT

The present discussion is based on work sponsored by the Swiss National Science Foundation, Grant No. 2100 050739. The discussers are grateful for the support.


Iman Kamehkhosh
email: iman.kamehkhosh@tu-dortmund.de

Dietmar Jannach
email: dietmar.jannach@aau.at

Geoffray Bonnin
email: bonnin@loria.fr

How Automated Recommendations Affect the Playlist Creation Behavior of Users

Keywords: Information Retrieval : Recommender Systems, Human Computer Interaction (HCI) : User studies Music Recommendation, Playlist Creation, User Study

Modern music platforms like Spotify support users to create new playlists through interactive tools. Given an empty or initial playlist, these tools often recommend additional songs, which could be included in the playlist based, e.g., on the title of the playlist or the set of tracks that are already in the playlist.

In this work, we analyze in which ways the recommendations of such playlist construction support tools influence the behavior of users and the characteristics of the resulting playlists. We report the results of a between-subjects user study involving 123 subjects. Our analysis shows that users provided with recommendation support were more engaged and explored more alternatives than the control group. Presumably influenced by the recommender, they also picked significantly less popular items, which leads to a higher potential for discovery. The effort required to browse the additional alternatives, however, increased the users' perceived difficulty of the process.

INTRODUCTION

Creating and sharing playlists is a common feature of most of today's music platforms. The manual construction of playlists by users can however be a comparably complicated and time consuming task [START_REF] Jo | More of an Art than a Science': Supporting the Creation of Playlists and Mixes[END_REF]. One way to help playlist creators is to provide them with automated suggestions for additional tracks while they are creating a playlist. Such a functionality can be found on some of today's music platforms, like Spotify and Pandora. In the literature, a number of algorithms were proposed in the past to determine a set of suitable tracks given, e.g., a partial playlist [START_REF] Bonnin | Automated Generation of Music Playlists: Survey and Experiments[END_REF]. The evaluation of such algorithms is mostly based on offline experiments with a focus on prediction accuracy, i.e., on the algorithms' capability of predicting the next tracks that were picked by users.

©2018. Copyright for the individual papers remains with the authors.

Copying permitted for private and academic purposes.

MILC '18, March 11, 2018, Tokyo, Japan While accurate track relevance predictions are important, such experiments cannot inform us about whether users will actually adopt the recommendation functionality and how the recommendations influence their behavior. With this work, we aim to obtain a better understanding of these aspects which, to our knowledge, have not been studied in that form in the literature before. We conducted a between-subjects user study involving 123 subjects, where the participants' task was to create playlists for a given topic using a web application that was developed for the study. One group of participants was provided with additional recommendation functionality, whereas the control group could only rely on the provided search functionality. A post-task questionnaire was used to assess additional aspects like the perceived difficulty of the task.

Our analyses revealed that the recommendation support was well accepted by the participants and that recommendations therefore represent a valuable tool for users. Almost half of the users that were provided with recommendations picked at least one of the recommended tracks, which is an unusually high proportion for the domain of recommender systems. Furthermore, we could observe that the selection of tracks was seemingly influenced by the recommendations even if no track was actually included in the playlists, i.e., the recommendations served as inspirations for the participants. Participants with recommendation support also significantly explored more tracks within the same period of time, and they more often chose less popular tracks based on the suggestions by the recommender. As users are more engaged and explore more options from the long tail, the recommender therefore measurably increases the potential of discovering new tracks or artists. Also, this increased user engagement can help providers to collect more information about the users' preferences.

Finally, the post-task questionnaire, to some surprise, revealed that the participants found the playlist creation task to be more difficult when recommendations were provided. This was the case even when they actually did not pick any of the tracks. This observation suggests that the user interface (UI) has to be carefully designed, since parts of this might be caused by the increased complexity of the web application that included a recommendation component.

PREVIOUS USER STUDIES

Compared to the number of papers that report results of offline experiments, the number of user studies is limited. Many of these previous user studies in the music domain focus on how users search for music and on the social or contextual aspects of listening to music [START_REF] Jo | More of an Art than a Science': Supporting the Creation of Playlists and Mixes[END_REF][START_REF] Sally | Finding New Music: A Diary Study of Everyday Encounters with Novel Songs[END_REF][START_REF] Downie | Music Information Retrieval[END_REF][START_REF] Lamont | Short-and Long-term Musical Preferences: What Makes a Favourite Piece of Music?[END_REF][START_REF] Ha | Survey Of Music Information Needs, Uses, And Seeking Behaviours: Preliminary Findings[END_REF]. In [START_REF] Jo | More of an Art than a Science': Supporting the Creation of Playlists and Mixes[END_REF], for example, interviews and web forum posts related to the construction of playlists were analyzed. The authors found that playlists (mixes) are often created with a theme or topic in mind, e.g., a genre, an event or a mood. In our study, we therefore also ask participants to create a playlist for one of several given topics. According to the suggestion in [START_REF] Jo | More of an Art than a Science': Supporting the Creation of Playlists and Mixes[END_REF], our online application does not automatically include additional tracks, but presents recommendations as a side information.

Questions of the UI design were also in the focus in [START_REF] Baur | Rush: Repeated Recommendations on Mobile Devices[END_REF] and [START_REF] Swearingen | Interaction Design for Recommender Systems[END_REF]. For instance, [START_REF] Baur | Rush: Repeated Recommendations on Mobile Devices[END_REF] proposed an interactive track recommendation service called rush, optimized for touchscreen devices and, among other aspects, analyzed its usability for lefthanded and right-handed users. In the application we used for our study, recommendations were positioned as a horizontal list on the bottom of the screen, which is common also for e-commerce sites like Amazon.com.

A few studies explore the users' interactions with music recommender systems [START_REF] Andjelkovic | Moodplay: Interactive Mood-based Music Discovery and Recommendation[END_REF][START_REF] Kamalzadeh | TagFlip: Active Mobile Music Discovery with Social Tags[END_REF] and the quality perception of music recommendations [START_REF] Barrington | Smarter than Genius? Human Evaluation of Music Recommender Systems[END_REF][START_REF] Kamehkhosh | User Perception of Next-Track Music Recommendations[END_REF]. The recent work of [START_REF] Kamehkhosh | User Perception of Next-Track Music Recommendations[END_REF] provided evidence that users prefer recommendations that are coherent with the recently played tracks in different dimensions. Their results also indicated that the participants tend to evaluate recommendations better when they know the track or the artist. In contrast to [START_REF] Kamehkhosh | User Perception of Next-Track Music Recommendations[END_REF], we do not compare different recommendation algorithms but the effects of the existence of a recommender.

Finally, questions related to the factors that influence which tracks users choose for inclusion in a playlist in different situations were analyzed in [START_REF] Kamalzadeh | A Survey on Music Listening and Management Behaviours[END_REF][START_REF] Stumpf | When Users Generate Music Playlists: When Words Leave Off, Music Begins?[END_REF] and [START_REF] Tintarev | Sequences of Diverse Song Recommendations: An Exploratory Study in a Commercial System[END_REF]. The study in [START_REF] Kamalzadeh | A Survey on Music Listening and Management Behaviours[END_REF], for example, showed that mood, genre, and artists are the most important factors for users when selecting the tracks, which is in line with the outcomes of the study of [START_REF] Jo | More of an Art than a Science': Supporting the Creation of Playlists and Mixes[END_REF]. Similar to their work, we explicitly asked users about their decision factors after the task and report the results in Section 5.

STUDY DESIGN

A general limitation of laboratory studies is that users might behave differently in a "simulated" situation than when they normally listen to music, e.g., at home. To alleviate this problem, we provided an online application to enable users to participate in the study when and where they wanted to.

All participants were asked to create a playlist -using the developed application -for one of the following pre-defined and randomly ordered themes: rock night, road trip, chill out, dance party, hip-hop club. 1 The participants were randomly assigned to one of two groups. One group (called Rec) received additional recommendations, as shown at the bottom of Figure 1. The control group (NoRec) was shown the same interface but without the recommendation bar at the bottom.

All participants could use the provided search functionality. When a user of the Rec group added an item to the playlist, the 20 provided recommendations were immediately updated. 2When the playlist contained at least six tracks, the participants could proceed to the post-task questionnaire.

The first part of the questionnaire was presented to all participants and contained a list of quality factors for playlists mentioned in the literature, which were either related to individual tracks (e.g., popularity or freshness) or to the list as a whole (e.g., artist homogeneity). The participants were asked to rank these quality factors or mark them as irrelevant.

In the next step, participants of the Rec group were asked if they had looked at the recommendations during the task and if so, how they assessed their quality in the following dimensions: relevance, novelty, accuracy, diversity (in terms of genre and artist), familiarity, popularity, and freshness. Participants could express their agreement with our provided statements, e.g., "The recommendations were novel", on a 7-point Likert scale item or state that they could not tell.

In the final step, all participants were asked (a) how often they create playlists, (b) about their musical expertise, and (c) how difficult they found the playlist creation task, again using 7-point Likert scale items. Free text form fields were provided for users to specify which part of the process they considered the most difficult one and for general comments and feedback. Finally, the participants could specify their age group.

SPOTIFY'S RECOMMENDATION ALGORITHMS

Both the search and the recommendation functionality in the study were implemented using the public API of Spotify. 3This allowed us to rely on industry-strength search and recommendation technology. As of 2017, Spotify is a leader across most streaming-service markets 4 and the recommendations produced by the service result from several years of A/B testing [START_REF] Sarrafi | A/B testing at Spotify[END_REF].

Companies do not usually reveal the details of the algorithms that they use. According to the documentation of the Web API of Spotify, recommendations are aimed to create a playliststyle listening experience based on seed artists, tracks and genres. In this context, the available information for a given seed item is used to find similar artists and tracks. Some additional information can be obtained from what the company presents in public presentations and industry reports. Public presentations around the year 2014, such as [START_REF] Johnson | Algorithmic Music Discovery at Spotify[END_REF] and [START_REF] Johnson | From Idea to Execution: Spotify's Discover Weekly[END_REF] indicate that Spotify relied at that time mainly on collaborative filtering techniques for generating recommendations. The latest presentation was made right after Spotify acquired The Echo Nest, a music intelligence platform that focused on the analysis of audio content, and Spotify announced they were going to also utilize content-based techniques. In more recent presentations, such as [START_REF] Steck | Interactive Recommender Systems with Netflix and Spotify[END_REF], the authors report that Spotify uses an ensemble of different techniques including NLP models and Recurrent Neural Networks as well as explicit feedback signals (e.g., thumbs-up / thumbs-down), and also audio features for certain recommendation tasks, but it is unclear from the presentations which techniques are used for which types of recommendations (radios, weekly recommendations, playlists, etc.).

RESULTS AND OBSERVATIONS

General Statistics

Participants. Most of the 123 participants who completed the study were students of universities in Germany and Brazil; a smaller part was via invitations on social network sites. Most (84%) of the participants were aged between 20 and 40. On a scale between 1 and 7, the median of the selfreported experience with music was 5, i.e., the majority of the participants considered themselves experienced or interested in music. Most of the participants, however, do not create playlists regularly; about 25% answered the question with a 5 or greater (median=3). The median of the perceived difficulty of the playlist creation task was 4, i.e., the majority of the subjects found the task comparably difficult. 5Topics and playlist length. Rock night and road trip were the most often selected themes. Each of them was selected in about 30% of all trials. Chill out (20%), dance party (15%) and hip-hop club (5%) were less frequently chosen. The average time for the participants to create one playlist (with at least 6 tracks) was 7.29 minutes and a created playlists contained, on average, 8.44 tracks.

Recommendation use. 57% of the participants were assigned to the Rec group (with recommendations). Almost half of these participants (49%) drag-and-dropped at least one of the recommended tracks to their playlists. We denote this group as RecUsed. The other half will be denoted as RecNotUsed.

Study Outcomes Impact of Recommendations on Users and their Behavior

Adoption of recommendations. When users actively used the recommendations, they relied on them to a significant extent. At the end, about 38% of the tracks of the playlists that were created by users of the RecUsed group were taken from the recommendations (mean=3.2 recommendations per playlist). This is a strong indicator of the general usefulness of a recommendation component in that domain, considering that in general e-commerce settings sometimes only every 100th click of a user is on a recommendation list and recommendations are used only in about 8% of the shopping sessions [START_REF] Jannach | Session-based Item Recommendation in E-Commerce: On Short-Term Intents, Reminders, Trends, and Discounts. User-Modeling and User-Adapted Interaction[END_REF].

Increased exploration. The participants that received recommendations played significantly6 more tracks when creating their playlist than the participants of the NoRec group (mean value 14.4 and 9.8, respectively). This value is even higher for the RecUsed subgroup, i.e., those who actively used the recommendations, with an average of 20.3. However, the participants of the Rec group, on average, only needed 30 seconds more to create the playlists 7 , i.e., the recommender helped them explore, and potentially discover, many more options in about the same time.

Such a high exploration is interesting as it not only increases the chances of music discovery but also allows the service provider to gather more information about users' preferences. This can in turn lead to providing better recommendations and a better listening experience for users, as well as higher customer retention and business value for the provider [START_REF] Jannach | Recommendations with a Purpose[END_REF].

Difficulty of playlist creation. To some surprise, the recommendation component did not make the playlist creation task easier for users but slightly added to the perceived complexity, with a median of 4 (higher complexity) for the Rec group and 3 for the NoRec group. This slight but not statistically significant difference could be caused by the more complex recommendation UI, as well as by the fact that users explored more options in the Rec condition as mentioned above (see, for example, [START_REF] Bollen | Understanding Choice Overload in Recommender Systems[END_REF] for a discussion on "choice overload" in recommender systems).

Impact of Recommendations on the Resulting Playlists

To analyze the impact of the provided recommendations on the resulting playlists, we queried the musical features of the tracks contained in the created playlists through Spotify's API.

Table 1 shows a list of these features. The average values and standard deviations of the features for each of the study groups are shown in Table 2. Several differences can be observed. We limit the discussion here to the most pronounced and statistically significant effects, which can be found with respect to the popularity and the freshness of the tracks that were used by the participants.

Popularity effect -promoting less popular items. 8 Users of the recommendation service added significantly less popular tracks to their playlists. This is in line with the observations from [START_REF] Jannach | Biases in Automated Music Playlist Generation: A Comparison of Next-Track Recommending Techniques[END_REF] where the recommendations by a commercial service were less popular (in terms of play-counts) than the tracks that users selected manually.

Recency effect -promoting newer tracks. Using recommendations also slightly but statistically significantly increased the freshness (release year) of the selected tracks. About 50% of the tracks of the playlists created by the participants who used the recommendations (RecUsed) were released in the last five years. This value is 40% for the RecNotUsed group and 34% for the NoRec group.

"Mere-Presence" Effect of Recommendations

When we compare the average musical features of the tracks recommended to the Rec group and those that were manually selected by the control group (NoRec), we can observe several significant differences regarding, e.g., danceablity, energy, popularity, freshness (release year), speechiness, or tempo i.e., the recommender often picks quite different tracks than users would do, see [START_REF] Jannach | Biases in Automated Music Playlist Generation: A Comparison of Next-Track Recommending Techniques[END_REF]. On the other hand, when comparing the recommendations to what the participants in the RecUsed group selected, the differences are no longer statistically significant, except for the popularity aspect. This is not surprising, as these participants often accepted the recommendations. The somewhat surprising aspect, however, is that the differences between the recommendations and the tracks of the RecNotUsed participants are also no longer significant. It is thus possible that the subjects in this group were biased (or inspired) by the presence of the recommendations. One indicator in favor of that possibility is that an overlap of 34% could be observed in terms of the artists that appeared in the recommendations and that were selected manually by the subjects. This means that the recommendations presumably influenced what users selected. This effect was previously investigated in a user study [START_REF] Köcher | Investigating Mere-Presence Effects of Recommendations on the Consumer Choice Process[END_REF], where the participants exhibited a tendency to select items that were content-wise similar to a (random) recommendation.

A possible explanation why participants still selected tracks that are on average more popular than the recommended ones might lie in the intended use or context of the playlist to be created. If participants, for example, assumed that the playlist is designed to be played for a group of listeners (e.g., at a dance party), they might prefer to pick tracks that are presumably known by many people. In fact, the tracks selected for "dance" playlists were significantly more popular than those for "chill out" playlists.

Investigating Quality Criteria for Playlists

In order to better understand which quality characteristics one should consider when designing algorithms for playlist construction support, we analyzed the rankings that were provided by the participants in the post-task questionnaire, see the section on the study design. To determine the overall ranking, we used the Modified Borda Count method [START_REF] Emerson | The Original Borda Count and Partial Voting[END_REF], which can be applied when some rankings are only partial, i.e., when not all items are ranked. The results are shown in Table 3.

The results indicate that, overall, the participants consider the homogeneity of musical features, e.g., tempo, energy or loudness along with the artist diversity of the resulting playlist as the most relevant quality criteria for playlists. On the other end of the spectrum, the order of the tracks in a playlist and their freshness were the least relevant aspects for the participants.

When looking at the different study groups, some smaller and not statistically significant differences in the rankings can be observed. Participants who used the recommendations considered track transitions less relevant than participants of the other groups. One potential explanation could be that the recommendations by the system (and, likewise, the created playlists) were perceived by users to be comparably coherent, e.g., in terms of the tempo, and the participants of the Rec group therefore paid less attention to the transitions. Another explanation is that the RecNotUsed participants are in general more demanding, and did not use some of the recommendations because they did not allow to make satisfying transitions.

Other differences with respect to the quality criteria rankings were observed for the "lyrics" aspect, which was considered more important for road-trip and hip-hop playlists. Finally, popularity was considered a very important criterion almost exclusively for dance playlists for the above-discussed reasons.

General Quality Perception of the Recommendations About 75% of the participants in the Rec group stated that they have looked at the recommendations when constructing the playlists. In the post-task questionnaire, these users were asked about their quality perception regarding the recommendations (as provided by Spotify's service). Different quality dimensions were assessed and users could express their agreement with related statements using 7-point Likert scale items.

In our analysis, we considered answers that were greater than 4 as positive responses.

Specifically, we asked users to what extent the recommendations (1) matched the given topic; (2) helped them discover novel tracks or artists; (3) matched their general interests; (4) were already known to them; [START_REF] Bonnin | Automated Generation of Music Playlists: Survey and Experiments[END_REF] were diverse in terms of genre and artist; (6) were generally popular or mainstream; and [START_REF] Sally | Finding New Music: A Diary Study of Everyday Encounters with Novel Songs[END_REF] were from trending music.

The results show that 62% of the respondents perceived the recommendation as topic-related (e.g., they match the playlist topic). More than half of the respondents also found the recommendations to match their interests and diverse in terms of genre and artist. With respect to novelty and familiarity, the results were mixed. In both cases many participants provided feedback with values above 4 (e.g., 45% in case of the novelty); a substantial amount of the participants however also found the recommendations of very limited novelty and familiarity.

We further compared the quality perception of the recommendations for the two groups of RecUsed and RecNotUsed. The only statistically significant difference we observed was in terms of novelty. In other words, those participants who used the recommendations in their playlists believed that the recommendations were significantly more novel (median=5) than the other participants who looked at the recommendations but did not accept them (median=3). One plausible explanation of this result could be that the participants who did not use the recommendations had a higher music expertise and were more demanding. However, this explanation is not in line with what the participants claimed in their answers, as the RecUsed group had a generally higher music expertise score than the RecNotUsed group (although the difference was not statistically significant). Another explanation is that participants were generally willing to discover novel music and the recommendations were mainly adopted by those participants who found the recommendations to be so.

Overall, the results indicate that Spotify's algorithms were quite successful in determining suitable tracks for recommendation. At the same time, the mixed results for some other factors are not necessarily negative either, because the impor-tance of factors like the appropriate popularity level of the tracks -according to our analysis -can depend on the topic and intended purpose of the playlist.

Qualitative Analysis

We finally asked the participants to specify what they found to be the most difficult aspect when creating a playlist. To remember the right tracks was most often mentioned as a difficult task by the participants (e.g., "Keeping up with the options. I often can't remember all the songs and performers I know").

To find good music for everyone was another difficult task indicated by many (e.g., "To meet the taste of all, as playlists are often created for occasions in which several people listen to them"). Some participants also emphasized on the smoothness of the track transitions (e.g., "The most difficult thing is to think of tracks that are connected on several dimensions at the same time such as their energy, the transitions, their level of sophistication, etc.").

Recommendation tools have the potential to help users in many of these dimensions, e.g., remind them of tracks or artists they know and liked in the past, to find tracks that are popular in a certain community, or to find tracks that are coherent with the current playlist with respect to musical features. In fact, several participants that had no recommendation support have written about the potential advantages of having such a system: "Sometimes I do not know what to search for and I expect the system suggest something. I prefer to choose just some songs and then listen/receive recommendations related to my initial choices".

CONCLUSIONS

With our study we aim to shed more light on the perception and adoption of recommender systems that are designed to support users in the playlist construction process. The results of the study not only showed a high general adoption ratealmost half of the participants picked at least one track from the provided recommendations, while the proportion in the domain of recommender systems is often lower than 10% -but also that the presentation of recommendations can indirectly influence the choices of the users. Our study, therefore, provides additional evidence that recommenders can be a valuable tool to steer consumer behavior.

Generally, participants who received recommendations also explored more alternatives during the process (without needing much more time), which can be interpreted as a higher engagement in the task. As a result, a recommender system has the potential to increase the involvement of the user with the service, which from a business perspective in the best case leads to higher customer retention rates.

Finally, the study also showed that including a recommendation component can lead to a higher perceived difficulty of the task for users, which emphasizes the importance of focusing on intuitive UI designs.

Figure 1 .

 1 Figure 1. Web application used in the study.

Table 1 .

 1 Description of the collected information for the tracks, as provided by Spotify: https://developer.spotify.com/web-api/ get-audio-features/.

	Information	Description
	Acousticness	Absence of electrical modifications in a track.
	Danceability	Suitability of a track for dancing, based on various
		information including the beat strength, tempo, and
		the stability of the rhythm.
	Energy	Intensity released throughout a track, based on
		various information including the loudness and
		segment durations.
	Instrumentalness Absence of vocal content in a track.
	Liveness	Presence of an audience in the recording.
	Loudness	Overall loudness of a track in decibels (dB).
	Popularity	Popularity of a track, based on the its total number of
		plays and the recency of those plays.
	Release year	Year of release of a track.
	Speechiness	Presence of spoken words in a track.
	Tempo	Speed of a track estimated in beats per minute
		(BPM).
	Valence	Musical positiveness conveyed by a track.

Table 2 .

 2 Average (Avg) and standard deviation (Std) of the musical features of the resulting playlists in different groups. * indicates statistical significance in comparison with the RecUsed group.

		RecUsed	RecNotUsed	NoRec
		(34 participants)	(36 participants)	(53 participants)
	Feature	Avg	Std	Avg	Std	Avg	Std
	Acousticness	0.22	0.28	0.17*	0.23	0.17*	0.26
	Danceability	0.56	0.18	0.59*	0.17	0.54	0.17
	Energy	0.68	0.24	0.70	0.19	0.73*	0.23
	Instrumentalness	0.16	0.32	0.12	0.28	0.12	0.27
	Liveness	0.20	0.17	0.21	0.18	0.21	0.17
	Loudness (dB)	-7.68	4.59	-7.60	3.67	-7.52	4.72
	Popularity	50.7	21.9	55.7*	17.1	54.3*	21.3
	Release year	2005	12.47 2002* 15.73 2003* 13.08
	Speechiness	0.07	0.07	0.08	0.08	0.08	0.08
	Tempo (BPM)	123.0	28.7	122.5	27.9	122.9	28.6
	Valence	0.50	0.26	0.53	0.25	0.49	0.24

Table 3 .

 3 Modified Borda Count: Ranking of Playlist Quality Criteria.

	Criteria	All	RecUsed RecNotUsed NoRec
	Homogeneity of musical	250	68	79	103
	features, e.g., tempo
	Artist diversity	195	55	62	78
	Transition	122	30	46	46
	Popularity	106	39	34	33
	Lyrics	95	32	34	29
	Order	74	12	33	29
	Freshness	32	12	11	9

Note that such themes actually also convey an intended use or purpose. The best example is the "road trip" theme, whose corresponding playlists are supposed to be played when driving a car for a prolonged amount of time.

Recommendations were displayed after the first track was added.

https://developer.spotify.com/web-api/

http://www.midiaresearch.com/downloads/ midia-streaming-music-metrics-bundle/

The collected data is ordinal, i.e., a ranking of the response levels is possible. However, we cannot assume equidistance between the response levels, and reporting mean and standard deviation values is often considered questionable in the literature.

To test for statistical significance for the ordinal data, we use the Mann-Whitney U test and for the interval data we use the Student's t-test, both with p < 0.05.

The differences were not statistically significant.

The popularity of the tracks were also determined using the Spotify's API, with values between 0 and 100 (lowest to highest popularity).

ACKNOWLEDGEMENTS

We would like to thank Adrian Skirzynski who contributed to this work with the results of his master's thesis.