Isabelle Mme

Calmet

M Gilles Bergametti

Je Remercie

Directeur M Philippe Blanc-Benon

M Michel Lmfa

Lance D'

Jérôme, Sébastien Jean-Marc Gilles

Srikanth Merci

Philippe Bernard

Aurélien, Emmanuel, Gherardo, Johannes, Zlatko, Quentin Antoine Que J'oublie

Dan-Gabriel Calugaru

Huang Gang

G Huang

S Simoëns

I Vinkovic

C Le Ribault

S Dupont

G Bergametti

Law-of-the-wall in a boundary-layer over regularly distributed roughness elements

Keywords: atmospheric boundary-layer, Gaussian hill, roughness, trapping of solid particles, recirculation zone variables Wall-bounded flow; roughness sublayer;

remarques pertinentes ayant contribué à améliorer les analyses présentées dans

Résumé

Afin de mieux comprendre les mécanismes liés à l'érosion du sol sous l'effet du vent, le transport de particules solides dans un écoulement de couche limite turbulente à l'échelle d'une soufflerie est étudié à l'aide de simulations numériques. La présence d'une ou plusieurs collines Gaussiennes au sol permet d'étudier les effets de la topographie sur le transport, le dépôt et la réémission de particules solides. L'écoulement du fluide porteur est résolu par la Simulation des Grandes Échelles (SGE). Des modèles de paroi pour la vitesse du fluide sont implémentés afin de mieux représenter l'écoulement proche d'une colline. Le mouvement des particules est pris en compte par un suivi Lagrangien. Des modèles d'envol et de rebond sont développés et utilisés pour prendre en compte l'émission et l'impact au sol des particules.

Dans la première partie, l'écoulement au dessus de collines transversales est simulé et validé par des comparaisons avec différentes expériences. Selon Oke [1988], l'écoulement dans la canopée urbaine peut être schématiquement caractérisé par différents régimes en fonction du positionnement relatif des obstacles. Ce concept est appliqué au cas des dunes, assimilées à des collines dans notre étude. L'accent est mis sur la zone de recirculation (ZR) formée derrière ces collines. Les variations de la ZR sont examinées en fonction de différents paramètres dont la configuration des collines et le nombre de Reynolds. De plus, une étude portant sur la souscouche rugueuse est effectuée de façon à déterminer l'effet de la rugosité due à la couche de particules solides au sol.

La seconde partie du travail porte sur la simulation des particules au dessus des collines. L'objectif est l'amélioration des modélisations concernant l'envol, le rebond et le couplage entre le fluide et les particules. Un premier travail de validation est réalisé en utilisant le modèle complet de transport des particules solides. En particulier, l'évolution du flux d'émission des particules, estimé par le modèle d'envol, en fonction du nombre de Shields, donne des résultats comparables aux modèles classiques de saltation et aux expériences de la littérature. Au-dessus des collines, le transport des particules solides est étudié par des profils de concentration et de vitesse moyenne. Pour analyser les résultats, deux cartographies sont réalisées. La première donne l'intensité des événements locaux et instantanés qui seraient à l'origine de l'évacuation des particules piégées au sein de la ZR. La seconde montre la distribution des particules déposées au sol. Ces résultats permettent d'identifier des zones sujettes à l'érosion et à l'accumulation autour des collines. Enfin, les flux des particules piégées et déposées à l'intérieur de la ZR sont quantifiés et comparés aux flux des particules émises en amont. Ces flux, bien que faibles par rapport au flux entrant, contribueraient aux migrations des dunes et à l'avancée des déserts.

Mots clés : couche-limite atmosphérique, colline Gaussienne, rugosité, piégeage de particules solides, zone de recirculation

Table of Contents

List of Figures xiii

Nomenclature xix 2.2 A computational grid cell depicting the staggering of the coordinates and dependent variables. Figure taken from the ARPS User's Guide [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF].

2.3 Illustration of a computational domain generated using coordinate transformation (Equation 2.30) and vertical grid stretching (Equation 2.31). Note that only every second grid line is shown for clarity. . . . . . . . . . . . . . . . . . . . .

2.4

Illustration of the spatial domain decomposition in ARPS from [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF].

xiv List of Figures 2.5 Rigid wall condition on a staggered grid. Two layers of the near-wall cells and one layer of boundary cells (under the physical wall) are pictured. . . . . . . .

2.6

Errors in the determination of the near-wall velocity gradient. Blue line: true velocity gradient. Red dashed line: approximation by Equation 2.41. . . . . . .

2.7

Comparison between the base state profile (red line, Equation 2.51 in Section 2.5.1)), the law-of-the-wall of von [START_REF] Von Kármán | The analogy between fluid friction and heat transfer[END_REF] (black, dotted line) and of Spalding [1961] (blue, dashed line), and the DNS results of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] (circles). . . . . Solid and dotted lines: DNS results of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] with Re D 660, for u 0 C and w 0 C , respectively. Red, LES. Blue, experiments. . . . . . . . 3.6 Particle behaviour based on the take-off curve of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] using the Stokes number and the crossing-trajectories parameter.

Introduction

In this chapter, an overview of the problem of wind erosion and sand transport is presented.

Then, the problem is formulated and decomposed into three key issues, using recent experimental and numerical results from literature. Finally, the approach and the objectives of the current study are addressed. The structure of the manuscript is given at the end.

Context

The transport of solid particles by wind is a typical issue in the study of aeolian processes.

It is a dynamic process of complex and multi-scale nature. Whereas aerosols interact with turbulence structures at a local scale, large-scale aeolian transport events widely present in nature, such as sand storms or pollen dispersion, take place across regional or even continental distances. Over the last decade, studies on these atmospheric phenomena often involve the use of remote sensing by satellites in order to obtain large amount of data over large scales. However, the physics behind these natural phenomena are not yet fully understood. Notably, two difficulties are involved in this subject: the turbulent nature of the large-scale atmospheric processes, and the behaviour of wind-driven particles, which act and interact at a smaller scale at which fundamental questions on turbulence remain unanswered.

Atmospheric boundary-layer

Planet Earth is covered by a layer of air, the atmosphere, that separates us from the outer space. Compared with the mean radius of Earth, 6370 kilometres, the atmosphere is very thin: 50 % of the atmosphere mass is within 5.5 kilometres and 99.9 % is below 49 kilometres above the sea level. Despite its apparent two-dimensional character, flow motions inside the atmosphere are largely 3D and present a multitude of scales, ranging from a few millimetres to continental scales. The relevant time scales thus range from a fraction of a second to several months or even years. In the meteorological domain, a three-level classification is often used to categorise these scales into micro-, meso-, and macro-scales. More loosely, terms such as local, regional, and global are also used for the same objective.

Notwithstanding regional and diurnal variations, a 1D approximation is often adopted to decompose the vertical structure of the atmosphere into four layers: the troposphere, the stratosphere, the mesosphere and the thermosphere, with increasing distance over the ground. Among them, the troposphere is the layer where most of the weather phenomena on Earth occur. At its lower end, about 1 km above the sea level, is the atmospheric boundary layer (ABL). In mesoscale studies, the scope is often limited to the so-called atmospheric surface layer, which corresponds to the lowest part of the ABL in which the effect of the planet surface is more relevant compared to the Coriolis force due to Earth's rotation. Variations of physical quantities such as wind speed, temperature and aerosol concentration with height are the sharpest inside this surface layer. These high gradients often induce high levels of momentum, energy and mass exchanges across the ABL. Moreover, owing to its proximity to the ground, this layer directly interacts with Earth's ecosystem and is thus of great relevance to human beings. This layer has often been studied from the surface by means of micrometeorological masts and towers.

The understanding of the surface layer is important because it is ultimately wind that powers the processes of solid particle transport by the atmosphere. The wind-driven entrainment, transport, and deposition of sand and mineral dust as well as their effects are generally referred to as "aeolian processes" (after the Greek god Aeolus, the keeper of the winds). These processes basically involve the interaction of the atmosphere with the lithosphere, or the solid surfaces of the planets. In the next section, we will present one process of aeolian nature that is of particular importance to human population: wind erosion.

Wind erosion and sand dunes

Wind erosion is a process of wind-forced movement of soil particles. In principle the process includes a whole sequence of particle emission, transport and deposition. From a global perspective, wind erosion is a geological and climatic phenomenon which takes place over long periods of time across the globe, albeit more frequent in arid areas. As a consequence of this process, large quantities of minerals and organic matter are carried with dust and redistributed around the world. This action is often beneficial, especially in the creation of the so-called "loess", which are wind-driven sediments brought from faraway sources and has contributed to the richness of farmlands over the world. It was recently found that 40 million tons of mineral dust is swept across the Atlantic ocean from the Sahara to the Amazon each year and greatly fertilizes the soil of the rain forests [START_REF] Koren | The Bodélé depression: a single spot in the Sahara that provides most of the mineral dust to the Amazon forest[END_REF]. Desertification is a term coined for the conversion of arable land to deserts. At present, deserts cover about a fourth of the world's land area. Besides, thousands of square kilometres are converted to deserts annually [START_REF] Greeley | Wind as a geological process: on Earth, Mars, Venus and Titan[END_REF]. Although there is still controversy in the role played by humans in the causes of deserts, it is of no doubt that wind erosion in dry regions, often aggravated by human activities, contributes directly to desertification. One example is the "desert pavement", coarser grains left behind by wind erosion forming a hard crust, which prevent any vegetation from growing. In order to improve characterisation and prediction of the wind erosion, it is necessary to understand the basic principles of wind-driven motion of particles in the ABL. This is presented in the next section.

Wind erosion is

Wind-driven particle motion

Sand drift and dust transport are the predominant manifestations of wind erosion. Particles are commonly transported by the following three modes, as illustrated in Figure 1.2:

1. suspension: particle smaller than 60 µm, often denoted as "dust".

2. saltation: particle between 60 µm and 2 mm, mostly sand.

3. traction, or creep: massive particle larger than 2 mm, often denoted as "granules" or "pebbles".

As illustrated in Figure 1.2, most of the wind-driven particle motions occur within a limited layer close to the surface, except for dust particles in suspension with a diameter often smaller than 60 µm [START_REF] Greeley | Wind as a geological process: on Earth, Mars, Venus and Titan[END_REF]. Based on in-situ measurements, the mean saltation height is found to be generally around 20 cm [START_REF] Pye | Aeolian Sand and Sand Dunes[END_REF]. The limiting factor for these particles is the gravity force, which constantly pulls back particles towards the wall, creating an impact.

During its stay in the air, since the wind speed prevails in the longitudinal direction and presents strong shear in the ABL, airborne particles often acquire larger streamwise velocities than the wall-normal ones. This disproportionality is reflected in the grazing angles of saltation trajectories. After rebound, the two velocity components become comparable with each other.

As the experimental results of [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF] show, a great part of horizontal momentum is transferred to the vertical one during the particle-bed interaction at impact. This is probably an effect of particle rotation. The impact of saltating particles also allows dust to be splashed into the atmosphere and subsequently remain suspended due to weak mass. On the other hand, larger grains, too heavy to be entrained by wind, can temporally engage in a slow, creeping motion after collisions with saltating particles. This could lead to an exposure of erodible elements, previously sheltered by the granule, and also a potential fragmentation of non-erodible elements, both contributing to further wind erosion.

Formulation of the physical problem

In the previous sections, we have provided some background information on wind erosion.

In this section, physical aspects of this process are analysed, with a review of relevant studies in literature from different domains. In Figure 1.3, a schematic view on the intimate relationship between flow inside the ABL, topography and particle transport is presented. The interactions between these three factors are discussed in the following sections.

Wind and particle

From the pioneering work of [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF], wind transport of solid particles has received continuous attention. Classical approaches adopting the steady state assumptions of [START_REF] Owen | Saltation of uniform grains in air[END_REF] provide the foundation of our understanding of this aeolian process. Several numerical studies also shed light on the particle saltation [START_REF] Anderson | Wind modification and bed response during saltation of sand in air[END_REF][START_REF] Shao | Numerical modelling of saltation in the atmospheric surface layer[END_REF].

Suspended dust in the atmosphere has been measured and modelled in Marticorena et al. [1997] in the context of a global dust cycle. Experimental studies of a saltation layer in wind tunnel were performed by [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF]; [START_REF] Ho | Aeolian sand transport: Length and height distributions of saltation trajectories[END_REF]; [START_REF] Nalpanis | Saltating particles over flat beds[END_REF] and [START_REF] Taniere | On the behaviour of solid particles in a horizontal boundary layer with turbulence and saltation effects[END_REF]. The objective of these works was to analyse established and stabilised saltation layer and its dependence on particle and flow characteristics.

Whereas the aforementioned studies were focused on the bulk aspects of particle saltation, in recent years researches were aimed at the fundamental nature of particle-laden flows, thanks to the rapid development of computational fluid mechanics (CFD). [START_REF] Marchioli | Mechanisms for particle transfer and segregation in a turbulent boundary layer[END_REF] studied a dilute dispersion of heavy particles using a pseudo-spectral direct numerical simulation (DNS) coupled with Lagrangian tracking of particles. In DNS, all relevant scales going from turbulence production to dissipation are resolved. Particle motion was assumed to be dependent on drag, gravity and lift forces, without feedback effects on the carrier flow. The objective was to investigate experimental observations on particle accumulation in the wall region. A link between the transfer of particles to, and away from the wall and the coherent structures present in the boundary-layer was established in this study. Recently, Vinkovic et al. [2011] used DNS and found that particles moving away from the wall are surrounded by ejection-like flow structures. The PhD work of [START_REF] Yu | Finite size particle transport in turbulent channel flow[END_REF] extended the study to include finite-size effects of the particle and found similar conclusions in the case of turbulent channel flow.

Despite its increasing role as a powerful research tool, the high computational cost of DNS makes it impractical for the simulation of complex flows at high Reynolds numbers, which is the case for the geophysical flows involved in aeolian processes. In recent years, large-eddy simulation (LES) for turbulent flow has been applied to the study of aeolian transport. Different from DNS, LES relies on both a resolution of large-scale, energetic turbulent motions and a modelling of smaller scales. Among the first LES of solid particle transport, Vinkovic et al. [2006b] developed a stochastic model for the simulation of passive dispersion, and later Vinkovic et al.

[2006a] applied the model to the simulation of sand saltation over a flat bed. The simulations were compared with the experimental results of [START_REF] Nalpanis | Saltating particles over flat beds[END_REF] and [START_REF] Taniere | On the behaviour of solid particles in a horizontal boundary layer with turbulence and saltation effects[END_REF] and satisfactory results were obtained.

We note two recent works that study the influence of turbulent structures on the spatial distribution of wind-driven particles, one on sand saltation [START_REF] Dupont | Modeling saltation intermittency[END_REF] and another on snow drift [START_REF] Groot Zwaaftink | Modelling Small-Scale Drifting Snow with a Lagrangian Stochastic Model Based on Large-Eddy Simulations[END_REF]. Using a Lagrangian stochastic model, fluctuations in the drifting snow flux are captured by LES in Groot Zwaaftink et al. [2014], and qualitatively compared with field and wind-tunnel measurements. [START_REF] Dupont | Modeling saltation intermittency[END_REF] conducted LES of a time-dependent sand saltation process and focused on its intermittent aspects. This work successfully reproduces "streamers", streamwisely elongated parcels of grains meandering near the surface, which have been observed by in-situ studies [START_REF] Baas | Formation and behavior of aeolian streamers[END_REF].

The omnipresence of topography in geophysical flows is of primary importance to the aeolian process. In this work, a first approximation is made in order to represent the topographic effects by an idealised terrain configuration. In the following sections, terrain effects on the flow are discussed using relevant researches related to this configuration.

Wind and topography

The topography studied here can be generally qualified as "complex terrain", which consists of irregularities of different scales compared to a flat surface. All these different scales are involved in the determination of local wind speed and wall shear. We point out the two scales involved in the study: a micro-and another meso-scale. At the micro-scale, the term "roughness" is often used. In our study, the influence of the roughness on the flow is represented in an average way. At the meso-scale, obstacles and hills exist, with dimensions comparable to a fraction of the ABL. These obstructions greatly modify flow characteristics of the incoming boundarylayer and in some occasions create large recirculation zones where particles can be trapped and deposited.

A bidimensional hill is widely used as an idealised topographical element in studies of flow over wall-mounted obstructions. This particular type of flow is a common backdrop for various physical phenomena, from pollutant dispersion over complex terrain [START_REF] Gong | A wind tunnel study of turbulent dispersion over two-and three-dimensional gentle hills from upwind point sources in neutral flow[END_REF] to the formation of sand dunes on Earth [START_REF] Charru | Subaqueous barchan dunes in turbulent shear flow. Part 2. Fluid flow[END_REF] and on Mars [START_REF] Araújo | Numerical modeling of the wind flow over a transverse dune[END_REF].

Among the numerous applications, the calculation of wind loads induced by the hilly terrain is an influential factor both in the planning of wind farms [START_REF] Carpenter | Investigation of wind speeds over multiple two-dimensional hills[END_REF] and buildings [START_REF] Bitsuamlak | Effects of upstream two-dimensional hills on design wind loads: a computational approach[END_REF].

Flows around various kinds of 2D hills without separation have been studied in the past using analytical, numerical, laboratory and in-situ investigations. [START_REF] Jackson | Turbulent wind flow over a low hill[END_REF] studied theoretically and analytically the disturbances generated by gentle hills on the mean flow by dividing the flow from the canopy to the upper BL into a series of layers with distinct dynamics.

In the proposed linearised theory, important simplifications were made to achieve an analytic solution. Using a rough hill of Gaussian shape, [START_REF] Hunt | Stably stratified shear flow over low hills[END_REF] extended the theory to account for different kinds of ABL stratification. [START_REF] Gong | A wind tunnel study of turbulent flow over model hills[END_REF] conducted wind-tunnel experiments and obtained measurements of the mean flow and turbulence over a cosine-square shaped hill. Their results confirmed the inviscid dynamics of the outer flow perturbations. The authors found comparable speed-up ratios as in the theoretical predictions of [START_REF] Jackson | Turbulent wind flow over a low hill[END_REF]. Using laser Doppler anemometry (LDA) measurements, [START_REF] Almeida | Wake flows behind two-dimensional model hills[END_REF] studied the recirculating flow established in the wake of two-dimensional, polynomial-shaped hills mounted in a water channel. The originality of this experiment lies in the configuration of a hill array, which allowed an original comparison between the flow around a single obstacle with that resulting from multiple hills. Using windtunnel measurements, [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] highlighted the lee-side flow separation behind a steep hill. They found that the flow behaviour is in stark contrast with the one behind a gentle hill, studied in Cao & Tamura [2007]. Moreover, by covering the surface with small cubes, the results of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] show that the separation zone grew as a result of the wall roughness. [START_REF] Chapman | Reynolds stress and sand transport over a foredune[END_REF] conducted in-situ measurements of the Reynolds stress and sand transport on the windward side of a coastal dune. They found characteristic correlations between sand transport and the different quadrants of the Reynolds stress.

Using LES, [START_REF] Dupont | Large-eddy simulation of turbulent flow over a forested hill: Validation and coherent structure identification[END_REF] studied a turbulent flow over a forested hill. Due to the joint influence of the hill and the vegetation canopy, a recirculation zone developed on the lee side. The focus of the authors was on the wake region near the canopy in which the turbulent structures were studied with care through the vorticity analysis and two-point velocity correlations. A weak correlation was found between the wake flow and the upstream flow. Starting from a Gaussian shaped hill, [START_REF] Araújo | Numerical modeling of the wind flow over a transverse dune[END_REF] used a morphodynamic model to achieve an asymmetric sand dune-form, and later conducted simulations by solving Reynolds Averaged Navier-Stokes (RANS) equations based on this dune shape. They found that the length of the separation bubble displayed a surprisingly strong proportionality to the wind intensity. Yet this finding is to be taken with caution, because of the strong role played by the instantaneous turbulent eddies, that are not accounted for by a RANS model, on the flow reattachment [START_REF] Tamura | Numerical study of unsteady wake flows over a hill for the oncoming boundary-layer turbulence[END_REF].

The interaction between wind and topography is an active subject of research in the domain of urban wind-engineering studies as well [Grimmond & Oke, 1999]. Clusters of buildings constructed in cities represent a distinct topographic feature to the overlying ABL. Over these buildings, an urban boundary-layer develops inside the ABL, often with distinct aerodynamic , thermodynamic properties and micro-climate. Readers are referred to [START_REF] Oke | Boundary layer climates[END_REF] for a complete review on this subject.

Depending on the building configuration, several studies have shown that flow over these obstacles can be classified into three regimes: isolated, wake and skimming flows. This is illustrated in Figure 1.4, Oke [1992]. For sparse elements (isolated flow), the interaction between the flow and the roughness elements is at its full extent since individual elements are completely exposed to the flow. As the roughness density increases, the momentum exchange between the roughness and the outer flow diminishes due to the sheltering effect between the elements [START_REF] Raupach | Conditional statistics of Reynolds stress in rough-wall and smooth-wall turbulent boundary layers[END_REF]. For extreme roughness density (skimming flow), the interaction between the flow and the roughness elements vanishes as a new smooth wall emerges. Grimmond & Oke [1999] delimited these three regimes by roughness density, based on the wind tunnel experiments of Hussain & Lee [1980]. These flow regimes have been reproduced experimentally by Simoëns et al. [2007] inside a wind-tunnel-scale boundary-layer using squared obstacles. Wind effect on obstacle is also an important subject in the study of dune migration. This could be achieved by coupling a LES code with the immersed boundary-method (IBM) [START_REF] Le Ribault | Immersed Boundary Method for Large Eddy Simulation and Lagrangian Stochastic Modeling of Passive Scalar Dispersion Downstream of an Obstacle[END_REF], not presented in this thesis, which allows to simulate migrating sand dunes with morphology features dependant on both time and wind intensities.

Topography and particles

As a product of aeolian processes, sand dunes represent one of the most fast-changing landscapes on Earth. Researches on the emergence of these bed forms have been conducted by [START_REF] Claudin | Field Evidence for the Upwind Velocity Shift at the Crest of Low Dunes[END_REF]; [START_REF] Franklin | Subaqueous barchan dunes in turbulent shear flow. Part 1. Dune motion[END_REF]. Readers are referred to [START_REF] Charru | Sand Ripples and Dunes[END_REF] for a review on this subject.

As illustrated in Figure 1.5, the sand mass moves forward through a continuous series of migration by saltation and creep over the stoss side. When sand grains carried by the wind reach the brink of the dune, the downward slope of the dune forces a fraction of them to either participate or trigger an avalanche on the slip face. Fed constantly by sand grains, the slip side of the dune is always on the verge of sliding. As a result, the slip side has generally a larger slope (about 30°) than the gentle windward one (about 18°).

As pointed out by [START_REF] Lancaster | Desert dune processes and dynamics[END_REF], most of field researches is focused on the stoss side of sand dunes, whereas the lee side has received little attention, despite the fact that the flow dynamics in the lee are equally important and present distinct features compared to the upstream et al. [2007] discussed the erosion potentials and geomorphological significance of the interdune region based on measurements of flow dynamics. However, no measurements of particle transport have been conducted by the authors and the dune spacing was not systematically varied. In our opinion, applying the flow classification proposed by [START_REF] Oke | Street design and urban canopy layer climate[END_REF], as described in

Our approach and goals

In the context of this thesis, an idealised topography, consisting of multiple 2D hills facing perpendicularly to the wind direction, located inside a laboratory-scale wind tunnel is adopted.

The hill geometry is symmetric in the streamwise direction and the slope tends smoothly to zero both at crest and feet without any sharp brinks. This is an idealised configuration compared with natural terrain. Nevertheless, symmetrical dune shapes do exist in nature, such as star dunes, in areas of complex wind regimes with strong seasonal changes in wind direction. On the other hand, sand dunes in nature are constantly under the joint influence of wind erosion and deposition of incoming sand grains and thus, present ever-changing shapes. After long periods of time, some large scale characteristic forms do emerge into the classic barchan, parabolic and star dunes, etc. A Gaussian shape thus represents, theoretically, a final geometric outcome at the end of a random overlapping of sand dunes of different shapes. Therefore, here we focus on the reproduction of the aerodynamics of such Gaussian shape using LES with the aim of studying schematic, but representative, recirculation zones in which solid particles transported in the upper layer could be trapped. Results obtained in this work could thus provide insights on to what extent such process participate to dunes migration and evolution.

In the next sections, we present, respectively, the experimental project related to this work and the numerical strategies adopted in this thesis.

Wind-tunnel studies

The experimental campaign presented in this section was conducted in the frame of the NFSC/ANR Sino-French program PEDO-COTESOF, "Particle Emission and Deposition Over Complex Terrain for Soil Fixation". During the experiments, denoted henceforth as PC09, the flow field around one isolated or several Gaussian hill(s) with various spacings was studied inside a wind tunnel [Simoëns et al., 2007] using particle image velocimetry (PIV) measurements. Solid particles were injected from an upstream source in order to study the transport over downstream hill(s). Particle concentration and velocity profiles were obtained using digital image treatment. The locations of the particle source and the Gaussian hills are shown in Figure 1.7. The spacing between two Gaussian hills is varied in order to study the influence of different flow regimes on the transport of particles. The influence of the spacing between obstacles are studied in the context of pollutant dispersion in Simoëns et al. [2007] and [START_REF] Simoëns | The flow across a street canyon of variable width-Part 2:. Scalar dispersion from a street level line source[END_REF], where passive scalar dispersion between squared obstacles was studied. Within this campaign, two objectives were fixed: first, to obtain some information on the way the hilly terrain modifies the concentration levels of solid particles and second, to provide a data basis for numerical validation. These points will be tackled through this manuscript. For more details on the experimental set-up, readers are referred to [START_REF] Simoëns | Influence of Gaussian hill on concentration of solid particles in suspension inside turbulent boudary layer[END_REF].

LES

In this thesis, LES is used as the numerical tool to investigate solid particle transport. Instantaneous characteristics of the flow are numerically represented in LES, since large, energetic fluid structures are resolved contrary to other "eddy-modelling" approaches such as RANS. This is an important aspect since the instantaneous influence of the flow is crucial in particle entrainment, as well as particle transport and trapping. From an aerodynamic perspective, since flow separation is largely instantaneous and intermittent, LES can better capture the large eddies generated behind the hill that has presumably a non-negligible influence on particle trapping.

Moreover, the general benefits of numerical simulation also apply to our work, such as the ease of conducting parametric studies by varying systematically the relevant non-dimensional numbers.

The Lagrangian particle-following method is adopted for the calculation of particle trajectories at each time step. The motion equation of particle is based on a point-wise hypothesis, which assumes essentially that eddies with sizes smaller than the grain are negligible. A twoway modelling is used to account for the feedback from particles to the flow. The numerical representation of the upstream, specially devised particle source allows a constant entrainment flux of grains into the flow.

The LES code ARPS, originally developed at the Center for Analysis and Prediction of Storms (CAPS), University of Oklahoma [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF], is used in this thesis for the study of particle transport. Previously, ARPS has been used for the development of a stochastic subgrid model for the study of passive scalar dispersion in Vinkovic et al. [2006b] as well as sand saltation over a flat surface in Vinkovic et al. [2006a] or [START_REF] Dupont | Modeling saltation intermittency[END_REF]. Recently, the code has been completely parallelised, notably for the calculation of solid particle motion, thanks to the work of S. Dupont at INRA, Bordeaux. This thesis is part of the continuous development of this numerical tool.

Objectives and outline

The aim of this thesis is to understand processes related to the problem of wind erosion and particle deposition over an idealised 2D, transverse hill. Concretely, the objectives of this thesis are 1. to conduct a numerical study and to compare with the experiments. To this aim, relevant physical models, such as particle lift-off, grain interaction with the surface and the flow are implemented in ARPS. Simulations of the complete sequence of a developing saltation process, including aerodynamic entrainment and particle rebound are conducted over a 2D-hill topography.

2. to conduct parametric studies using relevant non-dimensional numbers on flow characteristics and on particle trapping inside recirculation zones.

The structure of this manuscript is as follow. In Chapter

Chapter 2

Governing equations and numerical methods

The aim of this thesis is to simulate a spatially-developing turbulent boundary-layer (TBL) at the scale of a laboratory wind tunnel in order to study particle transport over Gaussian hills. The effects of the topography and wall roughness are considered in this study. An illustration of the three-dimensional domain used in the simulation of the TBL is shown in Figure 2.1. x denotes the streamwise direction in which the TBL develops, y the spanwise and z the wall-normal directions, respectively.

In this chapter, governing equations of fluid motion inside the TBL and their LES formulations are presented in Sections 2.1 and 2.2, respectively. Details of the numerical implementation such as the space discretization, the time integration and the boundary conditions are presented in Section 2.3. A wall model is employed in order to compensate for unresolved physics due to the coarseness of the near-wall grids. The details are given in Section 2.4. Finally, the generation of inflow fluctuation data, crucial to the simulation of the TBL is discussed in Section 2.5.

Figure 2.1: Sketch of the computational domain defined in a Cartesian coordinate system .x; y; z/. L x , L y and L y are the streamwise, the spanwise and the wall-normal extents of the domain.

Equations and assumptions

The continuity equation is derived by applying the principle of the mass conservation to a fluid parcel passing through an infinitesimal, fixed control volume. This yields [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]]

1 D Dt C r u D 0 ; (2.1)
where is the fluid density and u is the fluid velocity vector. The material derivative operator D=Dt for a variable A is defined by

DA Dt Á @A @t C u rA : (2.2)
DA Dt gives the derivative of the variable A related to a fluid particle in motion.

The Navier-Stokes (NS) equations describe the balance of forces on a fluid parcel. Among them are the body forces, which combine the gravity, characterised by the gravitational acceleration g, and the Coriolis force due to the rotation of the earth of angular velocity Ω. On the other hand, the contact forces (pressure, viscous friction) can be expressed by the stress tensor . The NS equations write

Du Dt D .g 2Ω u/ C 1 r : (2.3)
The stress tensor is related to the strain tensor by means of a constitutive equation. The strain-rate tensor S is the symmetric part of the velocity gradient tensor @u i =@x j ,

S ij D 1 2  @u i @x j C @u j @x i à : (2.4)
Here we use u i D .u 1 ; u 2 ; u 3 / D .u; v; w/ and x i D .x; y; z/ for convenience. A repeated subscript implies summation.

The constitutive equation is of the simplest form for the Newtonian fluid, in which the viscous stress depends linearly on the rate-of-strain S . The stress tensor is expressed as a sum of an isotropic part (I) related to the pressure p, and another deviatoric part (II):

ij D pı ij " ƒ' … I C .2S ij 2 3 @u j @x j ı ij / " ƒ' … II ; (2.5)
where D is the dynamic viscosity of the fluid, being the kinematic viscosity. In case of negligible spatial variation of temperature, can be assumed as homogeneous. The Kronecker tensor is defined by

