Isabelle Calmet

Philippe Drobinski

Richard Perkins

Jean-Marc Lacome

Chantal Bourdelas

Chantal Verger

Véronique Dehlinger

Cécile

Merci À Dominique

Yannick

Sébastien

Thierry Aurélien

Meissam, Bo, Ruiwei, Rémy, Shupeng Merci À Amélie

Youngseob

Jean Sylvain

Arièle Laurent Beaudet

Denis

Guillaume

Salma Paul-Antoine Lucie

Nicola Augustin

Laurent Sylvie

Thi Dubus

Thu

Christophe Bénédicte

Maya Stéphanie

SSL Shear Surface Layer TKE Turbulent Kinetic Energy VOC Volatile Organic Compound

Keywords: Atmospheric dispersion, atmospheric inhomogeneous flow, Code_Saturne, stable atmosphere, surface heterogeneity, turbulence. v Dispersion atmosphérique, écoulement atmosphérique inhomogène, Code_Saturne, atmosphère stable, hétérogénéité de surface, turbulence. SL Surface Layer

Étude expérimentale et numérique de la turbulence et de la dispersion atmosphériques en conditions stables et en champ proche sur un site complexe Remerciement

Introduction

The lowest layer of the atmosphere, known as the atmospheric boundary layer (ABL), is the part of the atmosphere in contact with the ground, the Earth's ecosystem and the life of human being. During the day time, the ABL is driven by surface heating that creates the convective boundary layer (CBL), whereas radiative cooling of the ground during the night leads to the stable boundary layer (SBL). This stable boundary layer is common during the night, then also called the nocturnal boundary layer (NBL), but can also occur during the winter time with little insolation and can persist for days in polar regions. Over the past half century, the progress in understanding the CBL has far out paced the SBL, since the much stronger forcing in the CBL makes measurement and modelling of turbulence much easier [START_REF] Fernando | Whither the stable boundary layer? A shift in the resarche agenda[END_REF]. Conversely, progress in understanding the SBL has been restrained, because it often evolves continuously with generally weak turbulence and several scaling regimes that are associated with different physical mechanisms [START_REF] Grachev | Sheba flux-profile relationships in the stable atmospheric boundary layer[END_REF].

The emissions from human activities have great impacts on the ecosystem and also on human's own living conditions, such as the global warming or the airpollution. Thus, it is interesting to understand how the pollutants disperse within the atmospheric boundary layer. [Mylne, 1992] suggested that under stable conditions, the initial dilution of the pollutant at the source is low due to the low wind speed, and the rate of turbulent dispersion is reduced because the turbulence is suppressed by the stable stratification. Figure 1 shows an example of smoke dispersion in a stable boundary layer where the plume is trapped close to the ground due to the low wind speed and little turbulence in the SBL. The pollutant dispersion in the SBL is much more complicated than in the CBL, as stable conditions can induce large fluctuations of pollutant concentrations with possible occurrence of very high values. Therefore, it is important to investigate the turbulence structures in the stable ABL and to extend our understanding of dispersion to include stable conditions. Since pollutant dispersion in a stable atmospheric boundary layer and in complex environment is still challenging to model and difficult to reproduce in a wind tunnel [Fernando andWeil, 2010, Grachev et al., 2007], this topic is becoming of major interest in the field of air-pollution from human activities, such as industrial risks and road transportation. tant effort dedicated to improving its performance and validation by experimental measurements. Meanwhile, another model developed at CEREA called Polair3D is applied for larger scale modelling. Also, the most challenging part of these works is to accurately simulating the flow and pollutant dispersion under stable conditions with complex geometry.

In order to study pollutant dispersion in SBL and to validate the performance of Code_Saturne in modelling such stable conditions, an experimental programme consisting in measuring pollutant dispersion in a stratified surface layer and in nearfield (less than 200 m) has been carried out at a site named SIRTA (Site Instrumental de Recherche par Télédétection Atmosphérique), on the campus of Ecole Polytechnique, about 20 km south of Paris. SIRTA is the atmospheric observatory of IPSL (Institut Pierre-Simon Laplace) regrouping several laboratories of the Parisian region (including CEREA). It is a complex site including building areas, open fields, a lake and several forest canopies. The research activities at SIRTA mainly involve the following themes:

-Cloud process: cloud and fog formation and dissipation, synergy of observation, numerical modelling, weather and fog forecast, etc. -Climate study: better understanding of local climate variability, anticipation of extreme temperature and energy demand, characterisation of evolution of different variables, climate models assessment and validation, etc. -Turbulence and dynamics: experimental and numerical studies of dynamic processes in the atmospheric boundary layer, characterisation of turbulence with varying weather conditions, measurement and micro-scale modelling of flow and pollutant dispersion, etc. -Reactive gases and aerosols: identification of sources, understanding of transport mechanisms, long-distance transport study, real-time study of pollution episodes, etc. -Renewable energies: wind resource estimation, dynamic modelling, measurement and modelling of photovoltaic panels behaviour under real operational conditions, etc.

Within the research activities of EDF R&D, CEREA and SIRTA, this thesis entitled "Experimental and numerical study of atmospheric turbulence and dispersion in stable condition and in near field at a complex site" has been carried out with two main objectives:

Introduction

-The first objective is to perform an experimental study of the flow and dispersion at the SIRTA. Through processing and analysis of wind and concentration measurements, one can characterise the wind, the turbulence and the pollutant dispersion at this complex site and investigate impacts of atmospheric stability and terrain heterogeneity on their properties. -The second objective is modelling and simulation of the SIRTA experiment using the CFD code, Code_Saturne. Taking into account real topography and meteorological condition of the site, the numerical study allows to validate and improve the performance of Code_Saturne in simulating flow and dispersion and comparing with field measurements. Chapter 1 introduces theories and characteristics of turbulent flow and pollutant dispersion in the atmospheric boundary layer. Chapter 2 describes in detail the SIRTA flow and dispersion experimental program. Then, wind and turbulence data processing and analysis are presented in Chapter 3 both for a specific tracer gas release period, and for two years of continuous measurements. Chapter 4 introduces different approaches in turbulent flow modelling and in concentration dispersion modelling, and summarises the equations applied in the atmospheric module of Code_Saturne for the numerical study. Next, results of numerical study for mean and turbulent flow are presented in Chapter 5 with comparison between simulation profiles and field measurements. Finally, concentration data analysis and simulation are presented in Chapter 6 for different tracer gas release periods.

Chapter 1

Atmospheric boundary layer and pollutant dispersion

Before studying the SIRTA wind and dispersion field measurement campaign, the background and theoretical understanding of turbulence and dispersion in atmosphere is indispensable. We are particularly interested in the phenomenon in the lowest layer of the atmosphere -the atmospheric boundary layer (ABL) which plays an important role in many fields closely related to human activities, for example, air pollution and dispersal of pollutants, agricultural meteorology, weather forecasting and climate. This chapter will first briefly introduce the definition of the atmospheric boundary layer and its structure. Then, equations, theories and characteristics of mean flows and turbulence in the ABL will be described in details. Finally, some notions about air pollution and pollutant dispersion will be presented.

Definitions and structure

The Earth is surrounded by a gaseous envelope -the atmosphere. It protects life on the Earth by absorbing ultraviolet solar radiation, warming the surface through heat retention and reducing the temperature differences between day and night. It can be divided into several layers, and the limits are set by the discontinuities of temperature variations with altitude (Fig. 1.1): the troposphere (from 0 to 8-15 km) where the temperature decreases with altitude, the stratosphere (about 8-15 km to 50 km) where the temperature increases with altitude till 0 • C (273 K) , the mesosphere (from 50 km to 80 km) where the temperature decreases again with altitude down to 180 K, the thermosphere (from 80 km to 350-800 km) where the temperature increases again with altitude, and the exosphere (from 350-800 km to 50 000 km). The atmospheric boundary layer that we are interested is in the lowest layer of the atmosphere, the troposphere. Figure 1.1 -Atmospheric layers and the typical vertical temperature profile through the atmosphere [START_REF] Wallace | Atmospheric science: an introductory servey[END_REF] In the atmospheric context, it has never been easy to define precisely what the boundary layer is [Garratt, 1992]. Nevertheless, [Stull, 1988] gave a useful definition which identifies that the boundary layer is the part of the troposphere that is directly influenced by the presence of the earth's surface and responds to surface forcings with a time-scale of about an hour or less. These forcings include frictional drag, evaporation and transpiration, heat transfer, pollutant emission, and terrain induced flow modification.

In addition, there is a diurnal variation of temperature in the lower troposphere which is not evident for greater altitudes. Such variation is a key feature of the ABL. This diurnal variation is not driven by direct forcing of solar radiation, since most of radiation is transmitted to the ground and absorbed by it. So, it is the ground that warms and cools in response to radiation forcing, which in turn forces the changes in ABL via transport processes, such as wind flows convection or turbulence. [Stull, 1988] has also summarised that air flow in the ABL can be divided into 1.1. Definitions and structure three broad categories:

-Mean wind : responsible for rapid horizontal transport (advection). The mean wind speed is smaller near the ground due to friction. Horizontal winds are usually on the order of 2 to 10 m/s in the boundary layer. Vertical mean winds are much smaller, usually on the order of centimetres per second. -Turbulence : responsible for the vertical transport and consists of eddies of different sizes. The high frequency of occurrence of turbulence near the ground is one of the characteristics that makes the boundary layer different from the rest of the atmosphere. -Waves : frequently observed in the nocturnal boundary layer and can be generated locally by mean-wind shear or by mean flow over obstacles. They transport little heat, humidity and other scalar such as pollutants, but they can transport momentum and energy.

Each category can exist separately in the ABL. Transport of quantities such as moisture, heat, momentum and pollutants is dominated in the horizontal by the mean wind, and in the vertical by the turbulence.

The boundary-layer thickness is quite variable in time and space, ranging from hundreds of meters to a few kilometres. The structure of ABL is presented in Fig. 1.2. During the diurnal evolution, the ABL can be distinguished in various forms. The three main components are: a mixed layer highly turbulent, a residual layer less turbulent and a stable boundary layer with sporadic turbulence. At the bottom of the ABL, we can define a surface layer (SL) which is in direct contact with the ground and is about 10 % of the thickness of the ABL. The evolution of this surface layer, whether it is part of a mixed layer or a stable boundary layer, is directly related to the insolation.

During the day, the boundary layer receives the heat transferred from the heated ground at the bottom and loses heat at the top by radiative cooling. It leads the hot air climbing and cool air descending. This is the phenomenon of thermal convection which leads to the formation of a mixed layer. This mixed layer grows above an unstable SL as the ground temperature remains higher than the air temperature above.

At the sunset (Fig. 1.2), the ground stops receiving solar radiative heating and the thermal convection ceases. There are two layers developed in place of the mixed layer. The upper layer is the residual layer, which does not directly touch the ground and has the same initial state as the former mixed layer. This is a neutral layer Figure 1.2 -Nocturnal evolution of atmospheric boundary layer structure [Stull, 1988] with a uniform turbulence intensity in all directions. The lower layer is the stable boundary layer, which is characterized by stable air with little turbulence. In this layer, the statically stable air tends to suppress turbulence and affects the bottom of the residual layer gradually overnight. Moreover, this stable boundary layer can also be formed during the day (usually in winter), provided that the ground temperature is cooler than the air above.

At the sunrise of the next day, the mixed layer regains the place of the stable boundary layer and the residual layer.

As for pollutant dispersion, the behaviour is also different in different layers. In a mixed layer, the pollutant concentration is usually homogeneous on the vertical because of the thermal convection. However, in a stable boundary layer, since there is no thermal convection and little turbulence, pollutants disperse little in the vertical direction but more quickly in the horizontal direction. In order to study the pollutant dispersion, it is essential to have a thorough understanding of flows and turbulence.

Mean flow and turbulence in the surface layer

The turbulence nature of the ABL is one of its most conspicuous and important features [Garratt, 1992]. However, turbulence in the lowest atmosphere layer is different from most turbulence studied in wind tunnels. In the ABL, turbulence associated with thermal convection coexists with turbulence generated by wind shear. In the study of SIRTA field measurements, we are particularly interested in the mean flow and turbulence in the surface layer.

Fundamental equations

Flows in the SL are described by the fundamental equations of fluid mechanics adapted to the scale of the lower part of the atmosphere [Stull, 1988, Garratt, 1992]. Four governing equations are applied: the mass-conservation or continuity equation, the momentum equation, the temperature equation, and an additional conservation equation for scalar quantities such as pollutant concentration.

Mass-conservation equation

For a general compressible fluid, the mass-conservation equation is written as follows:

∂ρ ∂t + ∂ρu i ∂x i = 0 (1.1)
where ρ is the fluid density, and u i is the fluid velocity for component i.

The atmosphere is often considered anelastic and the temporal variation term of air density is usually considered negligible (∂ρ/∂t = 0), which reduces the massconservation equation to:

∂ρu i ∂x i = 0 (1.2)

Momentum equation

The momentum equation given by Navier-Stokes equations for the flow in the atmospheric SL and with Coriolis force neglected is written as follows:

∂u i ∂t 1 + ∂ ∂x j (u i u j ) 2 = - 1 ρ ∂p ∂x i 3 + 1 ρ ∂τ ij ∂x j 4 -gδ i3 5 (1.3) 10
Chapter 1. Atmospheric boundary layer and pollutant dispersion with 1: storage of momentum (inertia) 2: advection of momentum 3: pressure-gradient forces 4: influence of viscous stress 5: term of gravity (g = 9.81, m.s -2 ) where p is the pressure, and τ ij is the stress tensor. The air in the atmosphere can be considered as a Newtonian fluid. Thus, the term 4 can be written as:

1 ρ ∂τ ij ∂x j = 1 ρ ∂ ∂x j µ ∂u i ∂x j + ∂u j ∂x i - 2 3 ∂u k ∂x k δ ij (1.4)
where µ is the dynamic viscosity (kg.m -1 .s -1 ). Considering that the dynamic viscosity is spatially homogeneous, one can deduce a kinematic viscosity ν = µ/ρ (m 2 .s -1 ).

With the incompressibility approximation (∂u i /∂x i = 0), the momentum equation (1.3) becomes:

∂u i ∂t + u j ∂u i ∂x j = - 1 ρ ∂p ∂x i + ∂ ∂x j ν ∂u i ∂x j + ∂u j ∂x i -gδ i3 (1.5)

Temperature equation

The ideal gas law can adequately describe the state of gases in the boundary layer:

p = ρR * T v (1.6)
where p is the pressure, ρ is the air density, R * is the gas constant for dry air (R * = 287J.kg -1 .K -1 ), and T v is the virtual absolute temperature expressed as:

T v = T (1 + 0.61q) (1.7)
where q is the water-vapour mixing ratio.

A potential temperature Θ is usually used in the ABL. It is defined as the temperature of an air parcel brought adiabatically from its initial position to the reference pressure 100 kPa:

Θ = T p 0 p R * /Cp
(1.8)

with C p the heat capacity at constant pressure, and p 0 = 100 kPa.

Mean flow and turbulence in the surface layer

The expression is the same for the virtual potential temperature:

Θ v = T v p 0 p R * /Cp
(1.9)

It should be noted that, when significant moisture is present, the virtual temperature/potential temperature (T v / Θ v ) is applied in calculation in order to include the contribution of moisture to the air density. However, over land the difference between Θ and Θ v (T and T v ) is small and often ignored [START_REF] Kaimal | Atmospheric Boundary Layer Flows -Their Structure and Measurement[END_REF]. In this work, we follow the same practice, recognizing that wherever buoyancy effects are involved, Θ v (T v ) would replace Θ (T ).

Then, the temperature equation can be written as an equation for potential temperature Θ:

∂Θ ∂t + u j ∂Θ ∂x j = 1 ρC p ∂ ∂x j λ ∂Θ ∂x j + S Θ (1.10)
where λ is the thermal conductivity, and S Θ is an extra thermal source such as radiative flux, sensible heat flux or latent heat flux.

Conservation equation of a scalar quantity

Let c i be the concentration of a pollutant i. The advection-diffusion transport equation for pollutant concentration can be expressed as follows:

∂c i ∂t + u j ∂c i ∂x j = ∂ ∂x j D mi ∂c i ∂x j + S ci (1.11)
where D mi is the mass molecular diffusivity of pollutant i, and S ci is the body source term for the remaining processes not already in the equation, such as chemical reactions.

Reynolds decomposition

In the atmospheric surface layer, flows are usually very turbulent. The Reynolds decomposition is applied to study the turbulence. This common approach consists in splitting variables such as temperature and wind into a mean part and a perturbation part. The mean part shows the effect of the mean temperature and the mean wind, while the perturbation part represents the effect of turbulence. For variables in my work, the Reynolds decomposition gives:

u i = u i + u i (1.12a)
Chapter 1. Atmospheric boundary layer and pollutant dispersion p = p + p (1.12b)

Θ = Θ + Θ (1.12c) c i = c i + c i (1.12d)
It should be noticed here that for the air density, the Boussinesq approximation has been applied. As described in [Stull, 1988], we consider that the air density barely varies in the ABL and every occurrence of ρ can be directly replaced by its mean part ρ.

The rules for calculating the Reynolds average are summarized as follows:

f = 0
(1.13a)

f + g = f + g (1.13b) f g = f g + f g (1.13c) cf = cf (1.13d) ∂f ∂s = ∂f ∂s (1.13e) f ds = f ds (1.13f)
where f and g are two variables and c is a constant.

Governing equations for turbulent flow

Next, we replace all the turbulent variables in the fundamental equations by their Reynolds decompositions and make the Reynolds average for the entire equations. Then we obtain the governing equations for turbulent flow in the ABL as follows:

Mass-conservation equation

After equation (1.2) and (1.12a) , the Reynolds averaged mass-conservation equation with the anelastic approximation becomes:

∂ρ u i ∂x i = 0 (1.14)
1.2. Mean flow and turbulence in the surface layer 13

Momentum equation

With decomposition of u i (1.12a) and p (1.12b), the Reynolds averaged momentum equation for an incompressible flow (1.5) becomes:

∂u i ∂t + u j ∂u i ∂x j = - 1 ρ ∂p ∂x i + ∂ ∂x j ν ∂u i ∂x j + ∂u j ∂x i -gδ i3 - ∂ ∂x j u i u j (1.15)
The term u i u j corresponds to the Reynolds stress tensor which represents the turbulent fluctuations in fluid momentum.

Temperature equation

With decomposition of u i (1.12a) and Θ (1.12c), the Reynolds averaged temperature equation (1.10) becomes:

∂Θ ∂t + u j ∂Θ ∂x j = 1 ρC p ∂ ∂x j λ ∂Θ ∂x j - ∂ ∂x j Θ u j + S Θ (1.16)
where Θ u j represent the turbulent heat flux.

Conservation equation of a scalar quantity

In the same way, with decomposition of c i (1.12d), the transport equation of pollutant concentration (1.11) becomes:

∂c i ∂t + u j ∂c i ∂x j = ∂ ∂x j D mi ∂c i ∂x j - ∂ ∂x j c i u j + S ci (1.17)
where c i u j represent the turbulent concentration flux of pollutant i in the direction j.

Equations (1.14) -(1.17) are called together as Reynolds Averaged Navier-Stokes equations (RANS).

Turbulence closure

We can see that the use of Reynolds average reveals new second order terms for turbulent flux (u i u j , Θ u j and c i u j ) which are unknown. If we derive equations to eliminate these second statistical moments, we involve third order moment terms and so on. Thus, to close the equations system, it is necessary to find another approach.

In 1877, Boussinesq has introduced a turbulent viscosity hypothesis which is mathematically analogous to the gradient-diffusion hypothesis [Pope, 2000]. According to the hypothesis, the deviatoric Reynolds stress (-ρu i u j + 2 3 ρkδ ij ) is proportional to the mean rate of strain:

-ρu i u j + 2 3 ρkδ ij = ρν t ∂u i ∂x j + ∂u j ∂x i (1.18)
where ν t is the turbulent viscosity and k is the turbulent kinetic energy (TKE):

k = 1 2 u 2 + v 2 + w 2 (1.19)
Thus, for u i u j ,

u i u j = -ν t ∂u i ∂x j + ∂u j ∂x i + 2 3 kδ ij (1.20)
Similarly for Θ u j and c i u j :

Θ u j = - λ t C p ρ ∂Θ ∂x j (1.21) c i u j = - ν t Sc t ∂c i ∂x j (1.22)
where λ t is the turbulent thermal conductivity, λ t /C p ρ = ν t /P r t with P r t the turbulent Prandtl number, and Sc t is the turbulent Schmidt number.

Atmospheric thermal stratification

The variations of temperature and specific humidity with height lead to the variations of density in the vertical, that is the thermal stratification of the atmosphere. [Stull, 1988, Garratt, 1992, Arya, 1999] have all discussed the concept of stability in the ABL and several parameters allowing to define the stability condition.

A quantitative measure of static stability or instability, often used in meteorology, is the static stability parameter [Arya, 1999]:

N 2 = g Θ v ∂Θ v ∂z (1.23)
where N is called Brunt-Vaisala frequency or buoyancy frequency.

Mean flow and turbulence in the surface layer

Imagining there is a small parcel air displacing over a small vertical distance ∆z from its initial equilibrium position, the buoyancy force on the parcel can be expressed as

F b = - g Θ v ∂Θ v ∂z ∆z (1.24)
A stratified boundary layer is called statically stable if this buoyancy force acts to bring any vertically displaced parcel back to its equilibrium position, unstable if the buoyancy force acts to push the parcel farther away from its equilibrium position, and neutral if there is no force. The criteria can be summarized as:

-the ABL is stable if

F b < 0 and ∂Θ v ∂z > 0 -the ABL is unstable if F b > 0 and ∂Θ v ∂z < 0 -the ABL is neutral if F b = 0 and ∂Θ v ∂z = 0
This traditional view of stability is static and has limitation for application, especially when it is used as a measure of turbulent mixing and diffusion, for example, in air pollution and boundary layer meteorology. Because static stability is not sufficient, [Stull, 1988, Garratt, 1992] have proposed a more general interpretation of stability based on the turbulent kinetic energy. The TKE equation for the incompressible flow and the horizontally homogeneous case can be written as:

∂k ∂t + w ∂k ∂z = g Θ v w Θ v 1 -u w ∂u ∂z -v w ∂v ∂z 2 - ∂ ∂z w k 3 + ∂ ∂z ν ∂k ∂z 4 - 1 ρ ∂ ∂z w p 5 -ε 6 (1.25)
with 1: thermal production/destruction by buoyancy 2: dynamical production by shear 3: turbulent energy diffusion 4: molecular energy diffusion (can be neglected) 5: energy distribution by pressure fluctuations 6: dissipation

Chapter 1. Atmospheric boundary layer and pollutant dispersion

The terms describing shear production and buoyancy production (∂Θ v /∂z < 0) or destruction (∂Θ v /∂z > 0) are rather important in determining the intensity of turbulence [Garratt, 1992]. Term 2 is usually positive. If the ABL is stable, term 1 is negative and can compensate the production by shear; if the ABL is unstable, term 1 is positive and contribute to the production of turbulence in addition to the production by shear. The ratio between terms 1 and 2 is defined as the flux Richardson number:

R f = g Θ v w Θ v u w ∂u ∂z + v w ∂v ∂z (1.26)
Applying the Boussinesq's hypothesis expressions (1.20) and (1.21) for the horizontally homogeneous case, we can also define the gradient Richardson number:

R i = g Θ v ∂Θ v ∂z ∂u ∂z 2 + ∂v ∂z 2 (1.27)
Richardson numbers, especially the gradient one R i , have been widely used as a thermal stability parameter and the criteria are:

-the ABL is stable if R f , R i > 0 -the ABL is unstable if R f , R i < 0 -the ABL is neutral if R f , R i = 0
There is also another stability parameter, deduced historically by Obukhov using the TKE equation and used extensively as a thermal stability parameter in the surface layer, which is called the Monin-Obukhov length. By multiplying the TKE equation (1.25) by -κz/u 3 * , terms 1, 2 and 6 become dimensionless. The dimensionless term 1 is usually assigned the symbol ζ and is defined as

ζ = z/L M O , where L M O is the Monin-obukhov length. Thus, ζ = z L M O = - κzgw Θ v u 3 * Θ v (1.28)
and

L M O = -u 3 * Θ v κgw Θ v (1.29)
where κ = 0.42 is the Von Karman constant and u * is the friction velocity:

u * = (u w 2 + v w 2 ) 1/4
(1.30)

Their criteria for the stability condition can be summarized as:

-the ABL is stable if ζ > 0, L M O > 0 -the ABL is unstable if ζ < 0, L M O < 0 -the ABL is neutral if ζ → 0, L M O → ∞
In my thesis, the Monin-Obukhov length is the parameter that I will mainly use to characterize the surface layer stability. Moreover, the non-dimensional height ζ is a primary parameter of the Monin-Obukhov similarity theory in the surface layer (see next section).

Monin-Obukhov similarity theory

For a variety of flow situations, particularly turbulent flows in the boundary layer, mixing layer and others, our knowledge of the governing physics is not complete enough to derive flux-profile relations based on first principles [Arya, 1999]. It is necessary to develop useful empirical relationships between variables of interest, for example, turbulent flux and mean gradients. These relationships are particularly needed for defining boundary conditions such as inlet profiles in the numerical simulations.

Based on dimensional analysis, [START_REF] Monin | Basic laws of turbulent mixing in the surface layer of the atmosphere[END_REF]] hypothesized that any dimensionless characteristic of the turbulence can depend only upon the surface friction velocity u * , the height z, the buoyancy variable g/Θ v and the heat flux w Θ v , which means upon ζ = z/L M O . Monin-Obukhov similarity theory gives flux-profile relationships for velocity and temperature in the surface layer by universal functions as follows:

κz u * ∂u ∂z = φ m (ζ) (1.31) κz Θ * ∂Θ ∂z = φ h (ζ) (1.32)
where tively dimensionless wind shear and dimensionless temperature gradient [Arya, 1999, Grachev et al., 2007]. Integral forms of equation (1.31) and (1.32) give profiles relationships as follows:

Θ * = -Q 0 /u * is
u u * = 1 κ ln z z 0 -ψ m (ζ) (1.33) Θ -Θ 0 Θ * = 1 κ ln z z 0 T -ψ h (ζ) (1.34)
where Θ 0 is the potential temperature at z = z 0 T , z 0 is the dynamic roughness length, z 0 T is the thermal roughness length, and ψ m (ζ) and ψ h (ζ) are uniquely related to φ m (ζ) and φ h (ζ) respectively as:

ψ m (ζ) = z/L M O z 0 /L M O [1 -φ m (ζ)] dζ ζ (1.35) ψ h (ζ) = z/L M O z 0 T /L M O [1 -φ h (ζ)] dζ ζ (1.36)
Many different approaches can be found in the literature for φ m (ζ) and φ h (ζ). In my thesis, we have retained the expression for neutral condition according to [Stull, 1988]. As for stratified situations, inspired by a technical report of German Engineers Association (report [VDI-3783-Part-8, 2002]), expressions suggested by [START_REF] Dyer | Flux gradent relationships in the constant flux layer[END_REF], Businger et al., 1971, Hicks, 1976] are applied. The deduction of velocity and temperature profiles, and other turbulent parameters are summarized in the following.

Neutral surface layer

For wind profile in neutral conditions, we have:

κz u * ∂u ∂z = 1 (1.37) κz Θ * ∂Θ ∂z = 0 (1.38)
After integrations from roughness length z 0 to a height z for velocity profile, we obtain: Under neutral condition, the potential temperature profile is constant:

u(z) = u * 0 κ ln z z 0 (1.
Θ(z) = Θ 0 (1.40) With the relationship dΘ dz = dT dz + g C p ,
we can deduce the profile of T as:

T (z) = T 0 - g C p z (1.41)
with g/C p = 9.81/1005 = 0.0098 K/m and T 0 the temperature at the ground surface (In practice, we use temperature measurement at the lowest level, which is 2 m above the ground level in the SIRTA experiment).

Also, for turbulent kinetic energy k, dissipation rate ε and turbulent viscosity ν t , we have

k(z) = u 2 * 0 C µ (1.42) ε(z) = u 3 * 0 κz (1.43) ν t (z) = C µ k 2 (z) ε(z) (1.44)
with C µ = 0.03 for atmospheric flows according to the work of [Duynkerke, 1988, Katul et al., 2004, Zaidi et al., 2013].

As for the upper limit of the boundary layer z i , we suppose that z i = ∞ for the neutral case.

Stratified surface layer

In this case, a large range of stability condition are taken into account for wind and temperature profiles reconstruction [START_REF] Dyer | Flux gradent relationships in the constant flux layer[END_REF], Businger et al., 1971, Hicks, 1976]:

Chapter 1. Atmospheric boundary layer and pollutant dispersion

κz u * ∂u ∂z =              (1 -15ζ) -1/4 for ζ < 0 1 + 5ζ for 0 ≤ ζ < 0.5 8 -4.25 1 ζ + 1 ζ 2 for 0.5 ≤ ζ < 10 0.7585ζ for ζ ≥ 10 (1.45) κz Θ * ∂Θ ∂z =      0.74 (1 -9ζ) -1/2 for ζ < 0 0 for ζ = 0 0.74 + 5ζ for ζ > 0 (1.46)
According to [VDI-3783-Part-8, 2002] and also [START_REF] Musson-Genon | Reconstruction of the surface-layer vertical structure from measurements of wind, temperature and humidity at two levels[END_REF], after integrations from z 0 to a height z for the velocity profile, and from z 0 T to z for the temperature profile, we obtain:

u(z) =                                                u * κ ln z z 0 -2 ln 1 + X 1 + X 0 -ln 1 + X 2 1 + X 2 0 +2 arctan(X) -2 arctan(X 0 )] for ζ < 0 u * κ ln z z 0 + 5(ζ -ζ 0 ) for 0 ≤ ζ < 0.5 u * κ 8 ln (2ζ) + 4.25 (ζ) -1 -0.5 (ζ) -2 -ln (2ζ 0 ) -5ζ 0 -4] for 0.5 ≤ ζ < 10 u * κ [0.7585ζ + 8 ln(20) -11.165 -ln (2ζ 0 ) -5ζ 0 ] for ζ ≥ 10 with X = (1 -15ζ) -1/4 , X 0 = (1 -15ζ 0 ) -1/4 , and ζ 0 = z 0 /L M O . T (z) =                        T 0 - g C p z + u 2 * (z) κ 2 gL M O T 0 0.74 ln z z 0 T + 2 ln 1 + Y 1 -Y 0 for ζ < 0 T 0 - g C p z for ζ = 0 T 0 - g C p z + u 2 * (z) κ 2 gL M O T 0 0.74 ln z z 0T + 5(ζ -ζ 0 T ) for ζ > 0 1.2.
Mean flow and turbulence in the surface layer

21 with Y = 0.74 (1 -9ζ) -1/2 , Y 0 = 0.74 (1 -9ζ 0 T ) -1/2 , and ζ 0 T = z 0 T /L M O .
Similar to the neutral case, T 0 is the temperature at the ground surface, and u * has been calculated from the velocity profile equations by using U ref and z ref which are velocity and height at a measurement point. These profiles are particularly applied to generate inlet meteorological conditions in the numerical simulations.

As for estimation of the upper limit of the boundary layer z i , [START_REF] Holtslag | Air Pollution Modeling and Its Application VIII, volume 15 of NATO -Challenges of Modern Society[END_REF] suggests the following simple approach which restricts the upper limit of the boundary layer to minimum heights for stable conditions and maximum heights for neutral conditions:

z i = max 250 m, min 800 m, C 2 u * L M O f c
with C 2 = 0.7 and Coriolis coefficient f c = 10 -4 s -1 for our latitudes.

According to [START_REF] Kerschgens | Comments on turbulence parameters for the calculation of dispersion in the atmospheric boundary layer[END_REF], standard deviations of wind speed components σ u , σ v , and σ w can be described in dependence on boundary layer parameters L M O , u * , w * and z i :

-for neutral and unstable stratification:

σ u = [(2.4u * ) 3 + (0.59w * ) 3 ] 1/3 exp -0.3 z z i σ v = [(2.0u * ) 3 + (0.59w * ) 3 ] 1/3 exp -0.3 z z i σ w =    (1.3u * exp -0.3 z z i 3 + 1.3 z z i 1/3 1 -0.8 z z i w * 3    1/3
where w * is a standard convective velocity introduced in case of an unstable stratification with strong ascending forces prevail,

w * = g Θ w Θ z i 1/3 = u * - z i κL M O 1/3
-for stable stratification:

σ u = 2.4u * exp -0.3 z z i σ v = 2.0u * exp -0.3 z z i σ w = 1.3u * exp -0.3 z z i

Chapter 1. Atmospheric boundary layer and pollutant dispersion

Then we can deduce turbulent kinetic energy k for corresponding stable condition through [START_REF] Kerschgens | Comments on turbulence parameters for the calculation of dispersion in the atmospheric boundary layer[END_REF] has also suggested expressions for dissipation rate ε:

k(z) = 1 2 u 2 + v 2 + w 2 = 1 2 σ 2 u + σ 2 v + σ 2 w [
-for neutral and unstable stratification:

ε(z) = max u 3 * κz 1 - z z i 2 + 2.5κ z z i + w 3 * z i 1.5 -1.3 z z i 1/3 , u 3 * κz
-for stable stratification:

ε(z) = u 3 * κz 1 + 4 z L M O
For turbulent viscosity ν t , we always have

ν t (z) = C µ k 2 (z) ε(z)
with C µ = 0.03 for atmospheric flows according to the work of [Duynkerke, 1988, Katul et al., 2004, Zaidi et al., 2013].

Turbulence statistics, autocorrelation and spectrum functions

As mentioned in the Sect.1.2.2 for the Reynolds decomposition, the statistical description of a turbulent flow can be separated into large scale variations and small scale fluctuations:

a = A + a
where a is any variable in the turbulent flow.

We can define the variance and the standard deviation as:

σ 2 a = 1 N N -1 i=0 a i -A 2 = 1 N N -1 i=0 a 2 i = a 2 (1.47)
and

σ a = a 2 1/2 (1.48)
where N is the number of data points.

Mean flow and turbulence in the surface layer

By analogy, we can define the covariance as

covar(a, b) = 1 N N -1 i=0 a i -A b i -B = a b (1.49)
We define also the probability density function for the fluctuation p(a ). This function describes the probability that a (t) is between a (t) et a (t) + ∆a (t):

P rob[a < a (t) ≤ a + ∆a ] p(a )∆a (1.50) and ∞ -∞ p(a )da = 1 (1.51)
And the first and second momentum for the variable a(t) can be written as:

a = ∞ -∞ a p(a )da = 0 a 2 = ∞ -∞ a 2 p(a )da = σ 2 a
For the first momentum, we recover the definition for the fluctuation part in the Reynolds decomposition. For the second momentum, we recover the definition of the variance.

Another concept that is widely used in the statistical description of turbulence is the autocorrelation function. An autocorrelation function of time describes the relationship between any variable at time t and the same variable at later time t + τ as [Arya, 1999] :

R aa (τ ) = a (t)a (t + τ ) (1.52)
where τ is called the time lag, and we suppose that the variable is stationary. The normalised form of the autocorrelation function, which is called the autocorrelation coefficient, is more commonly used:

R aa (τ ) = a (t)a (t + τ ) a 2 (1.53)
R(τ ) is an even function which depends only on the time lag τ :

R(τ ) = R(-τ ) R(0) = 1
Chapter 1. Atmospheric boundary layer and pollutant dispersion

R(τ ) ≤ R(0) R(τ ) → 0, when τ → ±∞
We can see that the correlation is expected to be perfect when there is no time lag (R(τ = 0) = 1), and to decrease with increasing time lag, and become insignificant for a long time lag. Then, an integral time scale widely used in data analysis can be deduced as follows:

T i = ∞ 0 R(τ )dτ (1.54)
which allows to characterise how long the variable remains significantly correlated with its previous values. This time scale, representing the area under the curve R(τ ) versus τ by definition, is often used as the characteristic time scale of large eddies in the turbulent flows. However, in practice, T i is usually difficult to define because it is very sensitive to the instationary effect of the time series. Thus, a more useful integral time scale has been define as :

T e = τe 0 R(τ )dτ (1.55)
where τ e is the time lag for which autocorrelation coefficient R(τ ) goes below 1/e = 0.37 for the first time. This approximation method has been applied in previous work such as [START_REF] Moore | Estimates of integral time scales from a 100-m meteorocogical tower at a plains site[END_REF], Mylne, 1992, Kaimal and Finnigan, 1994, Irvine et al., 1997].

Analogous to the time autocorrelation function, we can define the spatial autocorrelation function that describes correlate values of any variable at different points in space at a fixed time, and we suppose that the variable is homogeneous:

R aa (r) = a (x, t)a (x + r, t) a 2 (1.56)
We can see that R aa (r) depends both on the magnitude and direction of vector r. So R aa (r) can be expressed by its three scalar components R aa (r 1 ), R aa (r 2 ) and R aa (r 3 ), for example:

R r 1 aa (r 1 ) = a (x, t)a (x + e 1 r, t) a 2 (1.57)
Thus, each component has the same properties as the time autocorrelation function, and the spatial correlation functions allow to characterize the three-dimensional large eddy structure in the turbulent flows. Also, we can deduce the corresponding integral length scale for each component as:

L r 1 aa = ∞ 0 R r 1 aa (r 1 )dr 1 (1.58)
The integral length scale can be simply related to the integral time scale through Taylor's frozen turbulence hypothesis (see Sect 1.2.7).

Furthermore, any variable as function of time may be expressed as a combination of many sine waves of different frequencies. Spectrum functions are used to described the frequency decomposition of time series and to indicate the relative contribution of ranges of frequencies to the variance [Arya, 1999]. So, a spectral density function can be defined such that:

a 2 (ω, ω + ∆ω) = S aa (ω)∆ω (1.59)
where a (ω, ω + ∆ω) represents the filtered portion of a (t) in the frequency range between ω and ω + ∆ω, and

S aa (ω) = lim ∆ω→0 lim T →∞ 1 ∆ωT T 0 a 2 (t, ω, ω + ∆ω)dt
(1.60) Also, we have:

∞ 0 S aa (ω)dω = a 2 (1.61)
A normalized spectrum function is more commonly used, we can deduce:

S aa (ω) = S aa (ω) a 2 (1.62) and ∞ 0 S aa (ω)dω = 1 (1.63)
With ω = 2πf , the spectrum function can be expressed as:

S(f ) = 2πS(ω) S(f )df = S(ω)dω
Mathematically, the autocorrelation and spectrum functions are related by the following Fourier transform relations [Arya, 1999]:

S(ω) = 1 π ∞ -∞ R(τ )e -iωτ dτ = 2 π ∞ 0 R(τ )cos(ωτ )dτ (1.64a) R(τ ) = ∞ 0 S(ω)cos(ωτ )dω (1.64b)
or in terms of the linear frequency f :

S(f ) = 4 ∞ 0 R(τ )cos(2πf τ )dτ (1.65a) R(τ ) = ∞ 0 S(f )cos(2πf τ )df (1.65b)

Taylor's hypothesis

In the study of turbulent flow or pollutants dispersion, Taylor's hypothesis is widely applied because it offers a simple relationship between the spatial and temporal statistics for turbulent variables. It was originally proposed by Taylor [Taylor, 1938] as a convincing assumption:

If the velocity of air stream which carries the eddies is very much greater than the turbulent velocity, one may assume that the sequence of changes in u at the fixed point are simply due to the passage of an unchanging pattern of turbulent motion over the point.

It assumes that the turbulent field is frozen in time and transported horizontally past the observer. Thus, It is also named as the frozen turbulence hypothesis, which enables us to deduce a direct correspondence between the spatial changes in the mean flow direction and the temporal changes at a fixed point for any turbulence variable. Taking the wind speed as an example, we have [Arya, 1999]:

u (x, y, z, t) = u (x -ut, y, z, 0) (1.66) and ∂u ∂t = u ∂u ∂x (1.67)
With this hypothesis, we can establish a simple relationship between the integral length scale and time scale (T i (1.54) or T e (1.55)) as follows:

L = uT i (1.68)
or with approximatio:

L = uT e (1.69)
In practice, instead of measuring in a large region of space at a given time, we can simply measure at a fixed point in space during a long time period. From these time series, we can deduce the length scales that characterize the eddies structure through equation (1.53), (1.55) and (1.69).

However, we have to bear in mind that turbulence is not really frozen. Taylor's simplification is useful only for those cases where the turbulent eddies evolve with a time-scale longer than the time it takes the eddy to be advected past a sensor [START_REF] Powell | An Investigation of the Application of Taylor's Hypothesis to Atmospheric Boundary Layer Turbulence[END_REF].

Kolmogorov's hypothesis and spectral characteristics

According to [Arya, 1999, Pope, 2000], small-scale motions in all turbulent flows have some universal characteristics at sufficiently high Reynolds numbers. These characteristics are described by a similarity theory proposed by Kolmogorov [Kolmogorov, 1941]. He assumed that small-scale motions are statistically isotropic at sufficiently high Reynolds numbers, and proposed two similarity hypotheses for locally isotropic turbulence.

The first hypothesis states that, at sufficiently high Reynolds number, there is a range of small scales or high wave numbers for which turbulence structure is in statistical equilibrium and is uniquely determined by the energy dissipation rate ε and kinematic viscosity ν [Arya, 1999]. Here, the statistical equilibrium refers to the stationarity of small-scale turbulence, which can exist even if the large-scale turbulence is not stationary. This equilibrium is justified because the characteristic time-scale of small-scale motions is much smaller than that of the most energetic large-scale turbulent motions. This hypothesis is also called the equilibrium range hypothesis, from which we can deduce the following Kolmogorov's microscales:

Length: η = ν 3/4 ε -1/4
(1.70)

Velocity: v = ν 1/4 ε 1/4 (1.71)
Chapter 1. Atmospheric boundary layer and pollutant dispersion Kolmogorov's second similarity hypothesis is concerned with the inertial subrange of eddies or wave numbers in which the energy is transferred from larger eddies to smaller eddies at a fixed rate determined by the rate of energy dissipation ε. It states that at sufficiently high Reynolds numbers, there is an inertial subrange in which turbulent kinetic energy is neither producted nor dissipated and the turbulence structure depends only on ε [Arya, 1999]. An important consequence of this hypothesis is that in the inertial subrange, the energy spectrum is given by

E(κ) = αε 2/3 κ -5/3 (1.72)
where α is a universal Kolmogorov's constant, and κ represent the circular wave number (rad.m -1 ). This equation is also called the Kolmogorov's law or the -5/3 power law. Many observations of turbulence spectra in the atmosphere and in laboratory flows at high Reynolds number have confirmed this law.

Figure 1.3 -Schematic of energy spectrum in the atmospheric boundary layer, where L is the integral length scale of turbulence and η is the Kolmogorov microscale (inspired by [START_REF] Kaimal | Atmospheric Boundary Layer Flows -Their Structure and Measurement[END_REF] )

Taking into account these Kolmogorov's hypotheses, the energy spectrum in the ABL can be divided into several regions and with typical characteristics of each. The schematic is shown in Fig. 1.3 and three major spectral ranges can be identified:

The energy-containing range: It contains the bulk of the turbulent energy which is produced by buoyancy and shear. The characteristic length scale is the Eulerian integral length scale L which represents the size of the most energetic eddies in the turbulent flows. E(κ) reaches its maximum at a wave number corresponding roughly to the integral length scale (that is κ ∼ 1/L).

However,we cannot be more precise because E(κ) is not easily measured in the atmosphere [START_REF] Kaimal | Atmospheric Boundary Layer Flows -Their Structure and Measurement[END_REF]] . The inertial subrange: Energy is neither produced nor dissipated but transferred to successively smaller and smaller eddies. This energy-transfer processes is named as energy cascade. The spectrum form in this range respects well the -5/3 power law (1.72). The dissipation range: Turbulent kinetic energy is converted into internal energy. Its characteristic length scale is the Kolmogorov micro-scale which characterizes the size and dynamics of the smallest eddies in turbulence.

The universal equilibrium range contains the inertial subrange and the dissipation range.

In practice, one-dimensional spectra S 11 (κ 1 ), S 22 (κ 1 ) and S 33 (κ 1 ) (subscript 1 for component in the streamwise direction) of the three wind velocity components are usually deduced separately since the spectral forms tend to be different with different corresponding integral scales. And their forms are somewhat different from E(κ) but are predictable in both the inertial subrange and the energy containing range. At the end of the spectrum on low frequency side, S(κ 1 ) is seen to level off to a constant value as κ → 0, which is a consequence of the one-dimensional representation of the three dimensional turbulence spectrum [START_REF] Kaimal | Atmospheric Boundary Layer Flows -Their Structure and Measurement[END_REF]. In order to identify energy peaks in this type of spectrum, meteorologists use the form κ 1 S(κ 1 ). Thus, the spectrum falls off as κ 1 S(κ 1 ) ∝ κ +1 1 on the low frequency side and κ 1 S(κ 1 ) ∝ κ -2/3 1 on the high frequency side.

Moreover, [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF], Drobinski, 2005] have summarised different spectrum forms for the three velocity components in the surface layer (Fig. 1.4).

Very close to the ground, there is a layer dominated both by shear and blocking effect of the ground which is called the eddy surface layer (ESL). In this layer, eddies coming from upper layers are impinging onto the ground and stretched Figure 1.4 -Sketch of velocity spectra form in the surface layer [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF]: S 11 and S 33 represent the spectra of longitudinal and vertical velocity components; f l and f u are the frequency limits with change of slope in the spectra.

along the wind direction. The elongated eddies become larger than the distance to the surface while the vertical component is reduced by blocking. In this topdown process, turbulent flows lose their isotropy, and a f S(f ) ∝ f 0 subrange at intermediate-frequency range can be observed in the velocity spectra only for the horizontal components. The height of the ESL is estimated to be on the order of 10 m in the near neutral condition. The spectra form in the ESL can be written as follows:

     f S ii (f ) ∝ f -2/3 for f ≥ f u f S ii (f ) ∝ f 0 for f u ≥ f ≥ f l for i = (1, 2) f S ii (f ) ∝ f +1 for f ≤ f l (1.73) and f S 33 (f ) ∝ f -2/3 for f ≥ f u f S 33 (f ) ∝ f +1 for f ≤ f u (1.74)
where f l is the lower limit of the f S(f ) ∝ f 0 subrange, and f u is the upper limit of the f S(f ) ∝ f 0 subrange or the lower limit of the inertial subrange for the spectrum of vertical velocity components.

For the sublayer above the ESL, the shear affects the isotropy of turbulence without the blocking effect of the ground. The spectra form in this sublayer, which is named as the shear surface layer (SSL), can be written as follows:

     f S ii (f ) ∝ f -2/3 for f ≥ f u f S ii (f ) ∝ f 0 for f u ≥ f ≥ f l for i = (1, 2, 3) f S ii (f ) ∝ f +1 for f ≤ f l (1.75)
which means that the subrange with 0 slope can be found in both horizontal and vertical velocity spectra.

