#### Abstract

Hazardous weather, turbulence, wind, and thermals pose a ubiquitous challenge to Unmanned Aircraft Systems (UAS) and electric-Vertical Take-Off and Landing (e-VTOL) aircrafts, and the safe integration of UAS into urban area requires accurate high-granularity wind data especially during landing and takeoff phases. Two models, namely, Open-Source Field Operation and Manipulation (OpenFOAM) software package and Weather Research and Forecasting (WRF) model, are used in the present study to simulate airflow over Downtown Oklahoma City, Oklahoma, United States. Results show that computational fluid dynamics wind simulation driven by the atmospheric simulation significantly improves the simulated wind speed because the accurate modeling of the buildings affects wind patterns. The evaluation of different simulations against six Micronet stations shows that WRF-CFD numerical evaluation is a reliable method to understand the complicated wind flow within built-up areas. The comparison of wind distributions of simulations at different resolutions shows better description of wind variability and gusts generated by the urban flows. Simulations assuming anisotropy and isotropy of turbulence show small differences in the predicted wind speeds over Downtown Oklahoma City given the stable atmospheric stratification showing that turbulent eddy scales at the evaluation locations are within the inertial subrange and confirming that turbulence is locally isotropic.

#### 1. Introduction

Hazardous winds and turbulence are key barriers to the safe and autonomous navigation of small Unmanned Aircraft Systems (UAS) and future Urban Air Mobility (UAM) systems. The risk of damage or loss of aircraft due to weather is significantly high when operating in an urban environment where wind speeds and directions can change unpredictably due to the urban canyon effect. These challenges are exacerbated in urban settings due to the Venturi effect, urban heat islands, wind shear, microbursts, microclimates, and surface forcings (roughness, displacement height, built fraction). As UAS become more autonomous, airspace will become more congested with electric Vertical Take-Off and Landing (e-VTOL) aircrafts and air taxis, the prediction of winds within the Planetary Boundary Layer (PBL) must continuously evolve to quickly identify and respond to off-nominal situations that are needed for a UAS navigation.

Simulating winds within the PBL is challenging because the turbulent exchanges of momentum and energy between the atmosphere and both natural and human-made structures happen within the low layers of the PBL. The Surface Layer (SL) occupies the lowest 10% of the PBL, where the shear-driven turbulence is dominant. In the SL, the subgrid-scale influences of the turbulent eddies in the PBL are parameterized using Monin–Obukhov Similarity Theory (MOST) of Obukhov [1, 2].

The SL has two parts: Inertial Sublayer (ISL) and Roughness Sub-Layer (RSL). The ISL is the upper part of the SL, where MOST is valid and vertical variation of the turbulent fluxes is negligible. The RSL is the layer near and within the surface roughness elements (e.g., trees and buildings) where routine UAS operations will take place, and urban corridors and vertiports are installed. However, the MOST is not applicable within the RSL [3–5]. Moreover, the RSL has not been explicitly considered in global and mesoscale atmospheric models because the PBL is coarsely resolved. In fact, mesoscale atmospheric models, such as the WRF model are based on a terrain-following coordinate system and smoothing processes for terrain data to deal with the complex terrain [6], and great terrain differences hence exist between reality and the mesoscale models. Some studies applied the Weather Research and Forecasting (WRF) in the Large Eddy Simulation (LES) mode (WRF-LES) to increase model resolution [7, 8], but the low vertical resolution and the need for a smoothing process for terrain data still exist. Moreover, the computational cost of a WRF-LES simulation is high. References [9, 10] developed the Immersed Boundary Method (IBM) within the WRF model (WRF-IBM) which is a new WRF module that can solve flows over steep terrains and overcome the grid cells skewness using the IBM of Iaccarino and Verzicco [11] and Mittal and Iaccarino [12]. However, this method requires very high resolution and high computational power to handle complex urban terrain and limit interpolation errors.

