Severine Tomas
email: severine.tomas@cemagref.fr

Olivier Eiff

Valery Masson

Experimental investigation of the turbulent momentum transfers in a neutral boundary layer over a rough surface

Keywords: Functional lengths, Integral scales, Neutral boundary layer, Particle image velocimetry measurements, Rough walls, Turbulent budgets

The turbulent characteristics of the neutral boundary layer developing over rough surfaces are not well predicted with operational weather-forecasting models. The problem is attributed to inadequate mixing-length models, to the anisotropy of the ow and to a lack of controlled experimental data against which to validate numerical studies. Therefore, in order to address directly the modelling diculties of the development of a neutral boundary layer over rough surfaces and to investigate the turbulent momentum transfers of such a layer, a set of hydraulic ume experiments were carried out. In the experiments, the mean and turbulent quantities were measured by a particle image velocimetry (PIV) technique. The measured velocity variances and uxes (u i u j ) in longitudinal vertical planes allowed the vertical and longitudinal gradients (∂/∂z and ∂/∂x) of the mean and turbulent quantities (uxes, variances and third-order moments) to be evaluated and the terms of the evolution equations for ∂e/∂t, ∂u 2 /∂t, ∂w 2 /∂t and ∂u w /∂t to be quantied, where e is the turbulent kinetic energy. The results show that the pressure-correlation terms allow the turbulent energy to be transferred equitably from u 2 to w 2 . It appears that the repartition between the constitutive terms of the budget of e, u 2 , w 2 and u w is not signicantly aected by the development of the rough neutral boundary layer. For the whole evolution, the transfers of energy are governed by the same terms which are also very similar to the smooth-wall case. The PIV measurements also allowed the spatial integral scales to be computed directly and to be compared with the dissipative and mixing lengthscales, which were also computed from the data.

Introduction

Roughness eects on turbulent ow over surfaces have been studied since the mid-19th century. For atmospheric boundary layers, which are usually rough, roughness eects have also received particular attention; see for example the review of [START_REF] Raupach | Rough-wall boundary layers[END_REF]. When the roughness elements are high enough, the oset constant in the smooth-wall logarithmic law changes, as discussed in most textbooks (e.g., [START_REF] Schlichting | Boundary layer theory[END_REF]. The roughness eects are accounted for in the logarithmic law via dierently dened but equivalent roughness parameters: (i) the roughness length, used mostly in the atmospheric context, (ii) the equivalent sand roughness, ornally (iii) the roughness function (see [START_REF] Jiménez | Turbulent ows over rough walls[END_REF]). Yet, an important and long-standing question is whether the outer-layer ow (z + = u * z/ν > 50, where u * is the friction velocity and ν the kinematic viscosity (e.g. [START_REF] Pope | Turbulent Flows[END_REF] depends on the wall roughness. Some studies, investigating the effect of roughness on the turbulence structure [START_REF] Antonia | The response of a turbulent boundary layer to a step change in surface roughness. part 1. smooth to rough[END_REF][START_REF] Bandyopadhyay | Rough-wall turbulent boundary layers in the transition regime[END_REF][START_REF] Raupach | Rough-wall boundary layers[END_REF], conclude that there is an outerlayer similarity over smooth and rough walls. In this context, [START_REF] Townsend | The structure of turbulent shear ow[END_REF] similarity hypothesis is often invoked as an implicit "wall-similarity" [START_REF] Raupach | Rough-wall boundary layers[END_REF]) deduced from the general Reynolds number similarity [START_REF] Townsend | The structure of turbulent shear ow[END_REF]. The recent experimental studies of [START_REF] Flack | Experimental support for towsend's reynolds number similarity hypothesis on rough walls[END_REF] and [START_REF] Schultz | The rough-wall turbulent boundary layer from the hydraulically smooth to the fully rough regime[END_REF] support the outer-layer similarity of smooth-and rough-wall boundary layers in terms of both the mean ow and the Reynolds stresses. Similarly, the combined experimental and numerical investigation of [START_REF] Krogstad | An experimental and numerical study of channel ow with rough walls[END_REF] suggests an outer-layer very little aected by the roughness elements. However, an outer-layer structure unaected by the wall condition is not always found. For example, the experimental study of [START_REF] Krogstad | Comparison between roughand smooth-wall turbulent boundary layers[END_REF] showed that the roughness at the wall inuences the mean velocity and the turbulent stresses into the outer-layer.

In a recent review, [START_REF] Jiménez | Turbulent ows over rough walls[END_REF] readdressed this issue. He stressed that in addition to the roughness Reynolds number usually considered, i.e. Re * = k s u * /ν, where k s is the equivalent sand roughness, care needs to be taken to account for the fact that the roughness sublayer will increasingly aect a signicant portion of the logarithmic zone for low blockage ratios, δ/z h , where δ is the boundary-layer height and z h the height of the rough elements. [START_REF] Jiménez | Turbulent ows over rough walls[END_REF] estimated that to be free from direct roughness eects, δ/z h must be greater than at least 40. Thus, in order to attain fully rough regimes, the roughness Reynolds number needs to be greater than about 80. He concludes that the Reynolds number Re δ + = δu * /ν (also called Kármán number δ + ) should be greater than 4000. These are conditions that are usually satised in real atmospheric boundary-layer ows, with the possible exception of urban boundary layers.

Owing to the development of remote sensing technologies, in-situ measurements have become accurate enough to investigate the turbulence of the atmospheric boundary layer (e.g. [START_REF] Poulos | Cases-99: A comprehensive investigation of the stable nocturnal boundary layer[END_REF][START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF]Kunkel and[START_REF] Kunkel | Study of the near wall turbulent region of the high reynolds number boundary layer using an atmospheric ow[END_REF]. The experiments of [START_REF] Kunkel | Study of the near wall turbulent region of the high reynolds number boundary layer using an atmospheric ow[END_REF] over the western Utah Great Salt Lake Desert covered very high Kármán numbers (δ + ∼ = 10 6 ) with very small relative roughness heights (k s /δ ∼ = 10 -4 ). The results show good agreement with the similarity formulations [START_REF] Marusic | Streamwise turbulence intensity formulation for at-plate boundary layers[END_REF] based on Townsend's attached eddy model (1976) and thus support Townsend's outer-layer similarity. As the authors suggest, some of the observed outer-layer dierences in comparison with other studies might be attributed to blockage eects and/or low Reynolds numbers.

Yet, most boundary-layer studies, especially those concerned with roughness eects, only consider the fully developed boundary layer at equilibrium. The development of boundary layers is usually only considered in the context of a transition of well-developed boundary layer subject to a sudden change of wall roughness, also referred to as internal boundary layer, recent examples including those of Cheng and Castro (2002a;2002b) and [START_REF] Belcher | Adjustment of a turbulent boundary layer to a canopy of roughness elements[END_REF]. One exception is the study of [START_REF] Castro | Rough-wall boundar layers: mean ow universality[END_REF] who examined boundary layers developing in wind tunnels over various roughnesses at relatively high Reynolds numbers. However, he focused only on the mean ow for which he concluded that outer-layer ow similarity holds even for surprisingly low blockage ratios, down to δ/z h > 5, (i.e. signicantly smaller than the criterium suggested by [START_REF] Jiménez | Turbulent ows over rough walls[END_REF]). He did not examine the turbulence structure nor the functional and integral lengthscales.

A better knowledge of the turbulence structure including the functional and integral lengthscales is necessary to validate and improve numerical simulations. This includes 1D boundary-layer prediction models which are computationally inexpensive but need to be improved. For example, some 1D models still use simple expressions for the mixing length in neutral boundary layers, such as the Blackadar mixing length [START_REF] Blackadar | The vertuical distribution of wind and turbulent exchange in a neutral atmosphere[END_REF]) l = κzl 0 /(κz + l 0 ) where l 0 is an asymptotic value of the order of a few hundred metres which are not appropriate. This mixing length only depends on height, and not on the ow (it is equal to κz near the ground and to l 0 high above). The large values in the free troposphere predicted by this model are not valid, however, because the stable stratication tends to strongly limit the size of the eddies. Some mixing lengths tend to correct this with stability modulating the mixing length (Bougeault and Lacarrere, 1989), as in the mesoscale meteorological model Meso-NH [START_REF] Lafore | The Meso-NH atmospheric simulation system. Part I: adiabatic formulation and control simulations[END_REF] or the operational model AROME [START_REF] Bouttier | Seminar on Recent Developments in Data Assimilation for Atmosphere and Ocean, 8-12 September, ECMWF, European Center for Medium-Range Weather Forecasts Bouttier F[END_REF][START_REF] Ducrocq | Le projet arome[END_REF][START_REF] Bouttier | Seminar on Recent Developments in Data Assimilation for Atmosphere and Ocean, 8-12 September, ECMWF, European Center for Medium-Range Weather Forecasts Bouttier F[END_REF][START_REF] Seity | The arome-france convective scale operational model. Monthly Weather Review In press Spalart PR[END_REF]. Unfortunately, such mixing length parameterization does not correctly take into account the neutral case. It is not dened if there is no thermal inversion above the boundary layer. Indeed, at the top of the atmospheric boundary layer, there is usually a zone of high stability which separates the boundary layer from the free troposphere above. This is the reason why the mixing lengths are often based on a criterion of stability: in stable conditions, vortices are locally inhibited and mixing lengths are respectively small. In neutral or in unstable conditions, the mixing lengths are relatively large based on the assumption that vortices reach the boundarylayer top as dened by the area of high stability. While this assumption is true in unstable conditions, it is more questionable in neutral conditions where the turbulent structure is only determined by the surface friction and not related to the temperature gradient. Consequently, in a neutral boundary layer, the model overestimates the mixing length. Thus, since no current formulation of the mixing length takes into account the specic features of the neutral boundary layer and its streamwise vortices, except near the surface [START_REF] Redelsperger | A simple general subgrid model suitable both for surface layer and free-stream turbulence[END_REF][START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF]), there appears to be no formulation which satisfactorily predicts the complete structure of neutral boundary layers. This is, in part, due to the lack of data under known and controlled conditions. [START_REF] Hess | Evaluating models of the neutral, barotropic planetary boundary layer using integral measures. part ii. modelling observed conditions[END_REF] used in-situ observations to evaluate the rst-order models of [START_REF] Lettau | Theoretical wind spirals in the boundary layer of a barotropic atmosphere[END_REF] and [START_REF] Blackadar | The vertuical distribution of wind and turbulent exchange in a neutral atmosphere[END_REF][START_REF] Blackadar | A single layer theory of the vertical distribution of wind in a baroclinic neutral atmospheric boundary layer[END_REF] and the higher-order models of [START_REF] Freeman | Transport-dissipation analytical solutions of the e -turbulence model and their role in predictions of the neutral abl[END_REF] and Xu and Taylor (1997a;b) in the neutral case. While the rst-order models generally performed well, by tting a free parameter, the second-order models performed poorly. [START_REF] Hess | Evaluating models of the neutral, barotropic planetary boundary layer using integral measures. part ii. modelling observed conditions[END_REF] attribute this behaviour to deviations from the idealized conditions the higher-order models are based upon. In summary, it can be concluded that there is a need for better experimental data verication under controlled conditions, with access to the turbulent structure of the boundary layers.

For this purpose, it is proposed here to experimentally investigate the development of a neutral boundary layer under a potential free-stream ow over a rough surface. The study aims to provide an idealized and simple case for analysis and validation, for which no data-set under controlled and known conditions seems to be available. In order to attain high enough Reynolds numbers while retaining small enough relative roughness heights [START_REF] Jiménez | Turbulent ows over rough walls[END_REF]), the ow was simulated in the large water ume of the French National Center of Meteorological Research Center in Toulouse. Vertical two-dimensional velocity elds, spanning the development of the boundary layer, were measured via a particle image velocimetry (PIV) technique. The resulting experimental data describe the inertial and outer-layer, i.e. above the surface layer, of a fully rough and developing neutral boundary layer, in terms of the mean and turbulent quantities under known conditions with a high level of accuracy. The measured velocity variances and uxes (u 2 , w 2 and u w ) in a longitudinal vertical plane allow the vertical and longitudinal gradients (∂/∂z and ∂/∂x) of the mean and turbulent quantities (uxes, variances and third-order moments) to be evaluated. They also enable the advection, the dynamical production, the turbulent transport, the dissipation and pressure-correlation terms of the ∂e/∂t, ∂u 2 /∂t, ∂w 2 /∂t and ∂u w /∂t evolution equations to be estimated. Spatial correlation analysis permits the longitudinal and vertical integral lengths to be estimated and to be compared to dissipation lengths and mixing lengths evaluated directly from the data. The intent is to improve the parameterizations of functional lengths when equilibrium is reached. The data should also help asses LES simulations or other closure models simulations of developing atmospherictype rough wall neutral boundary layers.

