Darshan Venkatrayappa
email: darshan.venkatrayappa@mines-ales.fr

Philippe Montesinos
email: philippe.montesinos@mines-ales.fr

Daniel Diep
email: daniel.diep@mines-ales.fr

Baptiste Magnier
email: baptiste.magnier@mines-ales.fr

A novel image descriptor based on anisotropic filtering

Keywords: Anisotropic half derivative Gaussian kernel, Anisotropic half derivative exponential kernel, image descriptor, HoG, image matching

In this paper, we present a new image patch descriptor for object detection and image matching. The descriptor is based on the standard HoG pipeline. The descriptor is generated in a novel way, by embedding the response of an oriented anisotropic derivative half Gaussian kernel in the Histogram of Orientation Gradient (HoG) framework. By doing so, we are able to bin more curvature information. As a result, our descriptor performs better than the state of art descriptors such as SIFT, GLOH and DAISY. In addition to this, we repeat the same procedure by replacing the anisotropic derivative half Gaussian kernel with a computationally less complex anisotropic derivative half exponential kernel and achieve similar results. The proposed image descriptors using both the kernels are very robust and shows promising results for variations in brightness, scale, rotation, view point, blur and compression. We have extensively evaluated the effectiveness of the devised method with various challenging image pairs acquired under varying circumstances.

Introduction

Currently, most of the promising method for image content description use local features as their basis. Thus, image feature extraction has become an important and active research topic in the field of computer vision. Local image features form the foundation for many real time applications such as object detection and tracking, panorama stitching, video surveillance and image based retrieval. In addition to feature description these applications require invariance to rotation, scale, viewpoint and illumination changes. Many image description algorithms have been proposed to incorporate all these characteristics. To name a few: SIFT [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF], GLOH [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF], SURF [START_REF] Bay | Speeded-Up Robust Features[END_REF], LIOP [START_REF] Fan | Rotationally Invariant Descriptors Using Intensity Order Pooling[END_REF], ASIFT [START_REF] Morel | ASIFT: A New Framework for Fully Affine Invariant Image Comparison[END_REF], DAISY [START_REF] Tola | DAISY: An Efficient Dense Descriptor Applied to Wide Baseline Stereo[END_REF] and many more. Most of these descriptors are based on Histogram of Oriented Gradients (HoG) [START_REF] Dalal | Histograms of Oriented Gradients for Human Detection[END_REF]. SIFT, SURF and GLOH are robust to rotation and scale invariance characteristics. ASIFT improves on these descriptors by exhibiting good affine invariance properties. LIOP has succeeded in handling both the geometric distortions and monotonic changes in illumination. Some of these descriptors and their variants are compared and reviewed in detail in [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF].

Another category of feature descriptors use filter response as their basis. Schmid and Mohr [START_REF] Schmid | Local Gray-value Invariants for Image Retrieval[END_REF] use differential invariant responses to compute new local image descriptors. Differential invariant responses are obtained from a combination of Gaussian derivatives of different orders, which are invariant to 2dimensional rigid transformations. Koenderink et al. [START_REF] Koenderink | Representation of local geometry in the visual system[END_REF] formulate a methodology based on Gaussian function and its derivatives to capture the local geometry of the image. Larsen et al. [START_REF] Larsen | Jet-Based Local Image Descriptors[END_REF] follow a new approach for the construction of an image descriptor based on local k-jet, which uses filter bank responses for feature description. Palomares et al. [START_REF] Palomares | A new affine invariant method for image matching[END_REF] have come up with a local image descriptor issued from a filtering stage made of oriented anisotropic half-Gaussian smoothing convolution kernels. They achieve euclidean invariance in the matching stage. However, most of the descriptors based on filter response vaguely capture the geometry of the region around key-points, as a result they fail to compete with the HoG based descriptors.

Some researchers capture the geometry of the image patch using curvature information, and various edge related information. Monroy et al. [START_REF] Monroy | Beyond straight lines -Object detection using curvature[END_REF] use descriptor based on curvature histograms for object detection. This descriptor gives good results for images with strong line segments. Authors in [START_REF] Zitnick | Binary Coherent Edge Descriptors[END_REF] use location, orientation and length of the edges to construct an image patch descriptor. Here, the descriptor encodes the presence or absence of edges using a binary value for a range of possible edge positions and orientations. Eigenstetter et al. [START_REF] Eigenstetter | Visual Recognition using Embedded Feature Selection for Curvature Self-Similarity[END_REF] have proposed an object representation framework based on curvature self-similarity. This method goes beyond the popular approximation of objects using straight lines. However, like all descriptors using second order statistics, this approach also exhibits a very high dimensionality.

In our method rotation and affine normalized gray level image patch is extracted as in [START_REF] Mikolajczyk | Scale & Affine Invariant Interest Point Detectors[END_REF]. Then, we construct our descriptor on this image patch. Thus obtained descriptor is completely different than the above mentioned descriptors. It can be seen as a combination of all the above three approaches. In the remaining of this paper for simplicity we refer to anisotropic half Gaussian derivative kernel as AHGDK and anisotropic half exponential derivative kernel as AHEDK. In our descriptor :

1. The response of a rotating AHGDK or AHEDK around pixels is used to construct signatures. 2. The orientation of the edges and the anisotropic gradient directions are extracted from the signatures, thus capturing the geometry/curvature of the image patch. 3. These orientations are binned separately in different ways as in HoG to form different descriptors. 4. Additionally, we construct a faster variant of our descriptor by replacing the AHGDK with a computationally inexpensive AHEDK.

FILTERING STAGE

Anisotropic half Gaussian derivative kernel (AHGDK)

