Atmospheric boundary layer over urban roughness: Validation of large-eddy simulation Open Access

Cities have long been regarded as centers of economic and demographic growth, resulting in a significant concentration of population in their cores. According to the United Nations, by 2030, two-thirds of the global population will reside in urban areas. Thus, the study and understanding of urban flows to improve forecasting and develop accurate prediction methods is the focus of urban sustainability. Urban climate studies can be conducted at different scales, including the meteorological mesoscale, the meteorological microscale (up to about 2 km), the building scale (up to a few hundred meters), the indoor environment (⁠ $≈10 m$⁠), etc.

A large number of these studies have been conducted in generic or scaled-down geometries, from two-dimensional (2D) street canyons to 3D generic geometries, in order to understand the fundamental fluid dynamics associated with urban flows. In these studies, either Reynolds-averaged Navier–Stokes (RANS) equations or large-eddy simulations (LES) have been used.

Studies of flows over simplified 2D and 3D urban-like geometries have provided insight into the complex flow physics within street canyons due to their similarity to the topology of urban environments.^1–5 Among these studies, the complex flow pattern, turbulence statistics, dispersion of passive scalars, and wind loading are of particular interest.

Cheng and Castro⁶ considered a boundary layer over an urban-like roughness in multiple configurations. Both staggered and in-line matrices of square prisms with identical heights were investigated in the experiment. In addition, the study also included measurements of a staggered matrix comprising prisms of random heights. Prisms were equally spaced throughout the investigations where the Reynolds number considered ranged from $ReH=0.5×104$ to $1.2×104$ (⁠ $ReH$ here is based upon the prism height, H, and freestream velocity, $U∞$⁠). The study focused on the characteristics of the inertial and roughness sublayers. The roughness sublayer was found to have a depth of $1.80–1.85H$ for prism matrices of identical H. Dispersive stresses arising from the inhomogeneity were negligibly small as compared to spatially averaged Reynolds shear stresses. On the other hand, the thickness of the roughness sublayer increased substantially for the case with random heights. Castro et al.⁷ extended the investigations and focused on the turbulence statistics in the roughness layer for a staggered matrix of an identical height. A two-point correlation revealed the dominant scales of turbulence were of the same order as the H of the prisms. It also suggested the existence of a scale much smaller near the top of the canopy region. The results of stress anisotropy agreed with the previous findings in that the roughness increased the level of isotropy. This effect was most pronounced within the canopy layer, where sweep events dominated the momentum transport.

Xie and Castro³ performed LES and studied the turbulent flow over staggered wall-mounted cubes with both fixed and random heights. Reynolds numbers $ReH$ between $5×103$ and $5×106$ were considered. Numerical results were extensively validated against data from both direct numerical simulations and experiments, indicating the good predictive capability of LES for such flow scenarios. This partially flows over urban-like obstacles exhibit a wide inertial subrange, suggesting that turbulence reaches a quasi-isotropic state at a relatively lower frequency than non-vortex-shedding flows at similar Reynolds numbers. Moreover, steady RANS were also solved and proved inadequate, especially in the canopy region, due to the inherent unsteadiness of the flow.

Inagaki and Kanda⁸ conducted field experiments using scaled geometry to investigate the turbulence similarity over cubical obstacles. The roughness Reynolds number ranged between $104$ and $105$⁠. The study examined turbulence similarity within the inertial sublayer using inner-layer scaling variables. It was concluded that inner-layer scaling effectively described the wall-normal velocity fluctuations and Reynolds stresses, regardless of surface roughness or outer-layer conditions, under near-neutral stratification in atmospheric flow. However, inner-layer scaling did not hold for horizontal velocity fluctuations due to the influence of outer-layer disturbances.

While significant progress has been made using simplified geometries—providing insights into airflow and turbulence within cities, as well as improving the representation of drag forces and turbulent mixing—certain limitations remain, in particular when using scaled-down geometries. These simplifications often fail to capture the full complexity of real urban environments.⁹ In fact, scaled geometries may overlook the heterogeneity of building shapes, heights, and arrangements, leading to inaccuracies in representing wind flow patterns. Therefore, there is a need for improved parametrization of urban canopy turbulence that incorporates greater geometric detail.¹⁰ 

Studies of urban flows over semi-idealized or realistic urban geometries have been mostly performed using either RANS or unsteady RANS models. This preference is mainly attributed to the complexity of the flow, which spans a wide range of spatial and temporal scales, combined with the large dimensions of the domain under investigation, making scale-resolving simulations computationally challenging. However, RANS models are inherently limited in their ability to capture the unsteady features of turbulent flows within urban canopy layers, such as vortex shedding, separations, and reattachment. Consequently, the majority of such studies often suffer from reduced accuracy, particularly in narrow street canyons and downwind regions of buildings, where flow accelerations are more pronounced.^11,12

