## Abstract

The validity of omitting stability considerations when simulating transport and dispersion in the urban environment is explored using observations from the Joint Urban 2003 field experiment and computational fluid dynamics simulations of that experiment. Four releases of sulfur hexafluoride, during two daytime and two nighttime intensive observing periods (IOPs), are simulated using the building-resolving computational fluid dynamics model called the Finite Element Model in 3-Dimensions and Massively Parallelized (FEM3MP) to solve the Reynolds-averaged Navier–Stokes equations with two options of turbulence parameterizations. One option omits stability effects but has a superior turbulence parameterization using a nonlinear eddy viscosity (NEV) approach, and the other considers buoyancy effects with a simple linear eddy viscosity approach for turbulence parameterization. Model performance metrics are calculated by comparison with observed winds and tracer data in the downtown area and with observed winds and turbulence kinetic energy (TKE) profiles at a location immediately downwind of the central business district in the area labeled as the urban shadow. Model predictions of winds, concentrations, profiles of wind speed, wind direction, and friction velocity are generally consistent with and compare reasonably well to the field observations. Simulations using the NEV turbulence parameterization generally exhibit better agreement with observations. To explore further the assumption of a neutrally stable atmosphere within the urban area, TKE budget profiles slightly downwind of the urban wake region in the urban shadow are examined. Dissipation and shear production are the largest terms that may be calculated directly. The advection of TKE is calculated as a residual; as would be expected downwind of an urban area, the advection of TKE produced within the urban area is a very large term. Buoyancy effects may be neglected in favor of advection, shear production, and dissipation. For three of the IOPs, buoyancy production may be neglected entirely; for one IOP, buoyancy production contributes approximately 25% of the total TKE at this location. For both nighttime releases, the contribution of buoyancy to the total TKE budget is always negligible though positive. Results from the simulations provide estimates of the average TKE values in the upwind, downtown, downtown shadow, and urban wake zones of the computational domain. These values suggest that building-induced turbulence can cause the average turbulence intensity in the urban area to increase by as much as 7 times average upwind values, explaining the minimal role of buoyant forcing in the downtown region. The downtown shadow exhibits an exponential decay in average TKE, whereas the distant downwind wake region approaches the average upwind values. For long-duration releases in downtown and downtown shadow areas, the assumption of neutral stability is valid because building-induced turbulence dominates the budget. However, farther downwind in the urban wake region, which is found to be approximately 1500 m beyond the perimeter of downtown Oklahoma City, Oklahoma, the levels of building-induced turbulence greatly subside, and therefore the assumption of neutral stability is less valid.

## 1. Introduction

Because of large populations and difficult evacuation logistics, urban areas are the most consequential locations for an atmospheric release of hazardous material, whether the release is due to an industrial accident or to an act of terrorism. To protect urban populations, the use of observational and/or modeling tools is necessary in order to track and forecast the transport and dispersion of the hazardous material from such releases. The need for robust modeling tools was stressed in a recent report by the National Research Council of the National Academies (2003). Specifically, the report states that two categories of models are needed: 1) advanced rapid-execution models for emergency response applications, and 2) slower but more accurate computational fluid dynamics (CFD) models for situations such as emergency preparedness and postevent analyses.

Several CFD models are available and are continually being tested with data from actual releases of tracer material in urban areas such as Salt Lake City, Utah (e.g., Chan and Leach 2004), Oklahoma City, Oklahoma, and New York City, New York. These releases have taken place at various times of day and in all seasons. Yet many modeling tools often assume that nearly neutral conditions are appropriate to use for the entire urban area being simulated, perhaps because of the lack of representative field measurements or sufficiently sophisticated turbulence models. A strong argument supports the application of neutral stability to urban environments: atmospheric stability (as defined by the Richardson number) is determined by both mechanical stresses and buoyant forcing. For an urban setting, building-induced mechanical stresses are thought to dominate the production of turbulence, overriding any stability effects and creating nearly neutral atmospheric conditions at least in the most densely built-up urban areas and immediately downwind of these areas.

Results from our recent simulations of two releases of sulfur hexafluoride gas in the Oklahoma City urban area using a CFD model—the Finite Element Model in 3-Dimensions and Massively Parallelized (FEM3MP)—appear to support partially the assumption that urban areas tend toward neutral stability. Based on a model–data comparison for winds and concentration in the near field and velocity and turbulence profiles in the urban wake region, Chan and Leach (2007) observed that the neutral stability assumption appears to be valid for intensive operation period (IOP) 9 (a nighttime release with moderate winds) and also appears to be valid for IOP 3 (a daytime release with strong buoyant forcing) in the urban core area but is less valid in the urban wake region.

In this study, a more definitive characterization of stability conditions associated with the simulated releases is presented. We document the neutral stability caused by the enhancement of mechanical mixing in the urban core region and the decay of that enhanced turbulence downwind. We first describe our CFD model and the field experiments simulated. We present results from a model–data comparison for both the wind and concentration fields from four releases, using two different turbulence parameterizations, one considering buoyancy effects with a simple linear eddy viscosity (LEV) approach for modeling mechanical production of turbulence and one omitting stability effects but with a superior nonlinear eddy viscosity (NEV) approach for modeling turbulence processes. We next examine observed turbulence kinetic energy (TKE) budgets from a profile located immediately downwind of the urban area, with special attention to the role of observed buoyant forcing relative to the total TKE in the urban wake. Last, we compare the predicted TKE fields within the central business district urban zone with upwind and downwind zones to delineate the regions in which urban-generated turbulence dominates other means of turbulence production.