ı ij D ( 1 if i D j ; 0 otherwise : (2.6)
Combining Equations 2.3 and 2.5 gives the momentum equation for a compressible fluid in a rotational frame

@u i @t C u j @u i @u j D 1 @p @x i 2" ij k ˝j u k ı i3 g C @ @x j  2S ij 2 3 @u j @x j ı ij à : (2.7)
The second term on the RHS of Equation 2.7 is the Coriolis force, defined using " ij k which is the Levi-Civita symbol By replacing the pressure and the gravity terms in Equation 2.7 with the RHS of Equation 2.12, the momentum equation can be simplified into:

" ij k D 8
@u i @t C u j @u i @x j D 1 r @ @x i p r ı i3 g 2" ij k ˝j u k C @ x j .2 S ij / ; (2.13)
We note that the effect of density variation on the fluid under the Boussinesq approximation is only modelled through the term . = r / g, namely the buoyancy effect.

Among the three thermodynamic variables, , and p, two should be predicted and the other derived from the equation of state. In ARPS, the pressure p is directly solved due to the crucial role of the pressure gradient in the momentum equation (Equation 2.13). Using the Boussinesq approximation (Equation 2.12), this equation is obtained by taking the material derivative of the equation of state and replacing the time derivative of density by the velocity divergence using the mass continuity equation (Equation 2.1)

@ p @t C u j @ p @x j D r c 2 a  1 @ @t @u j @x j à C r gw ; (2.14)
where c a is the speed of sound in the air. The atmospheric boundary-layer (ABL) studied here being of neutral state, the potential temperature variations are neglected in this thesis.

Modelling approach of LES

Filtered equations

In LES, grid-scale filtering is applied to the flow velocity and the thermodynamic variables in order to obtain the filtered momentum equations. where D . x y z/ 1=3 is a measure of the filter width determined by the mesh sizes x, y and z in each direction.

A.x; t/ can thus be decomposed into a large-scale value and a small scale fluctuation

A D e A C A 00 :

(2.17)

By applying the filter to the momentum equation, LES aims to solve the filtered variables instead of the "true", unfiltered variables. The spatially filtered NS equations then write @e u i @t C e u j @e u i @x j D

1 e r @ @x i . e p ˛d e r Q u j @x j / e r ı i3 g 2" ij k ˝j e u k C @ @x j .2 f S ij R sgs ij / ; (2.18)
where ˛d corresponds to the artificial "divergence damping" terms designed to attenuate acoustic waves [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF].

The subgrid stress tensor R sgs ij arises from the filtering of the non-linear convective term A u j @u i @x j on the LHS of Equation 2.7:

R sgs ij D e u i u j Q u i Q u j : (2.19)
As a subgrid variable acting on the filtered flow field, R sgs ij characterises the momentum exchange between the larger scales and the subgrid ones. Since R sgs ij cannot be expressed using resolvedscale variables, it is obtained by turbulence closure models, presented in Section 2.2.2.

The filtered pressure equation (Equation 2.14) writes

@ f p @t C e u j @ f p @x j D r c 2 a 1 f r @ f r @t @e u j @x j ! C r ge w ; (2.20)
with the omission of several terms such as A u j @ p @x j e u j @ e p @x j .

Subgrid-scale modelling

Subgrid viscosity hypothesis

One common type of SGS modelling is the subgrid viscosity model. The concept of the subgrid viscosity is similar to the eddy-viscosity assumption used in the mixing-length theory proposed by [START_REF] Prandtl | Uber die ausgebildete Turbulenz[END_REF]. The mixing-length theory claims that the effect of the turbulence on the mean flow may be obtained by replacing the laminar viscosity with an "eddy viscosity".

In the framework of LES, the SGS viscosity model postulates that the energy transfer from the large scales to the small ones is similar to the molecular diffusion, which involves the concept of the SGS eddy viscosity.

In a similar form to the constitutive equation for the stress tensor in Equation 2.5, the deviatoric part of the SGS stress R sgs ij is modelled using the subgrid viscosity sgs as

R d ij D sgs .2 f S ij 2 3 @ Q u j @x j ı ij / ; (2.21) where R d is the deviator of R sgs R d ij D R sgs ij 1 3 R sgs kk ı ij : (2.22)
The isotropic part of the SGS stress 1 3 R sgs kk ı ij is proportional to the identity tensor. This term can be generally combined with the pressure term in Equation 2.18 and thus requires no explicit modelling [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF].

Contrary to the molecular viscosity, the eddy viscosity is not a propriety of the fluid but of the flow. In order to evaluate sgs , the Smagorinsky model is based on the local equilibrium hypothesis, which states that the flow is in constant spectral equilibrium and there is no accumulation of energy at any frequency. The Smagorinsky model assumes that this relationship is valid instantaneously and locally in the flow. The local use of this relationship is not theoretically justified, since it only ensures that the energy transfers through the cutoff are expressed correctly on average, and not locally. In the Smagorinsky model the subgrid viscosity is given by [START_REF] Le Ribault | LES of the compressed Taylor vortex flow using a finite volume/finite element method on unstructured grids[END_REF] sgs D .C s / 2 j Q Sj ;

(2.23)

where j Q S j D .2 f S d ij f S d ij / 1=2
and f S d is the deviatoric part of Q S on the RHS of Equation 2.21.

For C s , namely the Smagorinsky coefficient, [START_REF] Lilly | The representation of small scale turbulence in numerical simulation experiments[END_REF] proposed C s D 0:17 by a study of homogeneous isotropic turbulence. For inhomogeneous flows such as the TBL, it is generally accepted that this value is too dissipative. Thus in wall flow simulations, the Smagorinsky coefficient is generally decreased, e.g., to C s D 0:1 [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF][START_REF] Piomelli | Model consistency in large eddy simulation of turbulent channel flows[END_REF].

One of the problems with the Smagorinsky model is that the value of C s is not universal. In order to adapt the model to the local structure of the flow, [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] proposed a dynamic procedure that adjusts C s locally at each time step. By this dynamic procedure C s becomes a spatially and temporally varying flow parameter.

order turbulent kinetic energy transport equation

Unlike the Smagorinsky model, in the 1.5-order turbulent kinetic energy (TKE) closure [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a three-dimensional model[END_REF][START_REF] Moeng | A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence[END_REF], the subgrid viscosity sgs is not related to the filtered strainrate tensor. By considering one additional variable, the subgrid kinetic energy k sgs , the 1.5-order TKE model provides more information on the subgrid modes than the traditional subgrid viscosity model, bypassing the use of local equilibrium hypothesis. The filtered subgrid kinetic energy k sgs is defined as

k sgs D A u 00 k u 00 k 2 : (2.24)
k sgs needs to be resolved in addition to the filtered NS equations by a transport equation. The exact expression of k sgs is obtained from a second filtering of the product of u 00 i and the subgrid momentum equation, which itself is obtained by subtracting the filtered momentum equation (Equation 2.18) from the unfiltered NS equation (Equation 2.7). The transport equation of k sgs in its simplified form [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF] writes

@k sgs @t D Q u j @k sgs @x j " ƒ' … I R d ij e S ij " ƒ' … II C @ @x i  2. C sgs / @k sgs @x i à " ƒ' … III "ƒ'… IV : (2.25)
The terms on the right-hand side of Equation 2.25 represent, respectively, advection (I), production by shear stress (II), diffusion due to viscous as well as subgrid viscosities (III) and subgrid dissipation due to viscosity (IV).

The subgrid turbulent kinetic energy dissipation in IV writes, by definition

D B @u i @x j @u i @x j @ Q u i @x j @ Q u i @x j ! : (2.26)
Using dimensional analysis, is evaluated by

D C k 3 2
sgs ;

( 2.27) where the constant C takes the value of 3:9 in the wall-adjacent grids and 0:93 above, according to [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a three-dimensional model[END_REF] and [START_REF] Moeng | A Large-Eddy-Simulation Model for the Study of Planetary Boundary-Layer Turbulence[END_REF].

Finally, using the mixing length hypothesis, the subgrid eddy viscosity sgs is evaluated as the product of the characteristic velocity fluctuation ( p k sgs ) and a characteristic length of the subgrid turbulence l m sgs D C s l m q k sgs :

(2.28)

Since the mesh size is a good representation of the smallest scales present in the resolved flow, the mixing length l m is defined as

l m D (p
x y for horizontal motion, z for vertical motion.

(2.29)

According to [START_REF] Moeng | Evaluation of turbulent transport and dissipation closures in second-order modeling[END_REF], C s D 0:1. Equation 2.28 combined with Equations 2.21 and 2.25 completes the 1.5 order TKE closure.

Numerical method

In order to simulate a TBL flow with or without obstacles, we use the Advanced Regional Prediction System (ARPS), originally developed by the Center for Analysis and Prediction of Storms (CAPS) at the University of Oklahoma. This model was developed in the framework of numerical weather prediction of meso-scale meteorology and for general computational fluid dynamics applications. A detailed description of the standard version of ARPS and its validation cases are available in the ARPS User's Guide [START_REF] Xue | ARPS version 4.0 Users Guide[END_REF] and in [START_REF] Xue | The Advanced Regional Prediction System (ARPS) -A multi-scale nonhydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification[END_REF] and [START_REF] Xue | The Advanced Regional Prediction System (ARPS) -A multi-scale nonhydrostatic atmospheric simulation and prediction tool. part II: Model physics and applications[END_REF]. The following sections give several details of the numerical implementation in the context of the current thesis.

Discretization, mesh and parallelisation

Time discretization

Since ARPS resolves the NS equation with the Boussinesq approximation (Equation 2.18), acoustic waves are numerically present in the simulation, which severely limits the time step size of explicit time integration schemes. It is clear that high-frequency acoustic waves are not of major importance to the study of aeolian particle transport. To improve the model efficiency, a mode-splitting time integration technique presented in [START_REF] Klemp | The Simulation of Three-Dimensional Convective Storm Dynamics[END_REF] is employed. This technique divides the big integration time step into a number of computationally inexpensive small time steps and updates the acoustically active terms every small time step while computing all the other terms only once every big time step. Consequently, only a small part of the numerical resolution governed by the small time step is limited by the acoustic wave speed. The large time-step integration, using a centred leapfrog time differencing scheme [START_REF] Fletcher | Computational Techniques for Fluid Dynamics 1[END_REF], concerns mainly the potential temperature and the horizontal velocities. For small time steps, a second-order centred implicit "Crank-Nicholson" scheme is used.

This concerns notably the pressure and the vertical velocity, the "acoustically active" terms, in Equation 2.20.

Spatial discretization

The partial differential equations described in Sections 2.2.1 and 2.2.2 are spatially discretised using a fourth order quadratically-conservative scheme for the advective terms, and a second order finite-difference scheme for the other terms.

The discretised variables are defined on an Arakawa-C staggered grid [START_REF] Arakawa | Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow[END_REF] in ARPS. This particular type of grid configuration imposes that variables are not necessarily defined at identical spatial locations even if they share the same grid index .i; j; k/. Scalar variables, such as temperature and pressure, are defined at the centre of the cell. Locations of the coordinate variables x, y and z are staggered, as they are defined at the centre of the corresponding cell faces (Figure 2.2). The velocity components u,v and w are defined according to x, y and z. It follows that, e.g., the streamwise velocity gradient @u=@z is evaluated on the z plane, half a grid interval below the u point. et al., 1995].

Terrain-following mesh and vertical grid stretching

In ARPS, the mesh is constructed using a terrain-modified coordinate system. The key advantage of this coordinate system is that it adapts to the geometry of the underlying terrain near the ground, in the same way as the wind aligns to the surface, and reduces to a rectangular Cartesian system at a sufficient height above the hill, where the mean flow is largely horizontal.

Hence, it retains the advantages of both terrain-following and rectangular coordinate systems in the appropriate regions.

The transformed coordinate system ( , Á, ) is related to the Cartesian one by

D x Á D y (2.30) D .x; y; z/:
We note that the constant and Á surfaces are the same as those of constant x and y. The vertical coordinate of the terrain-modified coordinates, , is given by

D 8 < :.z flat z ref / z h z flat h C z ref if z ref Ä z Ä z flat ; z if z > z flat ;
( 2.31) z flat is the altitude from which coincides with z, and z ref is the reference altitude for , equal to 0 in this work. h.x; y/ is the elevation of the terrain. Note that the grid line D 0 is both the bottom boundary of the computational domain and a physical boundary.

In ARPS, the generated mesh is uniform in x and y directions. It is possible to refine the mesh in the wall-normal (z) direction. A vertical grid stretching is applied to the lower half of the domain in order to refine the mesh points. Grid size is varied according to a hyperbolic tangent function, which accommodates the coarseness of the upper domain to the finer grids in the near-wall region. An illustration of the computational domain is given in Figure 2. 2.30) and vertical grid stretching (Equation 2.31). Note that only every second grid line is shown for clarity.

Parallelisation

The implementation of ARPS on parallel computing architectures is based on the strategy of domain decomposition. The method consists in assigning subdomains of the full computational grid to separate processors, while ensuring the balance of work load between processors. In principle, a minimum amount of global information is required at each grid point and interprocessor communications are established at the boundaries of the subdomains.

In ARPS, the domain decomposition is employed in the horizontal directions, in order to comply with the direction of the mean flow inside the TBL. As shown in Figure 2.4, grid points located on the border areas of a subdomain require information from the adjacent subdomain at each time-step. This communication is necessary, e.g. for the calculation of the spatial finite differences. Values at the outer border are supplied by the neighbouring processor using the messages passing interface (MPI) between processors. In order to avoid redundancy in the exchange of information between the outer and the inner borders, the outer border data are stored in the local memories of the corresponding subdomain after each communication step.

Boundary conditions

At the boundaries of the computational domain, it is necessary to impose physically meaningful and numerically stable values to the discretised variables in order to enable the approximation of flux, gradients, etc., using the finite difference technique. In ARPS, the boundary conditions are enforced through the use of extra grid points defined outside the physical boundary of the domain, which are often referred to as the "fake" points. In the simulation of an ABL, only the lower boundary of the domain is physical. For the top and lateral boundaries, different types of boundary conditions are available in ARPS. The ones adopted in this thesis are presented here.

Bottom wall: rigid wall condition

On a staggered grid cell (Figure 2.5) , the rigid wall condition is reinforced by u 0 D u 1 for the streamwise velocity; v 0 D v 1 for the spanwise velocity; w 0 D w 2 and w 1 D 0 for the wall-normal velocity.

The detailed treatment of the wall boundary condition is further presented in Section 2.4.

Lateral wall: periodic condition

This choice is related to our assumption that the wind flow is statistically homogeneous in the spanwise direction.

Top boundary: zero-normal gradient

This is generally done by imposing cells at the top boundary to take the value of ones immediately interior to them.

Inlet

At the inlet, the inflow generation algorithm provides boundary values at each time step. This technique is further described in Section 2.5.

Outlet: radiative condition

This condition is described at the end of the section. For the outlet, an open boundary condition is used. This "radiative" condition is designed to allow waves created inside the domain to exit freely through the outlet with minimal repercussion [START_REF] Durran | The effects of moisture on trapped mountain lee waves[END_REF]. Concretely, radiation boundary conditions typically employ a simplified wave propagation equation in order to update accordingly the predicted boundary values at the outlet. In ARPS, the formulation of [START_REF] Orlanski | A simple boundary condition for unbounded hyperbolic flows[END_REF] is used for this condition. Due to the configuration of the physical problem, (a solid wall plus an upstream TBL at the beginning of the domain), the wall and the inlet boundaries need specific treatments, detailed in the following sections.

Near-wall treatment

The near-wall treatment is a set of specific procedures available in numerical simulations such as RANS and LES. Its main function is to tackle with the numerical errors that arise from the coarseness of the mesh near the physical wall. Through the use of wall models, the nearwall treatment relies on the modelling of the inner layer of the TBL in a Reynolds-averaged sense. It is a trade-off strategy between the computational cost of the simulation and the quality of the results. The wall modelling is especially necessary for complex flow situations at high Reynolds number, such as a TBL over wall-mounted obstacles, where the computational cost is usually prohibitively high. Here, the near-wall dynamics are of particular importance, since the aerodynamic entrainment, the rebound, and the splash of particles all take place near the wall.

Effects of the near-wall flow, e.g., the wall shear stress, need to be correctly predicted in order to support the study of solid particle transport in the TBL.

Evaluation of the wall shear stress

One particularity of the boundary-layer (BL) flow in comparison to other types of flow such as free shear flow is the existence of the inner layer, which is subject to the direct influence of the wall. The depth of this layer is commonly given by z < 0:1 ı [Pope, 2000], where ı is the 99% boundary-layer thickness. Outer-layer parameters such as ı and the free-stream velocity U 1 have minor influence here, whereas the viscosity reigns in this region because of the large velocity derivatives due to the non-slip condition at the wall. The second parameter in the inner layer is the wall shear stress w . We derive the expression of w from the momentum equations of the turbulent boundary-layer (TBL) flow. For simplicity, the flow field u i is henceforth noted as .u; v; w/.

Consider a two-dimensional turbulent boundary-layer flow .u; w/, after assuming that all spanwise variations are negligible compared to the variations in the other two directions. The stationary momentum equation for the mean streamwise velocity writes hui @hui @x C hwi

@hui @z D 1 @hpi @x C @ 2 hui @z 2 @ @z hu 0 w 0 i : (2.32)
The instantaneous velocity u has been divided into the Reynolds average hui and the corresponding fluctuation u 0 using the Reynolds decomposition u D hui C u 0 ;

( 2.33) where hui D

2T

Z T T udt :

(2.34)

The averaging period T is chosen to be sufficiently large compared to the relevant turbulent time-scales.

While the first term on the right hand side of Equation 2.32 corresponds to the pressure gradient, the second and the third terms can be combined into a total shear stress

Á @hui @z " ƒ' … I hu 0 w 0 i " ƒ' … II ; (2.35)
which is the sum of the viscous shear stress (I) plus an additional stress related to turbulent fluctuations, denoted as the Reynolds stress (II). For laminar flow, e.g., the total shear stress is equal to the viscous shear stress since the Reynolds stress is negligible.

Substituting Equation 2.35 into Equation 2.32 gives hui @hui @x C hwi @hui @z D 1 @hpi @x C

1 @ @z :

(2.36)

Towards the wall, the Reynolds stress tends to zero due to the non-slip boundary condition u D 0. At the wall, Equation 2.35 gives w D @hui @z ˇw :

(2.37)

The resulting wall shear stress, w Á j zD0 , is the viscous friction exerted on the fluid by the wall. w fixes the boundary value of the total shear stress across the inner layer. Besides, it determines the velocity gradient at the wall, which is non-zero and maximum for attached boundary-layer flow.

A velocity scale derived from w , appropriately named as the friction velocity u , writes

u Á p w = ; (2.38)
which is the velocity scale of the inner layer. Using and u , an appropriate length scale can be as well constructed as ı D =u , denoted as the viscous length scale.

Based on u and ı , two non-dimensional parameters can be obtained:

the distance from the wall measured in viscous length scale, or in "wall units"

z C Á z ı D u z ; (2.39)
the streamwise velocity scaled by the friction velocity

u C Á u u : (2.40)
The universal behaviour of the inner layer could thus be described by a law-of-the-wall generally of the form hu C i D f w .z C /, described in more details in Section 2.4.2.

If the computing of inner layer dynamics is set as an objective, this type of LES is qualified as "wall-resolved large-eddy simulation " (WRLES). In this case, no near-wall treatment is involved since no a priori knowledge of the flow is used in the simulation. It is acceptable to impose numerically the non-slip condition u D 0 in the wall-adjacent grids. Note that by construction z C is a local Reynolds number that estimates the balance between viscous and inertial effects. As the Reynolds number decreases towards the wall, the size of the energetic eddies is also reduced in the inner layer. In order to properly capture the near-wall dynamics of the inner layer over a smooth wall, a sufficiently fine mesh becomes thus necessary. [START_REF] Piomelli | Large-Eddy Simulations: Theory and Applications[END_REF] gave the following recommendations: z C 1 < 2, x C ' 50 150 and y C ' 15 40. Generally, the wall adjacent cells should lie in z C ' 1 in order to locate the first computation point inside the viscous sublayer.

In case of the staggered grids, introduced in Section 2.3.1, the non-permeable condition is applied to the bottom wall, which fixes the wall-normal velocity to zero: w D 0. For the other two velocity components, u and v, since they are not located directly at the wall boundary, the non-slip condition implies that the velocity gradient can be approximated by the following expression @u @z ˇw D u 1 z 1 ;

( 2.41) where z 1 corresponds to the centre of the wall-adjacent mesh, as illustrated in Figure 2.6. Substituting Equation 2.41 into Equation 2.37 gives the value of w .

While Equation 2.41 is widely used in DNS and wall-resolved LES for the determination of the wall velocity gradient and w , its range of application is limited. As shown in Figure 2.6, the error induced by Equation 2.41 in the determination of w is proportional to the grid size near the wall. This becomes problematic if the wall flow is subject to high Reynolds number, which is often the case in the majority of engineering applications. In fact, to fulfil z C 1 ' 1 at high Re requires that z 1 decreases. The resulting grid refinement in turn increases the number of computational points: [START_REF] Chapman | Computatisnal Aerodynamics Development and Outlook[END_REF] estimated that the number of points required to resolve the wall layer scales with Re 1:8 . This constraint makes the computational cost of DNS, even WRLES, prohibitively high in the computation of wall-bounded flow at large Reynolds number.

In order to circumvent this limitation, another widely used approach in LES is to model the near-wall dynamics instead of resolving them directly. In the next two sections we discuss the details of this strategy.

Wall stress model based on law-of-the-wall

WRLES of the TBL at high Reynolds number is particularly costly due to the stringent requirement on the near-wall grids. An alternative consists in modelling the near-wall dynamics instead of resolving them. This type of LES is generally denoted as the "wall-modelled largeeddy simulation" (WMLES). The main advantage of this approach is that the WMLES allows the first grid point to be placed inside the logarithmic layer (in practice, 30 Ä z C Ä 200), which is much less demanding than WRLES.

The first problem that arises in WMLES is that the non-slip condition and Equation 2.41 are no longer applicable due to the coarseness of the grid. As pointed earlier, the underestimation of w by Equation 2.41, illustrated by the difference between the blue and the red dashed line in Figure 2.6, is all the more significant when the mesh size is coarse.

In practice, the WMLES employs the wall stress models in order to determine, in a more sophisticated way than Equation 2.41, the wall shear stress w from the computed values provided by the LES. This is generally achieved by applying the non-dimensional law-of-the-wall that predicts the mean flow dynamics of the flow in the inner layer.

One classic formulation of the law-of-the-wall is proposed by von [START_REF] Von Kármán | The analogy between fluid friction and heat transfer[END_REF], often cited as the three-layer model

hui C D 8 < :z C if 0 < z C < 5 ; 5 ln z C 3:05 if 5 < z C < 30 ; 2:5 ln z C C 5:5 if z C > 30 (2.42)
Each line in Equation 2.42 corresponds to one of the three separate sublayers present in the inner layer of the TBL over a smooth wall:

1. the viscous sublayer in which the viscous effects reign;

2. the buffer layer;

3. the inertial sublayer in which neither the viscosity nor the large-scale motion effects are preponderant.

Alternatively, Spalding [1961] proposed a composite velocity profile which is a power-series interpolation scheme joining the viscous sublayer to the logarithmic region

z C Dhui C C exp. ÄB/ Ä exp.Ähui C / 1 Ähui C .Ähui C / 2 2 .Ähui C / 3 6
.Ähui C / 4 24 ;

( 2.43) where Ä is the von Kármán constant. We note that Spalding's law could not be transformed analytically into the form f w .z C /. Nevertheless, Equation 2.43 allows a continuous transition inside the buffer layer, between the viscous sublayer and the logarithmic layer, as illustrated in Based on the law-of-the-wall, the wall stress models approximate the instantaneous dependence between u , w and u i in a Reynolds-averaged sense. The simplification makes the wall stress models relatively straightforward to implement. It has been widely adopted since the pioneering LES study of [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF] on plane channel flow.

Calculation of the surface flux in ARPS

In ARPS, the evaluation of the wall shear stress is designated as the "parameterization of the surface flux". This is in accordance with the reference of the viscous flow study as a problem in momentum transport in the field of geophysical and chemical engineering [START_REF] White | Viscous Fluid Flow[END_REF].

The wall shear stress, w D u 2 , is regarded as a momentum flux through the inner layer and absorption by wall friction. In ARPS, the following relationships are used in order to evaluate the streamwise and the spanwise components of the wall shear stresses: 13 j w and 23 j w

13 j w D C dm U s u ; 23 j w D C dm U s v ;
( 2.44) where u and v are the two horizontal velocity components evaluated at the lowest grid level above the physical boundary and

U s D p u 2 C v 2 ˇz1
, namely the total horizontal wind speed at the lowest grid point. The bulk aerodynamic drag coefficient, C dm D .u =U s / 2 , is a nondimensional parameter. From Equation 2.44, the predicted w .x; t/ is proportional to the square of the local velocity U s .x; t /.

The formulation of [START_REF] Byun | On the analytical solutions of flux-profile relationships for the atmospheric surface layer[END_REF], valid in the ABL under neutral stratification, is used to determined the value of C dm in ARPS:

C dm D Â Ä ln .z 1 =z 0 // Ã 2 ; (2.45)
where z 0 is the roughness length. Equation 2.45 is based on the assumption that the wind speed in wall-adjacent cells, U s , follows the log-law instantaneously. In fact, by writing w Á q

2 13 C 2 23 ˇw D u 2 ; (2.46)
it is not difficult to derive from Equations 2.44 and 2.45

U s u D 1 Ä ln  z 1 z 0 à ; (2.47)
which is another expression of the log-law as in Equation 2.42.

One drawback of surface flux models is the ad hoc coefficients such as C dm and z 0 used as input parameters to the LES. Their values are in general fixed a priori from tabulated data depending on the surface type: water, soil, forest, etc. Although such approach is commonly used in the simulation of geophysical flows [START_REF] Byun | On the analytical solutions of flux-profile relationships for the atmospheric surface layer[END_REF], for the simulation of wind tunnel-scale TBL, more sophisticated wall models are necessary in order to achieve satisfactory results.

The wall model of Duprat et al. [2010] implemented in ARPS

Temmerman et al. [2003] reported that near-wall treatment is more influential than the subgrid modelling on the quality of the simulation results of flow separation on periodical hills, such as the locations of the separation and the reattachment points. It is not surprising that the wall models presented in the previous section do not perform well in predicting separated flows, since they are based on the equilibrium flow assumptions that break down when flow separates.

This also underlines the difficulty arising from the lack of universal scaling law for separated flow.

Steep, hilly surface with rounded shape, as the configuration studied in this thesis, further complicate the matter. In this scenario the flow separation, if present at all, is hardly stable, largely due to the irregularity of the separation and reattachment lines in space and time. As found in the experiments of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], the near-wall flow constantly detaches and reattaches even at the centre of the separation bubble formed behind a two-dimensional hill. This is in contrast with configurations where the separation point is imposed and stable, like in the backward-facing step [START_REF] Cabot | Near-wall models in large-eddy simulations of flow behind a backwardfacing step[END_REF].

Recently, there has been an increasing interest in the improvement of the classical wall stress models on separated flows [START_REF] Chen | Wall Modeling for Implicit Large-Eddy Simulation[END_REF][START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF][START_REF] Manhart | Near-wall scaling for turbulent boundary layers with adverse pressure gradient[END_REF]. A modified wall model was proposed in [START_REF] Manhart | Near-wall scaling for turbulent boundary layers with adverse pressure gradient[END_REF] by integrating the pressure into the estimation of the wall shear stress w . Depending on the nature of the pressure gradient, favourable or adverse, the boundary-layer flow is stabilised or, on the contrary, destabilised and probably runs into separation. In the formulations of the previous wall stress models presented in Section 2.4.2, the pressure gradient is not considered despite its essential role in the separation process, which could be problematic.

More recently, [START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF] extended the formulation of [START_REF] Manhart | Near-wall scaling for turbulent boundary layers with adverse pressure gradient[END_REF] and proposed a wall model capable of estimating w based on the informations of the inertial sublayer, instead of the viscous sublayer. This is of great practical interest, since, as we discussed in Section 2.4.1, the wall-adjacent grids in WMLES are generally located outside the viscous sublayer (z C > 5). In the following paragraphs we briefly present the wall model of [START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF].

By neglecting the convective terms close to the wall [Wang & Moin, 2002], Equation 2.36 is simplified into @ @z D @hpi @x :

(2.48)

By further assuming that the pressure gradient varies little in the wall-normal direction inside the wall-adjacent cells, an integration of Equation 2.48 with respect to z gives w D @hpi @x z 1 ;

D . @hui @z hu 0 w 0 i/ ˇz1 @hpi @x z 1 ;
where z 1 is related to the centre of wall-adjacent cells.

In [START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF], a turbulent eddy coefficient t , coupled with a damping function, is proposed. This function relies on a modified van Driest formula that takes into account the pressure gradient. Using t , the Reynolds stress can be related to the mean velocity gradient, in a similar way as to the viscous shear stress:

w D . C t / @hui @z ˇz1 @hpi @x z 1 : (2.49)
To apply Equation 2.49 in LES, one assumes that the filtered velocity Q u is equivalent to the averaged velocity hui close to the wall. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] showed that this assumption holds if cells are coarse enough close to the wall to contain sufficient near-wall eddies and if the time step is much larger than their time-scales, which is generally satisfied in WMLES.

The wall model of Duprat was first validated by a priori comparison to DNS data of channel and separated flows [START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF]. Subsequent LES showed that it allows for a good prediction of the wall shear stress both with and without streamwise pressure gradient, even when a very coarse grid is used. Readers are referred to the thesis of [START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF] for further details.

Initialisation and generation of turbulent inflow data

A spatially inhomogeneous and time-dependent inflow is necessary for the simulation of a TBL, since a non-periodic boundary condition should be imposed due to the physical configuration of the problem [START_REF] Lee | Direct numerical simulation of the turbulent boundary layer over a cube-roughened wall[END_REF][START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF][START_REF] Tabor | Inlet conditions for large eddy simulation: A review[END_REF]. The minimum requirement necessary to this aim is that the generated inflow velocity should satisfy the prescribed mean flow evolution as well as one-point, second-order statistics such as the Reynolds stress. In this thesis, a three-step procedure is adopted in order to fulfil this requirement. First, the computational field is initialised using a non-dimensional law-of-the-wall. The details are given in Section 2.5.1. A random fluctuation generation is then activated in order to inject flow variations at the inlet during the initial phase of the simulation. This method is presented in Section 2.5.2. Later, the inflow generation is taken over by the extraction/rescaling technique developed by [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF], presented in Section 2.5.3. The whole procedure allows the development of a realistic turbulent flow behaviour inside the whole computational domain starting from a constant, base flow state.

Initialisation

The time-dependent variables present in the discretised equations in ARPS needs to be initialised. The values of the thermodynamic variables p, and are determined by prescribed base state relationships. For the three velocity components, the initialised profile should be compatible with the mean flow behaviour across the depth of the BL.

According to Pope [2000], the mean velocity profile over the whole boundary-layer can be well represented by the sum of two functions, the law-of-the-wall f w .z C /, and a law of the wake, W .z=ı/. The latter represents physically the influence of outer structures on the boundary-layer.

Thus u C .z C ; z=ı/ D f w .z C / C W .z=ı/ : (2.50)
In this thesis, a power law approximation of Spalding's law, introduced in Section 2.4.2, is used to initialise the velocity field

f w .z C / D 8 < : U C 1 if z C > ı C ; 1 Ä ln z C C 5:5 if 90 < z C < ı C ; 32:22 z C 0:3147 C 24:52 if 4 < z C < 90 ; z C otherwise : (2.51)
A comparison between Equation 2.51 and Spalding's law is given in Figure 2.7. 2.51 in Section 2.5.1)), the law-of-the-wall of von [START_REF] Von Kármán | The analogy between fluid friction and heat transfer[END_REF] (black, dotted line) and of Spalding [1961] (blue, dashed line), and the DNS results of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] (circles). Equation 2.51 applies to the TBL over a smooth wall. Over a rough surface, protrusions of roughness elements alter the nature of the inner layer. In this thesis, a work has been done in order to propose a universal velocity profile for the rough-wall flows [START_REF] Huang | Law-of-the-wall in a boundary-layer over regularly distributed roughness elements[END_REF].

Details of this new law-of-the-wall are given in Appendix A.

From an extensive examination of experimental data, [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF] tabulated the value of W . Based on the data of Coles, Pope [2000] proposed the following approximation of

W W D ˘sin 2 z ı Á 2 : (2.52)
The constant ˘, the wake strength, can be determined using the values of the TBL thickness ı, the friction velocity u and the external velocity U 1 by applying Equation 2.50 at

z C D ı C .
The local value of the boundary-layer thickness of the base state velocity profile, ı.x/, is determined by the following formula given by [START_REF] Schlichting | Boundary-Layer Theory[END_REF]

U 1 ı D 0:21  Re x log Re x à ; (2.53)
with the Reynolds number defined as Re x D U 1 x= , between 10 5 and 10 6 . This condition is satisfied for the typical computational domain used in this work: L x U 1 = < 10 6 , L x being the streamwise extent of the domain.

The aforementioned initialisation procedure, used in all simulations presented in this thesis, has the sole aim of imposing stationary base state values to the time-dependent variables at the beginning of the computation. Additional treatments intending to accelerate the convergence of the simulation, e.g., by introducing fluctuations at the bottom wall [START_REF] De Villiers | The potential of large eddy simulation for the modelling of wall bounded flows[END_REF][START_REF] Duprat | Simulation numérique instationnaire des écoulements turbulent dans les diffuseur de centrales hydrauliques en vue de l'amélioration des performances[END_REF], are not used here, since the current implementation provides satisfactory results with reasonable efficiency.

Random fluctuation generation

In our simulations, the random fluctuation generation method is active at the initial phase of the calculation, since the domain is initiated using only the base state profiles. At each instant, the velocity fluctuations of the inlet velocity, u 0 i , are generated using the random number series j . Thus

u i D u r;i C u 0 i ; D u r;i C a ij j ; (2.54)
where the subscript r denotes the base state values. The random number series, j , are both derived from a Gaussian distribution with zero mean (h i i D 0) and unit variance ( i D 1). We note that the brackets denote a time average. These three sequences are statically independent from one another: h i j i D 0 for i ¤ j . Their contributions to the inflow velocity are regulated by the constant coefficients a ij of the matrix a.

As stated earlier, one necessary condition is that the generated inflow velocity verifies a prescribed Reynolds stress tensor, R ij . To this end, [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF] gave the following expression of a

a D 0 B B @ p R 11 0 0 R 21 =a 11 q R 22 a 2 21 0 R 31 =a 11 .R 32 a 21 a 31 /=a 22 q R 33 a 2 31 a 2 32 1 C C A (2.55)
Substituting Equation 2.55 into Equation 2.54 leads to hu 0

i u 0 j i D R ij .
In practice, inflow data are generated at each time step during the simulation, for each grid point located at the inlet, according to Equations 2.54 and 2.55. The prescribed Reynolds stress tensor R ij .z/ is determined from the classical profiles for a canonical boundary-layer (e.g. [START_REF] Schlichting | Boundary-Layer Theory[END_REF]) or specific experimental data.

In principle, the three sequences of random numbers should exhibit non-zero correlations, both temporally and spatially. This could be achieved through a fast Fourier transform, by imposing a spectrum of the fluctuations in the frequency space [START_REF] Le | Direct numerical simulation of turbulent flow over a backward-facing step[END_REF], which, unfortunately, adds complexity to the algorithm. In ARPS, a simpler procedure has been employed in order to generate temporally correlated random numbers at each grid point, according to the method of [START_REF] Deserno | How to generate exponentially correlated Gaussian random numbers[END_REF].

Consider a sequence of independent Gaussian random numbers with zero mean and unit variance g n , thus

prob.g n D x/ D 1 p 2 exp. x 2 =2/ 8n 2 N : (2.56)
Elements of g n are a priori uncorrelated with each other. From g n , we define a second sequence of random numbers n by

( 1 D g 1 ; nC1 D f n C p 1 f 2 g nC1 ; (2.57)
where f is a constant coefficient. It can be deduced that n is also a Gaussian sequence with zero mean and unit variance. More importantly, it can be shown that, for n , the autocorrelation coefficient c.nI m/ D f n . By defining f D exp. 1= 0 /, the autocorrelation function of n decays exponentially with a predefined correlation time 0 , which can be physically related to the Lagrangian integral time-scale from experimental results [START_REF] Swamy | Auto-correlation measurements and integral time scales in three-dimensional turbulent boundary layers[END_REF].

The randomly generated inflow data are intended to mimic the basic behaviour of a turbulent field, especially at the early stages of the simulation. In practice, after several flow-through times, meaningful flow characteristics begin to emerge inside the downstream flow field thanks to the random inflow generation at the inlet. From this moment, the development and sustainment of realistic turbulent structures is taken over by the extraction/rescaling technique, presented in the next section.

Inflow generation technique of Lund et al. [1998]

The inherent shortcoming of the random fluctuation inflow generation lies in the lack of realistic coherent structures. Since only the local and second-order statistics are satisfied for the velocity fluctuations, the resulting velocity field does not provide higher order statistics such as the skewness and the kurtosis. In such case, a long development domain is needed in order to remedy this lack of physical characteristics of the inflow.

Alternatively, [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF] proposed an approach of inflow generation for boundarylayers, in which the inflow informations are deduced from that located at a downstream station, far from the inlet. This makes sense since the downstream flow is resolved using the discretised NS equations, thus physically more realistic. Since then [START_REF] Tamura | LES study of turbulent boundary layer over a smooth and a rough 2d hill model[END_REF] has extended the method to a rough-wall case. Essentially, the inflow generation has two major procedures: a proper assessment and then an extrapolation of the downstream flow field to the inlet plane.

Before extraction, the velocity at the extraction plane is first decomposed into a mean and a fluctuating part. The decomposition is achieved by defining the mean as an average in the spanwise direction and in time. The velocity fluctuations are then defined as

u 0 D u u (2.58) w 0 D w w : (2.59)
The mean spanwise velocity, v, is assumed to be zero. In fact, subsequent simulation results show that the maximum of v barely reaches 1 % of U 1 , lending support to this assumption.

For the mean flow, the universal law-of-the-wall of the inner layer of the boundary-layer, as well as the velocity defect law in the outer layer, gives for the streamwise velocity

u inner .x; z C / D u .x/f 1 .z C /; U 1 u outer .x; Á/ D u .x/f 2 .Á/; (2.60)
where Á D z=ı is the outer layer non-dimensional coordinate and z C is the inner layer one. Functions f 1 and f 2 are assumed to be universal. The term u .x/ in Equation 2.60 accounts for the spatial development of the TBL in the streamwise direction.

We use the subscript "in" to denote the inlet flow and "ex" the extraction flow, as illustrated in Figure 2.8. From Equation 2.60 we have

u inner in .z C in / D u ex .z C in / ;
(2.61)

u outer in .Á in / D u ex .Á in / C .1 /U 1 ; (2.62)
where D u ;in =u ;ex . Using Equations 2.61 and 2.62, the flow field at a downstream position u ex is appropriately extracted and rescaled into the inlet field u in . The mean vertical velocity w, for its part, is assumed to scale as

w inner .z C / D U 1 g 1 .z C /; w outer .Á/ D U 1 g 2 .Á/;
(2.63) according to [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF]. Thus

w inner in .z C in / D w ex .z C in / ;
(2.64)

w outer in .Á in / D w ex .Á in / : (2.65)
v being assumed to be zero, no scaling is needed for the mean spanwise velocity.

The proper scaling of the velocity fluctuations is all the more crucial for the inflow genera- The exact form of the functions h 1 and h 2 is a priori unknown. Nevertheless, further approximation can be made on the streamwise homogeneity of these functions. In this case, the functions h 1 and h 2 are assumed to be periodic such that h 1;in h 1;ex and h 2;in h 2;ex . A similar expression is adopted for the spanwise fluctuation v 0 .

Finally, a composite velocity profile valid over the entire boundary-layer can be expressed using the weighted average of the inner and the outer profiles

u in D h u inner in C u 0 inner in i .1 W .Á in // C h u outer in C u 0 outer in i W .Á in / ; (2.72)
where the weighting function W is defined as

W .Á/ D 1 2 1 C tanh Ä a.Á b/ .1 2b/Á C b = tanh.˛/ : (2.73)
The constants are given by a D 4, b D 0:2 [START_REF] Lund | Generation of turbulent inflow data for spatiallydeveloping boundary layer simulations[END_REF]. A similar expression to Equation 2.72 can be written for v in and w in .

An illustration of the whole inflow generation procedure is given in Figure 2.9. Within one flow-through time, the downstream flow exhibits mainly a laminar behaviour, expressed by a weak velocity gradient near the wall (blue line in

Turbulent boundary-layer flow validation

In order to validate the performance of ARPS in simulating a spatially developing wind tunnel-scale boundary-layer, the following simulation case has been conducted and is presented here. The simulation is configured according to the experiments of Vinçont [1999] and Simoëns et al. [2007]. The characteristics of the boundary-layer are as follows [Simoëns et al., 2007;Vinçont, 1999] ı 0 is the value of ı prescribed at the inlet. The Reynolds number, defined as Re D ı 0 u = , equals to 500.

The computational domain is 18ı 1:7ı 2ı in the streamwise, spanwise and wall-normal directions, respectively. The mesh contains 1280 123 100 points. In wall units (using the wall shear evaluated at the inlet), the mesh resolution is x C D y C D z C wall D 7 and z C mean D 10, in accordance with the classical values used in the WMLES. The mesh is uniform in the streamwise and spanwise directions whereas a hyperbolic tangent stretching is applied in the wall-normal direction. The velocity field is initialised with the mean profile given by Equation 2.51. A wall model based on Spalding's law (Equation 2.43) is used in this simulation.

The extraction plane is located at 3ı downstream of the inlet. We assume that the extraction plane is located far enough from the inlet (i.e., at a distance exceeding the correlation length of the streamwise fluctuations) that artificial couplings are suppressed in the computed solution.

The total duration of the simulation is 17 flow-through times. The statistics are accumulated during the last 15 flow-through time in order to obtain statistically steady state data, which corresponds to 300 inertial time-scales ı=U 1 . Then, a spanwise average is applied to the data. . The simulation results are compared with the experiments of Simoëns et al. [2007] as well as the DNS results

of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] at Re D 660. The simulation results are closer to the experimental results at the same Reynolds number than to the DNS of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] at a higher Re . Overall, the results are satisfactory for the mean flow as well as for the RMS of velocity fluctuations.