[ [START_REF] Carlotti | Length scales in wall-bounded high-Reynolds-number turbulence[END_REF] has run LES simulations in neutrally stratified condition for a wall-bounded high-Reynolds-number turbulence and found evidence of the two sublayers of the surface layer which is consistent with Fig. 1.4. The spectra in [START_REF] Carlotti | Length scales in wall-bounded high-Reynolds-number turbulence[END_REF] showed that: very close to the ground (3 m height for S 11 and S 22 and 6 m height for S 33 ), a S(f ) ∝ f -1 subrange is visible on the longitudinal and transverse velocity fluctuation spectra only (eddy surface layer); at 22 m, a S(f ) ∝ f -1 subrange is found in the spectra of the three velocity fluctuation components and at 560 m, the turbulence displays isotropic characteristics.

Therefore, in the ABL, not only the spectra forms can be different between horizontal and vertical velocity components, they can also vary with height, environment homogeneity (presence of obstacles, forest, etc.) or even with stability conditions. This point will be discussed in the next Chapter with data analysis.

Atmospheric pollutant dispersion

Air pollution and air quality are problems that become more and more important nowadays. Air quality tends to be worse where air pollutants are emitted and transported. This phenomenon is particularly obvious in highly industrial areas or traffic-congested urban areas. So, local and urban air pollution is becoming a major problem since they are closely related to human activities and affect our living conditions. SIRTA field experiment is designed to have a better understanding of local pollutant dispersion in a semi-urban area under different meteorological conditions. In this section, notions about air pollution and pollutant dispersion will be introduced.

Air pollution problems

A condition of "air pollution" may be defined as a situation in which substances that result from anthropogenic activities are present at concentrations sufficiently high above their normal ambient levels to produce a measurable effect on human, animals, vegetation, or materials [START_REF] Seinfeld | Atmospheric Chemistry and Physics -From Air pollution to Climate Change[END_REF]. This definition could include any substance, natural or man-made chemical elements or compounds, whether noxious or benign. However, we are primarily concerned with substances that may cause significant undesirable effects.

Air pollution problems can be distinguished in three types according to the scales [Arya, 1999]: Type I is local/urban scale problems including indoor air pollution, local pollution around industrial sites (e.g. factories, power plants, waste disposal sites, etc.), and urban air pollution resulting from a variety of urban sources; Type II is regional/continental scale problems such as transport of sulphur and nitrogen oxides from major industrial areas to other regions, regional transport of ozone or particulates from forest fires, etc; Type III is hemispheric/global scale problems, which could lead to possible consequence on climate change, for example the increasing concentration of CO 2 and other greenhouse gases, the formation of the so-called ozone hole over Antarctica, and emissions from air-craft or volcanic eruptions.

All the air pollution problems have three main components [Arya, 1999]: (1) emission sources that produce air pollutants; (2) the atmosphere in which transport, diffusion, chemical transformations, and removal processes occur; (3) receptors (humans, animals, plants and materials) near the ground that respond to amounts of air pollutants reaching them. To fully understand and solve an air pollution problem, we need to study all the three components. Thus, air pollution is an interdisciplinary problem that requires efforts from specialists in vast domains such as chemistry, physics, meteorology, fluid mechanics, medicine, physiology, legislation and so on.

Origins and impact

The origins of air pollution are extremely varied. They can be divided into three broad categories as follows [Arya, 1999] Main concerns about the impact of air pollution are the effects on human health and on environment. Some of the effects are direct, specific and measurable, such as damages to vegetation and materials, many other effects are indirect and more difficult to measure, such as health effect on human beings and animals.

There are two main kinds of air pollution problems that have great impacts on human health, one is short-term, another is long-term. The short-term air pollution problems are usually due to accidental releases or extreme air pollution episodes, during which concentrations of pollutants in particular urban-industrial areas can reach excessively hight levels for several hours to several days. Exposure to relatively high air pollution levels over even short period can result in death and sickness. There are many examples in the history, for example, the 5-day great smog in London in the winter of 1952 caused nearly 4000 people died prematurely and made 100 000 more ill; the accidental release of heavy poisonous gas from a fertilizer plant in 1986 in Bhopal, India killed more than 3000 people and seriously injured many more. Nowadays, such acute exposure condition are rare due to awareness of government and air pollution control programs, and long-term air pollution problems become more common with growing urban population and increasing urbanization and industrialization worldwide. Long-term exposure to low ambient concentrations may not cause acute illness or deaths, but these exposures have been implicated in many other adversed health effects, such as respiratory and cardiovascular diseases; irritation of eyes, nose and throat; aggravation of an existing disease such as asthma; and impairment of performance of work and athletic activities [Arya, 1999]. Some epidemiological studies have already shown the relation between health effect and pollutants, for example, exposure to fine particulates is associated with increased incidence of respiratory illness, and eye irritation is associated with high concentrations of photochemical oxidants. How to reduce the pollutants emission in urban-industrial areas is still an issue that has been worked on by governments and many environmental agencies.

Impact of air pollution on environment can be divided into local and global effects. For local environment, air pollution can cause the reduction of visibility which is due to the absorption and scattering of light by airborne particles, for example urban smog. The gaseous as well as particles of pollutants in the atmosphere can significantly alter the radiation balance of the near surface layer in highly polluted urban atmosphere. On the one hand, the incident flux can be largely reduced by the absorption and scattering of particles before reaching the ground. On the other hand, the absorption and emission of long-wave radiation by the atmosphere may increase. Moreover, pollutants particles can contribute to the fog formation and precipitation since aerosols serve as cloud condensation nuclei. The primary pollutants such as SO 2 and NO 2 can transform into acidic species and turn to be the origin of acidic deposition (wet or dry), which could cause damage on vegetation and corrosion on building materials.

For the global impact of air pollution, the major concerns are stratospheric ozone depletion and climate change. Stratospheric ozone depletion, which is due to the extensive use of halogen compounds such as chlorofluorocarbons (CFCs), could possibly decrease the stratospheric temperatures and increase the amount of ultraviolet radiation reaching the ground. Meanwhile, the elevation of greenhouse gases concentrations in the air pollutants such as carbon dioxide, methane, nitrous oxide, CFCs and ozone is now recognized as the main cause of the global warming. Scientists have already made predictions for all kind of effects caused by the global warming, for example, the melting of glaciers and Antarctica ice cap, rising sea levels, flooding of low level coastal areas, destruction of estuaries, changes in regional temperature and rainfall patterns, and consequences on agriculture, forests and biological systems [Arya, 1999].

Dispersion in the atmospheric boundary layer

In the air pollution problem, my thesis concerns particularly the local part of atmospheric transport and diffusion of pollutants. Except for the exhaust emission from high-flying aircraft, almost all air pollutants from nature and anthropogenic sources are initially emitted within the ABL. Therefore, their short-range transport and diffusion are essentially determined by mean wind distribution and turbulence within the ABL [Arya, 1999]. The mean wind distribution essentially determines the average speeds and directions at which pollutants will be transported at different heights, while the fluctuations of turbulence determine the horizontal and vertical diffusion of pollutants around the cloud center.

Pollutants dispersion in the ABL depends strongly on the thermal stability condition. Figure 1.5 shows the schematics of some typical observed plume shapes and the corresponding wind speed and temperature profiles. For the stable case (a) with temperature inversion and small wind speed, vertical turbulence and diffusion tend to be suppressed by the stability. A highly meandering plume in the horizontal (not shown) appears due to variable horizontal wind direction. For the neutral case (d), the plume seems to have more or less equal spread in both vertical and lateral directions and looks like a cone. As for the fully convective case (c), the plume is taken by the large convective eddies and drifts upward and downward, causing the looping behaviour in the vertical. Overall, in a stable atmosphere with temperature inversion, turbulence vertical movements are attenuated and dispersion tends to be reduced; while in an unstable atmosphere, dispersion is usually amplified by the thermal convection.

Furthermore, pollutant dispersion also depends a lot on the turbulence mechanism and the ground surface characteristics. Thus, it is necessary to have precise description of turbulent flows in pollutant dispersion studies.

Figure 1.5 -Schematic of instantaneous plume patterns in the vertical and the corresponding wind speed and temperature profiles [Arya, 1999] Chapter 2

SIRTA flow and dispersion experimental program

In order to study pollutant dispersion in a stably stratified surface layer at small scales, a pollutant dispersion field experiment program has been conducted at the site SIRTA, "Site Instrumental de Recherche par Télédétection Atmosphérique" (Instrumented Research Site with Atmospheric Remote Sensing), on the campus of the Ecole Polytechnique in southern suburb of Paris. At this site, CEREA (Atmospheric Environment Teaching and Research Centre) has implemented many devices for meteorological measurements for several years in partnership with other laboratories in the Parisian region, such as LMD (Laboratoire de Météorologie Dynamique) and INERIS (Institut National de l'Environnement industriel et des RISques). In this chapter, background elements and objectives of the SIRTA experimental program will be introduced, measurement devices including their locations and operating principles will be presented, and several observation periods with gas releases will be described.

Context and objectives

Before planning this SIRTA experimental program, a literature review was conducted on previous campaigns performed by other research groups, including several campaigns organized by the Met Office in the UK [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF], Mylne, 1992, Mylne et al., 1996], the urban dispersion experiment titled Mock Urban Setting Test (MUST) [Biltoft, 2001], and the CASES-99 Field Experiment [START_REF] Poulos | Cases-99, a comprehensive investigation of the stable nocturnal boundary layer[END_REF], both held in the U.S.. It appears that the subject of pollutant dispersion in an urban environment has been fairly well documented and measurements were made on spatial scales ranging from several hundred meters to a few kilometers. There are campaigns such as CASES-99 in which the surface boundary layer has been studied in detail during strong thermal stability situations but without measuring the dispersion of pollutants. Some campaigns, such as those described by [Mylne, 1992], were well documented in terms of pollutant concentration measurements, but contained very few dynamic measurements. The MUST experiment was mainly dedicated to the study of flow dynamics and pollutant dispersion in an idealized urban field represented by rows of containers.

Taking into account the literature review, the objectives for the SIRTA turbulence and dispersion study have been defined to characterize at a complex site the fine structure of turbulence and associated dispersion through collocated high temporal and spatial resolution measurements of wind and concentration. These measurements are performed through an extensive network of ultrasonic anemometers measuring wind and turbulence and through photo-ionization detectors measuring concentration of a tracer gas. The instrumentation and the gas release mechanism in the MUST experiment have been adopted as a reference for the design of the SIRTA dispersion campaign.

The main features of the campaign are:

-Experiment in the near field (50 to 200 m); -Focus on stable thermal stratification, but it may include some neutral or slightly convective situations; -High-frequency measurements (about 10 Hz) to cover the entire frequency spectrum of fluctuations; -Large number of sensors measuring turbulence and concentration of tracer gas to document spatial inhomogeneities.

Figure 2.1 shows the different measurement areas in SIRTA domain. Our campaign is carried out in Zone 1 (Z1 on Fig. 2.1) which is limited in the north by a forest and in the south by a road. A lot of instruments are operated in routine mode at this site [START_REF] Haeffelin | Sirta, a ground-based atmospheric observatory for cloud and aerosol research[END_REF]. There were already some interesting measurements performed at this site which can inspire us on data analysis [Fesquet, 2008] and on modelling [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF]. In the future, it will be possible to study the effects of an idealized building by adding a container on the downstream side of the gas release point at the site.

Several specific meteorological conditions are required to perform this tracer release. After a climatological study based on observations between April 2007 and September 2009 on the site, we finally chose the following criteria:

-As the measurement zone is oriented along the east-west axis, wind direction Concerning the experimental devices (Fig. 2.2), dynamic measurements are mainly provided by ultrasonic anemometers, measuring three components of wind velocity (denoted u, v, w in meteorological coordinate) and speed of sound (from which the "sonic" air temperature T is derived, which is close to the virtual temperature) at 10 Hz. Propylene (C 3 H 6 ) is chosen as the tracer gas because of its low toxicity, low boiling point, low cost and low ionization potential. Photo-ionization detector (PID) is chosen to measure the gas concentration at 50 Hz, because of its good sensitivity to propylene. As indicate in [Biltoft, 2001] for the MUST campaign, the use of propylene/PIDs as the tracer/sampler combination is optimal for downwind sampling distances of several tens to several hundreds of meters because gas concentrations within these distances are expected to range from 0.5 to 500 ppmv, well within the PID operating range. Turbulence measurements have been recorded continuously since April 2012, while concentration measurements have been performed only during short releases of gas for specific meteorological conditions. These short gas releases last usually about an hour. Observation periods with gas releases are named as tracer tests or IOPs (Intensive observation periods) in the following.

Instruments and locations

As indicated above, the measurement devices include a large number of sensors documenting the flow dynamics and concentration of tracer gas at the same time. The source of tracer gas is situated at the "Dissemination" point. Around the source, there are four ultrasonic anemometers at four corners of a square. At 50 m from the source, there are five anemometers and five PIDs collocated and arranged in a circular arc. All the sensors above are at height of 3 m. In addition, there is a 10-m tower on the south-western part of the field on which is placed an anemometer at height of 10 m. And there is another 30-m tower on the south-eastern part of the field on which are placed two anemometers at height of 10 m and 30 m respectively.

One can notice in Fig. 2.3, that there are another five PIDs arranged in a circular arc at 100 m from the source. These PIDs will be used at a later stage of the dispersion campaign and are not included in the work of my thesis.

Moreover, as SIRTA is composed of an ensemble of state-of-the-art active and passive remote sensing instruments [START_REF] Haeffelin | Sirta, a ground-based atmospheric observatory for cloud and aerosol research[END_REF], there are many instruments operating routinely whose data might be helpful for us to have a better understanding of the atmosphere characteristics during measurements. For example, there were most of the time during the period considered for this work:

-one Lidar (Leosphere) measuring vertical profiles of three wind velocity components between about 40 m and 200 m, -one Sodar (Remtech) measuring vertical profiles of three wind velocity components and vertical velocity variance between about 100 m and 600 m, -one Radar UHF (Degreane) measuring vertical profiles of three wind velocity components between about 100 m and 3000 m, -one Microwave radiometer (RPG-HATPRO) which is a water vapour and oxygen multi-channel microwave profiler measuring temperature and humidity profiles from 0 m to about 10 km, -one Surface station including sensors measuring temperature, humidity, pressure, wind and precipitation, -etc.

In order to distinguish clearly every sensor in the following parts, a picture showing the locations and the labels of all the sensors involved in the measurements of the SIRTA dispersion campaign has been drawn in figure 2.4. represent ultrasonic anemometers. Four of them are placed at the corners of a 25m wide square, centred at the release point. They are named according to their positions: NE (north-east), NW (north-west), SE (south-east), and SW (southwest). Together, we denote them as "sonic square". There are five anemometers arranged in a circular arc with radius of 50 m, centred at the source. They are named according to the angle with respect to the east-west axis as 20N (20 • from this axis, on the north side), 10N, 0, 10S, and 20S. Together, we denote them as "sonic arc at 50 m". All the nine anemometers are at a height of 3 m. In addition, there is another ultrasonic anemometer at the top of a 10-m tower in the south-western part of the field, which is called "10mSW". Similarly, the other two ultrasonic anemometers, which are respectively at heights of 10 m and 30 m on a 30-m tower in the southeastern part of the field, are named "10mSE" and "30mSE". Green dots represent the PIDs. There are five of them collocated with the five ultrasonic anemometers at "sonic arc at 50 m". They are numbered from north to south by 1, 2, 3, 4, and 5. Another PID is positioned upstream of the source to measure background concentrations and is named "background".

In summary, turbulence is measured at three levels (3-m, 10-m, and 30-m height) with a dense network at 3 m, and concentrations are measured only at 3 m (measurements at a 10-m height have been obtained for the most recent IOP).

Furthermore, it should be noted that some shelters are present on the main platform at about 40 m to 50 m east from the release point.

Operating principles of main devices

The data that we are interested in the SIRTA dispersion campaign mainly comes from two types of sensors: ultrasonic anemometer and photo-ionization detector (PID). Therefore, it is useful to know their operating principles for a better data processing.

Ultrasonic Anemometer

Ultrasonic anemometer is a device that measures wind velocity in three directions and speed of sound. This instrument is used in many campaigns for direct measurement of small scale wind fluctuations and speed of sound from which one can deduce the heat flux. The operating principle of an ultrasonic anemometer is measuring the transit time of ultrasonic waves between three pairs of transducers carefully placed (Fig. 2.5 ). By calculating the difference between back and forth trajectory time, it allows to determine the various components of wind [START_REF] Gilbert | Un anemomètre ultra-sonique pour la mesure des paramètres moyens et turbulents de l'atmosphère[END_REF].

For example, the transit time for the back and forth trajectory between one pair of transducers are respectively:

t 1 = l c + V and t 2 = l c -V
where l is the distance between transducers, c is the speed of sound and V is the wind velocity along the axis of transducers.

It can be deduced as:

V = l 2 1 t 1 - 1 t 2 c = l 2 1 t 1 + 1 t 2
With the relationship between temperature T s and the sound speed for an ideal gas c = γRT s we have

T s = c 2 /γR
Here, T s is the so-called "sonic" temperature, which is close to the virtual temperature T v (see equation 1.7).

Anemometer uses three pairs of transducers oriented at a non-orthogonal coordinate system. Applying a method of passing the anemometer non-orthogonal reference to the meteorological reference, we can determine without ambiguity the components u, v, w of the wind in the meteorological reference. The effect of the orthogonal wind component to the measurements axis seems to have little impact on the accuracy of the measurement [START_REF] Gilbert | Un anemomètre ultra-sonique pour la mesure des paramètres moyens et turbulents de l'atmosphère[END_REF].

Then, from the instantaneous measurements at 10 Hz, we can calculate the mean values, variances and covariances, autocorrelation, etc. according to the needs of analysis.

Photo-ionisation detector (PID)

A photo-ionization detector is an ion sensor using photon energy to ionize the gas molecules. The gas is bombarded by photons, and electrons are extracted from the gas molecules which lose their electrons and become positively charged ions. So the ionized gas establishes an electrical current which is the signal output of the detector. The greater the gas concentration is, the more ions are produced, and the greater the current is. The current is then amplified and is proportional to the concentration. Next, the concentration is converted in ppm (parts per million) by using the calibration curve.

In our measurements, the PIDs measure the concentration value every 0.02 s (50 Hz). Before each measurement, they are all calibrated with an identical gain.

Description of IOPs

The campaign takes place when there is a forecasted "alert" in case of appropriate meteorological conditions. Intensive Observation Periods (IOPs) with gas releases have been performed since January 2012, informations about IOPs are summarized in table 2.1. The first phase took place between January and April 2012, mainly designed to test various components of the system and to verify the appropriateness of choice (instruments, meteorological conditions, etc.) and the ability to achieve the scientific objectives targeted under the experimental conditions. Next phases were planned from October 2012 with instrumental device set-up taking into account the feedback from the first phase, including improvement in gas flow control and sensors calibration.

In the first phase between January and April 2012, five IOPs have been carried out. The first IOP was conducted on 1st February 2012, which was the first experiment with full procedure of gas release, during which we focused on the technical aspects such as installation of the gas system and flow control of propylene. The second IOP (14 March 2012) was achieved with a gas release of approximately 22 min. But we had problems with PIDs measurements and difficulty in maintaining the gas flow. The third IOP (21 March 2012) lasted about 1.5 hours with a gas release of approximately 1 hour. The mass flow controller still failed to control gas flux properly, so we could not derive quantitative analysis from the concentration data. IOP-4 and 5 were planned on 27 March 2012 and 2 April 2012 respectively. However, there were no tracer test due to unsatisfactory wind conditions.

Second phase took placed since February 2013. After some modifications in the gas pipe, we succeeded to perform tracer tests with a better mass flow control. We tried to have tracer test as long as possible in a required meteorological condition, and turbulence measurements generally worked well. However, concentration measurements still have problems such as saturation or extreme high values due to inadequate gain in calibration. Also, during some IOPs, there were PIDs that simply did not work and needed some maintenance by the manufacture.

For data processing and analysis, the IOP-7 on 5 June 2013 has been chosen to present the results because it shows the best quality in PIDs measurements. Also, IOP-7 measured a relatively stabler surface layer than other IOPs, and its wind direction is quite appropriate to study the impact of the forest to the north of Zone 1. The most recent IOP-11 on 11 March 2015 was performed with an additional PID measuring concentration at height of 10 m over the position of PID-3. This IOP will be presented in Chapter 6, in the part of comparison between numerical simulation results and concentration measurements.

IOP

No. Table 2.1 -IOPs information ( * IOPs chosen for analysis)

Date

Chapter 3

Wind and turbulence data analysis

As mentioned in Chapter 1, a study of turbulence characteristics in SIRTA is necessary to understand the impact of terrain heterogeneity on three-dimensional flow, in order to prepare the study of dispersion and the numerical simulation. Thus, data processing and analysis of the SIRTA experiment are started with the measurements of wind. The raw data collected from ultrasonic anemometers includes three components of wind velocity and air temperature, which are taken every 0.1 s (10 Hz). This chapter will first present the procedures of data processing and the characterization of three-dimensional flow features by statistical analyses, correlations and spectral studies based on the measurements of IOP-7 (on 5 June 2013, from 18:48 to 20:17, see table 2.1). Then, through two years of continuous data, we will study the behaviour of the turbulence properties as function of surface layer stability condition and incident wind sector.

The same data analysis made for an earlier IOP on 21 March 2012 (IOP-3, see table 2.1) has been previously published in the article [START_REF] Wei | A preliminary analysis of measurements from a near-field pollutants dispersion campaign in a stratified surface layer[END_REF].

Data processing

Rotation of coordinate

The three components of wind velocity are given in raw data in the conventional meteorological reference noted as (u, v, w) (figure 3.1(a)). However, in order to highlight the theoretical characteristics of the atmospheric turbulence that we wish to observe, it is better to change the two horizontal wind components u and v into the longitudinal and transverse direction of the wind axis, which are denoted by a and b. So, the first step in sonic data processing consists of a reference frame rotation from meteorological frame of reference (u, v, w) to the one in the mean wind frame (a, b, w) (figure 3.1(b)).

By keeping the vertical wind direction, the reference rotation is carried horizontally in the following way: where dd is the instantaneous horizontal wind direction between 0 • and 360 • which is calculated previously from u and v :

dd = arctan -u -v
We define also the instantaneous horizontal wind speed as:

f f = √ u 2 + v 2 3.1.2 "Quick Look"

and sub-period selection

The instantaneous wind velocity (u, v, w) and temperature measurements provided by IOP-7 include approximately 53900 values. Before statistical analysis, we simply plot all the data to get an idea about their values and variation trends ("Quick Look").

The "Quick Look" figures of instantaneous three wind velocity components (u, v, w), temperature T s and horizontal wind direction dd are shown in appendix A.1. We can see that the raw data values are very fluctuating especially for the wind direction. Therefore, a data filter is needed to smooth out short-term fluctuations and highlight longer-term trends. Still, we can observe from the raw data that the fluctuation is around -1.5 m/s for u component, and around 0 for v and w components, which corresponds well to the meteorological conditions of low velocity and easterly wind required during tracer tests.

As for the data filter, a central moving average over 1 min has been chosen. That is to say, the value of nth min in the averaged series is calculated by averaging the values between the n-0.5 min and the n+0.5 min in the original series. Thus, there is no values for the first half minute and the last half minute in the average series. We obtained the averaged horizontal wind direction dd in figure 3.2. Other curves of averaged wind velocity and temperature are shown in appendix A.2.

From the curves of averaged data (figure 3.2 and appendix A.2) , we can observe that the 9 anemometers at height of 3 m show generally the same trend and variation with slight differences between four at "sonic square" and other five at "sonic arc at 50 m", which implies a proper functioning of all the anemometers. Moreover, we find that there is a strong variation in wind direction for several anemometers before 20th min of the IOP-7. And after 80 min measurements, wind direction dd tends to turn north and temperature T s begins to decrease while wind velocity u starts to increase (Fig. 3.2, A.7 and A.10).

In order to avoid the effects of these variations in analysis afterwards, we decided to choose a sub-period during which meteorological conditions are almost stationary for data processing. The choice is especially made with consideration of horizontal wind direction dd. It is found that dd is more stable for the measurements period between 20th min and 80th min during the IOP-7. Also, the wind velocity (u, v, w) and temperature T s are quite stable within this sub-period of 60 min. Therefore, we choose this sub-period for data processing and analysis of the IOP-7 in the following. It should be noticed that in practice, we take 36000 values (from 12000 to 48000 points of measurements) from the data of the IOP-7. These 36000 values correspond to the measurements from 19:08 to 20:08 on 5 June 2013.

Data analysis for IOP-7

Our turbulence data analysis is inspired by many previous studies.

[ [START_REF] Barthlott | Long-term Study of Coherent Sturctures in the Atmospheric Surface Layer[END_REF] have studied the coherent structures in the surface layer by Figure 3.2 -'Quick Look' at the horizontal wind direction dd filtered by central moving average over 1min of 9 anemometers at "sonic square" and "sonic arc at 50 m" during the IOP-7; from top to bottom: NE, NW, SE, SW, 20N, 10N, 0, 10S, 20S.

using wavelet analysis from SIRTA measurements. [Fesquet et al., 2009a] have investigated the impact of terrain heterogeneities on the average turbulence variables, by using SIRTA measurements, and shown that the terrain complexity has an important impact on turbulence variables such as wind direction, turbulent kinetic energy and momentum fluxes. [START_REF] Horst | HATS: Field Observations to Obtain Spatially Filtered Turbulence Fields from Crosswind Arrays of Sonic Anemometers in the Atmospheric Surface Layer[END_REF] have explored measurements from the Horizontal Array Turbulence Study (HATS) campaign through spatial cross correlations and examined the validity of Taylor's hypothesis for neutral stratification. [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF] have analysed the layered structure of the near-neutral surface layer especially through longitudinal and vertical velocity spectra by using the CASES-99 measurements, and shown that the turbulence anisotropy is due to a top-down mechanism (blocking of impinging eddies). Some previous studies for turbulence over canopy are also helpful to understand the impact of the forest at the site. [START_REF] Irvine | The evolution of turbulence across a forest edge[END_REF], Thomas and Foken, 2007, Chahine et al., 2014] have analysed field measurements for flows over a tall spruce canopy, a forest edge and a vineyard respectively, and have shown similar results for vertical profiles of turbulence variables and wind velocity spectra for flows over a vegetation zone. It has to be noticed that most studies have been limited to near-neutral cases, except [START_REF] Thomas | Organised motion in a tall spruce canopy: temporal scales, structure spacing and terrain effects[END_REF] who has included the mean wind profiles for stable and unstable stratifications.

Thus, wind and turbulence data analysis for the SIRTA experiment are carried out here through statistical analysis of turbulent variables, study of integral length scales, velocity correlation and velocity power spectra.

Statistical analysis

After selecting a data sub-period with relatively stationary weather conditions, the next step is to calculate statistical values for the time series such as mean, variance and other variables to characterize atmospheric turbulence.

First of all, we calculate mean values for (u, v, w) and T s , and we note them as u mean , v mean , w mean et T mean . We can deduce mean horizontal wind direction dd mean as:

dd mean = arctan -u mean -v mean
Next, we calculate instantaneous longitudinal and transverse wind velocity: For all the anemometers, b mean has been verified to be zero (∼ 10 -16 ) which is exactly the definition of transverse component of mean wind velocity. For a given averaging period, the longitudinal mean wind velocity a mean has been verified to have no difference with the horizontal mean wind speed f f mean . In the following comparison between the simulation results and the measurements (in Chapter 5), we compare directly the simulated horizontal mean wind speed f f mean with the measured longitudinal mean wind velocity a mean .

a b = -
Turbulence statistic such as mean wind direction dd mean and longitudinal velocity a mean , variances of the three wind components (σ 2 a , σ 2 b , σ 2 w ), turbulent kinetic energy k (1.19) , friction velocity u * (1.30) , vertical heat flux Q 0 = w Θ and Monin-Obukhov length L M O (1.29) are reported in table 3.1. They have been calculated over the selected sub-period of 60 min in the IOP-7 during which meteorological conditions are almost stationary.

We can observe that the statistical values for anemometers at same level generally agreed well with each other. In IOP-7, the wind is slightly south-east with wind velocity around 1.5 m.s -1 at height of 3 m. The positive values of L M O at all levels indicate that the measurement period is clearly stable. It has also been checked with the vertical temperature gradient. So this IOP corresponds well to the specific meteorological conditions required to perform the tracer experiment. The variances of the three wind components (σ 2 a , σ 2 b , σ 2 w ) have different order of magnitude which shows that the flow measured is strongly anisotropic. The ratios σ 2 a /u 2 * , σ 2 b /u 2 * and σ 2 w /u 2 * are consistent with the values found in [Stull, 1988, Arya, 1999] for similarity relationships in stable and neutral boundary layer. They also have different order of magnitude between three velocity components and their values are almost constant for anemometers at the same height.

If we look at the measurements at higher levels in Table 3.1, wind directions measured at 10 m and 30 m are much more north-east than that at height of 3 m, which is due to the effect of the forest to the north of the instrumented area, and will be discussed in Sect. 5.1 of Chapter 5 through a study of horizontal wind direction. The ratios σ 2 a /σ 2 w and σ 2 b /σ 2 w show relatively smaller values at 10 m and 30 m than From the measurements of IOP-7 at 3 m (Table 3.1), we can also observe the impact of the forest. The mean wind velocity measured by anemometers in the south is always greater than that measured in the north, and the mean wind direction also has a deviation between northern and southern anemometers. These are because the forest slows down and rotates the wind. Variances of three velocity components and friction velocities all have greater values for southern anemometers. These differences are small for σ 2 w but significant for σ 2 a and especially for σ 2 b , which can be explained by the subsidence and vertical diffusion of turbulence generated at the top of the forest by the north-easterly wind at higher levels (see simulation results in Sect.5.4 of Chapter 5). Similar characteristics due to the impact of the forest are found for all the IOPs with north-easterly wind at higher levels (10 m and 30 m). When the wind at higher levels is closer to the easterly wind, the differences between northern and southern anemometers are smaller. Also, the different magnitude between the variance of the three velocity components can be found for all the IOPs studied, which underscores the strong anisotropy of turbulence under stable conditions, very close to the ground.

Integral length scales

In order to understand more about the turbulence structure, we derive the integral time scales T e from the autocorrelation of three velocity components time series (see definition in Sect. 1.2.6 eq. 1.55), with approximation that it equals to the time lag where the velocity autocorrelation coefficient goes below 1/e for the first time (see example in Fig. 3.3). Then, we deduced the integral length scale L, which allow to characterize the dimension of most energetic eddies in the turbulent flows. Using Taylor's hypothesis, we can obtain the integral length scale from the simple relationship L = a mean T e as explained in Sect. 1.2.7 eq. 1.69, where a mean is the streamwise mean wind velocity.

The deduced integral length scales of three wind velocity components L aa , L bb and L ww can be interpreted as the spatial autocorrelation in the streamwise direction. Their values are reported in table 3.1. They have different order of magnitude with L aa >L bb >L ww , which indicates that, for the turbulence near the ground, espe- cially under stable conditions, there is a mix of large-scale structures and small-scale structures, and the small-scale structures seem to be dominant in the vertical direction. Integral length scales show also greater value with increasing altitude which implies that turbulence has larger eddies at higher levels. More discussion about variation of integral length scales as function of stability condition, measurement height and incident wind sector will be presented in Sect.3.3.

Velocity correlation

By analogy with the autocorrelation coefficient, one can define the crosscorrelation coefficient as

R a 1 a 2 (τ ) = a 1 (t)a 2 (t + τ ) a 2 1 a 2 2 (3.1)
with a 1 et a 2 representing the same turbulent variable at two different measuring points. In the data analysis, the cross-correlation coefficient is particularly calculated between pairs of anemometers (NE, NW), (SE, SW), (NW, 10N), (NW, 20N), (SW 10S) and (SW 20S), whose axis are close enough to the direction of the mean wind (easterly wind). As stated in [START_REF] Powell | An Investigation of the Application of Taylor's Hypothesis to Atmospheric Boundary Layer Turbulence[END_REF], the spatial cross correlation of velocity is a useful tool for examining the validity of Taylor's hypothesis of "frozen turbulence" and for determining the eddy advection velocity. Therefore, the calculation of velocity cross-correlation allows to demonstrate the transport of eddies between anemometers and the coefficient R a 1 a 2 becomes large when τ is in the same order of magnitude as the eddies advection time between two anemometers.

τ norm = τ t th = τ U dx cosθ (3.2)
where U is the average mean longitudinal wind speed between two correlated anemometers, dx is the distance between them, θ is the deviation of the mean wind direction from a direction parallel to their separation, and t th represents a theoretical time lag for eddy advection defined as t th = dx cosθ/U . The cross-correlation peaks reach up to 0.5 for the streamwise components, whereas the vertical components are poorly correlated (not shown). Theoretically, the cross-correlation peak should centre at τ norm = 1. However, they are all on the left of the vertical line at τ U/dx cosθ = 1. Defining the eddy advection velocity as

U adv = dx cosθ τ max (3.3)
where τ max is the time lag at the maximum correlation. The curves in Fig. 3.4 imply that U adv is higher than the mean wind speed U measured directly by the anemometers at the same level. 3.2 compares the mean wind speed and eddy advection speed deducted from spatial cross-correlation. A ratio of the eddy advection velocity to the mean wind speed is estimated as:

θ ( • ) U (m.s -1 ) U adv,a (m.s -1 ) U adv,b (m.s -1 )
r = U adv /U (3.4)
We note U adv,a and U adv,b as the eddy advection speeds deduced from the crosscorrelation of longitudinal and transverse velocity components respectively. Thus, we have r a = U adv,a /U and r b = U adv,b /U . Only anemometer couples with | θ |< 15 • have been considered to be useful to deduce the ratio r. Clearly, if the wind direction is not parallel to the line between the two anemometers, they do not sample exactly the same air parcels, thus decreasing the spatial cross correlation [START_REF] Horst | HATS: Field Observations to Obtain Spatially Filtered Turbulence Fields from Crosswind Arrays of Sonic Anemometers in the Atmospheric Surface Layer[END_REF]. So cross-correlations are not shown in Fig. 3.4 for anemometer couples (NE, NW) and (NW, 10N), because for them | θ |> 15 • . For other anemometer couples, U adv is always greater than U . Similar results were found in HATS (Horizontal Array Turbulence Study) field program presented in [START_REF] Horst | HATS: Field Observations to Obtain Spatially Filtered Turbulence Fields from Crosswind Arrays of Sonic Anemometers in the Atmospheric Surface Layer[END_REF], which found that the advection velocity may be up to 10%-20% larger than the mean wind speed obtained from measurements at height of 4 m and 8.3 m. [START_REF] Powell | An Investigation of the Application of Taylor's Hypothesis to Atmospheric Boundary Layer Turbulence[END_REF]] also deduced ratios with maximum values r = 1.15 for the three wind components. These values of the ratio r tend to show that Taylor's hypothesis of "frozen turbulence" might be not valid here, which can be due to the strong vertical velocity gradient in the surface layer very close to the ground: the eddy advection is probably affected by the flow at higher level where velocity is larger than that at height of 3 m. Moreover, measurements in [Powell andElderkin, 1974, Horst et al., 2004] are all under near-neutral stability conditions whose vertical shear is usually smaller than under a stable condition, which can explain the lower values of this ratio r compared to the values found here. The velocity cross-correlation at higher level for anemometers (10mSE, 10mSW) is very noisy for IOP-7, probably due to the north-easterly wind at higher altitude. More discussion about variation of ratios r a and r b as function of stability condition and measurement height will be presented in Sect.3.3.

Velocity power spectra

Power spectrum is another way to study turbulence structure. As first step, we compared the TKE power spectra with Kolmogorov's hypothesis (1.72) which implies the existence of an inertial subrange in the spectra. However, instead of following the -5/3 power law, a slope between -1 and -5/3 is found if we consider the whole spectral interval between the low frequency peak and 1 Hz. As mentioned in Sect.1.2.8, spectra studies are usually made on one-dimensional spectra of the three velocity components because of the different corresponding length scales between them. Therefore, we plotted here the power spectra of three velocity components for the measurement at heights of 3 m (Fig. 3.5(a)), 10 m (Fig. 3.5(b)) and 30 m (Fig. 3.5(c)). The spectra are calculated with a linear detrending procedure and premultiplied by frequency. The spectra shown in Fig. 3.5 are the average spectra of anemometers at the same level, so the spectra at height of 30 m and 10 m are much more fluctuating than those at height of 3 m because fewer sensors are available at these levels.

We observe that the spectra shapes correspond well to the description of [START_REF] Kaimal | Atmospheric Boundary Layer Flows -Their Structure and Measurement[END_REF], Drobinski et al., 2004, Carlotti and Drobinski, 2004] mentioned in Sect.1.2.8. The energy peaks can be identified in all the spectra. For the higher frequency side, there is the expected inertial subrange with f S(f ) ∝ f -2/3 for all the spectra. For the lower frequency side, there is the energy-containing range with f S(f ) ∝ f +1 which is less evident for the spectra of two horizontal velocity components at higher levels because they don't extent enough in low frequency region. Moreover, there is a frequency lag in spectrum shape between vertical velocity and horizontal ones showing strong anisotropy of turbulence. The vertical velocity spectrum is increasingly close to the horizontal ones with increasing height, which indicates a more isotropic turbulence at higher levels.

In our horizontal velocity spectra, the plateau at lower frequency range to the inertial subrange seems to have roughly a 0 slope which shows the existence of the f S(f ) ∝ f 0 subrange. [Richards et al., 1997, Hunt andMorrison, 2000, Figure 3.5 -Velocity power spectra f S(f ) at heights of (a) 3 m, (b) 10 m and (c) 30 m, calculated by averaging spectra of anemometers at the same level for the IOP-7 [START_REF] Högström | Theory and measurements for turbulence spectra and variances in the atmospheric neutral surface layer[END_REF] have all studied the velocity spectra shape from different field data and illustrated the existence of f S(f ) ∝ f 0 subrange for horizontal velocity spectra. [START_REF] Richards | Anisotropic turbulence in the atmospheric surface layer[END_REF], Högström et al., 2002] have found that the extension of this self-similar range decreases with increasing height. [START_REF] Högström | Theory and measurements for turbulence spectra and variances in the atmospheric neutral surface layer[END_REF] have also suggested that the upper limit of the 0 slope range f u decreases with increasing height. Unfortunately, our spectra at 10 m and 30 m are too noisy and don't extent low enough in frequency to see clearly the lower limit of the f S(f ) ∝ f 0 subrange for the horizontal velocity components. We can still find that the upper limit of this range roughly decreases with increasing height as pointed out by [START_REF] Högström | Theory and measurements for turbulence spectra and variances in the atmospheric neutral surface layer[END_REF].

Vertical velocity spectra have been shown to have a different form than horizontal velocity spectra, in all the literature mentioned above. At lower frequency, a range such that S ww (f ) ∝ u 2 * z ∼ constant has been proposed by [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF]. We find that our vertical velocity spectra respect well a slope of 1 corresponding to a range f S ww (f ) ∝ f prior to the power spectra peak, which is consistent with the typical form of vertical velocity spectra in the surface layer (Fig. 1.4 in Sect.1.2.8). The low frequency part of spectra seems to increase with height, as indicated in [START_REF] Högström | Theory and measurements for turbulence spectra and variances in the atmospheric neutral surface layer[END_REF] in the surface layer. However, some evidence of a f S(f ) ∝ f 0 subrange can also be found at an intermediate frequency range in vertical velocity spectra even at a height of 3 m. The origin of this phenomenon is still unknown, as all the theories have been made under near-neutral conditions, while our measurements are made under stable stratification which might cause a much thinner eddy surface layer. Thus, it is possible that our measurements are taken in a transition zone between the eddy surface layer and the shear surface layer above. Also, the forest and the shelters around might change in complex ways the flow by creating other shear layers and mixing them with the eddy surface layer.

Data analysis for two years of continuous measurements

Taking advantage of continuous measurements from ultrasonic anemometers, turbulence characteristics have been investigated with different stability conditions, incident wind sectors and measurement heights in order to shed light on the impact of terrain heterogeneity and the boundary-layer stability on flow characteristics such as wind velocity profiles, turbulence eddies length-scales and structures. Data analysis are made on two years of continuous measurements (April 2012 -March 2014) by calculating average statistical values and length scales for every 30min period. The statistical values calculated are the same as for an IOP (Table 3.1). Anemometers taken into account are the four anemometers at "sonic square" at 3-m height, and the three others on the two towers at height of 10 m and 30 m. Measurements from "sonic arc at 50 m" are not included because several anemometers were absent for months during these two years. After statistical analysis, variables values are selected for common periods of 30 min for all the anemometers and arranged according to measurement heights, different stability conditions and incident wind directions.

Influence of stability condition

A frequently used characterization of the atmospheric turbulence stability is the turbulence classification scheme originally developed by [Pasquill, 1961]. It categorized the atmospheric turbulence into six stability classes named A, B, C, D, E and F with class A corresponding to the extremely unstable cases, and class F to the moderately stable cases (Table 3.3). [Golder, 1972] then proposed a dependence relationship of Pasquill's stability classes on inverse Monin-Obukhov length 1/L M O and surface roughness. After [Golder, 1972], we deduce the range of Monin-Obukhov length for each stability class with a roughness length z 0 = 0.3 m which corresponds to the value estimated for the SIRTA site. [Golder, 1972] Before sorting 30-min periods into different stability classes, periods with northerly wind ( dd< 50 • or dd>330 • ) measured at height of 30 m are eliminated to avoid the impact of the forest. The impact of the forest with different wind sectors will be discussed in Sect. 3.3.2. We eliminate also those periods without valid values of length scales L aa , L bb and L ww calculated with the approximation method (equation 1.55 and 1.69). Then, we apply the selection ranges of L M O in table 3.3 for measurements at 10 m. L M O selection ranges for measurement at 3 m and 30 m are then deducted by keeping z/L M O constant. The number of final selected 30-min periods for each Pasquill's stability class and its percentage over the total period for anemometers SE, 10mSE and 30mSE are reported in table 3.4. The total selected periods number is 8465 and is identical for each anemometer. To facilitate the comparison between three measurement heights, we choose anemometers SE, 10mSE and 30mSE to present the results, because SE is relatively further from the forest to the north than other anemometers at 3 m, and 10mSE and 30mSE are collocated and thus easier to compare than 10mSW and 30mSE. Since stability classes A and B are less frequent, we combine them in order to have more representative sample compared with other stability classes. Table 3.4 -The final number of selected 30-min periods for each Pasquill's stability class and its percentage over the total number of selected periods for anemometers SE, 10mSE and 30mSE; The total number of selected periods is 8465 and is identical for each anemometer.

Figure 3.6 shows the averaged mean wind speed vertical profiles for different stability conditions. We can see that the profiles correspond well to the wind speed profile of flows in the surface-layer which varies approximately logarithmically with height. If we plot this profiles on semi-log graph (not shown), as suggested in [Stull, 1988], we obtain an almost straight line for the neutral class, and concave downward lines for more stable classes. The wind speed is greatest for the neutral condition (class D) and becomes smaller for more unstable classes and more stable ones. The velocity shear is stronger with increasing stability, which is consistent with the velocity profiles deduced from Monin-Obukhov's similarity relationships for different stability regimes [Garratt, 1992].

Figure 3.7 shows the variance of longitudinal and vertical velocity σ 2 a and σ 2 w profiles for different stability classes. We observe that σ 2 w is much smaller than σ 2 a which illustrates the strong anisotropic nature of the turbulence in the surface layer. σ 2 a and σ 2 w all increased with increasing altitude for all the stability classes due to stronger wind at higher altitude. However, it should be noted that the difference between σ 2 a and σ 2 w for different stability classes shown in Fig. 3.7 is influenced by the mean wind speed. In order to compare velocity fluctuation between stability classes, we discuss velocity standard deviation normalised by the median of the mean wind speed measured at 30 m for the same stability class σ a /U 30 and σ w /U 30 (Fig. 3.8). The result shows that the velocity fluctuation is smallest for the most stable class F due to the small wind speed and weak turbulence, and it is most important for the most unstable classes A-B due to the strong convection. Next, we also study the velocity fluctuation normalised by the friction velocity measured close to the ground (at 3 m) (Fig. 3.9). We found that σ a /u * 0 and σ w /u * 0 for the neutral class D are on the same order as the typical values proposed in [Stull, 1988, Arya, 1999, Garratt, 1992]. σ a /u * 0 and σ w /u * 0 have smallest values for the neutral condition and have greater values for more stable or more unstable conditions, which corresponds well to the curve of normalised wind fluctuation with stability parameter in [Garratt, 1992] (pp. 72 Fig. 3.16).

Figure 3.10 shows integral length scales of longitudinal and vertical wind components as function of measurement height for different stability classes. Integral length scales represent the size of largest eddies in a turbulent flow and are proved to be useful to characterize the anisotropy of turbulence. We can observe that integral length scales are increasing with height between 3 m and 30 m, which implies that turbulence structures are larger at higher level, as expected theoretically [START_REF] Carlotti | Length scales in wall-bounded high-Reynolds-number turbulence[END_REF] and previously observed in SIRTA field measurement and simulation results in [START_REF] Fesquet | Impact of terrain heterogeneity on coherent structure properties: numerical approach[END_REF].

Comparing integral length scales between wind components, we found L aa L ww for all the measurements heights and for all the stability conditions. Similar observations are also given in [Kaimal, 1973]. This difference between L aa and L ww can be linked to the mechanism suggested in [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF], which states that turbulence close to the ground is affected by a top-down mechanism. Large eddies impinge onto the ground where they generate internal boundary layers due to blocking (i.e. zero tangential velocity at the ground) in which smaller eddies develop. It seems that small structures are dominant for the vertical scales, while large structures are dominant for the longitudinal scales. This observation is also consistent with the simulation results in [START_REF] Fesquet | Impact of terrain heterogeneity on coherent structure properties: numerical approach[END_REF], which argued that the streamwise velocity is representative of large structures that impinge onto the ground through a top-down mechanism while the vertical velocity is more representative of small structures resulting from the impingement of the large structures.

Stability condition in the surface layer also has an impact on the turbulence structure. We observe that the length scales increase with increasing instability, which is consistent with the profiles presented in [Kaimal, 1973, Pena et al., 2010]. The increase is more remarkable for the measurements at higher levels (10 m and 30 m) where there are less effect of blocking by the ground. This indicates the inhibition of large-scale structure under stable stratification. Similar conclusions can be found in [START_REF] Barthlott | Long-term Study of Coherent Sturctures in the Atmospheric Surface Layer[END_REF], Fesquet et al., 2009a], where they all report the high probability of smaller scale eddies under stable conditions for their measurements at height of 30 m. Next, we would like to investigate the impact of stability on the validity of Taylor's hypothesis of "frozen turbulence" through the ratio r between eddy advection velocity and mean wind speed (equation 3.4 and 3.3). Variables such as the deviation of mean wind direction from anemometer axis θ, velocity correlations coefficient R and eddy advection velocity U adv are deduced for the common 30-min periods for anemometer couples (NE, NW), (SE, SW) and (10mSE, 10mSW). Then, periods with | θ |< 15 • and (R max -R min )/R max > 0.4 are selected to ensure a good quality in values of eddy advection velocity U adv and ratio r. The number of final selected 30-min periods for each Pasquill stability class and its percentage over the total selected period number for anemometer couples are reported in table 3.5. The total number of selected periods is 627 and is identical for each anemometer couple. In order to have more homogeneous data numbers, stability classes are merged to present de results as: unstable condition (classes A, B and C), neutral condition (class D) and stable condition (classes E and F). Table 3.5 -The final number of selected 30-min periods for each Pasquill's stability class and its percentage over the total number of selected periods for anemometer couples (NE, NW), (SE, SW) and (10mSE, 10mSW); The total number of selected periods is 627 and is identical for each anemometer couple.