## 2. The CFD model

Our model, developed under the sponsorship of the U.S. Department of Energy and Department of Homeland Security (DHS), is based on solving the three-dimensional, time-dependent, incompressible Navier–Stokes equations on massively parallel computer platforms. The numerical algorithm is based on finite-element discretization for effective treatment of complex building geometries and variable terrain, together with a semi-implicit projection scheme and modern iterative solvers developed by Gresho and Chan (1998) for efficient time integration. Physical processes treated in our code include turbulence modeling via Reynolds-averaged Navier–Stokes (RANS) and large-eddy simulation approaches (Calhoun et al. 2005), atmospheric stability, aerosols, UV radiation decay, surface energy budgets, and vegetative canopies, and so on. For convenience of code parallelization and computational speed, our current version of FEM3MP employs a structured mesh (but graded and distorted mesh is allowed) and explicitly resolved buildings within the computational domain are represented as solid blocks with velocity, pressure, and diffusivities set equal to 0.

While high-resolution CFD models are very useful for emergency planning, vulnerability analyses, and postevent analyses, such models usually require excessive computer resources and long turnaround times, thus rendering them unwieldy. As a means to improve the efficiency of FEM3MP, a simplified CFD approach, in which only targeted and important buildings are explicitly resolved using fine grid resolution and the remaining buildings are represented as drag forces (virtual buildings) using much coarser grid resolution, was implemented and demonstrated by Chan et al. (2004). The virtual building approach is discussed in detail in Chan and Leach (2007).

The present simulations resolve certain buildings with fine grid resolution and treat the remaining buildings as drag elements (virtual buildings) with coarser grid resolution. Buildings in the near-field, within 500 m of tracer releases, and buildings taller than 40 m are typically defined explicitly. Because the mixed explicit–virtual building approach reduces the number of grid points by approximately a factor of 10 (depending on the user’s choice of which buildings to make virtual), computational savings are also approximately an order of magnitude. However, the solution in the vicinity of virtual buildings is slightly less accurate because of coarser grid resolution and the presence of small values of velocity and diffusivities at the grid points associated with virtual buildings.

Turbulence within FEM3MP is typically parameterized with the NEV turbulence parameterization, which solves three turbulence equations as discussed in Gresho and Chan (1998). However, one limitation of our present version of the NEV turbulence model is that it assumes that the atmosphere being simulated is neutrally stable, that neither strong convection nor strong stratification is present. The present study seeks to determine the impact of such assumption on the ability of FEM3MP to simulate the continuous releases of sulfur hexafluoride (SF_{6}) during Joint Urban 2003 IOPs 2, 3, 8, and 9. For comparison, simulations utilizing a simpler (LEV) turbulence parameterization are also considered. The LEV parameterization (Chan et al. 1987; Calhoun et al. 2004) does consider the effects of atmospheric stability via a Monin–Obukhov-based stability function, but it is based on *K* theory for turbulence parameterization and is therefore less sophisticated than the NEV parameterization.

## 3. The Joint Urban 2003 field study

To provide quality-assured, high-resolution meteorological and tracer datasets for evaluation and validation of indoor and outdoor urban dispersion models, the U.S. DHS and Department of Defense–Defense Threat Reduction Agency cosponsored a series of dispersion experiments, named Joint Urban 2003 (JU2003), in Oklahoma City (OKC) during July 2003 (Allwine et al. 2004). These experiments took place during daytime and nighttime to include both convective and stable atmospheric conditions. A total of 10 IOPs were conducted and SF_{6} in the form of puffs or continuous sources was released over 6 daytime and 4 nighttime episodes.

Many wind and concentration sensors were used to collect wind and SF_{6} data over both long (up to 30 min) and short (0.1 s) time-averaging periods. In addition to measurements near the surface, wind and concentration profiles adjacent to the outside walls of several buildings were also taken. In one nocturnal case, balloons were deployed close to the tracer release area. Many of the released balloons exhibited quick ascents from ground level to the top of buildings, implying localized convection. For the present study, model predictions of SF_{6} concentrations are compared with 15-min averages taken from “Blue Boxes,” which collected ten 20-cm^{3} “sips” of air into a 200-cm^{3} bag over 15-min intervals from an intake valve about 1.5 m above the ground surface. The air samples from each bag were then analyzed using a Hewlett-Packard 5890 Series II gas chromatograph (GC) equipped with a model G1223A electron capture detector. The GC was calibrated using 49 SF_{6} standards ranging in concentration from 9.3 parts per trillion (ppt) to 980 parts per billion (ppb). Five separate contiguous calibration curves were constructed to quantify samples throughout the entire range. Extrapolation of the low-level calibration curve with intercept forced through the origin was used to quantify samples below 9.3 ppt. Samples above 980 ppb were diluted to provide SF_{6} concentrations within the calibration range. A set of 14 calibration check standards was used to monitor system stability and selected samples were run in replicate.