High-fidelity simulations have been proven to be more accurate in resolving the flow features in urban environments.^13–15 Xie and Castro¹⁵ studied the flow and dispersion characteristics in a 1:200 urban geometry corresponding to central London; they also concluded that a resolution of $1 m$ at full-scale size would be appropriate to describe the mean turbulent flow behavior. Other examples of this kind of study on scaled-down urban areas include Gousseau et al.¹⁶ and Hassan et al.¹⁷ While the former conducted RANS and LES simulations in a 1:200 geometric model to study pollutant dispersion in a neighborhood of Montreal, comparing their results against a wind tunnel prototype, the latter performed LES to examine traffic exhaust dispersion in a 1:150 scaled design of a separated street canyon. The critical role of LES in capturing unsteady flow features, higher-order turbulence statistics, and complex flow patterns was highlighted in the study of Tolias et al.,¹⁸ who also emphasized the importance of conducting a grid resolution study to achieve accurate results.

Kirkil and Lin¹⁹ conducted an LES over a 1.10 × 0.76 km² region of downtown Oklahoma City to analyze the effect of the building height on the flow configurations. Similarities to flow fields around cylinders or inside cavities were observed. Auvinen et al.²⁰ used LES to study the effect of several modeling parameters on the predictability of urban flows under realistic conditions, emphasizing the importance of mesh resolution at the pedestrian level and the height of the domain when it comes to the simulation of the transition of a fully developed atmospheric boundary layer (ABL) into an urban boundary layer. More recently, the impacts of urban geometry heterogeneity on flow statistics were studied by Cheng and Yang,²¹ who conducted simulations on realistic urban geometries covering an area of about 1.0 × 0.5 km². The study demonstrated the importance of turbulence statistics to improve the prediction of current urban canopy models.

Oh et al.²² performed an LES using seasonal atmospheric data, with a 3 m resolution at the pedestrian level, to assess wind and thermal comfort in urban environments and evaluate the impacts of high-rise buildings. The study, which considered a canopy-layer-based Reynolds number of $5.4×104$⁠, focused on the thermal wind environment around the Yonsei University campus, featuring two high-rise buildings. The results revealed that downdraft winds and large-scale circulations developed around realistic high-rise structures. These induced flows enhanced building-level wind speeds via the Venturi effect and reduced pedestrian-level temperatures by drawing cooler air downward.

A recent study by Shaukat and Giljarhus²³ has demonstrated that commonly used validation metrics for evaluating simulation quality should not be the sole criteria for assessment. Similarly, research by Tian et al.²⁴ has once again highlighted the necessity of using full-scale geometries to effectively analyze and capture the complexities of urban flows. Together, these findings emphasize the importance of conducting high-accuracy simulations on real-world urban geometries.

In a review of the published literature, Toparlar et al.²⁵ reported that less than 3% of investigations employed LES for modeling the urban microclimate in full-scale scenarios. With the increase in computational capacity, LES is increasingly used to simulate various phenomena such as wind patterns, pollutant dispersion,^26,27 and thermal effects²⁷ in full-scale urban environments in the presence of the atmospheric boundary layer (ABL).

In the present study, wall-modeled LES of the ABL over an urban roughness is implemented. Two geometries are considered, namely, a wall-mounted array of cubic prisms and the “Michel-Stadt” model. Simulations of both cases play an important role in the validation and understanding of the microclimate in full-scale urban areas for their resemblance to flow over complex terrain and in the vicinity of full-scale buildings. In most of the aforementioned LES studies, the inflow boundary conditions implemented fall into one of three categories: periodic boundary condition in the streamwise direction,^3,14,26 prescribed velocity profile at inflow,^13,17,19,27 or artificially generated turbulent inflow.^15,16,18,22 In the current study, a realistic turbulent ABL is implemented at the inlet via an online precursor simulation. This approach ensures a fully developed ABL with stochastic fluctuations that accurately represent the turbulent approaching flow.

The main objectives of the present work are to assess the performance of the methodology used to model the turbulent flow in complex urban scenarios, as well as to demonstrate the potential of this type of simulation for accurately modeling urban flows. A thorough validation, including a direct comparison of turbulence statistics with wind tunnel data, point-to-point analysis of various quantities, and evaluation via figures of merit, as well as energy spectra, gives credibility to the current setup. In addition, visualization of instantaneous flow structures is presented, providing insight into the flow evolution. In the following, we first present the definition of the study cases in Sec. II and the numerical methodology in Sec. III. The results are compared with wind tunnel data in Sec. IV. The main findings and the accuracy of the numerical setup are discussed in Sec. V, and conclusions are provided in Sec. VI.

In order to demonstrate the potential of using high-fidelity simulations to study turbulent flow in full-scale urban geometries, two different cases with available wind tunnel measurements have been considered. The selected cases are a 3D cubic prism array^28,29—an arrangement of idealized building blocks—and the “Michel-Stadt” BL3-3 case, which represents urban areas in central European cities.³⁰ 

The first case studied is the wind tunnel experiment conducted by Brown and Lawson,²⁹ which also serves to validate the numerical framework. This experiment investigated a neutral ABL approaching a 1:200 scaled array of 3D buildings arranged in seven rows. Each building has dimensions $L=W=H=0.15 m$⁠, where L and W represent the building's length and width. The streamwise and spanwise distances between blocks were set to H.