When compared to isotropic filters (Fig. 1(a)), anisotropic filters (Fig. 1(b)) have an advantage in detecting large linear structures. For anisotropic filters, at the corners the gradient magnitude decreases as the edge information under the scope of the filter decreases (Fig. 2(a)). Consequently, robustness to noise, when dealing with tiny geometric structures weakens. This weakness can be nullified by using AHGDK [START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF]. In addition to this, thanks to the elongated and oriented half kernels (Fig. 1(c), Fig. 1(e)), we are able to estimate the two edge directions as in Fig. 2(b). As illustrated in Fig. 1(f), at the pixel coordinate (x, y), a derivative kernel is applied to obtain a derivative information Q(x, y, θ) in a function of orientation θ ∈ [0; 360 • [:

Q(x, y, θ) = I θ * C • H (-y) • x • e - x 2 2λ 2 + y 2 2µ 2 , (1)
where I θ corresponds to a rotated image of orientation θ, C is a normalization coefficient [START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF] and (µ, λ) represent the standard deviation of AHGDK. Only the causal part of this filter along the Y axis is used, this is obtained by cutting the kernel in the middle, in an operation that corresponds to the Heaviside function H. As in [START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF], we have chosen to rotate the image instead of filter, thereby reducing the algorithmic complexity by making use of a recursive Gaussian filter [START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF].

Anisotropic half exponential derivative kernel (AHEDK)

We modify the Shen Castan filter [START_REF] Shen | An optimal linear operator for step edge detection[END_REF] to imitate the AHGDK. We call this modified filter as anisotropic half exponential derivative kernel (AHEDK). AHEDK shows similar characteristics and produces similar results as that of the above explained AHGDK. We use the recursive implementation of half exponential kernel which is of order 1. So, it is approximately 5 times faster than the recursive implementation of anisotropic half derivative kernel. By construction AHEDK exhibits derivative characteristics along the X direction and smoothing characteristics along the Y direction.

E(x, y, θ) = I θ * C 1 • H (y) • e (-ᵕy) • sign(x) • e (-α λ •|x|) (2)
The derivative information E(x, y, θ) in Eq.2 is obtained by spinning the AHEDK around a pixel (x, y), C 1 is a normalization coefficient and (α µ , α λ ) the height and width of the anisotropic half exponential kernel (see Fig. 1(d)). Similar to Eq.1 only the causal part of this filter along the Y axis is used by cutting the kernel in the middle, in an operation that corresponds to the Heaviside function H.

Anisotropic gradient magnitude and direction estimation

We construct the signature Q(x p , y p , θ) by considering the response of the above described rotating filter at a pixel location (x p , y p ) as in Fig. 1(f). Fig. 1(g) shows a sample signature obtained by applying the AHGDK at the the pixel location (x p , y p ) in steps of 2 • . Anisotropic gradient magnitude ∇I a and its associate direction η for the key-point at location (x p , y p ) are obtained by considering the global extrema G max and G min of the function Q(x p , y p , θ) along with θ 1 and θ 2 . The two angles θ 1 and θ 2 define a curve crossing the pixel (an incoming and outgoing direction) thus representing the two edge directions. Two of these global extrema are combined to maximize the gradient ∇I a , i.e:

               G max = max θ∈[0,360 • [ Q(x p , y p , θ) and θ 1 = arg max θ∈[0,360 • [ (Q(x p , y p , θ)) G min = min θ∈[0,360 • [ Q(x p , y p , θ) and θ 2 = arg min θ∈[0,360 • [ (Q(x p , y p , θ)) ∇I a = G max -G min η = θ 1 + θ2 2 (3)
By construction, our filter naturally detects curvature information with θ 1 and θ 2 using only 1st order derivatives. The arrows on the circle in Fig. 2 and on the square in Fig. 1 illustrates this property. We follow the same procedure to obtain the signature E(x p , y p , θ) and extract the edge and gradient directions using the AHEDK.

DESCRIPTOR CONSTRUCTION

Descriptor construction process is shown in Fig. 3. As in [START_REF] Mikolajczyk | Scale & Affine Invariant Interest Point Detectors[END_REF], we follow the standard procedure to obtain the rotation and affine normalized gray level image patch. This standard procedure is followed in the construction of all the descriptors (SIFT, DAISY, GLOH) used in our experiments. As in Fig. 3, for each pixel in the image patch, we spin the AHGDK and obtain a signature (for simplicity and proper viewing, in Fig. 3 signatures are not shown). From this signature we extract : i) Angle at the maxima, θ 1 and the response G max at θ 1 . ii) Angle at the minima, θ 2 and the response G min at θ 2 . iii) Anisotropic gradient angle η and its magnitude ∇I a . The next step is to construct 3 intermediate descriptors by forming the HoG of all the 3 angles separately. The angle θ 1 is weighed by G max , θ 2 by G min , η by ∇I a and binned as in Eq.4. 256 dimension descriptor obtained by concatenating HoG θ1 and HoG θ2 . iv) DESCT4-theta1-theta2-eta : A 384 dimension descriptor obtained by concatenating HoG θ1 , HoG θ2 and HoG η . We follow the same procedure to construct the descriptors using AHEDK and call them EXP-DESCT1-theta1-eta, EXP-DESCT2-theta2-eta, EXP-DESCT3-theta1-theta2 and EXP-DESCT4-theta1-theta2eta.

   HoG θ1 = {θ 1 bin1 , θ 1 bin2 , θ 1 bin3 , θ 1 bin4 .....θ 1 bin128 } HoG θ2 = {θ 2 bin1 , θ 2 bin2 , θ 2 bin3 , θ 2 bin4 .....θ 2 bin128 } HoG η = {η bin1 , η bin2 , η bin3 , η bin4 .........η bin128 } (4)

Experiments and Results

Dataset and Evaluation

The entire code was implemented on Matlab platform. Harris affine key points [START_REF] Mikolajczyk | Scale & Affine Invariant Interest Point Detectors[END_REF] were used for image patch extraction. Key points obtained from other detectors can also be used for extracting the image patches. We evaluate and compare the performance of our descriptor against the state of the art descriptors on the standard dataset using the standard protocol provided by Oxford group. The binaries and dataset are obtained from website linked to [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF] (http://www.robots.ox.ac.uk/ vgg/research/affine/). The dataset used in our experiments has different geometric and photometric transformations such as change of scale and image rotation (boat), viewpoint change (graf), image blur (bike), JPEG compression (compression) and illumination change (Leuven). For each type of image transformation there is a set of six images with established ground truth homographies.

We use the evaluation criterion as proposed by [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF]. The evaluation criteria is based on the number of correspondences, number of correct matches and the number of false matches between two images. Here, we test the descriptors using Similarity based method and Nearest neighbour method. Due to lack of space we restrain from going in to the details of these methods. A detailed description of these methods can be found in [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF]. The results are presented using the recall vs 1-precision curves. As in eq.5, recall is defined as the total number of correctly matched regions over the number of corresponding regions between two images of the same scene. From eq.6, 1-precision is represented by the number of false matches relative to the total number of matches. In all our experiments, Euclidean distance is used as the distance measure. AHGDK height µ and width λ is fixed to 6 and 1 respectively. AHEDK height α µ and width α λ is fixed to 1 and 0.2 respectively. Width and height parameters are chosen to have a ratio sharpness length suitable for robust edge detection [START_REF] Montesinos | A New Perceptual Edge Detector in Color Images[END_REF], which generally gives good results in most cases. This ratio is compatible with the angle filtering step.

Descriptor Performance

The performance of descriptors obtained using both AHGDK and AHEDK are compared against SIFT, GLOH and DAISY. For SIFT and GLOH, the descriptors are extracted from the binaries provided by Oxford group (http: //www.robots.ox.ac.uk/ ~vgg/research/affine/) [START_REF] Mikolajczyk | A Performance Evaluation of Local Descriptors[END_REF]. DAISY descriptor for patches is extracted from the code provided by [START_REF] Tola | DAISY: An Efficient Dense Descriptor Applied to Wide Baseline Stereo[END_REF]. We have compared our descriptor on all the images in the dataset. Due to lack of space, we omit showing the quantitative results for the image pair (1-6). Qualitative results using nearest neighbour approach for the image pair (1-4) using descriptors obtained from both AHGDK and AHEDK is shown in Fig. 7.

Rotation changes (boat).

Using similarity matching strategy, all of our descriptors obtained from both AHGDK and AHEDK have a clear advantage over SIFT, GLOH and DAISY. This can be seen from the 1st row of Fig. 4 and Fig. 5. When the descriptors obtained using AHGDK and AHEDK are compared against each other they exhibit similar results as shown in the 1st row of Fig. 6. The same can be said, when experimented with nearest neighbour matching strategy. Due to lack of space, for all the dataset we omit showing the quantitative results using nearest neighbour matching method.

Viewpoint changes (graff).

For image pair (1-2), (1-3) and (1-4), using similarity approach our descriptors and its variants performs similar to or better than that of SIFT, DAISY and GLOH. Image pair (1-5) is very challenging and all the descriptors including SIFT, DAISY and GLOH fail. The results can be seen from the

-precision

Recall SIFT GLOH DAISY DESCT1-teta1-eta DESCT2-teta2-eta DESCT3-teta1-teta2 DESCT4-teta1-teta2-eta (i) bike 1-2 (k) bike 1-3 (k) bike 1-4 (l)

-precision

Recall SIFT GLOH DAISY DESCT1-teta1-eta DESCT2-teta2-eta DESCT3-teta1-teta2 DESCT4-teta1-teta2-eta (m) compression 1-2 (n) compression 1-3 (o) compression 1-4 (p)

-precision

Recall

SIFT GLOH DAISY DESCT1-teta1-eta DESCT2-teta2-eta DESCT3-teta1-teta2 DESCT4-teta1-teta2-eta (q) Leuven 1-2 (r) Leuven 1-3 (s) Leuven 1-4 (t) Leuven 1-5
Fig. 4. Recall vs 1-Precision curves for SIFT, GLOH, DAISY and 4 descriptors obtained using AHGDK. Similarity matching is used for evaluation.

graphs in the 2nd row of Fig. 4 and Fig. 5. When the descriptors obtained using AHGDK and AHEDK are compared against each other they exhibit similar results as shown in the 2nd row of Fig. 6.

Variations in blur, compression and brightness.

For the blur (bike), compression, and brightness (Leuven) changes, using similarity approach all our descriptors outperform SIFT, DAISY and GLOH. The performance can be seen in the third, fourth and fifth rows of Fig. 4 and Fig. 5 respectively. When the descriptors obtained using AHGDK and AHEDK are compared against each other they exhibit similar results as shown in the 3rd, 4th and 5th row of Fig. 6.

Conclusion

This paper proposes novel image patch descriptors based on anisotropic half Gaussian derivative kernel (AHGDK) and anisotropic half exponential derivative kernel (AHEDK). The originality of our method is that it captures the

Fig. 1 .

 1 Fig. 1. (a) Isotropic gaussian kernel (b) Anisotropic gaussian kernel (c) AHGDK (d) AHEDK with αµ = 1 and α λ = .2 (e) A thin Gaussian derivative half-filter. (f) AHGDK applied to a keypoint (xp, yp) in the image patch. (g) Extrema of a function Q(xp, yp, θ). and ∆θ = 2 o . Note that the initial orientation of the filter is vertical, upwardly directed and steerable clockwise. ((h),(i),(j)) Synthetic square illustrating the image geometry captured by AHGDK. All the figures in this paper are best when viewed in color.

Fig. 2 .

 2 Fig. 2. Synthetic quarter circle used to illustrate the image geometry captured by AHEDK. Enlarge the image to get a proper view of the orientations.

Finally, weFig. 3 .

 3 Fig.3. Histogram construction process using 3 directions θ1, θ2 and η obtained using anisotropic half Gaussian kernel.

6 )

 6 matches No of correct matches + No of false matches (Our descriptors depends on 4 different parameters: ∆θ, N o-of-bins, height and width (µ, λ for AHGDK and α µ , α λ for AHEDK). The rotation step ∆θ is fixed to 10 • . Increasing the rotation step results in loss of information. The image patch is divided into 16 blocks. All blocks are of size 10x10(Since we are using a patch of size 41x41 the blocks in the extreme right and bottom have 11x11 size). The number of bins (N o-of-bins) is fixed to 8 per block, resulting in a 8 * 16 = 128 bins for 16 blocks. Increasing the number of bins results in almost the same performance but, increases the dimensionality of the descriptor.

5 Fig. 5 .

 55 Fig.5. Recall vs 1-Precision curves for SIFT, GLOH, DAISY and 4 descriptors obtained using AHEDK. Similarity matching is used for evaluation.

  (a)DESCT4-theta1-theta2-eta (boat 1-4) (b)EXP-DESCT4-theta1-theta2-eta (boat 1-4) (c)DESCT4-theta1-theta2-eta (Graf 1-4) (d)EXP-DESCT4-theta1-theta2-eta (Graf 1-4) (e)DESCT4-theta1-theta2-eta (bike 1-4) (f)EXP-DESCT4-theta1-theta2-eta (bike 1-4) (g)DESCT4-theta1-theta2-eta (comp 1-4) (h)EXP-DESCT4-theta1-theta2-eta (comp 1-4)(i)DESCT4-theta1-theta2-eta (leuven 1-4) (j)EXP-DESCT4-theta1-theta2-eta (leuven 1-4)

Fig. 7 .

 7 Fig. 7. Qualitative matching results for the Oxford dataset using the nearest neighbour matching method. Matches in the Left column are obtained using AHGDK and the matches in the right column are obtained using AHEDK.

Acknowledgements

This work is funded by L'institut mediterraneen des metiers de la longevite (I2ML), Nimes, France.

Authors Suppressed Due to Excessive Length Fig. 6. Recall vs 1-Precision curves for all descriptors obtained using both AHGDK and AHEDK. Similarity matching is used for evaluation.

: precision

Recall

: precision

: precision

geometry of the image patch by embedding the response of the AHGDK or AHEDK in a HoG framework. Our method incorporates edge direction as well as anisotropic gradient direction for generating the descriptor and its variations.

On the standard dataset provided by the Oxford group our descriptor and its variants outperform SIFT, GLOH and DAISY. In the future, we would like to learn the parameters by introducing a learning stage. We would also like to use this approach for object classification and image retrieval purpose. The speed of the descriptor generation can be boosted by parallel programming.


J Jacob

P Sagaut

Wind comfort assessment by means of large eddy simulation with lattice Boltzmann method in full scale city area

Keywords:

Large-eddy simulations based on the Lattice-Boltzmann method of the flow in a realistic, full scale urban area are performed to compare several wind comfort criteria. It is observed that popular criteria for pedestrian comfort lead to very different conclusions, due to the access to high spatio-temporal resolution data. Different mixed strategies based on the combination of several criteria are proposed and compared to enhance pedestrian wind comfort assessment in practical cases.

Introduction

Pedestrian comfort is a global field of urban physics dealing with wind comfort, pollutant dispersion and thermal comfort close to the ground of cities. It can be addressed using either wind tunnel or in situ measurements that provide data at specific locations or using Computational Fluid Dynamics (CFD) which gives access to data over wide areas with very high spatial resolution. All these approaches are complementary considering that in situ and wind tunnel measurements permit to generate complete database that will be used to validate numerical model and to define guidelines for CFD studies which permit to obtain plenty of data on full scale geometries to assess pedestrian comfort.

Pedestrian wind comfort is sensitive to several kind of parameters such as the local wind (mean velocity, turbulent intensity), the location of the city (atmospheric conditions, building density) or the peoples (age, weight) so it is necessary to find or define universal rules to study it. In the literature those studies are mainly based on a mixing of meteorological data, aerodynamic data and comfort criterion to address local wind comfort in cities. Different criteria have been proposed in the literature, which are observed to a significant dispersion of results in some cases. Ratcliff and Peterka [START_REF] Ratcliff | Comparison of pedestrian wind acceptability criteria[END_REF], Ohba et al. [START_REF] Ohba | Study on the assessment of environmental wind conditions at ground level in a built-up area -based on long-term measurements using portable 3-cup anemometers[END_REF], Bottema [START_REF] Bottema | A method for optimisation of wind discomfort criteria[END_REF] and Koss [START_REF] Koss | On differences and similarities of applied wind comfort criteria[END_REF] listed and compared several wind comfort criteria such as those discussed in Davenport [START_REF] Davenport | An approach to human comfort criteria for environmental wind conditions[END_REF], Gandemer [START_REF] Gandemer | Wind environment around buildings: aerodynamic concept[END_REF][START_REF] Gandemer | Aerodynamic studies of built-up areas made by C[END_REF], Isyumov and Davenport [START_REF] Isyumov | The ground level wind environment in built-up areas[END_REF][START_REF] Isyumov | Evaluation of the effects of tall buildings on pedestrian level wind environment[END_REF], Lawson and Penwarden [START_REF] Lawson | The effects of wind on people in the vicinity of buildings[END_REF], Melbourne [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] and Hunt et al. [START_REF] Hunt | The effects of wind on people; new criteria based on wind tunnel experiments[END_REF]. This variability of pedestrian comfort criteria is significantly impacted by the nature of their input data: time-averaged velocity, turbulent intensity and averaging period length, from a few seconds to a few hours. In order to reduce the uncertainty induced by this high sensitivity of wind comfort criteria, comparing them on the same case can be a good way to evaluate pedestrian comfort quality but many kind of aerodynamic data can be necessary to this end. CFD is then an interesting tool since it permits to assess different data at many locations with a moderate effort.

The use of CFD for urban flow simulation is more and more widespread. Many studies are available in the literature and best practice guidelines [START_REF] Franke | Best practice guideline for the CFD simulation of flows in the urban environment[END_REF][START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF][START_REF] Blocken | Computational Fluid Dynamics for urban physics: importance, scales, possibilities, limitations and ten tips and tricks towards accurate and reliable simulations[END_REF] on the use of CFD for that kind of application have been proposed. Most of existing CFD simulations have been performed using the Reynolds Average Navier Stokes (RANS) approach, which resolves only the mean, time-averaged flow while the turbulent motion is modelled, on simplified [START_REF] Blocken | CFD evaluation of wind speed conditions in passages between parallel buildingseffect of wall-function roughness modifications for the atmospheric boundary layer flow[END_REF][START_REF] Du | Effects of lift-up design on pedestrian level wind comfort in different building configurations under three wind directions[END_REF][START_REF] Yoshie | Cooperative project for CFD prediction of pedestrian wind environment in the Architectural Institute of Japan[END_REF][START_REF] Kim | A numerical study of the effects of ambient wind direction on flow and dispersion in urban street canyons using the RNG k-ε turbulence model[END_REF] or realistic [START_REF] Yoshie | Cooperative project for CFD prediction of pedestrian wind environment in the Architectural Institute of Japan[END_REF][START_REF] Blocken | Modification of pedestrian wind comfort in the Silvertop Tower passages by an automatic control system[END_REF][START_REF] Blocken | Pedestrian wind comfort around a large football stadium in an urban environment: CFD simulation, validation and application of the new Dutch wind nuisance standard[END_REF][START_REF] Blocken | CFD simulation for pedestrian wind comfort and wind safety in urban areas: general decision framework and case study for the Eindhoven university campus[END_REF][START_REF] Janssen | Pedestrian wind comfort around buildings: comparison of wind comfort criteria based on whole-flow field data for a complex case study[END_REF][START_REF] Kang | Development of a computational fluid dynamics model with tree drag parameterizations: application to pedestrian wind comfort in an urban area[END_REF][START_REF] Liu | CFD simulations of wind distribution in an urban community with a full-scale geometrical model[END_REF][START_REF] Montazeri | CFD evaluation of new secondskin facade concept for wind comfort on building balconies: case study for the Park Tower in Antwerp[END_REF][START_REF] Tominaga | Cross comparisons of CFD prediction for wind environment at pedestrian level around buildings. comparison of results for flowfield around building complex in actual urban area[END_REF][START_REF] Van Hooff | Coupled urban wind flow and indoor natural ventilation modelling on a high-resolution grid: a case study for the Amsterdam ArenA stadium[END_REF][START_REF] Zheng | Pedestrian-level wind environment on outdoor platforms of a thousand-meter-scale megatall building: sub-configuration experiment and wind comfort assessment[END_REF] geometries. However time-resolved approaches such as Detached Eddy Simulation (DES, see Ref. [START_REF] Sagaut | Multiscale and Multiresolution Approaches in Turbulence[END_REF]) have been used recently on both simplified isolated buildings geometries [START_REF] Liu | Detached eddy simulation of pedestrian-level wind and gust around an elevated building[END_REF][START_REF] Dadioti | Validation of a datached eddy simulation (DES) model in complex urban environment[END_REF] and complex urban area like Shinjuku AIJ test case [START_REF] Dadioti | Validation of a datached eddy simulation (DES) model in complex urban environment[END_REF]. Steady RANS simulations are widely used because of their low computational cost, but they do not permit to assess unsteady data, and therefore do not allow to consider all existing wind comfort criteria. Another commonly reported weakness is that they over-predict the turbulent energy dissipation, leading to a bad prediction of recirculation bubbles observed at the top or in the wake of buildings. To cure that problem Large Eddy Simulation (LES, see Refs. [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows: an Introduction[END_REF][START_REF] Garnier | Large-eddy Simulation for Compressible Flows[END_REF]) which allows for the direct resolution of a wide range of turbulent frequencies, can be used. However, the computational cost increase associated to LES due to its high spatio-temporal resolution can become a problem when simulations of flow over complex realistic geometries is targeted. These methods were first applied to simplified geometries or reduced area of cities, e.g. Refs. [START_REF] He | Evaluation of pedestrian winds in urban area by numerical approach[END_REF][START_REF] Xie | Large-eddy simulation of flows over random urbanlike obstacles[END_REF][START_REF] Razak | Analysis of airflow over building arrays for assessment of urban wind environment[END_REF][START_REF] Shi | Large eddy simulation of wind field and plume dispersion in building array[END_REF][START_REF] Letzel | High resolution urban large-eddy simulation studies from street canyon to neighbourhood scale[END_REF], but also for simulations of pollutant dispersion over complex urban areas [START_REF] Gousseau | Near-field pollutant dispersion in an actual urban area: analysis of the mass transport mechanism by highresolution Large Eddy Simulations[END_REF], wind loads on buildings [START_REF] Elshaer | LES evaluation of wind-induced responses for an isolated and a surrounded tall building[END_REF] and pedestrian wind comfort assessment [START_REF] Letzel | LES case study on pedestrian level ventilation in two neighbourhoods in Hong Kong[END_REF][START_REF] Adamek | Pedestrian level wind assessment through city development: a study of the financial district in Toronto[END_REF]. A review of urban CFD simulations at pedestrian level is given in Table 1. In order to reduce computational cost of LES simulation it is possible to switch from the classical CFD approach base on the Navier-Stokes equations to Lattice Boltzmann method (LBM) solvers which are very efficient for parallel computations of separated low Mach number flows. This is illustrated in Ahmad et al. [START_REF] Ahmad | Large-eddy simulation of the gust index in an urban area using the lattice Boltzmann method[END_REF], who performed simulations to assess pedestrian level gust index in a × × km km km 19.2 4.8 1.0 area of Tokyo with a m 2 finest grid resolution taking advantage of high efficiency of lattice Boltzmann method for parallel simulation. The LBM method is very interesting for CFD because it is fully local which avoid the use of complex and time consuming numerical methods and the complete algorithm needs only to access the data of the first order neighbors which increase the performance for parallel simulations. Furthermore the computational grid is based on a hierarchy of embedded uniform meshes with a ratio of 2 for the grid step between two successive refinement levels. The use of immersed boundary conditions allows to handle complex geometries such as city in a very easy and automatic way. This is also interesting considering that the ratio dx dt / is kept identical at all grid refinement levels, which means that only the nodes at the finest refinement level are computed every time step reducing the number of floating point operations during the simulation.

The present study deals with the application of an LES-LBM solver to wind comfort assessment at pedestrian level in full scale urban geometry using different wind comfort criteria. The aim is to compare different existing criteria thanks to the high space-time resolution data provided by LES in a realistic configuration, and to check their coherency and robustness with respect to the accuracy of input data. In Section 2 key features of the Lattice Boltzmann method used in this paper are presented. Section 3 presents the validation of the present method on a realistic urban configuration, namely the case F of the Architectural Institute of Japan open database [START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF][START_REF] Tominaga | Cross comparisons of CFD prediction for wind environment at pedestrian level around buildings. comparison of results for flowfield around building complex in actual urban area[END_REF]. Section 4 is devoted to the results obtained dealing with pedestrian wind comfort assessment at a height of 2 m from the ground. Conclusions are given in Section 5.

Numerical method

All the CFD simulations presented here have been carried out using a research version of ProLB [49] that use the Lattice Boltzmann Method [START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and beyond[END_REF][START_REF] Guo | The Lattice Boltzmann Method and its Applications in Engineering[END_REF][START_REF] Kruüger | The Lattice Boltzmann Method. Principles and Practice[END_REF][START_REF] He | Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation[END_REF][START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] to solve fluid dynamics equations. It is based on the resolution of Boltzmann equation (Eq. ( 1)) that describes the evolution of a particle distribution function = → → f f x c t ( , , ) which is related to the prob- ability density of particles with velocity → c at time t and position → x . This equation is solved on a DdQq (d dimensions, q discrete velocities) lattice.

∂ ∂ + → ⋅∇ → = f t c f f Ω( ) (1)
The lattice Boltzmann equation (Eq. ( 2)) with time step t Δ and space step ⎯ → ⎯⎯ = ⎯→ ⎯ x c t Δ Δ α is given by Ref. [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF].

→ + ⎯→ ⎯ + -→ = → f x c t t t f x t x t ( Δ , Δ ) ( , ) Ω ( , ) α α α α (2)
with f α the distribution function following the α direction and Ω α the collision operator. In this work the lattice use is the D3Q19 lattice given by,

→ = ⎧ ⎨ ⎩ = ± ± ± = - ± ± ± ± ± ± = - c α α α
(0,0,0) 0 ( 1,0,0), (0, 1,0), (0,0, 1) 1 6 ( 1, 1,0), ( 1,0, 1), (0, 1, 1) 7 18

α ( 3
)
The solution of Eq. ( 2) is usually computed in two steps using a Strang-type splitting. . The macroscopic quantities (density ρ, momentum ρu i and momentum flux tensor Π ij ) are then reconstructed computing the moments of the distribution functions, i.e.

∑ = = - ρ f α q α 0 1 (4a) ∑ = = - ρu f c i α q α αi 0 1 (4b) ∑ = = - f c c Π ij α q α αi αj 0 1 (4c)
The collision operator, which is interpreted as a relaxation towards an equilibrium state, may involve either a single relaxation time (BGK collision) [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] or several ones, leading to a multiple relaxation time (MRT) method [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. MRT models are expected to be numerically more stable than BGK models, because they allow for a finer tuning of the effective viscosity. Nevertheless, different methods such as collision kernel regularization [START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF] allow to increase the BGK model stability with a very small extra computational cost. Such a regularized collision operator with a single relaxation time expressed in Eq. ( 5) is used in the present work,

= - - = - τ f f τ f Ω 1 ( ) 1 α α α
eq α neq [START_REF] Davenport | An approach to human comfort criteria for environmental wind conditions[END_REF] where = + τ ν 0.5 3 is the dimensionless relaxation time, ν is the kine- matic viscosity, f α neq is the non equilibrium distribution function and f α eq is the local equilibrium distribution function given by Table 3 Different criteria for wind effect on human activities.

+ U σ 3.5 H U 1 m s 10 / 0.022% Isyumov [8] U H 1 m s 5.45 / 1.5% Soligo [48] U H 1 m s 2.5 / 20% B (Standing) NEN 8100 [45] U H 1 m s 5 / 5% Lawson [46] U H 1 m s 5.45 / 2% U s 3 m s 9.3 / 2% Melbourne [11] + U σ 3.5 H U 1 m s 13 / 0.022% Isyumov [8] U H 1 m s 7.95 / 1.5% Soligo [48] U H 1 m s 3.9 / 20% C (Walking) NEN 8100 [45] U H 1 m s 5 / 10% Lawson [46] U H 1 m s 7.95 / 2% U s 3 m s 13.6 / 2% Melbourne [11] + U σ 3.5 H U 1 m s 16 / 0.022% Isyumov [8] U H 1 m s 10.75 / 1.5% Soligo [48] U H 1 m s 5 / 20%
Fig. 6. Localization of the East and West points.

= ⎡ ⎣ ⎢ + + ⎤ ⎦ ⎥ f ρω c u c c Q u u 1 1 2 α eq α αi i s s αij i j 2 4 (6)
with c s the speed of sound, ω α the weighting coefficients of the lattice and

= - Q c c c δ αij αi αj s ij 2 .
Following the Chapman-Enskog expansion the non equilibrium function can be approximated by Refs. [START_REF] Latt | Lattice Boltzmann method with regularized pre-collision distribution functions[END_REF][START_REF] Latt | Hydrodynamic Limit of Lattice Boltzmann Equations[END_REF].

∑ ⎜ ⎟ ≈ - ⎛ ⎝ ∂ ∂ + ∂ ∂ ⎞ ⎠ = - f τω c Q ρu x ρu x Q c c c f f 2 2 ( ) α neq α s αij j i i j αij s α αi αj α α eq 2 4 (7)
It is seen from Eq. ( 7) that the non equilibrium function can be computed using either the macroscopic quantities and their space derivatives or the local distribution functions. In the following, it is chosen to denote f α neq FD , the non equilibrium function computed using finite differences and f α neq LBM , the non-equilibrium function computed from the local distribution functions, respectively. Using only f α neq LBM , can lead to numerical instability whereas considering only f α neq FD , can lead to strong numerical dissipation. To avoid these problems a mixed reconstruction is defined as

= + - f σf σ f (1 ) α neq α neq LBM α neq FD , , (8)
with ∈ σ [0,1] and,

∑ = - f Q c c c f f 2 ( ) α neq LB αij s α αi αj α α eq , 4 (9a)
⎜ ⎟ = - ⎛ ⎝ ∂ ∂ + ∂ ∂ ⎞ ⎠ f τω c Q ρu x ρu x 2 α neq FD α s αij j i i j , 2 (9b)
The regularized collision operator used in this paper is obtained combining Eq. ( 5) and Eq. ( 8), yielding

= - + - τ σf σ f Ω 1 [ ( 1 )
] with = σ 0.97, which is observed to be well suited for Implicit Large- Eddy Simulation. Therefore, no explicit subgrid model is added for the simulations presented hereafter.

α α neq LBM α neq FD , , (10)

Validation of the method

Description of the test case

The present lattice Boltzmann model is first validated with simulations of the Shinjuku area in Tokyo corresponding to the case F of the AIJ open source database [START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF][START_REF] Tominaga | Cross comparisons of CFD prediction for wind environment at pedestrian level around buildings. comparison of results for flowfield around building complex in actual urban area[END_REF] for the north wind condition. This area is a km 1 2 part of Tokyo with both tall (up to 225 m height) and small buildings. Velocity measurements (Fig. 1) are available for that case at several locations, including both field measurements carried out from December 1975 to November 1983 at the beginning of the development of that area and wind tunnel measurements carried out based on the 1977 field measurements conditions. These data are available online (https://www.aij.or.jp/jpn/publish/cfdguide/index_e. htm) or can be found in Refs. [START_REF] Yoshie | Cooperative project for CFD prediction of pedestrian wind environment in the Architectural Institute of Japan[END_REF][START_REF] Tominaga | Cross comparisons of CFD prediction for wind environment at pedestrian level around buildings. comparison of results for flowfield around building complex in actual urban area[END_REF].

It must be noticed here that no CAD data of the 1977 Shinjuku area were available so it was made thanks to photographs and old drawings as mentioned by Yoshie et al. [START_REF] Yoshie | Cooperative project for CFD prediction of pedestrian wind environment in the Architectural Institute of Japan[END_REF] and Tominaga et al. [START_REF] Tominaga | Cross comparisons of CFD prediction for wind environment at pedestrian level around buildings. comparison of results for flowfield around building complex in actual urban area[END_REF]. Therefore differences can exist between the real urban area at the field measurement time and the geometry used for the CFD computations.

Computational settings

As the Shinjuku area is about

× km km 1 1
with a height of = H 225 meters for the highest building (Shinjuku Mitsui Building, point D in Fig. 1) a computational domain of × × m m m 4600 5000 1500 was de- signed. This domain ensures an upstream length of H 4.5 , a downstream length of H 13 and a height of H 6.5 with lateral boundaries placed H 8 from the Shinjuku area in agreement with Franke et al. [START_REF] Franke | Best practice guideline for the CFD simulation of flows in the urban environment[END_REF] and Tominaga et al. [START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF] best practice guidelines. As explained in Section 1 the computational grid is composed of several embedded volumes with uniform mesh. The computations were performed considering 3 different finest grid resolutions summarized in Table 2. A view of the coarsest mesh is shown in Fig. 2 At the inflow, the mean velocity profile is reconstructed from the values given in the AIJ data base, leading to = U z z ( ) 0.895 0.276 . The lateral and top boundary conditions were defined as frictionless walls following the best practice guidelines [START_REF] Franke | Best practice guideline for the CFD simulation of flows in the urban environment[END_REF][START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF] and a constant static pressure condition was prescribed at the outlet. Standard wall law [START_REF] Afzal | Turbulent flow in a machine honed rough pipe for large Reynolds numbers: general roughness scaling laws[END_REF] was used for the floor and the different buildings walls.

Comparison with literature results

For all of the 3 designed grid 1 h of physical time was simulated to install the flow then another one hour time was used to collect samples and compute the flow statistics. The velocity at the different measurements points was normalized by the value of the velocity at the top of the Shinjuku Mitsui building (point D in Fig. 1). We notice that the grid size near the top of that building was the same as the grid size on the floor for all cases to ensure the same level of accuracy for all the measurements points.

We observe from Fig. 3 that the coarse grid do not allow a good recovery of experimental results close to the floor for a majority of measurement points. Grid refinement lead to a clear improvement of the prediction of normalized velocity and the normalized velocities computed from fine grid are in good agreement with both the field and wind tunnel measurements. Grid sensitivity analysis was carried out on the previous measurement points and on other points in Shinjuku area between the coarse and basic grid and between the basic and fine grid. Good agreement between basic and fine grid (Fig. 4) is observed mainly for normalized velocity higher than 0.35 and basic grid could be a good tradeoff between accuracy and computational cost. However in the next part of this work we will address pedestrian wind comfort at a height of 2 m from the ground and the fine grid seems to be a better choice for that kind of study to ensure accuracy of the results close to the floor.

Application to pedestrian wind comfort

Computational settings

The domain of study (Fig. 5) for wind comfort assessment is based . The use of the finest grid is in agreement with Tominaga et al. [START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF] best practices guidelines as we study pedestrian wind comfort on an area higher than the third grid point from the floor.

Wind comfort criteria

Different kind of criteria can be found in the literature, some of them are based on the direct effect of wind on people activities [START_REF] Hunt | The effects of wind on people; new criteria based on wind tunnel experiments[END_REF][START_REF] Murakami | New criteria for wind effects on pedestrians[END_REF] and are related to threshold values for the local effective velocity U e from which the realization of these activities become more difficult as summarized in Table 3.

Others criteria are based on two different parameters, namely the definition of the local effective velocity U e as previously, along with the estimation of the probability P that this local effective velocity exceed a threshold value U THR . This probability function is given by Eq. ( 11) [START_REF] Blocken | Modification of pedestrian wind comfort in the Silvertop Tower passages by an automatic control system[END_REF][START_REF] Sanz-Andres | Pedestrian wind comfort: feasibility study of criteria homogenisation[END_REF].

⎜ ⎟ > = ⋅ ⎡ ⎣ ⎢ -⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ P U x U A θ U G c θ ( ( ) ) 100 ( )exp ( ) θ T H R THR θ x k θ ( ) (11)
where A θ ( ), c θ ( ) and k θ ( ) are the Weibull parameters of the direction θ. Sanz-Andres and Cuerva [START_REF] Sanz-Andres | Pedestrian wind comfort: feasibility study of criteria homogenisation[END_REF] define the amplification factor G θ x at the point x influenced by a building by the ratio between the velocity U x ( ) at the point x and the velocity at the reference height of the meteorological station U ref .

= G U x U ( ) θ x ref ( 12
)
The definition of the effective velocity depends on the authors. The relation = + U U ασ e U is often used with different values of the parameter α: Melbourne [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] used a value = α 3.5, He and Song [START_REF] He | Evaluation of pedestrian winds in urban area by numerical approach[END_REF] used = α 3 and averages are computed over 10 min whereas Murakami and Deguchi [START_REF] Murakami | New criteria for wind effects on pedestrians[END_REF] use = α 0 with average values computed on 3 s. Lawson [START_REF] Lawson | The wind content of the built environment[END_REF] and Isyumov and Davenport [START_REF] Isyumov | The ground level wind environment in built-up areas[END_REF] define the effective velocity following the Beaufort scale which permit to adapt the method to the average period of the velocity. In the next part of this paper we will consider three different wind comfort criteria, Lawson [START_REF] Lawson | The wind content of the built environment[END_REF] criteria using velocity average over 3 s and 1 h, Melbourne [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] using the average velocity and the standard deviation of the velocity and the Dutch standard [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF] criteria using the mean velocity. The different pedestrian comfort criteria are summarized in Table 4 and classified following the Dutch standard [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF] quality class. Class A means that pedestrian comfort is good for long sitting, class B means that the comfort is good for short exposure and class C means that the pedestrian comfort is good enough for walking. It is worth noting here that Melbourne [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF], Isyumov and Davenport [START_REF] Isyumov | The ground level wind environment in built-up areas[END_REF] and Lawson [START_REF] Lawson | The wind content of the built environment[END_REF] defined quality classes with different comments and sometimes other kind of activities than the Dutch standard [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF] used here but the main ideas they use to create their classification are the same. In the author opinion all the different criteria presented in Table 3 and Table 4 are complementary, considering that the effect of wind is different following peoples. Children, adults or elderly peoples do not have the same definition of wind comfort so it seems important to use different criteria to assess pedestrian wind comfort in urban areas.

Wind comfort assessment

The computations were performed during 1 h of physical time to get a physically established unsteady turbulent flow and an additional 1 h of physical time is computed to collect samples needed to evaluate the various statistical quantities required for that study. During that time the hourly averaged velocity U , the standard deviation of the velocity σ U and the 1200 samples of 3 s averaged velocities were stored at every grid point of the area, then the local maximum value of the 3 s average velocity U max s 3 was extracted at each point. Fig. 7 and 8 display the time history of the velocity fluctuations for 1 h of physical time at two locations, one at the east side and one at the west side of the central building (Fig. 6). These two locations were chosen because they are located in areas with similar average velocity level ( = U ms 2.65 / for west point and = U ms 2.54 / at east point) but with very different values for maximal velocity and standard deviation of the velocity (Fig. 10). This means that different physical behaviors are observed on these two locations, mainly on the amplitudes of the velocity fluctuations which are more important on the east point (Fig. 8) than on the west point (Fig. 7).

It is also observed that the sampling time is long enough to get statistically converged quantities, since it is much larger than the observed energetic low frequencies. Fig. 9 shows the velocity histograms for these two points by step of m s 0.1 / and the associated Gaussian distributions. In both cases, the main range is located around m s 2.7 / (between 2.7 and 2.8 for west point and between 2.6 and 2.7 for east point) but the velocity distributions exhibit some significant differences, which indicate that the flow physics is not the same at the two locations. At the east point, which is in the wake of a tall building, one can observe that the histogram is flatter than at the west point, due to the strong coherent vortical events advected in the wake of the upstream building. Therefore, the flatness of the velocity pdf (probability density function) is not the same at the two test points, for the west point the kurtosis excess (

= - - k U U σ ( )/ 3 U 4 4
) is negative ( = -k 0.326) and positive for the east point ( = k 0.123). It is also ob- served that the two histograms are not symmetric with respect to their peak considering that the skewness (

= - s U U σ ( )/ U 3 3
) is negative at these two locations ( = -s 0.179 for west point and = - s 0.125 for east point). Another important observation is that the velocity pdf is not Gaussian at the two locations since the skewness and kurtosis excess of Gaussian pdf are equal to 0, showing that simplified wind models based on a Gaussian velocity model may lead to unsatisfactory results. Pedestrian wind comfort is linked to wind velocity magnitude but also to spatial and temporal velocity fluctuations. Wind speed variability is illustrated in Fig. 10 which displays the 1-hour averaged velocity field, an example of 3-seconds averaged velocity field and the 1hour standard deviation of velocity at 2m from the ground. A high spatial variability is observed on the three maps, as expected in dense urban areas in which many building wakes interact and/or impinge on other buildings and that channeling phenomenon is involved between close buildings. As expected, the 3-second averaging yields a less smooth field than the 1-hour averaging, with much wider high-velocity domains. The east test point is observed to be located in a high standard deviation/variability region (as large as 30% of the mean velocity) influenced by the wakes of three buildings and a strong channeling mechanism. Such a region can therefore be considered as not comfortable for pedestrians due to high velocity fluctuations following Murakami and Deguchi [START_REF] Murakami | New criteria for wind effects on pedestrians[END_REF] or Hunt et al. [START_REF] Hunt | The effects of wind on people; new criteria based on wind tunnel experiments[END_REF].

Murakami and Deguchi [START_REF] Murakami | New criteria for wind effects on pedestrians[END_REF] have shown that 3-seconds averaged velocity higher than m s 10 / can lead to serious effects on pedestrians or very serious effects when it exceeds m s 15 / . From Fig. 10 we can see that such velocity levels do not occur in the simulation with = U ms 1.69 / ref but Fig. 11 shows that increasing the inflow velocity can lead to exceed these threshold velocities. From these results we can say that depending on the inflow velocity it can be very difficult to walk on this area The Dutch standard criteria [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF] differs from the others ones [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] [46] since it is based on the probability of exceedance of the threshold velocity of m s 5 / whereas the others criteria define several velocities with the same maximal probability exceedance to define different pedestrian wind comfort quality class. Some of the classes defined in Table 4 are illustrated in Figs. [START_REF] Franke | Best practice guideline for the CFD simulation of flows in the urban environment[END_REF][START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF][START_REF] Blocken | Computational Fluid Dynamics for urban physics: importance, scales, possibilities, limitations and ten tips and tricks towards accurate and reliable simulations[END_REF] following the Dutch standard [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF], Melbourne [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] and Lawson [START_REF] Lawson | The wind content of the built environment[END_REF] criteria. As expected Fig. 13 shows that Lawson criterion computed with U H 1 is more restrictive for long sitting class than the Dutch standard since the threshold velocity and the maximum probability of exceedance are lower.

This figure also shows that building wakes are mostly considered as acceptable for long sitting exposure following the Dutch standard [START_REF] Nen | Wind comfort and wind danger in the built environment[END_REF] because of their low average velocity whereas they are not following Melbourne criteria [START_REF] Melbourne | Criteria for environmental wind conditions[END_REF] because of the important unsteady effects. It seems from Figs. 13-15 that lowering the averaging period on Lawson criterion defined with the Beaufort scale makes it more restrictive: class A area is reduced with computation from U s 3 and class C area is increased with U s 3 -based criterion. This comes from the definition of the effective velocity: using 1 h averaged velocity we account for steady effects of the wind probably more important for long exposure whereas using 3-seconds averaged velocity we take into account wind gust which can affect walking peoples. Fig. 15d) should not be analyzed in the same way since the main part of the study area is not suitable for standing but it shows that wind gust can affect the main part of this area and that pedestrian balance could be momentarily affected. From those figures (13)-( 15) Melbourne criteria taking account the local mean flow and the local turbulent intensity seems more restrictive than the Dutch standard, Ohba et al. [START_REF] Ohba | Study on the assessment of environmental wind conditions at ground level in a built-up area -based on long-term measurements using portable 3-cup anemometers[END_REF] while comparing several criteria concluded his assessment results in a severe evaluation of pedestrian comfort compared to Murakami et al. [START_REF] Murakami | New criteria for wind effects on pedestrians[END_REF] and Isyumov and Davenport [START_REF] Isyumov | The ground level wind environment in built-up areas[END_REF] criteria. Fig. 16 shows locations that do not belong to classes A, B or C, i.e. locations where pedestrian comfort is not acceptable for sitting or normal walking according to the different criteria. It seems that such locations are almost not present in the test domain according to Dutch standard whereas few locations, mostly in the wake of the top right building, are found following the other criteria. These locations are subjected to more important unsteady effects (see Fig. 3) than the others so it's not surprising to have few differences in the results of the different criteria depending on the definition of their effective velocity or the value of the maximum probability of the threshold velocity. From these results it seems that the use of different criteria for wind comfort assessment is interesting considering they are taking into account different parameter such as mean velocity, turbulent intensity or wind gust. These parameters are important parameters to define pedestrian wind comfort criteria.

The differences between the existing criteria is further emphasized looking at Table 5, which displays the percentage = R N N 100 / i (with N i the number of grid points in class i and N the total number of grid points in the area) of the total test area for each class obtained with the various criteria. The last line deals with the overlap of the point of the same class obtained with the four criteria. A striking conclusion is that the overlap is almost null, indicating that criteria are not equivalent at all and that they lead to very different conclusions in practical cases (see Fig. 17). The duration of the time-averaging period is also seen to be of crucial importance when using Lawson's criterion. As a matter of fact, a long averaging time results in a strong smoothing of the extrema of the instantaneous signal, yielding the prediction of large areas in zone A (about 1/3 of the total area in the present test case). On the opposite, considering a short time averaging does not eliminate the strong wind gusts, yielding a much more restricted comfort map.

Wind comfort criteria combinations

Since the different criteria proposed in the literature lead to different pedestrian comfort map as shown in Fig. 17, the use of a new criterion combining the previously defined criteria is necessary. In this study we proposed 4 new criteria:

• Criterion 1: for every grid point the less restrictive criterion assess is kept,

• Criterion 2: for every grid point the most restrictive criterion assess is kept,

• Criterion 3: for every grid point the most often predicted criterion is used, in case of a tie between two criteria the less restrictive is kept,

• Criterion 4: for every grid point the most often predicted criterion is used, in case of a tie between two criteria the more restrictive is kept. Fig. 18 shows the comfort maps obtained with these new criteria. The criterion 1 using the less restrictive assessment lead to comfort map quite similar to the one obtain with the Dutch standard criteria. The criterion 2 gives results very close to Lawson criterion that was more restrictive than the other ones. These two criteria are probably not the most representative since they are based on the extreme values assessed by the standard criteria. On the opposite the comfort map obtained using the proposed criteria 3 and 4 based on the most predicted class for each point seems to be a good compromise for wind comfort assessment. The main differences between these both criteria appears between class B and C around the two north buildings but the global aspect of the comfort map is very similar for these two criteria. In authors opinion, the use of CFD velocity fields combined with the proposed criteria 3 or 4 estimated from many different commonly used criteria seems an interesting way for wind comfort assessment since it will lead to more representative comfort map by removing extreme values.

Conclusions

Several LBM-based LES of urban flow in a realistic full scale case have been carried out, with application to wind comfort. A first observation is that LES provides a deep insight into flow physics and gives access to a huge amount of physically relevant and accurate data that can be used as input for wind comfort criteria. An observation is that wind velocity fluctuations, which originate in the combination of many physical phenomena (wakes, boundary layers, channeling …) are not Gaussian fluctuations. Using LES unsteady data, it is observed that wind comfort criteria that include wind unsteadiness as an input parameter can lead to conclusions that are very different from those based on mean flow velocity amplitude. This illustrates the fact that unsteady CFD approaches such as LES or DES should be privileged for wind comfort studies. Another observation is that the various existing criteria lead to significantly different results on the same case which complicates their practical use. In this study several methods have been proposed to combine the different criteria available in the literature. In order to get a robust conclusion in practical applications the use of a large amount of these classical wind comfort criteria is recommended to keep the more often predicted criteria.

Fig. 1 .

 1 Fig. 1. Positions of measurements points.

Fig. 2 .

 2 Fig. 2. Grid visualization for the coarse case.

Fig. 3 .

 3 Fig.3. Comparisons of the present simulations with measurement results: Field measurements * with standard deviation, wind tunnel measurements , coarse grid , basic grid and fine grid .

Fig. 4 .

 4 Fig. 4. Grid sensitivity analysis, • represent values on previous measurement points and * represent values on other points in the domain.

Fig. 5 .

 5 Fig. 5. Shinjuku area with new buildings (blue) and the area of interest (red), the arrow indicates the inflow direction. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Visualization of the velocity fluctuations on the West point in the case = U ms 1.69 / ref

Δ 32 .Δ 4 )Δ 8 )

 3248 . The coarse grid is composed of several layers of refinement levels close to the ground, = x m Δ 2 up to 70 m high (35 grid points) then x Δ is increased (by a ratio of 2) by layers of 12 points up to the maximal value = x m The areas around the tall buildings were refined ( = x m Δ 2 close to buildings) using individual refinement volumes to ensure that the different wakes will be correctly captured. Additional refinement levels were added around tall buildings, one volume ( = x m including the different individual areas and another one ( = x m surrounding the previously defined volume. The basic grid is based on the coarse one with an additional layer up to 40 m from the ground with = x m Δ 1 and the same procedure was applied from basic to fine grid with a layer of 20 m ( = x m Δ 0.5 ) from the ground. For all the 3 tested grids the measurements points of Fig. 1 are located in the finest resolution level.

Fig. 10 .

 10 Fig. 10. Visualization of a) 1 h averaged velocity, b) 3 s averaged velocity and c) 1 h standard deviation of the velocity computed with = U ms 1.69 / ref . The black points indicate the position of previous east and west points.

Fig. 11 .

 11 Fig. 11. Visualization of the maximum 3 s average velocity (m/s) for = U ms 5.07 / ref

Fig. 12 .

 12 Fig. 12. Visualization of the different quality class assessed with Dutch standard criteria using different reference values of the inflow velocity.

  on the Shinjuku area previously described. An area of interest of × m m 400 400 ( ∈ -x [ 70,330] and ∈ - y [ 130,270]) was defined and a building of x [100,160] and ∈ y[START_REF] Gousseau | Near-field pollutant dispersion in an actual urban area: analysis of the mass transport mechanism by highresolution Large Eddy Simulations[END_REF]100]) was added in that area. Since 1977 lot of new buildings were built in the Shinjuku area and our added building is quite similar to the the Shinjuku Center Building built in 1979. We notice that another building was added in the geometry ( another study coupled with this one to reduced the computational cost. The area of interest is a uniform grid with = x m Δ 0.5 containing 466,375 points 2 m far from the floor. The boundary conditions are the same as previously and 3 inflow velocity profiles were used, the velocity profile used in Sec. 3, 3 times and 6 times this profile corresponding to =

Fig. 13 .

 13 Fig. 13. Visualization of the class A (in blue) with a) NEN8100, b) Melbourne, c) Lawson with 1 h average velocity and d) Lawson with 3 s average velocity criteria. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 14 .

 14 Fig. 14. Visualization of the class B (in green) with a) NEN8100, b) Melbourne, c) Lawson with 1 h average velocity and d) Lawson with 3 s average velocity criteria. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 15 .

 15 Fig. 15. Visualization of the class C (in yellow) with a) NEN8100, b) Melbourne, c) Lawson with 1 h average velocity and d) Lawson with 3 s average velocity criteria. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 16 .

 16 Fig. 16. Visualization (in red) of the locations not belonging to class A, B or C with a) NEN8100, b) Melbourne, c) Lawson with 1 h average velocity and d) Lawson with 3 s average velocity criteria. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 17 .

 17 Fig. 17. Visualization of the locations belonging to a) Class A (blue), b) Class B (green), c) Class C (yellow) according to the four criteria and d) where the four criteria do not agree (black). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 18 .

 18 Fig. 18. Visualization of the wind comfort quality class assessed following the four criteria combinations proposed in this paper, a) Criterion 1, b) Criterion 2, c) Criterion 3 and d) Criterion 4. Quality class A is in blue, class B is in green, class C is in yellow and the other classes are in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1

 1 Example of existing CFD simulations at pedestrian level.

		Turbulence model	Computational domain	Area of interest	Pedestrian level	Wind comfort criteria	Grid points
	He [35]	LES	× 700 700 2 m			m 1.7
	Blocken [20]	RANS	× 900 700 2 m			m 1.75	Bottema [3]	× 2.9 10 6
	Blocken [21]	RANS	× 3000 3000 2 m	Amsterdam Arena	m 2	NEN8100 [45]	× 2.8 10 6
	Letzel [42]	LES	km 1.6 2	≃ km 1 2	m 2.5		× 7.2 10 8 &	× 1.6 10 9
	Jansen [23]	RANS	× 2077 1838 2 m	× 1918 1430 2 m	m 1.75	NEN8100 [45],	× 7.5 10 6
								Isyumov [8],
								Melbourne [11],
								Lawson [46]
	Montazeri [26]	RANS	× 2076 1963 2 m	Antwerp tower	m 1.7	NEN8100 [45]	≃	× 16 10 6
					balconies
	Blocken [22]	RANS	≃	× 2700 2300 2 m	× 1600 1100 2 m	m 1.75	NEN8100 [45]	× 7.5 10 6
	Shi [47]	RANS	× 3000 3000 2 m			m 1.5		× 1.14 10 6
	Zheng [29]	RANS	× 7950 7650 2 m	Outdoor platforms of	m 2	NEN8100 [45],	× 7.18 10 6
					megatall building		Lawson [46]
	Kang [24]	RANS	× 1000 1000 2 m			m 1.75	Isyumov [8]	≃	× 14.2 10 6
	Adamek [43]	LES	× 4270 2440 2 m	≃	× 600 600 2 m	m 1.5	Soligo [48]	× 1.1 10 6
	Present study	LES	× 4600 5000 2 m	Shinjuku area	m 2		× 22 10	6	-	× 136 10	6