Many of the in situ meteorological observations during JU2003 were within the urban canopy layer, at or below the mean height of the buildings in the immediate vicinity of the observations. Particularly useful for comparison with CFD predictions in the urban area are the Dugway Proving Ground (DPG) Portable Weather Information Display Stations (PWIDS) wind observations, collected by anemometers mounted on streetlights.

Measurements immediately downwind of the central business district (CBD) were obtained using a pseudotower (a ladder under tension, anchored by a construction crane) located about 750 m north of downtown OKC. At eight levels (7.8, 14.6, 21.5, 28.3, 42.5, 55.8, 69.7, and 83.2 m) along this crane, R. M. Young Co. Model 81000 ultrasonic anemometers were mounted, providing wind speed, wind direction, and virtual temperature measurements at those levels. The highest observation level was well above the mean building height of the Oklahoma City CBD (approximately 35 m) but below the height of the two tallest single buildings (approximately 120 m). Details on the construction of this pseudotower are presented in Gouveia et al. (2007). Virtual temperature profiles from this platform are not suitable for estimating atmospheric stability.

Each 30-min time series of 10-Hz data from the crane sonic anemometers is processed in the following fashion. First, to correct for any tilting in the sonic anemometer as mounted, the planar-fit tilt correction algorithm (Wilczak et al. 2001) is applied to the data; each 30-min time series is corrected independently. After the tilt correction is applied, the mean wind direction over the 30-min segment is calculated. In the event that the wind direction is northwesterly to northeasterly (315°–360° or 0°–45°), the data are rejected from further inspection because of expected contamination by the supporting crane. Northerly winds did not occur during any of the IOPs discussed here. The horizontal velocity data are rotated into a right-handed natural coordinate system, so that the streamwise coordinate *u* is aligned with the mean wind and the transverse component *υ* is perpendicular to *u* in the horizontal plane; *w* remains perpendicular to *u* in the vertical plane. For momentum and heat flux calculations, the mean of each time series over 30 min is removed. Measurements of turbulent kinetic energy (TKE) dissipation rate are also obtained from these data, using the inertial estimation method as described in Piper and Lundquist (2004).

Most terms of the TKE budget can be calculated from the crane dataset. The TKE budget is written

as in Stull (1988), or

where *T* represents the local storage or tendency of TKE, *A* represents the advection of TKE, *B* represents buoyant production or destruction of TKE, *S* represents shear production of TKE, TT represents the turbulent transport of TKE, *P* represents the pressure transport term, and *D* represents the dissipation rate of TKE. From the profile of sonic anemometers at the crane, *D,* TT*, S,* and *B* may be calculated, either uniquely from observations at each level, as for dissipation and buoyant production, or from spline fits to the profile of measurements, as for turbulent transport and shear production. The local storage or tendency term may also be estimated from the change in TKE at each level from one time series to the next time series; for the IOPs studied here, the storage term calculated in that fashion is 1–2 orders of magnitude smaller than *S*, the shear production term. The estimation of the pressure transport term would require sensitive pressure fluctuation instrumentation, and this term would include the transport or destruction of TKE by wave motions. The advection term *A* is expected to be the largest term at a location in an urban shadow. Without similar observations upwind and downwind of the crane, this advection term must be calculated as a residual, which lessens the confidence in this estimate as it now includes all errors from the calculations of the other terms, as well as any influence from pressure transport *P*.

To estimate the contribution of buoyant forcing to the total amount of TKE observed, the rate *B* of buoyant production of TKE is multiplied by a turbulent time scale *τ*, which is determined from the quotient of mean observed TKE over mean observed dissipation rate, following the model of Zeierman and Wolfshtein (1986). As discussed below, observed values of this time scale range from 60 to 170 s for these IOPs. This turbulent time scale is smaller than a typical convective turbulence time scale (approximated by the quotient of boundary layer height *h* over a convective velocity *w**) of 600 s. Shorter turbulent time scales are expected in an urban boundary layer given the local circulations (both thermal and mechanical) generated in such an environment.

## 4. Model–data comparison

### a. Model simulation parameters

Airflow and dispersion simulations for the first continuous release of IOPs 2, 3, 8, and 9 were performed. The first two of these IOPs occurred during the day, at 1100 local time (LT), which was central daylight time. The latter two occurred at night, at 2300 LT. In each case, SF_{6} was released near the ground as a point source for 30 min, with varying release rates as summarized in Table 1.

For all numerical simulations, a domain size of 1030 m × 3010 m × 425 m (in lateral, longitudinal, and vertical directions) was employed. A graded mesh consisting of 201 × 303 × 45 grid points, with a minimal grid spacing of ∼1 m near the ground surface, was used. Most of the buildings within 500 m of the release point were explicitly resolved, and the remaining buildings were treated as virtual buildings. The model domain and building representations appear in Fig. 1.