Moreover, using mesoscale models at the LES mode requires running intermediate simulations at horizontal resolutions that are within the grey zone of turbulence or Terra incognita [13] which refers to resolutions at which turbulent eddies in the PBL are partially resolved and partially parameterized. Turbulent eddies within the turbulence grey zone evolve in a grid-dependent way [14] and neither 1D PBL schemes nor LES formulations (such as Turbulent Kinetic Energy (TKE) of resolved eddies that is much higher than subgrid-scale modeled eddies) are tailored to model turbulence within this regime. Moreover, because of the small grid spacing, Reynolds Averaged Navier–Stokes (RANS) modeling is no longer valid and RANS assumptions (such as vertical fluxes that are fully parameterized and vertical velocity scales that are much smaller than horizontal ones) are not verified. Therefore, there is a growing necessity to find new numerical methods to represent these processes.

Fortunately, dedicated microscale models can partially compensate for the shortcomings of mesoscale models. First, Computational Fluid Dynamics (CFD) models can simulate the wind field with higher spatial resolutions (a few meters to tens of meters) than mesoscale models [15, 16]. Moreover, most microscale models are based on the finite volume method, which can improve their ability to depict realistic terrain [17].

Microscale prediction models can be classified into three broad categories based on the turbulence closure schemes: RANS, LES, and Direct Numerical Simulation (DNS) [18]. The choice among the three approaches is a balance between the cost and solution fidelity. The LES technique gained attention in both research and industry because it has the potential to provide more reliable results than RANS with lower cost than DNS. In addition, the ability of LES to resolve idealized urban canopy flows has been demonstrated using comparison wind tunnel and DNS results [19, 20]. The LES technique also has the advantage to have a weak dependency on Reynolds number and the choice of the Subgrid Scale (SGS) turbulence model [21].

However, these models have a shortcoming in coping with the boundary conditions and usually use simple wind profiles as the boundary condition. Therefore, coupling the mesoscale models that can handle atmospheric boundary layer physics, land surface-radiation interactions, and large-scale mean wind direction/speed variations with microscale models is one way of simulating the microscale wind field over complex terrain by providing more realistic boundary conditions to the CFD model that can provide a wind field simulation with much higher spatial resolution. Various meteorological microscale models in the community have been applied to simulate urban flows. For instance, the PArallel Large-eddy simulation Model (PALM) LES solver [22] was used by [21] to simulate local wind conditions in a real urban environment (Downtown Helsinki, Finland) under realistic atmospheric boundary layer conditions.

In recent years, a large amount of wind engineering and research on wind modeling has been conducted by using the atmospheric and CFD coupled system as computer hardware and computational resources are rapidly advancing. Reference [23] used WRF-LES to drive Open-Source Field Operation and Manipulation (OpenFOAM) simulations with an Unsteady Reynolds-Averaged Navier–Stokes (URANS) turbulence model to simulate the passage of a cold front over a wind-energy test site over a hill. The evaluation of wind speed against measurements from tower and UAS data showed that this coupling was able to simulate reasonably various flow micro features. Similar results are found in [24] who used the one-way offline coupling system to simulate wind fields in the Chough Mountain in China. Reference [18] coupled the CFD software OpenFOAM with the WRF model to investigate the atmospheric boundary layer structure, ambient wind, and pollutant dispersion within the built-up areas of Shenyang, China, and found that the OpenFOAM simulation driven by the mesoscale WRF simulation is a reliable method to understand the complex flows and dispersion around buildings. In this paper, we will use one-way offline coupling of LES using OpenFOAM and WRF-LES atmospheric simulation to simulate urban winds within the RSL.

This paper is structured as follows: in Section 2, used observations, the setup of the numerical simulations using atmospheric model and CFD model, and the coupling method were described. And the results and discussions of the numerical experiments were then presented in Section 3, and finally the conclusions were given in Section 4.

#### 2. Methods and Materials

##### 2.1. Observations

The Oklahoma City Micronet (OKCNET) is an operational surface observing network designed to improve monitoring across the Oklahoma City, Oklahoma, and shed insights on atmospheric processes such as urban heat island and thunderstorms evolution. Oklahoma City (OKC) is a rapidly growing urban area and its population increased by more than 10% between 2000 and 2010. In addition, most quality assured data are collected at airports and OKC is one of the rare urban areas where a measurement campaign, namely the OKCNET, was conducted by the Oklahoma Climatological Survey and between 2008 and 2013 within the urbanized center of the city in order to understand the complex nature of atmospheric processes within urban boundary layer, roughness sublayer in particular. Therefore, the Downtown OKC is a valuable area to study urban wind processes within the roughness sublayer. Figure1 shows the location of the state of Oklahoma, the study area of Downtown OKC.