The experimental procedure is presented in Sect. 2 and the boundarylayer development is discussed in Sect. 3. Then, the study focuses on the turbulent energy budgets and the transfers between the Reynolds stresses. The spatial integral scales as well as the dissipative and mixing lengths are discussed and compared in Sect. 5. The main conclusions are summarized and discussed in Sect. 6.

2

Experimental procedure

Boundary-layer ow facility

The experiments were carried out in a large, wide and high hydraulic ume (22 m long, 3 m wide and 1 m high) (Fig. 1a) to obtain high Reynolds numbers, to generate a two dimensional ow and to minimize other boundary eects. The ume was equipped with a series of grids to reduce the incident turbulence intensity and homogenize the incident ow, such that it can be considered uniform and potential. A one-metre long smooth surface behind the last grid was installed to allow the uniform and potential ow with a small smooth and laminar boundary layer to approach the roughness elements. A turbulent and rough boundary layer was then allowed to develop over a 12.3 m long and 3 m wide horizontal rough surface. The temperature of the water was maintained within ± 0.2 K, assuring neutral ow.

A Laser-Doppler-Velocimetry (LDV) technique was used for an initial estimation of the rough-wall boundary-layer thickness and for the validation of the ume characteristics. Vertical proles of the mean velocity and turbulent intensity were obtained at several abscissaes (x = [-0.16; 0.1; 6; 10.7; 12.4] m) for three transversal y-positions (y = [-0.5; 0; 0.5] m). For the incident free-stream velocity chosen (U ∞,0 = 0.35 m s -1 ), the turbulent intensity and the velocity proles for the three y-positions are essentially superimposed for each x-position. The mean deviation of these proles is about 1 % for both the mean velocity and for the intensity suggesting that the mean ow is essentially two dimensional. This is later conrmed by a very small residue of 5 the divergence ∂U /∂x + ∂W /∂z obtained from the PIV results (about 0.001 s -1 ). Also, the LDV measurements reveal that the free-stream turbulence intensity, dened by 1/2( u 2 /U + w 2 /W ), is less than 1 %, such that freestream ow can be regarded as potential.

Roughness elements

The roughness elements are constituted by LEGO c blocks 31.7 mm long, 15.7 mm wide and 9.6 mm high (Fig. 1b) to create a fully rough regime (Re * = u * z h /ν ∼ = 200). The arrangement (Fig. 1c) was chosen to represent a medium roughness length, z 0 , the roughness parameter usually used in the atmospheric context and dened by the logarithmic law:

U u * = 1 κ ln z -z d z 0 (1)
where κ is the Von Kármán's constant, xed at 0.4, and z d the zero-plane displacement length. The roughness length is a hydrodynamic quantity related to the drag on the roughness elements, which is usually estimated via purely geometrical parameters in terms of the plane area density (λ p = A p /A t ) and the frontal area density (λ f = A f /A t ) (e.g. [START_REF] Raupach | A wind-tunnel study of turbulent ow close to regularly arrayed rough surfaces[END_REF][START_REF] Coceal | A canopy model of mean winds through urban areas[END_REF]). A t is the area of the reference domain (here D x × D y ), A p is the plane area of the obstacles viewed from above and contained in the domain of reference and A f is the frontal area of each obstacle exposed to the wind (see Fig. 1b, Fig. 1c and Table 1 for the dimensions).

For later comparison with the measured data and for the initial design of the surface layout, several models to estimate the roughness length z 0 as well as the displacement length z d , based on the above mentioned geometrical densities, were evaluated. The results, based on the models of [START_REF] Kondo | Aerodynamic roughness over an inhomogeneous ground surface[END_REF]Yamazawa (1986), Bottema (1995;1997), [START_REF] Lettau | Note on aerodynamic roughness parameter estimation on the basis of roughness element description[END_REF] and [START_REF] Macdonald | An improved method for estimation of surface roughness of obstacle arrays[END_REF], are given in Table 2. The estimations of z d , proposed by [START_REF] Bottema | Parameterisation of aerodynamic roughness parameters in relation to air pollutant removal eciency of streets[END_REF]1997) and [START_REF] Macdonald | An improved method for estimation of surface roughness of obstacle arrays[END_REF], where an empirical drag coecient is also taken into account, are similar. The purely geometrical methods, developed by [START_REF] Kondo | Aerodynamic roughness over an inhomogeneous ground surface[END_REF] and [START_REF] Lettau | Note on aerodynamic roughness parameter estimation on the basis of roughness element description[END_REF] seem to over-estimate the modelling of z 0 in comparison with those of [START_REF] Bottema | Parameterisation of aerodynamic roughness parameters in relation to air pollutant removal eciency of streets[END_REF]1997) and [START_REF] Macdonald | An improved method for estimation of surface roughness of obstacle arrays[END_REF]. The latter values of z 0 , are close to the desired value of z 0 /δ for a medium roughness surface.