Steady logarithmic velocity profiles provided inflow boundary conditions. These profiles were created based on the 15-min averaged wind speeds and directions from the Pacific Northwest National Laboratory (PNNL) sodar located approximately 2 km SSW of downtown OKC and the hourly averaged data from the weather station on the rooftop of St. Anthony’s hospital at ∼1.5 km NW of downtown OKC. The 50-m wind speeds and wind direction for these four releases are summarized in Table 1; wind speeds varied from 5 to 7.2 m s^{−1}, and wind directions varied from 155° to 215°. IOP 2 was simulated twice, changing the wind direction only. The logarithmic inflow wind profiles appear with the results of the simulations in Figs. 2 –5.

For each simulated release, a quasi-steady-state flow field was established after ∼15 min of simulated time prior to the start of the dispersion simulation. The release of SF_{6} was modeled as a continuous source over a small area (covered by 2 × 2 cells on the ground surface) at a constant release rate and dispersion results indicate steady state was reached in about 20 min of simulated time. For all cases, the RANS approach with the three-equation NEV turbulence parameterization (Gresho and Chan 1998) was used and neutral atmospheric stability was assumed. Surface-based convection was not imposed because of the lack of appropriate and representative observations to provide boundary conditions. A second run using the LEV turbulence parameterization was also carried out for each IOP, as summarized in Table 1. The LEV parameterization does include the effects of atmospheric stability, but its modeling of mechanically produced turbulence is less sophisticated than that available in the NEV parameterization, as discussed above.

FEM3MP’s predictions of flow and tracer concentrations in the near and intermediate regions of the release point are presented and compared with observations. For brevity, only major results are presented and compared herein. More detailed model–data comparisons can be found in Chan and Leach (2007) for IOPs 3 and 9. Several of the statistical performance measures recommended by Hanna et al. (1993, 2003) are used to indicate the performance of our model. They are the fractional bias (FB), the geometric mean bias (MG), the normalized mean square error (NMSE), and the fraction of predictions within a factor of 2 or 5 (FAC2 or FAC5). For weighted differences in angles between predicted and measured wind vectors, the formula of scaled angle differences (SAA) of Calhoun et al. (2004) is used. As summarized in Warner et al. (2004), FB and MG measure the systematic bias of a model in terms of absolute differences and ratios; a perfect model would result in a FB of 0 and a MG of 1. NMSE measures the scatter associated with the predictions relative to observations, and for a perfect model would equal 0, as would SAA. To our knowledge, no community consensus defines acceptable ranges of these metrics.

The equations for these metrics are defined as

where *C* is the observation of interest (e.g., wind speed or tracer concentration), *C _{p}* is the model prediction, and

*C*is the observation. Overbars denote averages. The SAA, a model performance metric calculated from wind speeds and wind directions, is given by

_{o}where *ϕ _{i}* is the angle between predicted and observed velocity vectors and

*N*is the number of samples being averaged. The angle difference is scaled by the magnitude of the predicted velocity vector |

*U*| and then is normalized by the average of the magnitudes over all samples. By scaling the angles by the magnitudes, this metric weights the angles of the larger vectors more strongly to minimize the relatively unimportant errors in wind direction associated with small wind speeds.

_{i}### b. Model performance in the near field

#### 1) Flow field

Airflow in urban areas is extremely complex, with features such as flow separations, local stagnation regions, eddies of various size, and high velocity jets in street canyons. These features were all observed in our model simulations, as documented in Chan and Leach (2007). Visualization of the simulated flow fields are presented in Chan and Leach (2007) for IOPs 3 and 9. For a quantitative model–data comparison of wind vectors at a number of locations for IOPs 2, 3, 8, and 9, we use observations from the 20 DPG PWID sonic anemometers described above, located within the central business district zone. The observed wind speeds and wind directions at these locations for each 30-min period corresponding to the releases simulated with FEM3MP have been compared with the FEM3MP predictions. These statistical performance measures, based on the 20 PWID data points, are summarized in Table 1.

In all cases, the absolute value of the model’s fractional bias (FB) is reasonably good, less than 0.25. A tendency toward overprediction of wind speeds is seen, with the fractional bias regularly between 0 and −0.25, regardless of whether the sophisticated neutral-stability turbulence parameterization (NEV) or the simpler varying-stability turbulence parameterization (LEV) is utilized. Only for IOP 8 is the FB greater than zero. FEM3MP’s systematic overprediction of winds for these IOPs is also seen in the geometric mean bias (MG), the ratio between model wind speed and observed wind speed. Only for IOP 8 is the MG greater than 1. The scatter in the dataset of wind speed observations, as expressed in the NMSE, is small, less than 0.61 (with an average of 0.37 for the NEV simulations). The scaled angle differences are typically less than 30° and for IOP 2 are less than 20°. Last, at least 50% (or 72% on average) of FEM3MP simulated data agree with reported observations within a factor of 2 (FAC2).

Data from two NEV simulations of IOP 2 are reported here because of changing wind direction during the simulation period; according to these metrics based on near-field wind flow measurements, the IOP 2 simulation with a wind direction of 215° is more consistent with observations based on all the statistical metrics discussed above.