Figure 3.11 shows variation of the ratio between eddy advection velocity estimated with the correlation coefficient and mean wind speed r a and r b as function of stability and measurement height. We observe that the ratio decreases with increasing height, and r a > r b for a given stability class and a given height, which is in agreement with [START_REF] Powell | An Investigation of the Application of Taylor's Hypothesis to Atmospheric Boundary Layer Turbulence[END_REF] and [START_REF] Horst | HATS: Field Observations to Obtain Spatially Filtered Turbulence Fields from Crosswind Arrays of Sonic Anemometers in the Atmospheric Surface Layer[END_REF]. The ratio is smaller with increasing instability, which can be explained by the smaller wind shear in more unstable condition, thus U adv deduced from correlation is closer to the measured wind speed U . This is also the reason why the ratio r found in the IOP-7 has greater values than those in neutral conditions in the literature. In general, we The results show that Taylor's hypothesis is not strictly valid, especially in the stable surface layer with strong wind shear. The truth is that the eddies change shape continuously as they pass by the instrument. Theoretically, Taylor's hypothesis can be expected to hold only for eddies small enough so that velocity gradients across them are negligible compare to the advection velocity. [Kaimal, 1973] has proposed a simple relationship which considers the changes of eddies small when:

σ a /U 2.0 (3.5)
However, it seems that for our measurements in stable conditions and very close to the ground, this criteria of [Kaimal, 1973] is not sufficient for the validity of Taylor's hypothesis.

Influence of wind sector

In order to study the impact of the forest to the north of the instrumented area, we decide to investigate two years of data classified by incident wind sectors.

According to the topography of the SIRTA site (Fig. 2.1) , we choose a north sector with impact of the forest, a west sector with open terrain on the upstream side, and a east sector with different small obstacles (shelters, barriers, etc.) on the upstream side. The range of wind direction for each sector are shown in table 3.6. Data selection is made on the wind direction measured at 30 m (by 30mSE) averaged over 30-min period.

Incident wind sector

Description

Range of wind direction ( Table 3.7 -The final number of selected 30-min periods for each incident wind sector in neutral condition (class D) and its percentage over the total number of selected periods for anemometer SE, 10mSE, 30mSE.

After the wind sector selection, 30-min periods without valid values of length scales L aa , L bb and L ww are eliminated, and then divided into six Pasquill stability classes (table 3.3) like what we did in the previous section. For the comparison between different wind sectors, only selected 30-min periods in neutral condition (class D) are taken into account. The number of final selected 30-min periods for each incident wind sector in class D and its percentage over the total selected periods number for anemometer SE, 10mSE and 30mSE are reported in table 3.7. We obtain much less easterly wind than other wind sectors due to the weather condition of north-western France. Figure 3.12 shows the mean wind speed profiles for different incident wind sectors in neutral condition. As expected in north-western France, westerly wind is usually stronger than wind from other sectors. We observe that the wind profiles of east and west sectors are similar with the stronger shear close to the ground, while wind profile of north sector shows a much stronger deceleration in wind speed than others. Since the forest to the north of the instrumented area has a height of 15 m, it is clearly the forest which slows down the wind from the north sector. This deceleration of the flow across a vegetation zone has been also reported from data analyses in [START_REF] Irvine | The evolution of turbulence across a forest edge[END_REF], Thomas and Foken, 2007, Chahine et al., 2014].

Figure 3.13 and 3.14 present respectively the profiles of variances and integral length scales of longitudinal and vertical wind components for different incident wind sectors in neutral condition (class D). Like the wind profile, the behaviour of northerly wind is quite different than easterly and westerly wind due to the impact of the forest. The velocity variances σ 2 a and σ 2 w both have much greater values for the north sector than the others, while the integral length scales L aa and L ww both have smaller values for northerly wind. The difference is particularly obvious for the measurements at height of 10 m and 30 m, which is due to the large wind shear at the top of the forest. A strong turbulent kinetic energy with smaller scale eddies is created at the canopy top and advected downstream by the northerly wind to the instrumented area. Similar profiles of velocity fluctuations and length scales can be found in [START_REF] Irvine | The evolution of turbulence across a forest edge[END_REF], which shows the increase of velocity standard deviation and decrease of velocity length scales for the flow across a forest edge. More discussion about the forest effect on the turbulent flows in SIRTA Zone 1 will be presented in the numerical study in Chapter 5.

Moreover, the behaviours between easterly and westerly wind sectors are quite similar. Thus, from the measurements, it is difficult to identify the impact of the shelters and barriers on the east side of the instrumented area.

Conclusions

This Chapter has analysed the wind and turbulence data of the SIRTA field campaign. The analyses are made for a specific IOP, and for two years of continuous measurements.

The IOP-7 is an IOP with stable stratification and a north-easterly wind. Turbulence strong anisotropy in a stable surface layer is quantified by different order of magnitude between variances σ 2 a , σ 2 b and σ 2 w . Integral length scales deduced from the spatial autocorrelation in the streamwise direction also show different order of magnitude as L aa > L bb > L ww . The turbulent structures advection speed obtained by velocity cross-correlation between sensors at 3 m above the ground is higher than the measured average wind speed at this height, which confirms that the turbulence in the stable surface layer is affected by strong wind shear. Different spectral forms have been found between vertical and horizontal velocity components for this stable IOP. The spectral form seems coherent with the typical spectra in the eddy surface layer described in [START_REF] Drobinski | The structure of the near-neutral atmospheric surface layer[END_REF] with a f S(f ) ∝ f 0 subrange at intermediate frequency range for the horizontal velocity components. As for the vertical velocity spectra, evidence of this self-similar subrange can also be found, which is maybe due to a much thinner eddy surface layer under the stable condition or to perturbations from the forest and the shelters around.

Turbulence characteristics for different stability conditions and incident wind sectors have been investigated through analysis of two years of data. Mean wind speeds are larger under neutral conditions, while velocity fluctuations are larger under unstable conditions. Velocity fluctuations are very weak under the most stable conditions due to small wind speed and stratification. Integral length scales have higher values in a less stable surface layer and are increasing with height, which shows that turbulence in the surface layer exhibits larger structures with increasing height and instability. The different order of magnitude between integral length scales of different velocity components shown as L aa L ww can be found for all stability classes due to the blocking effect of the ground. This difference indicates that flows in the SIRTA Zone 1 are affected by a mixture of small scale and large scale structures, which are dominant for the vertical and longitudinal length scales respectively. Analysis of the velocity cross-correlation shows that the ratios of the eddy advection velocity to the mean wind speed r = U adv /U have slightly higher values under more stable conditions and for measurements closer to the ground. These results indicate that the ratios are affected by 2 factors: the wind vertical shear and the eddies vertical size, which are linked to atmospheric stability and to height above the ground level. The wind shear is more important under a more stable surface layer, and eddies have smaller size with decreasing height. Moreover, the study of different wind sectors illustrates that the forest to the north has a strong impact on the flow in Zone 1. It slows down the wind below its height (at 3 m and 10 m). The strong shear at the forest top induces larger wind velocity fluctuation and smaller eddy structures which are transported by the northerly wind and detected by the anemometers at 10 m and 30 m.

Chapter 4

Flow and dispersion modelling

The second part of my work is to perform numerical simulations with a CFD code, Code_Saturne, co-developed by CEREA and Electricité de France (EDF), and to verify if the code is able to reproduce the characteristics of flows and dispersion at the SIRTA site with many different meteorological conditions, by comparing with experimental data set.

Modelling of pollutants dispersion includes modelling of turbulent flow. This chapter will first introduce different approaches and models for turbulent flow modelling and dispersion modelling. Next, it will give a general presentation of the CFD code, Code_Saturne, and with more detail, a summary of the equations applied in the simulations.

Turbulent flow modelling

Different approaches

The particular properties of turbulent flows make it difficult to develop an accurate and tractable theory or model [Pope, 2000]. The velocity field is threedimensional, time dependent and quasi-random. There is a large range of time scales and length scales with the integral length scale L representative of the largest turbulent motions and the Kolmogorov's microscal η representative of the smallest motions down to order of 0.001 m. Moreover, as shown in Sect.1.2.3, the equations system is non-linear. Thus, there are many different numerical methods which allow the modelling of turbulent flows with supplementary hypotheses.

In general, three broad strategies are commonly used: Direct Numerical Simulation (DNS). In DNS, Navier-Stokes equations are solved for a complete time-dependent velocity field U (x, t), which means from fine-scale dissipative motions to large-scales overall flow. According to [Kolmogorov, 1941], the ratio between large-eddy scale and Kolmogorov's microscale L/η is proportional to Re 3/4 , where Re is the Reynolds number of the flow. Since all length-scales and time-scales have to be solved in DNS, for a three-dimensional resolution, the computational cost increases as Re 9/4 . In the atmosphere, the commonly encountered range of Re is from 10 5 to 10 8 , so Re 9/4 is approximately from 10 11 to 10 18 , which means that DNS is computationally very expensive. Nowadays, appearance of 'supercomputer' makes calculation in DNS feasible in some cases. However, this approach is still considered as too expensive and time consuming comparing with others methods especially for simulations with large domains. Thus, DNS is usually restricted to flows with low to moderate Reynolds number [Pope, 2000]. Large-Eddy Simulation (LES). In LES, the largest three-dimensional unsteady turbulent motions are directly represented, whereas the effects of the smallest-scale motions are modelled. The velocity field U (x, t) is decomposed into a filtered component U (x, t) which is representative of the largest-scale turbulent motions, and a residual component u (x, t) which is representative of the smallest-scale (or subgrid-scale (SGS) ) motions. The equations for the evolution of the filtered velocity field are derived from the Navier-Stokes equations with the momentum equation containing the residual-stress tensor (or SGS stress tensor). Closure for the equation system consists of modelling the residual-stress tensor with subgrid-scale models. The simplest and most commonly used SGS models are eddy-viscosity models, for example, the Smagorinski Model (see detailed expression in [Arya, 1999, Pope, 2000[START_REF] Dall'ozzo | Modélisation d'écoulements atmosphérique stratifié par Large-Eddy Simulation à l'aide de Code_Saturne[END_REF]). Finally, the model filtered equations are solved numerically for U (x, t), which provides an approximation to the large-scale motions in one realization of the turbulent flow. LES lies between DNS and RANS with dynamics of the largest-scale motions computed explicitly and the influence of the smallest-scales represented by simple SGS models. It is expected to be more accurate than RANS and less expensive than DNS, and is currently the best available approach to atmospheric turbulence and diffusion modelling [Arya, 1999]. However, the grid size and the SGS model should be chosen carefully according to the expected size or scale of the most energetic eddies in the flow in order to obtain desired accuracy in results. Also, there is an important difficulty in imposing the proper boundary conditions, which must include large eddy structure.

Reynolds Averaged Navier-Stokes (RANS). As presented in Chapter 1

Sect.1.2.3, RANS approach is about to decompose variables into an ensemble mean part and a perturbation part and deduce a system of Reynolds-averaged equations (1.14) to (1.17). This system of equations is solved for the mean velocity field but with the Reynolds stress u i u j and other similar terms Θ u j and c i u j appeared as unknowns. These unknown terms can be determined by a turbulence model: the most widely employed in RANS are turbulentviscosity models, such as the turbulence closures (1.20) -(1.21) -(1.22) for u i u j , Θ u j and c i u j presented in Sect.1.2.3; another model that is often used is a second-moment closure derived directly from Reynolds-stress transport equations [Pope, 2000, Hanjalić andLaunder, 2011] which is usually called the 'R ij ' model.

Although RANS is expected to give less accurate results than LES and DNS, it is still an approach that is commonly used in atmospheric science and many other industrial domains because of its lower computing requirements.

In the simulation part of my work, since it is the first numerical study for the SIRTA flow and dispersion campaign, we are particularly interested in RANS approach with turbulent-viscosity model as the turbulence closure. As this approach is still largely used for industrial applications, it is important to first evaluate its performances. In the future, RANS approach with 'R ij ' model and LES will be applied to give a more accurate comparison between measurements and simulation results.

Turbulent-viscosity models

Turbulent-viscosity models are based on the turbulent-viscosity hypothesis (Sect.1.2.3). According to this hypothesis, second order terms u i u j , Θ u j and c i u j can be expressed as (1.20) -(1.21) -(1.22) with appearance of a new term -turbulent viscosity ν t . To solve the system of equations, the remaining question is to determine an appropriate specification of ν t , which can be written as the product of a velocity scale u and a length scale l.

The simplest turbulence model for ν t is the mixing-length model, which is also called zero-equation model. The expression was first introduced by Prandtl in 1925 which specifies ν t as a function of geometry and flow parameter:

ν t = l 2 m | ∂u ∂z | (4.1)
where l m is the mixing length representing the length scale. In the log-law region, l m is usually estimated as l m = κz. The velocity scale is presented as a simple shear flow l m |∂u/∂z|. The major drawback of this simple model is that the mixing length has to be specified, but it is dependent on the geometry of the flow. The specification of l m requires a large measure of guesswork, and consequently there is little confidence in the accuracy of the results [Pope, 2000].

Another less empirical model for ν t is the k -l model, also called one-equation model or turbulent-kinetic-energy model. [Kolmogorov, 1942] suggested that it is better to base the velocity scale on the turbulent kinetic energy: u = C kl k 1/2 . The length scale is again taken to be the mixing length, and the turbulent viscosity becomes

ν t = C kl l m k 1/2 (4.2)
where C kl is a constant. The value of k is obtained by solving a transport model equation for k. This one-equation model has a modest advantage in accuracy over mixing-length model [Pope, 2000]. However, the major disadvantage remains that the length-scale l m must be specified.

The most widely used complete turbulent model and integrated in many CFD codes is the k -ε model, also called two-equation model. In this model, two turbulence quantities k and ε are solved from transport model equations, and the flow-dependent specification of l m is not required. From these two quantities, one can form a length scale proportional to k 3/2 /ε, and a time scale τ = k/ε. Thus, according to the standard k -ε model proposed in [START_REF] Launder | The numerical computation of turbulent flows[END_REF], the turbulent viscosity can be expressed as:

ν t = C µ k 2 /ε (4.3)
where C µ is one of the five model constants.

The standard model transport equations for k and ε are [START_REF] Launder | The numerical computation of turbulent flows[END_REF]:

∂k ∂t + u j ∂k ∂x j = ∂ ∂x j ν t σ k ∂k ∂x j + ν t ∂u i ∂x j + ∂u j ∂x i ∂u i ∂x j -ε (4.4) ∂ε ∂t + u j ∂ε ∂x j = ∂ ∂x j ν t σ ε ∂ε ∂x j + C ε1 ε k ν t ∂u i ∂x j + ∂u j ∂x i ∂u i ∂x j -C ε2 ε 2 k (4.5)
with the values of the constants as follows:

4.2. Dispersion modelling 83

C µ C ε1 C ε2 σ k σ ε 0.09 1.44 1.92 1.0 1.3 Table 4
.1 -Values of the constants in the standard k -ε model [START_REF] Launder | The numerical computation of turbulent flows[END_REF] The k -ε model is arguably the simplest complete turbulence model and has the broadest range of applicability. We have also chosen this model for our numerical study of the SIRTA campaign. It should be noted that the values of the standard k -ε model constants represent a compromise chosen to give the 'best' performance for a range of flows. For any particular flow, it seems that the accuracy of the model calculations can be improved by adjusting these constants [Pope, 2000]. The modified k -ε model and the values of its constants adapted for the atmospheric boundary-layer flows will be presented later in Sect.4.3.3.

Dispersion modelling

Different approaches

Like the turbulent flow modelling, there are many different approaches for the dispersion modelling. Here, I will introduce briefly several approaches that are usually employed for dispersion within the atmospheric boundary layer.

Gaussian models are simple dispersion models which consider that pollutants mean concentration has a Gaussian distribution during dispersion. They can model dispersion from a continuous or an instantaneous point source. Assuming a continuous point source in a uniform flow with horizontal homogeneous turbulence, the mean diffusion equation can be expressed as:

u ∂c ∂x = K y ∂ 2 c ∂y 2 + K z ∂ 2 c ∂z 2 + Q c δ( - → x -- → x 0 ) (4.6)
with σ y = 2K y x/u and σ z = 2K z x/u, where K y and K z are the turbulent diffusivities in the y and z directions respectively, Q c is the source emission rate, u is the mean transport velocity across the plume, and σ y and σ z are the standarddeviations of mean concentration distribution along the y and z directions.

The Gaussian plume formula can be resolved as:

c(x, y, z) = Q c 2πuσ y σ z exp - y 2 2σ 2 y - z 2 2σ 2 z (4.7)
If the continuous source is at height of H above the ground, assuming a perfectly reflecting surface, the Gaussian plume formula becomes:

c(x, y, z; H) = Q c 2πuσ y σ z exp - y 2 2σ 2 y exp - (z -H) 2 2σ 2 z + exp - (z + H) 2 2σ 2 z (4.8)
σ y and σ z , which are the model's parameters, can be seen as the representative length scales of plume lateral and vertical diffusion, and are dependent on characteristics of turbulent flows. The empirical estimation of σ y and σ z was first made by [Pasquill, 1961, Gifford, 1961] for the various Pasquill stability classes. Then, [Briggs, 1973] proposed empirical values adapted for open country and urban areas. More recently, [START_REF] Hanna | A baseline urban dispersion model evaluated with salt lake city and los angeles tracer data[END_REF] extended the empirical scheme of Briggs for situations of low wind speeds and large turbulence intensities in urban areas, and validated the model's performance with data set of two dispersion experiments (Salt Lake City Urban 2000 and Los Angeles 2001). Since Gaussian models are conceptually appealing and computationally cheap to use, they are extensively applied in assessing the impact of existing and planned sources of air pollution on local and urban air quality, particularly for regulatory application [Arya, 1999]. However, due to the approximation involved, they are limited to idealized situations such as relatively flat and homogeneous surfaces, relatively short distances (< 50 km), strong enough wind and moderately stable or unstable or neutral conditions. Moreover, Gaussian models have large uncertainties due to natural variability and simplified physics. Thus, more sophisticated numerical models should be used to simulate dispersion in complex flow conditions, for example, strong wind shear, inhomogeneous turbulence, presence of obstacles (e.g. buildings, forest) and very stable or convective conditions. By applying a Lagrangian approach in dispersion modelling, one can obtain better performance in pollutants dispersion prediction. The Lagrangian approach, which is to study the evolution of the fluid particles, consists in computing particles trajectories and properties from their previous positions and velocities. Imagine a particle in a turbulent flow at position (x , y , z ) at instant t , the probability of finding this particle at position (x, y, z) at a later instant t is φ(x, y, z, t)dxdydz. Thus, for N particles with different probability density function φ i (x, y, z, t), mean concentration at position (x, y, z) at instant t given by the Lagrangian approach can be expressed as [Milliez, 2006]:

c(x, y, z, t) = N i=1 φ i (x, y, z, t) (4.9)
Lagrangian particle-dispersion models are usually coupled with an Eulerian model or a stochastic model which provides the mean flow and turbulent velocity fields for calculation of particles trajectories. Although the Lagrangian approach can provide excellent description of each particle dispersion, three-dimensional mesoscale Lagrangian modelling generally require large computer resources, because a small grid size is needed to adequately resolve complex flows and a large number of particles must be considered to simulate dispersion from each source.

Eulerian CFD models are commonly employed in local scale dispersion modelling because they are accurate at small scale due to their ability to represent in detail the effects of relief and obstacles. [Carissimo andMacdonald, 2002, Milliez andCarissimo, 2008] have modelled dispersion in an idealized urban area with CFD models in RANS k -ε approach. [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF] have also compared the simulation results with the MUST data and shown an overall good agreement. [START_REF] Dejoan | Comparison between large-eddy simulation and reynoldsaveraged navier-stokes computations for the must field experiment. part ii: effects of incident wind angle deviation on the mean flow and plume dispersion[END_REF] have modelled the MUST field experiment with both LES and RANS and shown that the mean concentration profiles predicted by LES and RANS are both in good qualitative agreement with the measurements, and the fluctuation intensity of the concentration obtained from LES gives a satisfactory agreement with the experiment as well. [Gamel, 2015] has evaluated the validity of RANS k -ε approach by comparing with a wind tunnel experiment of turbulence flow and dispersion around an bi-dimensional obstacle. He has particularly analysed the impact of turbulent parameters on the simulation results especially for the values of C µ , σ k and Sc t .

The CFD approach is also the method that I will apply for the dispersion modelling here. The pollutant concentration will be modelled using the equation of conservation of a scalar quantity (1.11). Mean concentration will be calculated with the transport equation (1.17) deduced in Sect. 1.2.3. Concentration fluctuations modelling will be introduced in the next section.

Concentration fluctuation modelling

Many pollutants dispersion models, for example the Gaussian models, are mainly concerned with the mean concentration field because they may be directly compared with observed time-mean concentrations. However, as [Csanady, 1967] had pointed out, such a mean field is not sufficient to supply adequate information about the possible effects of the diffusing substance on living organisms, because instantaneous maxima, the persistence of given concentration and its occurrence frequency are all important problems when we study the hazardous effects of some pollutants or define the uncertainty in air quality models. Also, the concentration fluctuations can be usually as large as the mean concentration. Thus, concentration fluctuation is an important factor in the pollutants dispersion modelling.

There is already a wide range of approaches that have been developed for the concentration fluctuation modelling: (1) Models based on a transport equation for the concentration variance (K models) [Csanady, 1967, Andronopoulos et al., 2002, Sato and Sada, 2002, Hsieh et al., 2007, Milliez and Carissimo, 2008]; (2) Empirical Gaussian models [START_REF] Wilson | Predicting the spatial distribution of concentration fluctuations from a ground level source[END_REF]; (3) Similarity models [Mole, 2001]; (4) Meandering plume models [Gifford, 1959, Fackrell and Robins, 1982a, Fackrell and Robins, 1982b, Yee and Wilson, 2000];

(5) Probability density function models (PDF) [Mylne andMason, 1991, Lewis andChatwin, 1997]; (6) Statistical models (Lagrangian stochastic models) [START_REF] Cohen | Parameterization of the two-point velocity correlation function in a two-particle lagrangian stochastic model and its effect on the prediction of concentration variances due to a line source[END_REF]]; (7) Large-eddy simulation models [START_REF] Xie | Large-eddy simulation of dispersion: comparison between elevated and ground-level sources[END_REF], Dejoan et al., 2010].

The concentration fluctuations measured at a receptor can be characterized by the standard deviation of concentration: σ c . Here, we apply the type (1) for the concentration fluctuation modelling in our numerical simulations, which is also suitable for the CFD approach. A discussion about the performance of the Eulerian CFD model in simulating the near-source dispersion was made by [START_REF] Demael | Comparative evaluation of an eulerian cfd and gaussian plume models based on prairie grass dispersion experiment[END_REF].

The transport equation for the concentration fluctuations σ c is derived from the transport equation for the instantaneous concentrations c (1.11) and the transport equation for the mean concentrations (1.17). With the Reynolds decomposition of u (1.12a) and c (1.12d), multiplying the transport equations (1.11) and (1.17) by 2c , subtracting two equations and taking the average, we can obtain the transport equation for the concentration fluctuations σ 2 c = c 2 = (c -c) 2 :

∂c 2 ∂t 1 + u j ∂c 2 ∂x j 2 = ∂ ∂x j D m ∂c 2 ∂x j - ∂ ∂x j c 2 u j 3 -2c u j ∂c ∂x j 4 -2D m ∂c ∂x j ∂c ∂x j 5 (4.10)
with 1: temporal variation 2: advection by the mean velocity 3: molecular and turbulent diffusion 4: production term 5: dissipation term, also defined as ε c

As presented in Sect.1.2.3 for the other turbulent flux u i u j (1.20) , Θ u j (1.21) and c u j (1.22), the turbulent flux of concentration variance can also be modelled by using Boussinesq's hypothesis such that:

c 2 u j = - ν t Sc t ∂c 2 ∂x j (4.11)
With (1.22), the term 4 becomes:

-2c u j ∂c ∂x j = 2 ν t Sc t ∂c ∂x j ∂c ∂x j (4.12)
As for the term 5, [Csanady, 1967] proposed a simple expression :

ε c = c 2 T c (4.13)
where T c is dissipation time scale that is characteristic of the decay time of the concentration fluctuation scalar field. Different expressions for T c can be found in previous papers. For example, [Fackrell and Robins, 1982a, Andronopoulos et al., 2002, Hsieh et al., 2007] used description as follows:

T c = L c k 1/2 (4.14)
where L c is the dissipation length scale (or, the characteristic length scale for the concentration fluctuations) that is a function of the local characteristics of the plume. [START_REF] Sato | Estimation of gas concentration fluctuations using a numerical model and comparison with those of wind tunnel experiments[END_REF], Hsieh et al., 2007, Milliez and Carissimo, 2008, Gamel, 2015] employed another expression which assumes that the ratio of the dissipation time scale T k of the turbulent kinetic energy k to the dissipation time scale T c of the concentration variance c 2 is a constant R f :

T c T k = c 2 /ε c k/ε = R f (4.15)
[ [START_REF] Hsieh | Numerical modeling of passive scalar dispersion in an urban canopy layer[END_REF] also compared simulation results between these two models and found that the first model (4.14) provided better agreement with the experimental data than the second model (4.15). However, in the first model, the dissipation length scale L c is explicitly linked with the local plume scale, which is physically associated with the characteristic length scale of eddies that are comparable to the instantaneous plume size and are responsible for the internal mixing in the plume. These links make the first model much more complicated than the second one.

In my thesis, we choose the simpler model (4.15) for the dissipation of concentration fluctuations like [Milliez andCarissimo, 2008, Gamel, 2015], who had performed numerical simulations by using the same CDF code. Thus, the term 5 is expressed as:

2D m ∂c ∂x j ∂c ∂x j = ε c = 1 R f ε k c 2 (4.16)
Finally, with equations (1.22) (4.11) (4.12) and (4.16), the transport equation for the concentration fluctuations (4.10) becomes:

∂c 2 ∂t 1 + u j ∂c 2 ∂x j 2 = ∂ ∂x j D m + ν t Sc t ∂c 2 ∂x j 3 + 2 ν t Sc t ∂c ∂x j ∂c ∂x j 4 - 1 R f ε k c 2 5
(4.17)

Code_Saturne

General description

Numerical simulations are performed by using an open source CDF code, Code_Saturne, developed within Electricité de France (EDF) since 2001. The code is designed for laminar or turbulent flows in complex two-and three-dimensional geometries with complex physics. The set of equations considered consists of the Navier-Stokes equations for incompressible flows completed with equations for turbulence modelling and additional scalars (e.g. temperature, enthalpy, concentration of species, etc.). The code is based on the finite volume co-located method. The system of governing equations is discretized at the level of the volume cells throughout the flow domain. Then this system of non-linear algebraic equations obtained after the discretization is solved in an iterative way. As for the turbulent flow modelling, two approaches are available in Code_Saturne: RANS and LES (See Sect.4.1). And for RANS approach, many turbulence closures can be used: mixing-length model,

k -ε model, R ij -ε model, etc.
Code_Saturne has been used for industrial applications and research activities in many fields related to energy production, such as nuclear power thermal-hydraulics, gas and coal combustion, turbomachinery, heating, ventilation and air-conditioning, etc.

The numerical method and some examples of application are described in detail in [START_REF] Archambeau | Code_saturne: a finite volume code for the computation of turbulent incompressible flows -industrial applications[END_REF].

Atmospheric module

An atmospheric module within Code_Saturne has been co-developed by CEREA and EDF R&D, and is particularity adapted for the simulations of flows and pollutants dispersion within the atmospheric boundary layer. The main industrial applications are in wind-energy potential assessment, and pollutant dispersion and deposition around power plants. In this module, we add a thermodynamic variable, the potential temperature Θ. Code_Saturne then solves the conservation equation for Θ (1.10) together with all the other governing equations. The vertical variation of air density is also calculated according to the variations of pressure and potential temperature. Moreover, the analytical profiles deduced from Monin-Obukhov similarity theory (Sect.1.2.5) are used to define the boundary conditions.

There are many studies carried out with the atmospheric module of Code_Saturne. [START_REF] Demael | Comparative evaluation of an eulerian cfd and gaussian plume models based on prairie grass dispersion experiment[END_REF] have compared the performance of CFD simulation with two Gaussian plume models on the dataset of the Prairie Grass dispersion field experiment. [Milliez and[START_REF] Milliez | [END_REF]Carissimo, 2008] have studied the turbulent flows and pollutants dispersion in neutral and stable ABL for the MUST campaign. [START_REF] Qu | Micrometeorological modeling of radiative and convective effects with a buildingresolving code[END_REF], Qu et al., 2012] have evaluated the thermal effects of buildings on the local atmospheric flow with a coupled dynamic-radiative model. [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF] have simulated the flow over a complex semi-urban terrain (whole domain of Ecole Polytechnique) including all the SIRTA measurement zones for different wind sectors and with a canopy model integrated to simulate forest effect. [START_REF] Chahine | Modelling atmospheric effects on performance and plume dispersal from natural draft wet cooling towers[END_REF] have assessed the performance of a cooling tower by modelling the plume formation and its dispersal under realistic atmospheric conditions. All these numerical studies have been compared with different field measurements and have shown in general satisfactory agreement between them.

Equations

In this section, I would like to summarize all the equations applied for the numerical study of SIRTA experiment in my thesis. Simulations are run in RANS mode with the standard k -ε turbulence model modified for atmospheric flows. We assume that the flow is anelastic, and the Coriolis force is neglected. The governing equations with Boussinesq's hypothesis (Sect.1.2.3) for the mass, the momentum, the heat and the pollutants concentration become: ∂ρ u i ∂x i = 0 (4.18)

ρ ∂u i ∂t + u j ∂u i ∂x j = - ∂p ∂x i + ∂ ∂x j (µ + µ t ) ∂u i ∂x j + ∂u j ∂x i - 2 3 ρ ∂ ∂x j k -ρg i + S u,i (4.19) ρ ∂Θ ∂t + u j ∂Θ ∂x j = ∂ ∂x j λ C p + µ t P r t ∂Θ ∂x j (4.20) ρ ∂c ∂t + u j ∂c ∂x j = ∂ ∂x j D + µ t Sc t ∂c ∂x j (4.21)
where S u,i is an additional momentum sink term modelling the effect of the forest (See Sect.5.2.3), D is the molecular diffusion coefficient, and the extra sources terms for the heat and the concentration are neglected.

For the concentration fluctuations, we have:

ρ ∂c 2 ∂t + u j ∂c 2 ∂x j = ∂ ∂x j D + µ t Sc t ∂c 2 ∂x j + 2 µ t Sc t ∂c ∂x j ∂c ∂x j - ρ R f ε k c 2 (4.22)
With the k -ε turbulence model, the turbulent viscosity µ t is expressed as:

µ t = ρC µ k 2 /ε (4.23)
where we choose C µ = 0.03 for atmospheric flows according to the work of [Duynkerke, 1988, Katul et al., 2004, Zaidi et al., 2013]. The transport equations for the turbulent kinetic energy k (4.4) and the dissipation rate ε (4.5) become:

ρ ∂k ∂t + ρu j ∂k ∂x j = ∂ ∂x j µ t σ k ∂k ∂x j + P + G -ρε + S k (4.24) ρ ∂ε ∂t + ρu j ∂ε ∂x j = ∂ ∂x j µ t σ ε ∂ε ∂x j + C ε1 ε k (P + C ε3 G) -C ε2 ρ ε 2 k + S ε (4.25)
where S k and S ε are extra source terms representing the turbulence generation and destruction due to interaction of the mean flow motion with the vegetation elements [START_REF] Katul | Oneand two-equation models for canopy turbulence[END_REF], Zaidi et al., 2013] (See Sect.5.2.3), P and G are the production of k by the mean wind shear and the production or destruction rate due to buoyancy respectively:

P = µ t ∂u i ∂x j + ∂u j ∂x i ∂u i ∂x j (4.26) G = - µ t P r t 1 Θ ∂Θ ∂x i g i (4.27) C µ C ε1 C ε2 σ k σ ε 0.03 1.44 1.92 1.0 1.88 Table 4
.2 -Values of the constants in the k -ε model modified for atmospheric flows [START_REF] Detering | Application of the E-ε turbulence model to the atmospheric boundary layer[END_REF] The model constants also take the values modified for atmospheric flows following [START_REF] Detering | Application of the E-ε turbulence model to the atmospheric boundary layer[END_REF] (Table 4.2). And the value of C ε3 is taken after [Viollet, 1988]: C ε3 = 0 for a stably stratified atmosphere (G < 0), C ε3 = 1 for an unstably stratified atmosphere (G > 0).

Chapter 5

Numerical study for mean and turbulent flow

After a detailed study of turbulence characteristics from field data analysis (Chapter 3), we have discussed the impact of terrain heterogeneity and stability condition on the three-dimensional flow at SIRTA. In this, chapter, we are interested in simulating the inhomogeneous flow and turbulence over the SIRTA site and comparing with measurements. The simulation is performed with the CFD code, Code_Saturne (Sect.4.3). The goal is to analyse the ability of the atmospheric module in Code_Saturne to accurately reproduce the effects of the micro-scale heterogeneities, especially the impact of forest, on the mean flow and turbulence field at SIRTA, for neutral and stable conditions.

Many previous numerical studies can inspire our work in modelling and simulation. [START_REF] Milliez | Numerical simulations of pollutant dispersion in an idealized urban area, for different meteorological conditions[END_REF] have performed simulations with Code_Saturne in RANS mode, with k -ε turbulence closure in neutral and stable conditions and showed good overall agreement for dynamics comparison with the MUST data. [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF] have made simulations of flow over the whole domain of École Polytechnique (including SIRTA measurement zones) with Code_Saturne with k -ε turbulence closure and a canopy model [Sanz, 2003, Katul et al., 2004] for the forest area in neutral condition. Comparing with measurements averaged by direction from 3 sonic anemometers on 10-and 30-m masts, simulation results are consistent with measurements for normalized horizontal velocity and turbulence kinetic energy. [START_REF] Dalpé | Numerical simulation of wind flow near a forest edge[END_REF] have also performed simulations with a k -ε turbulence model and the canopy model of [Sanz, 2003] for different forest density distributions and found significant impact of the drag coefficient and forest density on simulation results within and above the forest. Large-Eddy Simulation (LES) has also been applied to study canopy flow under neutral stratification by [START_REF] Shaw | Large-eddy simulation of turbulent flow above and within a forest[END_REF][START_REF] Dupont | Influence of foliar density profile on canopy flow: a large-eddy simulation study[END_REF], Dupont and Brunet, 2008a[START_REF] Fesquet | Impact of terrain heterogeneity on coherent structure properties: numerical approach[END_REF]. [Shaw andSchumann, 1992, Dupont and[START_REF] Dupont | Influence of foliar density profile on canopy flow: a large-eddy simulation study[END_REF] studied the sensitivity of turbulent flow to the canopy morphology and showed that the typical features of the canopy flow become more pronounced as canopy density increases. [Dupont and Brunet, 2008a] illustrated different regions of turbulent flows over forest edges and showed that the flow adjusts faster with a denser canopy. [START_REF] Fesquet | Impact of terrain heterogeneity on coherent structure properties: numerical approach[END_REF] have characterised turbulence coherent structures close to the surface and their dependence on surface heterogeneity by using LES simulations compared to SIRTA measurements.

Sections 5.1, 5.2, 5.3 and 5.4 of this chapter are extracted from the article [START_REF] Wei | Experimental and numerical study of wind and turbulence in a near-field dispersion campaign at an inhomogeneous site. Boundary-Layer Meteorology[END_REF] submitted in Boundary-Layer Meteorology. A wind rotation effect due to the impact of terrain heterogeneity will be first investigated from 2 years of data. Then, after an introduction of numerical procedures and forest model used in simulations, discussions of the simulation results and comparison between simulations and field measurements will be presented. At last, a section (Sect. 5.5) will be added for an additional sensitivity study for the impact of the forest with a simplified geometry.

Impact of terrain heterogeneity

Site heterogeneity plays an important role for mean flow and turbulence in Zone 1. Since the measurement zone is oriented along the east-west axis with a forest canopy on the north side, wind flows are particularly affected in a large northerly sector by this forest [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF]. Figure 5.1 plots the frequency distribution of mean wind direction for continuous measurements over two years (April 2012 -March 2014) at three levels: 3 m (NE and SE), 10 m (10mSW), and 30 m (30mSE). The mean wind direction is deduced from 10-min averages of 10 Hz measurements. The forest height h is about 15 m. The frequency distribution of anemometer 30mSE shows the prevalence of south-westerly winds (around 230 • ) at the site, while this distribution is almost constant for all the other directions. As for lower levels, some new peaks appear around 90 • and 270 • . These new peaks are rather small for anemometer 10mSW but are much more pronounced for the measurements at 3-m height. Moreover, frequency distributions of anemometers NE and SE show very few winds for a large northerly sector (0 • -60 • and 310 • -360 • ). It appears that northerly winds have turned to easterly and westerly winds at a height of 3 m. This phenomenon is clearly due to a wind channelling effect from the forest to the north of the instrumented area.

In order to see how the wind turns at different levels, measurements over two years are grouped according to wind direction measured by the anemometer 30mSE. Figure 5.2 plots the mean wind direction of anemometers NE, SE, 10mSW and 30mSE for the groups defined by direction (30mSE) between 40 • -50 • (Fig. 5.2a) and 315 • -325 • (Fig. 5.2b). We can see clearly that for a north-easterly wind at 30 m, wind direction turns to easterly at 3 m (Fig. 5.2a), while for a north-westerly wind at 30 m, it turns to westerly at 3 m (Fig. 5.2b). There is also a slight difference between two anemometers at 3-m height. Wind direction given by the anemometer NE turns to become more parallel to the east-west axis than the SE and shows a larger spread, which is due to the NE's smaller distance to the forest edge. Also, with wind slowed down by the forest, direction variability becomes larger. As for measurements at 10-m height, wind direction begins to turn parallel to the east-west axis, but not as much as those at 3 m. Since the forest height is h = 15 m, the wind channelling effect is found mostly for the measurements below this height. The direction at 30 m remains unaffected.

Inside Zone 1, there are some shelters that are less than 3 m high on the east side of the instrumented area. However, the impact of these shelters on the turbulence has not been identified in the measurements (i.e. data from the closest anemometers in the "sonic square" that have different positions relative to the shelters). Moreover, according to [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF], who have made simulations of the flow over the whole domain of École Polytechnique (including other SIRTA measurement zones) for all the wind sectors (every 10 • ) and compared with the 10-m measurement in Zone 1, it is found that the production of turbulence is dominated by the effect of the northern canopy. Turbulent kinetic energy has higher values for northerly sectors which correspond to the forest and is nearly constant for the rest of wind directions. Therefore, in the present study, we suppose that the building to the south-west of Zone 1 and the urban area on the other side of the lake are sufficiently far away that they do not affect the instrumented area of Zone 1.

Numerical modelling

Numerical simulations are performed with Code_Saturne, an open source CFD code co-developed by CEREA and Électricité de France (EDF). The code is designed for laminar or turbulent flows with complex geometry and complex physics and is based on the finite volume co-located method. The set of equations consists of Navier-Stokes equations for incompressible flows completed with equations for turbulence modelling and additional scalars. The numerical method is described in detail by [START_REF] Archambeau | Code_saturne: a finite volume code for the computation of turbulent incompressible flows -industrial applications[END_REF]. Atmospheric specificities are described by [START_REF] Demael | Comparative evaluation of an eulerian cfd and gaussian plume models based on prairie grass dispersion experiment[END_REF] and [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF].

Simulation domain

The simulation domain is an area of dimension 1600 m (west to east) × 700 m (north to south) × 200 m (vertical) in Zone 1 (see modelling area in Fig. 5.3), with a progressive three-dimensional mesh refined in the instrumented area and near the ground. Figure 5.4 shows the mesh for the modelling area. The horizontal resolution ranges from 1 m in the instrumented area (180 m × 100 m) to 5 m for the rest of the modelling area. The vertical grid resolution varies from 0.5 m near the ground up to 10 m at 200-m height. A land-use file provided by the French National Institute of Geographic and Forest Information (IGN, "Institut national de l'information géographique et forestière") and manually corrected on the basis of satellite photographs is used to identify the different land cover types (forest, low vegetation, water, etc.) and assign to each type a corresponding roughness length z 0 . Shelters located eastward of the instrumented area are also taken into account explicitly in the mesh.

Boundary conditions

The boundary conditions in the present numerical study are constant in time and defined as follows:

-Inlet condition: Dirichlet type with analytical profiles. We apply Monin-Obukhov similarity theory and expressions proposed by [START_REF] Dyer | Flux gradent relationships in the constant flux layer[END_REF], Businger et al., 1971, Hicks, 1976] for the wind speed and temperature profiles and [START_REF] Kerschgens | Comments on turbulence parameters for the calculation of dispersion in the atmospheric boundary layer[END_REF] for the turbulence parameter profiles such as TKE and dissipation rate (see introduction in Sect. 1.2.5). During our simulations, analytical profiles are generated from a given wind speed at a reference height, a given wind direction, a given Monin-Obukhov length L M O and an aerodynamic roughness length z 0 . For the simulations of IOPs, due to the lack of measurements at the entrance of the simulated domain, we use wind direction and speed measured at 30mSE to generate the inlet profiles because winds receive less modification from the environment at this height and are more alike the winds at the entrance. Moreover, as we are interested in the layer very close to the ground, we use the value of L M O measured at 10mSE during the IOP for the generation of inlet profiles.

-Outlet condition: free outflow.

-The ground and shelters surfaces: according to the land use cover map, a constant roughness z 0 = 0.0001 m is selected for the lake and z 0 = 0.3 m is selected for the low vegetation area. This relatively high value is justified by the presence of various small obstacles like barriers. The forest is modelled by the canopy model introduced in the next section. Since there are large instruments on the roofs of the various shelters, a constant roughness z 0 = 0.3 m is also selected for all the shelter surfaces.

Canopy model

Concerning the modelling of the forest area, a canopy model described by [Sanz, 2003, Katul et al., 2004, Dalpé and Masson, 2009, Zaidi et al., 2013] has been applied for the present simulation. The model uses a porous media analogy combined with a modified k -ε turbulence model to simulate momentum losses and turbulence generation within the forest. The additional source terms in the Navier-Stokes momentum equation and the k and ε equations create a decrease in wind speed and modify turbulence.

The added source term in the Navier-Stokes momentum equation S u,i (See eq. 4.19) represents the momentum absorbed within the forest, while outside the forest S u,i = 0 [START_REF] Katul | Oneand two-equation models for canopy turbulence[END_REF]. It can be expressed as

* S u,i = -ραC d |U |U i (5.1)
where ρ is the air density, C d is the drag coefficient of the forest, α is the leaf area density, |U | is the local mean velocity magnitude, U i is the wind velocity for coordinate i. In the present study, the drag coefficient C d is set to 0.2 as proposed by [Katul et al., 2004, Dalpé and[START_REF] Dalpé | [END_REF]. The leaf area density α represents square meters of leaf surface per unit volume. Taking into account the type of trees in Zone 1 of SIRTA site (leaves are present from bottom to the top of the trees), the leaf area density is assumed constant with height. After [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF], the average value of α over the whole year is 0.5 m -1 , and α = 0.9 m -1 during summer and α = 0.25 m -1 during winter for the forest at SIRTA. The effect of the forest on TKE is modelled with additional terms S k and S ε in the k and ε equations (See eq. 4.24 and 4.25). They both consist of two terms, the first one corresponding to the turbulence generation due to breakage of the mean flow motion by vegetation, and the second one corresponding to the short-circuiting of the turbulence cascade. They are expressed as follows * [Sanz, 2003, Katul et al., 2004]:

S k = ραC d β p |U | 3 -ραC d β d k|U |
(5.2)

S ε = ραC d C ε4 β p ε k |U | 3 -ραC d C ε5 β d ε|U | (5.3)
where C ε4 , C ε5 , β p and β d are the model constants. The values of these constants are, according to propositions by [Sanz, 2003, Katul et al., 2004, Dalpé and Masson, 2009]:

β p = 1 (5.4) β d = C µ 2 0.05 2/3 β p + 3 σ k = 5.03 (5.5) C ε4 = C ε5 = σ k 2 σ ε -C µ /6 2 0.05 2/3 (C ε2 -C ε1 ) = 0.9
(5.6) * . As pointed out by S. Dupont, rapporteur of the thesis, Equations 5.1, 5.2 and 5.3 differ from the formula found in [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF]] by a 0.5 factor which was erroneously included in that paper. After verification we have found that this 0.5 has also been used in the numerical simulations presented here, giving an effective leaf area density α which is too small by a factor of 2.

Simulations under neutral conditions and comparison with 2-year measurements 101

A sensitivity study performed by [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF] showed that the values of C ε4 , C ε5 , β p and β d have little impact on simulation results.

Simulations under neutral conditions and comparison with 2-year measurements

Simulations are first conducted in RANS mode with a standard k -ε turbulence model under neutral conditions. The goal is to verify that Code_Saturne is able to reproduce correctly the ensemble mean flow in the experimental area. Simulations are performed with 4 typical wind directions: north-east (45 • ), north-west (320 • ), east (90 • ), and west (270 • ), in order to see how well the code can reproduce the impact of the forest on the wind rotation at low level and on the turbulence. The inlet conditions are generated as explained in Sect. 5.2.2 with wind speed ff = 2 m.s -1 at 10 m, with roughness length z 0 = 0.3 m and with turbulence at equilibrium. The forest is simulated with the average leaf area density α = 0.5 m -1 .

Figures 5.5 and 5.6 show the simulated vertical profiles of direction deviation dddd 30 , normalized horizontal wind speed f f /f f 30 and normalized turbulent kinetic energy T KE/f f 2 30 . The subscript '30' means values of variables at 30 m. The wind rotation is due to the impact of the forest to the north, which has been discussed for the measurements in Sect. 5.1 and can be observed clearly in the vertical profiles of dd -dd 30 for the 45 • and 320 • cases. The wind turns to be more parallel to the forest edge with decreasing height. For the simulation with easterly (90 • ) and westerly (270 • ) wind, there is no such directional shear because the wind is already parallel to the forest edge. The wind direction remains unchanged above the forest for all the cases, which is consistent with the discussion in Sect. 5.1. The normalized horizontal wind profiles (f f /f f 30 ) show wind deceleration below the forest height (h = 15 m) especially for the 45 • and 320 • cases. From the normalized turbulent kinetic energy profiles (T KE/f f 2 30 ) of the north-easterly and north-westerly wind cases, we notice that turbulence reaches a maximum slightly above the top of the canopy (around 20 m) due to the large vertical wind shear, and then decreases with decreasing height. Similar effects for wind speed and TKE have been discussed by [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF].