Conclusions

In this chapter, the governing equations and the assumptions used in our study of a turbulent . Solid and dotted lines: DNS results of Spalart [1988] with Re D 660, for u 0 C and w 0 C , respectively. Red, LES. Blue, experiments. velocity fluctuations at the inlet at each time step. Lastly, the simulation results of a TBL on the mean flow and the RMS velocities are presented and compared to the experimental ones at the same Reynolds number, which show a good agreement between LES and the experiments.

Chapter 3

Physics and modelling of wind blow particles

In order to study the transport of solid particles in the boundary-layer, the equation that governs the motion of individual airborne particles is given in this chapter. The forces acting on the particles are identified and their relative importance is evaluated. By assuming that particles are small with high density, a simplified equation of motion is obtained and described in Section 3.1. Special attention is devoted to the initiation of particle motion, described in Section 3.2. Particles resting on the floor are subject to aerodynamic entrainment, which is related to the turbulent structures in the near wall region. The modelling of particle collisions with the floor, which is recurrent in the saltation process, is presented in Section 3.3. At last, several important parameters that determine the nature of aeolian particle transport such as the Stokes number and the gravity parameter are presented in Section 3.4. These parameters are used in the discussion of the results presented in the last two chapters.

Motion of airborne particles

In this section, the motion equation of individual particles is presented. The relevant forces acting on particles carried by the flow are identified. Since the movement of each particle is individually simulated, a number of hypothesis on the kinematics and dynamics of particles are used. These simplifications, based on the characteristic properties of the grains studied here, are a reasonable compromise to the computational cost. They are described in the following paragraphs.

Point-wise particle approximation

The transport of solid particles discussed in the context of this study belongs to a specific case of dispersed two-phase flows, in which a continuous fluid phase (wind) carries a discrete phase (solid particles). First, we assume that the solid particles are spherical and rigid. Moreover, the particle diameter D p is smaller than the characteristic length scale of the carrier flow L, and the solid particle density p is much larger than the fluid density . Thus we deal with small and heavy particles: D p =L 1 and p = 1.

By assuming that the particle diameter is sufficiently small so that the local Reynolds number tends to zero, [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF] derived the equation for the motion of a solid particle considered as "point-wise"

dX p dt D U p m p dU p dt D V p DU Dt C 1 2 V p d.U U p / dt " ƒ' … I C 3 2 D 2 p . / 1=2 Z t 0 1 p t  DU D dU p d à d " ƒ' … II C V p . p / g C F L C F D " ƒ' … III : (3.1)
Here, X p is the particle position vector and U p is the particle velocity vector. V p is the particle volume, m p is the particle mass, F L is the lift force and F D is the drag. In Equation 3.1, forces acting on the particle are separated into three subgroups. Forces regrouped in I include, firstly, forces due to the acceleration of the undisturbed surrounding flow and secondly, the added-mass force. The added-mass force arises from the displacement of a virtual fluid parcel at the position of the solid particle X p .t/, due to the acceleration of the solid particle relative to the fluid phase.

Term II is the Basset force, which is a historical effect of anterior accelerations weighted by the inverse of elapsed time. Term III contains the gravity, the buoyancy force and the aerodynamic forces. The forces in term III will be discussed in the next subsection.

From Equation 3.1, it can be seen that the 'added-mass' is only relevant when p .

Besides, the force due to fluid acceleration is proportional to the fluid density as well. Both forces in term I are thus negligible at p . The Basset force in term II is associated with the history of particle motion. As argued by [START_REF] Minier | The pdf approach to turbulent polydispersed two-phase flows[END_REF], the Basset force arises from the flow unsteadiness near the particle. For heavy particles with p , the particle relaxation time ( p , defined in Section 3.1.2) is much larger than the viscous diffusion time of the flow.

In this case, the local flow unsteadiness is not relevant and the Basset force is several orders of magnitude smaller than the drag force in term III [START_REF] Greeley | Wind as a geological process: on Earth, Mars, Venus and Titan[END_REF]. Therefore, here the Basset force is neglected.

In the following subsection we discuss the forces regrouped in term III of Equation 3.1, namely the gravity force, the lift force F L and the drag force F D .

Equation of motion

Since the particle is carried by the flow, the primary force involved is the aerodynamic force F due to both tangential (viscous) stress and normal stress (pressure). The value of F averaged over the surface of the particle is expressed as

F D Z :n d Z pn d ; (3.2)
where n is a unit normal vector, the viscous stress tensor and p the pressure.

Generally, F is decomposed into the drag force F D and the lift force F L . By definition, F D is the projection of F in the direction of the particle velocity relative to the flow, while F L is the other projection orthogonal to F D . These forces are illustrated in Figure 3. The lift force originates from the flow velocity difference between the upwind and downwind side of the grain. Due to the Bernoulli effect, a difference in velocity generates a pressure gradient that points towards the side with the smaller velocity. The pressure gradient gives rise to the lift. Due to the spherical shape of the solid particles studied in this thesis, the aerodynamic lift is mainly due to the inhomogeneity of the flow. This lift force is commonly named as the Saffmann force [START_REF] Saffman | The lift on a small sphere in a slow shear flow[END_REF]. Following [START_REF] Zheng | Mechanics of wind-blown sand movements[END_REF] F L can be expressed as

F L D 1 8 a D 2 p C L .u 2 up u 2 down / : (3.
3) Equation 3.3 predicts that the lift force is generated if the fluid is accelerated over the upper side of the grain: u up > u down . Due to the small size of the solid particles considered here compared with the characteristic length scale of the flow structures, the local velocity gradient of the flow is generally insufficient to generate appreciable lift.

Another factor that favours the generation of lift forces is related to the rotation of the particle.

Experiments have shown that particles in saltation predominantly have topspin [START_REF] White | Magnus effect in saltation[END_REF], which accelerates u up and attenuates u down . This lift force due to rotation is commonly named the Magnus force. However, [START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF] found that with the typical rotation velocity (300 revolutions per second reported in [START_REF] White | Magnus effect in saltation[END_REF]), the average Magnus lift force is one order of magnitude lower than the drag force F D or the gravity force. We thus neglect contributions of the lift force to airborne particles.

The drag force opposes the relative motion of the particle to the flow. Through dimensional analysis, the drag force is found to be proportional to a characteristic area of the particles, as well as the square of the particle velocity and the fluid density. Using the drag coefficient C D , the drag force is often expressed as

F D D 1 2 C D SjU r jU r ; (3.4)
where S D 1 6 D 2 p is the windward, projected area of the spheric particle. U r is the particle velocity relative to the fluid, expressed as where P D m p g. We note that the vertical and the horizontal motions of the particle are uncoupled in Equation 3.9.

U r D U p .t
Equation 3.9 can be modified into p is denoted as the relaxation, or the response time of the solid particle, and is widely used in the literature [START_REF] Dupont | Modeling saltation intermittency[END_REF][START_REF] Grigoriadis | Lagrangian particle dispersion in turbulent flow over a wall mounted obstacle[END_REF]Vinkovic et al., 2006a].

d U p dt D U r 0 p C g : (3.
Nevertheless, it is necessary to point out that the true particle relaxation time, 0 p , is dependent on the particle Reynolds number. According to Equation 3.8, 0 p is one order of magnitude lower than p for Re p > 400. This value of Re p can occur, albeit not often, in case of particle saltation. For example, after a collision with the wall, the rebounding particle changes abruptly its direction of motion, increasing a priori its relative velocity to the flow and decreasing 0 p . Thus, p represents an upper limit of the particle relaxation time under general flow conditions, since particles are always quicker to adapt their motion to the local flow with 0 p Ä p .

Two-way coupling

The two-way coupling is the influence of solid particles on the carrier fluid. Trailing wakes formed behind the solid particles act as a sink of fluid momentum. In order to account for the momentum transfer from grains to fluid, an additional drag force to the fluid momentum equation (Equation 2.18) is introduced. The drag force due to the presence of particles averaged over a grid cell is given by [Vinkovic, 2005]

f twoway D 1 V grid N p X pD1 m p U .X p .t/; t/ U p .t/ p f .Re p / ; (3.13)
where V grid is the grid cell volume and N p is the number of resolved particles within the cell.

An additional term is also introduced in the transport equation of the subgrid turbulent kinetic energy (Equation 2.25)

A u 00 i f 00 i D p ˚p 2k sgs p C T L f .Re p / ; (3.14)
where f 00 i is the fluctuation component of the force from particles to fluid and ˚p is the volume fraction in the grid cell occupied by the particles. T L is the Lagrangian correlation time scale of the flow, estimated by

T L D 4k sgs 3C 0 e " : (3.15)
C 0 being the Kolmogorov constant. ˚p is the volume fraction given by:

˚p D P N p pD1 V p V grid : (3.16)
Elghobashi [1994] claimed that starting from ˚p 10 6 , the two-way interaction between particles and the fluid can no longer be omitted. Inter-particle collisions become non-negligible for ˚p 10 3 and a 'four-way' coupling is thus necessary. For cases considered in this thesis, the mean volume fraction is of order 10 5 inside the TBL, while peak values of ˚p can reach values as high as 10 4 in the saltation layer, as will be shown in Table 5.4. The effect of two- way coupling is particularly important near the source of solid particles, close to the particle emission source in our simulations. In this area, a large number of slow-motion grains just lifted off remain close to the wall. Although the volume fraction ˚p could reach 10 3 near the surface, no 'four-way' coupling has been used in ARPS, which could be further implemented in the future.

Stochastic model of Lagrangian solid particle tracking

In the Lagrangian tracking of solid particles, one difficulty lies in the determination of the fluid velocity along the solid particle trajectories. The fluid velocity as 'seen' by the grain, U .X p .t/; t / (Equation 3.5), is neither Eulerian nor fluid Lagrangian, since both the position of the solid particle and the fluid element evolve with time. In the case of inertial particles, the trajectory deviation between the two phases is pronounced. Yet this velocity bears resemblance to the Lagrangian velocity of a fluid particle especially in the limiting case where light solid particles nearly follow the motion of the flow. In this section we present succinctly the subgrid stochastic model for fluid particles developed during the thesis of Vinkovic [2005] and Aguirre [2005], before introducing the stochastic model for Lagrangian solid particle tracking [Vinkovic et al., 2006a] which is more relevant to this work.

By construction, the Lagrangian velocity of a fluid particle can be expressed as

u L .t/ D e u.x p .t/; t / C u 00 L .t/ ; (3.17)
where e u is the large-eddy contribution to the fluid velocity at the position occupied by the solid particle. It can be obtained from an interpolation of the LES velocity field to the location of the "fluid parcel" occupied by the solid particle, Eddies smaller than the filter size are 'ironed out' in e u. u 00 L is the small-scale Lagrangian velocity fluctuation around e u and has to be modelled.

In Vinkovic et al. [2006a], the following Langevin-type, stochastic differential equation governing the subgrid-scale fluid particle velocity is used

du 00 L D  1 T L C 1 2k sgs dk sgs dt à u 00 L dt C s 4k sgs 3T L dÁ Á Á ; (3.18)
In ARPS k sgs is resolved by the 1.5 order transport equation introduced in Section 2.2.2. dÁ Á Á is the increment of a Wiener process vector with zero mean and delta-correlated in time hdÁ i dÁ j i D ı ij dt .

As mentioned at the beginning of the section, solid particles tend to deviate from the fluid parcel that originally contained them, mainly due to their inertia and gravity, creating a decorrelation between the two phases. To account for this fact, the fluid velocity seen by the solid particle is estimated by a modified version of Equation 3.18 in which T p L , a Lagrangian decorrelation time scale of the fluid velocity replaces T L :

T p L D T L ˛grav C ˛inert ; (3.19)
here ˛grav and ˛inert are coefficients related to gravity and inertia. For further details of the subgrid solid particle formulation, readers are referred to the theses of [START_REF] Aguirre | Dispersion et mélange atmosphérique euléro-lagrangien de particules fluides réactives[END_REF] and Vinkovic [2005].

Aerodynamic entrainment

In this section we discuss the aerodynamic entrainment of particles at rest on the bed. This process allows particles to change from a static state on the floor to a mobile state by a liftoff. Another type of particle entrainment is related to the splash process, caused by violent collisions of particles impacting on the sand bed, ejecting a number of grains into the air. This process can be interpreted as an homogenisation of mechanical energy between particles already in saltation with a momentum surplus and the ones immobilised on the wall. However, owing to the configuration of a "sandbox-type" particle feeding in the PC09 experiments, instead of the commonly used sand feed by downward injection, spurious splash movement at the sand bed is avoided by design. Thus the initiation of particle transport by splash is thus not considered in the simulations. Thanks to this configuration, particles are only removed from the bed if the aerodynamic forces exerted by the wind exceed a critical value that holds the grains to the floor.

Different approaches to model the aerodynamic entrainment are discussed in this section.

We introduce, successively, the steady-state approach of [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF] and [START_REF] Shields | Anwendung der Ähnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung[END_REF] (Section 3.2.1), the empirical take-off curve of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] based on the experimental results (Section 3.2.2) and the take-off model based on the balance of forces (Section 3.2.3). The section is concluded by the presentation of the take-off model implemented in the numerical simulations conducted in this thesis (Section 3.2.4).

Threshold friction velocity

The driving forces of the aerodynamic entrainment are the drag and lift forces. The zone of interest is in the vicinity of the wall. The shear stress of the flow is a useful parameter in this regard, since it is related to the strength of the flow momentum sink, created by the wall and the particles carried by the flow. In the steady-state saltation model of [START_REF] Owen | Saltation of uniform grains in air[END_REF], a threshold friction velocity is defined using the critical shear stress t for particle motion initiation

u t D p t = (3.20)
By definition, the threshold friction velocity u t is an important flow parameter in the study of wind erosion. Several saltation models take u t as the main input to their predictions on the saltation flux [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF][START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF][START_REF] Lettau | Experimental and micrometeorological field studies of dune migration[END_REF][START_REF] Owen | Saltation of uniform grains in air[END_REF]White, 1979].

Based on the initial work of [START_REF] Shields | Anwendung der Ähnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung[END_REF], [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF] derived a take-off criterion using the following expression of the threshold friction velocity:

u t D C p p gD p ; (3.21)
where the density ratio coefficient p D p = 1. The constant C is between 0.1 and 0.2 for Re p > 3:5 according to [START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF]. [START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF] proposed C D 0:08 0:09. Using a regression method on the transport rate measurements, the authors obtained the lowest friction speeds from which the rate increases steeply for sand grains.

The squared ratio between u and u t , given by Equation 3.21, is proportional to the socalled Shields parameter

Sh D u 2 p gD p : (3.22)
Sh is a measure of the transport capacity of solid particles by the flow. Although results from laboratory and in-situ measurements have shown more than an order of magnitude variability in the value of this criterion [START_REF] Diplas | Coherent Flow Structures, Initiation of Motion, Sediment Transport and Morphological Feedbacks in Rivers[END_REF], the use of Equation 3.21 and the Shields number (Equation 3.22) remain the standard approach in the global models of wind erosion.

The take-off curve of Foucaut & Stanislas [1996]

Part of the discrepancies observed in the use of Bagnold's formula (Equation 3.21) for the threshold friction velocity (Equation 3.21) may be accounted for by the uncertainties related to the in-situ measurements and the resulting value of C . In order to obtain a finer characterization of this parameter, Foucaut & Stanislas [1996] conducted direct measurements of u t in a wind tunnel. Special care was taken by increasing the wind tunnel velocity gradually and linearly in order to obtain an accurate detection of particle emission.

In the subsequent analysis, [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] introduced an alternative set of nondimensionalized parameters to Sh and

D C p D D p u = , 8 < : u ? t D u u ref D u .g p / 1=3 D ? p D D p D ref D D p g p 2 Á 1=3 ; (3.23)
The main advantage of Foucaut's formulation lies in the definition of the non-dimensional parameters D ? p and u ? t . In fact, we deduce from D ref and u ref a Reynolds number

Re ref D D ref u ref D 1 ;
which states an equilibrium between the aerodynamic and viscous effects acting on the solid particle. Also we have

u 2 ref g p D ref D 1 ;
that relates to a state of balance between the aerodynamic effect and gravity. Thus, by separating the particle diameter and the friction velocity in the definitions of u ? t and D ? p , [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] argue that the relationship between u ? t and D ? p is more likely to be universal than the one between Sh and D C p .

Based on experimental results, the following expression of u ? t is given by Foucaut & Stanislas [1996]: This criterion can be used to evaluate, a priori, the erodibility conditions of the surface for a given set of flow conditions and particle characteristics. Nevertheless, Equation 3.24 remains empirical. Here, another modelling approach is used that takes into account the instantaneity of the flow as well as the lift force acting on the solid particle. In the case of Lagrangian particle tracking, this is physically more meaningful in our opinion.

u ? t D 22:71 D ?

Emission model based on the balance of forces

In the emission model proposed by [START_REF] Descamps | Érosion éolienne d'un lit de particules à large spectre granulométrique[END_REF], the aerodynamic entrainment is based on the balance of the forces exerted on a particle at the surface:

F L C F 0 L F adh C F 0 adh C P ; (3.25)
where the gravity force of the particle is P D m p g. The lift F L as well as the adhesion force F adh are considered as the sum of a mean value plus a fluctuating one, in order to account for the turbulent nature of the flow. Based on the experimental measurements of [START_REF] Mollinger | Measurement of the lift force on a particle fixed to the wall in the viscous sublayer of a fully developed turbulent boundary layer[END_REF], the mean lift force writes

F L D 15:5 2  u D p à 1:87 : (3.26)
The mean adhesion force, following [START_REF] Zimon | Adhesion of dust and powder[END_REF], is evaluated as

F adh D c 1 D p (3.27)
where c 1 D 1:43 10 5 N m 1 . Note that this value could be highly controversial and very dependent on the nature of the grain. Nevertheless, there has been evidence that for particles of diameter larger than 100 ţm with negligible humidity, gravity dominates over adhesion [START_REF] Zimon | Adhesion of dust and powder[END_REF].

The model of Descamps represent turbulent events such as ejections with predetermined statistical relationships. The criterion for take-off, Equation 3.25, is tested locally at each time step over a surface covered by solid particles. The most important drawback of the model developed by [START_REF] Descamps | Érosion éolienne d'un lit de particules à large spectre granulométrique[END_REF] is related to the evaluation of the lift: the determination of its value is based on statistical relationships decoupled from the instantaneous flow parameters such as the local fluid velocity and wall shear stress that inherit distinct characteristics from the turbulent flow.

Dupont [2012] (personal communication) proposed an instantaneous evaluation of the aerodynamic lift as a function of the local surface wind speed components

F L D c 2 D 2 p  1 Ä ln  30z 1 D p Ãà 2 C c 3 ! 1 .u 2 C v 2 / ; (3.28)
The value of c 2 D 47 is adjusted according to the take-off curve of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF].

The value of the constant c 3 D 10:6 is fixed to u D 2:7u and v D 1:8u at z 1 . By further assuming that the local wind speed follows a logarithmic profile, the mean lift force equals to

F L D c 2 D 2 p u 2 ; (3.29)
We note that the value of F L given by Equation 3.28 is 3 to 4 times larger than the one given by Equation 3.26. This underlines the high importance of correctly setting the coefficients in the take-off models.

By averaging Equation 3.25, we obtain the following mean-force balance for a particle resting on the bed

F L F adh C P (3.30)
Using Equations 3.29 and 3.30, the threshold friction velocity writes

u t D s 1 c 2 D 2 p 6 . p /gD 3 p C c 1 D p Á ; (3.31)
For a typical sand grain with p D 2500 kg m 3 , Equation 3.31 gives comparable values of u t to the take-off curve of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF].

Both Equations 3.24 and 3.31 provide an amplitude criterion for the threshold conditions of aerodynamic entrainment. However, these criteria are not sufficient given the objective of this thesis. The threshold friction velocity u t , obtained by these criteria, is only a time-averaged measure of the flow intensity, while the take-off process is an instantaneous process that may be influenced by extreme turbulent events.

Take-off model implemented in ARPS

In this section, a new take-off criterion is introduced, based on the instantaneous evaluation of the lift force and the impulse of the net forces exerted on the particle.

White [1940] claimed that it is the peak values, instead of the mean ones, of flow velocity and aerodynamic forces that are predominantly responsible for particle entrainment. This claim has been supported by experimental results showing that the turbulence level together with the Shields number enhance greatly sediment transport in water [START_REF] Schmeeckle | Forces on stationary particles in nearbed turbulent flows[END_REF][START_REF] Sumer | Influence of turbulence on bed load sediment transport[END_REF]. The role of turbulent structures on the initiation of particle motion is thus nonnegligible in the TBL. The use of a time averaged approach and of the threshold friction velocity is thus physically less meaningful for the cases studied here.

In the temporal signal of flow velocity, the turbulent events are manifested by rapid departures from the local mean conditions. For example, flow zones in the near wall region exhibit a pattern of low and high velocity streaks elongated in the streamwise direction, as illustrated in Figure 3.2. One parameter that is closely related to turbulent structures in particular is the instantaneous value of the Reynolds stress , defined in Equation 2.35. This link is highlighted by expressing as a sum of four conditionally averaged terms through the so-called quadrant analysis,

u 0 w 0 D u 0 ˚w0 " ƒ' … Q1 C u 0 « w 0 " ƒ' … Q2 C u 0 « w 0 « " ƒ' … Q3 C u 0 ˚w0 « " ƒ' … Q4 (3.32)
The terms from Q1 to Q4, each associated with a distinct type of flow structure, are conditionally averaged according to the indicated signs. Since in wall-bounded flows, the value of u 0 w 0 is largely negative, Equation 3.32 shows that the main contributions to the Reynolds stress result from the Q2 and Q4 events (Figure 3.3). At the same time, the physical interpretation of Q2 and Q4 is more straightforward than the other two, since Q4 events are generally related to high-speed fluid that swipes the surface from the upper part of the TBL, and Q2 events are generally related to the ejection of low momentum fluid moving away from the wall. The strongest among these turbulent events, commonly designated as gusts, contribute greatly to the instantaneous value of the Reynolds stress and enhance momentum exchange across the TBL. We can thus assume that these events have an impact on the solid particle entrainment from the surface. There is much experimental evidence in the field of sediment transport in rivers on the presence of Q4 events near or at the time of grain entrainment and a consistently high, positive drag force fluctuations during these events [START_REF] Dwivedi | Flow structures and hydrodynamic force during sediment entrainment[END_REF][START_REF] Hofland | Measurement of Fluctuating Pressures on Coarse Bed Material[END_REF]. These fluctuations enhance the incipient motion of solid particles. The dependence of particle lift force on the turbulent events is less clear. While it seems logical that a downward sweep of fluid parcel decreases the lift, [START_REF] Dwivedi | Hydrodynamic Forces Generated on a Spherical Sediment Particle during Entrainment[END_REF] identified high pressures beneath the grain and positive fluctuations in drag and lift during Q4 sweeps. On the other hand, in situ measurements of aeolian also revealed the substantial role played by Q4 events. [START_REF] Wiggs | Turbulent flow structures and aeolian sediment transport over a barchan sand dune[END_REF] found that turbulent structures with positive streamwise fluctuations (Q1 and Q4) have a larger influence on sand transport on dunes. Among the two, the Q4 sweep events are more representative both in terms of intensity and frequency, attested by the generally negative sign of the Reynolds stress. Quantitatively, it is reported in [START_REF] Wiggs | Turbulent flow structures and aeolian sediment transport over a barchan sand dune[END_REF] that up to 83 % and 95 % of transporting events at the toe and crest, respectively, can be accounted for by sweeps, based on high-frequency in-situ measurement data over dunes.

Here, we attempt to introduce an evaluation of the particle lift force that takes into account the non-stationary aspect of the flow while at the same time conforms with the experimental result of [START_REF] Mollinger | Measurement of the lift force on a particle fixed to the wall in the viscous sublayer of a fully developed turbulent boundary layer[END_REF] on the mean lift. The model is adapted in order to evaluate the instantaneous lift force and to test the particle entrainment in a local and instantaneous manner. This is made possible by LES since the large scales of the turbulent flow are resolved at each time step.

Firstly, we propose the following relationship by assuming the proportionality between the instantaneous value of the lift F L .X p ; t/ and the mean lift F

L F L .X p ; t / F L D u 0 ˚w0 « .X p ; t/ hu 0 ˚w0 « i y ; (3.33)
the bracket hi y denotes a spatial average along the transverse direction. The average lift force F L is based on Equation 3.26 derived from the experiments of [START_REF] Mollinger | Measurement of the lift force on a particle fixed to the wall in the viscous sublayer of a fully developed turbulent boundary layer[END_REF].

The velocity fluctuations are evaluated in relationship to their spatial average along the spanwise direction in wall-adjacent cells. Using Equation 3.33, we assume the predominant role of Q4 events in the dislodgement of stationary particles.

The main conclusion from the experimental results of [START_REF] Diplas | The Role of Impulse on the Initiation of Particle Movement Under Turbulent Flow Conditions[END_REF] is that not only the magnitude, but also the duration of energetic near-bed turbulent events is relevant in predicting grain entrainment. Thus, the product of force and its duration, or impulse, is a more appropriate parameter for identifying flow conditions favouring the particle take-off. This point was not considered in the classical emission models such as [START_REF] Descamps | Érosion éolienne d'un lit de particules à large spectre granulométrique[END_REF].

Based on these arguments, the following entrainment criterion is used in ARPS:

For w p .t 1 / D 0, IF 8 < : F L .t / > F adh C P for t 2 OEt 1 ; t 2 And t 2 X t 1 .F L F adh P / t > m p w 0 Then w p .t 2 / D t 2 X t 1 .F L F adh P / t=m p
In practice, we impose further t 2 t 1 < p .

The threshold lift-off velocity is w 0 D p 2gD p . This value corresponds to the amount of potential energy that needs to be overcome for a resting particle to be lifted to a height of one diameter over the surface. In this way, we only consider a take-off event as valid if the sand grain is effectively removed from the ground.

Particle-bed interaction

Statistical relationships on particle rebound

Submitted to gravity, all grains that are ejected into the saltation layer ultimately fall back to the ground and a collision with the surface follows. The consequence of the collision is that during a relatively short time period, an exchange of momentum and kinetic energy takes place between the incoming particle and the floor. Whether the particle rebounds or remains on the floor depends on the characteristics of the particle before impact, as well as the properties of the floor, such as the humidity, the packing density of the wall roughness, etc.

This collision process has to be modelled due to the complexity of the phenomenon. The multiple factors involved can be classified into three categories: the particle (morphology, size, material, etc.), the floor (humidity, components, etc.) and the flow (wind intensity, turbulence, etc.). All these factors are likely to interact with each other during a very short period of time due to the nature of the collision.

Experiments conducted in artificial conditions aim to study the details of the collision process, by propelling solid particles into a static bed of similar particles [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF][START_REF] Mitha | The grain-bed impact process in aeolian saltation[END_REF][START_REF] Rioual | Experimental study of the collision process of a grain on a two-dimensional granular bed[END_REF]Werner, 1990]. However, due to the complexity of the microtopography, the particle-bed collisions taking place in natural conditions are generally considered as a stochastic process. Particle-bed collisions are characterized by the size of the impacting grain, its angle and speed, and the nature of the local bed (grain size distribution, angle with respect to horizontal). These parameters are generally characterised from a statistical point of view.

In order to simplify the modelling, we assume here that the collision are characterized by the prediction of the rebound velocity V r based on the value of the impact velocity V i . As mentioned above, due to the complexity of the problem, the collision process is often considered as a stochastic one [START_REF] Anderson | Wind modification and bed response during saltation of sand in air[END_REF]. Based on the simulation results of [START_REF] Anderson | Wind modification and bed response during saltation of sand in air[END_REF], the probability that a particle in saltation rebounds upon impact can be approximated by P r D 0:95 .1 exp. ˛vi // ; (3.34) where v i is the magnitude of the impact velocity and ˛is an empirical parameter of order 2 s m -1 . Note that Equation 3.34 is a statistical expression deduced from numerical simulation results.

Particles that fail to rebound from the wall, as predicted by Equation 3.34, are not necessarily immobilized at the floor. Their motion may be better described as "traction" or "creep" (Section 1. 1.3). As they roll and slide along the wall, they can more readily return to the flow and resume saltation than particles immobilised on the ground. This type of particle motion is judged to be secondary to particle saltation and thus not considered in this thesis. Instead, particles losing much of their kinetic energy due to collision are immobilised temporarily at the floor.

Ceasing to be airborne, these particles are subsequently subject to the aerodynamic entrainment by the take-off model presented in Section 3.2.4.

The rebound velocity is often characterized by the mean rebound angle  r and its norm v r . For a typical particle-bed collision, since the incoming particle has been preferentially accelerated by the flow in the streamwise direction, the collision often takes place with a small grazing angle. Experimental values for the impact angle range from 5 ı to 15 ı , while the rebound angle is usually greater than 20 ı [START_REF] Nalpanis | Saltating particles over flat beds[END_REF][START_REF] White | Magnus effect in saltation[END_REF]. In the numerical simulation of [START_REF] Anderson | Wind modification and bed response during saltation of sand in air[END_REF], the mean rebound angle is found to be 35°to 45°and the rebound speed 50 % to 60 % of the impact velocity, while the impact angles and velocities are fixed a priori. In LES, Vinkovic et al. [2006a] used the following values for  r and v r :

v r D 0:3v i ˙0:25v i ; Â r D 30 ı ˙15 ı : (3.35)
These two stochastic variables are assumed to follow a uniform probability distribution. Coefficients used in Equation 3.35 are mainly based on the experimental results of [START_REF] Nalpanis | Saltating particles over flat beds[END_REF].

An alternative approach for modelling particle-wall collision is by parametrizing the energy restitution coefficients, defined as the ratio of the particle/bed velocity after and before the collision. In the numerical model of [START_REF] Kok | A comprehensive numerical model of steady state saltation (COMSALT)[END_REF], the fraction of kinetic energy retained by the rebounding particle is described by a Gaussian distribution and the rebound angle by an exponential distribution. The kinetic energy of the rebounding particles is 45 ˙22 % of the impacting kinetic energy. The rebound angle is given by an exponential distribution with a mean of 40°from horizontal.

Rebound model based on the experimental results of Beladjine et al. [2007]

When aeolian saltation is established over a flat surface, the variation of  i is rather limited as the process reaches a steady state. In fact, the mean rebounding particle velocity is assumed to be constant in the models described in Section 3.3.1. They are well suited for this scenario due to their simplicity. Yet, over a hilly terrain, collision scenarios multiply due to topography. Bed slope can both enhance collisions with certain impact angles and inhibit others. In this thesis, the effect of local bed slope is taken into account in the modelling of particle-wall collision.

In [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF], a collision experiment between an incident bead and a threedimensional granular packing made of particles identical to the impacting one was carried out.

The impact angle and the impacting speed were varied in a wide range: 10 ı < Â i < 90 ı and 50 < v i = p gD p < 200. Using a fast video camera, the authors captured and analysed the trajectories of the rebounding and ejected particles projected onto the incident plane. Although more attention was given to the velocity distribution of the ejected particles due to the splash, the rebound angle and the projected rebound velocity in the impact plane were also reported.

One drawback of applying the experimental results of [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF] to the simulation of particle saltation is that the experiment was conducted without a mean flow. Nevertheless, during the extremely short period of time and in the vicinity of the wall where collisions take place, we assume that the local flow has minor influence on the outcome of the collision.

One of the main findings of [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF] is that the mean restitution coefficients decrease with increasing  i . Two restitution coefficients, e z and e xz , are defined, respectively, in the wall-normal direction and inside the horizontal plane. That is, particles lose more energy while colliding at angles normal to the ground than at small grazing angles. The dependence of other parameters such as the impact velocity on e z and e xz was shown to be secondary. The authors proposed the following relationships for e z and e xz as a function of  i :

e z D A z sin  i B z e xz D A B sin  i (3.36)
with A z D 0:30, B z D 0:15, A D 0:87, and B D 0:72.

From Equation 3.36, we propose the following expression of  r (Figure 3.4):

 r D asin A z B z sin  i A B sin  i ! C ˛0 ; (3.37)
where ˛0 is the local bed slope. This parameter is not taken into account in the experiments of [START_REF] Beladjine | Collision process between an incident bead and a three-dimensional granular packing[END_REF] since it was conducted over a flat bed. In our work, the effect topography to particle rebound is accounted for by ˛0 to the first order. In Figure 3.4, Â xy is the horizontal deviation of the rebounding particle inside the plane (x; y). In ARPS, Â xy is assumed to follow a predefined Gaussian distribution with zero mean and of standard deviation around 10°according to values from literature [START_REF] Dupont | Modeling saltation intermittency[END_REF][START_REF] Kok | A comprehensive numerical model of steady state saltation (COMSALT)[END_REF].

Using Equation 3.37, two parameters among e xz , e z ,  r and v r are enough for modelling the impact/rebound process. In ARPS, e xz and  r are adopted and their values are assumed to follow a Gaussian distribution with the mean value determined by Equations 3.36 and 3.37. The following standard deviations are chosen:

 r D  r e xz D 2  r  r Á e xz (3.38)
The idea behind Equation 3.38 is that particles with higher rebound angles are more likely to lose more kinetic energy. This is supported by the observation of Anderson et al. [1991] based on experimental results, that when the impact angle increases, the rebound angle increases and the ratio of rebound speed to incident speed decreases.

Modes of particle transport and non-dimensional parameters

As discussed in Section 3.1.4, a heavy particle tends to deviate from the motion of the carrier fluid mainly due to the inertial and the gravitational settling effects. These two influence particle motion, or the 'trajectory-crossing' between the solid and the fluid particles as illustrated in

Inertia and relaxation time

The relaxation time p , defined in Equation 3.12, is commonly considered as the characteristic time scale of the solid particle carried by a fluid phase. Higher is p , better the solid particle maintains its proper motion due to its inertia. The relative importance of p can be characterized by the dimensionless Stokes number St, defined as

St D p = f ;
where f is a characteristic time-scale of the turbulence. Using the large eddy turn-over time inside the TBL, T L D ı=u , we define

St L D p =T L ;
as a Stokes number related to the large turbulent structures. If St L 1, even the large eddies can hardly influence the trajectories of solid particles. Their motion is in general ballistic. Alternatively, using T Á D p ı =u 3 , the Kolmogorov time scale characterising the smallest turbulent structures of the TBL, we define

St Á D p =T Á If St Á 1
, particles are responsive even to smallest eddies. In this case, their motion is mainly controlled by the carrier flow and its turbulence.

Gravitational settling and terminal velocity

The effect of gravity on particle motion is characterized by another parameter, the terminal velocity. By definition, the particle terminal velocity, w t , corresponds to the relative velocity u r at which the particle experiences zero acceleration, i.e., du p =dt D 0. In order to derive its expression, we rewrite Equation 3.10 into

du p;i dt D u r;i 0 p ı i3 g : (3.39)
We deduce from Equation 3.39 that the horizontal component of the terminal velocity is zero.

Thus w t D g 0 p :

(3.40)

w t provides an asymptotic value of particle velocity transported by the flow. For small values of Re p (Stokes regime), we have w t D g p . Now whether an airborne particle can remain suspended in the air depends on the sign and the amplitude of the terminal velocity relative to the ground. If the vertical component of fluid velocity is noted as w f , the absolute terminal velocity of the particle is

w p D w f C w t : (3.41)
The absolute terminal velocity is thus determined by p and w t . The discrepancy between the behaviours of different types of solid particles can be explained by Equation 3.41. Due to a relatively small size, dust particles usually have a small response time p and jw t j is closer to zero, therefore they tend to remain suspended in air (w p 0). In contrast, sand particles have a much larger response time and thus tend to fall back to the surface very quickly, participating probably to the saltation process.

The relative effect of gravity to the intensity of the carrier eddies is measured by the 'crossingtrajectories' parameter [START_REF] Yudine | Physical Considerations on Heavy-Particle Diffusion[END_REF]. In fact, if the turbulent eddy is sufficiently strong, the movement of the solid particles is continuously sustained by the flow in suspension. The net effect of gravity, which constantly drives all airborne particles in saltation to the floor, is thus suppressed. The crossing-trajectories parameter is defined as [START_REF] Taniere | On the behaviour of solid particles in a horizontal boundary layer with turbulence and saltation effects[END_REF] g D w t =u :

g is also denoted as the gravity parameter [Vinkovic et al., 2011]. Here u is considered as the relevant velocity scale of turbulent motions in the vertical direction.

For the mode of airborne particle motion, the value of the Stokes number generally deter-

Conclusions