To replicate the ABL in the wind tunnel experiments, spires and floor roughness elements were used to create an equivalent ABL, characterized by a reference velocity $Uref=3 m/s$ at the reference height, $Href$ (⁠ $Href=H$⁠) and a friction velocity $u*=0.24 m/s$⁠. In the experiments, vertical profiles of streamwise and spanwise velocities, together with turbulent kinetic energy (TKE), were measured at various streamwise positions along the centerline of the array. Although the $ReH$ considered is larger than that in the reference experiment by approximately two orders of magnitude, both values are well above the critical threshold of $ReH≈1.1×104$⁠, beyond which the flow becomes effectively independent of Reynolds number effects.^31,32 Therefore, any interference due to $ReH$ that could potentially compromise the flow similarity is highly limited.

The second case considered is part of the CEDVAL-LES database—a collection of datasets designed for validating LES.³⁰ The experiments were conducted at the wind tunnel facility of the Meteorological Institute at the University of Hamburg. The specific urban geometry under study, referred to as the “Michel-Stadt” case (reference BL3-3), represents a semi-idealized urban layout typical of residential areas in central European cities. The scaled model, at 1:225, consists of 60 flat-roof building blocks with courtyards and roof heights of $0.625H (15 m)$⁠, $0.75H (18 m)$⁠, and $H (24 m)$⁠, covering a total area of $1320×830 m2$ at full scale.

In the experiments, the approaching ABL was modeled using roughness elements to replicate a highly rough flow environment. This setup was characterized by a surface roughness length $z0=1.53 m$ and a friction velocity $u*=0.596 m/s$⁠. The reference velocity was defined as $Uref=6.1 m/s$ at a reference height of $Href=100 m$⁠.

Laser Doppler anemometry (LDA) measurements were used to obtain time-series data for two velocity components (streamwise and lateral) across 40 vertical profiles and 5 horizontal levels within the central part of the city. These levels correspond to heights of 2, 9, 18, 27, and 30 m, as described by Hertwig et al.³³ Additionally, first- and second-order flow statistics are available for comparison. For the remainder of the paper, the aforementioned cases will be referred to as cases 1 and 2, respectively.

The filtered incompressible Navier–Stokes equations, Eqs. (1) and (2), are solved using SOD2D (Spectral high-Order coDe 2 solve partial Differential equations).³⁵ It is an open-source code based upon a spectral element version of Galerkin's finite element method continuous model.³⁶ Projection stabilization precludes numerical oscillations due to dominant convection while introducing very low numerical dissipation. The aliasing effects of the reduced order integration caused by employing spectral element method integration for convective terms are countered with a skew-symmetric splitting.³⁷ For temporal discretization, a BDF-EXT3 high-order operator splitting approach is used to solve the velocity–pressure coupling, allowing for equal order interpolation of velocity and pressure.³⁸ Regarding spatial discretization, in the present study, the mesh is based on fourth-order hexahedral elements.

SOD2D has been developed to efficiently utilize the growing computational power available globally, with a particular emphasis on leveraging GPUs, which play a central role in modern high-performance computing (HPC). To ensure versatility, the code is designed to run on both GPU and CPU architectures. Specifically, SOD2D is written in Fortran and employs Message Passing Interface (MPI) for coarse-grained parallelism and OpenACC for fine-grained parallelism, enabling efficient execution on heterogeneous systems. Additionally, the code uses HDF5 for input/output operations, a robust and widely validated library for HPC applications. For further details, readers are referred to Folk et al.³⁹ 

Figure 1 shows the numerical setup and computational domains for both cases considered. For case 1, shown in Fig. 1(a), the computational domain has a size of $x×y×z=60H×33H×10H$ (where $H=30 m$ denotes the height of the square prism). In accordance with best practice guidelines for simulating urban flows,⁴⁰ the 3D block array is placed 6.7H from the inlet and 10H from the lateral boundaries of the domain. The outlet of the domain is placed 40.3H downstream of the array.

For case 2, as seen in Fig. 1(b), the domain has a size of $x×y×z=159.4H×91.5H×30H$ (H denotes maximum building height $H=24 m$⁠) incorporating a heterogeneous geometry representative of a semi-idealized neighborhood. In the study domain, the geometry is placed along the centerline of the horizontal plane, 20H downstream of the inlet. The lateral boundaries of the domain are placed at a distance of 28.5H from the city, while the outlet of the domain is located at 84.4H downstream of the city.

In order to simulate the flow around the complex urban geometries, unstructured hexahedron-element meshes are constructed. Within each element, the solution is represented by tensor products of fourth-order polynomials in each direction.

For both cases, three levels of refinement have been considered. Details about the grid resolutions used for both cases are given in Fig. 2 and Table I. The computational meshes are designed to resolve an important part of the inertial subrange in the region of interest. To achieve this, three zones of refinement with progressively increasing mesh resolution are implemented. Within the city, the finest meshes are specified with resolutions below $1 m$ (see Z1 in Fig. 2). This results in meshes of $22.8×106$⁠, $62.2×106$⁠, and $207.07×106$ nodes for case 1 and $22.6×106$⁠, $78.11×106$⁠, and $217.84×106$ nodes for case 2. In wall units, i.e., $Δ+=u*/ν$⁠, the grid resolution within the region of interest (zone 1) for case 1 is $Δx+×Δy+×Δzground+=5.38×3.68×5.38 (×103)$ near the ground floor and $Δx+×Δy+×Δzwall+=5.38×3.15×5.38 (×103)$ in proximity to the buildings. For case 2, the respective resolutions are $Δx+×Δy+×Δzground+=4.24×3.01×3.96 (×103)$ and $Δx+×Δy+×Δzwall+=3.51×2.36×3.83 (×103)$⁠, respectively. Despite the fact that the wall-unit values are larger than the conventional ones in wall-modeled LES by approximately an order of magnitude, the grid resolution is proven to be fine enough to accurately capture flow features by the resolved Kolmogorov $−2/3$ scaling in energy spectra (discussed in Sec. V). It is also noteworthy that typical grid resolutions for LES in full-scale geometries reported in the literature are generally on the order of $2.5–5 m$⁠,^41–44 while resolutions down to $1 m$ are rarely reported.^15,45 According to observations of Xie and Castro,¹⁵ a resolution of about $1 m$ is sufficient to achieve reasonable accuracy in the predictions of turbulence statistics in full-scale urban geometries. Thus, in the present work, within the zone of interest—i.e., the urban geometry identified as Z1—grid mesh with resolution below $1 m$ is employed, which is achieved at an affordable cost.