Differences between northern and southern anemometers at 3 m can also be distinguished in the simulation profiles. The simulations with 45 • and 320 • winds show a larger effect of the forest on the profiles of anemometers NE and NW, with a larger deviation in wind direction, a stronger wind deceleration and a smaller turbulent kinetic energy below forest height than other anemometers. These smaller TKE values can be related to the vertical advection and diffusion processes downstream of the forest (See Sect. 5.4). Even for the simulation with easterly wind (90 • ), the vertical profiles corresponding to NE and NW anemometer locations show smaller wind speed than others below forest height. Thus, these two locations are still affected by the forest, which can be explained by the shape of its edge. For the same reason, they exhibit greater values of the normalized TKE above 3 m than at south anemometer locations (SE and SW), where little difference with the inlet profile is seen.

Simulation profiles are compared with continuous measurements over 2 years for selected wind directions: 40 Overall, we find that the simulations are able to reproduce well the main characteristics of the flow observed in the measurements. The change in wind direction below the forest height is underestimated by simulations for both northwesterly (320 • ) and north-easterly (45 • ) winds. However, for the case of easterly (90 • ) winds, simulations and measurements agree well. In the case of westerly winds (270 • ), there is a slight wind rotation in the measurements, but not in the simulation results. The reason is that the forest geometry has been extracted from a satellite map and the land-use map with uncertainties of at least 5 m. [START_REF] Zaidi | Numerical simulations of the microscale heterogeneties of turbulence observed on a complex site[END_REF] already showed that the comparison of calculated TKE with measurements was strongly improved at this site by increasing the resolution of the land-use map. However, an uncertainty of 5 m on the position of the forest edge can still have a significant impact on the results. Since wind rotation effects are very sensitive to the shape of the forest and the distance to it, a difference of 5 m could induce such a difference between simulations and measurements. Also, it is possible that the average leaf area density α = 0.5 m -1 is not sufficient to model the dense forest to the north; this point will be discussed further in Sect. 5.4. The wind speed decrease is clearly overestimated by the model, especially for the 45 • and 320 • cases, where there is a strong impact of the forest. This could be due to the fact that the boundary conditions are generated from analytical profiles, which might underestimate the wind speed very close to the ground. Nevertheless, the difference between northern and southern anemometers is correctly reproduced. The normalized TKE is well reproduced for the north-easterly (45 • ) case, and is slightly overestimated for the 5.4. Simulations under stable conditions and comparison with IOP-7105 north-westerly (320 • ) case, but remains within the measurement variability. The overestimation shown for the easterly (90 • ) and westerly (270 • ) cases might again be due to the uncertainty in the forest shape and the analytical inlet profiles.

Simulations under stable conditions and comparison with IOP-7

In order to confirm that Code_Saturne can also simulate realistically a stably stratified atmosphere, simulations are performed in RANS mode with standard k -ε turbulence model and stable thermal stratification. Inlet conditions are generated with measurements of anemometers 10mSE and 30mSE during the IOP on 5 June 2013: wind direction dd = 58 • , wind speed f f = 3.5 m.s -1 , L M O = 130 m and z 0 = 0.3 m.

Figures 5.7 and 5.8 show simulation results of pressure, wind speed and TKE in horizontal and vertical cross-sections for Zone 1 and a zoom on the instrumented area. The simulation is run with a leaf area density α = 0.9 m -1 , because the IOP was in summer. Simulation results are consistent with the phenomena observed in measurements during the IOP and discussed in Sect. 3.2. Figure 5.7b shows the strong heterogeneity of the wind speed at 3-m height and in particular the decrease of wind speed due to the forest. We can see more clearly in Fig. 5.7d the effect on the speed in the instrumented area and the change of direction, which leads to a divergent flow in this area. The flow tends to follow the forest edge, which might be partly explained by the pressure gradient due to the low pressure region within the forest area (Fig. 5.7a). Figures 5.7c and 5.7e show that TKE is very low inside and close to the forest, but then increases towards the south because of the wind shear. It leads to a strong TKE gradient in the instrumented area. This is in agreement with the relatively high TKE values given by anemometers SE, SW and 20S compared to the northern anemometers (see Table 3.1 in Chapter 3). The impact of the shelters on TKE seems to be dominated by this effect of the forest area.

Figure 5.8 exhibits the vertical cross-section for TKE and vertical velocity w. A positive vertical velocity is found in front of the forest edge (on the northern side), and a negative one is found behind the forest (on the southern side). A region of strong turbulence above the canopy is observed, which is induced by the large wind shear at the canopy top. Maximum TKE occurs around 1.5 to 2 times the tree height and decays rapidly inside the forest as previously observed by [ [START_REF] Shaw | Large-eddy simulation of turbulent flow above and within a forest[END_REF], Dupont et al., 2006, Dupont and Brunet, 2008a]. The development of this turbulent region coincides with the decrease in vertical velocity as observed by [Dupont and Brunet, 2008a]. These high values of TKE are advected downstream of the forest and impact the instrumented area. Due to the vertical advection and diffusion, these high values reach the layer close to the ground, which explains the higher TKE measured and calculated for the southern anemometers.

Figure 5.9 compares vertical profiles between measurements and simulations with two leaf area densities: α = 0.9 and 2 m -1 , for the case of IOP-7. This comparison is presented only for two anemometers at 3 m (20N and 20S), which are the closest and furthest anemometers from the forest edge, and for another 2 anemometers at 10 m and 30 m (10mSE and 30mSE) on the tower in Zone 1. Measurements used for the comparison are statistical values averaged over the sub-period of 60 min during IOP-7 given in Table 3.1 in Chapter 3. Similar characteristics can be found when comparing with the simulation of north-easterly winds in Sect. 5.3. The rotation in wind direction and the wind deceleration are stronger with decreasing altitude and with decreasing distance from the forest edge. A maximum TKE can also be seen around the forest canopy height in the simulated profiles. TKE has smaller values with decreasing distance to the forest. Statistics of the simulations performance are shown in Table 5.1. The errors are calculated with the measurements at 3 m and are considered as the indicator of the Figure 5.9 -Comparison of vertical profiles of (a) horizontal wind direction, (b) horizontal wind speed, (c) TKE and (d) vertical velocity between simulations and measurements for IOP-7 ; 'dashed line': inlet condition; 'solid line': simulations with α = 0.9; 'dash-dot line': simulation with α = 2; 'colour circles': measurements. simulation performance for the individual IOP. Comparing with the measurements of the IOP on 5 June (IOP-7), the simulation under stable conditions reproduces well the characteristics of the flow. However, as for the neutral case in Sect. 5.3, it underestimates the wind rotation and overestimates the wind deceleration below the forest canopy height. TKE has been underestimated in the simulation of this IOP. The underestimation is about 0.18 m 2 .s -2 at 3 m, but is much greater at 30 m. As for the vertical velocity w, since the instrumented area is behind the forest with a north-easterly wind, negative vertical velocities can be found for the simulated profiles and the measurements which are consistent with Fig. 5.8b. The simulation underestimates the vertical velocity by 0.02 to 0.06 m.s -1 below the forest canopy height and by 0.08 m.s -1 at 30 m. The measured vertical velocity shows greater 5.4. Simulations under stable conditions and comparison with IOP-7109 values near the forest edge at 3 m. The simulation does not show this difference at 3 m but at higher levels, around the forest canopy height. As for the neutral cases, the underestimation of the wind rotation and the vertical velocity might be due to the uncertainty in the geometry of the forest. The overestimation of the wind deceleration and the underestimation of the TKE might be due to the inlet profiles which are calculated analytically using the measurements inside the domain. Particularly, these inlet profiles are generated with a vertically uniform stability. However, in reality, the stability condition can vary with height. During the IOPs performed at the beginning of the night, the measurements at 30 m usually show a less stable stratification than that at 3 m. Thus, the analytical inlet profiles might underestimate the wind velocity very close to the ground and overestimate the effect of stability at 30 m.

Performance statistics for simulations of two other IOPs are also shown in Table 5.1. The IOP on 20 November 2014 (IOP-10) is under near-neutral conditions with north-easterly winds. The wind rotation and deceleration are comparable to the IOP on 5 June but are less important, because the wind is more easterly. The simulation shows better agreement with the measurements for TKE and less difference in wind rotation. The IOP on 11 March 2015 (IOP-11) is under stable conditions and with slightly south-easterly winds. Since there is almost no impact of the forest, the difference between the northern and southern anemometers in wind direction and speed are less important in the measurements. TKE and the vertical velocity have also smaller variations between 3, 10 and 30 m. The simulation agrees better with the measurements in the horizontal wind speed and direction compared to the simulation for the IOP on 5 June (IOP-7). Thus, the IOP on 5 June 2013, which is under stable conditions and with a strong impact of the forest, is the most complex and challenging situation for the numerical study.

Figure 5.9 also shows simulation profiles with two different leaf area densities α = 0.9 and 2 m -1 . A leaf area density of 2 m -1 is used as the maximum value in the leaf area density vertical profile in the simulations of [START_REF] Dalpé | Numerical simulation of wind flow near a forest edge[END_REF]. As discussed previously, a constant value of α is taken here which is justified by the tree variety. Even if this value remains on the high side, it allows us to test the sensitivity of simulation results to this parameter. While the wind direction above the forest remains similar between α = 0.9 and 2 m -1 , the wind rotation below the forest height increases with increasing leaf area density, as mentioned by 5.1 -Performance statistics for the simulations of three IOPs all with z 0 = 0.3 m and α = 0.9 m -1 ; calculations are made with 9 anemometers at 3 m; Bias =

IOP Inlet condition Bias RMSE MRE dd( • ) at 30 m f f (m/s) at 30 m L M O (m) at 10 m dd( • ) f f (m/s) TKE (m 2 /s 2 ) w (m/s) dd( • ) f f (m/s) TKE (m 2 /s 2 ) w (m/s) f f (m/s) TKE (m 2 /s 2 ) w (m/
1 n (x s -x o ), RMSE = 1 n (x s -x o ) 2 , and MRE = 1 n 2|x s -x o |/|x s + x o |,
with x s and x o representing respectively the simulation results and the measurements, and n = 9.

Simulations under stable conditions and comparison with IOP-7111

Sensitivity study of IOP-7

Bias RMSE MRE dd( • ) f f (m/s) TKE (m 2 /s 2 ) w (m/s) dd( • ) f f (m/s) TKE (m 2 /s 2 ) w (m/s) f f (m/s) TKE (m 2 /s 2 ) w (m/s)
Inlet condition with analytical profiles , z 0 =0.3 m -14.7 -0.34 0.17 Table 5.2 -Same as Table 5.1 for the sensitivity studies of the IOP on 5 June 2013 (IOP-7).

a. Note: It might appear that the simulation under neutral condition gives better results in f f and TKE. However, it gives worse result in wind direction dd, which is a very important factor for dispersion simulation.

[ Dupont and Brunet, 2008b] † . TKE profiles show little difference between both leaf area densities. The vertical velocity increases with increasing leaf area density as also shown by [Dupont and Brunet, 2008a]. In general, simulation results improve with increasing leaf area density, especially for the wind direction at 3 m.

Besides the leaf area density, we have also studied the impact of the inlet stability condition and the roughness length (z 0 ) on the simulation results for the IOP on 5 June (IOP-7). Performance statistics for all the sensitivity studies are shown in Table 5.2. With a neutral condition, simulation results are particularly improved for the wind speed and TKE. When decreasing the roughness length to z 0 = 0.1 m over the entire modelling domain, the simulation shows stronger wind rotation but underestimates even more TKE. Nevertheless, we found that these impacts are limited compared to those of the leaf area density especially for the wind rotation effect. We believe that the impact of the forest is the most important factor in accurately reproducing the flows at this site.

Lastly, as mentioned previously, differences between calculated and measured values might be due at least partly to the inlet profiles, which could not be measured and, therefore, must be built with theoretical profiles and measurements inside the simulation domain. Even at 30 m, these measurements at the site might be different from their upwind values. An attempt has been made to run additional simulations using a linear interpolation of the measurements at the site (wind, temperature and TKE obtained with radiometer, LIDAR and SODAR) as inlet profiles. The measurements taken into account are from 10-m and higher levels, since we would like to let the simulations generate the flow close to the ground. Figure 5.10 compares vertical profiles between measurements and simulations with two leaf area densities: α = 0.9 and 2 m -1 . Statistical errors between the simulation results and the measurements are also shown in Table 5.2. Regarding the underestimation of the wind rotation effect, we have closer values for the wind horizontal and vertical velocity, and even an overestimation of TKE around and below the forest canopy height. We conclude that an adequate inlet profile is another key factor in simulating the flow over a complex site, especially under stable stratification. In the future, this question about the inlet profiles might be addressed by performing data assimilation of on-site measurements as proposed by [START_REF] Zajaczkowski | A preliminary study of assimilating numerical weather prediction data into computational fluid dynamics models for wind prediction[END_REF].

Figure 5.10 -Same as Fig. 5.9 for simulations with the inlet profiles obtained from the interpolation of the measurements at the site.

Sensitivity study for the impact of the forest

Since we are aware of the strong impact of the forest on the wind measurement and on the simulation results, a sensitivity study is performed next. First, more wind rotation effect studies are added like in Sect.5.1 by selecting different periods of data. Then, the impact of leaf area density on different turbulence variables are investigated from simulations with a simplified forest geometry.

Study on wind directions

Like in Fig. 5.1, we plot the mean wind direction frequency distribution from two years of data but separating the neutral cases (Fig. 5.11(a)) from the stable cases (Fig. 5.11(b)). From the distribution of anemometer 30mSE, we observed that the prevailing wind direction (230 • ) has higher frequency in neutral conditions, while easterly and westerly winds are more frequent in stable cases. Thus, for anemometers at lower levels (10mSW, NE and SE), the peaks around 90 • and 270 • seem to have greater frequency in stable conditions. As for the large northerly sector (0 • -60 • and 310 • -360 • ), the frequency distribution of measurements at 3 m is almost zero for neutral cases, but has a small value for stable ones. Similarly to Fig. 5.2 for the neutral cases, we plot here for the stable cases the mean wind direction of anemometers NE, SE, 10mSW and 30mSE for the groups defined by mean wind direction (measured by 30mSE) between 40 • -50 • (Fig. 5.12(a)) and 315 • -325 • (Fig. 5.12(b)). The wind directional shear for anemometers below the forest height and the slight difference between two anemometers at 3-m height can be also observed in stable conditions (Fig. 5.12) like previously in the neutral ones in Sect.5.1 (Fig. 5.2). The median values of mean wind direction at 10 m (10mSW) show little difference between stable and neutral cases. Meanwhile, the rotation effect seems to be stronger at height of 3 m for the case of 315 • -325 • (Fig. 5.12(b)). Moreover, because of small wind speed in stable conditions, wind directions show larger spread in Fig. 5.12 than in Fig. 5.2, particularly for the anemometer NE with the deceleration of wind speed due to the forest. This also explains the non zero frequency for the large northerly sector in stable conditions (Fig. 5.11(b)). There- fore, it seems that the wind rotation effect might be slightly stronger under stable conditions for the measurements very close to the ground. However, it is difficult to draw a global conclusion whether the stability really plays a role in the wind rotation effect due to the relatively small data number of stable cases and their large spread in wind direction.

Next, we would like to see the impact of forest density on the wind rotation from measurements, just like we have varied the leaf area density α in the numerical simulation. The mean wind direction frequency distributions from two years of data are plotted for the months of summer (June, July, August) and the months of winter (Dec., Jan., Feb) in Fig. 5.13(a) and 5.13(b). The big difference between distributions during summer and winter is the peaks around 90 • and 270 • for the anemometers below the forest height. For the two anemometers at 3 m, those peaks are very narrow in winter but have much wider range in summer, which indicates that the wind rotation affects more periods in summer than in winter. Those peaks are rather hard to find for the anemometer 10mSW during winter which means the wind rotation effect is small at this height. Figure 5.14 shows the mean wind direction of anemometers NE, SE, 10mSW and 30mSE for the groups of mean wind direction (measured by 30mSE) between 40 • -50 • and 315 • -325 • during the months of summer and winter respectively. As expected, wind rotation is more important during summer than winter, since the forest density is higher in summer. The difference between anemometer NE and SE is more obvious for the summer months because of the stronger impact of the forest. These findings from field measurements are consistent with the simulation results of different leaf area densities.

Simulation with simplified geometry

In order to simplify the complex terrain and to isolate the forest effect, the sensitivity study is made on a simplified geometry. The simulation domain has the same dimension as in the modelling area presented in Sect.5.2.1: 1600 m (west to east) × 700 m (north to south) × 200 m (vertical), with a progressive three-dimensional mesh refined near the ground (from 0.5 m near the ground until 10 m at 200 m). Instead of modelling the real topography, we simply model a parallelepipedic forest area in the center of the domain. The parallelepipedic forest area has approximately the same dimension as the forest to the north of Zone 1: 1050 m (west to east) × 165 m (north to south) × 15 m (vertical). As shown in Fig. 5.15, the horizontal The main objective of the sensitivity study is to investigate the impact of forest density on wind and turbulence. Therefore, simulations are performed with various values of the leaf area density α = 0.2, 0.5, 0.9, 2 and 10 m -1 , from realistic values during winter and summer or averaged over the year, to the extreme values available in the literature. Moreover, in order to observe better the differences between a flow across a canopy and an obstacle, we also replace the parallelepipedic forest area by a solid obstacle of the identical dimension.

The boundary conditions applied are as follows:

-Inlet condition is deduced from the analytical profiles presented previously in Chapter 1 Sect.1.2.5. The incident wind is the same as the IOP-7 with wind direction = 58 • and wind velocity = 3.5 m.s -1 at height of 30 m. Simulations for the comparison between obstacle and forest with various values of α are all made in neutral condition. -Outlet condition is free outflow as usual.

-The ground and obstacle surfaces are simulated with an identical roughness length z 0 = 0.1 m which is a typical value for low vegetation ground. To study the impact on turbulent flows, vertical profiles of different variables are extracted before, inside and after the simplified forest or obstacle. Positions and names of the extracted points are also shown in Fig. 5.15.

Results and discussions

Vertical profiles of horizontal mean wind direction dd, mean wind speed f f , TKE and vertical velocity w are shown in Fig. 5.16, 5.17, 5.18 and 5.19 respectively, for various values of the leaf area density and the solid obstacle case, and for several extracted points. The vertical cross-sections of TKE, vertical velocity w and pressure fields are exhibited in Fig 5.21,5.22 and 5.23 for two different forest density and the obstacle case. The horizontal cross-section at 3-m height of pressure filed is also presented in Fig. 5.24.

The wind directional vertical shear can be observed in front of (Fig. 5.16(a)) and behind (Fig. 5.16(b)) the center of the forest canopy especially for the flow below its height. Consistently with the simulation results presented in [Dupont and Brunet, 2008b] ‡ , the wind rotation increases with decreasing height and increasing canopy density. The rotation is more important behind the forest.

For the horizontal mean wind speed, the deceleration is stronger with a higher forest density in front of the canopy (Fig. 5.17(a)) and within the canopy (Fig. 5.17(c) and 5.17(d)) (no wind below the canopy height for the obstacle case). For the profiles at point 'entrance' and 'exit', a strong gradient around the forest top can be observed. This gradient is enhanced with increasing canopy density as stated in [START_REF] Dupont | Influence of foliar density profile on canopy flow: a large-eddy simulation study[END_REF]. Within the canopy, the deceleration of the flow is stronger with increasing density as observed in [Dupont and Brunet, 2008a]. At 'exit', the wind speed within the canopy becomes almost zero for the extreme value case of α = 10 m -1 as in the obstacle case. Behind the center of the canopy (Fig. 5.17(b)), the trend is difficult to observe since there is little difference between these cases. However, it is interesting to point out that below 10 m, the wind seems to become slightly stronger with increasing density, which might be due to the fact that the wind circumvents the dense forest, like around a solid obstacle.

The TKE has different profiles before and after the encounter with the forest edge. At point 'ext 2', with increasing forest density, the TKE becomes stronger above the forest height and weaker below the forest height (Fig. 5.18(a)). At point 'center' and 'ext 3', the strong TKE region induced by the wind shear above the forest top can be observed (Fig. 5.18(b) and 5.18(c)). The TKE peak value seems to decrease with increasing forest density for the cases with leaf area density α = 0.2, 0,5 and 0.9 m -1 . Within the forest, the TKE decreases more rapidly with increasing density as observed in the simulations of [Dupont and Brunet, 2008a].

The vertical velocity is positive in front of the forest zone (Fig. 5.19(a)) and negative behind it (Fig. 5.19(b)). The vertical velocity is stronger with higher forest density due to the stronger flow distortion, and is maximal with the obstacle. At 'ext 3', w changes sign for the flow very close to the ground because of the small recirculation created behind the forest, which can be observed in Fig. 5.22. These observations are consistent with the previous results in Sect.5.4 and in [Dupont and Brunet, 2008a].

The dissipation rate ε in front of the forest zone (at point 'ext 2') has little difference between different forest densities (Fig. 5.20(a)). After the encounter with the forest, at point 'center' and 'ext 3', the dissipation rate has greater value with increasing forest density above the forest height (Fig. 5.20(b) and 5.20(c)), which can explain the decrease of the TKE peak value with increasing density. Then, the dissipation rate decreases rapidly below the forest height. We would also like to study the different behaviour between flows across a canopy and across a solid obstacle by replacing the forest area by an obstacle with the identical dimension. We observe that the flow across an obstacle seems like a flow across an extremely dense canopy. The horizontal wind directional shear is stronger for the obstacle case than for the forest case with α = 10 m -1 (Fig. 5.16). The horizontal wind deceleration is more important in front of the obstacle area (Fig. 5.17(a)), and the gradient is greater at the top of the area (Fig. 5.17(c) and 5.17(d)) because of zero wind speed inside the obstacle. At 'ext 3', the obstacle gives a greater wind speed than the forest cases (Fig, 5.17(b)), which can be related to the horizontal pressure gradient (Fig. 5.24(c)). As for the TKE, the obstacle case shows the strongest TKE around the obstacle height at point ' ext 2' (Fig. 5.18(a)), which This sensitivity study has investigated in particular the wind rotation effect as a function of forest density. We can observe from the vertical profiles (Fig. 5.16) and from the cross-sections (Fig. 5.23 and 5.24) that the greater is the leaf area density, the stronger the wind flow tends to follow the forest, with the obstacle as the extreme case. The change in wind direction can be considered as a change in velocity components. As shown in the Navier-Stokes momentum equation (4.19), we have verified that: outside the forest or obstacle zone, the wind advection term ρu j (∂u i )/(∂x j ) varies with the same magnitude as the pressure-gradient force term ∂p/∂x i ; within the forest, the term ρu j (∂u i )/(∂x j ) varies with the same magnitude as the sum of the pressure-gradient force term ∂p/∂x i and the forest aerodynamic drag term S u,i . The low pressure region is more remarkable with increasing forest density, which gives a stronger pressure-gradient force. Also, the aerodynamic drag S u,i increases with increasing α. The obstacle can be considered as the forest with infinite aerodynamic drag which prevents the wind flow from entering. Therefore, the wind rotation is more important with increasing forest density which is consistent with the observations from field data and simulation results.

In addition, we have compared the difference between stable and neutral conditions. For simulation under stable condition, we take L M O = 131 m as for the IOP-7. Similarly to the measurements, we find little difference in wind rotation between stable and neutral condition for identical forest density.

In summary, the obstacle can be considered as an extremely dense forest with α → ∞. The same as the simulation of the IOP-7 with two different leaf area densities, we can conclude that the features of canopy flow become more pronounced as canopy density increases. The extracted point 'ext 3' is approximately on the same location as the measurement area in Zone 1 to the forest. The profiles of variables at point 'ext 3' generally have the same appearance as the simulation results in Sect. 5.3 and 5.4 especially for the equivalent leaf area density values.

Conclusions

Experimental and numerical studies of mean flow and turbulence characteristics were performed for the near-field pollutant dispersion campaign carried out at the SIRTA site. Prior to the dispersion modelling, data processing and analysis were performed for the dynamic measurements from ultrasonic anemometers at heights of 3 m, 10 m, and 30 m. First, two years of continuous measurements were used to investigate the effect of the forest on the mean flow for different wind direction sectors. Second, analysis of measurements of a specific IOP (60 min) was conducted to study the turbulence characteristics under stable stratified condition during a pollutant dispersion test. Then, numerical simulations were performed with Code_Saturne applying standard k -ε turbulence model and a canopy model for neutral and stable conditions. Comparisons were made between simulation profiles and measurements at different levels.

Global analysis of continuous measurements over two years highlights the impact of terrain heterogeneity on turbulence in Zone 1 of the site. This impact is dominated by the effect of the northern forest canopy (15 m high), which rotates upper level winds (at 30 m) from a large northerly sector into lower level easterly and westerly winds (at 10 and 3 m). This effect is found mostly for the measurements below the forest height. The directional shear is stronger with decreasing height and decreasing distance from the forest edge.

Numerical simulations are able to reproduce the characteristics of the mean flow at the measurement site, especially the impact of the forest for different wind directions, under both neutral and stable atmospheric conditions. Simulation results illustrate also that the forest changes wind direction and slows down the wind below its canopy height and generates higher TKE above it. For the IOP on 5 June 2013 (IOP-7) which is characterized by a north-easterly wind at 30 m, simulations show a divergent flow in the instrumented area and reproduce the differences in wind speed and direction and TKE between the northern and southern anemometers, as observed in the measurements. They also show that this flow divergence is due partly to the pressure gradient on the forest edge, which leads to the wind becoming tangent to the edge of the forest. Overall, they are in satisfactory agreement with the measurements. They underestimate wind speed and the wind rotation (15 • in the worst case studied) for both neutral and stable conditions, which might be due to the uncertainty of about 5 m in the position of the forest edge. They also tend to underestimate TKE for the stable IOP case, which could be due to an overestimation of the stability in the inlet profiles at upper levels.

Sensitivity studies were performed for the impact of leaf area density α on the simulation results. They showed that the typical features of canopy flow become more pronounced as canopy density increases, as previously observed in many other simulations. With higher values of leaf area density, one can have a larger wind rotation below the forest height and a stronger vertical wind around the height of the forest edge. Comparing with other sensitivity studies conducted for different stability conditions and roughness lengths, it is found that the impact of the forest is the most important factor in accurately reproducing the flows at the site. In the future, a more accurate geometry of the forest would be required to improve further the simulations quality.

Then, simulations have been made with the inlet profiles generated directly from the interpolation of the measurements. The underestimation of wind speed and TKE was reduced, but the wind rotation was still underestimated. We hope to improve those results in the future with data assimilation.

Finally, an additional sensitivity study has been made for the impact of the forest. The 2 years of data has been rearranged for periods in stable and neutral conditions, and for periods in months of summer and winter. Simulations have been performed with a simplified forest with various values of leaf area density. Measurements and simulation results show that the impact of the forest becomes more important with increasing the canopy density. The wind rotation is more pronounced with greater value of α, and has little difference between neutral and stable conditions. Moreover, by replacing the simplified forest by a solid obstacle of the identical dimension, we found that the flow across an obstacle seems like a flow across an extremely dense forest with α → ∞. The impact of the obstacle on the wind directional shear, the horizontal mean wind speed, the vertical velocity and the pressure field are observed to be even more pronounced than the forest case with α = 10 m -1 . The only difference with the forest cases is the absence of the strong turbulent zone above the obstacle top. In the obstacle case, this strong turbulent zone is inhibited by the large values of dissipation rate.

Chapter 6

Concentration data analysis and numerical simulation

The second part of the experimental and numerical study of the SIRTA experiment is the pollutants dispersion. The data concerns the measurements of propylene gas concentration during the IOPs. In this chapter, concentration data processing and analysis for the tracer test will be presented first. The results mainly concern the IOP-7 (on 5 June 2013, from 18:48 to 20:17 see table 2.1). The IOP-11 (on 11 March 2015, from 18:41 to 21:03) with an additional PID at height of 10 m over the PID-3 will be also discussed. The wind and concentration data analysis of the IOP-11 is introduced briefly in the Appendix B. Next, simulation results will be compared with the measurements especially for gas concentration and its fluctuations for these two IOPs. Impact of parameters in the concentration dissipation model will be also discussed.

Data processing

The raw data retrieved from PIDs is treated using calibration coefficients in order to calculate the concentration of propylene gas in ppm by volume. It is recorded every 0.02 s, that is at 50 Hz. During the IOP-7, the tracer test lasts about 76 min, from 18:48 to 20:04. The concentration data series include approximately 228 000 values. Like wind data processing, we first simply plot the "Quick Look" of concentration data. The whole "Quick Look" for IOP-7 are included in Appendix A.3, here we show only 10-min plot to explain the data processing.

From Fig. 6.1 (10-min plot) and A.11 (the whole IOP-7), we observe that the PID-1, 2, 3 and 5 can detect the tracer gas concentration correctly because they all show significant peaks from about 7th minute, while PID-background detects concentration very close to zero all the time. Also, we can find out that the PID-4 does not work correctly in the IOP-7 since it records very few peaks which are small and hard to see. These raw data contains invalid negative values (for example PID-3 Figure 6.1 -'Quick Look' at raw concentration data for the POI-7 from 40th to 50th minutes; from top to bottom: PID-1, -2, -3, -4, -5 and -background; see the whole IOP in Fig. A.11. in Fig. 6.1), and is not ready for analysis. First, we eliminate these invalid data by linear interpolation between neighbours that have non-negative values. Next, in order to remove sensor drift and background concentration from measurement of each PID, two methods have been tried: 1st method : baseline

The baseline method is inspired by the data processing described in [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF] and [Mylne, 1992]. They proposed that a representation of the zero baseline can be derived by dividing the series into small sections and fitting a complex curve through the near-minimum points in each section. These near-minimum points in each section are taken to be the values with approximately 90% of the data in the section lying above it. In practice, the final fitted curve is not very different from a straight line. Thus, our baseline is deduced from the average of the 200 smallest values every 5 min (15 000 values). The baseline obtained from each PID is applied to remove sensor drift and background concentration from its own signal. The baseline removal results are quite satisfactory, since all the small values out of the concentration peaks are brought to near zero. This method introduces inevitable negative values, about 100 negative values over 15 000 values every 5 minutes, but these negative values are very close to zero (between 0 and -0.1 ppm) and have little impact on the data analysis. The advantage of this method is that it allows, for each sensor, to remove both its own drift and the background concentration at the same time. However, the measurements of the background concentration seems not useful here.

2nd method : filtered background concentration

We seek another method which can use the measurements of background concentration. Firstly, the instantaneous background concentration measured by the PID-background is fitted by a 4th order polynomial. Then, we remove this filtered background concentration from the measurements of all the other PIDs. Unfortunately, the results of this method is less satisfactory than the 1st one. It is found that the sensor drifts of different PID are not identical and can not be eliminated simply by the removal of filtered background concentration given by the PID-background. At the same time, we introduce many negative values whose number varies for different PID, making the sensor drift removal of each PID more difficult.

Finally, we choose the baseline method and we obtain the concentration data ready for statistical analysis, plotted in Fig. 6.3 (10-min plot) and A.13 (the whole IOP-7). We observe that PIDs 1, 2, and 3 detected most of the concentration peaks, which is consistent with the slightly south-easterly wind measured by the anemometers at 3 m during the IOP-7. In order to have more accurate comparison with the simulation results, concentration data analysis are made for measurements within the selected sub-period of 60 min (from 19:08 to 20:08, see Sect. 3.1.2) during which meteorological conditions are almost stationary. Thus, data analysis are made for the concentration measurements from 19:08 to 20:04 during the IOP-7.

6.2 Data analysis [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF] have studied a series of tracer experiments in nearneutral and slightly convective condition with measurements at range of 50 m to 1000 m from the source. [Mylne, 1992] has analysed a set of tracer experiments near the surface at short range from the source and in a stably stratified nocturnal boundary layer. He has also compared the results with those obtained in [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF] in near-neutral stability. The concentration data analysis of the SIRTA dispersion experiment is mainly based on the methods applied in these two papers. Their discussions and conclusions also help us to better understand our own concentration data.

Development of dispersing plume

Before showing the results of data analysis, it is interesting to discuss a conceptual model of plume dispersion proposed by [Mylne andMason, 1991, Mylne, 1992]. Figure 6.4 shows a typical form of an instantaneous dispersing plume. The model states that the dispersing plume is under the influence of turbulent eddies on a wide range of scales. Thus, it suggests three different regimes according to the relative size of the eddies length scale λ and the instantaneous plume width σ i as follows: Figure 6.4 -Illustration of the conceptual model of plume dispersion [Mylne, 1992]. It shows the three different effects of eddies of length scale λ according to the relative size of the crosswind length scale of the instantaneous plume σ i .

I. σ i

λ: Close to the source, where turbulent eddies are large compared to the instantaneous transverse plume, these eddies act to move the whole plume around thus cause the plume to meander. This meandering induces the intermittency effect in concentration measurements by detecting periods of non-zero and zero concentration whose durations are determined by the frequency of the meander. [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF]] also showed that, for the plume dominated by meandering, the intermittent concentrations are better fitted by the exponential PDF (probability density function). II. σ i ∼ λ: As the plume travels further and grows into a scale of order of λ, the eddies of scale λ act to break up the plume and entrain pocket of clean air into it. This process generates and maintains internal concentration variance within the instantaneous plume. III. σ i λ: At much longer distances, when the plumes has grown much wider than λ, eddies of scale λ act to mix the pollutant within the plume, smoothing out the concentration fluctuations and transferring the variance to smaller scales, which is the process characterised by the inertial subrange in the concentration spectrum. At this regime, the concentrations are less intermittent, and the instantaneous plume becomes less different from the mean plume.