The OKCNET consists of 40 stations including 4 MESOscale NETwork (MESONET) stations and 36 Micronet stations mounted on traffic signals. Six stations in Downtown OKC are retained for our study. Table 1 shows the characteristics of the six stations that will be used for comparison and evaluation. The data are collected using WINDCAP sensors with an accuracy of 0.3 m s^{−1} [25] and quality assured in near real-time at an interval of 1 min. Figure 1 displays the location of the six Micronet stations inside Downtown OKC.

##### 2.2. WRF Model Setup

The WRF model was used to provide the boundary conditions for the CFD model. In recent years, WRF was widely used in both academic research and industry [26, 27]. A fully compressible and nonhydrostatic dynamic framework is adopted in the ARW module. The simulated three domains (D1, D2, and D3) are shown in Figure 1. The horizontal resolutions of D1, D2, and D3 are 3 km, 1 km, and 400 m. The three domains are centered on 97.51529º W and 35.47239º N. Table 2 shows the three domains’ settings. The outermost D1, intermediate D2, and innermost D3 domains have 74 × 61, 112 × 97, and 150 × 130 grid points in the south-north and east-west directions. The WRF model contains 80 vertical levels and the lowest 30 levels are below 1 km. The adaptive time stepping is used to guarantee the numerical stability of the WRF model. Table 3 shows the configuration and the parameterizations used in the simulations over D1 and D2. For the simulation over D3, similar parameterizations are used except that no PBL parameterization is used, the subgrid-scale parameterization for small eddies modeling used is the 1.5 TKE closure. The LES assumptions such as TKE of resolved eddies is very much higher than the TKE of modeled eddies which are not verified for D3, and the simulation with a horizontal resolution of 400 m was used to test the impact of terrain smoothing on the dynamic of the flow. Two PBL schemes are tested for wind simulation over D2; one is nonlocal and the other is local: the Yonsei University (YSU) and the Mellor–Yamada–Janjic (MYJ) schemes. The YSU scheme is a first-order closure, nonlocal eddy-diffusivity (K theory) approach with explicit treatment of the PBL top [31]. This scheme provides an updated countergradient flux term originally included in the medium range forecast (MRF) scheme. Here, the PBL top is determined by the level at which minimum turbulent flux exists (heat, momentum, and moisture). The YSU scheme was updated in WRF version 3.4 to account for a bug in the code in the calculation of the momentum diffusion coefficient [32]. The author in reference [32] claims this improves wind speed estimation in stable regimes by decreasing vertical mixing and is one of the only PBL schemes to focus improvements in wind. The MYJ scheme is a local mixing, 1.5-order closure, prognostic TKE scheme, used operationally in the North American Mesoscale Forecast System (NAM). The TKE vertical distribution is controlled via a master mixing length, which is used in the TKE budget equation [33]. Local closure implies only small turbulent eddy contributions to the TKE distribution which are accounted for. A comparison between simulations over D2 using the two different cumulus parameterizations of different types is shown in Figure 2. The first one is Kain–Fritsch scheme [30] which is a mass flux or moisture convergence scheme where the different cloud processes including the convective updraughts, downdraughts, entrainment, and detrainment in clouds are modeled. The second one is Betts–Miller–Janjic cumulus scheme [33–35] which is an adjustment scheme that relaxes the large-scale environment towards reference thermodynamic profiles. This comparison over domain D2 shows a very small difference over the domain, and that difference between the two simulations can locally be as high as 4% because the atmosphere is mainly stable during the winter season. The Kain-Fritsch scheme is used in the WRF simulations setup.

##### 2.3. CFD Model Setup

As CFD model, the OpenFOAM software package version 6 which is a set of C++ libraries is used to solve complex problems in fluid mechanics. In this study, we use the pisoFoam solver which is a transient solver for incompressible and turbulent flow using the Pressure Implicit with Splitting of Operators (PISO) algorithm. The PISO is an efficient method to solve the Navier–Stokes equations in unsteady problems by resolving LES filtered equations.