To impose a realistic turbulent ABL inlet condition, a precursor simulation is implemented simultaneously with the downstream simulation of the study domain (online precursor). Typically, two methodologies are commonly used to impose inflow conditions: synthetic data generation (see, for instance, Xie and Castro⁴⁶) or a database storing the spatial and temporal variation of the inflow conditions from a precursor simulation (e.g., Shaukat and Giljarhus²³). While the former may not fully replicate the turbulent characteristics of a real ABL, the latter requires storing a substantial amount of high-frequency data, which should match the time step and the integration time of the main simulation. In this study, we utilize two domains, i.e., the precursor and the main domain, which are solved simultaneously. This approach eliminates the need to store inflow data and seamlessly adapts to the spatial dimensions of the domain and the time step variations of the city simulation. To appropriately set up the inflow conditions, periodic boundary conditions are applied in the x- and y-directions within the precursor domain (see Fig. 1), where the turbulent ABL is driven by an imposed pressure gradient in the x-direction. The outlet of the precursor simulation serves to “drive” the inlet of the study domain: the inflow of the study domain is updated at each time step with data from the outlet of the precursor simulation.

To ensure that the inflow aligns with the conditions documented in the experiments, Fig. 3 compares the numerical results for the fine mesh in both cases with the measurements described in Sec. II. Both profiles of streamwise mean velocity and TKE are plotted. As shown in Fig. 3, the precursor simulation accurately reproduces these statistics in both cases.

In addition, the top and lateral boundaries of the computational domain are modeled as free-slip walls, assuming zero normal gradients for all variables. At the outlet, zero static pressure is imposed. Within urban canyons, non-equilibrium 3D flows are expected around buildings. Consequently, applying an equilibrium wall model on the wall of the buildings (e.g., a log-law near the wall) is inadequate. Therefore, no-slip conditions are applied at the walls, while for the ground, a wall model is employed.

The flow field takes approximately 10 flow-throughs (FT) to pass the transient phase and reach the quasi-steady state. One FT is defined based on the length of the city and the velocity at the edge of the atmospheric boundary layer: $FT=Lcity/Uedge$⁠. This corresponds to 48T and 125T, respectively, for each case considered (T denotes one eddy turnover time, $T=H/u*$⁠). Following the initial transient phase, flow statistics are collected over approximately 40FTs for both cases. This corresponds to 192T and 500T for cases 1 and 2, respectively. For further details regarding the temporal convergence of the flow, please refer to the discussion in Sec. V B.

Finally, it is important to highlight the computational effort it takes for the present simulations. Both cases were executed on Marenostrum V, utilizing eight nodes from the accelerated partition, each equipped with 4 H100 GPUs. For case 1, the setup achieved a performance of 0.128 s per iteration, corresponding to 0.62 ns/dof. Statistical data collection on the finest grid (207 × 10⁶ grid points) over 40FTs required a wall time of 0.4 days, resulting in a computational cost of 404.9 GPU-h. Similarly, for case 2, simulations on the finest grid required a wall time of 6.7 days using 32 GPUs, with a total computational cost of 4744 GPU-h. These results emphasize the efficiency of the current setup in performing accurate and computationally efficient simulations for complex urban environments, achieving a balance between precision and resource utilization.

Extensive comparisons of the results have been performed for both cases. Both profiles of streamwise mean velocity and TKE have been compared to experimental measurements obtained from a wind tunnel. Moreover, for case 2, validation metrics have been used, along with energy spectra obtained at different stations.

Validation of the numerical results begins with grid-convergence studies. For each case, three grid resolutions are employed within the study domain to assess the influence of grid refinement on the numerical results. For case 1, Figs. 4–6 compare profiles of streamwise mean velocity, $u¯/Uref$ (where $·¯$ represents temporal averaging and $Uref$ denotes the reference velocity), wall-normal mean velocity, $w¯/Uref$⁠, and TKE, $12ui′ui′¯/Uref2$ (⁠ $ui′=ui¯−ui$ represents the fluctuating components) with experimental results at selected locations along the plane $y=0$⁠. For clarity, results for the medium grid resolution are omitted. In general, discrepancies between the results obtained using the coarse and fine meshes are minor for both first- and second-order statistics, indicating that the fine mesh resolution is sufficient to achieve good agreement with experimental data. Therefore, including the medium grid results would not alter the conclusions; instead, it would compromise the clarity of the comparison and distract from the observation of key discrepancies.