In this chapter, the governing equation of particle motion relative to the carrier flow is introduced. The particle drag and gravity forces are considered as the dominant ones, by accounting for the small particle size relative to the length scales of the flow, and the high density compared to the fluid. The aerodynamic entrainment of solid particles on the wall is related to the turbulent structures in the near wall region by the take-off model based on a force balance approach. A model of particle rebound at the wall is derived from the experimental results on particle-wall collision. The rebound angle and the restitution coefficient are modelled stochastically as a function of the characteristics of the impacting particle. Lastly, several non-dimensional parameters, such as the Stokes number and the gravity parameter, that determine the type of particle motion inside the ABL are presented.

Chapter 4

LES of turbulent boundary-layer flow over 2D hills

As mentioned in the Introduction, the recirculation zone (RZ) formed due to flow separation behind a 2D hill is the key aerodynamic feature studied in this thesis. It results from the interaction between the flow and the topography. The flow inside the RZ presents generally distinct characteristics compared with the outer flow. These characteristics are of special interest to the study of particle transport over obstacles, since air parcels as well as particles carried by the flow spend on average a longer residence time there. This leads to an increase of, e.g., chemical reactions and particle deposition. An example of the hill-induced flow separation is given in Figure 4.1, based on the LES results obtained by ARPS, previously described in Chapter 2.

The image shows notably the appearance of pockets of backflow formed behind the hill and structures of intense vorticity emanating from the hill crest. This chapter is structured as follows. A validation by comparison with a reference experimental case [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], denoted as CT06, of a turbulent flow over a steep sinusoidal hill is presented in Section 4.1. Then, the experimental campaign (PC09) to which this thesis is related is detailed in Section 4.2. Further numerical tests are first conducted over a single, Gaussian hill. Then, the simulation results with successive Gaussian hills according to the PC09 configuration are presented and discussed. In Section 4.3, flow separation behind 2D hills with different shapes is studied. Using LES, the critical angle of flow separation, the dependence of the RZ characteristics on the hill slope, the Reynolds number and wall conditions are investigated and discussed.

Aerodynamic case CT06

Description

In [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], wind tunnel experiments using split-fibre and cross-wire probes have been carried out in order to study a turbulent flow over a 2D steep hill. The experiment showed that a steady, relatively large separation occurs on the lee side that modifies the downstream flow characteristics. Small cubes are placed on the hill surface in order to further study the effects of wall roughness on the RZ. We note that this study corresponds one of the rare experimental cases with a 2D transversal hill, carried out in both smooth-and rough-wall conditions, for which detailed comparisons on second-order flow statistics with numerical simulations exist in literature [START_REF] Cao | Numerical study on turbulent boundary layers over two-dimensional hills -Effects of surface roughness and slope[END_REF][START_REF] Tamura | LES study of turbulent boundary layer over a smooth and a rough 2d hill model[END_REF]. In this section, we validate the LES by comparing the computed flow characteristics with the experimental data, in both smooth-and rough-wall conditions.

The geometry of the sinusoidal hill is sinusoidal and defined as

h.x/ D 8 < : H cos 2 x 2L Á if jxj < L ; 0 otherwise : (4.1)
The height of the hill, H, is 40 mm and the length L, which corresponds to the streamwise extent of the hill, is 100 mm. The maximal slope on the lee side is 32°. A rough wall covered with 5 mm-high cubes exhibits a roughness length z 0 of 0:2 mm. The characteristics of the flow are given in Table 4.1. The depth of the undisturbed TBL at the location of the sinusoidal hill, ı 0 , gives a nominal blockage ratio H=ı 0 of 0:16. At this ratio, the effect of blockage was assumed negligible by [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF]. Two Reynolds numbers are defined here: Re D ı 0 u = , also noted as ı

Validation of mean flow and turbulence characteristics

The CT06 case is simulated by ARPS. Since the height of the sinusoidal hill and the TBL thickness are relatively large, grid spacings can be optimised, decreasing computational cost.

Numerical details of the corresponding simulations are given in Table 4.2. The grid for the smooth case is more refined near the wall compared to the rough case in the vertical direction, in order to tackle the higher flow velocities in this region. The surface flux wall model of ARPS, described in Section 2.4.2 is used in the simulations. For the rough-wall case, the roughness effect is taken into account by the use of the roughness length z 0 .

Case In what follows, profiles of mean velocity and RMS of velocity fluctuations are presented in In Figure 4.2, the mean streamwise velocity profiles around the sinusoidal hill are presented.

N x N y N z L x =H L y =H L z =H
Coordinates are scaled by the hill height. The origin of the longitudinal coordinate is at the centre of the hill. The velocity is normalised by the free-stream velocity U 1 . For each velocity profile, U=U 1 D 0 is indicated by a dashed line passing through the corresponding X=H, from which the profile is taken. The simulation results (full line) are compared to the experimental data (points). Globally, the mean velocity profiles on the lee side of the hill differ greatly from those on the upstream side mainly due to the presence of the RZ. For both smooth-and rough-wall conditions, similar flow behaviours around the 2D hill can be identified: the flow decelerates on arriving at the upwind hill foot before accelerating gradually on the windward side, culminating in the remarkable velocity speed-up on the hill top. Although the flow reattaches slightly after 5H (smooth wall) or 6H (rough wall) after the hill, a longer distance is needed for the flow to recover its upstream TBL state before the hill. This is due to the development of an internal boundary layer that starts as soon as the flow reattaches. Globally, the agreement between the LES and the experiments is quite good, especially for the smooth-wall case. For the rough-wall case, although the wall roughness is only accounted for by a roughness length in the simulation, the mean flow results are very close to the experimental ones. Behind the hill, the increase in the streamwise RMS is remarkable. This is mainly due to the separation of the shear flow shortly after the hill crest. The peak of the RMS in this region reaches the double of the windward value. The location of the peak is between 3H and 4H for both smooth-and rough-wall conditions. On the smooth wall, the persistence of a wake region with higher intensities along the downstream zone hints at a quicker reattachment than on the rough wall. An active mixing due to this wake region allows the flow to reattach more rapidly. This is supported by both LES and experiments (Figure 4.2).

U/U ∞ 1 1 1 1 1 1 1 1 1 Z / H X / H (b) Rough wall
It is shown here that in the CT06 case, the LES performs fairly well, for both smooth-and rough-wall conditions. The main characteristics of the turbulent flow behind the hill are duly reproduced, in the upper shear layer and inside the recirculation zone. The main discrepancies lie in the wake region behind the rough hill, where the LES overestimates the turbulence strength.

This is possibly due to the lack of explicit simulation of cubic roughness, which, according to [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], apparently suppress the turbulence level behind the hill. Nevertheless, these discrepancies remain small. Besides, measurements of turbulence levels could be underestimated in this region due to both the presence of backflows and the spatial resolution of the measurements.

Here, the code has been validated for a single-hill case. In the next section, a case with multiple hills is presented. This experimental case has the advantage of containing solid particles and, to the best of our knowledge, represents the first experimental results on the solid particle transport over 2D hills.

Successive-Gaussian-hill case PC09

In this section, the experimental campaign (PC09) and the relevant LES results are presented.

Aerodynamic aspects of the experimental configuration are given in Section 4.2.1. In order to find the optimum numerical configuration regarding in particular the reproduction of RZ in the PC09 case, a single-hill configuration is devised. Different grids and wall models are tested and the main results are presented in Section 4.2.2. Simulation results on the successive-hill case are discussed in Section 4.2.3.

Experimental configuration

In the PC09 configuration, flow regimes were chosen as a compromise between the aerodynamic entrainment of upstream particles and the appearance of a notable recirculation zone on the lee side or between the hills. Three regimes are studied, each with different free-stream velocity and friction velocity. The set-up is depicted in At the beginning of the experimental campaign, the following Gaussian hill shape was used

Mnemonic ı 0 (cm) U 1 (m s -1 ) Re Re H H (cm) z 0 (µm) z C 0 V30 7
("AC case") h.x/ D H exp  0:35 x H Á 2 à ; (4.2)
where the height is H D 10 mm. However, this hill geometry did not generate a noticeable recirculation zone, necessary for the study of a potential solid particle trapping, on the lee side at large flow regimes.

The second 2D Gaussian hill type, hereafter referred to as "NC case", was then devised, it has the following geometry

h.x/ D H exp  0:7 x H Á 2 à : (4.
3)

The maximal slope on the lee side is 31°, which is nearly 2 times steeper than the AC case. The hill height gives a blockage ratio of 1=7 at the hill top. With a model/field scale of 1:10 000, the experimental configuration is consistent with an ABL of 150 m and an equivalent dune height of 20 m.

In the following sections, the objective is to simulate the PC09 cases with the NC hill geometry to study the flow characteristics and particle transport. Special attention is given to the V90 case with the highest wind speed, which is susceptible of transporting substantial numbers of solid particles. The NC hill geometry, with a steep slope close to that of the CT06 sinusoidal hill, could generate large RZ on the lee side, which increases the potential of solid particle trapping in this region.

Numerical tests on grid and wall model 4.2.2.1 The single-hill numerical configuration

In this section, the numerical details of the tests on the choice of the grid and the wall model are given. We point out the extra difficulties in the numerical implementation of the PC09 case compared to the CT06 case, due to the smaller boundary-layer depth and the lower hill. The dimensions of the experimental set-up impose a stringent requirement on grid sizes used in the numerical study. The vertical grid spacing z, proportional to the hill height, is four times finer than the CT06 case. The lateral grid spacing, limited by the narrow width of the steep hill (around 3H for NC), is five times smaller than in the CT06 case. In this case, the time step is also reduced in order to ensure numerical stability. One criterion is given by the Courant-Friedrichs-Lewy (CFL) condition, which states that C D u dt=dx Ä 1, indicating that the time step should be shorter than the time that a fluid particle takes to pass through one grid cell. A proportionality is thus established between the time step and the smallest mesh size. For the simulations presented here, the CFL number is close to 0:3 for the minimum grid cell near the hill crest for the highest flow regime. This results in a reduction of the time step by a factor of 6 compared to the CT06 case. All these factors lead to a substantial increase of the computational cost compared to the CT06 case.

Three grids with different mesh sizes have been tested, noted as A, B and C. For each grid, the streamwise and the spanwise grid cells are of equal size, whereas the vertical grid spacing is more refined near the surface aiming a better resolution of the near-wall flow. From A, B to C, the grid are gradually refined near the wall. Cells of the coarsest grid A are cubic in the near-wall region with a dimension a tenth of the hill height. For B, the resolution of the first computational point from the wall is doubled compared to A, with an adjustment of the grid aspect ratio. The total number of grid points, the streamwise and the spanwise spacings remain the same between A and B. In the finest grid C, both the resolution and the number of grid points are doubled in all the directions compared to grid B. This has been arranged in order to maintain both the same grid aspect ratio and domain dimensions between B and C. The domain size and mesh resolution of each case are summarized in Table 4.4. size after a 2D hill reported in the literature [START_REF] Almeida | Wake flows behind two-dimensional model hills[END_REF][START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF][START_REF] Kanda | Particle-Image Velocimetry Measurements of Separation and Re-attachment of Airflow over Two-Dimensional Hills with Various Slope Angles and Approach-Flow Characteristics[END_REF]. The whole length of the domain is thus 64H. The hill has a height of 1 cm, and the undisturbed boundary-layer thickness is ı 0 D 7 cm at the location of the Gaussian hill. For the simulation of rough-wall flow, the roughness length z 0 is fixed at 8.2 10 -6 µm that is the reference value obtained from the PC09 experiments. Numerical parameters and domain characteristics expressed in wall units for U 1 D 11:2 m s 1 are given in Table 4.5. We note that for A, the coarsest grid, the vertical spacing of wall-bounded grids is z C min D 30, which is a minimum value and falls in the range of the logarithmic region. For B, the vertical grid spacings are refined and the their centres are located at z C min D 15 inside the buffer layer. For C, the most fine grid, z C min D 7:5 reaches the end of the viscous sublayer. These values are in accordance with the classical range of nondimensionalised grid resolutions for WMLES [START_REF] Chen | Wall Modeling for Implicit Large-Eddy Simulation[END_REF][START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF]. The grid resolution in wall units depends on local flow characteristics.The wall shear stress increases and the grid resolution becomes coarser in wall units in regions where the flow accelerates, generally on the windward side. By using the LES data, the friction velocity can be evaluated a posteriori according to u .x/ D p j w .x/j= , based on the values of the wall shear stress w given by the wall model and averaged both temporally and in the transverse direction. In

N x N y N z L x =H L y =H L z =H
x C y C z C mean z C min A 30

Comparisons of grids

As mentioned in Chapter 2, the choice of grid spacing and wall model has a non negligible effect on the simulation. The objective of the first part of this section is thus to study the influence of these factors on the simulation results. The three grids presented in Section 4.2.2.1, A, B and C, are tested for the same case. For the smooth-wall cases, two wall models, SPD based on the wall function of Spalding [1961] and DBM of [START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF], are used. The rough-wall cases are simulated using a modified wall model (RSL) that takes into account the effect of the roughness sublayer. This is presented in details in Appendix A.

From the LES data, a separation line can be identified as the line segment passing through all points where the mean velocity is zero, separating the RZ into seemingly congruent halves. This line is obtained by a linear interpolation of the flow field in the wall-adjacent meshes. By linearly extrapolating the separation line and identifying its intersections with the wall, the positions at which the flow separates and reattaches are obtained, with an uncertainty of one mesh ( x). The separation point, x sep , is defined as the upstream point at which the streamwise velocity reverse its sign compared to the prevailing wind direction, and the reattachment point, x reat , is identified downstream by the opposite sign change.

Comparisons of wall models

Comparing the two wall models, SPD fails to predict flow separation using the coarsest grid A and gives a smaller separation length on grid B. Based on the experimental data of the PC09 case [START_REF] Simoëns | Influence of Gaussian hill on concentration of solid particles in suspension inside turbulent boudary layer[END_REF], the underestimation of the separation bubble on grid B by SPD is reflected in both a delayed detachment and an advanced reattachment. As pointed out by [START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF], this is mainly due to the assumption of an equilibrium boundary-layer related to the wall function of Spalding [1961] on which SPD is constructed. This assumption is not adapted for separated flows. In the following sections, simulations are carried out using the DBM model for the smooth wall cases.

In Figure 4.7, we compare side-by-side mean streamlines obtained from the simulations to the experimental data [START_REF] Simoëns | Influence of Gaussian hill on concentration of solid particles in suspension inside turbulent boudary layer[END_REF]. The separation and reattachment points obtained by LES and by the experiments over the smooth and the rough walls are give in Table 4.7.

Qualitatively, the two approaches agree rather well. We note that even with the finest grid C and the DBM wall model, the simulation underestimates the length of the RZ. The results over the rough wall are in better agreement. Overall, the experimental data give a larger RZ over the smooth wall, contrary to the simulations. The difference between the smooth-and rough-wall experimental cases is less than 10%. We note that in literature, there is still an ongoing debate on whether wall roughness enhances or suppress flow separation [START_REF] Kanda | Particle-Image Velocimetry Measurements of Separation and Re-attachment of Airflow over Two-Dimensional Hills with Various Slope Angles and Approach-Flow Characteristics[END_REF]. Due to the prohibitive cost related to the use of grid C (at least one order of magnitude higher than the other cases), grid B is used in the simulations presented in the following sections.

PC09 simulations 4.2.3.1 Numerical configuration

In order to simulate the PC09 case, a successive-hill numerical configuration is devised based on an extension of grid B in the longitudinal direction (Table 4.8). The domain includes an upstream Gaussian hill, 80H after the inlet, and two Gaussian hills separated by a distance of 3 or 8H, depending on the experimental configuration, as illustrated in Figure 4.8. The centre of the valley between the double hills is located at 71H after the first hill. All three Gaussian hills are of NC shape. The simulation is initiated using a TBL field at U 1 D 11:2 m s 1 from a precursor simulation, and after 100 inertial time scales, flow statistics are collected and averaged over 10 000 inertial time scales and spatially along the transverse direction.

N x N y N z L x =H L y =H L z =H

Simulation results of the successive-hill cases

In this section, the numerical results of the successive-hill cases (3H and 8H) are presented.

As an illustration, the mean streamlines are plotted for these two cases in Figure 4.9. The most noticeable feature is the formation of two recirculation zones between and behind the two Gaussian hills, respectively. The RZ behind the double hills is smaller than the one behind a single NC hill (Figure 4.7). This is probably due to the shear layer formed from the first hill top that induces high velocity fluctuations and enhances exchange, allowing a quicker reattachment behind the second hill. For the 3H case, the outer flow skirts the hills and the two separation bubbles, as if they formed a single envelope. A similar behaviour was already noted in [START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF] for the case around a single obstacle. Comparison on the lee-side flow between the 3H case and a single-hill case will be studied in details using velocity and Reynolds stress profiles in the following paragraphs. As the hill spacing increases from 3H to 8H, there is an expansion of the inter-hill vortex and a deflation of the lee-side RZ (Figure 4.9). For the 3H case, whereas the small size of the first RZ is apparently related to the limited valley size, the second RZ, larger than in the 8H case, implies that the inter-hill wake fluctuations are also less intense. For the 8H case, mean streamlines point downwards into the inter-hill valley at the rear of the inter-hill RZ, before moving upwards along the windward slope of the second hill. This penetration could lead to better momentum exchanges between the valley and the outer flow and a potential solid particle trapping inside the RZ.

In Figure 4.10, we compare the simulation results between the smooth-and rough-wall conditions for the 3H case. We note that for the mean horizontal velocity profiles (Figure 4.10a) and Reynolds stress (Figure 4.10b), the two wall conditions give comparably the same evolution before, between and behind the Gaussian hills. The main difference lies in the mean velocity speed-up at the first hill, in which the smooth-wall value is slightly larger than the rough-wall one. Another difference lies in the small deficiency of the rough-wall Reynolds stress in the wake regions.

The same conclusions are obtained for the 8H case. The corresponding results are shown in single-hill profiles are shifted in such a way that the hill centre coincides with that of the first hill of the successive-hill case, in order to emphasise the extra effects of the double hills.

In Figure 4.11a, we note that although an additional vortex bubble is formed between the crests and the flow separates and reattaches on the second hill, it does not have much consequences on the overall evolution of the outer flow. Starting from X=H D 4, the successive-hill profiles collapse with the upstream-hill ones, albeit with a slightly higher deceleration across the boundary-layer. The same similarity seems to be present in the Reynolds stress profiles, as shown in Figure 4.11b. From X=H D 4 to 1:5, the influence of the downstream configuration is hardly noticeable on u 0 w 0 , whereas the local windward slope seems to be the only influential factor on the flow. Inside the vortex bubble at X=H D 0, the Reynolds stress falls to nearly zero, hinting at a lack of exchange with the outside flow, in conformity with the nature of the skimming flow, described in Introduction. The subsequent wind-ward side of the second hill suppresses apparently flow fluctuations as in the CT06 case, which results in a peak of u 0 w 0 that is 60% lower than the one behind the upstream single hill at X=H D 1:5. Nevertheless, as the flow separates once more over the lee, which contributes to generation of the Reynolds stress hereafter, the two u 0 w 0 profiles collapse at X=H D 6.

In summary, the double hill case with 3H spacing behaves as much as a whole to the upcoming flow. Effects of the double-hill configuration are negligible starting from X=H D 6. The similarity in the mean flow and turbulence levels at this location is quite surprising, if we take into account the differences in sizes and locations of the RZ behind the single hill and between double hills.

For the 8H case presented in Figure 4.12, profiles behind the first hill (between X=H D 4 and 1) are replotted with a 8H shift to the right, in order to compare with profiles behind the second hill. Similarly to Figure 4.11, the upstream single-hill profiles are also plotted in For the 8H case, the mean flow evolution before the first hill does not deviate much from that of the 3H case. We note the resemblance between the mean velocity profile at X=H D 4 in Figure 4.12a and the one at X=H D 1:5 in Figure 4.11a. Over the second hill top, the flow speed up is also reduced from the one at the first hill top, simmilar to the 3H case. The main difference from the 3H case lies in the lee side of the second hill, where the flow separates but reattaches much more quickly, depsite the reduced speed-up at the hill top, compared to the flow seperation behind the first hill. Figure 4.12a shows that, on the lee side of the first double-hill, the separated shear layer seems to be less developed than behind the upstream single-hill, shown by the smaller velocity gradient between X=H D 2 and 2. At X=H D 2 in Figure 4.12b, the peak of u 0 w 0 at X=H D 2 starts with a value 10% higher than its counterpart after the upstream single-hill.

Yet, this u 0 w 0 peak is quickly overtaken at X=H D 0, probably due to the presence of the second hill. At the same time, an upwards spreading of u 0 w 0 to the higher portion of the flow appears, even before the encounter with the second hill . This implies an extension of vertical mixing thanks to the double hills, in conformity with the nature of the wake flow (Introduction). The second hill allows the flow to reattach more rapidly, at X=H D 2, than in the upstream-hill case. On the other hand, the second hill seems to have little effect on the shear layer above Z=H D 1. Behind the crest, the flow separates but reattaches shortly afterwards. We note a lowering of the u 0 w 0 peak at X=H D 5:5 in relationship to the upstream-hill case, roughly at the centre of the second RZ. This hints at the smaller size of the second separation bubble. Another false "peak" seems to be overshadowed at Z=H D 2, which is probably a history effect from the upstream u 0 w 0 peak at X=H D 4. As the flow reattaches shortly after X=H D 7, the peak of u 0 w 0 falls to Z=H D 1 and its lintensity, albeit smaller than its counterpart at Z=H D 1 (open squares), remains stronger than the one after the upstream single-hill (filled squares).

The richness in the variety of flow characteristics in the successive-hill configuration motivates us to use this geometry to study the particle transport, as presented in the next chapter.

Such richness could be linked to the conceptual flow regimes over obstacles proposed by Oke [1988]. The results obtained in this section confirm that the 3H case presents skimming flow behaviours and the 8H case the wake flow ones.

Study of the mean flow separation and recirculation behind a 2D hill

In this section, a parametric study on the recirculation zone behind a 2D hill with sinusoidal or Gaussian geometry is presented. First, the simulated case are presented. Section 4.3.2 is focused on the influences of the hill slope and flow parameters on the separation and reattachment points. The effects of the Reynolds number, the wall condition and the blockage ratio are investigated. In the second part (Section 4.3.3), the aim is to characterise the separation bubble under different flow/hill conditions using adequate geometric and kinematic parameters. The dependence of these attributes on the hill geometry, the Reynolds number and the wall condition are presented and discussed.

Description of the simulation cases

A series of sinusoidal hills as well as the two Gaussian hill shapes, AC and NC, from the PC09 cases are used in this study. The maximum hill slope angle varies from 9°to 32°. The hill height H is 1 cm for all cases. These hill geometries are illustrated in Figure 4.13, where the streamwise and vertical coordinates are normalised by H. In order to investigate the effects of the Reynolds number and the blockage ratio, the free-stream velocity U 1 is varied from 2:3 to 11.2 m s -1 , and the thickness of the boundary-layer from 5 to 10 cm, covering a Reynolds number range of Re D 500--3000, and a blockage ratio variation from 0:1 to 0:2. For all simulations presented in this section, the single-hill configuration and grid B presented in Section 4.2.2.1 are used. The simulation is initiated using a TBL field obtained from a precursor simulation after the convergence of flow statistics. For the highest flow regime, statistics are collected and averaged over a period of 1000 ı 0 =U 1 and spatially along the transverse direction.

Flow separation and reattachment

Critical angle of flow separation

The aim of this section is to examine the critical hill slope that leads to the flow separation on the lee side. In literature, the experimental value of this critical slope angle varies around 14°and 16° [START_REF] Cao | Numerical study on turbulent boundary layers over two-dimensional hills -Effects of surface roughness and slope[END_REF][START_REF] Kanda | Particle-Image Velocimetry Measurements of Separation and Re-attachment of Airflow over Two-Dimensional Hills with Various Slope Angles and Approach-Flow Characteristics[END_REF]. The problem is not trivial. Whereas obstacles with sharp edges, e.g., backward-facing step [START_REF] Cabot | Near-wall models in large-eddy simulations of flow behind a backwardfacing step[END_REF] and triangular ridges [START_REF] Finnigan | Air Flow Over Complex Terrain[END_REF] generally force the flow to separate at the edge, the curly surface of the hill does not induce flow separation in a systematic manner. The hill geometry is thus a determinant factor on the flow behaviour behind it.

To investigate this problem, the generic sinusoidal hill shape defined by Equation 4.1 is used.

This hill geometry has been widely adopted in previous studies [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], 2007;[START_REF] Gong | A wind tunnel study of turbulent flow over model hills[END_REF][START_REF] Tamura | LES study of turbulent boundary layer over a smooth and a rough 2d hill model[END_REF] with various hill slopes. The ratio L=H determines the geometry of the hill. If L=H ! 1, the sinusoidal hill tends towards a flat surface. Starting from L=H D 10, we gradually decrease this ratio, which in turn increases the maximum slope of the hill. The free-stream velocity is fixed to 11.2 m s -1 , a flow regime sufficiently strong to initiate and sustain solid particle saltation according to the PC09 experiments.

In Figure 4.14, the gradual appearance of a downhill recirculation zone is captured by the mean flow streamlines around a 2D hill with growing slope angle from 9°(L=H D 10) to 32°( L=H D 2:5). The ratio L=H clearly has a significant influence on the flow dynamics especially after the hill. At L=H D 10 (9°), the flow before and after the hill is overall symmetric, except very near the wall on the lee side. Apparently, the flow inertia is strong enough to resist the unfavourable pressure gradient created by the hill. This attached-flow regime seems to break down at L=H D 7:5 (12°). The LES results show that the flow separates but rapidly reattaches on the lee side. At L=H D 5 (18°), the flow separation becomes quite visible. A small separation bubble appears on the lee side of the wall, apparently squeezed by the outer streamlines. Finally, a large RZ is established firmly behind the hill at L=H D 2:5 (32°), occupying a significant area with extents comparable to the dimensions of the hill. Globally, these results are in accordance with the experimental results from literature showing that the critical angle of flow separation behind a 2D hill is between 14°and 16° [START_REF] Finnigan | Air Flow Over Complex Terrain[END_REF][START_REF] Kanda | Particle-Image Velocimetry Measurements of Separation and Re-attachment of Airflow over Two-Dimensional Hills with Various Slope Angles and Approach-Flow Characteristics[END_REF].

For the mean flow, streamlines adhere to the surface without separation on the L=H D 10 sinusoidal hill. Yet, a more careful look into the lee-side flow field reveals that weak flow separations do occur intermittently over this slightly inclined surface. The duration of these "separations" is in general very short (barely 0:5% of T L ), which relates to the absence of a mean separation bubble. This intermittency in flow separation behind a 2D hill has already been studied by [START_REF] Fröhlich | Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions[END_REF] and [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] on cases in which a large flow separation occurs. The main conclusion of the authors is that a single separation "streamline" is not meaningful for the instantaneous flow due to the high irregularity of the flow. In the following, results from the time-averaged flow separation line are presented.

The streamwise locations of the separation and reattachment points and distance in-between (denoted as the separation length) are summarised in Table 4.10 for each case. The experimental results of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] and Cao & Tamura [2007], albeit at a higher Reynolds number, are added for comparison. Globally, the results reveal that flow separation is more likely to occur when the hill slope is high. In fact, for the maximal slope angles between 12°and 32°, the length of RZ increases approximately linearly with hill slope, as shown in Figure 4

Influence of the Reynolds number and of the wall condition

From the PC09 simulations, we have seen that the NC hill is able to generate a large RZ at U 1 D 11:2 m s 1 . In order to study the Reynolds number dependency of the flow separation on this hill, flow velocity of the upstream TBL is systematically varied from U 1 D 2:3 to 11.2 m s -1 , for both smooth-and rough-wall conditions, as shown in Table 4.9. The undisturbed boundary-layer depth is fixed (ı 0 D 0:07 cm). 4.10.

The LES results show that on a smooth wall, an increasing Reynolds number shortens the separation length on the NC hill in a substantial way (Table 4.11). Interestingly, we find that the separation point is more sensible to Re H than the reattachment point at high Reynolds numbers. For the rough wall, an increasing Reynolds number increases the flow separation length, although its influence is much weaker than on a smooth wall (Table 4.12). This increase is in accordance with the numerical result of [START_REF] Araújo | Numerical modeling of the wind flow over a transverse dune[END_REF] for a RZ behind a 2D sand dune.

Influence of the blockage ratio

In [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF], the authors judged that the blockage ratio of H=ı 0 D 1=6 is small enough that the effects of the TBL depth can be omitted. In this section we set out to examine the blockage effect by a 2D hill on the TBL, while maintaining the other flow characteristics (U 1 and u ) and the Gaussian hill geometry constant. Both AC and NC hill shapes are tested, with the flow velocity U 1 fixed at 11.2 m s -1 . The results are summarised in Table 4.13.

For the smooth wall, a higher value of H=ı 0 delays the separation and advances the reattachment. This could be explained by the fact that a higher blockage produces a Venturi effect that allows the flow to better resist separation. As the blockage ratio increases from 1=10 to 1=5, the separation region in the AC smooth cases grows by 20%. In the NC smooth cases, the increase is around 7%. Moreover, the blockage ratio has negligible effects if H=ı 0 < 1=7.

For the NC rough cases, a change in the blockage ratio does not generate significant variations in the flow separation. Their values are too close to the uncertainty in the determination of separation/reattachment points (0:1H).

In conclusion, it seems that in the range between 1=10 and 1=5, the blockage ratio has a nonnegligible inhibiting effect on flow separation behind a mild, smooth hill. Over a steep hill, this effect seems to be only relevant between 1=7 to 1=5, whereas over the rough wall, its influence is weak over the whole range.

Parametric study of the separation bubble

Identification of the geometrical and kinematic parameters

In the previous sections, the flow separation was studied using the separation/reattachment locations and the separation length, which give an overall description of the separation bubble.

Nevertheless, whereas the separation length provides an estimation of the streamwise extent of the RZ, the extent of the RZ in the vertical direction, equally crucial to the potential particle trapping, cannot be fully appreciated using this parameter.

Before proposing other parameters more suitable to the characterisation of the RZ geometry, a proper border of the RZ needs to be defined. We define this border using the recirculation streamline, which corresponds to the longest streamline forming a clockwise loop behind the hill. Streamlines outside this loop only skirt over the RZ. Then, the bulb-shape RZ can be characterised by a major axe (L CD ) aligned with the hill surface, and a minor one (L AB ) perpendicular to the hill surface.

The major axe is the longest segment connecting the two extremities of the RZ along the wall surface. The minor axe is defined using the maximum distance between the recirculation streamline and the hill surface in the vertical direction. This distance quantifies the extent to which the outer streamlines need to deviate from the hill to circumvent the RZ, thus giving a measure of the RZ intensity complementary to the major axe. The point at which the two axes intersects is denoted as the centre of the RZ. As shown in Figure 4.16, the two axes L AB and L CD quantitatively describe the 2D dimensions of the RZ: when the hill slope steepens, the RZ expands along the two axes as both L AB and L CD react by increasing their values. Finally, the circumference of the RZ, L R , is another measure of the RZ size. A characteristic time scale R can be constructed as

R D Z L R dl u l : (4.4)
Here u l is the component of the fluid velocity tangential to the recirculation streamline. R corresponds to the time duration that a fluid parcel takes to complete a clockwise loop around the RZ. This kinematic parameter can be used to estimate the intensity of the RZ: smaller the value of R , more efficient the RZ will be in the transfer of mass and momentum. This time scale is later used to compare with the trapping time durations of the solid particles in Section 5.4.

The aforementioned parameters are derived in the cases with ı 0 D 7 cm, from Table 4.9, and summarised in Table 4.14. In several cases, the hill slope is too mild and the flow separation is too weak to provide measures of all the geometric parameters such as the minor axe length. In these cases, only the RZ centre and the major axe length are given. Cases without mean flow separation are marked by a minus sign.

Two observations can be made on the influence of the hill slope:

1. By increasing the hill slope, the centre of the RZ at .x c ; z c / moves both upwards and towards the hill top. The former is purely an aerodynamic behaviour, since the decrease in L results in a priori a downward and backward shift of .x c ; z c /. The lengths of the minor and the major axes respond both proportionally to the hill slope increase. Among the two, the growth of the minor axe seems to be more significant, since the ratio L AB =L CD grows with the hill slope.

Re

Angle 2. As expected, the circumference of the RZ grows with the hill slope. The RZ time scale R decreases with it. This may be a sign of the enhancement of the backflow inside the RZ strong enough to overcome the increase of L R .

x c =H z c =H L AB =H L CD =H L AB =L CD L R =H U 1 R =L R 500
The dependences of the RZ characteristics on the Reynolds number and the hill geometry are presented and discussed in the following sections.

Reynolds number effects in smooth-and rough-wall cases

In Figure 4.17, we compare side-by-side the effects of the wall roughness on the size of the RZ behind the NC hill. The free-stream velocity increases from 2.3 to 11.2 m s -1 . First, qualitative characteristics of the RZ are essentially the same between 7.9 to 11.2 for both wall conditions. This implies an independence of the Reynolds number effects. Secondly, the rough wall hill gives larger RZ than their smooth-wall counterpart in general. This can be attributed to the smaller velocity near the wall, which makes the flow more inclined to be separated on the lee side.

In Figure 4.18, we note that in smooth-wall condition, the Reynolds number effects on the RZ length, L R , are strong at small velocity, with a 40% decrease for cases from 2.3 to 11.2 m s -1 . At high velocity, the abrupt increase of the RZ size at the largest Reynolds number could not be easily explained. This suggests, that the current LES results on the smooth wall may indeed underestimate the RZ dimensions, as suggested in Section 4.2.2.2. For the rough wall condition, the RZ length L R hardly varies with Re H , except for the case of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] with a Reynolds number of order 10 4 . We note that the rough wall in this case presents a relative roughness length z 0 =H six time higher than in other cases. On the other hand, the form ratio between the two axes of the RZ evolves slowly with Re H , with a maximal increase of less than 20%. In general, whereas the shape of the RZ seems to be not sensible to the wall condition, the size of the RZ varies greatly between rough and smooth walls, and this discrepancy increases with the flow velocity. These tendencies are in accordance with the evolution of the separation length with Re H , presented in Tables 4.11 and 4.12.

Comparison between Gaussian and sinusoidal geometries

A priori, the maximal slope of the hill is not the only factor that determines the flow behaviours after the hill. Here we compare the numerical results obtained on two similar hill geometries, a Gaussian hill (NC) and another sinusoidal hill with L=H D 2:5, both in rough wall conditions, for three flow regimes. The main difference between the two hill shapes lies near the foot where the hilly terrain meets the flat surface. As shown in Figure 4.13, the transition is longer and smoother for the NC hill. Since their maximal slope angles differ only by 1°, we expect that any difference in the RZ is mainly caused by the detailed geometric shape of the hill.

As shown in Figure 4.19, the longer hill foot of the NC hill does not seem to induce noticeable difference in the RZ. For all cases, the recirculation zone generally reaches a little less than 6H behind the hill top. The RZ shapes are also quite similar and implies a Reynolds number independence for these two hill types as early as 7.9 m s -1 (Re H D 4000). Both the RZ length and the form ratio present a slight increase with the Reynolds number before Re H D 10 4 , as shown in Figure 4.20. The higher values of L R at Re H D 10 4 are likely to be caused by the cube roughness used in the CT06 case, which is larger than in the other cases relative to the hill height. Globally, in the range of Reynolds number studied here, the difference between Gaussian and sinusoidal hill is quite small in terms of the RZ characteristics. The aforementioned results obtained from sinusoidal hills could thus be comparable to Gaussian hills with similar hill slope.

As the objective of the thesis is the solid particle transport and trapping over 2D hills, we have not considered the effect of other roughness parameters, such as the roughness density. We note that this roughness effect could also be important and needs to be studied in the future.

A priori study of particle lift over terrain

The hill-induced modifications on the wind are captured by the instantaneous flow field from LES. These data, especially in the near-wall region, allow to evaluate the particle lift force according to the take-off model presented in Section 3.2.4. Maps of particle lift statistics can thus be plotted, and patterns of wind erosion and particle deposition over the Gaussian hills can be estimated a priori. To this aim, we define a fictive particle of diameter 200 µm and density 1000 kg m -3 , according to the mean diameter characteristics in the PC09 experiments. Later, the particles are assumed to be homogeneously distributed at each grid point on the wall. Exposed to the local flow and fixed to the wall, they are subject to the forces of lift, gravity and adhesion.

The statistics of the lift force are reported in this section.

Based on Equation 3.33, the formulation relating particle lift to events of strong turbulent structures writes,