Table 2

 2 Grid parameters and computational time for the Shinjuku area test case.

	Grid	x Δ min (m)	t Δ min (s)	Grid refinement level	Grid points (10 6 )	Number of processors	Computational time for 1 h
	Coarse	2	0.03	5	22	120	9h
	Basic	1	0.015	6	54	240	20h
	Fine	0.5	0.0075	7	136	504	50h

Table 4

 4 Different wind comfort criteria.

	Class	Criteria	U e	U THR	P max
	A (Sitting)	NEN 8100 [45]	U H 1	m s 5 /	2.5%
		Lawson [46]	U H 1	m s 3.35 /	2%
			U s 3	m s 5.7 /	2%
		Melbourne [11]

Table 5

 5 Fraction (in %) of points of the Shinjuku area belonging to the different quality class according to several criteria and overlap between them.

	Criteria	Zone A	Zone B	Zone C	Other zone
	NEN8100 [45]	57.5%	29%	13.5%	0%
	Melbourne [11]	41.5%	30.8%	20.8%	6.9%
	Lawson [46] (U H 1 )	34%	24.6%	33.7%	7.7%
	Lawson [46] (U s 3 )	13.8%	29.7%	51.1%	5.4%
	According to all criteria	13.6%	0.14%	2.9%	0%

Acknowledgments

This work was supported by the French project CLIMB, with the financial support of BPIFrance (Project No. P3543-24000), in the framework of the program "Investissement d'Avenir: Calcul Intensif et Simulation Numérique". This work was performed using HPC resources from GENCI-TGCC/CINES (Grant 2017-A0012A07679).