The $u¯/Uref$ profiles (see Fig. 4) demonstrate good agreement across all stations. Notably, the reverse flow in the recirculation region behind the blocks is reasonably predicted [see Figs. 4(c) and 4(e)]. However, a slight under-prediction of TKE is observed in these regions (Fig. 6), which can likely be attributed to an insufficient resolution close to the ground. Meanwhile, a wall function is not advisable in this region due to the lack of local stress equilibrium in the areas around buildings where the flow exhibits significant separation. The simulations capture the TKE peak above building heights well, corresponding to the shear layer development along rooftops.

Despite the overall good agreement in $u¯/Uref$ profiles, discrepancies are observed in $w¯/Uref$ profiles across all stations (see Fig. 5). They persist even with grid refinement, suggesting that they may stem from experimental uncertainties rather than insufficient resolution. It should be pointed out that reproducing an ABL in a wind tunnel is a challenging task, as it depends not only on the turbulence generation method but also on factors such as Reynolds number mismatch and wall effects. These factors can influence turbulence structures and introduce discrepancies between experimental results and full-scale atmospheric simulations.

For case 2, experimental data include $u¯/Uref$⁠, $v¯/Uref$⁠, and root mean square (rms) of the streamwise fluctuating component, $urms′$⁠, at various locations. For clarity, only the comparison of these statistics in the $x–z$ plane along $y=0$ is presented in Figs. 7–9, while the remaining data are provided in  Appendix A. Similar to case 1, good agreement is observed between the results of different mesh resolutions. $u¯/Uref$ profiles are accurately predicted at all locations; discrepancies are evident in $v¯/Uref$ profiles though. This is particularly the case near building corners and cross-road areas [e.g., Figs. 8(a), 8(d), and 8(f)]. Above rooftops, the flow becomes nearly one-dimensional, and the agreement with experimental data is excellent, reflecting the accuracy of the inflow conditions generated by the upstream precursor simulation.

As in case 1, the simulations accurately capture the $urms′/Uref$ peaks above building heights (see Fig. 9). These peaks are indicative of shear layers formed near the rooftop level of the buildings. However, numerical predictions fail to fully resolve the fluctuations within the reverse-flow regions between buildings [e.g., Figs. 9(e) and 9(f)]. This under-prediction may be attributed to the inherent challenges in modeling separated flows. Overall, the grid-convergence study confirms that the fine meshes with resolution down to 0.75 m are capable of producing reliable results.

Figures 10 and 11 visualize the instantaneous streamwise velocity fields, $u/Uref$⁠, for both cases, providing an insight into the flow features at different planes. For case 1, see Fig. 10(a), a large-scale vortex dominates the spaces between the square prisms, while an extensive wake region develops downstream of the trailing prism. For the current aspect ratio $S/H=2$ (i.e., a ratio of the face-to-face distance between the prisms, S, to the height of a cube), according to Oke,⁴⁷ a skimming flow configuration is expected. That is, the flow above the rooftop bypasses the canyon, and a stable vortex is formed within the canyon. At the pedestrian height $[z/H=0.067 (2 m)]$⁠, see Fig. 10(b), the velocity field is characteristic of the reverse flow in the wake of each square prism. Due to the conservation of mass, the flow accelerates through the aisles as it enters the matrix. This results in the formation of a stable recirculation within the canyon between two consecutive buildings (not shown here), as predicted by Oke.⁴⁷ This phenomenon is well pronounced in the mean-flow field, as shown in Fig. 10(d) where a spanwise symmetry about the center plane (⁠ $y=0$⁠) and a quasi-periodicity in the x-direction within the array of prisms are also observed. As the distance increases beyond the prism height $[z/H=1.25 (37.5 m)]$⁠, see Figs. 10(c) and 10(e), the impact of the square prisms on the flow field is primarily observed in the reduction of velocity magnitude within the matrix, where the reverse flow in the wake region is no longer perceptible.

Figure 11 presents the instantaneous flow structures of case 2. One distinctive feature of a semi-idealized urban configuration is its geometrical complexity, shown in Fig. 11(a). A large-scale vortex occupies compact inter-building spaces, while an extensive wake region forms in larger open areas. Close to the ground, at $z/H=0.083 (2.0 m)$ [see Fig. 11(b)], the flow field within the enclosed building block is mostly dominated by the reverse flow. The wake region downstream of the trailing column of the buildings, see Fig. 11(d), does not appear particularly extensive. This is likely a result of spatial inhomogeneity upstream. As the distance increases beyond the buildings and reaches $z/H=1.25 (30 m)$⁠, the reverse flow is barely observed within and downstream of the high-rise buildings $(z=24 m)$⁠. The impact of the low- and high-rise buildings on the flow field reduces as distance increases away from the ground, as seen in Figs. 11(c) and 11(e). The specific distance at which the urban geometry considered no longer influences the flow field can be evaluated by identifying the roughness sublayer.⁶ This question will be addressed in detail via mean statistics in further studies.

In Sec. IV, we provided a detailed comparison of the results obtained at different locations with those measured in wind tunnel experiments, based on first- and second-order statistics. Although the analysis of these statistics provides critical insights by highlighting areas where predicted values differ from experimental measurements, such as recirculating flows in particular, close to the ground, these discrepancies are inherently challenging to quantify.