The algorithm can be summed up as follows:(1)Set the boundary conditions.(2)Solve the discretized momentum equation to compute an intermediate velocity field.(3)Compute the mass fluxes at the cells faces.(4)Solve the pressure equation.(5)Correct the mass fluxes at the cell faces.(6)Correct the velocities based on the new pressure field.(7)Update the boundary conditions.(8)Repeat from 3 for the prescribed number of times.(9)Increase the time step and repeat from 1.

Steps 4 and 5 can be repeated for a prescribed number of times to correct nonorthogonality.

The schemes used to model various terms of the Navier–Stokes equations are summarized in Table 4.

Figure 1 shows the horizontal buildings configuration at the surface with the red boxes showing real buildings that are modeled. The real buildings and their complex shapes are simplified. The heights are found in https://skyscraperpage.com/diagrams/?stateID=64. The average height is 51 m with the tallest and lowest buildings that are 257.4 m and 30 m, respectively. A nonuniform mesh is 3-dimensional mesh containing hexahedral and split-hexahedral cells. The WRF data from one grid that was nearest to the location of the CFD domain was used to provide the initial and boundary conditions to OpenFOAM. Moreover, the boundary condition on the ground is set to “Wall” (Dirichlet), and the top boundary was set to “symmetry” (Neumann). The four lateral boundary conditions were set to “patch” boundary conditions. “Patch” allows a different type of boundary conditions depending on the prevailing wind direction. For example, if the prevailing wind direction is northeast, the northern and eastern domain boundaries are Dirichlet BC and the southern and western boundaries are Neumann for Velocity.

The mesh was generated using the snappyHexMesh utility. This option consists of conforming to the surface by iteratively, refining a starting mesh and morphing the resulting split-hex mesh to the surface.

##### 2.4. Coupling Method

The horizontal spacings of WRF model ranging from hundreds of meters to kilometers are different from those of the CFD model of several meters. Therefore, the coupling of WRF and CFD models is an extremely challenging task.

To solve this problem of resolution differences, we used Cressman interpolation [36].

The wind speed at a CFD grid point in the boundary faces only depends on the distances from the nearby WRF grid points as formulated in (1) where is the wind speed at a CFD grid point and N_{WRF} represents the number of WRF grid points that have influence on the CFD grid point. is the wind speed at a WRF grid point and represents the weighting function which depends on the distances between a P (x_{cfd}, y_{cfd}, z_{cfd}) and nearby WRF grid points (x_{i}, y_{i}, z_{i}) as shown in Figure 3. Two weighting functions are used here: the Cressman and Barnes schemes as explained in sections {5.1} and {5.2}, respectively:

The communication from the parent to the child domains is ensured by an hourly data flow through boundary conditions. The interpolation in space is acting as a coupling from WRF grid to the CFD grid through the initial and boundary conditions and we simulate the airflow over buildings in time.

The fields communicated to OpenFOAM are the three velocity components. The interpolation method consists of two steps: (1) the coordinates of the cell’s centers of the four lateral faces and the top face are transformed from Cartesian form to geographic form. (2) 3D interpolation in space of the fields from instantaneous WRF produce outputs every hour.

In addition, the CFD model is run for an hour and is initialized from the resolved fields of the previous run. A spin-up of one hour is not sufficient to reach a balance of spine-up effects. Reference [37] found that a spin-up of several days is needed until equilibrium is reached. References [37, 38] used a 24-hour spin-up time to balance the spin-up effects. In this study, given the computational capacity and project timeline, a spin-up of one hour is used.

###### 2.4.1. Cressman Method