F L D 15:5 2 p w = D p ! 1:87 u 0 ˚w0 « hu 0 ˚w0 « i y ; (4.5)
the overline denotes a temporal average and the bracket denotes a spanwise average. Quantities such as w , u 0 ˚w0 « and u 0 w 0 are evaluated locally and instantaneously from the LES-resolved flow. We point out that in Equation 4.5, u and w are projections of the near-wall wind speed in the tangential and wall-normal directions to the surface, respectively. The second relationship relates the mean lift to the mean wall shear stress using the formula given by Equation 3.26:

F L D 15:5 2  u D p à 1:87 : (4.6)
In this estimation, only the temporal average of a term proportional to u D p j w j= is needed to obtain the mean lift. As opposed to Equation 4.5, Equation 4.6 is readily applicable to eddymodelling simulations such as RANS.

Figure 4.21 shows that the mean particle lift over the flat surface obtained using event-based method (Equation 4.5) is almost doubled from the one obtained from the mean wall shear stress (Equation 4.6). The discrepancy between the two estimations becomes smaller on the windward side of the hill. Inside the RZ, the second evaluation apparently overestimates the lift force.

After X=H D 6, the event-based evaluation generally gives larger, and sometimes doubled, values compared to the mean one.

Qualitatively, both estimations identify the windward sides of the hills as highly susceptible of wind erosion, whereas both the valley between the double hills and zones behind the hill have high potential for particle deposition, due to the weak particle lift in these areas. These results confirm that the RZ occupying large spaces between and behind the hills is likely to trap incoming particles. In the next chapter, solid particles are released from an upstream position and the particle trapping and deposition are duly investigated.

Conclusions

In this chapter, the aerodynamic features of the flow over one, or several, 2D hills are studied.

The formation of recirculation zones on the lee side of the hills has received most attention.

Simulations described in this chapter not only represent an application of the LES code ARPS to flows over complex terrain, but also lay groundwork for the study of particle transport around the Gaussian hills, presented in the next chapter.

The beginning of the chapter presents a validation case by using a reference experimental case CT06 of a TBL flow over a steep, sinusoidal hill. Velocity and second-order flow statistics are compared with the experimental data of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF]. Later, the experimental project PC09 related to the current thesis is presented. The single-hill configuration is first used to test different wall models and grid configurations. The wall model of [START_REF] Duprat | A wall-layer model for large-eddy simulations of turbulent flows with/out pressure gradient[END_REF] is shown to give a RZ size closer to the experimental value for the smooth-wall case, than the wall model based on the wall function proposed by Spalding [1961] for a canonical TBL flow. Nevertheless, a sufficiently high grid resolution seems to be the primary factor affecting the RZ size predicted by LES, especially at high Reynolds numbers. The rough-wall case is proved to be less severe in this regard, as simulations with a medium grid resolution give comparable results with the experiments. By adopting these numerical configurations, the successive-hill cases are simulated. The results of the two cases, 3H and 8H, are first compared with the upstream-hill profiles. Differences in the mean flow and the Reynolds stress fields between the two cases are also discussed. It is shown that globally, the mean flow behaviour of the 3H case can be characterised by the skimming flow, and the 8H case by the wake flow, according to the terminology of [START_REF] Oke | Street design and urban canopy layer climate[END_REF].

The third part of chapter is devoted to a parametric study of flow separation behind a 2D hill, of Gaussian or sinusoidal geometry. We first study the hill slope at which the flow separates from its lee side. At Re H D 6000, the critical angle of flow separation is found to be 12°over a smooth sinusoidal hill. Then, we study the influences of the Reynolds number, wall conditions and the blockage ratio on the locations of flow separation and reattachment. On the smooth hill, the larger Reynolds number seems to postpone flow separation and advance reattachment, whereas the opposite is true on the rough hill. We further identify several parameters describing the geometric and aerodynamic properties of the separation bubble. The dependence of these parameters on flow and wall conditions is subsequently studied. To the best of our knowledge, this is the first time that this type of analysis are carried out in the study of flow separation behind 2D hills. Finally, values of particle lift are deduced from the simulation data for the successive-hill cases. Two methods of estimation are tested. The first method takes full advantage of the instantaneous, local flow field from LES. The second one relies solely on the value of the mean wall shear stress, which has the merit of being widely available in eddy-modelling simulations such as RANS. The discrepancy between the two estimations is shown to be as large as 100% over the flat surface, owing to eddies of intense events that are not captured by the mean wall shear stress in RANS models. The windward side of the hill is identified as susceptible of wind erosion, and the lee side as regions of potential particle deposition, based on the relative intensity of particle lift to gravity and adhesion. This a priori estimation reveals the potential zones for particle trapping and deposition, which will be looked in details in the next chapter.

Chapter 5 LES study of solid particle transport

In Chapter 4, some aerodynamic aspects of Gaussian hills for isolated and grouped configurations have been studied. In particular, the properties of the recirculation zones due to the terrain have been thoroughly studied. It allows to obtain some geometric properties of RZ as a function of flow regimes and hill spacings. In the same vein as in the experimental campaign of PC09, the various Gaussian-hill configurations are maintained as solid particles are introduced from an upstream position inside the TBL. The objective is now to determine whether there are some particularities of ejection or deposition of solid particles due to RZ, and if so, how does these patterns vary with flow regimes and hill spacings.

The experimental and numerical configurations are duly presented in Section 5.1. Based on the modelling approaches presented in Chapter 3, the particle transport over a flat wall is first studied in Section 5.2. Later, in Section 5.3, particle saltation over Gaussian hills are studied using concentration and velocity profiles. An example of the simulation results is illustrated qualitatively in Figure 5.1. An evaluation of particle forces, used in our force-balance model of particle take-off, is also made around the Gaussian hills. Finally, we tackle the problem of particle trapping in the context of this thesis in Section 5.4. Based particle deposition maps, preferential regions of deposition are first identified. The particle trapping is formulated in the context of Gaussian hills configuration. The outcome of this particle trapping and its dependence on the Shields number and hill-spacing are subsequently analysed using simulation results.

Configuration

We first present the experimental case that is used to validate our complete model for solid particle transport. At our best knowledge this is the only experimental case that presents quantitative studies on both aerodynamic aspects of the RZ and solid particle concentration. After will be presented the numerical configuration to mimic this case.

Experimental configuration

In the experimental campaign, solid particles are supplied continuously from a feeding device on the wall upstream of the Gaussian hills. These grains have a mean diameter of 200 µm, with variations between 170 and 250 µm, and a density around 1000 kg m -3 . We point out that particles of same characteristics are glued all along the wall and over the hill surface in order to obtain rough wall conditions close to a static sand bed. The particle density is smaller than the classical value of sand density (between 2500 and 2750 kg m -3 ). This is to approach values in ABL particle transport, in terms of the Stokes and Shields numbers as given in Table 5.1, because of the reduced length-and time-scales inside a wind tunnel compared to natural conditions. We note that for all flow regimes, the solid particle motion in the PC09 case belongs to the "modified saltation" type, according to the condition St L < 1 and g > 1 [START_REF] Taniere | On the behaviour of solid particles in a horizontal boundary layer with turbulence and saltation effects[END_REF]. In this thesis, the name "sandbox" denotes the particular particle-feeding device that was used experimentally in order to supply solid grains from a tank located beneath the working section of the wind tunnel. Powered by an upward-moving piston at a fixed pace, particles are constantly introduced into the incoming TBL flow via a 20 cm 10 cm slot, as shown in Figure 5.2. We point out that the width of the wind tunnel is 20 cm. We note the relatively long extent of the sandbox compared to the hill dimensions. This allows a substantial and continuous supply of solid particles ready to be set into motion, if the entrainment condition is fulfilled.

D P (µm) p;0 (s) u (m s -1 ) T L (s) T Á (
Particular attention was paid to avoid the gradual accumulation of grains at the sandbox during all recordings of experimental data. This in turn regulates the regime of the electric motor to which the piston is connected according to the upstream TBL flow. At equilibrium, the piston velocity has been gradually increased to a stable value, in order to preserve a constant shape of sand bed at the box. The only intrusion to the flow induced by the sandbox is the continuous introduction of solid particles. Moreover, a direct visualisation of the incipient particle motion is possible, which cannot be done if the grains are released/bombarded from a source located high up of the working section. In summary, the sandbox feeding assures a stabilised supply of grains and avoids spurious splash phenomenon in an unnatural way, without altering other TBL flow characteristics except the introduction of solid particles.

N // D S box p 3D 2 P =4 3 6 D 2 p 3 S box D 2 P ;
The coefficient 3=6 accounts for the ratio between the number of circles (projection of spherical sand grain) to the number of equilateral triangles,(marked by dashed lines in Figure 5.3a). Let we define box as the plane roughness density of this particle layer: where N // is the number of trains with this layer. Inside a vertical view, as shown in Figure 5.3b, the number of layers sent into the TBL per unit time due to the upward motion of the piston is

N ? D P Q v dt H p ;
H p being the distance between successive layers of grains, which equals to p 2=3D p if particles are closely packed. The volume fraction of solid particles contained in the unit volume is thus

Á D N ? N // V P P Q v dt :
The mass flow rate at the sandbox is

P m 0 D S box P z Á p : (5.1)

Numerical implementation

Owing to the special design of the sandbox configuration, its numerical implementation of the particle source is quite straightforward. Due to the slow upward speed at the sandbox (of the order of several millimetres per minute) particles supplied at the sandbox has a negligible initial velocity when first exposed to the flow. This is reflected in the LES that particle motions at the sandbox are uniquely initiated through aerodynamic entrainment, without any arbitrary, numerical initiation of particle velocity or altitude. Instead of imposing a fixed rate of particle emission, after the instant at which a grain is set into motion, another one is systematically added for potential entrainment at future instants. We assume that this implementation conforms with the experimental set-up, in which the supply rate of solid particles is only dictated by the momentum of the flow, and no net erosion occurs and modifies the surface at the sandbox.

In this chapter, two numerical configurations have been used for the study of solid particle transport, denoted as the NC and PC09 configuration. The mesh size of NC in the same as the grid B, presented in Table 4.4, and that of PC09 have been given in Table 4.8. Both the NC and PC09 configurations have been originally used for the TBL simulations in Chapter 4.

We point out that the NC configuration is also adapted into a single-hill case in Section 5.4, with an isolated Gaussian hill located 50H downstream of the inlet, similar to the single-hill configuration presented in Section 4.2.2 without particle emission. The PC09 configuration includes an upstream, isolated hill not far from the sandbox location and two downstream hills relatively close to each other, with a spacing of 3H or 8H. These two cases are further identified as the 3H and the 8H cases. Numerical studies of the PC09 cases are mainly conducted within the following zones:

1. Z01 : this zone includes the sandbox of dimensions 20H 6H (longitudinal spanwise). Its origin is located at the upstream edge of the sandbox, 23H after the inlet.

2. Z02 (the upstream Gaussian hill) : this zone is centred longitudinally at the upstream hill crest, 36H after the sandbox.

3. Z03 (the double Gaussian hills) : the origin of this zone is at the midpoint between the two hill tops, located 71H downstream of the upstream hill top.

The closely-packed configuration implies a high particle volume fraction, most likely larger than 1 , near the sandbox. The inter-particle actions is probably not negligible in this location according to [START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF]. Similarly, due to the constant supply of resting grains on the wall, the splash could be of importance as well. Nevertheless, both four-way coupling and splash of particles are not included in the simulation results presented in this chapter. To our best knowledge, no well-validated model exists in the literature on the inter-particle interactions and splash in the case of incipient particle motion. Yet, it has to be noted that at least in our cases, these two processes tend to cancel out each other, since the inter-particle friction and collision hinder all relative motions between grains and thus, should mitigate effects due to the "bombardment" of impacting particles, such as particle splash.

In the following sections, the simulation results of particle transport within a TBL case without Gaussian hill are presented and compared to analytical relationships in literature. The results related to the PC09 cases within Z03 are presented subsequently.

Particle transport inside a turbulent boundary-layer

Description

In order to validate the simulation of particle transport, we study a developing saltation layer inside a TBL. The domain is of the same dimensions as the NC configuration. It starts from the aerodynamic entrainment at the sandbox and end at X=H D 64, the end of the domain (Figure 5.4a). Particles reaching the outlet are evacuated from the simulation. In order to study the influence of the Shields parameter, several flow regimes have been used, with U 1 ranging from 2.3 m s -1 to 20 m s -1 .

In the next sections we will focus on one particular parameter, the saltation mass flux, by comparing the simulations results to the experiments as well as predictions of empirical formulas from the literature.

Saltation mass flux

By definition, the streamwise mass flux q x .x; z/ corresponds to the mass of solid particles passing through a unit area, perpendicular to the longitudinal direction, per unit time. The mass flux has thus the dimensions of [kg m -2 s -1 ]. In our case, the unit area is equal to L y z p at given x and z. L y is the spanwise length of the computational domain, which is equal to the width of the numeric sandbox. A longitudinal averaging interval of x p D 0:1H is further used in the calculation of the mass flux. The flux is thus averaged inside a volume V, between x and x C x p , z and z C z p , and across L y . The calculated mass flux writes

q x .x; z/ D 1 x p L y z p X V m p u p ; (5.2)
where m p and u p are the mass and the streamwise velocity of individual particles, respectively.

In the results presented here, x p and z p are taken as 0:1H. The unit time related to the calculated mass flux in Equation 5.2 is thus equal to an ensemble average , within V, of the time that each grain takes to travel through x p .

It is often useful to examine the saltation mass flux across the whole saltation layer. By definition, the saltation mass flux, Q, is the vertical integration of the streamwise mass flux q x that is function of the altitude z.

Q.x/ D Z 1 0 q x .x; z/dz : (5.3)
This quantity represents the total amount of sand grains in movement measured at a fixed streamwise position. In The saltation flux is a crucial parameter for the study of wind erosion, as it can be obtained with relative ease from both in-situ and laboratory experiments. Since the seminal works of [START_REF] Shields | Anwendung der Ähnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung[END_REF] and [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF], much effort has been devoted to formulating equations that effectively predict the saltation mass flux [START_REF] Bagnold | The physics of wind blown sand and desert dunes[END_REF][START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF][START_REF] Kawamura | Study on sand movement by wind[END_REF][START_REF] Lettau | Experimental and micrometeorological field studies of dune migration[END_REF][START_REF] Owen | Saltation of uniform grains in air[END_REF][START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF][START_REF] Sørensen | On the rate of aeolian sand transport[END_REF]. As a result, various saltation models are available in the literature that relates the saltation flux to the characteristics of the flow and solid particles.

Several of these models are compared to the simulation results in Figure 5.5. Further details of these formulas are given in Table 5.2. We note that these formulas generally establish a powerlaw evolution with the friction velocity at high Shields numbers. Relationships of Bagnold and Owen [1964] predict cubic dependence of Q on the friction velocity u , whereas the models of Sørensen [1991] and of [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF] give Q s / u 2 . At low Shields numbers, a transitional regime exists starting from the threshold range of particle saltation. According to the aforementioned saltation models, this regime is apparently regulated by a miscellaneous set of nondimensional ratios such as W t =u t and u t = p gD P , noted collectively as ˇin Table 5.2.

Several coefficients of these models are determined from experiments with specific particle and flow characteristics. The relevant details are given in Table 5.3.

Source Q=Q 0 Constants Bagnold [1941] ˛V 3 ˛D 1:8 p D p =D ref ; D ref D 250 ţm Owen [1964] .1 V 2 / .˛C ˇV 1 / V 3 ˛D 0:25 , ˇD 0:33W t =u t Sørensen [1991] .1 V 2 / .˛C ˇV 1 / V 2 ˛D 3 , ˇD 3:9 Creyssels [2009] ˛.1 V 2 / V 2 ˛D 28 = p 2=3 , D .u t = p gD P / 3

Note

Obtained from the take-off curve of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] Table 5.3: Flow and particle characteristics of the cases used to compare with our simulation results and experiments.

For Sh high enough (Sh > 0:1), our results exhibit a behaviour similar to that predicted in [START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF]. It has to be noted that the splash is not present in our simulations, corre-sponding to our experimental mimicking. Besides, our values of friction velocity u are taken from the undisturbed TBL, which are probably greater than ones determined inside a developed saltation layer, as in other experimental studies. This could explain the apparent shift between our results and the saltation models shown in Figure 5.5.

Threshold regime of particle transport

Using ARPS, the simulation case with u D 0:188 m s 1 gives a small but non-negligible fraction of grains entrained and transported by the flow, with Q D 0:28 g m 1 s 1 . At u D 0:111 m s 1 , the computed mass flux drops to an insignificant value of 2 10 5 g m 1 s 1 (not shown in Figure 5.5). Apparently, the wind intensity is too weak to entrain grains into the central region of the flow. According to the take-off curve of of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] (Equation 3.24), the threshold friction velocity at which solid particle transport can be sustained is 0.157 m s -1 for PVC particles of 200 µm. The LES results are thus in accordance with the prediction of the take-off curve and with the experiments [START_REF] Simoëns | Rapport ANR PEDO-COTESOF[END_REF].

For sand particle with 242 µm and of 2500 kg m -3 , used in the experiments of [START_REF] Sørensen | On the rate of aeolian sand transport[END_REF] and [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF], the threshold friction velocity is estimated to be 0.258 m s -1 using the take-off curve. Yet both authors reported smaller values in their experiments (Table 5.3). We remind that in our case, particles are entrained into motion from the sandbox on the wall contrary to the sand feed that injects particles from an upper location, commonly used in other wind-tunnel studies such as [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF] and [START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF]. In our opinion, this sand feed runs the risk of increasing, at low u , spurious splash phenomena compared to in-situ and real conditions. For high u , this effect is less important as injected particles tend to be more rapidly evacuated by the flow without reaching the ground. Since particle splash is absent in our configuration, this artefact does not exist in our cases, which could explain the accordance between our results and the prediction of the take-off curve.

Results of particle transport over Gaussian hills

In this section, we simulate particle motions as in the PC09 experimental configurations in order to study the particle transport over successive Gaussian hills, at the highest Shields number of Sh D 0:12 available in the PC09 experiments (U 1 D 11:2 m s 1 ). Cases are first initiated using relevant aerodynamic cases without particles as presented in Section 4.2.3. Particles are later released from the source once the aerodynamic convergence is achieved.

The attention is given to the key zones identified in Section 5.1, namely Z01, Z02 and Z03 (Figure 5.4). We first validate the simulation results on particle concentration using the PC09 experiments. Then, particle velocity profiles and the evolution of forces exerted on particles on the wall are given for the same zones. In the end, we present the wall shear stress around the Gaussian hills and relate it to the particle deposition map.

Particle concentration profiles

The particle concentration presented in this section is defined as the mean volume fraction occupied by all the solid particles inside the unit volume V as defined in Section 5.2.2

c.x; z/ D 1 x p L y z p X V V P ; (5.4)
where V P is the volume occupied by individual solid particles comprised in the unit volume. In Figure 5.6, concentration profiles from locations upstream of the double hills are plotted. The coordinate origin is fixed at the centre of the upstream Gaussian hill in this figure. Simulation results (full line) are compared with the experimental ones (points). The concentration profiles, obtained from the simulation results, are normalised by their local maximum, given in Table 5.4. Each profile is delimited using two dashed lines marking zero and maximum, respectively. The coordinate origin is set at the top of the upstream hill.

X=H

-25.00 -3.50 -1.50 0.00 1.50 3.00 c max 1.12 10 -3 2.40 10 -4 5.87 10 -4 2.18 10 -4 1.23 10 -4 9.58 10 -5 Z=H .c max / 0.00 0.00 0.33 1.00 1.73 2.11

Table 5.4: Peak volume fraction and vertical location of the peak around the upstream hill from the simulation results.

The first profile in Figure 5.6, at X=H D 25, corresponds to the concentration profile near the downstream edge of the sandbox. A non-negligible quantity of solid particles are constantly entrained and set into saltation. The vertical extent of this developing saltation layer reaches roughly 1:5H at a distance of 5H away from the source.

The subsequent profiles in Figure 5.6 describe how the particles as a whole react to the change of terrain from a flat surface to a steep, transverse hill. Globally, the agreement between the simulation result and experimental one is reasonable: the locations of the concentration peak at various positions are well predicted by LES. Nevertheless, the simulation seems to overestimate the saltation height growth behind the Gaussian hill. In fact, a high percentage of solid particles are moved into the upper part of the TBL (z > H) behind the hill. This upward migration is less dramatic in the experimental results nevertheless. Besides, the simulation apparently underestimates the quantity of solid particles passing through the RZ compared to the experiments.

The predicted extension of the saltation layer over the Gaussian hill could be explained by the fact that solid particles, impacting on the windward side of the hill, are abruptly redirected upwards due high hill slopes. As an inertial effect, these rebounding grains maintain this upward motion before reaching high layers of rapid flow motion, which immediately evacuate them to the downstream of the hill region.

Figure 5.7 shows the concentration profiles around the two hills separated by 3H. A better collapse with the experimental profiles is achieved around the second hill. We note in particular that a small fraction of particles are concentrated inside the small vortex between the crests. The simulation apparently overestimates the quantity of incoming solid particles over the top of the first hill. This could be related to a potential error in the experimental data, due the important difference between the X D 0 profile in Figure 5.6 and the X D 1:5 one in Figure 5.7. Table 5.5: Peak volume fraction and vertical location of the peak around the 3H hill from the simulation results.

In Figure 5.8, the concentration profiles around the two hills spaced 8H apart are plotted. The agreement between simulation results and experiments are reasonable. From Table 5.6, we note that the overall particle concentration presents a continuous decrease compared to the one around the upstream hill. Similar to Figure 5.6, the simulation apparently underestimates the quantity of solid particles passing between the two hills according to the experimental results.

We note that the global evolution of concentration peaks is similar between the 3H and 8H cases by comparing Tables 5.5 and 5.6. This implies a weak influence of the hill spacing on the overall evolution of mean particle concentration due to Gaussian hills. We have shown that the modelling of the complete solid particle saltation is valid. We will now focus on mechanisms that could trap or bring to deposition solid particles. To this aim, we examine in the next section the mean particle motion around the hills.

0 1 2 3 4 5 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 C/C max 0 1 0 1 0 1 Z / H X /

Particle velocity profiles

As the simulation results are validated using particle concentration profiles, in this part we focus on the kinematic aspects of the particle transport inside Z03 for the 3H and 8H cases, respectively. We recall that as solid particles are carried by the wind, a feedback effect from particles to the wind is accounted for in LES using two-way modelling (Section 3.1.3). This effect hinders any relative motion between wind and particles. Consequently, a reduction of wind speed induced by particles is in place in the velocity profiles presented in this section, in contrast with the ones presented in Section 4.2.3, where no particles are simulated.

Particle velocities presented in this section represent an Eulerian average of particle motion inside a unit volume V of size 0:1H L y 0:1H. The horizontal velocity profiles (u p ) are presented in Figures 5.9a As a general tendency, particle motion in the streamwise direction follow the prevailing wind speed across the TBL. It is the result of a long saltation layer development, starting from the sandbox located at 100H upstream of the first profiles given in Figures 5.9 and 5.10. History effects of the upstream hill are totally unremarkable in the profiles. Before the double hills, due to the absence of mean flow in the vertical direction , the mean particle velocity is much weaker, as particles in saltation keeps a back-and-forth movement, in the vertical direction, due to gravity and rebound. The Gaussian hills later deviates both solid grains and fluid into an upward motion, which explains the uniformly positive value of w P at the hill crests.

As for the effects of particles on the flow, not surprisingly, fluid velocity suffers from a slight deceleration compared to the case without particles. The vertical fluid motion sees its absolute value suppressed, as illustrated notably by the X=H D 1:5 profile in Figure 5.9b as well as X=H D 4 in Figure 5.10b, between grey (without particles) and black lines (with particles). On the whole, this feedback effect is weak on the mean flow, due to the diluted particle concentration over the double hills. Indeed, even the particle concentration peaks in Z03, given in Tables 5.5

and 5.6, border on the lower end of the range in which solid particles affects turbulence structures in a sensible way, according to [START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF].

Figure 5.9 illustrates the evolution of particle velocity side-by-side with the streamwise wind speed at the same locations over the double hills with 3H spacing. At X=H D 4, the mean particle motion almost follows the fluid motion across the depth of the TBL until a small lag begins to form from Z=H D 4. This could be related to the effect of gravity that ultimately forbids high-energy solid particles from catching up with the flow in the upper part of the TBL.

At the same time, the vertical particle velocity is substantially small compared to the flow. The presence of the Gaussian hill is already felt in the wind field at X=H D 4, which is not the case for solid particles. Near the ground a higher percentage of solid particles exhibit an upward motion. This near-wall, positive w p is in accordance with experimental observations that particles impact the ground with a smaller angle and rebound with a higher angle over a flat rough surface. As a result, the vertical velocity of the particle not only reverses its sign but sees its magnitude increase.

At the first hill top at X=H D 1:5, the streamwise particle velocity presents an inflection point at around Z=H D 3, above which u p increases steadily with height, similar to the X=H D 4 profile before the hill. The plateau between Z=H D 1 and 3 may be due to the blocking effect of the windward side, which redirect solid particles towards upper heights. This mean upward motion of particles is illustrated by the strong value of w p , shown in Figure 5.9b, whose peak (Z=H D 3) even slightly exceeds the wind speed up over the hill crest.

Over the second hill between X=H D 0 and 1:5, u p remains out-of-phase with the flow without notable variation in the vertical direction for Z=H > 1. The high peak of w p persists and move upwards from Z=H D 4:5 to 5 between X=H D 0 and 1:5. Inside the RZ at X=H D 0, we note the weak particle velocity in both horizontal and vertical directions. -5

-4 -3 -2 -1 0 1 2 3 4 5 6 U/U ∞ 0 1 0 1 0 1 0 1 0 1 Z / H X / H
(a) Longitudinal velocity profiles.

0 1 2 3 4 5 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 W/U ∞ 0 0.2 0 0.2 0 0.2 0 0.2 0 0.2 Z / H X / H
(b) Vertical velocity profiles.

Figure 5.9: Comparison between particle velocity (symbols), fluid velocity undisturbed by particles (grey line) and fluid velocity modified by two-way effects (black line) inside Z03 for the 3H case.

Behind the double hills at X=H D 4, an apparent discontinuity exists in u p slightly below the hill top at around Z=H D 0:8. A separate layer seems to appear below this level, populated with slow particles with sometimes reversed particle motion, due to the presence of backflows inside the RZ. Above Z=H D 0:8, grain motion maintains the same behaviour as X=H D 1:5. This particular behaviour, possibly related to particle trapping, is discussed further in Section 5.4.

For w p , whereas the peak remains high above Z=H D 5, its absolute amplitude decreases to half the value at X=H D 1:5. We note in particular the mean downward motion of grains for Z=H < 2 : gravity regains the solid particles previously bouncing up at the windward side of the double-hill. No discontinuity in w p profiles could be observed around Z=H D 0:8. This is probably due to the fact that u p and w p are not directly coupled in the particle motion equation described in Equation 3.9. On the other hand, it also implies that the trapping effect of the RZ acts mainly in the longitudinal direction.

In Figure 5.10, we compare the evolutions of particle and fluid velocities along the 8H Gaussian hills. The behaviour is overall similar to the 3H case: the plateau in u p develops after the first hill top and persists until reaching second hill; the strong positive peak of w p appears at the first hill top and moves upwards with decreasing amplitude All along the valley, the mean vertical motion of solid particles is directed towards the surface increasing the potential for particle trapping.

The larger spacing of 8H has bring modifications to the particle motion inside the RZ. The discontinuity in u p appears right after the first hill from X=H D 2 and continues at least to the foot of the second hill at X=H D 2. This suggests a continuous pattern of particle motion, probably of the trapped particles, across the whole valley. For w p , we note that from X=H D 0 to 2, particle motion in the vertical direction is in phase with the wind between Z=H D 0 and 2. This agreement breaks down over the second hill top, where a second peak in w p appears around Z=H D 2.

0 1 2 3 4 5 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 U/U ∞ 0 1 0 1 0 1 0 1 0 1 Z / H X / H
(a) Longitudinal velocity profiles.

0 1 2 3 4 5 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 W/U ∞ 0 0.2 0 0.2 0 0.2 0 0.2 0 0.2 Z / H X / H
(b) Vertical velocity profiles. In this section, the mean particle motion inside the saltation layer is investigated and compared to the wind motion. In the next section, we look into the particles on the ground susceptible of being re-entrained into the flow by the particle lift force on the ground.

Particle forces at the wall

In the previous sections, the development of a saltation layer over the Gaussian hills is studied through particle concentration and velocity profiles. Along the way, a fraction of solid particles in saltation lose much of their kinetic energy and fail to rebound from the ground after the impact. For these grains temporally immobilised at the wall, whether they can be once more entrained into the flow depends on the balance between a lift force exerted by the wind, modelled by our take-off model described in Section 3.2.4, and the sum of gravity (P ) and adhesion (F adh ), forces resisting grain movement.

In this section, the values of this particle lift force are given for the PC09 cases at U 1 D 11:2 m s 1 , and compared to the a priori estimations of Section 4.4. In addition to the effect of topography, feedback effects from solid particles to the wind have also impacts on these values.

A "lift-dominant" (LD) event can be defined, according to the instantaneous value of the particle lift, as F L .x; t/ > P C F adh . Under the influence of these events, According to the frequency and the intensity of these events, solid particles on the wall can gain enough momentum and return to the flow.

In Figure 5.11a, the ratio between F L and the sum of P and F adh in zone Z01, at the sandbox location is plotted. The mean intensity of the LD events corresponds thus to regions with F L =.P C F adh / > 1. We remind that the sandbox starts from X=H D 0 and ends at 20 in this figure. Within the sandbox, particle lift force rapidly drops to a third of its upstream value and does not recover until the end of the sandbox at X=H D 20. The cause of this fall can be largely attributed to an sudden accumulation of low-speed solid particles, yet to be accelerated into a more rapid motion. As shown in Section 5.3.1, particle concentration near the wall is at its maximum near the sandbox. The effect of the two-way coupling that models the particle feedback on the flow is thus substantial. As a result, the flow over the sandbox is greatly decelerated, resulting apparently a fall of energetic LD events.

We note that within the sandbox, the dense concentration of the grains near the wall maintains the value of the mean lift at a level slightly lower than the sum of gravity and adhesion, bordering the limit of the LD events. This reveals that the modelled aerodynamic entrainment by the flow reaches a saturated level at the particle source. This echoes the second hypothesis of [START_REF] Owen | Saltation of uniform grains in air[END_REF] on the aeolian transport of particles that both the concentration of particles engaging in the saltation and the flow carrying them adjust themselves so that the shear stress exerted by the wind on the ground is just sufficient to onset particle motion at the wall. Likewise, our results

show that the current take-off model, as described in Section 3.2.4, displays a self-balancing mechanism for solid particles emission at the sandbox.

As for the mean particle lift, a similar reduction in the time percentage of the LD events, T LD , also occurs in Z01 as shown in Figure 5.11b. Upstream of the sandbox where no solid particles are available in the simulation, the total duration of the LD events occupies on average 28% of the time, for a fictive particle located at the wall. By introducing solid particles into the flow, the duration of the LD events quickly decreases by approaching the sandbox, before stabilising around 17%, with little variation within the sandbox as the mean lift force. Starting from the end of the sandbox at X=H D 20, the LD duration recovers gradually and tends to 25%, a little less than its upstream level.

In Figure 5.12, evolutions of the particle lift force and the LD time percentages are plotted for the Z02 and Z03 areas. Compared to the sandbox location, the presence of airborne particles slightly decreases the particle lift force, but their effect is much weaker, especially on the occurrence of the LD events. This is due to the fact that, a larger number of solid particles begin to gain more kinetic energy and exhibits higher hop height in an incipient saltation layer. This results in a spreading of the overall particle concentration across the TBL and subsequently a weakening of concentration near the wall, as shown in Tables 5.5 and5.6. Figure 5.12 (a) reveals that, as the encounter with the Gaussian hill is imminent, the LD duration gradually drops to, 20% without, and 19% with grains, at X=H D 2:5. A similar decrease also appears in the particle lift. As the flow accelerates over the windward side of the hill, the LD events become more frequent before reaching a maximum around 30 to 32% slightly ahead of first hill crest. The mean lift force also reaches a maximum value before the hill crest.

Both the particle lift force and the LD duration falls again abruptly to nearly zero at X=H D 0:8, near the location at which the downhill slope reaches its maximum. The small bumps are likely caused by the intermittent backflows inside the RZ. Downstream of the mean RZ, the LD duration begins to gradually recover towards a new stabilised value of 25% at X=H D 8.

Over the double hills, as shown in (b) and (c) in Figure 5.12, The evolutions of the LD durations are alike and not particularly sensitive to the hill spacing. Whereas the LD duration increases by 20% on the windward side compared to its upstream value, it decreases more dramatically behind and between the Gaussian hills. In the 3H case, a second peak appears slightly ahead of the second hill crest, with an intensity much higher than the upstream peak. For the mean particle lift, the 8H case leads to a more balanced distribution of the LD intensity over the two hills. Moreover, the mean lift remains close to zero inside the mean RZ, and it is difficult for trapped particles to escape from this region with enough momentum. Behind this region, the particle lift force recovers towards its upstream value after X=H D 4 for the 3H case and after X=H D 6 for the 8H case.

In conclusion, we confirm that in the presence of solid particles, the windward side of the hills is highly susceptible of wind erosion due to the high value of mean lift and long duration of the LD events. In the valley between the double hills and on the lee sides, the mean particle lift force and the LD duration are both weak, forming potential areas of particle deposition. In the next section, individual particles are tracked in order to study the particle trapping inside these regions in which the particle lift force is weak compared to gravity and adhesion.

Near-wall particle transport and deposition

As pointed out in Section 3.2.1, the wall shear stress is of primary importance for the entrainment of solid particles. Its absolute value, related to the velocity gradient at the wall, directly measures the potential of particle transport near the wall. Therefore, classical saltation equations, as described in Section 5.2, commonly use the wall shear stress as a key parameter to evaluate the flow capacity to bring particles into motion. In our case, the wall shear stress is subject to both the terrain-induced modifications and the particle feedback. As depicted in Fig- Both the absolute value of w and the sign of wall velocity gradient are given by LES.

In Figure 5.13, the friction velocity in regions where the flow attaches to the wall is marked in red. In these areas, the particle transport is carried on along the prevailing wind direction of the upper flow. In comparison, the friction velocity related to backflow regions inside the RZ are marked in blue. Since the flow changes its direction near the wall, it is reasonable to assume that the backflow transports grains backwards towards the lee side of the hill. This "reverse transport", as studied in a conceptual way in [START_REF] Araújo | Numerical modeling of the wind flow over a transverse dune[END_REF] without considering the involvement of solid particles, could cause particle trapping behaviour inside the RZ.

In order to further explore the link between particle trapping and flow characteristics inside the RZ, we set out to identify regions of preferential deposition in the Z02 and Z03 areas. The net deposition of particles are represented using the packing density formed due to their accumulation, D , scaled by the closely-packed density at the sandbox, box defined in Section 5.1.1.

Evolution of this value during a period of T D 30ı 0 =u is shown in Figure 5.14 for the Z02 and Z03 zones. We remind that ı 0 is the depth of an undisturbed TBL at the origin in the figure.

The sign of D is used to indicate, respectively, regions of net erosion if D < 0, and ones of net deposition if D > 0.

The results show that the sporadic areas of net erosion are weak and more likely to occur on the windward face of the hill, whereas the net deposition of particles is more remarkable behind and between hill crests. Although the accumulations due to these deposited grains are in general one, even two orders of magnitude sparser than that at the sandbox, their value is substantial at locations where preferential depositions occur. Not surprisingly, these regions of high deposition coincide with areas with weak mean particle lift force (Figure 5.12). More interestingly, comparing Figures 5.13 and 5.14, we note that these regions are also located at the junction of positive and reverse transport. This could be explained by the fact that, instantaneously, frequent alterations take place in the direction of particle transport in these areas. For particles already trapped by the RZ, this back-and-forth motion may lead to a final deposition. Due to the weak lift force, they spend a longer time immobilised on the wall and contribute to the net accumulation in these regions.

Particle trapping by recirculation zones

As pointed out at the beginning of Chapter 4, one particularity of the RZ is its ability to delay the motion of air parcels. As a result, particles travelling through the RZ presumably spend a longer time isolated from the outer flow. Due to the weak flow velocity and lift force in this region, a larger amount of particles are deposited and accumulated in these regions. The results in Sections 5.3.3 and 5.3.4 provide some evidence in support of this claim. In this section, a direct study on particle trapping is conducted in which particles passing through the RZ are individually followed in order to study their interactions with the RZ.

Definition of particle trapping

We define particle trapping in this study as a particular state of particle motion in which the aerodynamic effects of the RZ become predominant in the determination of particle movement.

In practice, the following criteria are used in the study of particle trapping:

1. a particle is considered to be trapped by the RZ if it passes through the RZ and decelerates due to the backflows inside the RZ. For the single-hill case, the reference plane is located at the hill top at X=H D 0 in Fig- ure 5.15. For the PC09 cases, it is set on the first hill top within Z03 at X=H D 1:5 and X=H D 4, for the 3H and 8H cases, respectively. For each flow regime, 3000 samples have been collected with a sampling frequency of 600 Hz, 100 times higher than the particle relaxation frequency. The total duration is thus 4:5 s, which corresponds to 12 to 50 integral eddy turn-over times, depending on the flow regime.

In what follows, we provide statistics of particle trapping and deposition behind a single Gaussian hills at different Shields number.

Q I Q D Q T Q Out x z -h

Influence of the Shields number

The Shields number is a measure of flow capacity to set solid particles into motion. Flow with a stronger Shields number entrain more particles into the flow, and also supply them with a higher amount of kinetic energy. These particles in rapid motion, once trapped by the RZ, could escape with relative ease due to their high momentum and inertia. On the other hand, as the Shields number decreases, there exist a threshold regime under which no particle transport is sustainable, as mentioned in Section 5.2, Therefore, it is likely that there exists an optimum Shields number at which large number of incoming particles are trapped and end up deposited by the RZ.

This question is investigated over a range of Shields numbers from 0:02 to 0:41 by changing flow regimes with the same solid particle characteristics as in the PC09 experiments. The results are summarised in Table 5 In Figure 5.17 ). These peaks confirm the existence of optimum flow regimes of particle trapping and deposition between Sh D 0:03 and 0:07. For Shields numbers higher than 0:1, the particle deposition inside the RZ quickly becomes insubstantial, whereas the number of trapped particle decreases more gradually.

In Figure 5.18, ratios between trapped, and deposited particles, to the incoming flux are plotted, respectively. Among incoming particles, we note that the fractions of both trapped and deposited particles decrease monotonically with the Shields number. The number of trapped particles represents in general less than 1% of the total incoming particles for Sh > 0:07. Among these trapped particles, at least half are able to escape the RZ during the sampling period for all regimes studied here. This percentage of escaped particles raises quickly to over 80% for Sh > 0:07. In Figure 5.20, we note that the Shields number Sh D 0:07 seems to be the optimum flow regime for particle trapping and deposition between the two hills. In the range of Shields numbers studied here, the number of trapped particles is significantly higher in 8H case than in 3H case. This can be attributed to the larger RZ size in the 8H case. On the other hand, the absolute number of deposited particles is also higher at Sh D 0:12 in 3H case. In relative terms, particle deposition is more frequent among trapped particles in the 3H case, except at the lowest flow regime, as shown by the evolution of Q D =Q T in Figure 5.21 This high efficiency in particle deposition could be explained by the fact that the wall shear as well as the mean lift force at the second hill crest is weaker in the 3H case, as shown in Figure 5.12.

Q T / Q 0 , Q D / Q 0 Q I / Q 0 Sh Q T / Q 0 Q D / Q 0 Q I / Q 0
Q T / Q I , Q D / Q I Q D / Q T Sh Q T / Q I Q D / Q I Q D / Q T
Q T / Q 0 , Q D / Q 0 Q I / Q 0 Sh Q T / Q 0 Q D / Q 0 Q I / Q 0
Q T / Q I , Q D / Q I Q D / Q T Sh Q T / Q I Q D / Q I Q D / Q T

Conclusions

In this chapter, the main objective is to study the solid particle transport around the Gaussian hills, based on the knowledge gained from the aerodynamic studies of the flow around the Gaussian hills, presented in the previous chapter.

The chapter begins by a presentation of experimental configurations, with particular attention given to the solid particle feeding device. Numerical configurations are accordingly devised in order to simulate the PC09 cases with particle transport.

In the second section, we present a validation case of particle transport over a flat surface.

The simulation results on the saltation flux are compared with classical analytical formulas from literature. It is shown that, with particle characteristics studied here, the particle flux evolves asymptotically as a square of the friction velocity for Shields number larger than 0:1. This result is in accordance with the saltation models proposed by [START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF] and [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF].

Another point concerns the threshold particle transport. We note that our simulation results

give a threshold Shields number greater than the experiments of [START_REF] Creyssels | Saltating particles in a turbulent boundary layer: experiment and theory[END_REF] and [START_REF] Sørensen | An analytic model of wind-blown sand transport[END_REF]. A possible explanation of this discrepancy is our use of a sandbox emission instead of the commonly adopted sand feed injection. The latter method could artificially create splash movement at the sand bed, and thus enhance particle saltation.

In the third section, simulation results of the PC09 cases are presented. The simulated particle concentration profiles are in good agreement with the experimental results. We note in particular a non-negligible particle concentration inside the RZ formed by the Gaussian hills in both 3H and 8H cases. Then, the particle velocity is compared with the local wind velocity modified by the presence of solid particles. Globally, particle motions are in phase with with the mean flow motion. The adaptation of the particles to the presence of Gaussian hills is more late-occurring than of the TBL flow. Over the hills, the mean particle motion experiences a deceleration in the longitudinal direction and an enhancement in the vertical direction. A large fraction of particles are reflected into the upper layer of the boundary-layer, due to their impact at the steep hill slope. As a result, the average saltation height nearly doubles before and after the hill, as shown by the particle concentration profiles.

From particle velocity profiles, it is revealed that inside the large RZ formed between Gaussian hills and near the wall, there is a distinct layer filled with slow particles with particular motion patterns. Their behaviour is reflected in the apparent discontinuity of the longitudinal particle velocity profiles inside the RZ. This leads us to claim that these particles, under the influence of the RZ, participate in a reverse transport [START_REF] Iversen | The effect of surface slope on saltation threshold[END_REF]] separated from the mean particle transport.

Next, an evaluation of particle forces near the wall under the effects of both particles and Gaussian hills is conducted. This leads to an identification of regions in which particle lift-up is less likely to occur, thus revealing potential regions of particle deposition. On the other hand, regions with intense, frequent particle lift-up events are related to zones susceptible of wind erosion. The result is compared with a priori estimations presented at the end of the last chapter.

It is found in particular that at the sandbox location, the particle emission , modelled by the take-off model in our simulations, self-regulates to the extent that the mean particle lift force is equal to the sum of gravity and adhesion.

As an important link between flow and solid particle motion, the wall shear stress is highlighted and studied by considering its intensity and direction around the Gaussian hills. We find that the backflow inside the RZ generates a substantial mean shear opposite to the prevailing wind direction. Later, regions with high particle deposition are identified from the simulation results. It is found that these deposited particles generally lie within the backflow regions behind the hill, and they strongly coincide with the regions of weak particle lift-up events.

We thus set out to characterise the particles trapped and involved in a reverse transport by the backflows in the RZ. Specific criteria of particle trapping and deposition inside RZ are proposed.

Based on the simulation data, a large quantity of solid particles are tracked in order to investigate their potential encounter and interaction with the backflow inside the RZ. First, the influence of the Shields number on particle trapping is examined over a single Gaussian hill. For Sh > 0:07, the percentages of trapped, and deposited, particles represent less than 1%, and 0:1%, of the total incoming particles, respectively. Next, in the PC09 cases, the valley between the two hill crests are shown to maintain the fluxes of particle trapping and deposition at large Shields numbers. In particular, the 8H case captures more particles inside the RZ than the 3H case, whereas particles are more likely to deposit shortly after the second hill in the 3H case.

Chapter 6

Conclusions

This thesis presents a numerical study of the physical problem of soil erosion by wind. By using LES and modelling particle emission, rebound and feedback to the wind, we study particle transport inside a TBL in the presence of obstacles, to improve our knowledge on wind erosion and particle deposition over complex terrains. The aeolian process of wind erosion is studied through the schematic representation of a set of interactions between wind, topography and solid particles. The schematic topography is chosen by accounting for the remarks from [START_REF] Oke | Street design and urban canopy layer climate[END_REF] for urban canopy. We adapt such a reasoning to solid particles and desert dunes.

The first part of the thesis is devoted to the general study of a turbulent boundary-layer in the presence of transverse hills. A good agreement is found in the computed mean and RMS velocities compared to the experiments of [START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF]. Later, the simulations results are compared and analysed in comparison with the PC09 experiments [START_REF] Simoëns | Rapport ANR PEDO-COTESOF[END_REF]. The flow between and behind the Gaussian hills is dominated by large recirculation zones. In the 3H case, the flow characteristics on the whole is similar to the flow behind a single hill, and thus to the skimming flow in general, based on the terminology of [START_REF] Oke | Street design and urban canopy layer climate[END_REF]. In the 8H case, strong interactions take place between the shear layer forming from the first hill top and the windward side of the second hill. The strong vertical mixing across the valley in-between exhibits flow characteristics similar to that of the wake flow (again using the terminology of Oke). Influences of the hill slope, Reynolds number and wall condition on the occurrence and characteristics of the lee-side flow separation over a 2D hill are investigated. Furthermore, in the evaluation of particle lift at the ground, the windward side of the hill is identified as an area susceptible of wind erosion, while on the lee side and inside the mean RZ, particle deposition is more likely to occur and persist for a longer time.

In the second part, mass flux inside a TBL carrying particles emitted from a particle source is compared to analytical relationships in literature and the PC09 experiments. Asymptotically, the saltation flux is proportional to the square of the friction velocity for Shields numbers larger than 0:1. The threshold saltation regime predicted by the LES lies between 0:111 and 0.188 m s -1 , in accordance with the value given by the analytical formula of [START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF] and the PC09 experimental value. At the sandbox, we find that the concentration of particles taking off and the local flow carrying them adjust themselves so that the mean particle lift exerted by the wind on the ground grains is just sufficient to initiate their motion, in support of the second hypothesis of [START_REF] Owen | Saltation of uniform grains in air[END_REF] concerning particle saltation. Inside the saltation layer developed from the sandbox and over the Gaussian hills, the predicted particle concentration is in overall agreement with the experiments. Then, the particle velocity profiles reveal a particular pattern of particle motion inside the RZ, distinct from the one within the outer flow. The relevant particles are decelerated by the backflows and engage in a reverse transport towards the upstream, lee side of the hill. In accordance, a mapping of particle accumulation on the wall reveals that preferential particle deposition occurs in areas where intense and frequent lift-dominant events occur, and often at the junction between regions with mean and reverse particle transport behaviours. Last but not least, a study of particle trapping inside the recirculation zones is conducted. The fluxes of the particle trapping and deposition inside the RZ are compared to the total flux of particles incoming from an upstream reference location. Their percentages are shown to be of the order of a thousandth of the incoming flux. An optimum flow regime for particle trapping and deposition is found between Sh D 0:03 and 0:07. From Sh D 0:1, particle trapping and deposition decrease with increasing Shields numbers. For the successive-Gaussian cases, the larger canyon in the 8H case is able to attract a larger amount of particles into its RZ compared to both the 3H and the single-hill cases. Nevertheless, the amount of deposited particles is similar between the 3H and 8H case.

In view of these conclusions, further work could contribute to the following points :

1. An implementation of inter-particle interactions, e.g., the inter-particle collision model developed in the thesis of Vinkovic [2005], will allow us to better characterise particle behaviours at emission, deposition, and those related to the near-wall, closely packed particles in general. A modelling of the rotational motion of the particle both at the surface ("rolling") and in-air is necessary, for example, to represent the avalanche of grains on the lee side of the dune, crucial to the formation and migration of sand dunes, 2. In the current thesis, the splash process is not considered in the simulations, as it is not necessary in regard to the PC09 experimental set-up. To better compare with saltation flux measurements from wind-tunnel over sand beds and in-situ experiments in deserts, another numerical configuration needs to be devised, in which the surface is fully covered by sand grains, ready to be emitted into the air by impacting grains. In this case, the splash process needs to be considered in order to simulate a saltation layer in equilibrium with the flow as in the work of [START_REF] Dupont | Modeling saltation intermittency[END_REF].

3. For the future study on the morphology of sand dunes, the immersed boundary-method needs to be applied to the problem of solid particle transport. The implementation of this method, coupled with the deposition and re-emission of particles predicted by the LES, can take into account the two-way interaction between particle and topography and further improve our understanding on the problem of wind erosion in general.

As opposed to the log-region, the roughness sublayer present above rough surfaces is still poorly understood due to the complex interaction between wakes developing behind roughness elements. To investigate the spatially averaged flow velocity in this region, a dataset has been collected from several direct numerical simulations and wind-tunnel experiments available in the literature. A generalised law-of-the-wall has been derived, applicable to a roughness sublayer present over regularly distributed roughness elements. The key roughness parameter of this new law is the effective height ε, which characterises the interaction between the roughness and the outer flow in a temporally and spatially averaged sense. A morphometric study reveals that ε is closely related to a new roughness density parameter, λ 2 , that accounts for the roughness element shape and the inter-element spacing. This allows ε to be a universal parameter on roughness characterisation. The derived values of the classical roughness length z 0 of the log-law compare well with previous experimental data and geometrical model predictions. Finally, the main properties of the roughness sublayer such as its height are discussed using the geometrical and the roughness parameters proposed in the study.

Glossary

Introduction

Widely present in nature and in industrial applications, turbulent rough-wall flows have received continuous interest since the twentieth century. [START_REF] Nikuradse | Strömungswiderstand in rauhen Rohren[END_REF][START_REF] Perry | Rough wall turbulent boundary layers[END_REF][START_REF] Raupach | Rough-wall turbulent boundary layers[END_REF][START_REF] Jiménez | Turbulent flows over rough walls[END_REF] Rough-wall flows are more complex than smooth-wall flows because of the complex interactions between the surface roughness elements and the outer flow. Surface roughness elements are known to increase the surface drag (with some exceptions such as the riblets), and to induce turbulent structures that may modify the outer-layer flow. [START_REF] Krogstad | Structure of turbulent boundary layers on smooth and rough walls[END_REF][START_REF] Orlandi | Direct numerical simulation of turbulent channel flow with wall velocity disturbances[END_REF][START_REF] Monty | A high Reynolds number turbulent boundary layer with regular 'Braille-Type' roughness[END_REF] Far away from the rough wall, wind tunnel experiments [START_REF] Raupach | A wind-tunnel study of turbulent flow close to regularly arrayed rough surfaces[END_REF] and in-situ measurements in rural [START_REF] Marticorena | Factors controlling threshold friction velocity in semiarid and arid areas of the United States[END_REF] and urban environments [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF][START_REF] Cheng | Near wall flow over urban-like roughness[END_REF] (U) exhibits a logarithmic form. Depending on the scaling parameter, two expressions of this logarithmic profile, also referred to as the log-law, have been proposed in the past.

-The first expression uses the roughness length z 0 of the surface, and writes as:

U + = κ -1 ln(z/z 0 ) , (1)
where U + = U/u * , u * is the friction velocity, κ the von Kármán constant (VKC), and z the wall-normal coordinate. Here, the origin of z is defined at the displacement height z d , located within the roughness height h, as illustrated in Figure 1. The roughness length z 0 is the height at which the velocity would be zero according to the log-law. In this way, z 0 characterises the absorption of the flow momentum by the underlying rough wall. At high Reynolds number, the flow dependence of z 0 becomes negligible and the resulting rough-wall flow is considered 'fully rough' as opposed to 'transitionally rough' at lower Reynolds number. [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] Only in this case, is z 0 a unique function of surface geometrical properties. However, it is difficult to estimate z 0 of an arbitrary rough surface from its geometrical properties such as h. This is especially true for realistic roughness. Contrary to what its name suggests, z 0 is more of a mathematical integration constant than a topographical scale of the roughness. [START_REF] Hanna | Wind flow and vapor cloud dispersion at industrial and urban sites[END_REF] Nevertheless, several empirical models of z 0 have been proposed in the literature as a function of geometrical parameters. [START_REF] Raupach | Drag and drag partition on rough surfaces[END_REF][START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF][START_REF] Kastner-Klein | Mean flow and turbulence characteristics in an urban roughness sublayer[END_REF][START_REF] Millward-Hopkins | Estimating aerodynamic parameters of urbanlike surfaces with heterogeneous building heights[END_REF] However, the empirical constants employed in these models lack from physical interpretations, and are generally tuned by a limited number of measurements. This severely limits the reliability of these models and may explain the large scatter in the predicted values of z 0 .[17,18] -The second expression of the log-law is an extension of the log-law over idealisedsmooth wall. The log-law over smooth wall is theoretically based on the existence of an intermediate region, named the inertial sublayer (IS), where the only relevant length scale is z, instead of a viscous length scale or an external length scale. Hence, the non-dimensionalised velocity gradient, dU + / dz + , can be expressed as a function of z + in wall units (z + = zu * /ν, where ν is the kinematic viscosity of the flow). After integration, U + follows a logarithmic profile U + = κ -1 ln z + + B, where B is an integration constant. Based on several experiments such as in [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF] and [START_REF] Hama | Boundary layer characteristics for smooth and rough surfaces[END_REF], it has been claimed that the presence of roughness elements modifies the mean flow through a downward shift U + :

U + = κ -1 ln z + + B -U + , (2)
where U + is named the roughness function. Below the log region, i.e. in the near-wall region, various formulations of a unified velocity profile have been proposed for smooth-wall flow, such as the formula of Spalding [START_REF] Spalding | A single formula for the "law of the wall[END_REF], where a viscous sublayer holds in the vicinity of the wall, a linear velocity profile applies up to z + ∼ 5, and the log profile is considered valid for z + > 30. On the other hand, the lower limit of the logarithmic region remains unknown to rough-wall flow, and no analytical velocity profile formulations have been proposed to our knowledge in the near-wall region.

The earliest works on rough-wall flow focused on the global effect of the wall roughness elements, without probing into the near-wall region: Nikuradse [START_REF] Nikuradse | Strömungswiderstand in rauhen Rohren[END_REF] studied experimentally the friction factor related to the pressure drop in the sand-grain roughened pipe; Moody [START_REF] Moody | Friction factors for pipe flow[END_REF] tabulated this factor as a function of the Reynolds number, covering both aerodynamically smooth and rough surfaces. Later, Townsend [START_REF] Townsend | The structure of turbulent shear flow[END_REF] proposed a wall-similarity hypothesis, where the roughness disturbances are confined to a zone near the wall, while outside this region, the roughness effect is only felt through a roughness-modified wall shear stress. Although this hypothesis provides some support to the validity of the log-law over a rough wall, it does not clarify the near-wall region.

In the context of drag reduction, Bechert and Bartenwerfer [START_REF] Bechert | The viscous flow on surfaces with longitudinal ribs[END_REF] studied the alterations of the mean longitudinal flow produced by surfaces with longitudinal ribs. Bechert and Bartenwerfer [25] defined the protrusion height by arguing that the velocity profile appears as if it originated from an equivalent plane wall located at a distance below the riblet tips. Luchini et al. [26] defined a similar protrusion height in the case of cross-flow over grooved surfaces. By analysing the Stokes flow of a fluid across the grooves, Luchini et al. [START_REF] Luchini | Resistance of a grooved surface to parallel flow and cross-flow[END_REF] characterised the differential effect of the wall on parallel flow and cross-flow by defining the difference of the two (parallel and cross-flow) protrusion heights. The authors' analysis lead to conclude that what makes the corrugated wall effective in retarding cross-flow, which presumably decreases the overall turbulent drag to the flow, is that the virtual plane wall seen by cross-flow is located deeper into the fluid than the one seen by parallel flow. Later, by shear-stress measurements over conventional riblet configurations, Bechert et al. [START_REF] Bechert | Experiments on drag-reducing surfaces and their optimization with an adjustable geometry[END_REF] supported the theoretical model proposed by Luchini [START_REF] Luchini | Effects of riblets on the growth of laminar and turbulent boundary layers[END_REF] based on the assumption that riblets impede the fluctuating turbulent cross-flow near the wall, and in this way reduce momentum transfer and shear stress. By direct numerical simulations (DNS) of flow over riblets with increasing sizes, García-Mayoral and Jiménez [START_REF] García-Mayoral | Hydrodynamic stability and breakdown of the viscous regime over riblets[END_REF] studied the breakdown of the viscous regime leading to drag increase. The authors found that the breakdown is caused by long spanwise rollers that develop from a two-dimensional (2D) Kelvin-Helmholtz-like instability of the mean streamwise flow. By a simplified linear instability model, García-Mayoral and Jiménez [START_REF] García-Mayoral | Hydrodynamic stability and breakdown of the viscous regime over riblets[END_REF] explained that the breakdown spacing scales with the Downloaded by [Gang Huang] at 11:23 01 March 2016 square root of the groove cross-section. As defined in these studies, [START_REF] Bechert | The viscous flow on surfaces with longitudinal ribs[END_REF][START_REF] Luchini | Effects of riblets on the growth of laminar and turbulent boundary layers[END_REF][START_REF] García-Mayoral | Hydrodynamic stability and breakdown of the viscous regime over riblets[END_REF] the ratio of the protrusion height to the period of corrugations is a purely geometrical parameter depending only on the shape of the wall corrugations and neither on their size nor on the actual speed of the driving fluid stream. The idea introduced is that the corrugated wall is equivalent to a plane wall located at a distance below the riblet tips.

With the development of measurement techniques and numerical methods, the nearwall region of a rough-wall flow in general has received growing interest. [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF][START_REF] Orlandi | On the generation of turbulent wall friction[END_REF][START_REF] Leonardi | Direct numerical simulations of turbulent channel flow with transverse square bars on one wall[END_REF][START_REF] Florens | Defining the roughness sublayer and its turbulence statistics[END_REF] The notion of the roughness sublayer (RS) was first proposed by [START_REF] Raupach | A wind-tunnel study of turbulent flow close to regularly arrayed rough surfaces[END_REF], which corresponds to the upper limit of a 'roughness wake effect' region where the spatially averaged velocity departures from the log-law. More generally, Jiménez [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] proposed that the RS should be considered as 'the region with direct roughness effect' . Since the near-wall flow is subject to the influence of individual roughness elements, it appears necessary to look at spatially averaged velocity profile, providing an alternative to encompass the effect of individual roughness elements. [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF] If the RS is seen as a mixing layer between roughness elements and upper boundary-layer, this is the only way to characterise such a flow part with a unique parameter such as the roughness density. Through spatial averaging, an equivalent representation of the RS emerges with spatially homogeneous properties. Furthermore, by delimiting the RS to the lower limit of the logarithmic region z r , the RS can be seen as a roughness inner layer in which the roughness effect modifies directly the mean flow dynamics. [START_REF] Florens | Defining the roughness sublayer and its turbulence statistics[END_REF] In this way, RS bears resemblance to the viscous wall region (z + < 50 according to Pope [START_REF] Pope | Turbulent flows[END_REF]) over a smooth wall in which the direct effect of molecular viscosity is dominant. The RS height is still a research issue. A multitude of criteria have been proposed in the literature to determine the RS height. In [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF], the RS height is deduced by comparing the defection of the velocity profiles between smooth and rough-wall flows. Ashrafian et al. [START_REF] Ashrafian | DNS of turbulent flow in a rod-roughened channel[END_REF] defined the RS height as the upper limit of spatial inhomogeneities in the Reynolds stress profiles.

The objective of this paper is to propose a universal profile of the mean wind velocity with a roughness scaling depending on the flow regime and on a roughness density parameter that could not be overlooked. To that purpose, DNS and wind-tunnel experiments of near-surface wind flows over various roughness element densities and configurations available in the literature are used. The novelty of the roughness scaling is that it is based on the effective height ε, which characterises the depth of the mixing layer in the upper roughness element layer participating directly to the turbulent exchanges with the outer flow (Figure 1).

The paper is structured as follows. First, the roughness density that characterises the surface geometrical properties, the roughness scaling based on ε, and the definition of the non-dimensional velocity are introduced in Section 2, followed by a description of the wind flow data-set collected from the literature. In Section 3, the generalised law-of-the-wall is deduced from a wind profile matching. The dependence of the effective height on surface geometrical properties is studied subsequently. From the obtained law-of-the-wall, we also deduce an estimation of the classical but here modified (roughness density-dependent) roughness length z 0 , and compare it to the experimental data of Cheng et al. [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF] and to the predictions of the geometrical model of Kastner-Klein and Rotach [START_REF] Kastner-Klein | Mean flow and turbulence characteristics in an urban roughness sublayer[END_REF] and MacDonald et al. [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF]. In Section 4, the flow regimes and the characteristics of the RS are discussed using the newly defined roughness parameters. Conclusions are given in Section 5. Downloaded by [Gang Huang] at 11:23 01 March 2016

Key parameters and data-set

Roughness density

The flow dynamics within the RS depend on the roughness element configuration. The roughness density of a surface is one of the most influential parameters characterising the friction induced by a rough surface on the flow. [START_REF] Lee | Turbulent boundary layers over rod-and cube-roughened walls[END_REF][START_REF] Ahn | Statistics of the turbulent boundary layers over 3D cube-roughened walls[END_REF] Several definitions exist for the roughness density, depending on the representative area of the roughness elements. For a surface area S covered by n roughness elements of length l, breadth b, and height h, the frontal area density is defined as λ f = nbh/S and the plane area density as λ p = nbl/S. By definition, λ p does not depend on the roughness height h, as opposed to λ f . Since the dataset used in this study is mainly composed of roughness elements with aspect ratio l/h close to unity, we hereafter assume λ p ∼ λ f .

Several studies showed that the roughness shape such as the windward slope of the roughness elements plays also a meaningful role in the roughness effects of the surface on the flow. [START_REF] Orlandi | Direct numerical simulation of three-dimensional turbulent rough channels: parameterization and flow physics[END_REF] Elements with sharp edges appear more effective in creating fully rough flow than rounded sand grains. [START_REF] Snyder | The critical Reynolds number for rough-wall boundary layers[END_REF] Since the shape factor is not directly accounted for neither in λ f nor in λ p , we define here a new roughness density λ by pondering λ p (or λ f ) with the ratio between the average ( h) and the maximum (h max ) heights of the roughness elements:

λ = λ p h h max . (3)
Hence, for square-shaped roughness elements, λ equals to λ p , while for cases with roundedshape elements, λ is smaller than λ p . It will be further shown in Section 3.2 that it is necessary to consider, in our work, the anisotropy of roughness elements to completely account for the geometric properties.

Flow regimes and effective height

As classical parameter used in the log-law, the displacement height z d characterises the roughness level of the wall. A first approach to determine z d consists in fitting the velocity profile to the log-law (Equation ( 1) or Equation ( 2)) by considering z d as an offset to the origin of z in the log-law. Although this fitting technique has proven to be accurate for aerodynamically smooth surfaces, it remains uncertain for rough-wall cases as neither z 0 nor U + are known a priori. [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF] Moreover, fitting the velocity profile blindly to the loglaw without knowing its range of validity beforehand could give aberrant values of flow parameters such as z 0 . With a more physical approach, Jackson [START_REF] Jackson | On the displacement height in the logarithmic velocity profile[END_REF] proposed to link z d to the height at which the total drag acts on the roughness elements. For λ p = 25%, Cheng et al. [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF] could not obtain a reasonably good fit of the log-law profiles without changing the roughness parameters such as the friction velocity. According to the authors, [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF] this suggests that although Jackson's theory may work well at low density, it breaks down when the flow begins to skim the roughness. After assuming κ = 0.41 (Equation ( 1)), Leonardi et al. [40] compared the log-law fitting of the mean velocity data to the approach of Jackson [START_REF] Jackson | On the displacement height in the logarithmic velocity profile[END_REF]. The authors [START_REF] Leonardi | Properties of d-and k-type roughness in a turbulent channel flow[END_REF] found good agreement for big values of λ p , while larger differences are observed when λ p is smaller than 20%.

Here, we argue that the inherent problem of z d is due to its definition that neglects the physics underlying the exchanges between the inside of roughness elements and the layer just above the roughness top. In fact, the value of z d relates to the bottom of the canopy layer not involved in the exchange process that is at the origin of the RS. In comparison, the effective height ε characterises the level of exposure of roughness elements to the incoming flow, as illustrated in Figure 1. It is this exposure that leads to the creation of a mixing-layer in which the interaction between the roughness and the outer flow takes place. We, therefore, claim that ε is a more appropriate roughness length scale that seems to characterise all kinds of flow regimes in terms of the exchange between the interior of roughness elements and the outer flow. Unlike in [START_REF] Perry | Rough wall turbulent boundary layers[END_REF], we define ε by using the same dependence to the roughness density as in [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF], allowing to note three different classes of flow regimes that are described in the following paragraph. To highlight the difference between our analysis with ε to the usual one with z d , we define z d = h -ε in Figure 1 that is not directly related to z d used classically as a coefficient of the log-law.

Depending on the roughness configuration, several studies have shown that the flow can be classified following three regimes: isolated, wake, and skimming flow regimes. [START_REF] Leonardi | Direct numerical simulations of turbulent channel flow with transverse square bars on one wall[END_REF][START_REF] Simoëns | The flow across a street canyon of variable width-Part 1: kinematic description[END_REF][START_REF] Dupont | Modeling aeolian erosion in presence of vegetation[END_REF] For sparse roughness elements (isolated flow regime), the interaction between the flow and the roughness elements is at its full extent since individual elements are completely exposed to the flow. As the roughness density grows (wake flow regime), the momentum exchange between the roughness and the outer flow diminishes due to the sheltering effect between the elements. [START_REF] Raupach | Drag and drag partition on rough surfaces[END_REF] For extreme roughness density (skimming flow regime), the interaction between the flow and the roughness elements vanishes as a new smooth wall emerges. Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF] delimited these three regimes by λ p , based on the wind tunnel experiment of [START_REF] Hussain | A wind tunnel study of the mean pressure forces acting on large groups of low-rise buildings[END_REF]. This led them to propose an evolution of z d as a function of λ p :

-Isolated flow with λ p < 0.15 and z d < 0.2h, in which individual roughness elements are completely exposed to the flow; -Wake flow with 0.15 < λ p < 0.35 and 0.2h < z d < 0.75h, marked by strong wake interferences between elements; -Skimming flow with λ p > 0.35 and 0.75h < z d < h, in which an isolated vortex is created in the roughness cavity with little interaction with the outer flow. We will modify this classification with the parameters λ 2 (defined in Section 3.2) and ε, allowing us to provide a new roughness density flow regimes completely defined by ε and λ 2 alone.

Definitions of the non-dimensional wall distance and velocity

To characterise the mean flow inside the RS, a roughness scaling is applied to the vertical coordinate z using the roughness length scale ε

z -= z/ε . ( 4
)
For the streamwise velocity, the following non-dimensionalisation is used: The slip velocity, U + w , corresponds to the spatially averaged streamwise velocity in wall units at the top of roughness elements.

U -= U + -U + w , ( 5
The value of U + w is generally non-zero depending on the roughness configuration. We chose to subtract U + w to U + because only the region starting from the top of roughness elements is of interest in this study. The scaled velocity U -has been previously used in [START_REF] Orlandi | Direct numerical simulation of three-dimensional turbulent rough channels: parameterization and flow physics[END_REF][START_REF] Orlandi | DNS of turbulent channel flows with two-and three-dimensional roughness[END_REF]. The authors defined their log-law using U -as a function of z + , while systematically considering the displacement height at the top of roughness elements (z d = h), which is different from our approach.

Data-set

Mean streamwise wind velocity profiles were collected from recent DNS and experimental results available in the literature. This data-set includes both channel [START_REF] Leonardi | Direct numerical simulations of turbulent channel flow with transverse square bars on one wall[END_REF][START_REF] Leonardi | Properties of d-and k-type roughness in a turbulent channel flow[END_REF][START_REF] Orlandi | DNS of turbulent channel flows with two-and three-dimensional roughness[END_REF][START_REF] Leonardi | Channel flow over large cube roughness: a direct numerical simulation study[END_REF] and boundary-layer [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF][START_REF] Lee | Turbulent boundary layers over rod-and cube-roughened walls[END_REF][START_REF] Ahn | Statistics of the turbulent boundary layers over 3D cube-roughened walls[END_REF] flows. Although this data-set includes mainly 2D and threedimensional (3D) roughness elements with cubical section, some roughness cases with circular and triangular section are also present, [START_REF] Orlandi | Direct numerical simulation of three-dimensional turbulent rough channels: parameterization and flow physics[END_REF] which justifies the need to account for the element edge shape in the roughness density formulation (Equation ( 3)). For 2D elements, both longitudinal and transverse configurations are considered. [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF][START_REF] Orlandi | Turbulent channel flow with either transverse or longitudinal roughness elements on one wall[END_REF] The plan area density λ p varies from 4% to 50%, covering all three flow regimes according to the ranges proposed by Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF]. To remove periodic fluctuations caused by individual roughness elements, the collected velocity profiles have been averaged both temporally and spatially in the original papers. The bulk Reynolds number Re b , based on the bulk velocity and the external length scale of the flow, varies from 7 × 10 3 [START_REF] Lee | Turbulent boundary layers over rod-and cube-roughened walls[END_REF] to 5 × 10 4 . [START_REF] Leonardi | Channel flow over large cube roughness: a direct numerical simulation study[END_REF] Both transitionally and fully rough regimes are covered. Table 1 gives further details on the data-set.

Figure 2 illustrates the wind velocity profiles of our data-set, classified following the three flow regimes defined according to the ranges proposed by Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF]. At this stage of the study, ε is unknown and was set to h. The following observations can be made from Figure 2: -For all regimes, the departure from the log-law is noticeable in the near-wall region. The logarithmic region is more recognisable in the isolated flow regime than in the other two regimes. An RS can be seen in all velocity profiles. -For the wake and isolated flow regimes, the velocity profiles exhibit an inflection point at the top of the roughness elements. -For the skimming flow regime, the departure from the log-law is the most pronounced and the size of the RS (z r /h) appears proportional to the roughness density. -Despite some similarities between velocity profiles of each flow regime, the profiles do not collapse as ε has been assumed constant (ε = h) at this stage for all roughness cases.

Results

The law-of-the-wall following the three flow regimes

We remind that, in the non-dimensional velocity profiles plotted in Figure 2, ε was not known and has been replaced by h. Using the non-dimensional velocity and wall distance defined in Section 2. the roughness top to the upper limit of the log region, for each flow regime (see Figure 3), by considering ε as a roughness density-dependent parameter. We obtained the following velocity profiles in the log and RS regions.

In the log-region, the velocity profile follows the same log-law for the three flow regimes:

U -= κ -1 ln z -+ W, [START_REF] Orlandi | Direct numerical simulation of turbulent channel flow with wall velocity disturbances[END_REF] with W = 1.2 for the wake and isolated flow regimes, and W = 4.5 for the skimming flow regime. By combining Equation ( 1) and Equation ( 6), we can relate W to z 0 as: Within the RS, the deviation of the velocity profile from the log-law is modelled as an exponentially decreasing function of z -for both isolated and wake flow regimes, from a maximum W = 1.2 at z -= 1 (the top of roughness elements) to roughly 0 at z -= 2.2 (Figure 3(a) and 3(b)), leading to:

W = κ -1 ln(ε/z 0 ) -U + w . ( 7
U -= κ -1 ln z -+ W -W exp( z --1 ) -4 . ( 8
)
Note that the case K1 from [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF] exhibits the most significant departure from the other profiles. In this case, we have W = 3.8 instead of 1.2. This discrepancy will be discussed in Section 4.3. Combining Equation [START_REF] Monty | A high Reynolds number turbulent boundary layer with regular 'Braille-Type' roughness[END_REF] with Equation ( 8) leads to

U + = κ -1 ln(z/z 0 ) -W exp( z/ε -1 ) -4 , ( 9
)
for z -ࣙ 1. Equation ( 9) constitutes a modified law-of-the-wall to Equation ( 1) using the roughness parameter ε, for wake and isolated flow regimes.

For the skimming flow regime, the following collapsed profile is obtained within the RS (Figure 3(c)):

U -= 4.8 ln z -. [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF] This equation resembles to the law-of-the-wall proposed by von Kármán [START_REF] Von Kármán | The analogy between fluid friction and heat transfer[END_REF] between z + = 5 and z + = 30 over a aerodynamically smooth wall.

Effective height

To study the dependence of ε on the roughness element properties, Figure 4 shows ε/h, obtained from the wind velocity profile matching as a function of the roughness density λ for 2D (Figure 4(b)) and 3D (Figure 4(a)) roughness elements. Although we do not share the same definition of ε or z d , the values of hz d obtained from the data of Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF] are also presented on the same figure for comparison. Qualitatively, the evolution of ε/h are in accordance with the one of hz d from [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF]. More importantly, Figure 4 suggests that ε/h is correlated with the flow regimes. Mostly notable within the wake flow regime, 2D roughness elements apparently give higher values of ε than 3D ones. This difference suggests that other geometrical parameters than λ may influence ε, such as the distance between roughness elements. We define below these new parameters that will replace λ as defined in Section 2.1.