Last, the flow results using the simpler LEV turbulence parameterization (which accounts for convection) are surprisingly good, often slightly better than the results using the more sophisticated NEV turbulence parameterization (which assumes neutral stability). This subtle improvement in wind field simulation cannot be construed as an endorsement of the LEV turbulence parameterization. As shown below, the NEV turbulence parameterization performs consistently better in predicting the concentration field in all cases, likely because the more sophisticated building-induced turbulence parameterization allows for more enhanced plume broadening as well as upwind transport of tracer concentrations.

#### 2) Concentration field

The predicted concentration patterns are also compared with observed 15-min averages of SF6 concentrations from Blue Boxes, of which 10–15 were deployed during each IOP. The locations of the Blue Boxes changed for each IOP, depending on the release location and expectations for vertical transport of SF6. The statistical performance measures of the simulations, in comparison with the observations, appear in Table 1.

For all IOPs, FEM3MP with the NEV turbulence parameterization could predict the significant upwind and lateral spread of material as indicated by the measured data. The fractional bias metric does indicate a systematic overprediction of concentrations. One factor contributing to this overprediction could be the difference between the steady inflow conditions used by the model, as compared with the unsteadiness of the actual winds. This variability would cause a plume to meander, thus reducing the peak values of concentration. Another possible influence could be the averaging time of the data (15 min), which is relatively long in comparison with the duration of the release (30 min) and likely does not reflect a “steady state” result as obtained by the model. Closer inspection of the differences between observed and predicted values indicates no clear pattern to this overprediction (not shown).

To delineate the stability effects on the dispersion results for these releases, simulations with the LEV parameterization, which includes stability effects but lacks a sophisticated parameterization for building-induced turbulence, were also performed; the resulting metrics for the LEV parameterization also appear in Table 1. Regardless of the metric used to compare observed concentrations with predicted concentrations, the NEV parameterization (which includes more sophisticated building-induced turbulence) performs much better than the LEV parameterization for both daytime and nighttime releases. These results suggest that building-induced turbulence dominates over any turbulence reduction due to the slightly stable conditions at night or turbulence increase due to convective conditions during the day, thus making the assumption of neutral stability in model simulations justifiable for long-duration releases under moderate wind conditions within a built-up urban area.

### c. Comparison of simulations with observations in the urban shadow

Wind speed, direction, and turbulence data collected on the crane pseudotower (in the urban shadow region, within 1000 m of the downtown area) are analyzed to construct winds and TKE profiles in the area for comparison with the FEM3MP simulations with the NEV turbulence parameterization. These results are considered together to assess the stability conditions associated with these simulations in the urban shadow, as well as to describe the fidelity of the model’s RANS simulations to the observed atmosphere.

In Fig. 2, a comparison of predicted versus observed profiles of four variables at the crane station for IOP 2 (simulated with an incoming wind direction of 205°) is presented. Included in the comparison are profiles of wind speed, wind direction, friction velocity, and TKE. As is seen in the figure, there is good agreement between the predicted and observed profiles for the wind speed (upper left panel). The FEM3MP simulation slightly underpredicts wind speed relative to all but the top level of the crane, and FEM3MP predicts wind directions that deviate from those measured at the crane by less than 10°.

The friction velocity *u*_{*} may be calculated two different ways from the crane data: first, by summing the momentum fluxes at each level independently (the “local” method), as in

and second, by assuming that the profile of wind speeds fit the relationship

where *k* is the von Kármán constant of 0.4, *U* is the mean wind speed at each height *z*, and surface roughness *z*_{0} is set to 0.5 m in the urban wake to achive optimal agreement between the local and profile methods. [Note that Gouveia et al. (2007) use a least squares fit method based on one month of data from the crane pseudotower and find *z*_{0} to be at least 0.5 for the wind directions considered here.] Because only this second calculation may be applied to FEM3MP output, only this “profile” method is presented in Figs. 2 –5. The simulated estimates of friction velocity are lower than the observed values in the 10–70-m levels, as is the case with the wind profiles.

In the lower-right panel, the TKE predicted by the model is shown, along with the total TKE observed at the crane site (line with diamond symbols). TKE is directly calculated by the NEV turbulence model in FEM3MP, and the data analysis method calculates TKE as 0.5( + + ). Also presented is the TKE profile with an estimate of the buoyant production contribution removed, against which the predicted profile should be compared. The two observed TKE profiles suggest that buoyant production contributes only between 5% and 25% of the total TKE budget for this IOP, one of the two daytime IOPs in which buoyant forcing is expected to be a significant contributor to the total turbulence budget.

In the lower-right panel of Fig. 3, a similar profile of simulated, observed, and observed minus buoyant TKE is presented for IOP 3. For this case, the shapes of the TKE profiles are very similar, though the model underpredicts the amount of TKE observed at the actual crane site. The estimate of buoyant contribution to the total observed TKE is much smaller for this daytime IOP, less than 10% at maximum. The crane observations do indicate higher wind speeds at the highest levels of the crane for IOP 3 than for IOP2, thus suggesting increased shear and therefore increased mechanical production of turbulence in the CBD area. As in the case of IOP 2, the simulations underpredict wind speed and friction velocity at the crane location. The wind direction profile indicates that the simulation wind direction was likely off by 5°–10°. With the shift in wind direction and the observations that buoyant production of turbulence is an order of magnitude smaller than shear production of turbulence, we assume that the omission of buoyancy effects in the NEV turbulence parameterization is not responsible for the deviations from the observations.