Figure 3 shows a sketch map of the Cressman interpolation function where P represents a CFD grid point and O_{1}, O_{2}, and O_{3} represent the WRF grid points. The black circles indicate the influence range of each WRF grid point. It can be seen in the example given in Figure 3 that the wind speed at point P is influenced only by O_{1} and O_{2} but not O_{3}. The weighting function is defined by (2) where r_{L} and r_{z} represent the horizontal and vertical influence radii of a WRF grid point respectively. The Cressman weighting function depends on the horizontal (d_{xy}) and vertical (d_{z}) distances between the CFD grid point of coordinates (x_{p}, y_{p}, z_{p}) and nearby WRF grid points of coordinates (x_{i}, y_{i}, z_{i}) as defined by (3) and (4), respectively. Since the WRF data are not sparse, both horizontal and vertical radii, r_{h} and r_{z}, of influence are set as 1 km assuming the isotropy of turbulence. Paper [39] used a radius of influence of 520 km which is about twice the mean station separations (250). That is why a value for the radius of influence of 1 km that is approximately twice the average spacing of the WRF data (400 m) tends to be a reasonable compromise between undersmoothing and oversmoothing. The assumption of turbulence isotropy will be discussed in Subsection 3.5 where two simulations will be compared: one with r_{h} = r_{z} = 1 km and the other with r_{h} = 1 km and r_{z} = 500 m to represent turbulence anisotropy.

###### 2.4.2. Barnes Method

Because Barnes scheme was initially designed for horizontal filters, we are using the generalized 3D Barnes scheme developed by Lu and Broning [40] as defined by (5) where *L*_{h} and are the horizontal and vertical length scale parameters:

##### 2.5. Simulations’ Summary

All performed simulations, their types, domains and configurations and used parameterizations are summarized in Table 5.

##### 2.6. Model-Measurements Comparison Method

In order to evaluate the accuracy and temporal correlation of numerical simulations, the statistical scores defined in Table 6 are calculated: the simulated mean (), the observed mean (), the Root Mean Square Error (RMSE), the correlation (R), and the Mean Bias (MB).

#### 3. Results and Discussions

##### 3.1. General Evaluation

The simulated wind speeds using different simulations are evaluated against observations over the six Micronet stations and their statistical comparison is shown in Table 7.

During the simulation period and because of the low surface heating, the atmosphere is mainly stable over D2 as shown in Figure 4 where the vertical gradient of the potential temperature between 20.0 and 2.5 m is positive. The atmosphere is neutrally stratified in some sparse areas mainly located in the north of D2.

The observed wind speeds slow low values lower than 3 m s^{−1} over the six stations. All simulations overestimate the wind speeds: the simulated means range from 2.17 to 4.36 m s^{−1} against the measured mean of 1.99 m s^{−1}. The bias of 16% of the observed mean is 0.3 m s^{−1} bias on average, which is at the minimal detectable difference of the WINDCAP sensor. The coupled simulations (S5 and S6) have better performance than WRF simulations (S1, S2, S3, and S4) with higher correlations (58%), lower errors (0.70–0.72 m s^{−1}), and biases (16–18%); the simulated means of S5 and S6 are 2.03 and 2.17 m s^{−1}, respectively, against the observed wind speed (1.99 m s^{−1}). This difference is mainly due to the very high resolution of the coupled simulations and a better representation of the microscale terrain features (buildings) and wind turbulence within the urban canopy. S5 was able to accurately predict wind speeds with a correlation of 58%, and the lowest mean bias of 16% and S1 simulation had the higher bias with 236% because of the coarse resolution (3 km) compared to all other simulations. All WRF simulations (S1, S2, S3, and S4) overestimate wind speeds with a mean bias exceeding 100% because of uncertainties on the terrain topography and BL parameterization. The use of LES model in the WRF simulation with a 400 m horizontal resolution shows a negligible difference as the simulated mean decreased from 3.42 to 3.39 m s^{−1} and the MB slightly decreased from 143 to 140%.

According to wind speed’ time series shown in Figure 5, all WRF simulations, namely, S2, S2, S3, and S4, show higher discrepancies during the second simulation day as the four simulations predict more than twice the observed hourly wind speeds. Figure 6 shows horizontal slices of the simulated wind speed and direction over the CFD domain at 2 m, 35 m, and 200 m using S5 in two different times: January 1 at 12 pm and January 2 at 6 pm.

On January 1 at 12 pm, the winds are blowing from the NW direction with an angle of 314 degree with 5.56 m s^{−1} in a time when winds are blowing from the SSE direction with an angle of 154 degree with a speed of 1.11 m s^{−1}. Figure 6 also displays how the air flow behaves around the building blocks and demonstrates microscale wind flow features such as vorticity behind and between buildings, Venturi effects observed between buildings, and corner acceleration in which wind speeds almost double due to sharp buildings shapes. Downwashing is also observed because of the difference in buildings heights creating turbulent recirculations.