For an identical roughness density λ, different arrangements of roughness elements are possible. These configurations are not equally sensible to the prevailing wind direction and also channel the wake flows in different ways between elements. Let us note L x as the streamwise spacing between individual roughness elements, and L y as the spanwise one. Ideally, a 2D roughness configuration is obtained by decreasing L x or L y between roughness elements, assumed here to be of the basic shape b = l (satisfied in all 3D roughness cases in our data-set), to the minimum value. For example, we have L y /L x = λ p and L y /L x = 1/λ p for transversal and longitudinal bars, respectively.

Since L x follows the prevailing wind direction along which most roughness wakes develop, increasing L x /L y while keeping unchanged λ may attenuate the roughness wake interference, and thus favour the transition from the wake to the isolated flow. This led us to hypothesise that the wind flow regimes could be delimited by a new roughness density λ 2 that depends on the ratio = L y /L x :

λ 2 = λ α , ( 11
)
where α is a constant to be determined empirically. The value of can be deduced for roughness configurations with regular distribution. For example, we have = 1 for uniformly distributed roughness, = 2 for aligned, and = 0.5 for staggered roughness configurations, respectively. These configurations are depicted in With this new roughness density, a good fit of the effective height as a function of λ 2 is obtained for both 2D and 3D elements with α = 0.2 (Figure 6). For λ 2 between 0 and 1, we obtained:

ε h = ⎧ ⎨ ⎩ 1 i f λ 2 < 0.09 0.035 λ -1.4 2 if 0.09 < λ 2 < 0.4 0.125 if λ 2 > 0.4 . ( 12
)
To determine ε from data, we have used a best fitting, as classically done to obtain z 0 in the log-region, of the RS. By definition, z -= z/ε, thus we have only one parameter, ε, to determine. Yet, as pointed out in the introduction, a best fitting is indeed subject to errors. Nevertheless, our objective was to obtain an analytical law for ε, as given by Equation ( 12). This equation can then be used to obtain ε for a new set of data instead of the best fitting, in the limit of validity and for the restrictions given in the discussion of Section 4.4.

Downloaded by [Gang

Roughness length

A new parameterisation of the roughness length (z 0 ) as a function of the effective height (ε) and the slip velocity (U + w ) can now be derived by rewriting Equation 7:

z 0 = ε exp(-κ (U + w + W )) . ( 13
)
To relate the roughness length to the surface properties, we assume that U + w depends only on surface roughness configurations at high roughness Reynolds number h + = hu * /ν. Perry et al. [START_REF] Perry | Rough wall turbulent boundary layers[END_REF] identified a near-wall velocity scale at the edge of the RS in rough-wall flow, and drew an analogy between this velocity scale and the velocity at the outer edge of the viscous sublayer over a smooth wall. Here, U + w bears resemblance to this near-wall velocity scale since U + w is defined at the top of the roughness elements, which separates the canopy layer and the RS. As argued by Perry et al. [START_REF] Perry | Rough wall turbulent boundary layers[END_REF], their near-wall velocity scale can be viewed as the inverse of a local drag coefficient u * /U r , which should be constant at large Reynolds number.

Figure 7 shows U + w as a function of λ 2 . Note that less data is available to define a parameterisation of U + w as this parameter is usually not reported in papers. [START_REF] Orlandi | Turbulent channel flow with either transverse or longitudinal roughness elements on one wall[END_REF] Generally speaking, U + w decreases with λ 2 , with larger values of U + w for 3D roughness elements than 2D ones. Several cases depart nevertheless from the average tendency shown in the figure. The K1 case shows a lower value of U + w that may be related to its smaller roughness Reynolds number (h + < 15). The parallel riblet case, O1, exhibits a high value of slip velocity, which is probably due to the drag-reduction effect of the riblet.

By excluding the K1 case and several cases with large roughness densities (λ 2 ࣙ 50%, O1 and OL4) , the following fitting parameterisations of U + w are obtained: -for 3D roughness cases (h + > 70): good agreement is found with the experimental data of Cheng et al. [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF]. The result does not compare well with the prediction of Kastner-Klein and Rotach [START_REF] Kastner-Klein | Mean flow and turbulence characteristics in an urban roughness sublayer[END_REF]. For isolated flow (λ 2 < 0.09), z 0 may be overestimated in our result. This may be due to a discrepancy of the parameterisation of U + w . In Figure 8(a), we show that the peak of z 0 for 2D roughness is about two times higher than for 3D roughness. This is supported by the experimental results collected and shown in [START_REF] Raupach | Rough-wall turbulent boundary layers[END_REF].

U + w = -
In Figure 8(b), we illustrate in particular the difference in roughness effects between aligned and staggered roughness configurations. The model of MacDonald et al. [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF] for uniform roughness configuration is plotted for comparison. Using λ 2 , the model predictions of MacDonald et al. [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF] and the experimental data do not collapse into a single curve, contrary to the prediction of our model. The discrepancy is especially large at small roughness densities for the aligned case (C4). For the wake flow regimes, our results are closer to the experimental cases of C1 (staggered) and C2 (aligned) from [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF] than the model of MacDonald et al. [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF].

Discussion

A generalised law-of-the-wall over rough surfaces has been obtained. In this section, we revisit the classification of flow regimes, and attempt to retrieve the main characteristics of the RS, such as z r and the velocity deviation W. Finally, the application range and the uncertainties of the proposed wall-law are discussed.

Flow regimes

We observed in Section 2.2 that the evolution of the effective height ratio ε/h is correlated with the flow regimes. Since a reasonable collapse of the wind velocity profiles has been obtained for each flow regime (Figure 6), it seems physically appealing and was shown above to be more general and needing less geometrical dependence to separate the three regimes by λ 2 and ε, instead of λ p and z d as in [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF]. The result of the velocity profile matching in Section 3.1 delimits the flow regimes as follows:

-Isolated flow: λ 2 < 0.09 and ε ∼ 1.

-Wake flow: 0.09 < λ 2 < 0.4 and 0.15 < ε < 1.

-Skimming flow: λ 2 > 0.4 and ε ∼ 0.15. This delimitation has the advantage of being applicable to both 2D and 3D roughness elements, whereas the formulation by Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF] could not. It takes into account the non-square shape and the distances between roughness elements through λ 2 , as opposed to the delimitation proposed by Grimmond and Oke [START_REF] Grimmond | Aerodynamic properties of urban areas derived from analysis of surface form[END_REF]. Furthermore, this parameterisation is supported by the physical behaviour in-between roughness elements and in the RS.

However, cases classified as skimming flow regime show a relatively large dispersion in the values of ε/h, with an average value of 0.15 (Figure 6). It appears to us that λ 2 > 0.4 is not a sufficient condition to define the skimming flow regime. For example, the parallel riblet case (noted as O1 in Figure 6) with λ 2 = 0.5 and the transversal square bar case (L1) with λ 2 = 0.43, have been both classified as skimming flow regime, while the transversal riblet case (O2) with λ 2 = 0.43 has been classified as a wake flow regime based on the velocity profile. This discrepancy may be due to the more favourable exchange with the outer Downloaded by [Gang Huang] at 11:23 01 March 2016 flow of the windward triangular shape of O2, and thus the establishment of a wake flow regime. This underlines the importance of roughness shape on the onset of skimming flow at large λ 2 . Further investigation is, therefore, needed to improve the parameterisation of the skimming flow regime.

Roughness sublayer height

The RS height z r informs on the lower limit of the log region (Figure 1). Using the generalised law-of-the-wall (Equations ( 8) and ( 10)), we find that z r is proportional to the effective height ε for the three wind flow regimes. For wake and isolated flow regimes, z r is about 2.2ε above the displacement height, a height at which the departure from the log-law is less than 1% of the maximum deviation (W). For the skimming flow regime, z r is about 7.4ε, which corresponds to the height where the wind profile within the RS equals that in the log region. Based on the evolution of ε with λ 2 in Figure 6, it is interesting to note that the minimum of z r falls within the range of wake flow around λ 2 = 0.3.

Compared to the mean roughness element height h, RS reaches about 0.75h for the skimming flow regime from the roughness element top, between 0.2h and 1.2h for the wake flow regime and approximately 1.2h for the isolated flow regime. This is consistent with the usual values reported in the literature [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF][START_REF] Florens | Defining the roughness sublayer and its turbulence statistics[END_REF][START_REF] Rotach | On the influence of the urban roughness sublayer on turbulence and dispersion[END_REF] where h < z r < 2h. Note that these previous studies set the roughness height as the upper limit of the spatial inhomogeneity, while here we use the spatially averaged velocity profiles to even out the disturbances caused by individual roughness elements. This could explain our lower estimates of z r .

In [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF], it is reported that the log-law can be extended to the roughness top for the spatially averaged velocity profiles. It implies that the height of the RS is zero in their cases. The discrepancy can be explained by the following arguments. First, we remark that the reported value of z d , determined by fitting the experimental data to the log-law in both IS and RS regions by the authors, results in a value of ε 50% lower than the one predicted by Equation ( 12), thus a lower value of z r than our prediction. Besides, with a roughness density of 25% (Table 1), z r is relatively low as stated earlier. The consequence is that it would have been necessary to examine the streamwise velocity within a distance of 0.3h above the roughness top in their experience. This is extremely difficult to achieve because the flow is largely chaotic at this wall distance. In our opinion, these factors combined could explain the oversight of the velocity departure from the log-law in the case of Cheng and Castro [11].

Fully rough vs. transitionally rough regimes

Here, we would like to highlight the fact that different parameters necessary to establish the present law, such as W, ε, are now well established for high enough Reynolds number. Yet, for transitionally rough regimes, the dependence on Reynolds number needs to be further studied.

For the isolated flow regime, we were not able to match the K1 [START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF] wind velocity profile with other cases of the same regime. In K1, we obtained W = 3.8 instead of 1.2 for other cases. The K1 case shares the same roughness configuration as the case of Leonardi et al. [START_REF] Leonardi | Direct numerical simulations of turbulent channel flow with transverse square bars on one wall[END_REF] but with a lower roughness Reynolds number, h + = 13 instead of 92. Hence, W seems to increase for low h + . Based on the fact that in average, the form drag is two orders of Downloaded by [Gang Huang] at 11:23 01 March 2016 magnitude higher than the friction drag, Krogstad et al. [33] qualified K1 as 'close to fully rough' . In our opinion, the relatively large value of W as well as the lower than average slip velocity reported in Section 3.3 are due to the viscous effect. This is supported by the locally large value of the friction drag in K1 (up to 15% of the form drag). [START_REF] Ashrafian | DNS of turbulent flow in a rod-roughened channel[END_REF] This nonnegligible viscous effect is the key characteristic of transitionally rough flow. Therefore, our result does not support the hypothesis of Bandyopadhyay [50] that h + = 10 marks the limit of fully rough flow, since 2D roughness elements with h + up to 13 probably still exhibits strong transitional roughness effects in terms of W.

The larger value of W for the transitionally rough case suggests that W relates to the aerodynamic smoothness of the surface, independently of the flow regime. Similarly, the roughness function U + defined in Equation ( 2) has been used as a measure of the aerodynamic roughness. [START_REF] Blocken | CFD simulation of the atmospheric boundary layer: wall function problems[END_REF] Yet, it is difficult to determine a critical value of U + for fully rough flow, not least because U + is generally proportional with h + based on experimental results. [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] From Equations ( 1), [START_REF] Perry | Rough wall turbulent boundary layers[END_REF], and ( 13), U + can be expressed as

U + = κ -1 ln h + + ln (ε/h) -(U + w + W ) + B, (16)
which confirms the proportionality between U + and h + for a given rough surface in the fully rough regime.

On the right-hand side of Equation ( 16), the first and the second terms have different signs and thus act in opposition to each other. In the first term, ε/h relates to the roughness exposure to the flow that contributes to the form drag, whereas the second term, proportional to U + w and W, is probably related to the effect of the viscous friction as claimed earlier. This reminds us the remark of Jiménez [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] on the double effects of the wall roughness: the creation of form drag and the weakening of viscous friction. The first effect is usually more pronounced than the second one since a roughened surface generally exhibits a higher drag. Yet, in particular conditions, for instance riblets, the second effect becomes disproportionately larger and the velocity is sufficiently enhanced in the viscous region to decrease the net drag. [START_REF] Choi | Direct numerical simulation of turbulent flow over riblets[END_REF] Similar phenomena appear to occur in the so-called 'd-type' roughness, [START_REF] Perry | Rough wall turbulent boundary layers[END_REF] for which the proportionality between U + and h + breaks down. [START_REF] Leonardi | Properties of d-and k-type roughness in a turbulent channel flow[END_REF] Further investigations on the balance between ε/h, U + w and W are indeed necessary in order to provide insights on the transitional roughness, the 'd-type' roughness as well as the drag reduction.

Application range and uncertainties of the study

Our study is limited to roughness elements distributed regularly and with simple geometry and equal height, missing variabilities in geometrical shapes and configurations inherent in realistic surfaces. [START_REF] Yuan | Estimation and prediction of the roughness function on realistic surfaces[END_REF] We have restricted intentionally the number of parameters describing roughness elements to two geometrical parameters: the roughness density λ and the spacing factor . Our study has demonstrated the crucial roles of both parameters in the parameterisation of the RS. The proposed law-of-the-wall, based on the geometrical and the roughness parameters, has the potential to be applied in large-eddy simulation (LES) models, in which the small-scale turbulent structures near the rough wall usually remain unsolved. [START_REF] Saito | Large eddy simulation of smooth-wall, transitional and fully rough-wall channel flow[END_REF] Several factors can influence the wall-law proposed in this study. The first one is related to the lack of information in some data sources. For example, values of h + is not reported Downloaded by [Gang Huang] at 11:23 01 March 2016 in [START_REF] Leonardi | Direct numerical simulations of turbulent channel flow with transverse square bars on one wall[END_REF] and an indirect estimation of h + based on the drag measurements is, therefore, done for this work. This may have impacted the determination of ε. Besides, errors in the original data due in particular to the spatial averaging process can affect the estimation of ε and z 0 . Additionally, the parameterisation of the slip velocity is sensitive to the friction velocity, whose estimation is still an open issue for rough-wall flow. [START_REF] Weber | Remarks on the definition and estimation of friction velocity[END_REF][START_REF] King | Investigations of the law-of-the-wall over sparse roughness elements[END_REF] For the roughness length, since z 0 is expressed using an exponential function (Equation ( 13)), its value is extremely sensible to the input variables such as the slip velocity. The uncertainty is thus relatively high. Despite these uncertainties, our work provides some insightful results on the RS. In our opinion, the aforementioned factors are unlikely to alter conclusions drawn in this study, based mainly on the relatively satisfactory results obtained in the formulation of the law-of-the-wall.

Furthermore, it will be interesting to see to what extent the universality of the log-law with a sublayer as introduced by von Kármán for smooth surface could still be a challenge to demonstrate, in case of rough surfaces with the concept of RS. For the present work, the VKC validity is not questioned as it is in the papers from [START_REF] Frenzen | On the magnitude and apparent range of variation of the von Karman constant in the atmospheric surface layer[END_REF] and [START_REF] Leonardi | Channel flow over large cube roughness: a direct numerical simulation study[END_REF] and more recently from [START_REF] Claffey | Evaluations of the von Kármán constant in the atmospheric surface layer[END_REF]. This last paper give a very good overview of the possible modellings of VKC leading to a value lower than 0.41 and questioning the universality of the VKC. With [START_REF] Yaglom | Data on turbulence characteristics in the atmospheric surface layer[END_REF] from [START_REF] Claffey | Evaluations of the von Kármán constant in the atmospheric surface layer[END_REF], we think that no data basis exists, for the moment, to reject the hypothesis that κ = 0.4 (without the second decimal place) for roughness surface. In the paper by Andreas et al. [START_REF] Claffey | Evaluations of the von Kármán constant in the atmospheric surface layer[END_REF], it is demonstrated how different authors try to reinterpret or to correct experimental in situ data with new value for VKC. They even tend to conclude for the universality of the VKC but for a lower value, of 0.39, than the originally accepted one. Nevertheless, the paper by Leonardi and Castro [45] reviews pertinent scaling that was introduced in the literature to assert the VKC with its classical value based on physical soundings. Even more, they attempt to introduce the variation of VKC with λ p (as defined in the present paper). As explained in their conclusions, Spalart admits that U + = f(z + and otherthings) is actually a 'frontal attack on the log-law' . They concluded with different scenarios as a trail for further research on the VKC universality issue. In some sense, we maintain Jackson's [START_REF] Jackson | On the displacement height in the logarithmic velocity profile[END_REF] objectives for 'VKC equal to 0.41 and fixed, over all λ p ' , but the displacement (ε here, z d for Jackson [START_REF] Jackson | On the displacement height in the logarithmic velocity profile[END_REF]) is not seen, here, as the direct height at which the axial drag acts. It takes into account the mixing layer above the roughness elements separately from the flow part inside the roughness elements (the inside part is not directly characterised by our proposal) that did not exchange directly with the upper RSL.

Conclusions

The RS over regularly distributed roughness elements of various configurations has been analysed. A universal non-dimensional velocity profile has been derived by matching velocity profiles from DNS and wind-tunnel experiments, for three flow regimes: isolated, wake, and skimming flow. This profile extends the classical logarithmic law-of-the-wall to the RS.

To achieve this result, we define a roughness scaling based on the effective height ε and a new roughness density parameter λ 2 accounting for the roughness element shape and inter-element spacing. It is shown that ε is the relevant length scale in the RS. The affinity of ε with the flow regimes renders its dependence on the surface properties straightforward. Compared to the usual frontal area density, λ 2 characterises 2D roughness separately from the 3D ones. Moreover, λ 2 benefits from a well-defined range (between 0 and 1), without an Downloaded by [Gang Huang] at 11:23 01 March 2016 explicit dependence on the roughness height h. We obtain a promising relationship between ε and λ 2 for a wide range of roughness configurations. Besides, the parameterisation of ε provides a systematic way to determine the displacement height z d in experimental and numerical studies.

From the generalised law-of-the-wall, the classical roughness length z 0 of the surface is reformulated using ε and λ 2 . It predicts successfully the peak of z 0 near the onset of the wake flow regime, and a higher z 0 for 2D roughness as well as staggered roughness, in agreement with previous experiments. This result emphasises the reliability of the new law-of-the-wall and the relevance of the geometrical parameters used in the parameterisation of the RS.

Finally, two important characteristics of the RS are analysed: the RS height (z r ) and the velocity departure from the log-law (W). Their dependence on the surface properties is studied for each flow regime. We find that z r decreases with the roughness density and reaches its minimum for the wake flow regime at about λ 2 ∼ 0.3 before the onset of the skimming flow regime. Our investigation reveals also that W, which is constant for each flow regime at fully rough state, could be viewed as a measure of the aerodynamic smoothness of the rough wall.

This work raised several issues that need further examination:

-The dependence of the slip velocity U + w on the Reynolds number needs to be clarified. Besides, since U + w is defined at the top of roughness elements, a proper investigation of the underlying canopy layer may also be necessary.

-The parameterisation of ε highlights the growing effect of other geometrical parameters such as the roughness shape with increasing roughness density. To better characterise the skimming flow regime, more sophisticated parameters based on the roughness topography may be appropriate, such as the root-mean-square roughness height, which is commonly used in the studies of irregular roughness. [START_REF] Mejia-Alvarez | Wall-parallel stereo particle-image velocimetry measurements in the roughness sublayer of turbulent flow overlying highly irregular roughness[END_REF][START_REF] Flack | The onset of roughness effects in the transitionally rough regime[END_REF] -Further, it will be interesting to see to what extent the effective height defined in our work, proposed on the concept of upper mixing layer exchanges with roughness cavities, could improve drag reduction profiles as a function of riblet geometry from the concept proposed by Bechert and Bartenwerfer,Luchini et al.,[START_REF] Luchini | Resistance of a grooved surface to parallel flow and cross-flow[END_REF][START_REF] García-Mayoral | Hydrodynamic stability and breakdown of the viscous regime over riblets[END_REF].

1 Introduction 1 . 1

 111 Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Atmospheric boundary-layer . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Wind erosion and sand dunes . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Wind-driven particle motion . . . . . . . . . . . . . . . . . . . . . . . 1.2 Formulation of the physical problem . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Wind and particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Wind and topography . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Topography and particles . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Our approach and goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Wind-tunnel studies . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Objectives and outline . . . . . . . . . . . . . . . . . . . . . . . . . . List of Figures 1.1 Mega-dunes of Badain Jaran desert. . . . . . . . . . . . . . . . . . . . . . . . 1.2 Diagram showing principal modes of particle motion. Figure taken from Greeley & Iversen [1987]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Interlink between wind flow, topography and particle transport.Figure adapted from Lancaster [2011]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The flow regimes associated with wind flow over building arrays of increasing spacing-to-height ratio. Figure taken from Oke [1988]. . . . . . . . . . . . . . 1.5 Sketch illustrating different zones around a transverse sand dune. Figure taken from Lancaster [2011]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Conceptional process of dune migration due to wind erosion and particle deposition events at preferential locations. Figure taken from Ernstsen et al. [2007]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Experimental set-up of PC09. Wind direction from right to left. . . . . . . . . 2.1 Sketch of the computational domain defined in a Cartesian coordinate system .x; y; z/. L x , L y and L y are the streamwise, the spanwise and the wall-normal extents of the domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 1

 1 Main forces acting on an airborne particle. . . . . . . . . . . . . . . . . . . . . 3.2 Different types of coherent structures. Figure taken from Guingo [2008]. . . . 3.3 Quadrant analysis defined by the relative values of fluctuating velocity components u 0 and w 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Representation of the collision process. Figure modified from Crassous et al. [2007] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The trajectory-crossing effect of a heavy particle (full circle) resulting from the inertial and the gravitational settling effects. The fluid parcel initially carrying the solid particle is denoted by an open circle. Figure taken from Shao [2009].