PIV measurements

The PIV experiments were performed at a free-stream velocity incident on the roughness element, U ∞,0 , of 0.35 m s -1 . In order to characterize the upstream ow and the boundary-layer development over the rough surface, measurements were performed at six x-positions. The ow elds, 0.5 m high and 0.6 m wide in order to capture and exceed the maximum boundary-layer depth, were centered at x = [0.2, 1, 3, 6, 8, 11] m with Re x = U ∞ x/ν = [5 × 10 4 , 3 × 10 5 , 1 × 10 6 , 2 × 10 6 , 3 × 10 6 , 4 × 10 6 ], respectively. The images were taken in the center of the channel (y = 0 m) with a X-STREAM VISION 10 bit CMOS camera with a 1260 × 1024 pixel resolution and equipped with a NIKKOR camera lens with a 50 mm focal length. Image acquisition was controlled by a frequency generator which allows the acquisition of N pairs of images at 1 Hz. A typical time interval between the two images of a pair, ∆t, was about 0.017 s. A high number of pairs was chosen to assure convergence of the turbulent statistics (N = 999). The water was seeded with 60 µm polyamide particles (Orgasol c 2000) illuminated with a laser light sheet 5 mm thick, created with an oscillating mirror and a 25 Watt Spectra Physics Beamlok Argon-Ion laser.

The PIV images were processed with the algorithms and interface developed by [START_REF] Fincham | Low cost, high resolution dpiv for measurement of turbulent uid ow[END_REF] and [START_REF] Fincham | Advanced optimization of correlation imaging velocimetry algorithms[END_REF], which permitted to eectively eliminate the peak-locking bias errors. A resolution of 2 pix mm -1 with a correlation box size of 20 pixels and a grid spacing of 10 pixels in the x-direction, and a correlation box size of 15 pixels and a grid spacing of 8 pixels in z-direction, yields a longitudinal resolution of about 5 mm and a vertical resolution of 4 mm, respectively. This resolution is roughly equivalent to the laser light thickness, thus giving an essentially isotropic measurement resolution (without signicant loss of out-of plane moving particles, estimated at less than 10 %, well below the limit of 30 % suggested by [START_REF] Fincham | Low cost, high resolution dpiv for measurement of turbulent uid ow[END_REF]).

Longitudinal spatial averaging of the measured velocity elds at a given x-position was performed to further increase statistical convergence. It is, however, only justied for the further downstream x-positions (x > 3 m). At the beginning of the rough surface (x = 0.2 m), the velocity eld is subdivided into three parts (smooth, transition, rough) and x-averaged for each one. The uctuations (u i ) in each case are the dierence between the raw eld (u i ) and the spatially and temporally averaged eld (U i ) which is given by:

U i (z) ≡ 1 N N p=1 1 n x nx q=1 u i (q, z, p) , ( 2
)
where N is the number of computed velocity elds (N = 999) and n x is the number of velocity vector columns used (n x = 90, for a complete eld).