The two daytime IOPs are expected to depict the strongest role of buoyant forcing, while nighttime IOPs are expected to be more classically neutral. In Fig. 4, a comparison of predicted versus observed profiles of wind speed, wind direction, friction velocity, and TKE, at the crane pseudotower for nocturnal IOP 8 is presented. Although the simulation again underpredicts wind speed, friction velocity, and TKE, the shapes of the profiles are generally consistent with observations. Model predictions for both wind speed and friction velocity at the crane top are only approximately 65% of the observed values. The discrepancies between observations and simulations for IOP 8 could be due to larger-scale flow processes not currently accounted for in our simulation, particularly the vertical transport of momentum and turbulence. Lundquist and Mirocha (2006, 2007) discuss the occurrence of a nocturnal low-level jet on the night of IOP 8 based on data from the pseudotower and the PNNL boundary layer wind profiler 2 km SSW of the OKC urban core. Shear generated by the low-level jet (seen in Fig. 15 of Lundquist and Mirocha 2006) could be responsible for the vertical transport of momentum from upper levels into the lower levels simulated here, thereby increasing wind speed, friction velocity, and TKE. Because this present simulation excludes the possibility of vertical transport of momentum from outside the simulation domain, such processes are not simulated here and therefore could explain some of the discrepancies between observations and simulations.

In the upper-right panel, the predicted and observed profiles for the wind direction are compared. The agreement between predicted and observed profiles is generally good except for *z* = 15–30 m, with the largest underprediction of the angle by ∼10° near *z* = 15 m AGL. In the lower-right panel, profiles of observed TKE, observed TKE without buoyant production, and predicted TKE are displayed. The two observed TKE profiles suggest that buoyant production, in this case, is negligible. Also, as seen in Table 2, observations indicate the rate of buoyant production of turbulence is 2 orders of magnitude smaller than shear production of turbulence (and observations show buoyant production of turbulence rather than suppression). The predicted and observed profiles have similar shapes, although the predicted values are only about 55% of those observed. Again, the possibility of a nocturnal jet present during the release is a plausible explanation for the higher TKE observed. Most important in the context of this study, buoyant forcing plays no consequential role, either in production or suppression of turbulence in the crane observations, underscoring the suitability of an assumption of neutral stability for urban nocturnal simulations of long-duration releases.

Last, Fig. 5 compares predicted and observed profiles for IOP 9. The wind speeds, as predicted by the model, are slight overpredictions relative to observations, while the wind direction is as much as 20° away from observed wind directions. Profiles of friction velocity are slightly greater than observed, which is consistent with the slightly elevated wind speeds observed. The TKE profiles predicted by the model agree very well with those observed at the crane location. As with IOP 8, the role of buoyancy is minimal during this nocturnal release, with buoyant effects creating a small amount of turbulence, an order of magnitude less than shear effects. The simulation of this IOP agrees well with the observations, perhaps because of the fact that, as discussed in Lundquist and Mirocha (2006, 2007), no vertical transport of momentum from outside the simulation domain affects the flow.

### d. Observed TKE budgets in the urban shadow

TKE budgets are estimated from the crane pseudotower data, with dissipation, shear production, buoyant production or destruction, and turbulent transport calculated directly. The residual term includes advected turbulence, pressure transport, and the accumulated errors from the directly calculated terms. The TKE budgets corresponding to the releases simulated from IOPs 2, 3, 8, and 9 appear in Fig. 6, presented dimensionally. The largest directly calculated term for all IOPs is the dissipation term, which is consistent with expectations for a site immediately downwind of a large source of TKE production. Values for dissipation rate range from 0.004 to 0.09 m^{2} s^{−3}, well within the range observed in disturbed boundary layers (Piper and Lundquist 2004). Shear production tends to be the next largest term, although for IOP 3, turbulent transport is very large at the lower levels of the crane. Last, buoyant production or destruction of TKE is small even for the daytime IOPs (IOP 2 and IOP 3), which implies that buoyant forcing is not a critical consideration for simulations of long-duration releases in the shadow of an urban area. The storage of TKE is 1–2 orders of magnitude smaller than the shear production term.

The shapes of the budget profiles merits comment. For all IOPs, shear production increases with height, as expected for canopy flows. Dissipation rates are higher near the surface than aloft for the daytime IOPs, which exhibit moderate average wind speeds and flow from the south-southwest. For the nighttime IOPs, with lighter average wind speeds and flow from the south-southeast, dissipation rates increase with height. For all IOPs, the turbulent transport of TKE tends to peak around 30 m, which is approximately the average height of buildings in the downtown area. The residual term, which includes the advection term as well as any pressure transport and the accumulated errors of the directly calculated terms, mirrors largely the dissipation term. Because of the large dissipative term, we hypothesize, in view of the average TKE values discussed in the next subsection, that most of the observed TKE is generated within the urban matrix and advected downwind to be dissipated rather than advected in from outside the computational domain. A complete test of that hypothesis would require two sets of similar observations, both upwind and far downwind of the urban domain, and we encourage future field experiments to deploy such instrumentation if possible.