In this sense, validation metrics can serve as a complementary quantitative assessment of the high-fidelity simulation results. These metrics, proposed by Franke et al.,⁴⁰ are frequently used in existing studies (for instance, Vranckx et al.,⁴⁸ Yu and Thé,⁴⁹ Tolias et al.,¹⁸ Ding et al.,⁵⁰ etc.) and here they are applied to case 2, “Michel-Stadt” BL3-3, as extensive measurements at more than 1000 locations are available in the database.

In this study, $Pi$ and $Oi$ denote the normalized numerical predictions and experimental values, respectively. The parameter D represents the allowable relative deviation, which accounts for the reproducibility of wind tunnel measurements in built environments and the permissible inaccuracies in simulation results. Based on the recommendations of Schatzmann et al.,^51, D is set to 0.25. The parameter W defines the allowable absolute deviation, reflecting uncertainties caused by interpolation and measurement repeatability. This value is case dependent, with reported ranges in the literature from 0.01 to 0.1.^18,49,52,53 In this study, W is assigned values of 0.0165 for streamwise velocity and 0.0288 for spanwise velocity, following the work of Tolias et al.¹⁸ These values align with the experimental measurements conducted for case 2. The results for these parameters are summarized in Table II.

Table II presents the validation metrics for data at different heights, along with the overall values for both streamwise and spanwise velocity components. Scatter plots comparing high-fidelity simulation results with experimental data, focusing on $u/Uref$ and $v/Uref$ and their fluctuations, can also provide a valuable metric for assessing the confidence level of the numerical measurements. These plots, shown in Fig. 12, include the 1-to-1 line along with the 1-to-2 and 2-to-1 reference lines to better contextualize the level of agreement.

The comparison between LES and wind tunnel experiments shows a generally good agreement in $u/Uref$⁠, with higher accuracy at greater heights. The FAC2 and HR values show that LES performs well, particularly above 18 m, where both metrics exceed 0.9. At lower heights, discrepancies become more evident, though FAC2 remains above 0.7, surpassing the 0.3 thresholds considered acceptable for urban flow simulations, as suggested by Hanna and Chang.⁵⁴ Similarly, HR values are mostly above 0.66, the admissible limit for this quantity.⁵⁵ These results suggest that LES provides a higher level of accuracy than is typically required for urban flow modeling, particularly in capturing the mean-flow field.

The scatter plots further illustrate this agreement [see Figs. 12(a) and 12(b)]. The LES predictions align well with the experimental data, with the majority of points clustering around the identity line (solid black line), especially for $u/Uref$ [Fig. 12(a)]. However, for $v/Uref$⁠, the scatter is noticeably larger [Fig. 12(b)], indicating a higher degree of variability between LES and experimental results. This discrepancy is consistent with the lower HR values for $v¯/Uref$⁠, suggesting that LES predicts larger lateral velocities than observed in the wind tunnel experiments. Given that LES employs a low-dissipation numerical scheme and achieves a high resolution of 0.75 m, it is unlikely that numerical errors are solely responsible for this difference. The wind tunnel measurements may also underestimate lateral velocity fluctuations due to physical constraints, such as probe alignment issues, wind tunnel sidewall effects, or Reynolds number mismatches. These effects can also be added to the fact that wind tunnel turbulence generation using spires and roughness elements may not fully replicate full-scale turbulence structures, particularly in the lateral direction. Consequently, the summation of numerical uncertainties along with those in the wind tunnel might be the reason for the larger deviations observed in this velocity component compared to those in the streamwise direction.

The rms velocity comparisons [Figs. 12(c) and 12(d)] provide additional insights into the representation of turbulence intensity. The LES captures the fluctuations in the streamwise component relatively well, especially at mid and upper heights (green, blue, and magenta points). However, at lower heights (red points, $h=2 m$⁠), there is a tendency for LES to slightly deviate from the rms values of the experiments. The velocity fluctuations for the lateral component tend to have a larger scattering along the identity line, with a slight underestimation compared to the experiments, especially at low and mid-heights. As discussed before, deviations at lower heights might, in general, be attributed to insufficient grid resolution necessary to represent small-scale turbulence near walls, but also to some limitations in the experimental setup. Given that urban flow turbulence is highly anisotropic, with complex wake interactions and shear layers, the Reynolds number mismatch between LES (full-scale) and the wind tunnel (1:200 scale) may lead to differences in how turbulence develops and is measured. While FAC2 and HR are useful for assessing LES accuracy, these metrics have limitations and should not be the sole basis for validation. FAC2 does not account for systematic biases; it simply evaluates whether LES predictions fall within a factor of two of experimental values. HR is more restrictive, but it still does not fully capture whether the LES reproduces the correct turbulence structures, spectral content, or flow dynamics. Given these limitations, energy spectra analysis is essential to assess whether LES captures the correct distribution of TKE across scales. This is particularly important in urban environments, where multi-scale turbulence, anisotropy, and wake interactions play a crucial role in flow development.