The integral time scale can be estimated to be larger than the advection time (i.e., Taylor's hypothesis with u /U << 1 where advection dominates). Therefore, a conservative estimate based on the maximum advection time and the sampling frequency of 1 Hz yields a roughly 50 % overlap of the large scales between consecutive instantaneous velocity elds, or about 2000 independent samples at large scales. This is larger than the 1000 large scales suggested by [START_REF] Tennekes | A rst course in turbulence[END_REF] for example to ensure statistical convergence. This sampling is higher than the Direct Numerical Simulations (DNS) of [START_REF] Coceal | Mean ow and turbulence statistics over groups of urban-like cubical obstacle[END_REF] who averaged over only 400 large scales. It should also be noted that the increased spatial dispersion as the top of the roughness elements is approached (i.e. the roughness sublayer), prevents representative mean statistics to be evaluated with single vertical plane measurements. If meaningful statistics in the roughness sublayer were desired, several vertical planes across a typical roughness area would need to be measured to allow spatial (double-) averaging in all horizontal planes (x, y). This was not the objective of this study, so the measurements and statistical evaluations were limited to above the roughness layer (taken as z > 2z h , e.g. [START_REF] Macdonald | Modelling the mean velocity prole in the urban canopy layer[END_REF].

3 Boundary-layer development

Mean velocity eld

In this section, the mean velocity elds as the boundary layer develops over the rough surface are presented. The log-law scaling parameters used to normalize the velocity and length quantities are also determined.

Fig. 2a-b show the temporally averaged longitudinal and vertical velocity elds, U (x, z) and W (x, z), respectively, around the beginning of the rough bed (x = 0 m). In the longitudinal velocity eld U (x, z) for x > 0, the development of the rough boundary layer can be seen. It can also be seen that the incident smooth boundary layer is relatively thin, as desired, of the order of the roughness elements (δ smooth /z h ∼ = 1). The ensuing ow can thus be considered as a new boundary layer which starts to develop on the roughness elements under a potential and uniform free stream, as opposed to under an existing boundary layer. The eect of the surface roughness discontinuity, at x = 0, can clearly be seen in the strong W perturbation centred around x = 0, with a strong positive increase in W from eectively 0 to about 0.03 m s -1 . This rise is symmetric with respect to the discontinuity at x = 0 m, with a radius of inuence of about 0.18 m or about 10z h . Also, since around the discontinuity the ow over the obstacles must rst accelerate as it is deviated, it causes an initial acceleration of U (∂U /∂x > 0) with an accompanying boundary-layer height decrease, as is indeed observed in Fig. 2a. Further downstream (x > 0.03 m), the new boundary-layer growth dominates due to the drag exerted, and ∂U /∂x < 0 while ∂δ/∂x > 0. Fig. 3 shows the mean longitudinal velocity proles obtained via time and spatial averaging at each measurement station (A-F). Because of the nite height of the water in the ume, the boundary-layer growth is associated with a slight increase of the free-stream velocity (U ∞ ) to balance the ow rate due to the increasing displacement thickness in the boundary layer. This increase in U ∞ is weak, about 5 % over all measurement stations, and in agreement with the preliminary LDV measurements (Sect. 2.1). The mean streamwise velocity prole in the log-law region given by Eq. 1, can be used to determine the roughness length z 0 and the friction velocity u * . To evaluate the displacement height z d independently, we chose [START_REF] Macdonald | Modelling the mean velocity prole in the urban canopy layer[END_REF] semi-empirical proposition:

z d = 1 + α -λp (λ p -1) z h ( 3
)
where α is an empirical coecient, C d the drag coecient xed at 1.2 with a corrective factor β. For cubical and staggered obstacles, [START_REF] Macdonald | An improved method for estimation of surface roughness of obstacle arrays[END_REF] give α = 4.43 and β = 1 which yields z d = 0.76 z h . [START_REF] Macdonald | Modelling the mean velocity prole in the urban canopy layer[END_REF] has veried the results of [START_REF] Raupach | A wind-tunnel study of turbulent ow close to regularly arrayed rough surfaces[END_REF] and concluded that the log-law region, characterized by Eq. 1, is valid from z = 2.3z h to z = 3.5z h . Thus, by linear regression in the region 2.3z h < z < 3.5z h , for x = [6, 8, 11] m, u * and z 0 was obtained (Table 3). The roughness length z 0 and the friction velocity u * , for all values of x > 6 m were found to be approximately constant: 0.25 mm ± 0.05 mm and 0.02 m s -1 ± 0.001 m s -1 , respectively. The value of z 0 agrees with the semi-empirical prediction established by [START_REF] Macdonald | An improved method for estimation of surface roughness of obstacle arrays[END_REF] reported earlier in Table 2. For all values of x < 6 m (positions A to C), estimating u * and z 0 with the log-law is not expected to work since the log-law region is too small for relatively large z h /δ, as discussed earlier.

Another way to determine u * , independent of the log-law, is by assuming u 2 * = -u w z=z d . To evaluate the Reynolds stress in the canopy at z = z d which was not measured at this height, we can either extrapolate the linear stress variation in the outer layer or the essentially constant-stress part below (z/δ < 0.2), as shown in a typical Reynolds stress prole in Fig. 4a. Here, we chose to take an average of the obtained minimum and maximum u * values, shown in Fig. 4b (squares) with those obtained from the log law (triangles). It can be seen that the friction velocity decreases with fetch until an equilibrium value is reached, as already observed by [START_REF] Rao | The structure of the two dimensional internal boundary layer over a sudden change of surface roughness[END_REF] for a boundary layer over a change of surface roughness. In the case of a newly developing boundary layer under a freestream ow, a detailed discussion of the evolution of u * is given by [START_REF] Castro | Rough-wall boundar layers: mean ow universality[END_REF]. It can be noted that the ow is fully rough as desired, since z + = z 0 u * ν ≥ 2, or equivalently k + s = ksu * ν ≥ 70 with k s = z 0 /0.033 (e.g. [START_REF] Jiménez | Turbulent ows over rough walls[END_REF]).

Vertical development of the boundary layer

The boundary-layer height, δ(x), has been determined via two methods. One is the height where U reaches 99 % of the free-stream velocity (U ∞ ). The other is the height where the ux -u w decreases to 5 % of its maximum value. The boundary-layer heights based on these two denitions and normalized by z 0 are shown in Fig. 5. They can be seen to be in close agreement although the height based on the mean velocity is consistently higher than the height based on the turbulent ux, a sign of both the robustness and the arbitrary nature of both criteria. This lead to the use of the average of both criteria to t the evolution of the boundary-layer height, yielding:

δ z 0 ( x z 0 ) = 0.18 ( x z 0 ) 0.79 -203. (4)
As expected, this law diers from the classic smooth turbulent boundarylayer growth rate which predicts a power law exponent of 0.85 and a coecient about 0.08 [START_REF] White | Viscous uid ow. 2nd ed[END_REF]. Our power law exponent is close to the well-known 4/5th power law, however, starting with Elliot's model (1958), to predict the growth of an internal boundary layer developing after a sudden change of roughness. Numerous experimental studies of developing internal boundary layers (e.g. [START_REF] Antonia | The response of a turbulent boundary layer to a step change in surface roughness. part 1. smooth to rough[END_REF][START_REF] Pendergrass | Dispersion in neutral boundary layer over a step change in surface roughness i. mean ow and turbulence structure[END_REF] as well as numerical ones (e.g. [START_REF] Lin | The eect of surface roughness on ow structures in a neutrally stratied planetary boundary layer ow[END_REF]) have yielded exponents close to 0.8. Towsend (1966) theoretical analysis shows that for a large change in friction velocity, very smooth to very rough, the ow behaves essentially as that of a boundary layer developing below a uniform free-stream velocity, i.e. the situation studied here. However, the coecient and depth of the boundary layer measured is less than predicted by Elliot's theory (1958) which is valid for small changes in friction velocity. Indeed, our measured coecient of 0.18 is less than the 0.68 predicted by Elliot with an estimated smooth roughness length of z 0 1 of 3 × 10 -5 m [START_REF] Perrier | Veine hydraulique du cnrm: fonctionnement en écoulement neutre. Note de travail de l'etablissement d'etudes et de recherches météorologique n 216[END_REF]. Also, our coecient is still smaller than the 0.35 predicted by [START_REF] Pendergrass | Dispersion in neutral boundary layer over a step change in surface roughness i. mean ow and turbulence structure[END_REF] for a roughness change parameter M = ln(

z 0 1 z 0 2 ) = -2.3 which is close to our M = -2.1.

Turbulent quantities

We now turn our attention to the evolution of the turbulent quantities. Fig. 6a-c show the proles of u 2 /u * 2 , w 2 /u * 2 and -u w /u * 2 , respectively, as a function of z/δ where δ is the local mean value of the boundary-layer height (Fig. 5) and u * is the local value obtained from the Reynolds stress prole, for all measured x-positions. The normalized intensity of these variances and covariances, at last initially, increases with x. According to the study of [START_REF] Antonia | The response of a turbulent boundary layer to a step change in surface roughness. part 1. smooth to rough[END_REF] again of an internal boundary-layer growth, the increase of u 2 /u * 2 is linked to the intensication of the vertical gradient of U which is accompanied with an increase of dynamical production. Here, the stresses' growth with x tends to decrease until an equilibrium is reached where the proles overlap. This equilibrium is reached between x = 3 m and x = 6 m (Re x = 1 × 10 6 and 2 × 10 6 ) for all three stresses. Both u 2 /u * 2 and -u w /u * 2 are maximum near the top of roughness elements. This feature can be related with the existence of turbulent structures which are generated by the roughness elements [START_REF] Kim | Turbulence statistics in fully developed channel ow at low reynolds number[END_REF]Castro et al. 2006). The existence of these structures is supported by the experimental work of Hommena and Adrian (2003). From these maxima, u 2 /u * 2 decreases linearly up to the top of the outer layer. The decrease of -u w /u * 2 is quasi-linear, as expected for a pressure driven ow. Yet, near the rough wall the shear stress remains roughly constant as usually observed for neutral rough boundary-layer measurements (Cheng andCastro 2002a and2002b;[START_REF] Chow | Explicit ltering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer ow[END_REF][START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF]). This might be attributed, as the roughness elements are approached, to a lack of sucient spatial resolution in this high intensity region and also some possible remaining spatial dispersion not accounted for in single plane measurements as the roughness layer is approached. Fig. 6b shows that the w 2 /u * 2 proles in equilibrium (for x > 6 m) have a maximum near z/δ ∼ = 0.3 but vary very weakly between 0 < z/δ < 0.5 around a value of about 1. The quasi-constancy of w 2 /u * 2 until 0.5δ agrees with the proles observed by [START_REF] Panofsky | The atmospheric boundary layer below 150 meters[END_REF], [START_REF] Yaglom | Similarity laws for wall turbulence ows[END_REF] and Castro et al. (2006) and simulated by [START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF] and [START_REF] Coceal | Mean ow and turbulence statistics over groups of urban-like cubical obstacle[END_REF]. Near the surface, (z/δ ∼ = 0.1), the normalized variances lie between 4 and 4.5 for u 2 /u 2 * and around 0.8 for w 2 /u 2 * . These values are slightly lower than those deduced from the in-situ measurements and the LES results of [START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF] which lie between 5 and 6 for u 2 /u 2 * and 1 to 2 for w 2 /u 2 * . Fig. 6a and 6b reveal that w 2 /u * 2 is about four times smaller than u 2 /u * 2 , i.e. the ow is strongly anisotropic, at large scales. It is also the ratio established by [START_REF] Moeng | A comparison of shear-and buoyancy-driven planetary boundary layer ows[END_REF] or the observations of [START_REF] Nicholls | Aircraft observations of the structure of the lower boundary layer over the sea[END_REF] and [START_REF] Grant | Observations of boundary layer structure made during the 1981 kontur experiment[END_REF]

who measured u 2 /w 2 ∼ = 4.
It can be concluded that this experiment captures the turbulent structure of a developing rough neutral boundary layer, which reaches equilibrium between Re x ∼ = 1×10 6 and 2×10 6 (3 < x < 6 m). In order to investigate the transfer processes governing such a developing boundary layer, the budgets of the turbulent kinetic energy e, u 2 , w 2 and u w are determined in the next section.

Normalized budgets of the second-order turbulent quantities This section examines the budgets of the second-order turbulent transport equations to help evaluate the numerical models and characterize the transfers between the turbulent quantities. The Reynolds-stress transport equation is given by Eq. 5 where Coriolis and thermal production eects are not considered,

∂ ∂t (u i u j ) = -U k ∂u i u j ∂x k ADV -u k u i ∂U j ∂x k + u k u j ∂U i ∂x k DP - 1 ρ 0 u i ∂p ∂x j + u j ∂p ∂x i P C
(5)

-2 ν ∂u i ∂x k ∂u j ∂x k DISS - ∂u k u i u j ∂x k T R + ν ∂ 2 u i u j ∂x k ∂x k DIF F
.

The index k ranges from 1 to 3 and refers, respectively, to the longitudinal, transverse and vertical component. The terms on the right hand side represent: advection (ADV ), dynamical production (DP ), pressurecorrelation (P C), dissipation (DISS), turbulent transport (T R) and molecular diusion(DIF F ). It may be noted that what is often referred to as the dissipation in Eq. 5 is in fact the pseudo-dissipation, ε, dened by:

ε = ν ∂u i ∂x j ∂u i ∂x j . ( 6
)
It can be related to the true dissipation, ε by:

ε = ε + ν ∂ 2 u i u j ∂x i ∂x j . (7)
The diusion term can be neglected since it is the product of the viscosity (10 -6 m 2 s -2 ) by second-order derivatives. As discussed in Sect. 2.1, under the present two-dimensional conditions, the transverse gradients (∂/∂y) of the time-averaged quantities and the transverse mean velocity (V ) can also be neglected. Also, for x-stations at 1 m and further, the vertical mean velocity (W ) is negligible as discussed in Sect. 3.1. Finally, since the ow is statistically stationary, ∂u i u j /∂t = 0.

Turbulent kinetic energy evolution equation

In addition to the above general approximations, in the budget of the turbulent kinetic energy (e = u i u i /2), the pressure-correlation term (P C) can be neglected as it is expected to be a relatively weak transport term as for smooth boundary layers [START_REF] Pope | Turbulent Flows[END_REF]. This implies:

P C u 2 + P C v 2 + P C w 2 ∼ = 0 . (8)
Eq. 5 applied to e, thus reduces to

0 = -U ∂e ∂x ADV -u 2 ∂U ∂x -u w ∂U ∂z DP -ε DISS - ∂u e ∂x - ∂w e ∂z T R . ( 9
)
ε cannot be estimated directly from the measurements, at least with reasonably accuracy, rst, because the transverse gradients have not been measured and second, because the spatial resolution, 5 mm is too large. It is estimated to be about 16 dissipation scales (η) via the turbulence intensities in Fig. 6 and the later integrale-scale computations (Sect. 5.1). This is not suciently small to include the peak of the dissipation spectrum, peaking at about 24η [START_REF] Pope | Turbulent Flows[END_REF]. However, since ε is the only unknown quantity of the measured budget Eq. 9, it can be estimated via the residual. Also, for evaluating the turbulent kinetic energy, e an assumption needs to be made for v 2 . The measurements of Cheng and Castro (2002b) and [START_REF] Macdonald | Physical modelling of urban roughness using arrays of regular roughness elements[END_REF] as well as the simulations of [START_REF] Castro | Turbulence over urban-type roughness: deductions from wind-tunnel measurements[END_REF] show that v 2 ∼ = 0.5(u 2 + w 2 ), while the in-situ experiments of [START_REF] Drobinski | The structure of near-neutral surface layer[END_REF][START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF] and the simulations of [START_REF] Moeng | A comparison of shear-and buoyancy-driven planetary boundary layer ows[END_REF] suggest that the coecient is closer to 0.4. Here, the coecient of 0.5 is retained since because the studies which support this approximation have a conguration close to our experiment, notably for the characteristics of the canopy (z 0 /z h and u * /U ∞ ). It can also be noted that this uncertainty on v 2 aects the estimation of the turbulent transport and the advection terms by a factor of 7 %, and the pseudo-dissipation term by about 4 %.

The turbulent kinetic energy budget is shown in Fig. 7a. Above z/δ ≥ 0.1, advection is a minor term, even smaller than the also minor vertical turbulent transport term. However, the vertical turbulent transport is not negligible at the top of the boundary layer, for 0.6 < z/δ < 1. Indeed, above z/δ > 0.6, although weak compared to the dissipation and the dynamical production in the lower layer, it becomes the only signicant source here. Advection also grows in importance, becoming a minor sink term peaking near z/δ ∼ = 0.8. This balance between the turbulent transport and the advection agrees with the simulations of [START_REF] Mason | Large-eddy simulation of the neutral static stability planetary boundary layer[END_REF] of the neutral boundary layer.

In the outer layer, but below z/δ ∼ = 0.6, the dynamical production essentially balances the dissipation as observed by [START_REF] Glendening | Large eddy simulation of internal boundary layers created by a change in surface roughness[END_REF] in LES of an internal boundary layer. The preponderant contribution of the dynamical production has also been assessed by [START_REF] Drobinski | The structure of near-neutral surface layer[END_REF] and more recently by [START_REF] Castro | Turbulence over urban-type roughness: deductions from wind-tunnel measurements[END_REF] and [START_REF] Burattini | Comparison between experiments and direct numerical simulations in a channel ow with roughness on one wall[END_REF]. It is the classic assumption which stipulates that the production rate equals the dissipation.

It also supports the outer layer similarity hypothesis, the budgets being very similar as for smooth boundary layers [START_REF] Pope | Turbulent Flows[END_REF]. Fig. 7b shows the longitudinal evolution of the dynamical production term of e, the most signicant term of the budget of e measured directly. It can be seen that for x = 1 m and x = 3 m (positions B and C) the dynamical production is higher than for the further x-positions in the range 0.3 < z/δ < 0.7. This source "surplus" is balanced by an increase in the sink of the vertical turbulent transport of e (not shown) which decays by x = 6 m. By x = 6 m at the latest, all the proles have converged (including the other terms, not shown) in accordance with the observation with regard to the turbulent stresses (Fig. 6a-c) that equilibrium is reached in the region 3 m < x < 6 m.

u 2 evolution equation

For stationary and two dimensional ows, Eq. 5 applied to u 2 reduces to

0 = -U ∂u 2 ∂x ADV -2 u 2 ∂U ∂x + u w ∂U ∂z DP - 2 ρ 0 u ∂p ∂x P C (10) -2 ν ( ∂u ∂x ) 2 + ( ∂u ∂z ) 2 ε u 2 - ∂u 3 ∂x - ∂w u 2 ∂z T R
.

The dissipation term (ε u 2 ) and the pressure-correlation (P C) term are the two unknown terms in the budget. Thus, in order to balance the budget, an additional hypothesis needs to be made. [START_REF] Pope | Turbulent Flows[END_REF], based on DNS data of Spalart (1988) of a boundary layer over a at smooth plate, shows that close to the wall, the anisotropy in the dissipation rate (ε u i u j ) of the Reynolds stresses, u i u j , is clearly large. However, for z/δ > 0.1, there is approximate isotropy. The small level of anisotropy in ε u i u j for z/δ > 0.1 can be attributed to the relatively low Reynolds number of the DNS of Spalart (1988). For rough-wall boundary layers, the experiments of [START_REF] Saddoughi | Local isotropy in turbulent boundary layer at high reynolds number[END_REF] and [START_REF] Saddoughi | Local isotropy in complex turbulent boundary layer at high reynolds number[END_REF] conrm that at this altitude there is also a local isotropy of ε u i u j at high Reynolds numbers. Some studies even conclude that the roughness increases the degree of isotropy closer to the wall and this impact is greater with three-dimensional roughness [START_REF] Castro | Turbulence over urban-type roughness: deductions from wind-tunnel measurements[END_REF]. Therefore, as the present experiment was performed at high Reynolds number (Re δ + = δu * /ν > 3 × 10 3 at x = 6 m) and with three-dimensional roughness, the dissipation rate (ε u i u j ) of u i u j can reasonably be expressed as

ε u i u j = 2 3 ε δ ij , (11)
where δ ij is the Kronecker function. Thus, taking ε u 2 = 2/3ε, with ε obtained from Eq. 9, the pressure-correlation is the only residual term of the budget of u 2 . The estimation of ε u 2 and consequently of the pressurecorrelation term is thus again essentially aected by the uncertainty on v 2 by a factor inferior to 4 %. The budget of u 2 (Fig. 7c) shows, as expected, that the intensity of all the terms decreases as z/δ increases. The turbulence acts principally in the lower part of the outer layer. Above the surface layer, the dynamical production is the only signicant source as shown by [START_REF] Pope | Turbulent Flows[END_REF] for a smooth boundary layer. The advection by mean wind or turbulent transport are not really important for the evolution of u 2 , except at the top of the outer layer where it peaks near z/δ ∼ = 0.8 and where the turbulent transport becomes the dominant source, as for the budget of e. The pressure-correlation term is signicant throughout. Indeed, it is the main sink for u 2 above the surface layer. Similarly to boundary layers over a smooth surface (Spalart 1988), the pressure uctuations appear to distribute the energy between the dierent components: from u 2 to v 2 and w 2 . This is conrmed in the next section. Fig. 7d shows that equilibrium of the most signicant term in the u 2 budget, the dynamical production term, is reached between x = 3 m and x = 6 m, as for the budget of e and in agreement with the longitudinal evolution of u /u * 2 (Sect. 3.3).

w 2 evolution equation

For stationary and two dimensional ows, Eq. 5 applied to w 2 reduces to

0 = -U ∂w 2 ∂x ADV - 2 ρ 0 w ∂p ∂z P C (12) -2 ν ( ∂w ∂x ) 2 + ( ∂w ∂z ) 2 ε w 2 - ∂u w 2 ∂x - ∂w 3 ∂z T R
.

Here, the dynamical production term depends on W and can thus be neglected. Also, as ε u 2 for the u 2 budget, ε w 2 is estimated by Eq. 11. The pressure-correlation term (P C) is then again the only residual term of the measured budget of w 2 . The w 2 budget (Fig. 7e) shows that the turbulent transport term (T R) is again weak, peaking near z/δ ∼ = 0.8 as a source similar to u 2 and e but smaller in magnitude. The advection, also weak, peaks near z/δ ∼ = 0.8 as a sink similar to u 2 and e. The main balance, below z/δ ∼ = 0.8, is between the pressure-correlation and the dissipation term. In the outer layer, around z/δ = 0.8, the main source for w 2 is thus the pressure-correlation term as seen in the recent computations of [START_REF] Ashraan | Roughness eects in turbulent channel ow[END_REF]. This observation conrms the previous hypothesis: the pressure uctuations re-distribute the energy from u 2 to w 2 . It can also be observed that the w 2 pressure-correlation term attains half the absolute value of the u 2 pressurecorrelation term, suggesting that the energy is equally partitioned between w 2 and v 2 . This is conrmed by considering that the sum of the pressurecorrelation terms of u 2 , v 2 and w 2 is zero (Eq. 8) since the pressurecorrelation term of e is assumed to be zero (see Sect. 4.1), which allowed P C v 2 to be evaluated. Yet, it is important to insist on the fact that the pressure-correlation terms of u 2 and w 2 are individually and independently estimated: each one is the residue of their respective budget for the ve downstream x-positions. Finally, Fig. 7f shows the longitudinal evolution of P C w 2 term. In agreement with all other terms it can be seen that equilib- rium is reached for x > 3 m. Fig. 8a shows all three pressure-correlation terms at x = 6 m while Fig. 8b shows the ratio -P C u 2 /P C w 2 . It can be seen in Fig. 8a that P C v 2 and P C w 2 are well superposed for all z/δ. The mean value of the ratio -P C u 2 /P C w 2 in Fig. 8b reveals a ratio close to 2 for 0 < z/δ < 0.8. For z/δ > 0.8, the ratio surprisingly decreases close to 1, but the magnitude of the terms is probably too small to get a signicant estimation of the ratio. Awaiting more precise measurements for z/δ > 0.8, it can thus be concluded that the pressure uctuations equitably distribute the energy contained in the longitudinal uctuations to the transverse and vertical uctuations:

P C v 2 ∼ = P C w 2 ∼ = - 1 2 P C u 2 . ( 13
)
This redistribution, observed for all x-positions, is essentially the same as has been observed over smooth boundary layers [START_REF] Pope | Turbulent Flows[END_REF], suggesting that the turbulent structure in the outer layer is not inuenced by the roughness elements.

u w evolution equation

For u w , under the stationarity and 2D assumptions, Eq. 5 reduces to

0 = -U ∂u w ∂x ADV -u 2 ∂U ∂x + u w ∂U ∂z DP (14) - 1 ρ 0 u ∂p ∂z + w ∂p ∂x P C - ∂u 2 w ∂x - ∂u w 2 ∂z T R
, where the dissipation term ε u w has been neglected (see Eq. 11). The pressure-correlation term is thus the residual term here. The budget of u w is presented in Fig. 7g. Fig. 7g shows that the advection and turbulent transport terms are again relatively weak but with reversed signs compared to u 2 , w 2 and e, and are still peaking near z/δ ∼ = 0.8. Below, the main balance is between the dynamical production (source) and the pressure-correlation (sink), similar to the budget of u 2 where dissipation contributes as a sink. This result agrees with the measurements of [START_REF] Wyngaard | Atmospheric turbulence[END_REF] and the simulations of [START_REF] Ashraan | Roughness eects in turbulent channel ow[END_REF]. [START_REF] Wyngaard | Atmospheric turbulence[END_REF] shows that the turbulent transport is negligible and that the pressure-correlation acts as a sink term which locally balances the dynamical production term. This repartition of energy was also observed by [START_REF] Mulhearn | A wind-tunnel boundary layer study of the eects of a surface roughness change: rough to smooth[END_REF] in wind-tunnel modelling of a rough-to-smooth internal boundary-layer type transition.

In conclusion, all the budgets reach equilibrium between stations at x = 3 m and x = 6 m while the distribution of the budget terms is very similar to smooth turbulent boundary layers and internal boundary layers.

5

Length scales

In this section we investigate the dissipative lengthscale l ε and the mixing length l m based on the measured statistics, as well as the integral lengthscales based on spatial correlations. An accurate quantication of these scales is essential as they drive the modelling of the ow. l ε ∼ = u 3 /ε and l m ∼ = (K/(∂U /∂z)) 1/2 are functional lengthscales, based on the ratio of statistical properties of the ow [START_REF] Hunt | Cross-correlation and length-scales in turbulent ows near surfaces[END_REF], the former often being used to estimate the latter which involves the sought-after Reynolds stress. Here, u is a characteristic turbulent velocity, and K is the turbulent eddy viscosity or exchange coecient given by K = -u w /(∂U /∂z) for horizontally homogeneous ow. First-order models are based on the so-called Prandtl-Kolomogorov relation, which assumes K ∼ = ul, where u is most often taken as e 1/2 and l is a lengthscale of the energy containing eddies, well described by the integral scale, but usually taken from a functional lengthscale. The aim of this section is to characterize the turbulent scales of a neutral boundarylayer ow and to improve parameterization of the functional lengths.

Integral lengthscales

Integral lengthscales can be obtained directly from spatial correlations with the present PIV measurements, without relying on the often necessary Taylor hypothesis. In particular, the spatial correlation function in the (x, z) plane,

R u i u j ( r, δr) = u i ( r)u j ( r + δr) u i ( r)u j ( r) , (15)
can be estimated directly, where r = x i + z k is the reference position and δr taken as either δx i or δz k is the displacement, which is either positive or negative. The correlation functions are computed in each of the N = 999 uctuating (x, z) velocity elds with r for all measured points in the velocity elds. The resulting correlation functions are then averaged over all N for convergence.

The integral lengthscales at each reference x-position are then computed via the integral of the correlation functions (Eq. 15):

L u i u j ,δ r ( r) = ∞ 0 R u i u j ( r, δr)d δr. ( 16
)
To assure systematic convergence, the integration for large r was stopped when the correlation function reaches 0.1. The measurements permit 16 dierent integral scales to be obtained, i.e. for the four stresses (u 2 , w 2 , u w and w u ), each in two directions (δx i and δz k), with δx and δz either positive or negative. As expected, it was observed that δx either positive or negative does not impact on L u i u j ,δx ( r), underlining that the ow is homogenous in x, at least at the scale of the measurement elds. However, for L u i u j ,δz ( r), the results revealed a general tendency of larger integral scales on the top half of the boundary layer for negative δz displacements than for positive ones, and vice-versa on the bottom half. Since [START_REF] Carlotti | Length scales in wall-bounded high-reynoldsnumber turbulence[END_REF] argued that the smallest vertical integral scales should be close to the mixing length, the minimum integral-scale value between the positive and negative δz displacements at each z-level (for all four stresses) was chosen. The resulting integral scales are regrouped in Fig. 9a-h for all x-positions. It can be seen that equilibrium of the integral scales is also reached between 3 m < x < 6 m, as for the turbulent statistics and budgets. Nevertheless, the rst x-positions (x = 0.2 m and x = 1 m) exhibit the same tendencies except for L u 2 ,δx . It can also be noted that L u 2 ,δx and L u w ,δx reveal a greater gap between x = 3 m and x = 6 m. This illustrates that the u uctuations are more aected by the growing of the boundary layer than w which are more associated with the local turbulence. Also, it can be observed that L u 2 ,δx dominates all other lengthscales, as expected. This dominance of L u 2 ,δx is due to the existence of coherent structures which are elongated in the longitudinal direction (e.g. [START_REF] Drobinski | Numerical and experimental investigation of the neutral atmospheric surface layer[END_REF]). These structures are localized in the bottom of the boundary layer with a maxima at z/δ = 0.2 which can also be observed in the L u w ,δx(z) proles. For 1D modelling purpose, it is more appropriate to consider the vertical scales, given by δz displacements. As seen in the left column of Fig. 9, all the vertical integral scales increase linearly until a maximum value is reached around L u i u j ,δz /δ = 0.4 in the range of 0.3 < z/δ < 0.6. Above, the vertical integral scales decrease to L w 2 ,δz /δ = 0.2, L u w ,δz /δ = 0.2 and L w u ,δz /δ = 0.1 at z/δ = 1. Only, L u 2 ,δz increases again above z/δ = 0.9, which underlines the fact that signicant exchanges of longitudinal structures occur between the boundary layer and the free stream. The vertical velocity uctuation L w 2 ,δz and L w 2 ,δx are very close however, suggesting that the lengthscales of w 2 are essentially isotropic. Fig. 10a compares directly the proles of the resulting vertical integral-scales of the four stresses. Here, it can be seen more clearly that the vertical integral lengthscales essentially have the same shape and intensity, increasing in the bottom half of the boundary layer and decreasing above. More precisely, below z/δ ≤ 0.2, in the inertial sublayer, all proles collapse very well, but above there are some dierences. Also, as expected due to the inhomogeneity in the z-direction, the integral scales of the u w and w u uxes are not equal, with a marked shift in the maximum scale.

To compare the integral lengthscales with l m and l ε used in 1D prediction models at equilibrium, a representative integral lengthscale (L z ) was taken as the average of the four vertical lengthscales in Fig. 10a. This mean integral lengthscale (L z ), shown in Fig. 10b, is linear in the lower part of the outer layer, as expected. Above, L z tends to be hyperbolic. The maximum value reached by integral scale is L z /δ ∼ = 0.3 at z/δ = 0.5. A function of the form az/(b + z k ) is thus expected to be adapted, as the dashed line in Fig. 10b shows. The tted normalized equation is:

L z δ = A z-z h δ -s B + z-z h δ -s p , (17)
with A = 0.26; s = 0.05; p = 2.4 and B = 0.24.

Dissipative lengths

In the following sections, we consider only the region where the ow has reached equilibrium, i.e. for x > 6 m. [START_REF] Taylor | Statistical theory of turbulence: Parts i-iii[END_REF] rst suggested that ε ∼ = u 3 /l, usually written as:

ε = C ε e 3/2 l ε . ( 18
)
where C ε is a constant largely xed at 0.845 (Schmidt and Schuman 1989;[START_REF] Cuijpers | Large eddy simulation of trade wind cumulus clouds[END_REF][START_REF] Canuto | New third-order moments for the convective boundary layer[END_REF]Cheng et al. 2002). Using the dissipation rate obtained from the budget of e (Eq. 9), l ε can thus be determined via Eq. 18. Fig. 11a shows the resulting vertical prole of l ε normalized by the experimental integral length L z of Fig. 10b. Clearly, except between 0.3 < z/δ < 0.7, where l ε ∼ = 2.4L z , l ε /L z is not constant. This is not surprising since the two scales should only be proportional in homogeneous turbulence.

Mixing length

The mixing length model,

u w = -C m l 2 m | ∂U ∂z | ∂U ∂z , (19)
where C m is a constant, is used to model the Reynolds stress when one lengthscale dominates. It can be used here to deduce the mixing length l m from the measured mean velocity gradient and Reynolds stress prole, also for x > 6 m. Fig. 11b shows the resulting normalized prole l m C

1/2 m /L z . In comparison with the dissipative length, l /L z , the mixing length l m matches L z better. Neglecting the strong nonlinearities near the surface for z ≤ 0.1δ, the scales are approximately proportional over a wider range, between 0.3 < z/δ < 0.9, with a ratio of proportionality (l m C 1/2 m /L z ) of about 0.37.

The standard mixing length in the logarithmic layer,

l m = κ(z -z d ), (20)
in combination with l m C 1/2 m determined from Eq. 19, permits C m to be determined in the log-law region via linear regression, yielding an average C m = 0.51 for the three latest x-stations. Thus, the ratio of proportionality l m /L z in the zone 0.3 ≤ z/δ ≤ 0.9 is about 0.52. Fig. 12 shows l m /δ with C m = 0.51 as well as the standard mixing length in the logarithmic layer (Eq. 20). The formulation of [START_REF] Blackadar | The vertuical distribution of wind and turbulent exchange in a neutral atmosphere[END_REF],

l m = l 0 κz κz + l 0 , (21)
has been also given in Fig. 12, by determining l 0 by a least squares t of l m for the last three x-stations. The resulting three l 0 /δ ratios are very close, with an average of 0.24. It can be seen that the normalized measurements for all three stations follow essentially the same curve, which might be described as a weak S. In the central outer-layer, between 0.3 ≤ z/δ ≤ 0.7, l m /δ is largely constant at about 0.13, before retreating slightly and then increasing as the top of the boundary layer is approached. Not surprisingly, Blackadar's relation used for all stability conditions does not t the data well. A parametrization based on either the dissipative scale or the vertical integral scale would also yield signicant dierences, so direct parametrization of the mixing length based on the measured l m proles in Fig. 12 might be most appropriate.

Conclusion

The main goal of this study was to investigate the turbulence of the neutral boundary layer to help validate numerical simulation models and to help improve the parameterizations of the turbulent processes for one-dimensional neutral boundary-layer models. The experiment presented in this paper reproduces the development of a neutral boundary layer in a water ume over a rough surface at high Reynolds numbers. The velocity measurements were performed via a two-dimensional particle image velocimetry (PIV) technique in the vertical symmetry plane of the ow. The turbulent statistics at several longitudinal positions were obtained by spatially and temporally averaging the PIV velocities elds, assuring statistical convergence for the higher-order statistics.

The boundary-layer growth was established on the basis of the mean longitudinal velocities and the turbulent shear stress, yielding a 0.8 power law in agreement with previous results for the growth of internal boundary layer.

However, it appears that our 0.18 coecient is smaller than the one of an internal boundary layer. The beginning of the roughness surface is well revealed by mean vertical velocity perturbation whose streamwise extend is about 0.18 m. An analysis of the budget of turbulent kinetic energy (e) shows that shear production is the main source of turbulence in the outer layer, as for turbulent boundary layers over smooth surface. The residual term of this budget yielded the dissipation rate. This allowed also the dissipation of the u 2 and w 2 budgets to be estimated. Finally, we were able to quantify the pressure uctuations and describe its role in transferring energy between the dierent variances. Again, these transfers were found to mirror the exchanges due to the pressure-correlation terms in smooth boundary layers, suggesting the turbulent structure of the outer layer over a rough surface is very similar to one over a smooth surface. The development of the rough neutral boundary layer does not signicantly modify the repartition between the constitutive terms of the budget of e, u 2 , w 2 and u w . Whatever the abscissa, even before equilibrium is reached, the same terms govern the transfer of energy, except for the vertical turbulent transport of e which is a minor term for 0.1 < z/δ < 0.6 at all stations except for x = 1 m where it is not negligible, and where it appears as a sink of energy which balances the "surplus" of dynamical production. All budgets as well the rst and second order statistics reach equilibrium between 12000 ≤ x/z 0 ≤ 24000.

Finally, the spatially resolved measurements allowed us to compute the spatial integral lengthscales directly, including the vertical ones, more appropriate for modelling purposes. These integral length were used to compare with the dissipative and the mixing lengths functional scales. It is shown that the ratio of the dissipative scale to vertical integral lengthscale, not expected to be constant in this inhomogeneous ow, indeed varies signicantly by a fac-tor up to 2. Direct calculation of the mixing length suggest a better match with the vertical integral scale, with a narrow constant of proportionality for z/δ ≥ 0.1 of about l m /L z ∼ = 0.52. Further numerical tests need to be done in order to validate the new proposal of L z and verify whether if it improve the modelling of neutral boundary layer.

Table 1: Dimensions of the reference domain (streamwise D x and crosswise D y ), the spacing between the roughness elements contained in the domain (streamwise W x and crosswise W y ), the dimensions of the roughness elements (height z h , length L t and width L s ), the dierent areas (the area of the reference domain A t , the plane area of the obstacles viewed from above and contained in the domain of reference A p and the frontal area of each obstacle exposed to the wind A f ) and densities (the plane area density λ p and the frontal area density λ f ) (all units in mm). * = -u w z=z d from the u w proles: the minimum value is obtained by averaging the rst points of the prole (for z/δ < 0.2) and the maximum value is given by the intersection between the extrapolation of the linear regression of u w for 0.5 < z/δ < 0.9 and z = z d . (b) Squares symbols: friction velocity u * estimated by method (a) with u * = (u * max + u * max )/2 via -u w shear stress proles and Triangle symbols: friction velocity u * estimated via the log-law ts in the region 2.3z h < z < 3.5z h (see Table 3). Fig. 5: Boundary-layer height evolution (δ) normalized by the roughness length (z 0 ), the square symbols: 5 % of u w max criterion; the X symbols: 99 % of U ∞ criterion; the triangle symbols: average of δ u w 5% and δ U 99% ; the solid line: the best t of the average of both criteria. * = -u w z=z d from the u w proles: the minimum value is obtained by averaging the rst points of the prole (for z/δ < 0.2) and the maximum value is given by the intersection between the extrapolation of the linear regression of u w for 0.5 < z/δ < 0.9 and z = z d . (b) Squares symbols: friction velocity u * estimated by method (a) with u * = (u * max + u * max )/2 via -u w shear stress proles and Triangle symbols: friction velocity u * estimated via the log-law ts in the region 2.3z h < z < 3.5z h (see Table 3).

DOMAIN

29

Tomas S., Eiff O. and Masson V. (2011) Experimental investigation of the turbulent momentum transfers in a neutral boundary layer over a rough surface, Boundary-Layer Meteorology, 138, 385-411 author-produced version of the final draft post-refeering the original publication is available at www.springerlink.com -doi : 10.1007/s10546-010-9566-0

z 0 δ / z x / z 0 δ (x / z ) = 0.18 (x / z ) -203 0 0 0 0.79 δ mean U u'w' 99% 0.5% δ δ
Figure 5: Boundary-layer height evolution (δ) normalized by the roughness length (z 0 ), the square symbols: 5 % of u w max criterion; the X symbols: 99 % of U ∞ criterion; the triangle symbols: average of δ u w 5% and δ U 99% ; the solid line: the best t of the average of both criteria.

36

Tomas S., Eiff O. and Masson V. (2011) Experimental investigation of the turbulent momentum transfers in a neutral boundary layer over a rough surface, 138, author-produced version of the final draft post-refeering the original publication is available at www.springerlink.com -doi : 10.1007/s10546-010-9566-0

FiguresFig. 1 :

 1 Figures Fig. 1: (a) Sketch of the experimental set-up. (b) Typical element, LEGO c . (c) Arrangement of the roughness elements.

Fig 2 :

 2 Fig 2: Mean velocity elds at position A (x = 0 m), near the beginning of the rough bed (a) U (x, z), (b) W (x, z).

Fig. 3 :

 3 Fig. 3: Mean longitudinal velocity proles U (z) at all measurement xpositions (A-F).

Fig. 4 :

 4 Fig. 4: (a) Method used to estimate the friction velocity u 2* = -u w z=z d from the u w proles: the minimum value is obtained by averaging the rst points of the prole (for z/δ < 0.2) and the maximum value is given by the intersection between the extrapolation of the linear regression of u w for 0.5 < z/δ < 0.9 and z = z d . (b) Squares symbols: friction velocity u * estimated by method (a) with u * = (u * max + u * max )/2 via -u w shear stress proles and Triangle symbols: friction velocity u * estimated via the log-law ts in the region 2.3z h < z < 3.5z h (see Table3).

Fig. 6 :

 6 Fig. 6: Locally normalized proles of (a) u 2 /u * 2 , (b) w 2 /u * 2 and (c) u w /u * 2 as a function of z/δ for all measured x-positions.

Fig. 7 :

 7 Fig.7: Locally normalized budget of e (a), u 2 (c), w 2 (e) and u w (g) as a function of z/δ for the D position (x = 6 m). Dominant source terms of each budget: the dynamical production term (DP ) for the budget of e (b), u 2 (d), u w 2 (h) and the pressure-correlation term (P C) for the budget of w 2 as a function of z/δ for all x-positions.

Fig. 8 :

 8 Fig. 8: (a) Locally normalized magnitude of the pressure-correlation terms P C of the budget of u 2 , w 2 and v 2 as a function of z/δ for the D position (x = 6 m). (b) Ratio of the normalized the pressure-correlation terms -P C u 2 /P C w 2 as a function of z/δ for the D position (x = 6 m).

Fig. 9 :

 9 Fig.9: Normalized minimum vertical and longitudinal integral lengthscale L u i u j ,δz (z), in the left column and L u i u j ,δx (z), in the right column, for all x-positions for u 2 (a and b), w 2 (c and d), u w (e and f) and w u (g and h).

Fig. 10 :

 10 Fig. 10: (a) Normalized minimum vertical integral lengthscales for u 2 , w 2 , u w and w u . (b) Mean vertical integral lengthscale (L z ) of (a) and its t with Lz δ = [A( z-z h δ -s)]/[B + ( z-z h δ -s) p ], with A = 0.26; s = 0.05; p = 2.4 and B = 0.24.

Fig. 11 :

 11 Fig. 11: (a) Ratio of the integral lengthscale and the dissipative length (l ε /L z ) at developed x-positions as a function of z/δ. (b) Ratio of the integral lengthscale and the mixing length (l m C 1/2 m /L z ) at developed x-positions as a function of z/δ.

Fig. 12 :Figure 1 :-

 121 Fig.12: Measured mixing length l m /δ deduced from the Eq. 19 (symbols) at developed x-positions as a function of z/δ, compared to the Eq. 20, dashed line, l m = κ(z -d) and Eq. solid line l m = l 0 κz κz+l 0(Blakadar, 1962).

Figure 2 :Figure 3 :

 23 Figure 2: Mean velocity elds at position A (x = 0 m), near the beginning of the rough bed (a) U (x, z), (b) W (x, z).

Figure 4 :

 4 Figure4: (a) Method used to estimate the friction velocity u 2 * = -u w z=z d from the u w proles: the minimum value is obtained by averaging the rst points of the prole (for z/δ < 0.2) and the maximum value is given by the intersection between the extrapolation of the linear regression of u w for 0.5 < z/δ < 0.9 and z = z d . (b) Squares symbols: friction velocity u * estimated by method (a) with u * = (u * max + u * max )/2 via -u w shear stress proles and Triangle symbols: friction velocity u * estimated via the log-law ts in the region 2.3z h < z < 3.5z h (see Table3).

Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 6789101112 Figure 6: Locally normalized proles of (a) u 2 /u * 2 , (b) w 2 /u * 2 and (c) u w /u * 2 as a function of z/δ for all measured x-positions.

Table 2 :

 2 Roughness parameters z d and z 0 and ratio z d /z h and z 0 /z h obtained from the predictions based on analytical theories.

		SPACINGS OBSTACLES	AREAS		DENSITIES
	D x D y	W x W y	z h L t L s	A t	A p	A f	λ p	λ f
	32 61.3	8	8	10 32 16 1962 1024 480 0.5	0.2
		Lettau Kondo et Yamazawa	Bottema	Macdonald et al.
		(1969)		(1986)	(1995, 1997)		(1998)
	z 0 (mm)	1.1		1.2	0.01		0.25
	z 0 / z h	0.11		0.12	0.001		0.025
	z d (mm)	NA		NA	6			7
	z d / z h	NA		NA	0.66		0.76

Table 3 :

 3 Estimation of the roughness parameters, z 0 and u * by linear regression of the log-law in the region 2.3z h < z < 3.5z h and for x > 6 m with z d = 7.6 mm, for U ∞,0 =0.35 m s -1 .

		x (m)
	6	8	11
	u * (m s -1 ) 0.020 0.021 0.019
	z 0 (mm) 0.253 0.346 0.223

Tomas S.,Eiff O. and Masson V. (2011) Experimental investigation of the turbulent momentum transfers in a neutral boundary layer over a rough surface, Boundary-Layer Meteorology, 138, 385-411 author-produced version of the final draft post-refeering the original publication is available at www.springerlink.com -doi : 10.1007/s10546-010-9566-0

Acknowledgements

This work study was supported by the MesoScale Meteorology and Experiment Meteorology groups (GMME and GMEI) of the National Center of Meteorological Research of Météo-France. The authors would like to thank B. Beaudoin, J-C. Boulay, J.-C. Canonici, M. Morera, S. Lassus Pigat and H. Schaner of the CNRM-GAME uid mechanics laboratory (SPEA) for their support during the experiments for providing us with valuable comments and availability.