The statistical comparisons shown in Tables 8–10 show that simulations S1, S2, and S3 are predicting similar wind speeds over the six stations as all stations are located at the same 3 km and 1 km grid cell. S4 shows a different wind speed over KCB105 because it is located in a different 400 m grid cell than the five other stations.

Moreover, the computational times between an OpenFOAM and WRF simulation are compared. Given the resolution difference, OpenFOAM simulation takes approximately 7 times more than a WRF simulation of 1 km resolution over nested domain D2. Table 2 shows the central processing unit (CPU) time of 1 simulated hour using S2 and a LES using OpenFOAM.

##### 3.2. Impact of the Horizontal Resolution

The comparison between the two mesoscale simulations S1 and S2 showed a slight improvement of the predicted wind speeds: the simulated means over all stations decreased from 4.36 to 3.42 m s^{−1}, respectively. This slight improvement is due to the better representation of topography and terrain height as shown in a grid cell-based representation in Figure 7: the terrain height of D1 ranges from 180 to 660 m.a.s.l. in a time when it ranges from 250 to 450 m.a.s.l. over D2. The comparison also shows that although a small improvement of simulated wind speeds over the six stations with the increasing horizontal resolution as the MB improved from 238% to 143%, discrepancies remain because of uncertainties on the built fraction, terrain height, and features. Although the LES simulation S4 shows a negligible difference as the simulated mean decreased from 3.42 to 3.39 m s^{−1}, S4 was not able to resolve turbulent motions of urban wind because of uncertainties on terrain topography indicating that smoothing the terrain from 1000 m to 400 m does not appear to affect flow dynamics. Besides, relevant eddies have scale below 400 m within both daytime and nighttime PBL and cannot be captured by this discretization setup. As mentioned by Cuxart [42], the minimum horizontal resolution to be used in LES that is falling into the inertial subrange of the turbulence spectrum is 100 m.

##### 3.3. Impact of PBL Scheme

The comparison of the S2 and S3 shows that S2 is closer to the observations and that YSU parameterization predicts slightly better the wind speeds than MYJ parameterization: the simulated mean using S2 is 3.42 m s^{−1} which is lower than the one of S3 which is 3.69 ms^{−1} and the mean bias of S2 (143%) is lower than the one with S3 (161%). This difference is mainly due to the fact that the YSU scheme takes into account the nonlocal momentum mixing and that the YSU treats vertical mixing more realistically because turbulent eddies can span. Figure 8 shows the relative difference of wind speed over D2 between S3 and S2. The wind speed at 10 m.a.g.l. differences is as high as 20% over D2 due to the impact of the PBL scheme.

##### 3.4. Impact of Coupling Scheme

The comparison between S5 and S6 shows that simulated wind speeds using S5 are slightly closer to observations 2.03 m·s^{−1} (S5) than 2.17 m·s^{−1} using S6. Therefore, Barnes interpolation slightly overestimates the wind speeds at boundaries of the domain leading to a slight overestimation of wind speeds. Figure 9 shows that this difference is less than 2% and is mainly located over turbulence vorticity areas, namely, buildings corners and wakes and between the buildings where Venturi effect is visible.

##### 3.5. Impact of Horizontal Resolution on Wind Gusts Prediction

Wind gusts are also an important safety challenge for unmanned aviation especially using small UAS. A comparison of wind gust predictability is shown here by comparing the upper tails of the probability density functions of the predicted wind speeds using different simulations: S1, S2, S4, and S5.

Figure 10 shows that the predicted wind variability explored using PDF is different among the four simulations and Table 11 shows the mean and standard deviation of the four PDFs. The wind PDFs of S1 and S4 have a Gaussian/lognormal shape in a time when wind PDFs of S2 and S5 show an anomaly with a second peak at a wind speed of 1.3 m s^{−1} approximately. The highest wind speed variability (2.60 m s^{−1}) is found using the coupled simulation S5 where increasing wind speeds and gusts are predicted. This difference is forced by the coupling method due to the fact that terrain features are better modeled in S5 instead of a flat terrain as modeled using the WRF model.