Figures 13 and 14 present the pre-multiplied energy spectra of the streamwise fluctuating velocity component for case 2. The energy spectra are non-dimensionalized by the standard deviation of the velocity fluctuations $σu=u′2$⁠, i.e., $Euu*=Euu/σu2$⁠, and the frequency is non-dimensionalized as $f*=f H/Uref$⁠. The selected locations are representative of the flow field in different areas, including open spaces enclosed by buildings, rooftops, streets, and road intersections. The exact location of each probe is given in Table III.

The energy spectra exhibit a well-defined inertial subrange with the expected $−2/3$ power-law decay when the energy density is multiplied by the frequency,⁵⁶ regardless of the numerical probe's location. For the probes located within the urban canopy layer (Fig. 13), the energy spectra decay spans just over one frequency decade. This behavior is consistent with the findings of Poggi and Katul,⁵⁷ who noted that within the urban canopy layer, the isotropic assumption becomes less valid due to wake production caused by interactions with buildings. As a result, the inertial subrange only emerges beyond the scales associated with wake production.

Above the urban canopy layer (Fig. 14), a clearer inertial subrange extending over approximately two decades is observed. This indicates a better separation between energy-containing scales and smaller dissipative scales, as turbulence becomes more isotropic with distance from the ground. Additionally, the energy spectra in both figures show a distinct filter cutoff around frequencies of $O(10)$⁠, where a noticeable change in slope occurs.

Several key conclusions can be drawn from the analysis of these results. First, an accurate capture of the inertial subrange highlights the high quality of the computations. The high-order methods used in the current implementation do not produce spurious oscillations at high frequencies or excessive numerical dissipation. The latter would result in a faster decay in the inertial subrange with a slope steeper than that predicted by the Kolmogorov hypothesis. Moreover, there is no energy pileup at high frequencies, which is typically indicative of poor resolution when using low-dissipation numerical schemes or insufficient dissipation of the SGS model.⁵⁸ 

Second, the resolution used in these simulations is sufficient to distinctly capture the separation between energy-containing scales and dissipative scales. This is particularly important in the context of LES, where the subgrid-scale model requires a computational mesh fine enough for the filter cutoff to occur well within the inertial subrange. Achieving this ensures that the simulations yield accurate results and effectively represent the underlying physics of turbulent flows.

In urban flow simulations, ensuring the convergence of turbulence statistics over time is essential for capturing the complex dynamics of wind interactions with buildings and street canyons. Due to the inherently unsteady and multi-scale nature of urban turbulence, statistical quantities such as mean velocity and turbulent fluctuations might require a relatively long averaging period to provide reliable results. Without proper temporal convergence, the assessment of flow structures, shear layers, and recirculation zones may be misleading, affecting the accuracy of numerical studies. Given the sensitivity of second-order statistics to turbulent fluctuations, a thorough evaluation of convergence is necessary to ensure that simulations yield robust and physically meaningful results. On the other hand, achieving temporal convergence in turbulence statistics can significantly increase computational costs. Thus, the integration time is a trade-off between statistical accuracy and computational feasibility, as longer simulations demand greater resources, particularly when high-resolution LES is used. As discussed in Sec. III C, statistics are collected over a period of 40FTs in both cases, which corresponds to 192T and 500T, respectively. It is worth noting that although this integration time is similar to those reported in the literature for comparable cases,^{13,23,59–61} a detailed assessment of time convergence is still necessary. In what follows, we provide an analysis of the time convergence of the computed statistics.

Figure 15 presents the temporal convergence study of case 1, examining both $u¯/Uref$ and TKE. Case 2 is not shown as it exhibits a similar convergence behavior. A comparison of three sampling lengths [see Figs. 15(c) and 15(d)] reveals that while 20FTs exhibit noticeable variations in both streamwise velocity and TKE, an averaging period of 40FTs is sufficient to achieve a temporal convergence of second-order statistics, i.e., TKE in this case. Profiles of $u¯/Uref$ converge substantially faster than those of TKE. This is an expected result, as TKE is affected by small-scale fluctuations which require a much longer integration to statistically capture the full range of fluctuations. This is, in particular, true in those zones featuring wake interactions, recirculations, separations, and reattachment, where flow complexity drives slower convergence of flow statistics.

The study has presented a wall-modeled large-eddy simulation (LES) of an atmospheric boundary layer (ABL) over idealized urban roughnesses. Two cases are considered: a three-dimensional cubic prism array and the “Michel-Stadt” urban model. These cases are selected to investigate the accuracy of LES in predicting urban turbulence and to validate the methodology against wind tunnel experiments. The simulations have been performed using a high-resolution, low-dissipation numerical scheme, with a spatial resolution below 0.75 m within the urban canopy, exceeding typical LES studies in the literature. A precursor simulation is used to generate realistic inflow conditions, ensuring that the simulated ABL matches the experimental setup.

The inflow of the study domain is driven by an online precursor simulation, ensuring a physically consistent ABL that matches the wind tunnel inflow conditions. This approach eliminates the need for synthetic turbulence generation and allows for a realistic turbulent inflow that dynamically adapts to the simulation of the study domain. A comparison of the inflow velocity profile with experimental measurements indicates excellent agreement, reinforcing the validity of the ABL development in the numerical model. The precursor simulation approach contributes to the accurate reproduction of mean-velocity profiles and turbulence statistics, further supporting the reliability of the LES setup.

