Under the sponsorship of the U.S. Department of Energy and U.S. Department of Homeland Security, a computational fluid dynamics (CFD) model for simulating airflow and dispersion of chemical/biological agents released in urban areas has recently been developed. This model, the Finite Element Model in 3-Dimensions and Massively Parallelized (FEM3MP), is based on solving the three-dimensional, time-dependent Navier–Stokes equations with appropriate physics submodels on massively parallel computer platforms. It employs finite-element discretization for effective treatment of complex geometries and a semi-implicit projection scheme for efficient time integration. A simplified CFD approach, using both explicitly resolved and virtual buildings, was implemented to improve further the model’s efficiency. Results from our model are continuously being verified against measured data from wind-tunnel and field studies. Herein, this model is further evaluated using observed data from intensive operation periods (IOP) 3 and 9 of the Joint Urban 2003 field study conducted in Oklahoma City, Oklahoma, in July 2003. The model simulations of wind and concentration fields in the near and intermediate regions, as well as profiles of wind speed, wind direction, friction velocity, and turbulent kinetic energy (TKE) in the urban wake region, are generally consistent with and compared reasonably well to field observations. In addition, this model was able to reproduce the observed split plume of IOP 3 and the end vortices along Park Avenue in IOP 9. The dispersion results and TKE profiles at the crane station indicate that the effects of convective mixing are relatively important for the daytime release of IOP 3 but that the stable effects are relatively unimportant for the nighttime release of IOP 9. Results of this study also suggest that the simplified CFD approach implemented in FEM3MP can be a cost-effective tool for simulating urban dispersion problems.
Urban areas are the mostly likely locations for atmospheric releases of hazardous material, whether resulting from industrial accidents or from terrorist acts. To protect the population effectively, there is a great need for observational and modeling tools to track and to forecast the transport and dispersion of the hazardous material from such releases. The need for robust modeling tools, among others, was stressed in a recent report by the Committee on the Atmospheric Dispersion of Hazardous Material Releases (2003). Among the recommended priorities for improving modeling capabilities, the report states “New dispersion modeling constructs need to be further explored and possibly adapted for operational use in urban settings. This includes advanced, short execution time models, slower but more accurate computational fluid dynamics and large-eddy simulation models, and models with adaptive grids.”
Under the sponsorship of the U.S. Department of Energy (DOE) and U.S. Department of Homeland Security (DHS), we have recently developed a computational fluid dynamics (CFD) model for simulating airflow and dispersion of chemical/biological agents released in the urban environment. Our model, the Finite Element Model in 3-Dimensions and Massively Parallelized (FEM3MP), 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 with the Reynolds-averaged Navier–Stokes (RANS) and large-eddy simulation (LES) approaches (Chan and Stevens 2000), aerosols (Chan et al. 1999), UV radiation decay, surface energy budgets (Lee and Brown 2001), and vegetative canopies (Chan et al. 2002).
Predictions from our model are continuously being verified against measured data from wind-tunnel and field studies. Some of the model evaluation studies using wind-tunnel data were reported in Chan and Stevens (2000) and Chan et al. (2001). Examples of model evaluations using observed data from various field experiments were documented in Calhoun et al. (2004, 2005), Chan et al. (2003), Lee et al. (2002), and Humphreys et al. (2003).
In addition to model evaluation studies, our model has also been used to investigate the effects of inflow turbulence on dispersion scenarios involving nighttime releases under light and highly variable winds, such as the case of intensive operation period (IOP) 7 of the Urban 2000 experiment (Allwine et al. 2002). Through a series of controlled numerical experiments with variations in time-dependent forcing and turbulence intensity from the inflow boundary, Chan and Leach (2004) demonstrated that to simulate urban dispersion scenarios successfully under light and highly variable winds it is necessary to use appropriate time-dependent forcing and turbulence from the larger-scale flow through the inflow boundary. Their results also indicate that inflow turbulence is as important, if not more so, than building-induced mechanical turbulence in dispersion scenarios under the above conditions.
Although 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 unsuitable for emergency-response applications. Our goal is to develop a sufficiently fast CFD urban dispersion model for integration into the DOE National Atmospheric Release Advisory Center operational modeling system. As a step toward reaching such a goal, a simplified CFD approach, with options to use either virtual buildings (represented as large drag forces) only or a combination of explicitly resolved and virtual buildings, has been developed (Chan et al. 2004). With the latter option, only targeted and important buildings are explicitly resolved using fine grid resolution, and the remaining buildings are modeled as virtual buildings using coarser grid resolution. The drag forces were modeled as a sink term in the mean momentum equations, similar to the canopy-drag term proposed by Yamada (1982) and Brown and Williams (1998).
In this paper, the FEM3MP model is further evaluated in a more detailed comparison with a subset of the data collected in the Joint Urban 2003 field study: our model using both explicitly resolved and virtual buildings was used to simulate the first continuous release of IOPs 3 and 9, and simulation results of wind and concentration in the near and intermediate regions (from the source locations), as well as profiles of wind speed, wind direction, friction velocity, and turbulent kinetic energy (TKE) at the pseudotower location in the urban wake region, are compared with the field observations. In the following, we first give a brief description of the simplified CFD approach implemented in FEM3MP and then discuss briefly the field experiments being simulated, present a comparison between model results and observed data, and then offer a few concluding remarks.
2. A simplified CFD approach
For convenience of code parallelization and computational speed, our current version of FEM3MP employs a structured mesh (but a graded and distorted mesh is allowed), and buildings within the computational domain are represented as solid blocks with velocity, pressure, and diffusivities set equal to 0. The crux of our simplified CFD approach (Chan et al. 2004) is to resolve targeted and important buildings explicitly with fine grid resolution and to treat the remaining buildings as virtual buildings with coarser grid resolution. An initial evaluation of the approach, using data from the Joint Urban 2003 field experiment (Allwine et al. 2004), can be found in Humphreys et al. (2004).
The virtual buildings are modeled as drag forces in the mean momentum equations by using a term that is similar (but not identical) to the canopy-drag term proposed by Yamada (1982) and Brown and Williams (1998): at grid points associated with the virtual buildings, a term in the form of Cd|U|Ui, where Cd is a drag coefficient (m−1), |U| is the local wind speed, and Ui is the ith velocity component, is added to the mean momentum equations. By using a large value of the drag coefficient, as will be seen in the example below, a virtual building can be made to act like an explicitly resolved (solid) building. The drag term is linearized and treated implicitly in the time-stepping algorithm. In addition, the diffusivities at the corresponding grid points are set equal to the molecular diffusivity of air, instead of zero, to avoid numerical instability.
In addition to computational efficiency, the virtual-building approach also offers the advantages of much easier grid generation and a direct means for evaluating the drag forces to aid the parameterization of the urban canopy in larger-scale CFD models. However, predicted results in the vicinity of the virtual buildings are, as expected, slightly less accurate, because coarser grids are used and small values of velocity and diffusivities (instead of zero) are present as a result of the drag force approximation.
The following simple example, using a relatively coarse grid, shows how the virtual-building approach performs in simulating the airflow and tracer dispersion around a cubical building of unit dimension (H = 1). Four simulations were performed. The cubical building was explicitly resolved (as a solid building) in one simulation and modeled as a virtual building (drag forces) with a drag coefficient of 15, 50, and 100, respectively in the remaining simulations. A computational domain of 8H × 6H × 2H in the longitudinal, lateral, and vertical directions, with a graded mesh consisting of 43 × 33 × 15 grid points was used in all simulations. For simplicity, a logarithmic velocity profile, with friction velocity u* = 0.0356 m s−1, roughness length z0 = 0.001 m, and U = 0.6 m s−1 at H, was specified on the inlet plane, and a similarity K-theory model (Dyer 1974) was used for turbulence parameterization. In each case, a steady-state wind field was established first (after 60 s of simulated time) and then followed by the dispersion simulation associated with a ground-level tracer released continuously over an area of 2 cells at 1.5H in front of the cube for a duration of 120 s.
In Fig. 1, the predicted steady-state velocity and normalized concentration on the vertical plane of symmetry with the cube modeled as a virtual building (top three panels) and as a solid building (bottom panel) are shown. The normalized concentration is defined as X = CUH2/Q, in which C is the calculated concentration (mass fraction), U is the reference velocity at H, and Q is the source rate. Predicted wind vectors obtained from the virtual-building simulations indicate that, as the drag coefficient is increased to 100, the wind speeds and airflow patterns, including flow separation on the rooftop, a small eddy and the indication of a dividing streamline on the windward side, and a large eddy on the leeward side of the cubical building, become very close to those obtained from the simulation using the solid-building approach. A similar trend of convergence is observed for the concentration field.
In Fig. 2, the corresponding velocity vectors and concentration field on the plane z/H = 0.2 obtained from the four simulations are compared. Again, as the drag coefficient is increased to 100, results from the virtual-building approach become very close to those obtained from the run using the solid-building approach. The similar features include a diverging flow in the front, separation on the sides, and a reverse flow with two counterrotating vortices in the wake region, as well as the horseshoe-shaped concentration patterns on the plane. These results suggest that a drag coefficient of O(100) is a reasonable parameter to use.
As alluded to earlier, because of the approximate nature of the predicted flow in the vicinity of the virtual building, less accurate predictions, such as spurious infiltration of the tracer (of lower concentrations) into the virtual building, are expected. It must be emphasized that the main idea of the simplified CFD approach is to use virtual buildings judicially in places where large gradients of the concentration field are absent and also at locations where only reasonably accurate results are warranted. In these respects, the example presented here is a severe test of the virtual-building approach, and the results from this example are considered to be reasonable. Explicitly resolved buildings should always be used wherever highly accurate predictions are important.
Using four 2.4-GHz “Xeon” processors, the central processing unit (CPU) times required to generate the steady-state flow field by the runs using the virtual-building approach (with Cd = 15, 50, and 100) are 131, 133, and 134 s, respectively, and the solid-building approach requires 136 s. Total CPU times required for the four dispersion simulations are basically the same. The savings in computational cost by using the virtual-building approach, although insignificant, are nevertheless noticeable. For practical urban dispersion simulations, which involve typically hundreds or more buildings, the savings in computational cost by modeling a large fraction of the buildings as virtual buildings can be significant. For instance, in simulating a hypothetical tracer gas release in the Salt Lake City, Utah, downtown area involving O(100) buildings with an all-virtual-building approach, together with a coarser grid (than the one used in an all-solid-building approach), Chan et al. (2004) were able to reproduce reasonably well the main features of concentration patterns predicted by the all-solid-building approach, but with an order-of-magnitude savings in computational cost.
The option of using only virtual buildings has been developed with a certain category of emergency-response applications in mind, for which airflow in street canyons still needs to be fairly well resolved but a compromise between accuracy and computational speed is definitely required. This type of application may become feasible with a simplified CFD approach using only virtual buildings, together with a simple turbulence model, a relatively coarse and graded grid, and modest computer resources. For these reasons, the value of the drag coefficient was made based on a simple K-theory turbulence model. The value of the drag coefficient should perhaps have been calibrated for the nonlinear eddy viscosity (NEV) turbulence model as well; however, computational time could increase greatly because of the necessity of solving three additional (turbulence) equations and the need for finer grid resolution to realize fully the higher accuracy offered by the NEV turbulence model. Therefore, the combination of virtual buildings and NEV turbulence model has not been considered for emergency-response applications. Nevertheless, some of our recent urban dispersion simulations have been made using virtual buildings exclusively, together with either linear eddy viscosity (LEV) or NEV turbulence models. Results from these simulations seem to indicate the same value of drag coefficient is appropriate for both turbulence models.
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 the Defense Threat Reduction Agency of the U.S. Department of Defense cosponsored a series of dispersion experiments, named Joint Urban 2003 (JU2003), in Oklahoma City (OKC), Oklahoma, during July of 2003 (Allwine et al. 2004). These experiments are complementary to the Urban 2000 field study (Allwine et al. 2002) in that they provide another comprehensive field dataset for the evaluation of CFD and other dispersion models. In contrast to the Urban 2000 experiments, which were conducted entirely at night, these experiments took place during daytime and nighttime to include both convective and stable atmospheric conditions. Prior to the field study, FEM3MP simulations were performed to provide guidance for the selection of release sites and the deployment of wind and concentration sensors.
A total of 10 IOPs were conducted, and sulfur hexafluoride (SF6) in the form of puffs or continuous sources was released over six daytime and four nighttime episodes. Many wind and concentration sensors were used to collect wind and SF6 data over both long and short time-averaging periods. In addition to surface measurements, wind and concentration profiles adjacent to the outside walls of several buildings were also taken. In some cases, balloons were deployed close to the tracer release area. Many of the released balloons exhibited quick ascents from ground level to the rooftop of buildings, implying highly convective conditions.
During the JU2003 experiment, a pseudotower, supported by a 90-m crane and fitted with sonic anemometers at eight levels, was also deployed for wind and turbulence observations. The pseudotower (crane station) was located in the urban wake region at approximately 750 m north-northwest from the OKC central business district. Winds, temperature, and turbulence data collected on the pseudotower have been analyzed by Lundquist and Chan (2007) to construct profiles of wind speed, wind direction, friction velocity, and TKE. These observed profiles are also used in this model evaluation study.
All of the data used in this study were downloaded from the JU2003 Internet site (https://ju2003-dpg.dpg.army.mil). The locations used were those reported in the README files or in the file metadata. All data locations that were reported in latitude–longitude were converted to universal transverse Mercator coordinates using standard conversion algorithms. The wind data were averaged to 30 min from the various sensors, either by averaging the original 10-Hz data or by averaging data that had previously been averaged over shorter time intervals. The tracer data used for comparison with model results had been obtained by using 15- or 30-min averaging time.
4. Model–data comparison
In this study, airflow and dispersion simulations associated with the first continuous release of IOPs 3 and 9, which are daytime and nighttime releases, respectively, were performed. In each case, SF6 was released near the ground as a point source for 30 min, with a release rate of 5.0 g s−1 for IOP 3 and 2.0 g s−1 for IOP 9. Shown in Fig. 3 are the footprints of buildings in the central business district of OKC, with the release locations indicated by S3 for IOP 3 and S9 for IOP 9. The tallest building in the area is approximately 120 m high and the average building height in the area is ∼30 m.
In the numerical simulations, a computational domain of 1030 m × 3010 m × 425 m (in the lateral, longitudinal, and vertical directions) was used. A graded mesh consisting of 201 × 303 × 45 grid points, with a minimal grid spacing of ∼1 m near the ground surface and certain explicitly resolved buildings, was used. Most of the buildings within ∼500 m of the release points were explicitly resolved, and the remaining buildings were modeled as virtual buildings.
Steady logarithmic velocity profiles were constructed from nearby sodar and weather station observations and used as the inflow conditions. The resulting wind speed at z = 50 m is 6.5 m s−1 for IOP 3 and 7.2 m s−1 for IOP 9; the estimated average wind direction is 185° for IOP 3 and 180° for IOP 9. A comparison of the inflow wind speed and direction profiles (solid lines) with sodar data (dashed lines) for the two IOPs is shown in Fig. 4, indicating that the constructed profiles are a reasonable representation of the observations. Based on the sodar and sonic anemometer data we examined, directional shear and vertical motion do not appear to be significant for the simulated releases. However, for IOP 8, because of the presence of a nocturnal low-level jet (Lundquist and Mirocha 2006), directional shear and vertical motion may be significant and may have to be considered in numerical simulations.
FEM3MP does have the capability to use a realistic profile (including wind shear and time variations) other than a logarithmic profile. However, appropriate measurements for defining detailed inflow profiles are not available; thus logarithmic profiles were constructed from nearby sodar and weather observations and were used in our simulations. The stable nocturnal conditions exhibited by the observed surface data in Fig. 4 do not appear to be representative of the nearly neutral atmospheric conditions observed elsewhere and are thus not incorporated into the inflow conditions. However, our model simulations (to be shown later) for winds and concentration field in OKC downtown and urban wake regions are generally consistent with observations, which seem to suggest that flow and concentration in the downtown and wake regions are only weakly affected by upwind surface winds (below the average building height). This is probably true, because winds below the average building height are often greatly altered by buildings in the urban area, in addition to the changes caused by the upwind fetch (which is about 600 m long in our domain of simulations). In general, appropriate inlet boundary conditions are important for accurate flow and dispersion simulations and should be further studied.
For IOP 3 (Fig. 4a), the actual wind and direction profiles are from a sodar deployed in the botanical gardens about 100 m to the south of the release point. For IOP 9 (Fig. 4b), the observed wind speed and direction profiles are from a sodar about 4 km downwind of OKC. We recognize that this is downwind but feel that it is more representative of the upwind conditions of the release, for which the crane data (to be shown later) indicate the atmospheric conditions for this nighttime release were nearly neutral. The sodar observations from the botanical gardens and another upwind site at the OKC maintenance yard exhibited characteristics of a nocturnal decoupling of the surface layer from the boundary layer above. Such stable nocturnal conditions do not appear to be representative of the nearly neutral atmospheric conditions observed elsewhere.
On the inflow boundary, the above velocity profile, together with values of TKE and epsilon (of the turbulence equations) that are consistent with the specified velocity profile, was imposed. On the walls of the explicitly resolved buildings and ground surface, no-slip boundary conditions (zero velocity) and zero values of TKE and epsilon were specified. No penetration (zero vertical velocity) was applied on the top boundary, and natural boundary conditions (zero normal and tangential stresses) were used on the remaining boundary.
For each simulated release, a quasi-steady-state flow field was first established after ∼10 min of simulated time before the tracer was released. The release of SF6 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 that the steady state was reached in about 20 min of simulated time. For both cases, the RANS approach with an NEV turbulence model (Gresho and Chan 1998) was used and neutral atmospheric stability was assumed.
In the following, model predictions of flow and concentration in the near and intermediate regions of the release point are presented and are compared with observed data. Several of the statistical performance measures recommended by Hanna et al. (2004) are used to assess the performance of our model. They are the fraction of predictions within a factor of 2 or 5 (FAC2 or FAC5, respectively), fractional bias (FB), geometric mean bias (MG), and normalized mean square error (NMSE). For differences in angles between predicted and measured velocity vectors, the formula of scaled average angle differences (SAA) devised by Calhoun et al. (2004), with larger vectors carrying more weights, is also used.
The equations for these metrics are defined as
in which C is the data type being evaluated (e.g., wind speed or tracer concentration), Cp is the model result, and Co is the observation, with overbars denoting averages. In the above metrics, FB and MG measure the systematic bias of a model in terms of differences and ratios and NMSE measures the scatter associated with the model output relative to observations. A perfect model would have FACx = 1.0, FB = 0, MG = 1.0, and NMSE = 0.
The SAA, a model performance metric calculated from wind speeds and wind directions, is given by
in which ϕ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 |Ui| 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.
Additionally predicted and observed profiles of wind speed, wind direction, friction velocity, and TKE at the crane station (at about 750 m north-northwest from the OKC central business district) are compared. The friction velocity u* is calculated by assuming that the profile of wind speeds fit the relationship
in which k is the von Kármán constant of 0.4, U is the mean wind speed at each height z, and surface roughness z0 is set to 0.5 m in the urban wake.
a. IOP 3
Airflow in urban areas is extremely complex, as is illustrated in Fig. 5 for IOP 3. In the figure, predicted wind vectors and corresponding wind speeds (grayscale contours) near the ground (z = 2 m) are displayed. Some of the complex flow features include flow separations, stagnation zones, various sizes of eddies, and high-velocity jets in street canyons. Also, it is interesting to see how the flow separates at the southwest corner of the building in the center, which is obviously the cause of a split plume as will be shown later.
In Fig. 6, the predicted wind vectors (without arrowheads) in the downtown area are compared with the 30-min-averaged data (with arrowheads) measured by Dugway Proving Ground (DPG) portable weather information display system (PWIDS) and sonic anemometers on three towers along Park Avenue (y = 430–460 m). Overall, the agreement between model predictions and field observations is good. The values of statistical performance measures are SAA = 25.3, FAC2 = 0.74, FB = −0.21, MG = 0.79, and NMSE = 0.30.
There were several anemometers at various levels on each of the towers along Park Avenue. Only the level from the tower nearest to 8 m is used to compare with the simulations from FEM3MP. For this IOP, there is reasonably good agreement, especially for the two westernmost towers. For the tower near the eastern end of the Park Avenue, the observed wind vector is much smaller in magnitude, with close to a 90° direction error. This tower was near a bank of trees that was not represented in the FEM3MP simulations, and this smaller vector may represent vegetation-canopy effects.
In Fig. 7, a comparison of simulated and observed profiles of wind speed, wind direction, friction velocity, and TKE is presented. In general, there is a reasonable agreement between the simulated and observed profiles for wind speed and friction velocity, with the modeled values in the range of 60%–75% of the observations. The wind direction profiles indicate that the simulated wind direction is likely off by 5°–15°. In the lower-right panel, profiles of simulated, observed, and observed minus buoyant TKE are presented. The shapes of the TKE profiles are fairly similar, with the simulated TKE values in the range of 50%–90% of the observed values. As is seen in the figure, the estimated buoyant contribution to the total observed TKE is only 10% at maximum.
A discussion of the TKE budget and evaluation of TKE and its dissipation rate from the crane data can be found in Lundquist and Chan (2007). To estimate the contribution of buoyant forcing to the total TKE, the rate of buoyant production is multiplied by a turbulent time scale τ, which is determined from the quotient of TKE over the dissipation rate, following the model of Zeierman and Wolfshtein (1986). The estimated turbulent time scales for the two simulated experiments are 40–145 s for IOP 3 (a daytime release) and 70–85 s for IOP 9 (a nighttime release).
The simulated concentrations (solid line) along Broadway Avenue (at abscissa x = 115–150 m in Fig. 6) are compared with the time-averaged data (circles) in Fig. 8. The model results are from x = 140 m, and the field data are a collection of observations from sensors located near x = 140 m, with most of the data averaged over t = 15–30 min. The exceptions are values at downwind y = 670, 775, and 890 m, which were averaged over t = 0–30 min. The agreement is good and is within a factor of 2 in the urban area (downwind distance y < 600 m), beyond which the modeled values are much higher than observed. The discrepancies are mainly attributable to the following: the observed values at y = 670, 775, and 890 m are most likely too low because they were obtained with a 30-min averaging time (for a 30-min release), the gas sensors might have been too far apart (∼200 m at the 2000-m arc) to capture the higher concentration values, and the assumption of neutral atmospheric conditions in our current NEV turbulence model may not be appropriate in the urban wake region and beyond. To evaluate the stability effects, another simulation using an LEV turbulence model based on the similarity K theory (Dyer 1974), together with an estimated Monin–Obukhov length of −200 m, was performed. A better agreement between model results (dashed line) and data in the far field was indeed observed. However, the agreement in the near field is not as good because of a less sophisticated representation of turbulence near buildings.
The LEV model is actually a simplified version, that is, without heavy-gas effects, of the modified K-theory turbulence model developed by Chan et al. (1987) for simulating atmospheric dispersion of heavy gases over variable terrain. The model accounts for turbulence resulting from atmospheric stability but has no explicit treatment for building-induced turbulence. Thus, the model is more appropriate for upwind and regions farther away from the urban area. Despite it being simpler and less accurate than other turbulence models such as NEV and LES, the LEV turbulence model is capable of capturing the major features of airflow around a complex building (Calhoun et al. 2004) and has also been observed to perform reasonably well in one of our recent urban flow and dispersion studies (Lundquist and Chan 2007). On the other hand, the NEV model (Gresho and Chan 1998) has treated building-induced turbulence explicitly through solving three additional equations for turbulence and is more appropriate for urban flow and dispersion simulations in general. However, our version of the NEV model can treat only airflow under nearly neutral atmospheric conditions. Without including the stability effects, results can be inaccurate in regions in which atmospheric stability plays an important role, such as the urban-wake and farther-downwind regions.
In Fig. 9, modeled concentration patterns in the source area are compared with the measured data (small squares with the same color scheme). Except for missing or underpredicting a few very low concentrations near the left edge of the plume, the predicted concentrations generally agree well with the observations. In addition, the simulation was able to predict a split plume in front of the building nearby. The values of statistical performance measures are FAC5 = 0.42, FB = −0.56, MG = 6.2, and NMSE = 14. The values of MG and NMSE are high because of a bias produced by the presence of two high concentrations (one is near the edge of the plume and the other is to the south of the building) within a small sampling population. Hanna et al. (2004) pointed out that the values of MG and NMSE could be overly influenced by infrequently occurring high observed and/or simulated data. When the two pairs of highest concentrations were excluded in the performance evaluation, the values of statistical performance measures became FAC5 = 0.50, FB = −0.35, MG = 4.82, and NMSE = 0.61.
b. IOP 9
In this section, sample flow and dispersion results from simulations of the IOP 9 release are presented and are compared with available data. In Fig. 10, the modeled wind vectors and speeds (grayscale contours) in the source area are depicted to illustrate again the complexity of airflow in an urban area, including stagnation zones in front of the buildings, flow separations on the sides, jetting in street canyons, and large wakes behind buildings. In addition, there are two counterrotating vortices behind the wide building on the south side of Park Avenue (y = 430–460 m). Such end vortices were also reported by Brown et al. (2004) in the field study. Another interesting feature of the wind field is the strong reverse flow (northerly winds within the northwest quadrant of the figure) in front of the 105-m-tall Kerr-McGee building (near the north edge of the figure).
In Fig. 11, simulated wind vectors in the downtown area are compared with the 30-min-averaged data measured by DPG PWIDS and sonic anemometers on four towers along Park Avenue (y = 430–460 m). Again, the overall agreement between model results and field measurements is good, especially in the source area (middle of the figure). The statistical performance measures are similar to those of the previous case: SAA = 34.2, FAC2 = 0.71, FB = −0.20, MG = 0.98, and NMSE = 0.48.
In this IOP, the observed winds along Park Avenue are more variable. Again, only the wind vectors from the towers nearest the 8-m level along Park Avenue are compared. The vectors agree very well at the more southern tower in midblock and at the westernmost tower, where both the simulation and observation show very light wind speeds. The wind vectors do not agree well at the more northern midblock tower or at the tower near the east end of Park Avenue. The observations hint at the existence of a counterclockwise eddy in the eastern half of the block, whereas the simulation has such an eddy but it does not penetrate as far west in the urban canyon (see Fig. 10). Some of the discrepancies may be due to the use of the 180° wind direction for the inflow conditions.
A close comparison between Figs. 6 and 11 reveals that some of the wind vectors are considerably different in speed and direction, even though the inflow wind directions differ merely by 5° and the wind speeds differ only by ∼10%.
In Fig. 12, the modeled versus observed profiles of wind speed, wind direction, friction velocity, and TKE are compared. In general, the shapes of the predicted profiles are in good agreement with those observed. The wind speeds are slightly overpredicted relative to observations, and the wind direction is about 10°–15° away from observed wind directions. Modeled values of friction velocity are slightly greater than are observed, which is consistent with the slightly overpredicted wind speeds. The TKE profiles given by the model agree very well with those observed at the crane location, and the role of buoyancy is minimal during this nocturnal release. Again, the contribution of buoyant forcing to the total TKE was obtained by the product of the rate of buoyant production multiplied by a turbulent time scale τ, which is determined from the quotient of TKE over the dissipation rate, following the model of Zeierman and Wolfshtein (1986). The estimated turbulent time scale for this nighttime release is 70–85 s.
Because of the southerly ambient winds and a lower source rate (2.0 instead of 5.0 g s−1) released in the middle of Park Avenue, observed concentrations along Broadway Avenue were much lower than those observed in IOP 3. In addition, some of the data were not usable because of missing samples, field sampling problems, or problems during laboratory analysis. For these reasons, it was decided to compare the downwind concentration along the plume “centerline” instead.
In Fig. 13, modeled concentrations (solid line) along the plume centerline (at abscissa x = 32 m in Fig. 11) are compared with the time-averaged data (circles). The field data are a collection of observations from sensors located near x = 32 m, with most of the data averaged over t = 15–30 min. Two exceptions are the values at downwind y = 775 and 2370 m, which were averaged over t = 0–30 min. Despite considerable underpredictions in the near field (downwind distance y < 800 m), the overall agreement between model results and observations is within a factor of 2. To evaluate the stability effects, a simulation using the LEV turbulence model based on the similarity K theory (Dyer 1974), together with an estimated Monin–Obukhov length of 300 m, was performed. In this case, the simulation using the LEV model, together with parameterization for stable stratification, was unable to yield the higher concentrations expected for all locations, because the LEV model is not sophisticated enough to model the complexity of turbulence near buildings and in the wake region.
Unlike the daytime release of IOP 3, dispersion results in the intermediate region and beyond are only slightly affected by the slightly stable atmospheric stability. These results are consistent with the observed TKE data at the crane station (lower-right panel of Fig. 12), which indicates any turbulence reduction resulting from the slightly stable conditions at night is minimal. Considering this fact and the finding by Lundquist and Chan (2007) that building-induced turbulence is dominant in an urban area, it is justifiable to use the neutral stability assumption for this nighttime dispersion simulation.
The modeled concentrations in the source area are compared with measured data in Fig. 14. Again, except for a few locations near the western edge of the plume, the simulated results generally agree well with the observed data. In addition, the model was able to simulate successfully the unusually high concentration of over 500 ppb near the east end and south side of Park Avenue. The counterclockwise vortex near the east end of the street is believed to be mainly responsible for producing such a surprisingly high concentration at the location. As quantitative measures of model performance, the statistics are FAC5 = 0.56, FB = −0.39, MG = 2.0, and NMSE = 0.96.
As mentioned earlier, balloons released as visible tracers during some of the experiments exhibited quick ascents from ground level to the top of buildings, implying significant updrafts during those experiments. In Fig. 15, updrafts in street canyons and their effects on tracer dispersion are illustrated. In Fig. 15a, the simulated and observed concentration profiles on the outside walls of the building at the northeast corner of Park Avenue in Fig. 14 are compared. In general, our model predictions are able to reproduce a lofting plume observed in the field and also match fairly well with the measured values, mostly within a factor of 3. The largest overpredictions occur at the south and east sides of the building and are very likely due to the inaccurate inflow direction used in the simulation. The estimated inflow direction may be off by 10°–15°, as suggested by results in the upper-right panel of Fig. 12. Because there is a large open space with some lower buildings on the south and east sides of the building, specifying an inflow direction of 165°–170° (rather than 180°) would have resulted in more airflow to reduce the modeled concentrations on those sides of the building to make the model output agree better with the observed data. In Fig. 15b, the simulated concentrations and wind vectors on the plane cutting through the source (x = 32 m) are shown. The wind vectors show clearly the strong updrafts in building wakes and intense vortex motions around the buildings. As a result, the plume rises and reaches above the rooftops of certain buildings.
In this paper, FEM3MP has been evaluated using wind and concentration data obtained from IOPs 3 and 9 of the Joint Urban 2003 experiment. Our model simulations for the two IOPs, for both wind and concentration fields in the near and intermediate regions, as well as profiles of wind speed, wind direction, friction velocity, and TKE at the crane station, are generally consistent with and compared reasonably well to field observations. In addition, our model was able to reproduce the split plume observed in IOP 3 and the end vortices and an unusually high concentration near the east end of Park Avenue observed in IOP 9.
Judging from the crane data and modeled dispersion results, the effects of convective mixing in the daytime release (IOP 3) appear to be relatively important in the intermediate (urban wake) region and should be considered appropriately. On the other hand, for the nighttime release of IOP 9, the crane data indicate any turbulence reduction resulting from the slightly stable conditions at night is minimal. Considering this fact and the finding by Lundquist and Chan (2007) that building-induced turbulence is dominant in the urban area, it is justifiable to use the neutral stability assumption in the urban dispersion simulation for IOP 9.
The overall results of this study suggest that the simplified CFD approach implemented in FEM3MP, with explicitly resolved and virtual buildings, can be a cost-effective tool for simulating urban dispersion problems. We will further evaluate and improve our model, with the goal to produce a sufficiently fast CFD model for integration into the DOE National Atmospheric Release Advisory Center operational modeling system.
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. The authors thank Drs. David Stevens, Tom Humphreys, and Bob Lee for their contributions to FEM3MP and in the early phase of model evaluation studies using the Joint Urban 2003 data.
Corresponding author address: S. T. Chan, Lawrence Livermore National Laboratory, P.O. Box 808, L-103, Livermore, CA 94551. Email: firstname.lastname@example.org
This article included in the Urban 2003 Experiment (JU2003) special collection.