The PDF of the coupled simulation S5 and S4 has the lowest variability (1.12 m s^{−1}) confirming that the 400 m resolution falls under the turbulence grey zone and LES assumptions are too coarse and not sufficient to provide accurate representation of vertical transport and wind gusts even in winter stable atmospheric conditions.

##### 3.6. Isotropy versus Anisotropy of Turbulence Modeling

Depending on the atmospheric stability conditions, turbulence could be anisotropic [43]. The anisotropy is taken into account by modifying the horizontal and vertical length scales in the coupling process. The assumption of turbulence isotropy is then examined here by comparing two simulations: S5 where Lh = Lv = 1 km and where Lh = 1 km and Lv = 500 m. Figure 11 shows the absolute difference between the two simulations on January 1 at 1 am at a height of 4 m. The simulation with anisotropic turbulence predicts higher wind speeds. The difference is mainly located where the flow cross-sections become constricted, over buildings edges, and can be as high as 0.7 m s^{−1} which is almost 12% of the wind speed predicted by WRF-1 km and WRF-LES-400 m. This result confirms the local isotropy hypothesis in Kolmogorov’s theory [44] that is related to turbulence scales within the inertial subrange. Similar conclusion was found in the work of Chan [45] by comparing eddy dissipation rate (EDR^{1/3}) from both a Doppler lidar and a 20-Hz anemometer at the Hong Kong International Airport (China).

#### 4. Conclusions

This study examines the two-day simulation of urban airflow over Downtown Oklahoma by using a WRF-OpenFOAM coupled model. Four WRF numerical simulations were performed using different horizontal resolutions and different boundary layer parameterization. Improved wind speeds were predicted with a 1 km WRF simulation. The WRF-LES simulation of 400 m resolution was not able to capture small scale wind eddies indicating that the representation of the terrain in the built-up area is of paramount importance to accurately predict urban winds which proves that the WRF setup and schemes are not tailored to resolve urban winds at the horizontal resolutions that fall within the turbulence grey zone. WRF simulations with two different urban canopy models are not shown here, but no significant difference was found between their predicted wind speeds. By comparing wind distributions predicted by simulations at different resolutions, the modeling of wind gusts was better represented using coupled simulations. Sensitivity study to turbulence isotropy and anisotropy showed that during the simulated two days of stable stratification, turbulence can be considered as isotropic and that dominant scales within the urban canopy layer are within the inertial subrange. More in-depth analysis of turbulence scales can be conducted in the future once higher-frequency data from anemometers or remote sensing instruments are available using a spectral analysis.

From this preliminary study, the simulations that are able to accurately predict the means of wind speed over the six urban stations are the coupled simulations where the atmospheric WRF-LES model is used to provide boundary conditions to the OpenFOAM model in which terrain topology is better represented. Future work should focus on sensitivity studies to run OpenFOAM without the buildings to validate the effect of resolution of input data and compare it with the effect of explicit resolution of turbulence.

These simulations can be the basis of decision-support tool for areas to avoid and dynamical geofencing for UAS operation based on vehicles tolerance to wind speeds. It was found that urban airflows were complicated in the presence of real building clusters and under varying atmospheric boundary structures. The slight difference was found between coupled simulations using Barnes and Cressman interpolation schemes. More robust comparison of winds aloft can be performed using data from remote sensing instruments such us radars or lidars.

Future work will focus on ingesting real observations from optimally placed anemometers into the WRF-CFD model and impact will be assessed in terms of improving the predicted wind speed and more accurately reproducing hyperlocal wind effects. In order to operationalize this model, a lot of sensitivity studies and tests mentioned in the paper need to be conducted in addition to the acceleration of the CFD simulations using machine learning or model reduction techniques.

#### Data Availability

The data used to support the findings of the study can be obtained from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was funded by University of North Dakota. Oklahoma Mesonet data have been provided through the courtesy of the Oklahoma Mesonet, which is jointly operated by Oklahoma State University and the University of Oklahoma. Continued funding for maintenance of the network was provided by the taxpayers of Oklahoma.