For each IOP, the portion of TKE due to buoyancy effects is estimated by multiplying the buoyant production rate *B* by a turbulent time scale *τ*, which is given by the total TKE over the dissipation rate ɛ. Values averaged over all crane levels for each IOP are summarized in Table 2. The turbulent time scales are all short, between 1 and 2 min, likely because of the small size of turbulent eddies at the crane pseudotower just downwind of the central business district in the urban shadow. (Recall that for convective boundary layers over homogeneous fetch, a turbulent time scale is typically calculated to be on the order of 10–20 min, for a boundary layer 1–2 km deep with convective eddies rising at 1–2 m s^{−1}.) The consistency between estimates of the turbulent time scale, for different IOPs, does suggest the possibility of coherent structures downwind of the urban area. In summary, only for IOP 2, which had a lower wind speed than the other daytime IOP, does a portion of TKE due to buoyant forcing become significant. At this location in the urban shadow, the rate of buoyant production of turbulence is always at least an order of magnitude smaller than the shear production of turbulence. It is also interesting to note that in all IOPs buoyancy effects are always observed to produce turbulence, rather than suppress it, even during the nighttime IOPs.

### e. Variability of TKE in the upwind urban–urban shadow–downwind wake regions

As argued above, neutral stability is often presumed in urban areas because of the enhanced mechanical mixing found in complexes of building geometries. In all of these simulations, significant TKE due to building-induced turbulence was generated. Because of the generation of strong turbulent shear stresses as a result of the flow impinging on building walls and flow separation from building edges, regions with the highest turbulence intensity are generally on the edges of the CBD in close alignment with the streamwise direction, as seen in Fig. 7, which depicts TKE contours at 32 m above the surface (average building height) for all four IOPs, focusing on the central business district. The regions of highest turbulence intensity coincide with regions of large horizontal shear, and so most of the TKE production in these regions can therefore be attributed to shear production.

Averaged TKE, integrated across the computational domain and for the lowest 200 m, offers insight into the generation and decay of mechanically produced TKE. These integrated values are presented in Fig. 8 as a function of north–south distance for all IOPs, which loosely approximates distance upwind and downwind. For comparison purposes, four zones may be defined, somewhat arbitrarily. The zones are defined as follows: upwind (*y* = from −400 to 0 m, in Fig. 1), urban CBD (*y* = 200–1000 m), downwind urban shadow (*y* = 1000–2000 m), and downwind wake (*y* = 2000–2600 m). (The location of the crane-observing tower, used above, is at *y* = 1200 m, within the urban shadow region.) Note that only the urban CBD zone appears in its entirety in Fig. 7. The same crosswind extent of 900 m, centered on the intersection of Reno and Robinson, and the same vertical extent, from the surface to 200 m, are used in all four zones. The averaged TKE values are presented in Table 3 and appear in Fig. 8 as horizontal dotted lines. The crane pseudotower is located within the urban shadow; values observed at that location also appear in Fig. 8 as asterisks. Substantial deviation between the crane observation at one location and the average across the east–west band is expected.

These values indicate a significant increase in turbulence intensity due to increased shear production of TKE in the urban CBD zone for all IOPs. The turbulence intensity in the urban area is at least 5 times that at upwind and is about three times (with the exception of a value of 1.6 for IOP 9, which will be discussed later) that in the urban wake. The urban shadow region, or the transition between the CBD and the downtown urban wake region, exhibits an exponential decay of turbulence that seems to be a function of wind speed, and so mean values are not particularly descriptive for this region in which turbulence generated within the CBD dissipates as it is advected downwind. Because of turbulence transport from the urban area, turbulence intensity in the far downwind wake area is slightly higher than that in the upwind region. Assuming the averaged TKE value at upwind is representative of the turbulence level in the rural area, the above estimates suggest that building-induced mechanical stresses have caused the turbulence intensity to increase by as much as seven times in the urban area. The two sets of values (for IOP 2) also indicate the incoming wind flow direction is not a strong factor in determining the averaged TKE values in different zones, at least for the slightly varying wind directions considered herein for Oklahoma City.

The ratio between CBD turbulence and turbulence in the urban wake exhibits a slight trend with incoming wind speed: IOP 9, with the highest incoming wind speed of 7.2 m s^{−1}, has the lowest ratio of CBD to wake turbulence. Because CBD-generated turbulence is advected downwind faster in this simulation, as compared with the other simulations, it is conceivable that less dissipation of CBD turbulence has taken place by the time that turbulence enters the “urban wake” zone.

The information in Fig. 8, particularly the length of the urban shadow exponential decay, could be used to indicate the regions in which neutral stability may be assumed appropriately. However, this calculation is highly dependent on wind speed and somewhat on wind direction, and should therefore be pursued in the context of creating a reference for locations of interest rather than general guidelines for all urban areas. The nature of the curves in Fig. 8 are indicative of turbulence generated in Oklahoma City and may not be representative of other locations, especially those with different urban geometries.

## 5. Conclusions