The concentration data is better fitted by the Gaussian PDF.

It should be noted that, during pollutants dispersion, these three processes act to the plume simultaneously with eddies of all scales. [Mylne, 1992] also suggested that, in stable conditions, the second process is reduced and the plume usually moves into the third regime closer to the source than in near-neutral or convective conditions.

In our study, the characterization of the dispersing plume from data analysis is based on this conceptual model.

Concentration time series

Under stable stratification, vertical motions are suppressed by the buoyancy forces , and the TKE production is reduced. Consequently, the rate of turbulent diffusion is small and the instantaneous plume grows relatively slowly [Mylne, 1992]. However, the relatively larger scale horizontal motions are not suppressed by the vertical stability forces, and may persist for long periods since the turbulent friction is reduced. Thus, the meandering of the plume dispersion might be more easily observed than under neutral or convective conditions. This meandering effect can be simply observed from instantaneous concentration measurements. For example, concentration data of experiments in [Mylne, 1992] shows continuous fluctuating non-zero concentrations for a substantial period of time (usually several minutes) and nothing during another long period. The concentration data from the SIRTA experiment shows also similar characteristic. From 10 min measurements in the IOP-7 (Fig. 6.3), we can see that each PID detects alternatively fluctuating non-zero concentrations during several minutes and near-zero concentrations during another several minutes. Except that PID-2 detects more frequently large concentration peaks. We suppose that PID-2 is located near the center line of the mean plume during the IOP-7, and detects more often the instantaneous meandering plume.

Since PID-4 did not work during the IOP-7, here I use the concentration time series of the IOP-11 (See introduction in the Appendix B) to illustrate more clearly the correlation between the plume meandering and the change of velocity v. Figure 6.5 plots the 10-min concentration data from 20th to 30th minute of the IOP-11. We observe that continuous fluctuating non-zero concentration is mostly detected by PID-5 for the first 5 minutes, then great peaks begin to appear in PID-2 and -3 for the second 5-minute period. Fig. 6.7 shows the velocity v filtered by central moving average over 1 min for the five anemometers collocated with PIDs. From the first part of Fig. 6.7, we find that the velocity v changes from -0.5 m.s -1 to near zero during this 10-minute period. Next, Figure 6.6 plots the 10 min concentration data from 70th to 80th minute of the IOP-11. We observe that, from 73th to 75th minute, PID-1 detects some concentration peaks while PID-3 and -4 detect almost nothing. If we look at the variation of the velocity v (the second part of Fig. 6.7), we find that v is positive during this 10-minute period and has a sudden increase around 73th minute and then a decrease around 75th minute. Thus, from measurement of the SIRTA experiment, we also find the evidence of the plume meandering caused by large-scale horizontal motions as described in [Mylne, 1992]. After the observation of time series, we plot concentration histogram. Figure 6.8 shows the concentration histograms for the IOP-7 within the 60 min sub-period and for the PID-1, -2, -3, -5 and -background. Keeping in mind the conceptual model of plume dispersion, if measurements are made in near-source region (regime I and II in Fig. 6.4), histograms should be close to an exponential distribution [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF]. Because sensors detect intermittent concentrations, with periods of near-zero concentration and periods of concentration peaks. Those peaks with high concentration values are detected in small numbers. On the other hand, if measurements are made at long distance from the source (regime III in Fig. 6.4), histogram should be close to a Gaussian distribution since concentration within the plume varies around its average value.

Concentration histogram

In the IOP-7 (Fig. 6.8), exponential distributions were obtained for all the PIDs except that the PID-background's histogram is much narrower than others. Most of the data detected by the PID-background is near-zero which is logical for the background concentration. PID-1 and PID-2 detect much more high concentrations than PID-3 and PID-5 which is consistent with the fact that they detect more concentration peaks during the IOP. These histograms confirm that in the SIRTA dispersion experiment, 50 m from the source is a near-source region where pollutant dispersion is influenced by the plume meandering effect.

Statistical analysis

Statistical values such as intermittency factor I (probability that the concentration is non-zero), mean concentration C, conditional mean C p , standard deviation of concentration σ C , fluctuation intensity σ C /C, conditional intensity (σ C /C) p are presented in Table 6.1, where the subscript p is used to indicate that the statistical parameter is calculated from non-zero concentrations only. These statistical values are calculated for the data of the selected 60-min sub-period in the IOP-7. In practice, non-zero concentrations are defined as those which exceed some threshold value C T , in order to ensure that real tracer gas concentration can be distinguished from background and instrumental noise. In the IOP-7, a threshold value C T = 2 ppm has been chosen to define the non-zero concentrations. With this specific value, we can have the intermittency factor I almost zero for the PID-background (I = 0.002) and I between 0.4 and 0.6 for the PIDs near the centreline of the plume, which are the typical values suggested in [Mylne, 1992]. Thus, those conditional statistical parameters are obtained based on concentrations larger than the threshold value C T .

From Table 6.1, we observe that intermittency factors and mean concentrations vary in the same way and have higher values for PID-1 and -2 which is coherent with the south-easterly wind measured at 3 m height. However, the conditional mean C p is near-constant and much higher than C, which indicates that the plume behaviour is dominated by meandering caused by turbulent eddies larger than the instantaneous plume. Similar results are reported in [Mylne, 1992] for measurement under stable conditions. It suggests that the mean concentration in stable condition seems to be determined primarily by the meandering of the plume, which has the effect of greatly reducing the concentration at a fix point averaged over half an hour or more.

The standard deviation of concentration σ C has higher values for PIDs who detect more concentration peaks. With larger intermittency factor I, the fluctuation intensity σ C /C becomes smaller. However, the conditional intensity (σ C /C) p is approximately constant between PID-1, -2 and -3, except PID-5 who detects few concentration peaks and is supposed to be usually out of or on the edge of the plume. These results show again the evidence of plume meandering, and indicate that the turbulent fluctuations within the plume are independent of its meandering position as stated in [Mylne, 1992].

Similar observations for concentration statistical analysis can be made for other IOPs. For example, the IOP-11 (See Appendix B Table B.2), which is performed under a more stable stratification than the IOP-7, shows much greater difference in values between non-conditional and conditional statistical parameters (C and C p , σ C /C and (σ C /C) p ).

Integral length scale

We deduce the integral time scales and length scales from the autocorrelation of the concentration time series in the same manner as for the velocity components (Sect. 3.2.2). As shown in Fig. 6.9, time scale of concentration T C is obtained with the approximation that it is equal to the time lag τ e where the autocorrelation coefficient goes below 1/e for the first time. The length scale of the plume structure is then calculated with relationship L C = a mean T C , where a mean is the longitudinal mean wind speed measured by the collocated anemometer. The values of L C deduced from the concentration measurements within the 60 min sub-period of the IOP-7 are reported in Table 6.1. Figure 6.9 -Autocorrelation of the concentration data from the measurement of PID-1 in the IOP-7; Integral time scale is defined approximately by the time lag τ e where the autocorrelation coefficient goes below 1/e for the first time.

Comparing with the turbulent length scales L aa , L bb and L ww during the same period (See Table 3.1), we find that the length scales of the plume L C are larger than the vertical turbulent length scale L ww but smaller than the horizontal ones L aa and L bb . Moreover, L C seems closer to the magnitude of transverse length scale L bb . This result implies that, at 50 m from the source in a stable condition as the IOP-7, turbulence eddies seem to have larger horizontal scales than the plume and they cause the plume to meander horizontally. The crosswind velocity might have more effect than other components on the plume meandering. [Mylne, 1992] has also compared the concentration autocorrelation function R(τ ) between a stable case and a near-neutral case. He comments that the time scales in stable conditions are about an order of magnitude larger than those in neutral conditions. R(τ ) in stable conditions can stay close to unity for time lags up to about 1 s, while R(τ ) in neutral conditions falls much more rapidly. Our functions R(τ ) deduced from different PIDs during the IOP-7 seem to fall less rapidly than the neutral example shown in [Mylne, 1992], and they usually stay larger than 0.7 before 1 s which is relatively high for an autocorrelation coefficient. However, since the plume scale depends on several factors, such as the distance from the source, the diffusion rate due to small-scale turbulence and the action of eddies of a similar size which act to break up the plume, we can not draw too detailed conclusion by comparing our plume scales with those in the literature.

Concentration power spectra

Normalized concentration spectra f.S(f )/σ 2 C are calculated for PID-1, -2, -3, -5 and -background and plotted in Fig. 6.10. Spectra generally follow well the -2/3 slope on the high frequency range showing the existence of the inertial subrange. Because of the intermittency effect, it is much more difficult to observe a clear inertial subrange in the measured concentration spectra in stable condition. [Mylne, 1992] suggested that in order to observe a clear inertial subrange, it is necessary to select time series with as large intermittency factor as possible. That is why the spectra of PID-5 and -background are much more fluctuating than that of PID-1, -2 and -3, since they are out of or on the edge of the plume and have an intermittency factor much smaller than others. However, we find that the spectrum of PID-background still shows somehow an inertial subrange, which implies that there might be other sources around during the IOP. We suppose that the concentration detected might come from the road on the other side of the forest to the north of the SIRTA site.

As stated in [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF], in the development of a plume, the three kind of processes shown in Fig. 6.4 act simultaneously for eddies of all scales, so that as the plume grows with increasing distance from the source, large eddies are incorporated into the processes generating the concentration inertial subrange. They also pointed out that, for measurements at a point close to the source where the plume is dominated by meandering motions, the inertial subrange only occurs at very 148 Chapter 6. Concentration data analysis and numerical simulation high frequencies, because only small eddies can contribute to the mixing within the plume. Comparing the concentration spectra of the PID-1, -2, -3 and -5 in the IOP-7 with spectra in [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF] measured at much larger distance (740 m) from the source, the inertial subrange in the IOP-7 can be found only for frequency higher than 0.1 Hz, while the inertial subrange shown in [START_REF] Mylne | Concentration fluctuation measurements in a dispersing plume at a range of up to 1000 m[END_REF] (Fig. 6.10(f)) can be identified since 0.01 Hz. This fact shows again the evidence that our measurements are strongly influenced by the plume meandering process. Moreover, the spectra of PID-background seems to follow the -2/3 slope from a lower frequency than others, which implies that it might measure the concentration far from a source as the process III described previously in Fig. 6.4. C . The inlet conditions are set to be zero for both C and σ 2 C (no flux at the ground). Their initial conditions are also defined as zero for all the modelling domain. At the source cell, a volume source term of concentration C is imposed, using the real flow rate during the corresponding IOP. However, the value of concentration fluctuation σ 2 C at the source is unknown during the experiment. According to simulation results in [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF] who used the same tool to simulate the dispersion of the MUST campaign, we have tried σ C =0%C, σ C =1%C and σ C =10%C at the source for the SIRTA dispersion.

Simulation results and discussions

It should be noted that simulation in Code_Saturne gives concentration in mass fraction such that:

C CS = M C 3 H 6 M C 3 H 6 + M air ≈ M C 3 H 6 M air (6.1)
while concentration measurements are in ppm by volume (noted as ppmv). Before analysing the results, we have to convert concentration in mass fraction into ppm by volume:

C CS • 1/ρ C 3 H 6 1/ρ air • 10 6 = C in ppmv (6.2)
with density of propylene ρ C 3 H 6 =1.81 kg.m 3 , and air density ρ air = 1.23 kg.m 3 .

IOP-7 with wind rotation effect

The gas release during the IOP-7 has a flow rate D=160 l.min -1 . The dispersion is simulated with two different leaf area densities α = 0.9 m -1 and 2 m -1 . Figure 6.11 shows simulation results for the concentration field in horizontal and vertical cross-sections. The horizontal cross-sections are at height of 3 m, with horizontal wind vectors and grey circles representing the devices location in the instrumented area. Vertical cross-sections are at the location of the source along a E-W line with grey circles representing the source (3 m), PID-3 (3 m) and the additional PID-6 at 10 m over PID-3. Although PID-6 has been only added at the most recent IOP-11, here I still mark it in the simulation results in order to have a reference for height. We can see in Fig. 6.11 that the plume shapes are different with different forest densities. Since the case with α = 2 m -1 induces a stronger wind rotation, we observe higher mean concentration in the north close to the forest. Figure 6.12 compares the simulated mean concentrations with measurements for the IOP-7. The peaks of mean concentration C in simulations appear to be more south than the measurements, which can be explained by the fact that the wind directional shear due to wind channelling effect from the forest has been underestimated by simulations (See simulation results in Sect. 5.4). However, the peak values of C between simulations and measurements are in good agreement. With increasing leaf area density α in the simulation, the wind turns more to the north thus brings concentration peak also more to the north, and the agreement is better with the measurements. Figure 6.13 shows the horizontal and vertical cross-sections of concentration fluctuations σ C for simulations with the dissipation model constant R f =0.5 (eq. 4.16) and the value of σ C =1%C at the source for α = 0.9 and 2 m -1 . We observe that the concentration fluctuations have larger spread than the mean concentration both horizontally and vertically. For the vertical cross-sections, the concentration fluctuations σ C reach their maximum levels at a greater height than the mean concentrations, which has also been observed in [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF]. ues of parameter R f in the dissipation model (eq. 4.16) for leaf area density α= 0.9 m -1 (Fig. 6.14(a)) and α= 2 m -1 (Fig. 6.14(b)). Like the mean concentration, the leaf area density affects the degree of wind rotation, thus changes the position of maximum values for σ C . Peaks of σ C all turn more to the north with α= 2 m -1 than with α= 0.9 m -1 . [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF] has studied the influence of the emission rate at the source. They have performed simulations with the values of σ C (source)/C(source)= 0, 1%, 10% and 100% and found that the difference between 0 and 1% is nearly non-existent and is small between 1% and 10%, and 100% is unrealistic and largely overestimated σ C . We have also performed simulations with the values of σ C (source)/C(source)= 0, 1% and 10%. We find that σ C =0%C gives almost the same concentration fluctuations as σ C =1%C with difference no more than 1 ppm (not shown in figures). However, σ C =10%C gives concentration fluctuations much larger (more than factor 2) than σ C =1%C. Comparing with the measurements, we find that σ C (source)=1%C(source) gives closer results, which implies that the concentration fluctuations at the source is very small in the SIRTA experiment.

The influence of parameter R f in the dissipation model has been also studied previously in [Hsieh et al., 2007, Milliez andCarissimo, 2008]. [ [START_REF] Milliez | Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area[END_REF] have shown the dependence of the results on the values of R f . Similarly in our simulations, we obtain smaller σ C with smaller R f . Values of R f proposed in the previous papers varied from 0.9 to 0.5. For the SIRTA experiment, simulations with all these values tend to overestimate the concentration fluctuations when compared with the measurements. As stated in [Hsieh et al., 2007, Milliez andCarissimo, 2008] and in Sect. 4.2.2, R f is a model constant that represents the ratio of the concentration dissipation time scale to the dissipation time scale of the turbulence (equation 4.15). Such results implies that during the SIRTA field measurement, the concentration dissipation time scale is much smaller than that of the turbulence, which is consistent with the situation of near-source region and in stable condition (regime I in Fig. 6.4). Also, the overestimation of σ C might be due to the fact that the concentration measurements at SIRTA are at a smaller range (50 m), thus with stronger intermittancy effect, than experiments reported in the literature. To have a better result, we could reduce more R f , but the value of 0.5 is already very small according to the literature. Actually, the simple dissipation model that we applied assume that the dissipation time scale is proportional to the turbulence time scale (T C ∝ k/ε). This model is strictly valid for homogeneous turbulent mixing in a homogeneous scalar field, as indicated in [START_REF] Hsieh | Numerical modeling of passive scalar dispersion in an urban canopy layer[END_REF]. That is to say, this model is valid for the case where the size of the instantaneous plume is much larger than the length scale of the background turbulence, which is not the case for the SIRTA experiment. To overcome the limit of this simple model, one should apply another more complex model with a dissipation length scale that varies with the local characteristics of the plume (See Sect. 4.2.2 eq. 4.14).

For now, the existing simulations show that, with the value of σ C (source)/C(source)= 1% and the parameter R f =0.5, we obtain the closest results to the measurements for the concentration fluctuations.

IOP-11 with additional measurement

A new PID has been added in this IOP at height of 10 m above PID-3, which allows us to have an idea about concentration dispersion at a higher level and compare with vertical profiles obtained in the simulations. As we can see from Table 6.2, the IOP-11 has almost an easterly wind at height of 3 m and a slightly south-easterly wind at higher levels, thus dispersion should be less affected by the wind rotation due to the forest to the north. The gas release during the IOP-11 has a flow rate 156 Chapter 6. Concentration data analysis and numerical simulation D=200 l.min -1 .

- Figure 6.15 shows the horizontal and vertical cross sections for the mean concentration (Fig. 6.15(a)) and the concentration fluctuations (Fig. 6.15(b)) simulated with the leaf area density α = 0.5 m -1 , R f =0.5 and σ C (source)=1%C(source). All the cross-sections are at the same locations as previously presented for the IOP-7. The mean concentration is also at the same scale as the IOP-7. We can observe that the dispersion of the plume is coherent with the wind direction. Since there is no wind rotation effect (and thus no flow divergence) and a higher flow rate, the plume appears to be narrower and the concentrations are higher at the "sonic arc at 50 m" than that shown in Fig. 6.11. The concentration fluctuations also have higher values than that in the IOP-7. Similarly, σ C has larger spread than the mean concentration and it reaches its maximum at a greater height than the mean concentration. Moreover, the IOP-11 is under a more stable stratification than the IOP-7, which might be another reason for the narrower concentration plume and higher values of C and σ C detected in the IOP-11 than in the IOP-7.

Simulations for the IOP-11 have been performed with 2 different leaf area densities α = 0.5 and 0.9 m -1 . Figure 6.16 shows the comparison of mean concentration with measurements at height of 3 m (Fig. 6.16(a)) and vertical profiles at location of ). From Fig. 6.16(a) we observe that the peak values of C are in good agreement between simulations and measurements. Without the wind rotation effect, the positions of maximum of C are also very close. The reason for the simulated concentration peaks slightly more to the north than in the measurements is that the meteorological boundary conditions are calculated with measurements at higher levels where the wind are detected to be slightly south-east. From Fig. 6.16(b), we can see that the simulations slightly underestimate the concentration at 3 m (PID-3) due to the shift of the concentration peak position, and slightly overestimate the concentration at 10 m (PID-6). Moreover, we find that there is little difference between simulations with α = 0.5 and 0.9 m -1 for the dispersion as well as for the wind and turbulence, which confirms that in this case the forest does not have a significant effect due to the wind direction. Figure 6.17 shows the comparison of concentration fluctuations with measurements at 3 m (Fig. 6.17(a)) and vertical profiles at location of PID-3 and PID-6 (Fig. 6.17(b)). The same sensitivity study is made for the value of the emission rate at the source and the constant R f in the dissipation model. The same conclusions are reached with the results of the IOP-7: it is with the value of σ C (source)/C(source)= 1% and the parameter R f =0.5 that simulation gives best results comparing with the measurements. However, at 10 m, all the simulations seem to largely overestimate the concentration fluctuation compared with the measurements of PID-6 which is only 4 ppm. This might be due to the overestimation of the TKE at height of 10 m in the simulations of the IOP-11.

Conclusions

Previous studies show that plume dispersion in a near-source region appears to be meander-dominated, and the meandering occurs more easily in stable conditions, due to large scale horizontal motions and small turbulent diffusion rate. Concentration data analysis shows that the measurements of the SIRTA dispersion experiment in stable condition at 50 m from the source seem to be under the influence of plume meandering. Alternative periods of zero and non-zero concentrations, each persisting for several minutes, are observed in the concentration time series. Concentration histograms are closer to exponential distributions which means that high concentration values are detected in small numbers. Great differences are found between non-conditional and conditional statistical parameters (C and C p , σ C /C and (σ C /C) p ), and conditional statistical values are almost constants between the PIDs across the plume, which indicate that at 50 m from the source, the plume behaviour is dominated by the meandering caused by turbulent eddies larger than the scale of instantaneous plume. The plume meandering effect has largely reduced the mean concentration and increased the fluctuation intensity when they are calculated for a long averaging period. Length scales of the plume are calculated to be smaller than the scales of turbulent horizontal structures. Concentration fluctuation spectra exhibit inertial subrange at very high frequency, which corresponds to the eddies of smaller scales than that of the plume.

Simulation results show that the mean concentrations and the concentration fluctuations are sensitive to the simulated flow features such as wind direction and turbulent kinetic energy which are difficult to reproduce exactly as in the real situations. The mean concentrations show a good agreement with measurements in values for all the IOPs studied, except that the position of the concentration peak depends on the accuracy of simulated wind rotation below the forest height. The concentration fluctuations obtained from the simulations seem to be affected significantly by the fluctuation emission rate at the source and the values of the constant R f in the dissipation model. From the existing simulations, we find that the value of σ C (source)/C(source)= 1% and the parameter R f =0.5 give the closest results to the measurements for the concentration fluctuations. However, the current simple dissipation model, which is valid for homogeneous turbulent mixing, might not be adequate for the SIRTA experiment, since we have shown from the data analysis that the measurements are under the meandering effect and the plume size is smaller than the large-scale turbulent structures. Another more complex dissipation model taking into account the instantaneous plume characteristics (Sect. 4.2.2 eq. 4.14) might be able to give better results. Furthermore, simulations with RANS approach can only give the average and fluctuation values. If we want to illustrate the meandering effect during the pollutants dispersion, unsteady models such as large-eddy simulations would be more appropriate.

Conclusions and Perspectives

Conclusions

The goal of this thesis is to study atmospheric flow and dispersion over a complex site with a focus on stable stratification. This study is based on both numerical simulations and analysis of measurements obtained at the SIRTA site. The first part is the experimental study with data analysis of wind and concentration measurements, with the aim of characterizing the wind, the turbulence and the pollutants dispersion in the atmospheric boundary layer under the impact of stability and field heterogeneity. The second part is the numerical study with modelling and simulations of the field experiment taking into account the real situations during the IOPs, which allows to validate and improve the performance of the atmospheric module of the CFD code, Code_Saturne, developed at EDF and CEREA. Data analysis have been made for two years of continuous measurements for wind, and during IOPs for wind and concentration. Turbulence at the SIRTA site has been shown to be strongly anisotropic with different orders of magnitude between the three velocity components variances, and with different spectra forms between vertical and horizontal velocity components. The turbulent structures advection speed obtained by velocity cross-correlation between sensors at 3 m above the ground is found to be higher than the measured average wind speed at this height. These results can be linked to the fact that turbulence in the stable surface layer is under the impact of a strong wind shear and the blocking effect due to the presence of the ground. These characteristics are found to be consistent with the results of other field or theoretical studies for the turbulent flows in the atmospheric surface-layer very close to the ground surface, such as [START_REF] Powell | An Investigation of the Application of Taylor's Hypothesis to Atmospheric Boundary Layer Turbulence[END_REF], Horst et al., 2004, Drobinski et al., 2004].