1 6

 16 Figure taken from Taniere et al. [1997]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Contours of the instantaneous streamwise velocity ("UU", filled contours) and the spanwise vorticity ("VOR", contour lines) in the mid-span plane obtained from our simulations. Velocity levels are scaled by the free-stream velocity U 1 , and the vorticity levels scaled by U 1 =ı. Coordinates are scaled using the hill height H. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Mean velocity profiles over smooth and rough hills. The velocity is scaled by the free-stream velocity U 1 . Points: experiments of Cao & Tamura [2006]. Full lines: LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 RMS profiles of the streamwise velocity fluctuation over smooth and rough hills. RMS velocities are scaled by the free-stream velocity U 1 . Points: experiments of Cao & Tamura [2006]. Full lines: LES. . . . . . . . . . . . . . . . . . . . . 4.4 Sketch of the wind tunnel set-up in the PC09 experiments. . . . . . . . . . . . 4.5 Numerical configuration (not to scale) of the single-hill configuration. . . . . . 4.6 Values of the non-dimensional spacings x C (blue line with squares) and z C (red line with circles) in the wall-adjacent cells for grid B along the streamwise direction around the Gaussian hill. . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Comparisons of mean streamlines over an isolated Gaussian hill between experimental data (left) and the LES data (right) for the V90 case. . . . . . . . . . . 4.8 Numerical configuration (not to scale) for the PC09 simulations. . . . . . . . .4.9 Mean streamlines around double Gaussian hills in rough-wall condition obtained by LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Comparisons between smooth-and rough-wall successive-hill cases with 3H spacing of mean longitudinal velocity and Reynolds stress profiles. Lines: smoothwall profiles. Open circles: rough-wall profiles. . . . . . . . . . . . . . . . . . 4.11 Mean longitudinal velocity and Reynolds stress profiles over Gaussian hills with 3H spacing in rough-wall condition. Lines: successive-hill profiles. Filled circles: upstream single-hill profiles. . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Mean longitudinal velocity and Reynolds stress profiles over Gaussian hills with 8H spacing in rough-wall condition. Lines: successive-hill profiles. Filled circles: upstream single-hill profiles. Open squares: successive-hill profiles between X=H D 4 and 1 shifted to the right by 8H. . . . . . . . . . . . . . . 4.13 The local hill height h scaled by the maximum height H, plotted as a function of the normalised streamwise position, X=H, relative to the hill top. Comparison between different 2D hill geometries. . . . . . . . . . . . . . . . . . . . . . . xvi List of Figures 4.14 Mean flow streamline plots showing the gradual appearance of a stable recirculation zone behind a 2D sinusoidal hill with increasing slope and decreasing ratio L=H. Flow is from left to right, with U 1 D 11:2 m s 1 . . . . . . . . . . 4.15 Separation length (L s D .x reat x sep /=H/ as a function of the maximum hill slope of a sinusoidal hill, obtained from the LES results in Table 4.10. . . . . . 4.16 Illustration of the key characteristics of the RZ behind a 2D hill. . . . . . . . . 4.17 Dependence of the recirculation zone behind the NC hill in smooth-wall (left) and rough-wall condition (right) on the flow regime. . . . . . . . . . . . . . . . 4.18 Evolutions of (a) the RZ circumference L R =H and (b) the RZ form ratio L AB =L CD as a function of the Reynolds number Re H . Comparison between rough (square) and smooth (circle) wall conditions. Filled symbol: sinusoidal hill. Open symbol: Gaussian hill (NC). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.19 RZ comparisons between the NC (left) and the L=H D 2:5 case (right) at different flow regime in the rough-wall condition. . . . . . . . . . . . . . . . . . 4.20 Comparison of the RZ characteristics, L R =H in (a), L cd =L AB in (b), between the NC (red, open square) and the L=H D 2:5 case (blue, filled square), at different Reynolds numbers in rough wall condition. CT06 simulation results at Re H D 10 4 are added for comparison. . . . . . . . . . . . . . . . . . . . . . . 4.21 Mean particle lift relative to the sum of gravity and adhesive forces as a function of the non-dimensionalised longitudinal coordinate X=H. Blue lines: values given by the mean wall shear stress. Red lines with circles: values from the event-based estimation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Snapshot showing that solid particles carried by the TBL are disturbed by the Gaussian hill, forming parcels of different concentrations on its lee side after hill crest. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Sandbox-type particle-feeding device. Figure adapted from Taniere et al. [1997]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Arrangement of closely packed grains: (a): a horizontal view of the particle layer. (b): a lateral view of two such layers separated by H p . . . . . . . . . . 5.4 Sketch of the two numerical configurations. . . . . . . . . . . . . . . . . . . . 5.5 Non-dimensionalized saltation mass flux as a function of the Shields parameter. Comparison between simulation results, analytical predictions and experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Concentration profiles at the end of the sandbox (first profile from left) and in Z02. Each profile is delimited using two dashed lines marking zero and maximum, respectively. The coordinate origin is set at the top of the upstream hill. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Concentration profiles in Z03 of the 3H case. The coordinate origin is set at the mid-point between the double hills. . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Concentration profiles in Z03 of the 8H case. The coordinate origin is set at the mid-point between the double hills. . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Comparison between particle velocity (symbols), fluid velocity undisturbed by particles (grey line) and fluid velocity modified by two-way effects (black line) inside Z03 for the 3H case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Comparison between particle velocity (symbols), fluid velocity undisturbed by particles (grey line) and fluid velocity modified by two-way effects (black line) inside Z03 for the 8H case. (a), longitudinal particle velocity u p ; (b), vertical particle velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 (a), Mean intensity and (b), time percentage of the lift-dominant events in zone Z01 as a function of the non-dimensionalised longitudinal coordinate X=H. Blue solid line: without particle. Yellow circle: with particle. . . . . . . . . . . 5.12 Mean intensity (left) and time percentage of the lift-dominant events (b) in zones Z02 and Z03 as a function of the non-dimensionalised longitudinal coordinate X=H. Blue solid line: without particle. Yellow line with circle: with particle. For each figure, the origin of the coordinated has been set, in (a), at the hill top and in (b) and (c), at the midpoint between the hill tops. . . . . . . . . . . . . . 5.13 Shear velocity u along the hill surface scaled by its upstream value u 0 in the TBL. U 1 D 11:2 m s 1 . Red colour corresponds to region where the flow attaches to the surface, blue colour to the backflow. Coloured arrows indicate a priori directions of particle transport in the near wall region. . . . . . . . . . . 5.14 Instantaneous map of the net accumulation of solid particles on the wall during a period of T D 30 ı 0 =u at U 1 D 11:2 m s 1 . Areas in erosion are marked in red and deposition is marked in yellow. Hill height in grey, not to scale. The deposition density D is scaled relative to "closely packed" roughness density at the sandbox box . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Instantaneous snapshot at Z02 showing particles trapped inside the recirculation region, bordered by separation streamline marked in white. Coloured contour for absolute fluid velocity (UF), grey-scaled contours for absolute particle velocity (UP). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii List of Figures 5.16 Definitions of the incoming particle flux Q I , trapping flux Q T , deposition flux Q D , escape flux Q out . Border of the recirculation zone is schematically represented by a dashed line. The dotted line indicates the reference plane through which incoming particles are tracked and analysed individually. h is the local altitude of the surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.17 Particle trapping flux Q T (blue square), deposition flux Q D (red circle) and incoming flux Q I (black triangle), scaled by Q 0 D p D p p gD p , as a function of the Shields parameter. Dashed line: threshold flow regime of particle transport. 125 5.18 Left-hand ordinate scale: percentage of trapped grains Q T (blue square) and deposited grains Q D (red circle) among total incoming particles Q I . Right-hand ordinate scale: ratio between deposited and trapped particles, Q D =Q T (black triangle). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.19 Instantaneous snapshot at Z03 showing particles trapped inside the recirculation regions between and behind Gaussian hills. Coloured contour for absolute fluid velocity (UF), grey-scaled contours for absolute particle velocity (UP). Both velocities are nondimensionalised by p gD p . . . . . . . . . . . . . . . . . . . 126 5.20 Particle trapping flux Q T (blue square), deposition flux Q D (red circle) and incoming flux Q I (black triangle), scaled by Q 0 D p D p p gD p , as a function of the Shields parameter. Comparison between 3H (filled symbols) and 8H cases (empty symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.21 Left-hand ordinate scale: percentage of trapped grains Q T (blue square) and deposited grains Q D (red circle) among total incoming particles Q I . Right-hand ordinate scale: ratio between deposited and trapped particles, Q D =Q T (black triangle). Comparison between 3H (filled symbols) and 8H cases (empty symbols).127 B.1 Comparisons between smooth-and rough-wall successive-hill cases with 8H spacing of mean longitudinal velocity and Reynolds stress profiles. Lines: smoothwall profiles. Open circles: rough-wall profiles. . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Mega-dunes of Badain Jaran desert.

Figure 1 . 2 :

 12 Figure 1.2: Diagram showing principal modes of particle motion. Figure taken from Greeley & Iversen [1987].

Figure 1 . 3 :

 13 Figure 1.3: Interlink between wind flow, topography and particle transport.Figure adapted from Lancaster [2011].

Figure 1 . 4 :

 14 Figure 1.4: The flow regimes associated with wind flow over building arrays of increasing spacing-to-height ratio. Figure taken from Oke [1988].

Figure 1 . 5 :

 15 Figure 1.5: Sketch illustrating different zones around a transverse sand dune. Figure taken from Lancaster [2011].

Figure 1 . 4 ,

 14 Figure 1.4, to the inter-dune flow dynamics, could potentially contribute to our understanding on the subject of solid particle transport around sand dunes.

Figure 1 . 6 :

 16 Figure 1.6: Conceptional process of dune migration due to wind erosion and particle deposition events at preferential locations. Figure taken from Ernstsen et al. [2007].

Figure 1 . 7 :

 17 Figure 1.7: Experimental set-up of PC09. Wind direction from right to left.

Figure 2 . 2 :

 22 Figure 2.2: A computational grid cell depicting the staggering of the coordinates and dependent variables. Figure taken from the ARPS User's Guide[START_REF] Xue | ARPS version 4.0 Users Guide[END_REF].

  3.

Figure 2 . 3 :

 23 Figure 2.3: Illustration of a computational domain generated using coordinate transformation (Equation2.30) and vertical grid stretching (Equation2.31). Note that only every second grid line is shown for clarity.

Figure 2 . 4 :

 24 Figure 2.4: Illustration of the spatial domain decomposition in ARPS from Xue et al. [1995].

Figure 2 . 5 :

 25 Figure 2.5: Rigid wall condition on a staggered grid. Two layers of the near-wall cells and one layer of boundary cells (under the physical wall) are pictured.

Figure 2 . 6 :

 26 Figure 2.6: Errors in the determination of the near-wall velocity gradient. Blue line: true velocity gradient. Red dashed line: approximation by Equation 2.41.

Figure 2

 2 Figure 2.7.

Figure 2 . 7 :

 27 Figure 2.7: Comparison between the base state profile (red line, Equation2.51 in Section 2.5.1)), the law-of-the-wall of von[START_REF] Von Kármán | The analogy between fluid friction and heat transfer[END_REF] (black, dotted line) and ofSpalding [1961] (blue, dashed line), and the DNS results of[START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Rj= 1410[END_REF] (circles).

Figure 2 . 8 :

 28 Figure 2.8: Illustration of the extraction/rescaling technique, applied to the extraction plane ("ex") and the inlet plane ("in").

Figure 2 .

 2 9a) compared to the base state TBL profile. The flow-through time is defined by L x =U 1 , L x being the streamwise extent of the domain (Figure 2.1). At this stage, the inflow velocity is mainly generated using random number sequences (red line in Figure 2.9a). After ten flow-through times, realistic eddies seem to appear in the extraction plane (blue line in Figure 2.8). Velocity fluctuations are thus assessed and recycled back to the inlet (red line in Figure 2.9b) according to the previously presented method.

  10 flow-through times.

Figure 2 . 9 :

 29 Figure 2.9: Two examples of instantaneous streamwise velocity profiles at the extraction plane and at the inlet, respectively. The base state profile is plotted for comparison.

Figure 2 .

 2 Figure 2.10 shows the mean velocity profile. The simulation results are compared with the experimental results of Vinçont et al. [2000], as well as the law-of-the-wall of Spalding [1961].The most noticeable difference lies within the buffer region. Satisfactory results are obtained for the mean flow inside the logarithmic region. Figure2.11 shows the profiles of the root-mean-

Figure 2 .

 2 11 shows the profiles of the root-meansquare (RMS) of the longitudinal and vertical velocity fluctuations, u 0 C and w 0 C

Figure 2 . 10 :

 210 Figure 2.10: Mean streamwise velocity profile. Blue squares: Experimental results of Vinçont et al. [2000] at Re D 500. Red circles: LES. Solid line, Spalding's law.

Figure 2 . 11 :

 211 Figure 2.11: RMS profiles. Blue squares: Experimental results of Simoëns et al. [2007] at Re D 500; Open symbols: streamwise RMS u 0 C. Filled symbols: wall-normal RMS w 0 C . Solid and dotted lines: DNS results ofSpalart [1988] with Re D 660, for u 0 C and w 0 C , respectively. Red, LES. Blue, experiments.

  1.

Figure 3 . 1 :

 31 Figure 3.1: Main forces acting on an airborne particle.

Figure 3 . 2 :

 32 Figure 3.2: Different types of coherent structures. Figure taken from Guingo [2008].

Figure 3 . 3 :

 33 Figure 3.3: Quadrant analysis defined by the relative values of fluctuating velocity components u 0 and w 0 .

Figure 3 . 4 :

 34 Figure 3.4: Representation of the collision process. Figure modified from Crassous et al. [2007]

Figure 3 . 5 .

 35 Figure 3.5. In this section relevant non-dimensional parameters related to these factors are introduced. Using these parameters, different modes of aeolian transport can be classified, based on both particle characteristics and flow conditions.

Figure 3 . 5 :

 35 Figure 3.5: The trajectory-crossing effect of a heavy particle (full circle) resulting from the inertial and the gravitational settling effects. The fluid parcel initially carrying the solid particle is denoted by an open circle. Figure taken from Shao [2009].

  mines whether particles participates to saltation or remain in suspension : particles are assumed be in pure suspension if St Á 1, and in pure saltation if St L 1. Yet a vast range of Stokes numbers spread between the two asymptotic regimes. In Taniere et al. [1997], g is used to further single out two intermediate regimes: "modified suspension" (St Á > 1, g < 1) and "modified saltation" (St L < 1, g > 1) as depicted in Figure 3.6.

Figure 3 . 6 :

 36 Figure 3.6: Particle behaviour based on the take-off curve of Foucaut & Stanislas [1996] using the Stokes number and the crossing-trajectories parameter. Figure taken from Taniere et al. [1997].

Figure 4 . 1 :

 41 Figure 4.1: Contours of the instantaneous streamwise velocity ("UU", filled contours) and the spanwise vorticity ("VOR", contour lines) in the mid-span plane obtained from our simulations. Velocity levels are scaled by the free-stream velocity U 1 , and the vorticity levels scaled by U 1 =ı. Coordinates are scaled using the hill height H.

Figures 4 .

 4 Figures 4.2 and 4.3 and discussed in this section. Locations at which the profiles are chosen are marked by a dashed line. The profiles are plotted up to Z=H D 4. In this layer the hill-induced modifications are most pronounced.

Figure 4 . 2 :

 42 Figure 4.2: Mean velocity profiles over smooth and rough hills. The velocity is scaled by the free-stream velocity U 1 . Points: experiments of Cao & Tamura [2006]. Full lines: LES.

Figure 4 . 3 :

 43 Figure 4.3: RMS profiles of the streamwise velocity fluctuation over smooth and rough hills. RMS velocities are scaled by the free-stream velocity U 1 . Points: experiments of Cao & Tamura [2006]. Full lines: LES.

Figure 4 .

 4 [START_REF] Jiménez | Turbulent flows over rough walls[END_REF] and flow characteristics are listed in Table4.3. Both smooth-and rough-wall conditions were studied. For the rough case, small PVC particles (with an average diameter of 200 µm), closely packed, are glued on the wall. We note the relatively fine size of wall roughness used in this study. In particular, the ratio z 0 =H D is 16% of the value used in[START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] with cube roughness. We note that the objectives are different between the two cases: sand transport for PC09 and urban canopy studies for CT06.

Figure 4 . 4 :

 44 Figure 4.4: Sketch of the wind tunnel set-up in the PC09 experiments.

Figure 4 .

 4 Figure 4.5 shows the numerical configuration of the domain, shared by all three grids, in a streamwise-vertical plane. The geometry is bidimensional and a periodical boundary condition is imposed in the spanwise direction. The extraction plane used for the onset of turbulent inflow is at 20H after the inlet. The hill centre is located downstream at 50H from the inlet. A remaining 14H extends from the hill top to the domain outlet, covering at least two times the classical RZ

Figure 4 . 5 :

 45 Figure 4.5: Numerical configuration (not to scale) of the single-hill configuration.

Figure 4 . 6 ,

 46 x C and z C are plotted as a function of the streamwise coordinate X=H, around a Gaussian hill at U 1 D 11:2 m s1 . It shows that the hilly terrain leads to a near threefold increase of x C up to 90 on the windward size, compared to the upstream values.

Figure 4 . 6 :

 46 Figure 4.6: Values of the non-dimensional spacings x C (blue line with squares) and z C (red line with circles) in the wall-adjacent cells for grid B along the streamwise direction around the Gaussian hill.

Figure 4 . 7 :

 47 Figure 4.7: Comparisons of mean streamlines over an isolated Gaussian hill between experimental data (left) and the LES data (right) for the V90 case.

Figure 4 . 8 :

 48 Figure 4.8: Numerical configuration (not to scale) for the PC09 simulations.

Figure 4 . 9 :

 49 Figure 4.9: Mean streamlines around double Gaussian hills in rough-wall condition obtained by LES.

Figure B. 1 , 1 Figure 4 . 10 :

 11410 Figure B.1, Appendix B.

1 Figure 4 . 11 :

 1411 Figure 4.11: Mean longitudinal velocity and Reynolds stress profiles over Gaussian hills with 3H spacing in rough-wall condition. Lines: successive-hill profiles. Filled circles: upstream single-hill profiles.

Figures

  Figures 4.12aand 4.12b.

  Figures 4.12aand 4.12b.

Figure 4 . 13 :

 413 Figure 4.13: The local hill height h scaled by the maximum height H, plotted as a function of the normalised streamwise position, X=H, relative to the hill top. Comparison between different 2D hill geometries.

  32°, L=H D 2:5

Figure 4 . 14 :

 414 Figure 4.14: Mean flow streamline plots showing the gradual appearance of a stable recirculation zone behind a 2D sinusoidal hill with increasing slope and decreasing ratio L=H. Flow is from left to right, with U 1 D 11:2 m s 1 .

Figure 4 . 15 :

 415 Figure 4.15: Separation length (L s D .x reat x sep /=H/ as a function of the maximum hill slope of a sinusoidal hill, obtained from the LES results in Table4.10.

  Hill with steep slope (32°).

Figure 4 . 16 :

 416 Figure 4.16: Illustration of the key characteristics of the RZ behind a 2D hill.

Figure 4 . 17 :

 417 Figure 4.17: Dependence of the recirculation zone behind the NC hill in smooth-wall (left) and rough-wall condition (right) on the flow regime.

Figure 4 . 18 :

 418 Figure 4.18: Evolutions of (a) the RZ circumference L R =H and (b) the RZ form ratio L AB =L CD as a function of the Reynolds number Re H . Comparison between rough (square) and smooth (circle) wall conditions. Filled symbol: sinusoidal hill. Open symbol: Gaussian hill (NC).

Figure 4 .

 4 Figure 4.19: RZ comparisons between the NC (left) and the L=H D 2:5 case (right) at different flow regime in the rough-wall condition.

Figure 4 . 20 :

 420 Figure 4.20: Comparison of the RZ characteristics, L R =H in (a), L cd =L AB in (b), between the NC (red, open square) and the L=H D 2:5 case (blue, filled square), at different Reynolds numbers in rough wall condition. CT06 simulation results at Re H D 10 4 are added for comparison.

Figure 4 . 21 :

 421 Figure 4.21: Mean particle lift relative to the sum of gravity and adhesive forces as a function of the non-dimensionalised longitudinal coordinate X=H. Blue lines: values given by the mean wall shear stress. Red lines with circles: values from the event-based estimation.

Figure 5 . 1 :

 51 Figure 5.1: Snapshot showing that solid particles carried by the TBL are disturbed by the Gaussian hill, forming parcels of different concentrations on its lee side after hill crest.

Figure 5 . 2 :

 52 Figure 5.2: Sandbox-type particle-feeding device. Figure adapted from Taniere et al. [1997].

Figure 5 . 3 :

 53 Figure 5.3: Arrangement of closely packed grains: (a): a horizontal view of the particle layer. (b): a lateral view of two such layers separated by H p .

Figure 5 .Figure 5 . 4 :

 554 Figure 5.4: Sketch of the two numerical configurations.

Figure 5 . 5 ,

 55 the saltation mass flux, scaled by Q 0 D p D p p gD p , is plotted as a function of the Shields parameter, Sh D u 2 =. p gD p /. The total mass flux is evaluated at the outlet of the domain.

Figure 5 . 5 :

 55 Figure 5.5: Non-dimensionalized saltation mass flux as a function of the Shields parameter. Comparison between simulation results, analytical predictions and experiments.

Figure 5 . 6 :

 56 Figure 5.6: Concentration profiles at the end of the sandbox (first profile from left) and in Z02. Each profile is delimited using two dashed lines marking zero and maximum, respectively. The coordinate origin is set at the top of the upstream hill.

Figure 5 . 7 :

 57 Figure 5.7: Concentration profiles in Z03 of the 3H case. The coordinate origin is set at the mid-point between the double hills.

  and 5.10a, and the vertical velocity profiles (w p ) are shown in Figures 5.9b and 5.10b. The profiles of wind velocity scaled by U 1 in the presence of and without solid particles are added for comparison, respectively.

Figure 5 . 10 :

 510 Figure 5.10: Comparison between particle velocity (symbols), fluid velocity undisturbed by particles (grey line) and fluid velocity modified by two-way effects (black line) inside Z03 for the 8H case. (a), longitudinal particle velocity u p ; (b), vertical particle velocity.

Figure 5 .

 5 Figure 5.11: (a), Mean intensity and (b), time percentage of the lift-dominant events in zone Z01 as a function of the non-dimensionalised longitudinal coordinate X=H. Blue solid line: without particle. Yellow circle: with particle.

Figure 5 . 12 :

 512 Figure 5.12: Mean intensity (left) and time percentage of the lift-dominant events (b) in zones Z02 and Z03 as a function of the non-dimensionalised longitudinal coordinate X=H. Blue solid line: without particle. Yellow line with circle: with particle. For each figure, the origin of the coordinated has been set, in (a), at the hill top and in (b) and (c), at the midpoint between the hill tops.

ure 5 .

 5 13, the irregular evolution of the local wall shear stress is represented by the ratio of the local friction velocity u to its upstream value u ref . By definition, u is given by

Figure 5 . 14 :

 514 Figure 5.14: Instantaneous map of the net accumulation of solid particles on the wall during a period of T D 30 ı 0 =u at U 1 D 11:2 m s 1 . Areas in erosion are marked in red and deposition is marked in yellow. Hill height in grey, not to scale. The deposition density D is scaled relative to "closely packed" roughness density at the sandbox box .

  2. a particle is considered to be deposited inside the RZ, if the trapped particle loses enough momentum to stay immobilised on the wall.At an upstream location of the relevant RZ, all particles passing through a reference plane are tracked individually during a fixed period. At each instant, their states of motion are updated based on their instantaneous locations using the principles proposed above. As an example, an instantaneous plot of particle positions and the absolute fluid velocity contours are given in Figure 5.15. Particles are coloured with their absolute velocity. Both velocities are nondimensionalised by p gD p . Particles are judged inside, or outside the RZ, by comparing its position with the instantaneous recirculation streamline, plotted in a white line in Figure 5.15. This streamline is a hypothesised border that separates the RZ from the outer flow. The knowledge of the instantaneous state of each incoming particle allows us to obtain statistics on the following particle fluxes: the incoming particle flux at the reference plane, Q I , the particle trapping flux through the border of the RZ, Q T and lastly, the deposition flux at the wall inside the RZ, Q D . These fluxes have the same dimensions as the saltation flux: [kg m -1 s -1 ]. An illustration depicting the relevant particle fluxes is given in Figure 5.16.

Figure 5 . 15 :

 515 Figure 5.15: Instantaneous snapshot at Z02 showing particles trapped inside the recirculation region, bordered by separation streamline marked in white. Coloured contour for absolute fluid velocity (UF), grey-scaled contours for absolute particle velocity (UP).

Figure 5 . 16 :

 516 Figure 5.16: Definitions of the incoming particle flux Q I , trapping flux Q T , deposition flux Q D , escape flux Q out . Border of the recirculation zone is schematically represented by a dashed line. The dotted line indicates the reference plane through which incoming particles are tracked and analysed individually. h is the local altitude of the surface.

  , the particle fluxes Q I , Q T and Q I , scaled by Q 0 D p D p p D p g, are plotted as a function of Sh. Q 0 is chosen as a good parameter to scale the particle flux since it is only dependent on particle characteristics fixed beforehand. The flux of incoming particles, Q I , increases with the Shields number in the same way as in the TBL case of Section 5.2. As for particles trapped inside the RZ (Q T ), and the deposited ones (Q D ), their flux reaches their maximum in the intermediate range of the Shields number. In fact, the number of trapped particles increases most rapidly at Sh D 0:07 (u D 0:34 m s 1 ), whereas the particle deposition most often occurs at Sh D 0:03 (u D 0:22 m s 1

Figure 5 . 17 :

 517 Figure 5.17: Particle trapping flux Q T (blue square), deposition flux Q D (red circle) and incoming flux Q I (black triangle), scaled by Q 0 D p D p p gD p , as a function of the Shields parameter. Dashed line: threshold flow regime of particle transport.

Figure 5 . 18 :

 518 Figure 5.18: Left-hand ordinate scale: percentage of trapped grains Q T (blue square) and deposited grains Q D (red circle) among total incoming particles Q I . Right-hand ordinate scale: ratio between deposited and trapped particles, Q D =Q T (black triangle).

5. 4 . 3 Figure 5 . 19 :

 43519 Figure 5.19: Instantaneous snapshot at Z03 showing particles trapped inside the recirculation regions between and behind Gaussian hills. Coloured contour for absolute fluid velocity (UF), grey-scaled contours for absolute particle velocity (UP). Both velocities are nondimensionalised by p gD p .

U

Figure 5 . 20 :

 520 Figure 5.20: Particle trapping flux Q T (blue square), deposition flux Q D (red circle) and incoming flux Q I (black triangle), scaled by Q 0 D p D p p gD p , as a function of the Shields parameter. Comparison between 3H (filled symbols) and 8H cases (empty symbols).

Figure 5 .

 5 Figure 5.21: Left-hand ordinate scale: percentage of trapped grains Q T (blue square) and deposited grains Q D (red circle) among total incoming particles Q I . Right-hand ordinate scale: ratio between deposited and trapped particles, Q D =Q T (black triangle). Comparison between 3H (filled symbols) and 8H cases (empty symbols).

  Figure . Schematic representation of the developed wall flow over an array of roughness elements with identical height h. The logarithmic profile of streamwise velocity (dashed line) is compared to the proposed law-of-the-wall (solid line). CL stands for canopy layer, RS for roughness sublayer with the upper limit z r , and IS for inertial sublayer.

[ 19 ]

 19 Downloaded by [Gang Huang] at 11:23 01 March 2016

)

  Downloaded by [Gang Huang] at 11:23 01 March 2016

Figure  .

  Figure . Non-dimensional velocity profiles U -as a function of z/h. Cases are classified into different flow regimes based on λ p according to Grimmond and Oke [], and are listed in order of increasing λ p (values on the right side of the case names). For reference, κ - ln (z/h) is plotted in a dotted line. (a) Isolated flow. (b) Wake flow. (c) Skimming flow.

Figure  .

  Figure . Non-dimensional velocity profiles U -as a function of z -. Cases are classified into different flow regimes according to their values of ε/h. Values of λ  (Equation ()) and ε/h are indicated on the rightside of the case name. Cases are listed in order of increasing λ  in the legend. (a) Isolated flow. (b) Wake flow. (c) Skimming flow.

Figure  .

  Figure . The effective height ε normalised by the roughness height h as a function of the roughness density λ, for D (a) and D (b) roughness elements. The solid line is an average fit to the 'reasonable envelop' given by Grimmond and Oke []. The dashed lines delimit the three flow regimes according to Grimmond and Oke []. Label I corresponds to isolated flow, II to wake flow, and III to skimming flow. Each symbol corresponds to one type of roughness element with one or several roughness densities. (a) D roughness. (b) D roughness.

Figure 5 .Figure  .Figure  .

 5 Figure . Arrangements of uniform, staggered, and aligned configurations composed of roughness elements of breadth b and length l. Dashed line denotes the minimal repeating unit of each roughness configuration. L x and L y correspond to the streamwise and spanwise dimensions of the unit, respectively. (a) Uniform. (b) Staggered. (c) Aligned.

Figure  .

  Figure . The slip velocity U + w as a function of the roughness density λ  for D roughness (solid symbols) and D roughness (open symbols). The solid line corresponds to Equation () and the dotted line to Equation ().

14 )

 14 1.55 ln λ 2 + 0.21, (Downloaded by [Gang Huang] at 11:23 01 March 2016

  2.8 Illustration of the extraction/rescaling technique, applied to the extraction plane ("ex") and the inlet plane ("in"). . . . . . . . . . . . . . . . . . . . . . . . . .

2.9 Two examples of instantaneous streamwise velocity profiles at the extraction plane and at the inlet, respectively. The base state profile is plotted for comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Mean streamwise velocity profile. Blue squares: Experimental results of Vinçont et al. [2000] at Re D 500. Red circles: LES. Solid line, Spalding's law. . . . 2.11 RMS profiles. Blue squares: Experimental results of Simoëns et al. [2007] at Re D 500; Open symbols: streamwise RMS u 0 C . Filled symbols: wall-normal RMS w 0 C

  . . . . . 161

	symbol symbol	unit unit	description description
	C dm u t	-m s -1	bulk aerodynamic drag coefficients for surface momentum fluxes threshold friction velocity
	C L u i	-m s -1	lift coefficient fluid velocity using Einstein notation
	C s U p	-m s -1	Smagorinsky coefficient particle velocity
	D p U r	m m s -1	particle diameter relative velocity of the particle with respect to the fluid
	Nomenclature D ? p -e xz -m s -1 hui u m s -1	nondimensional diameter defined by Foucaut & Stanislas [1996] rebound restitution coefficient in the impact plane Reynolds-averaged velocity Mean flow velocity
	e z U C w	-m s -1	restitution coefficient of the wall-normal velocity slip velocity
	F D V i	N m s -1	drag force exerted on particle impact velocity
	F L v i	N m s -1	lift force exerted on particle norm of the impact velocity
	Acronyms g V p	m s -2 m 3	acceleration of gravity particle volume
	H V r	m m s -1	external length scale of wall flow rebound velocity
	symbol h v r	m m s -1	description peak-to-valley roughness height norm of the rebound velocity
	ABL h C W	-m s -1	atmospheric boundary-layer roughness Reynolds number velocity deviation from the log-law in the RS
	ARPS h max x	m m	Advanced Regional Prediction System maximum height of the roughness element streamwise coordinate
	BL N h X p	m m	boundary-layer average height of the roughness element particle position vector
	CAPS k sgs z	m 2 s -2 m	Center for Analysis and Prediction of Storms subgrid-scale turbulent kinetic energy distance from the wall modified by the displacement height
	CFL l z 0	m m	Courant-Friedrichs-Lewy number streamwise length of roughness roughness length
	DBM l m z d	m m	wall model proposed by Duprat [2010] mixing length displacement height
	DNS L x z r	m m	direct numerical simulation streamwise extent of the domain upper limit of roughness sublayer
	IBM L y	m	immersed boundary-method spanwise extent of the domain
	LDA L y	m	laser Doppler anemometry wall-normal extent of the domain
	LES Re Greek symbols -	large-eddy simulation Reynolds number
	LMFA Re b	-	Laboratoire de Mécaniques des Fluides et d'Acoustique bulk Reynolds number U b H=
	PIV Re p symbol	-unit	particle image velocimetry particle Reynolds number description
	RANS Re "	-m	Reynolds Averaged Navier-Stokes friction Reynolds number (Kármán number) effective height
	RMS Re x Ä	--	root mean square flat plate Reynolds number von Kármán constant
	RSL S	m 2 -	wall model based on the formula of Huang et al. [2016] characteristic area of particle shape-dependent plan area density
	RSL S frontal 2	--	roughness sublayer frontal area of roughness spacing dependent roughness density
	SPD Sh f	--	wall model based on the formula of Spalding [1961] Shields number frontal area density
	TBL S plan p	--	tubulent boundary-layer plan area of roughness plan area density
	WMLES T L	s m 2 s -1	wall-modelled large-eddy simulation Large eddy turn-over time inside the TBL kinematic viscosity
	T L ˘-s	Lagrangian correlation time scale (fluid) roughness spacing ratio
	T p L p	s -	Lagrangian correlation time scale (grain) density ratio between solid particle and fluid
	Roman symbols u m s -1 Â i °impact angle fluid velocity vector
	u  r	m s -1 °rebound angle friction velocity
	symbol U b  xy	unit m s -1 °horizontal deviaton angle of rebound description bulk velocity of the flow
	B ı	U C	--m	additive constant in the log-law roughness function 99% boundary-layer thickness
	b C D ı ij U 1 u ? ˚p	m -m s -1 --	spanwise breadth of roughness drag coefficient free-stream velocity Kronecker tensor particle volume fraction

t

nondimensional friction velocity defined by

[START_REF] Foucaut | Take-off threshold velocity of solid particles lying under a turbulent boundary layer[END_REF]

  While Stokes law can be generally applied in the domain of biomechanics and microfluidic studies, experimental results show that the Stokes law begins to fail at about Re p D 10[START_REF] Shao | Physics and Modelling of Wind Erosion[END_REF]. For airborne particles transported inside the ABL, Re p generally exceeds this value.

						7, is used
	in this study
	C D	D	24 Re p	f .Re p / ;
	f .Re p / D	8 : 1 C 0:15 Re 0:687 0:0183 Re p p <	otherwise : if Re p < 1000 ;	(3.8)
	The gravity is the driving factors of particle saltation. By considering the gravity and ne-
	glecting the buoyancy effect for p		, the resulting equation of motion of airborne solid
	particles used in this thesis writes
		m p	dU p dt	D	2	C
				Re p D	jU r jD p ;	(3.6)
	The classical Stokes law is valid for small particle Reynolds numbers (Re p	1), in which
	viscous effects dominate inertial effects [Werner, 1990]. For this case, the Stokes law applies
	and writes
					C D D	24 Re p	:	(3.7)

/ U .X p .t/; t/ :

(3.5)

For solid particles, the particle Reynolds number can be constructed as For this range of particle Reynolds numbers, local flow dynamics around the particle depart greatly from the Stokes regime. Especially, the strong inertial effects of the particle bring the flow to separation at the rear of the particle, rendering Equation 3.7 inapplicable. For these high Reynolds-number flow regimes, empirical formulae exist in the literature for the evaluation of C D as the ones listed in

[START_REF] Clift | Bubbles, Drops and Particles Academic[END_REF]

. The one proposed by

[START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF]

, which accounts for higher order dependences of C D on Re p compared to Equation

3

.D SjU r jU r C P ;

(3.9)

Table 4 .

 4 1: Characteristics of the reference turbulent boundary-layer in the case of[START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF].

	Case	ı 0 (cm)	U 1 (m s -1 )	Re	Re H	H (cm)	z 0 (mm)	z C 0
	Smooth			3200	12 000		-	-
		25	5.84			4
	Rough			4900	10 500		0.2	3:97

C 0 , and Re H D U H H= , where U H is the undisturbed TBL velocity at the crest position.

Table 4 . 2

 42

	x=H	y=H	z mean =H	z min =H

: Domain sizes and mesh resolutions for the CT06 simulations.

Table 4 .

 4

			500	1200		-	-
		2.3
	V30 ?		500	1300		8.2	0.06
	V60	7.92	1600	4300	1	-	-
	V60 ?		1600	4000		8.2	0.18
	V90		2000	6000		-	-
		11.2
	V90 ?		2000	5400		8.2	0.25

3: Flow characteristics for the PC09 experimental cases. ?, rough-wall cases.

Table 4 .

 4

	x=H	y=H	z mean =H	z min =H

4: Domain sizes and mesh resolutions for the single-hill simulations.

Table 4 . 5

 45

	30	46	30

: Mesh resolutions in wall units for the single-hill simulation at U 1 D 11:2 m s 1 . The friction velocity corresponds to that of the undisturbed TBL.

Table 4 .

 4

	Wall model	Grid	x sep =H	x reat =H	.x reat x sep /=H	N P	CPU-H
		A	-	-	-	96	1423
	SPD
		B	0:73	3:55	2:82	96	1536
		A	0:47	3:82	2:35	96	1448
	DBM	B	0:51	4:59	4:08	96	1547
		C	0:30	5:45	5:15	192	16 939
	Number of processors.
	CPU hours necessary for 30 large eddy turn-over time scales (ı=u ).

6 summarises the results in terms of the streamwise locations of the separation and the reattachment points using different wall models (SPD and DBM) and grids. The reference data from the experiments are based on

[START_REF] Simoëns | Influence of Gaussian hill on concentration of solid particles in suspension inside turbulent boudary layer[END_REF]

. The amount of CPU-hours (CPU-H) and the number of processors (N P ) are also presented for each test case. The overall size of the flow separation is given by .x reat x sep /=H. The uncertainty of these results is x=H D 0:1 for cases A and B and x=H D 0:05 for case C. In general, a finer grid systematically predicts a larger separation length independently of the wall model used. Using SPD, the more refined grid B enables the LES to correctly predict flow separation. Using DBM, the simulation with grid A correctly predicts flow separation despite the coarseness of the grid. A substantial gain (70%) in terms of the separation length is achieved by using the more refined grid B instead of A. We note that this improvement is not a direct consequence of the grid refinement since A and B share the same x=H. By passing subsequently to the finest C, another gain of 25% is obtained using DBM.

Table 4 . 6 :

 46 Effect of the grid and wall model on the extents of the RZ at U 1 D 11:2 m s 1 over the smooth wall.

Table 4 .

 4 8: Domain size and mesh resolution for the PC09 simulations.

	x=H	y=H	z mean =H	z min =H

Table 4

 4

.3. All hill geometries are tested using smooth-wall condition. The rough wall condition is tested on most of the hill types with intermediate hill slopes. Details of flow parameters and hill geometries are given in Table

4

.9.

Table 4 .

 4

		Hill		U 1 (m s -1 )	ı 0 (m)	Re
	Mnemonic	Angle	H (m)
				11.2	0.05
	AC	17°0.01	11.2	0.07
				11.2	0.10
				2.3	0.07
				7.9	0.07
	NC	31°0.01	11.2	0.05
				11.2	0.07
				11.2	0.10
				2.3	0.07
				5.7	0.07
	NC ?	31°0.01	7.9	0.07
				11.2	0.05
				11.2	0.07
				11.2	0.10
	L/H=2.5	32°0.01	11.2	0.07
				2.3	0.07
	L/H=2.5 ?	32°0.01	7.9	0.07
				11.2	0.07
	L/H= 5	18°0.01	11.2	0.07
				2.3	0.07
	L/H=7.5	12°0.01
				11.2	0.07
	L/H=7.5 ?	12°0.01	11.2	0.07
				2.3	0.07
	L/H=10	9°0.01	7.9	0.07
				11.2	0.07
	L/H=10 ?	9°0.01	11.2	0.07

9: Summary of geometric characteristics of 2D hills and flow parameters. Angle corresponds to the maximum hill slope. ?, rough-wall condition. , Gaussian hill.

Table 4 .

 4 .15. 10: Influences of the hill slope on the lee-side flow separation at U 1 D 11:2 m s 1 behind a smooth sinusoidal hill. Data from[START_REF] Cao | Experimental study on roughness effects on turbulent boundary layer flow over a two-dimensional steep hill[END_REF] andCao & Tamura [2007] are added for comparison.

	Hill	Angle	Re H	x sep =H	x reat =H	.x reat x sep /=H
	L/H=10	9°6000
	L/H=7.5	12°6000	5:75	6:50	0:75
	L/H=5	18°6000	2:48	5:40	1:92
	L/H=2.5	32°6000	0:61	4:20	3:59
	Cao & Tamura [2007]	12°12 000
	Cao & Tamura [2006]	32°12 000	NA	5:40	NA

Table 4 .

 4 11: Effects of the Reynolds number on the extents of the RZ for NC cases over a smooth wall.

	The influence of the Reynolds number on geometric parameters of the RZ is further discussed
	in Section 4.3.3.
	Hill	Re H	x sep =H	x reat =H	.x reat x sep /=H	%
		1200	0:40	6:29	5:89	0.0
	NC	4300	0:42	4:78	4:37	-25.8
		6000	0:50	4:62	4:12	-30.0

Table 4 .

 4 14: Summary of the RZ characteristics, defined in Section 4.3.3.1, at different flow velocities and hill slopes. ?, rough-wall cases. , Gaussian hill geometry. By default, sinusoidal hill geometry and smooth wall condition.

		9°7.37	0.31	0.29	3.69	0.08	7.53	327.83
		12°5.98	0.34	0.50	4.53	0.11	9.20	283.10
		31°?	3.07	0.52	0.82	4.76	0.17	9.79	106.95
		32°?	2.45	0.59	0.77	4.50	0.17	9.29	108.24
		9°
	1600	31°	2.63	0.48	0.75	3.51	0.21	7.33	72.96
		31°?	2.98	0.61	0.94	4.68	0.20	9.80	96.11
		32°?	2.60	0.63	0.92	4.90	0.19	10.16	100.81
		9°
		9°?	7.10	0.20		1.00
		12°6.44	0.04		0.90
		12°?	5.03	0.36	0.22	2.57	0.08	5.20	135.91
	2000	18°3.95	0.26	0.24	1.77	0.14	3.60	71.11
		31°	2.62	0.46	0.71	3.43	0.21	7.14	70.13
		31°?	2.96	0.63	0.96	4.81	0.20	10.01	98.71
		32°2.35	0.46	0.72	2.82	0.25	5.94	62.92
		32°?	2.49	0.65	0.94	4.98	0.19	10.37	108.19

Table 5 .

 5 1: Particle and flow characteristics in the experimental configuration.

	millis)	St L	St Á	g	Sh

Table 5 .

 5 2: Analytical expressions of the nondimensionalised saltation mass flux Q=Q 0 as a function of the nondimensionalised friction velocity V D u =u t .

	Source u Sørensen [1991] D p (µm) p (kg m -3 ) u .m s 1 / 242 2500 OE0:18; 0:82	0:18	OE 7; 136	7
	Creyssels [2009]	242	2500	OE0:24; 0:67	0:20	OE11; 90	9
	Simoëns [2015]	200	1000	OE0:23; 0:46	0:23	OE32; 129	32
	ARPS	200	1000	OE0:11; 0:82	0:15	OE 7; 411	15

t .m s 1 / Sh (10 3 ) Sh t

(10 3 )

H

  Figure 5.8: Concentration profiles in Z03 of the 8H case. The coordinate origin is set at the mid-point between the double hills.

	X=H	-4.00	0.00	4.00
	c max	7.99 10 -5	4.08 10 -5	4.35 10 -5
	Z=H .c max /	1.05	1.30	1.05

Table 5 .

 5 6: Peak volume fraction and vertical location of the peak around the 8H hill from the simulation results.

Table 5 .

 5 .7. 7: Statistics on particle trapping and deposition in zone I behind the upstream Gaussian hill over 3000 instants.

	U 1	Sh	Incoming	Trapped Number	%	Deposited Number	%
	4.0	2.1 10 -2	17043	874	5.12	432	2.53
	5.7	3.1 10 -2	67331	2044	3.03	824	1.22
	7.9	7.0 10 -2	377306	3640	0.96	705	0.18
	11.2	1.2 10 -1	922440	2004	0.21	82	0.00
	20.0	4.1 10 -1	3383292	626	0.00	1	0.00

Table 5 .

 5 8: Statistics on particle trapping and deposition in zone I between two Gaussian hills over 3000 instants.

	1	Sh	Spacing	Incoming	Trapped Number	%	Deposited Number	%
	5.7	3.1 10 -2	3H	67537	104	0.15	11	0.00
			8H	67855	4297	6.24	1047	1.54
	7.9	7.0 10 -2	3H	373269	2191	0.58	840	0.22
			8H	372479	7159	1.92	1066	0.28
	11.2	1.2 10 -1	3H	899537	1511	0.16	382	0.00
			8H	901759	3865	0.42	237	0.00

  external length scale of wall flow h peak-to-valley roughness height h + roughness Reynolds number h max maximum height of the roughness element h average height of the roughness element

	l streamwise length of roughness
	L x streamwise spacing
	L y spanwise spacing
	Re b bulk Reynolds number U b H/ν
	S total area
	S frontal frontal area of roughness
	S plan plan area of roughness
	U spatially and temporally averaged streamwise velocity
	u * friction velocity
	U + roughness function
	U b bulk velocity of the flow
	U + w	slip velocity
	W velocity deviation from the log-law in the RS
	x streamwise coordinate
	z distance from the wall modified by the displacement height
	z 0 roughness length
	z d displacement height
	z r upper limit of roughness sublayer
	z d displacement height related to the effective height
	Greek symbols
	ε effective height
	κ von Kármán constant
	λ shape-dependent plan area density
	λ 2 spacing dependent roughness density
	λ f frontal area density
	λ p plan area density
	ν kinematic viscosity CL canopy layer DNS direct numerical simulation roughness spacing ratio
	IS inertial sublayer
	LES large-eddy simulation
	RS roughness sublayer
	VKC von Kármán constant
	Roman symbols
	B additive constant in the log-law
	b spanwise breadth of roughness

CONTACT G. Huang gang.huang@ecl.ec-lyon.fr H

Table  .

  Description of the data-set used in the study. z d is the displacement height reported by the authors. CH: channel flow. Stg.: staggered roughness spacings. Long.: longitudinal D cases. Sq.: square-shape cases. Tri.: triangular-shaped cases. Tran.: transversal D cases. EXP: experimental cases. BL: boundary-layer flow.

	Case	λ p		λ 	h +	ε/h	z d /h	z d /h	U + w	Description
	[]									DNS CH
	OL	.	.	.	.	.	.	.	.	Stg. cylinders
	OL	.	.	.	.	.	.	.	.	Stg. cubes
	OL	.	.	.	.	.	.	.	.	Uniform cubes
	OL	.	.	.	.	.	.	.	.	Long. sq. bars
	[]									DNS CH
	O	.	.	.	.	.	.	.	.	Long. tri. bars
	O	.	.	.	.	.	.	.	NA	Tran. tri. bars
	O	.	.	.	.	.	.	.	NA
	O	.	.	.	.	.	.	.	NA	Tran. round bars
	[]									DNS CH
	L	.	.	.	.	.	.	.	.	Tran. sq. bars
	L	.	.	.	.	.	.	.	.
	L	.	.	.	.	.	.	.	.
	L	.	.	.	.	.	.	.	.
	L	.	.	.	.	.	.	.	.
	L	.	.	.	.	.	.	.	.
	[]									DNS CH
	LC	.	.	.		.	.	.	.	Stg. cubes
	LC	.	.	.		.	.	.	.
	LC	.	.	.		.	.	.	.
	LC	.	.	.		.	.	.	.
	LC	.	.	.		.	.	.	.
	LC	.	.	.		.	.	.	.
	[]									DNS BL
	A	.	.	.	.	.	.	.	.	Stg. cubes
	A	.	.	.	.	.	.	.	.
	A	.	.	.	.	.	.	.	.
	A	.	.	.	.	.	.	.	.
	[]									Tran. sq. bars
	K	.	.	.	.	.	.	NA	.	DNS BL
	K	.	.	.	.	.	.	NA	.	EXP BL
	[,]									EXP BL
	C	.	.	.		NA	NA	.	NA	Stg. cubes
	C	.	.	.		NA	NA	.	NA	Aligned cubes
	C	.	.	.		NA	NA	.	NA	Stg. cubes
	C	.	.	.		NA	NA	.	NA	Aligned cubes

)

  Downloaded by [Gang Huang] at 11:23 01 March 2016

				L6 0.027 1.00	K2 0.082 1.00
				LC1 0.034 1.00	K1 0.082 1.00
				A4 0.055 1.00	A3 0.083 1.00
	U -			L5 0.063 1.00	Eq. 8 with W=1.2
				L4 0.082 1.00	Eq. 6 with W=3.8
	10 0	10 1	10 2
	z	-
		(a)
				A2 0.096 0.60	L2 0.18 0.35
				LC2 0.096 0.75	O3 0.19 0.30
				LC3 0.11 0.60	LC6 0.21 0.30
	U -			A1 0.15 0.55	OL2 0.21 0.30
				LC4 0.14 0.53	OL3 0.25 0.23
				L3 0.14 0.60	O4 0.26 0.13
				LC5 0.17 0.50	O2 0.43 0.15
	10 0	10 1	10 2	OL1 0.17 0.45	Eq. 8 with W=1.2
	z	-
		(b)

©  Taylor & Francis Downloaded by [Gang Huang] at 11:23 01 March 2016

Downloaded by [Gang Huang] at 11:23 01 March 2016

Remerciements

Appendix A

Mean velocity profile over regularly distributed roughness elements

In this appendix, we present an article entitled "A wall function for the mean velocity over regularly distributed roughness elements". This study relates the conceptional flow regimes proposed by [START_REF] Oke | Street design and urban canopy layer climate[END_REF]. to the roughness sublayer developed over over regularly based roughness elements, based on various DNS and experimental studies. A roughness length scale, the effective height, is derived and proved essential in achieving a universal velocity profile. A parametric study allows to relate this length scale to the geometric properties of the rough wall, in particular the roughness density. The paper has been submitted to the "Journal of Turbulence" and the status is currently under review. -for 2D roughness cases:

Combining Equation [START_REF] Raupach | Drag and drag partition on rough surfaces[END_REF] with Equation ( 14) or Equation ( 15) leads to z 0 /h as a function of λ 2 , as shown in Figure 8(a) and 8(b). On the same figure, z 0 /h is compared with the values obtained from several experiments [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF][START_REF] Krogstad | An experimental and numerical study of channel flow with rough walls[END_REF] and from the models of Kastner-Klein and Rotach [START_REF] Kastner-Klein | Mean flow and turbulence characteristics in an urban roughness sublayer[END_REF] and of MacDonald et al. [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF]. The expected peak of z 0 /h is well predicted by our parameterisation. This supports the claim that this peak is correlated with the onset of the wake flow regime. [START_REF] Jia | Effects of surface roughness element spacing on boundary-layer velocity profile parameters[END_REF] The value of λ 2 at the z 0 /h peak (0.1) is slightly lower than the value reported in [START_REF] Macdonald | An improved method for the estimation of surface roughness of obstacle arrays[END_REF] (0.15 ∼ 0.20). For the wake flow regimes (0.09 < λ 2 < 0.4),

Downloaded by [Gang Huang] at 11:23 01 March 2016

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix B Additional figures