In this paper, the validity of omitting stability considerations when simulating long-duration transport and dispersion in the urban environment is explored using observations from the Oklahoma City Joint Urban 2003 field experiment and computations fluid dynamics simulations of that experiment. Four IOPs, two daytime (IOP 2 and IOP 3) and two nighttime (IOP 8 and IOP 9), have been simulated using the building-resolving FEM3MP computational fluid dynamics model to solve the Reynolds-averaged Navier–Stokes equations using two possible turbulence parameterizations and a mixed virtual–explicit building representation. Model performance metrics are calculated for the downtown area based on wind fields and tracer concentrations; in the urban shadow, model wind and turbulence fields are compared with observed profiles. Model predictions for all four IOPs, regarding winds, concentrations, profiles of wind speed, wind direction, and friction velocity, are generally consistent with and compare reasonably well to the field observations. Simulations using the NEV turbulence parameterization, which omits stability considerations but includes complex building-induced turbulence production, generally agree better with observations than simulations that do account for stability at the expense of a less-complex mechanical turbulence model, although neither set of simulations accounts for forced convection. The shapes of the predicted TKE profiles are similar to those observed in the urban shadow for all cases. However, the predicted turbulence intensities are in general too low, except for the case of IOP 9, for which wind speeds were somewhat overpredicted as well.

Inspection of the observations, including the TKE budget as measured at the crane pseudotower in the urban shadow, indicates a minimal role of buoyant forcing. The largest term of the TKE budget is the dissipation rate, which is balanced by a residual term that includes the advection of turbulence from upwind as well as pressure transport of TKE and errors in other calculations. Considering the large urban matrix found upwind of the observing station, and that significant turbulence is expected to be generated by flow through that area, we assume that much of the residual is attributable to advection of TKE from the downtown CBD. The buoyant production of TKE is always found to be positive, even at night, but typically found to be at least an order of magnitude less than the shear production of TKE. We conclude that buoyancy plays a minimal role in the production of TKE even in the urban shadow. The FEM3MP simulations’ deviations from the observations therefore do not seem to be related to “missing” buoyant forcing due to the assumption of neutral stability in the NEV turbulence parameterization. These deviations are most likely due to other sources including inaccuracies in boundary conditions (e.g., time-varying inflow conditions), inadequate representation of large-scale forcing (e.g., low-level jets, in the case of IOP 8), and to a lesser extent the possible limitations of the turbulence parameterizations used in this study.

Turbulent time scales are calculated to quantify further the importance of buoyant forcing relative to the total TKE. For IOP 2 (a daytime release with low wind speeds), the buoyancy production contributes no more than 25% to the total TKE budget in the region. Buoyancy plays less of a role in the other daytime release, IOP 3, probably because of the stronger winds in that case and the increased role of mechanical forcing. For both IOP 8 and IOP 9, nighttime releases, the contribution of buoyancy production to the total TKE budget in the same region is negligible.

Estimates of the average TKE values in various zones of the computational domain suggest that building-induced turbulence can cause the average turbulence intensity in the urban area to increase by as much as 7 times typical upwind values. This increase persists in the urban shadow but decays exponentially, relaxing by approximately 1 km downwind. In a wake region defined to be 1.5 km downwind and beyond, simulations show most of the urban-produced TKE has decayed. The exact downwind distance required for decay seems to be a function of wind speed for two reasons. First, higher wind speeds involve more shear and more production of turbulence, and second, higher wind speeds advect turbulence downwind faster. The dissipation of a given amount of TKE would require more distance to occur in a case with higher wind speed. Average TKE values in the upwind, downtown, and urban wake regions appear independent of the wind directions considered in this study.

Considering that buoyancy effects are negligible during the nighttime releases (IOP 8 and IOP 9), the assumption of neutral stability for nighttime simulations in the downtown and downtown shadow is acceptable. For one daytime release (IOP 2), although the buoyancy production could contribute up to 25% of the TKE budget, the assumption of neutral stability in the downtown and downtown shadow is still valid because building-induced turbulence is the dominant source of TKE in the area. Because of higher wind speeds in the other daytime release (IOP 3), buoyancy effects are minimal and can be neglected. However, in areas approximately 1.5 km and farther downwind, the levels of building-induced turbulence subside rapidly. Therefore, in these downwind regions, the assumption of neutral stability is less valid. Considering that building-induced mechanical turbulence plays a subdued role in these regions, it is likely that buoyancy effects should be considered in flow and dispersion simulations for those regions, at least for cities like Oklahoma City with concentrated central urban areas.

These conclusions are based on RANS simulations and relatively long releases (30 min) of tracer gas, which average out the effects of localized buoyant plumes and differential heating. Those localized buoyant effects, which were observed even at night in Oklahoma City, could become consequential for shorter-duration or “instantaneous” releases even within the downtown area.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory, under Contract W-7405-Eng-48. This document’s tracking number is UCRL-JRNL-217489. The authors appreciate insights from discussions with B. Kosović.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* J. K. Lundquist, Lawrence Livermore National Laboratory, P.O. Box 808, L-103, Livermore, CA 94551. Email: lundquist1@llnl.gov