The mean-velocity profiles and velocity fluctuations show excellent agreement between LES and experiments at most heights, particularly for the streamwise component. The LES accurately captures the mean-flow patterns, including wake regions and the acceleration over rooftops, which are key characteristics of urban flows. The turbulence kinetic energy (TKE) has also been well predicted at most locations, with only slight underestimations near the surface, which could be attributed to SGS modeling effects or unresolved small-scale turbulence near walls. However, for the spanwise turbulence fluctuations, LES shows a general tendency to slightly underestimate root mean square values at lower heights, which might be related to measurement constraints in the wind tunnel, including probe misalignment or blockage effects limiting lateral velocity fluctuations together with some insufficient resolution near the wall.

The validation against wind tunnel experiments demonstrates a strong agreement in the streamwise velocity component, particularly at higher elevations where discrepancies observed are minimal. The validation metrics, FAC2 and hit rate (HR), exceed the commonly accepted thresholds of 0.3 and 0.66, respectively, confirming that LES predictions are well within standard accuracy levels for urban flow modeling. However, larger discrepancies have been observed in the spanwise velocity component, which could be attributed to a combination of numerical and experimental factors, including differences in turbulence generation, lateral boundary effects in the wind tunnel, and Reynolds number mismatches between full-scale and wind tunnel conditions.

While the validation metrics (FAC2 and HR) confirm the reliability of LES in predicting mean velocities, they do not fully characterize turbulence structure or spectral energy distribution. A detailed spectral analysis has therefore been conducted to evaluate whether LES properly captures the turbulent energy cascade across scales. The energy spectra plots confirmed that LES maintains low numerical dissipation, as there is no excessive damping of small-scale fluctuations, and the inertial subrange follows the expected turbulence decay slopes. This is a key finding, as it validates that the LES formulation does not introduce excessive numerical diffusion, which could otherwise affect turbulence dynamics in urban flows.

The LES resolution of 0.75 m within the urban canopy, which is finer than typical LES studies in urban environments, allows for an accurate representation of turbulent structures. The spectral resolution analysis demonstrated that LES resolves a significant portion of the inertial subrange, ensuring that most of the energy-containing turbulence scales are explicitly captured. The Kolmogorov $−2/3$ slope in the pre-multiplied spectra has been well reproduced below and above the urban canopy, indicating that the energy cascade is correctly represented. However, inside the urban canopy, some deviations have been observed due to the influence of building-induced turbulence and wake interactions. These findings reinforce the high quality of the LES predictions and the suitability of the numerical scheme for urban flow modeling.

Overall, this study highlights the effectiveness of the current methodology for accurately capturing the physics of urban flows, particularly when high-order, low-dissipation schemes and fine spatial resolution are used. Despite the limitations of FAC2 and HR as validation metrics, the findings surpass standard validation criteria, emphasizing the need to complement the analyses with spectral decomposition to fully assess turbulence dynamics. Furthermore, this work sets the foundation for a more comprehensive study that will encompass a more in-depth analysis of turbulent quantities and roughness characterization.

This work has been partially financially supported by “Agència de Gestió d'Ajuts Universitaris i de Recerca” under the call CLIMA 2023 (ref. 2023 CLIMA 00097) and the APPWIND project (ref. PLEC2021-007943). O. Lehmkuhl's work is financed by a Ramón y Cajal postdoctoral contract by the Ministerio de Economía y Competitividad, Secretaría de Estado de Investigación, Desarrollo e Innovación, Spain (RYC2018-025949-I). The authors acknowledge the support of the Departament de Recerca i Universitats de la Generalitat de Catalunya through the research group Large-scale Computational Fluid Dynamics (ref. 2021 SGR 00902) and the Turbulence and Aerodynamics Research Group (ref. 2021 SGR 01051).

The authors have no conflicts to disclose.

Ming Teng: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Josep M. Duró Diaz: Data curation (equal); Software (equal); Visualization (equal); Writing – original draft (equal). Ernest Mestres: Data curation (equal); Software (equal); Writing – original draft (equal). Jordi Muela Castro: Methodology (equal); Writing – original draft (equal). Oriol Lehmkuhl: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Ivette Rodriguez: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In addition to the data comparison presented in Sec. IV, additional comparisons with experimental measurements are shown here to further reinforce the comprehensive validation of case 2. Figures 16–18 present the comparisons of $u¯/Uref$⁠, $v¯/Uref$ and $u′rms$ in the $x–z$ plane along $y/Href=−9.375$⁠, while those along $y/Href=9.375$ in the $x–z$ plane are shown in Figs. 19–21.

The minor discrepancies observed in all profiles between the two grid resolutions indicate the grid independence of the numerical results. A satisfactory agreement with experimental measurements is obtained in $u¯/Uref$ profiles across all locations (see Figs. 16 and 19). In comparison of $v¯/Uref$⁠, discrepancies are largely pronounced at a pedestrian level around building corners and road intersections [e.g., Figs. 17(c), 17(g), and 20(e)]. Although the simulations accurately predict the peak values of $u′rms/Uref$ above the roofs, an under-prediction within the reverse-flow regions between buildings is observed [see Figs. 18(b), 21(b), and 21(e)]. Additionally, similar under-predictions are also clearly identified at certain locations far above the roof-level [e.g., Figs. 18(a)–18(d)].