By dividing 2 years of data into different stability classes, we have found that the velocity fluctuations normalised by the mean wind speed decrease continuously from very unstable to very stable stratification. Integral length scales have higher values in a less stable surface layer and are increasing with altitude, which shows that turbulence in the surface layer exhibits larger structures with increasing height and instability. The different order of magnitude between integral length scales of different velocity components shown as L aa L ww can be found for all the stability classes due to the blocking effect of the ground. This difference indicates that flows in the SIRTA Zone 1 are affected by a mixture of small scale and large scale structures, which are dominant for the vertical length scale and the longitudinal length scale respectively. Analysis of the velocity cross-correlation shows that the ratios of the eddy advection velocity to the mean wind speed r = U adv /U have slightly higher values under more stable conditions and for measurements closer to the ground. This is coherent with the fact that the wind shear is more important near the ground and under a more stable surface layer. By studying different wind direction sectors, we have illustrated the impact of the forest to the north on the flows of the measurements area. The forest changes wind velocity and direction for a large northerly sector. It also induces strong velocity fluctuations and creates smaller eddies above its top. These findings are consistent with the works on the studies of flows under different stability conditions [Kaimal, 1973, Barthlott et al., 2007, Fesquet et al., 2009a[START_REF] Fesquet | Impact of terrain heterogeneity on coherent structure properties: numerical approach[END_REF], Pena et al., 2010], and on the observations of flows across vegetation zones [START_REF] Irvine | The evolution of turbulence across a forest edge[END_REF], Thomas and Foken, 2007, Chahine et al., 2014].

Concentration data analysis has shown the evidence of the meandering effect on the dispersion.

The observations are consistent with the work of [Mylne andMason, 1991, Mylne, 1992] for the studies of pollutants dispersion near the source and under stable conditions. Evidence has been found from the concentration time series, histograms, length scales and power spectra, which shows that the dispersing plume is under the impact of meandering: Alternative periods of zero and non-zero concentrations, each persisting for several minutes, are observed in the concentration time series; Concentration histograms are close to exponential distributions; Great differences are found between non-conditional and conditional statistical parameters; Concentration fluctuation spectra exhibit inertial subrange at very high frequency. Thus, 50 m from the source is a near-source region where the plume scale is smaller than the large scale turbulent structures, which is the main cause of the plume meandering effect. Also, the great difference between conditional and non conditional statistical values has shown that the meandering effect is able to largely reduce the mean concentration and increase the fluctuation intensity when the analysis are made for a long averaging period. It is possible that, during the long periods with relatively low concentrations, the meandering plume has induced high concentrations at some places for short periods which could be already harmful for human health. Thus, it is necessary to study the dispersion under stable condition for the need of risk assessment.

The numerical study of flows and dispersions have been carried out by using the CFD simulations in RANS approach with the standard k -ε turbulent model adapted for the atmospheric flows and a canopy model for the forest. The pollutant dispersion has been modelled by transport equations for the mean concentrations and the concentration fluctuations.

The simulations have been shown to reproduce correctly the characteristics of mean flows on the measurement site, especially the impact of the forest for different wind directions, in both neutral and stable conditions. However, the exact wind rotation effect and TKE in the real situations are difficult to obtain in the simulations. This might be partly due to the inlet profiles which are not measured but are built with the theoretical profiles and the measurements inside the simulation domain. Also, the stability condition varies between 3 m, 10 m and 30 m in the reality, while during the simulations we supposed that the stability condition is uniform vertically. Another reason might be that there was a position uncertainty of at least 5 m when we extracted the forest geometry from the landuse file and the satellite map. Since the wind rotation and the TKE detected at the sensors locations are very sensitive to the shape of the forest and the distance to it, a difference of 5 m could be able to change the results. Still, the performance of the canopy model has been found to correspond to the simulation works of others for the canopy flows [START_REF] Shaw | Large-eddy simulation of turbulent flow above and within a forest[END_REF][START_REF] Dupont | Influence of foliar density profile on canopy flow: a large-eddy simulation study[END_REF], Dupont and Brunet, 2008a, Dalpé and Masson, 2009, Zaidi et al., 2013], such as the wind directional vertical shear, the pressure depression within the canopy and the strong TKE generated above. The sensitivity study on the impact of the forest show that the typical features of canopy flows become more pronounced with increasing canopy density, which is also the conclusion found in the work of [START_REF] Dupont | Influence of foliar density profile on canopy flow: a large-eddy simulation study[END_REF].

The accuracy of simulated pollutants dispersion depends on the accuracy of simulated flows. The maximum values of the mean concentrations have been correctly reproduced by the simulations for all the IOPs studied. However, the maximum concentration positions depend on the simulated wind rotation effect. The concentration fluctuations have been shown to be affected by the emission rate at the source and the parameter of the simple model of the dissipation term.

Overall, the performance of the atmospheric module of Code_Saturne in RANS mode is satisfactory after comparing the simulation results with the measurements.

Perspectives

This thesis is the first work on the SIRTA dispersion experimental program. New field measurements have been obtained recently that are not included in my work. In the future, new PIDs could be added in another circular arc at 100 m from the source at height of 3 m and also at higher levels. These measurements would allow us to have a more complete vision of the plume dispersion, and to identify the change of plume characteristics as function of distance to the source.

With the continuous measurements of wind since 2012 and the short-periods of concentration measurements during around a dozen of IOPs, the SIRTA dispersion experimental program contains already a very rich dataset. A lot of additional interesting data analysis could be carried out based on my work. For example, the behaviour of the velocity spectra can be studied in different stability conditions and in different wind sectors. The turbulent dissipation rate ε can be deduced from velocity power spectra of the measurements and then compared with the ε in the simulations. Moreover, correlation studies can be made between wind data and concentration data. In my work, we have only related them visually by the observations of time series. I believe that correlation studies under different stability conditions and wind sectors could be a great step in better understanding the actions of turbulence on pollutant dispersion.

To evaluate the performance of the atmospheric module in Code_Saturne, we have started with the RANS approach with the k -ε turbulent model. We made this choice because this kind of calculation is less time consuming, and it corresponds to what is performed currently for operational applications in the domain of atmospheric dispersion at EDF. However, we can only compare the averaged values with the measurements and only for variables such as the mean wind speed and direction, the TKE, the mean concentration and its fluctuation. In the future, it is planned to use another turbulent model R ij -ε in the RANS approach. This is a second moment closure derived directly from Reynolds-stress transport equations, which is able to take into account the turbulence anisotropy in simulations. With the R ij -ε model, we could compare the velocity fluctuations and other second order moments such as w T and w c between simulations and measurements. Furthermore, to simulate the meandering effect under stable conditions, simulations with the Large-Eddy Simulations would be useful. Simulations in LES could illustrate unsteady turbu- ) 0.17 0.17 0.14 0.12 0.17

Figure 1 -

 1 Figure 1 -Rising smoke in Lochcarron, Scotland, forms a ceiling over the valley due to a temperature inversion. (Photo by http://moonrise.ws/ S/V Moonrise, on January 2006)

  the convective temperature scale with Q 0 = w Θ the kinematic flux of sensible heat, φ m (ζ) and φ h (ζ) are some non-dimensional universal functions (stability profile functions) of ζ = z/L, which imply respec-Chapter 1. Atmospheric boundary layer and pollutant dispersion

  39) with κ = 0.42 and u * 0 = U ref κ/ ln(z ref /z 0 ) where U ref and z ref are velocity and height at a measurement point.

Figure 2

 2 Figure 2.1 -Measurement areas in the SIRTA field. The SIRTA dispersion experiment is carried out in Zone 1 (Z1).

Figure 2

 2 Figure 2.2 -(a) Tracer gas release point and (b) collocated ultrasonic anemometer (front) and photo ionisation detector (behind)

Figure 2 .

 2 Figure 2.3 shows a satellite map of spatial arrangements of various instruments used in SIRTA experimental program in Zone 1.

Figure 2

 2 Figure 2.3 -Satellite map of spatial arrangements of various instruments used in SIRTA experimental program

Figure 2

 2 Figure 2.4 -Sensors locations in Zone 1 with their labels: purple diamond -source; red triangle -ultrasonic anemometer; green dot -photo-ionization detector.

Figure 2 .

 2 Figure 2.4 shows the position of sensors used for the experiment. Red triangles

Figure 2

 2 Figure 2.5 -3D ultrasonic anemometer

  worked; stable condition; an additional PID measuring at 10-m height above PID-3

Figure 3

 3 Figure 3.1 -(a) Conventional meteorological reference (u, v, w) and (b) reference in the mean wind frame (a, b, w)

Figure 3

 3 Figure 3.3 -Illustration of definition of the time lag τ e for the calculation of the approximation integral time scale T e from data of IOP-7

Figure 3

 3 Figure 3.4 -Spatial cross-correlation of anemometer couples (SE, SW), (NW, 20N), (SW, 10S) and (SW, 20S) for longitudinal and transversal wind components as a function of normalized time lag for IOP-7

Figure 3

 3 Figure 3.6 -Mean wind speed as function of measurement height for different stability classes. Values shown represent the median of mean wind speed calculated for the selected 30-min periods over two years

Figure 3 Figure 3

 33 Figure 3.7 -Velocity variance σ 2 a and σ 2 w as function of measurement height for different stability classes. Values shown represent the median of velocity variance calculated for the selected 30-min periods over two years

Figure 3

 3 Figure 3.10 -Integral length scale L aa and L ww as function of measurement height for different stability classes. Values shown represent the median of integral length scales calculated for the selected 30-min periods over two years

Figure 3 .

 3 Figure 3.11 -Ratios between eddy advection velocity and mean wind speed r a and r b as function of measurement height for different stability classes. Values shown represent the median of ratios calculated for the selected 30-min period over two years.

Figure 3 .

 3 Figure 3.12 -Mean wind speed as function of measurement height for different incident wind sectors in neutral condition (class D). Values shown represent the median of mean wind speed calculated for the selected 30-min period over two years.

Figure 3

 3 Figure 3.13 -Velocity variance σ 2 a and σ 2 w as function of measurement height for different incident wind sectors in neutral condition (class D). Values shown represent the median of variances calculated for the selected 30-min periods over two years.

Figure 3 .

 3 Figure 3.14 -Integral length scale L aa and L ww as function of measurement height for different incident wind sectors in neutral condition (class D). Values shown represent the median of integral length scales calculated for the selected 30-min periods over two years.

Figure 5

 5 Figure 5.1 -Mean wind direction frequency distributions for measurements over two-year period (April 2012 -March 2014) at three levels: 3 m (NE and SE), 10 m (10mSW) and 30 m (30mSE). At 3 m, there is a decrease in frequency in the northerly sector and a corresponding increase in the easterly and westerly sectors, which is absent for measurement at 30 m.

Figure 5

 5 Figure 5.2 -Mean wind direction from measurements over two years (April 2012 -March 2014) at three levels: 3 m (NE and SE), 10 m (10mSW) and 30 m (30mSE), selected for the periods for which mean wind direction measured by 30mSE is (a) 40 • -50 • and (b) 315 • -325 • , and for |L M O | > 200 m measured by 10mSW. The dashed lines extend from the lowest to highest values of the data with "diamond" at the median. The box indicates the 25th and 75th percentiles of the data.

Figure 5

 5 Figure 5.3 -Measurement areas in SIRTA field: red rectangular -modelling area; yellow rectangular -instrumented area.

Figure 5

 5 Figure 5.4 -Mesh for the modelling area taking into account the land-use map and the eastward shelters: (a) horizontal cross-section at ground level; (b) zoom on the instrumented area at ground level; (c) vertical cross-section. Grey circles indicate the positions of measurement devices.

Figure 5 . 5 -

 55 Figure 5.5 -Comparison of wind direction deviation dd-dd 30 (top), normalized horizontal wind speed f f /f f 30 (middle) and normalized turbulent kinetic energy T KE/f f 2 30 (bottom) for classical logarithmic analytical profiles (dashed lines), simulations (solid lines) under neutral conditions of wind directions: (a) 45 • and (b) 90 • , and continuous measurements (box plots) over two years (April 2012 -March 2014) for selected wind directions: (a) 40 • -50 • and (b) 85 • -95 • at 30 m with |L M O | > 200 m at 10 m. Box plots represent measurements of anemometers, from bottom to top: NE, NW, SE, SW, 10mSW, 10mSE, and 30mSE (only for T KE/f f 2 30 ). The boxes extend from the lower to upper quartile values of the data with vertical red line at the median. The dashed lines extend from the lowest to the highest values of the data. The subscript '30' means values of variables at 30 m. The horizontal red line indicates the vertical level of 3 m.

5. 3 .

 3 Figure 5.6 -Same as Fig. 5.5 for simulations (solid lines) under neutral conditions of wind directions: (a) 320 • and (b) 270 • , and continuous measurements (box plots) over two years for selected wind directions: (a) 315 • -325 • and (b) 265 • -275 • at 30 m with |L M O | > 200 m at 10 m.

  • -50 • (Fig. 5.5a), 85 • -95 • (Fig. 5.5b), 315 • -325 • (Fig. 5.6a) and 265 • -275 • (Fig. 5.6b) at height of 30 m with |L M O | > 200 m at 10 m. The box plots in the figures represent the distribution of measurements for each direction.

Figure 5

 5 Figure 5.7 -Horizontal cross-section at 3 m of (a) pressure field with horizontal wind vectors for Zone 1, (b) velocity field for Zone 1, (c) TKE field for Zone 1, (d) velocity field with horizontal wind vectors for the instrumented area, and (e) TKE field for the instrumented area. Grey circles represent devices in the instrumented area.

Figure 5

 5 Figure 5.11 -Mean wind direction frequency distributions from measurements over two years (April 2012 -March 2014) at three levels: 3 m (NE and SE), 10 m (10mSW) and 30 m (30mSE), selected for the periods for which L M O measured by 10mSW is (a) neutral |L M O | > 200 m (10-min average data number: 36067) and (b) stable 0 < L M O < 200 m (10-min average data number: 22174).

Figure 5

 5 Figure 5.12 -Mean wind direction from measurements over two years (April 2012 -March 2014) at three levels: 3m (NE and SE), 10m (10mSW) and 30m (30mSE), selected for the periods for which mean wind direction measured by 30mSE is (a) 40 • -50 • and (b) 315 • -325 • , and for 0 <L M O < 200 m measured by 10mSW. The dashed lines extend from the lowest to highest values of the data with "diamond" at the median. The box indicates the 25th and 75th percentiles of the data.

Figure 5

 5 Figure 5.13 -Same as Fig.5.11, selected for the periods during (a) summer (June, July, August) (10-min average data number: 17902) and (b) winter (Dec., Jan., Feb) (10-min average data number: 22787).

Figure 5

 5 Figure 5.14 -Same as Fig.5.12, selected for the periods for which mean wind direction measured by 30mSE is (a) 40 • -50 • during summer, (b) 315 • -325 • during summer, (c) 40 • -50 • during winter, and (d) 315 • -325 • during winter.

Figure 5 .

 5 Figure 5.15 -Mesh for the simplified geometry, with a uniform horizontal resolution of 5 m and same vertical resolution as in Fig. 5.4. Circles indicate the extracted points for vertical profiles.

Figure 5

 5 Figure 5.16 -Horizontal mean wind direction vertical profiles for several extracted points: (a) ext 2 and (b) ext 3, obtained from simulations with simplified geometry for the forest of α = 0.2, 0.5, 0.9, 2 and 10 m -1 , and the solid obstacle.

Figure 5

 5 Figure 5.17 -Horizontal mean wind speed vertical profiles for several extracted points: (a) ext 2, (b) ext 3, (c) entrance and (b) exit, obtained from simulations with simplified geometry for the forest of α = 0.2, 0.5, 0.9, 2 and 10 m -1 , and the solid obstacle.

Figure 5

 5 Figure 5.18 -TKE vertical profiles for several extracted points: (a) ext 2, (b) ext 3 and (c) center, obtained from simulations with simplified geometry for the forest of α = 0.2, 0.5, 0.9, 2 and 10 m -1 , and the solid obstacle.

Figure 5 .

 5 Figure 5.19 -Vertical profiles of vertical velocity w for several extracted points: (a) ext 2, (b) ext 3 and (c) center, obtained from simulations with simplified geometry for the forest of α = 0.2, 0.5, 0.9, 2 and 10 m -1 , and the solid obstacle.

Figure 5

 5 Figure 5.20 -Dissipation rate ε vertical profiles for several extracted points: (a) ext 2, (b) ext 3 and (c) center. obtained from simulations with simplified geometry for the forest of α = 0.2, 0.5, 0.9, 2 and 10 m -1 , and the solid obstacle.

Figure 5 .

 5 Figure 5.21 -Vertical cross-sections of TKE (m 2 .s -2 ) at the center of the domain along the S-N line, obtained from simulations with simplified geometry for the forest of (a) α = 0.2 m -1 , (b) α = 2 m -1 , and (c) the solid obstacle.

Figure 5 .

 5 Figure 5.22 -Vertical cross-sections of vertical velocity (m.s -1 ) at the center of the domain along the S-N line, obtained from simulations with simplified geometry for the forest of (a) α = 0.2 m -1 , (b) α = 2 m -1 , and (c) the solid obstacle.

Figure 5 .

 5 Figure 5.23 -Vertical cross-sections of pressure (Pa) at the center of the domain along the S-N line, obtained from simulations with simplified geometry for the forest of (a) α = 0.2 m -1 , (b) α = 2 m -1 , and (c) the solid obstacle.

Figure 5

 5 Figure 5.24 -Horizontal cross-sections of pressure (Pa) at height of 3 m, obtained from simulations with simplified geometry for the forest of (a) α = 0.2 m -1 , (b) α = 2 m -1 , and (c) the solid obstacle.

Figure 6

 6 Figure 6.2 -'Quick Look' at concentration data after elimination of negative values and with the baseline for the IOP-7 from 40th to 50th minutes; from top to bottom: PID-1, -2, -3, -4, -5 and -background; see the whole IOP in Fig. A.12.

Figure 6

 6 Figure 6.3 -'Quick Look' at concentration data after the baseline removal for the IOP-7 from 40th to 50th minutes; from top to bottom: PID-1, -2, -3, -4, -5 and -background; see the whole IOP in Fig. A.13.

Figure 6

 6 Figure 6.5 -'Quick Look' at concentration data after the baseline removal for the IOP-11 from 20th to 30th minutes; from top to bottom: PID-1, -2, -3, -4, -5, -6 and -background; see the whole IOP in Fig. B.6.

Figure 6

 6 Figure 6.6 -'Quick Look' at concentration data after the baseline removal for the IOP-11 from 70th to 80th minutes; from top to bottom: PID-1, -2, -3, -4, -5, -6 and -background; see the whole IOP in Fig. B.6.

Figure 6

 6 Figure 6.7 -'Quick Look' at wind velocity v (m/s) filtered by central moving average over 1 min of 5 anemometers at "sonic arc at 50 m" during the IOP-11 for two 10min periods; from top to bottom: 20N, 10N, 0, 10S, 20S; see the whole IOP in Fig. B.3.

Figure 6

 6 Figure 6.8 -Concentration histogram for the IOP-7; from top to bottom: PID-1, -2, -3, -5 and -background; the axis-y indicates the percentage of data number.

Figure 6

 6 Figure 6.10 -Normalized concentration spectra f.S(f )/σ 2 C for PID-1 (a), -2 (b), -3 (c), -5 (d) and -background (e) during the IOP-7, and for (f) four sample concentration spectra (continuous lines) from different points on a cross-section through a plume at 740 m range shown in [Mylne and Mason, 1991]; dashed lines on the graph indicate the -2/3 slope representing the inertial subrange.

  Simulations for the pollutants dispersion are performed under the different condition of IOPs by using the transport equations for the mean concentration C (equation 4.21) and the concentration fluctuations c 2 = σ 2 C (equation 4.22) previously presented in the Sect. 4.3.3. To do so, extra scalars are added during simulations with Code_Saturne: a scalar representing the mean concentration C and a scalar representing the variance of the mean concentration σ 2

  Figure 6.11 -Horizontal and vertical cross-sections of mean concentration C for simulations of the IOP-7 with (a) α = 0.9 m -1 and (b) α = 2 m -1 ; horizontal crosssections (top) are at height of 3 m, with horizontal wind vectors and grey circles representing the devices; vertical cross-sections (bottom) are at the location of the source along a E-W line, with grey circles representing the source (3 m), PID-3 (3 m) and the future additional PID (10 m) over PID-3; color maps are in logarithmic scale.

Figure 6 .

 6 Figure6.12 -Comparison of measured mean concentration C with simulations for α = 0.9 and 2 m -1 for the IOP-7.

Figure 6 Figure 6

 66 Figure 6.13 -Horizontal and vertical cross-sections of concentration fluctuation σ C for simulations of the IOP-7 with R f =0.5 and σ C =1%C at the source for (a) α = 0.9 m -1 and (b) α = 2 m -1 ; cross-sections, vectors and grey circles are identical with Fig. 6.11; color maps are in logarithmic scale.

Figure 6 .Figure 6

 66 Figure 6.15 -Horizontal and vertical cross-sections of (a) mean concentration C and (b) concentration fluctuation σ C for simulations of the IOP-11 with α = 0.5 m -1 , R f =0.5 and σ C =1%C at the source; horizontal cross-sections (top) are at height of 3 m, with horizontal wind vectors and grey circles representing the devices; vertical cross-sections (bottom) are at the location of the source along a E-W line, with grey circles representing the source (3 m), PID-3 (3 m) and PID-6 (10 m) above PID-3; color maps are in logarithmic scale.

Figure 6

 6 Figure 6.17 -Comparison of concentration fluctuation σ C for simulations of the POI-11 with α = 0.5 m -1 ; Fig. (a) compares σ C at height of 3 m with measurements of PID-1, -2, -3, -4 and -5; Fig. (b) compares σ C vertical profiles with measurements of PID-3 (3 m) and PID-6 (10 m); blue color represents simulation with R f = 0.8; red color represents simulations with R f = 0.5.

Figure

  Figure A.2 -'Quick Look' at instantaneous wind velocity u (m/s) of 9 anemometers at "sonic square" and "sonic arc at 50 m" during the IOP-7; from top to bottom: NE, NW, SE, SW, 20N, 10N, 0, 10S, 20S.

Figure

  Figure A.4 -'Quick Look' at instantaneous wind velocity w (m/s) of 9 anemometers at "sonic square" and "sonic arc at 50 m" during the IOP-7; from top to bottom: NE, NW, SE, SW, 20N, 10N, 0, 10S, 20S.

Figure

  Figure A.8 -'Quick Look' at wind velocity v (m/s) filtered by central moving average over 1min of 9 anemometers at "sonic square" and "sonic arc at 50 m" during the IOP-7; from top to bottom: NE, NW, SE, SW, 20N, 10N, 0, 10S, 20S.

Figure

  Figure B.2 -'Quick Look' at wind velocity u (m/s) filtered by central moving average over 1min of 9 anemometers at "sonic square" and "sonic arc at 50 m" during the IOP-11; from top to bottom: NE, NW, SE, SW, 20N, 10N, 0, 10S, 20S.

  Also, we can deduce longitudinal and transverse velocity of mean wind from u mean , v mean and dd mean as:

	54		Chapter 3. Wind and turbulence data analysis
	a mean b mean	=	-sin(dd mean ) -cos(dd mean ) cos(dd mean ) -sin(dd mean )	u mean v mean
				u
				v

sin(dd mean ) -cos(dd mean ) cos(dd mean ) -sin(dd mean )

Table 3 .

 3 1 -Statistical values of 12 anemometers calculated from the 60 min sub-period data of the IOP-7 at 3 m, which indicate that the flows are less anisotropic at higher altitude. In addition, L M O increases with height. By calculating the values of z/L M O , we find that the surface layer is slightly less stable at 10 m and 30 m than at 3 m.

Table 3 .

 3 2 -Comparison between mean wind speed and eddy advection speed deduced from spatial cross-correlation, where r a and r b are ratios of the eddy advection velocity to the mean wind speed r = U adv /U for horizontal wind velocity components Table

	r a	r b

Table 3 .

 3 6 -Incident wind sector and direction range for data selection.

	• )

  Table 6.2 -Some statistical values of several anemometers calculated from the 90-min sub-period data of the IOP-11 (See all the statistical values of all the 12 anemometers in Table B.1 )

		NE NW SE SW 20N 20S 10mSW 10mSE 30mSE
	z (m)	3	3	3	3	3	3	10	10	30
	dd mean (deg) 84.2 85.3 90.4 89.9 84.7 95.4	94.1	94.4	97.4
	a mean (m.s -1 ) 1.33 1.42 1.36 1.55 1.44 1.49	2.33	2.30	3.39
	L M O (m)	8	18	10	10	9	17	9	10	47

  Table B.1 -Statistical values of 12 anemometers calculated from the 90-min sub-period data of the IOP-11

	0.6	6.05	5.29	2.08	2.91	2.55	39.71	19.35	7.81
	1.0	8.87	7.69	2.37	3.75	3.25	40.25	16.10	4.14
	1	19.75	13.38	3.70	5.33	3.67	39.82	19.56	3.03
		11.24 8.00 6.63 7.10 9.70 10.74 8.33 8.26 6.22	6.51 5.33 5.61 4.73 6.12 7.44 5.56 4.96 3.56	2.37 1.78 2.04 1.78 1.53 2.48 1.39 1.39 1.33	4.75 4.50 3.25 4.00 6.33 4.33 6.00 5.00 4.67	2.75 3.00 2.75 2.67 4.00 3.00 4.00 3.00 2.67	22.20 18.49 8.98 17.07 35.07 31.89 27.23 26.96 18.66	6.25 6.83 4.49 5.43 8.48 7.64 8.04 6.86 4.48	1.86 1.42 1.50 1.40 2.01 1.65 1.55 1.60 1.19
		σ 2 a /u 2 *	σ 2 b /u 2 *	σ 2 w /u 2 *	σ 2 a /σ 2 w	σ 2 b /σ 2 w	L aa (m)	L bb (m)	L ww (m)

(a) (b) (c)

(a) (b)

†. As pointed out by S. Dupont, rapporteur of the thesis, the origin of the wind rotation in the reference cited is related to the Coriolis force. While in our study, the Coriolis force is neglected.

‡. As pointed out by S. Dupont, rapporteur of the thesis, the origin of the wind rotation in the reference cited is related to the Coriolis force. While in our study, the Coriolis force is neglected.

lent motions and their actions on the plume dispersion. Unfortunately, it will be difficult to generate the accurate inlet conditions based on the measurements inside the modelling area and only at 3 levels (3 m, 10 m and 30 m).

Appendices

Appendix A "Quick Look" at IOP-7

A.1 "Quick Look" at raw wind data

The below figures show 'Quick Look' at raw data measured during IOP-7 by "sonic square" (NE, NW, SE, SW) and "sonic arc at 50 m" (20N, 10N, 0, 10S, 20S). A.2 "Quick Look" at filtered wind data

A.3 "Quick Look" at concentration data

Brief analysis for the IOP-11

The IOP-11 takes place on 11 March 2015, from 18:41 to 21:03. It lasts about 142 min with the gas release for about 122 min (from 18:41 to 20:43). It should be notice that a new PID has been added for this IOP. It is placed at height of 10 m above PID-3. Thus we can have an ideal about the vertical profile of concentration and compare with simulation results.

As for data analysis, a sub-period of 90 min, from 30th min to 120th min, has been chosen with relatively stationary meteorological conditions. This sub-period including 45 0000 values for wind data and 270 000 values for concentration data, corresponds to the measurements from 19:11 to 20:41 on 11 Mars 2015.

B.1 "Quick Look" at filtered wind data -2, -3, -4, -5, -6 and -background. which is consistent with the easterly wind measured by the anemometers at height of 3 m during the IOP-11.

B.3 Statistical analysis for wind and concentration data

Statistical values of wind and concentration are calculated for the 90-min subperiod of the IOP-11 in Table B.1 and B.2. The wind is east at height of 3 m and is slightly south-east for the higher levels. The impact of the forest to the measurements is less remarkable than in the IOP-7 with an north-easterly wind. The wind directional vertical shear is much less important. The difference between anemometers in the south and in the north is also less evident.

Unfortunately, the background concentration is very fluctuating, which make the threshold value difficult to choose and thus conditional concentration variables less accurate.

B.4 Simulation results for wind and turbulence

Figure B.7 compares the simulation results with the measurements for the IOP-11. Since the impact of the forest is much less important for this IOP with an easterly wind, the simulation agrees better with the measurements in the horizontal wind velocity and direction and also has less difference in TKE and vertical velocity comparing to the simulation for the , the wind rotation is correctly reproduced for anemometer 20S and 10mSE but not for 20N, which might be due to the uncertainty of the forest shape when extracting from the satellite map and the land use map. The difference between simulation and measurements in T KE and w might be again due to the inlet profiles which are built with the theoretical profiles and the measurements inside the simulation domain.