Abstract

Results are presented from a series of numerical studies designed to investigate the atmospheric boundary layer structure, ambient wind, and pollutant source location and their impacts on the wind field and pollutant distribution within the built-up areas of Shenyang, China. Two models, namely, Open Source Field Operation and Manipulation (OpenFOAM) software package and Weather Research and Forecasting (WRF) model, are used in the present study. Then the high resolution computational fluid dynamics (CFD) numerical experiments were performed under the typical simulated atmospheric boundary conditions. It was found that the atmospheric boundary structure played a crucial role in the pollution within the building cluster, which determined the potential turbulent diffusion ability of the atmospheric surface layer; the change of the ambient wind direction can significantly affect the dispersion pattern of pollutants, which was a more sensitive factor than the ambient wind speed; under a given atmospheric state, the location of the pollution sources would dramatically determine the pollution patterns within built-up areas. The WRF-CFD numerical evaluation is a reliable method to understand the complicated flow and dispersion within built-up areas.

1. Introduction

The air pollution within urban environments with dense population has become an important research issue in the past two decades, leading to numerous studies of the airflow and dispersion patterns around various building configurations. Field measurements and laboratory-scale physical modeling such as water-tank or wind-tunnel experiments are two traditional ways to study wind flow and dispersion process around buildings [18]. However, both these two methods are costly, which makes them difficult for practical application.

With rapid development in computer hardware and numerical algorithms, computational fluid dynamics (CFD) models have become common tools to simulating and predicting flow and pollutant dispersion patterns in real built-up areas. The CFD simulation method consists of solving the transport (advection and diffusion) equation of concentration based on the velocity field obtained from the Navier-Stokes equations. CFD prognostic models can be classified into three broad categories according to the respective turbulence closure schemes: Reynolds-averaged Navier–Stokes (RANS) approach [912], Large-eddy simulation (LES) [1315], and direct numerical simulation (DNS) approach. The choice among these three methods is a balance between the cost and the goal. In this study, the RANS model will be used because of its higher computational efficiency.

Furthermore, there are several studies that have numerically examined the air flow and dispersion in built-up areas by using CFD models with boundary conditions given by mesoscale models. Tewari et al. demonstrated that by using output of the WRF model as the initial and boundary conditions, the prediction ability of a CFD model employed over an urban area can be significantly improved [16]. Liu et al. investigated wind field and traffic pollutant dispersion at street level by coupling of LES in the local area with mesoscale atmospheric model WRF in the entire city [17]. Miao et al. coupled a CFD software package with a mesoscale weather model WRF to studying the airflow and dispersion of pollutant in a complex urban area of Beijing, China [18].

Shenyang located in Northeast China Plain is an important center of mineral and base of heavy industry of China. Pollution managements have been performed during recent years; however, air pollution remains a severe problem in Shenyang compared to other cities. Previous study performed by Miao et al. [18] made a probe into studying and predicting urban flow and dispersion in densely built-up areas using a CFD model OpenFOAM coupled with a mesoscale model WRF, and just one case was carried on under unstable atmospheric boundary condition in summer. In this study, a series of sensitivity numerical experiments have been performed to study the sensitivity of factors associated with flow and dispersion process within the built-up area of Shenyang with high resolution, such as the atmospheric boundary structure, ambient wind direction and speed, and the location of the pollution sources, using the coupled WRF-CFD model. Furthermore, an urban canopy parameterization scheme was used in WRF simulation to better simulate urban meteorological conditions. The study may support decision-making for pollution control strategies and traffic planning in Shenyang.

The remainder of the paper is organized as follows. In Section 2, the mesoscale model and CFD model were described. In Section 3, the setup of numerical experiments was given. And the results and discussions of the numerical experiments were then presented in Section 4, and finally the conclusions were given in Section 5.

2. Model Description and Coupling Method

2.1. Mesoscale Model

A community model, the WRF model, was used in this study as a mesoscale model. WRF is a nonhydrostatic mesoscale model with several options for dynamic cores, as well as various choices for physical parameterizations. In this study, we used the Advanced Research WRF Version 3.4, which was released in April 2012. The WRF model was used to simulate the wind and turbulence fields on urban scale, whose horizontal dimension is from tens of kilometers to hundreds of kilometers [19]. The WRF-simulated results can be used to provide the boundary conditions for the CFD model to recalculate flow fields for the subdomain scale simulation [18, 20], whose horizontal dimension is several kilometers [21].

2.2. CFD Model

As a CFD model, the Open Source Field Operation and Manipulation (OpenFOAM) software package version 2.1.1 which is a collection of C++ libraries was used to solve complex problems in fluid mechanics. OpenFOAM is completely free distributed and allows user to develop some specific solvers, which can be integrated with already existing tools. For further information about the model, please refer to the Programmer’s Guide available with OpenFOAM [22]. In this study, we used the solver simpleFoam, one of the standard solvers in OpenFOAM.

The simpleFoam solver, which is a steady-state solver for incompressible and turbulent flow, is used to solve the Reynolds-averaged Navier-Stokes (RANS) equations with the standard turbulence model by using the SIMPLE scheme [23]. The SIMPLE scheme is an iterative method to solve algebraic system of equations. A termination criterion of 10−4 is used for all field variables in this study. A brief introduction of iterative procedure of the SIMPLE method can be found in Miao et al. [20].

Step 1. Set up the initial guess field of pressure and velocity.

Step 2. Solve the momentum equation to compute the intermediate velocity field.

Step 3. Solve the pressure equation.

Step 4. Correct the velocities on the basis of the new pressure field.

Step 5. Update the guess fields of pressure and velocity.

Step 6. Repeat Steps 2 to 5 until convergence.

3. Numerical Experiments Setup

3.1. WRF Simulation Setup

The WRF model was configured with 3 two-way interactive, nested grids with horizontal grid spacing of 37.5, 7.5, and 1.5 km, and horizontal grid dimensions were 40 × 40, 56 × 46, and 51 × 51, respectively (Figure 1(a)). The innermost domain was set with the center located at 41.77°N, 123.41°E and covered most of Shenyang and its adjacent areas (Figure 1(b)). In the vertical direction we used 28 atmospheric eta levels, ranging from surface up to 50 hPa. There were 14 levels within the lowest 2 km in the atmosphere to better resolve the atmospheric boundary layer structure. The initial conditions were obtained from 1° × 1° Final Operational Global Analysis (FNL) data produced by the NCEP (National Centers for Environmental Prediction) every 6 hours, which was adopted as the background meteorological data fields. Since the typical features (i.e., boundary layer height and evolution process) of the atmospheric boundary layer structures of different seasons were different, two 42-h WRF numerical experiments (Summer-Run and Winter-Run) were designed. The details of the two WRF simulations are listed in Table 1. There were no obvious synoptic systems over northern China during the two WRF-simulated periods.

The land use dataset based on Moderate Resolution Imaging Spectroradiometer (MODIS) was used in this study, in which 20 land use categories were included, and the land use category of the innermost domain was shown in Figure 1(b). We used the Mellor-Yamada-Janjic PBL scheme [24, 25], which was capable of predicting TKE. Other physical options chosen were the Rapid Radiative Transfer Model (RRTM) long-wave scheme [26], the Dudhia shortwave scheme [27], the WRF Single-Moment (WSM) 3-class simple ice scheme [28], Kain-Fritsch (new eta) scheme [29], and the Noah land-surface scheme [30] with a single-layer UCM (Urban Canopy Model) [31, 32]. The first 16 hours of the simulations were used as spin-up, and WRF output interval was set to 1 hour. The coupling method of Miao et al. [18, 20] is used to couple WRF and OpenFOAM. The validations of OpenFOAM can be found at the study of Miao et al. [18, 20].

3.2. CFD Simulation Setup

The CFD model simulation domain covers a part of Northeastern University, China, located at the center of Shenyang. A Google Earth image and the CFD model domain with horizontal building configuration at the surface were shown in Figure 2. The actual buildings and their different shapes were simplified to cuboids, and the average height of the buildings was 24.5 m, with the tallest and lowest building being 81 m and 6 m. The horizontal grid interval was 10 m, and the horizontal grid dimension was 86 × 81. In the vertical, a nonuniform grid system with 50 levels was employed, in which the vertical grid interval was 5 m up to the 30th level and the grid interval of the topmost 20 levels was 10 m. The WRF data from one grid that was nearest to the location of the CFD domain was used to provide initial and boundary conditions to OpenFOAM. The velocity components and TKE at the inflow condition can be directly given by WRF simulation, while the momentum diffusion coefficient was indirectly calculated by Grisogono scheme [33, 34] and the TKE dissipation rate is parameterized by using these two parameters:where is the TKE, is its dissipation rate, is the turbulent exchange coefficient, and is an empirical constant. To simulate the wind and turbulent fields of different situation, the initial and boundary conditions of velocity components and TKE were set by the simulation of WRF of different moments (i.e., 1400 LST and 2200 LST). The zero-gradient boundary condition was employed at the outflow boundaries, including the top boundary. The building surface was set as being in a no-slip condition and the wall turbulence functions were employed to the grid points adjacent to the walls [35].

Three point sources located at different places were hypothetically set on the ground level in the built-up area, which were continuously emitted during the simulation period. The locations of the two point sources are shown in Figure 2(b) labeled by the crosses: one was set at the southern side of building-A, whose height was 81 m, the other one was set in the center of build-up area. The emission rate of each source was 10 ppm s−1, and the time step was set to 1 s. A series of sensitivity studies have been performed, resulting in that the concentration pattern was stable when the number of integral steps reached 5000, and thus the total number of the integrate steps was set to 7000 for all the dispersion numerical experiments.

4. Results and Discussions

In this section, the two series numerical experiments of the WRF-CFD coupled model on winter and summer would be presented.

4.1. Simulation Results in Winter
4.1.1. Simulation of WRF Winter Experiment

Figure 3(a) showed the observed sea level pressure field at 0800 LST 16 February 2012. At this time, the observed sea level pressure field exhibited a high pressure system over the Mongolia region and a low pressure system over the sea of Okhotsk. Shenyang was located in the uniform pressure field in front of Mongolia high circulation, without obvious synoptic system. Severe air pollution often happens in this situation over Shenyang area [36].

The measurements from 17 meteorological observatories (Figure 1(b)) located in Shenyang area were used to examine the accuracy of the WRF model. Figure 4 showed the diurnal variations of average observed and WRF-simulated 2-m temperature and 10-m speed. Both the observation and simulation are the mean value of the 17 stations or corresponding grids. It is noted that the observed temperature was higher than simulated one during the winter which was the heating season of Shenyang area (Figure 7(a)). This difference may owe to the underestimation of the anthropogenic heat in WRF. As noted in earlier studies, the anthropogenic heat has larger impact on surface air temperature in winter than in summer [37]. Despite this deficiency, the diurnal temperature cycle was replicated well by WRF model. The accuracy of wind speed simulation was lower than temperature; the simulated wind speed was higher than observed one. This is a general tendency by the model to overestimate wind speed found in many earlier studies [3840]. A new method to parameterize the effects that the unresolved topography exerts over the surface circulations has been implemented in WRFv3.4.1+, to correct the high wind speed bias [40]. In addition, the large anthropogenic heat in winter may also contribute to the difference between simulated and observed wind speed.

Figure 5 exhibited the 10-m wind vector fields and 2-m temperature fields of the innermost WRF domain on 16th of February 2012. During the studied period, the northwest wind dominated the Shenyang area, and the wind became weaker when the ground temperature gets higher. Furthermore, we found that the wind speed in the urban area of Shenyang was lower than the suburb areas, which was particularly obvious at 1000 LST and 1400 LST. At 2200 LST, the wind in the southeast off the city turned to the west.

Temperature stratification and boundary layer height are two important meteorological factors related to pollutant dispersion process. Temperature stratification affects vertical dispersion of air pollutant in atmospheric boundary layer, and low-level inversion which trapped humidity and pollutant is a key factor affecting regional air quality. The boundary layer height determines the volume of air which can be used to dilute pollutants. The WRF-simulated potential temperature profile and height-time section of potential temperature at the grid cell that was nearest to the location of the CFD domain in WRF were given in Figure 6. The MYJ scheme determines the PBL height using the TKE profile. Since the TKE is largest within the PBL, MYJ defines the top of the PBL to be the height where the TKE decreases to a prescribed low value [41], which leads to unreliable results during winter night. Instead, the diurnal variation of boundary layer height calculated using the 1.5-theta-increase method was shown in Figure 6. The 1.5-theta-increase method defines PBL heights as the level at which the potential temperature first exceeds the minimum potential temperature within the boundary layer by 1.5 K [42, 43]. The weak inversion can be seen at 0400 LST, and the surface layer began to turn into the mixed layer at 0800 LST; at 2000 LST, the inversion was formed again because of the cooling of surface (Figure 6(a)). As was shown in Figure 6(e), we can clearly see that the inversion was formed at around 1700 LST and lasted the whole night. Furthermore, the boundary layer height (mixed layer height) of the periods in which the inversion existed was lower than any other time of the day (Figure 6(c)), which was not good for the vertical diffusion of pollutant released from the ground.

4.1.2. Simulation of CFD Experiments in Winter

According to the analysis of WRF simulation results, two CFD numerical experiments with different atmospheric boundary conditions (unstable and stable) at different moments were set up, that is, 1400 LST and 2200 LST. At 1400 LST, an unstable mixed layer was fully developed and reached the highest boundary layer height (about 1.8 km) of the winter simulation period, and this unstable boundary structure was good for the vertical mixing process of the pollutant. On the contrary, a typical stable nocturnal boundary layer was established at 2200 LST, with a lower boundary layer height (about 300 m).

Figure 7 exhibited the simulated concentration (ppm) and wind vector fields for two heights at 1400 LST on 16th of February. At this moment, the ambient wind direction given by the WRF model was northwest, and the speed was about 3.5 m s−1 both at = 2.5 m and = 22.5 m. At = 2.5 m, wind vectors within the buildings were quite different from the free wind vectors on the boundary regions under the influences of the buildings. It was particularly obvious in the southern part of the build-up regions, where the distribution of the buildings was denser than the other areas. When the airflow passed through the building-A, two lee side vortexes were established; a strong one was on the southern lee side of the building-A and the other weak one was on the east. And wind speed around the edge of building-A was fastest in the build-up area. Since the building distribution was sparser at the height of 22.5 m, the wind field there was stronger and simpler. The two vortexes still existed on the lee side of building-A. Besides that, at the height of 22.5 m, most of wind vectors within the building followed the ambient wind direction.

At = 2.5 m, when the point source was located at P1, driven by the vortex on the lee side of building-A, a pollutant plume with concentrations higher than 5 ppm extended as an ellipse-pattern, and most of the pollutants were dispersed toward the southeast (Figure 7(a)). The point source P2 was located at a more open area; thus the pollutant concentration was high only at the southeast of the source due to the block effect of the building there (Figure 7(b)). Similar to the wind vector fields, the concentration field at = 22.5 m was simpler than that of lower level; most of the pollutants were transported by the ambient wind direction.

At 2200 LST, a typical nocturnal stable layer was formed near the surface (Figure 6(e)), which may inhibit the dispersion process of pollutants. Although the ambient wind speed of 2200 LST was higher than that of 1400 LST, the pollution could be more serious at 2200 LST with the same point source. For example, in the P1 point source experiment, the polluted area with concentration higher than 0.1 ppm was smaller at night (Figures 8(a) and 8(c)), which reflected that the diffusion ability was weaker at that moment, and more pollutants were trapped in the vortex on the leeward side of building-A.

For the P2 point source experiment, as the wind speed increased, the shape of pollution plume was changed a little, and at the height of 2.5 m more pollutants were transported along the channel formed by the buildings located on the north of P2. Unlike P1, placing on the leeward side building, the P2 point source was set up in a more open area; the pollution of P2 was dominated by the convection process rather than the turbulent diffusion process. Therefore, comparing to the daytime unstable situation, as wind speed increased, the convection of the pollutants strengthened, and then the pollution was relaxed a little at night (Figures 8(b) and 8(d)).

In short, it was found that buildings may alter the flow and concentration fields significantly. The position of the point source had significant influences on the concentration pattern. Comparing the flow and dispersion pattern at two moments, it was found that the convection and turbulent diffusion were the two important processes to determine the pollution serious level. Within the nocturnal stable atmospheric boundary structure, although the ambient wind speed was increased, the pollution within the buildings may be more serious because of the weaker turbulent diffusion.

4.2. Simulation Results in Summer
4.2.1. Simulation of WRF Summer Experiment

Figure 3(b) showed the observed sea level pressure field at 0800 LST 13 August 2012. Similar to the winter case, the Shenyang area was under the control of this high pressure system with weak synoptic system.

Figures 4(c) and 4(d) exhibited the diurnal variations of the observed and simulated 2-m temperature and 10-m wind speed, to examine the WRF output. The simulated temperature was a little lower than observed one at night, and the temperature during the daytime was simulated well. Despite the fact that the simulated wind speed was slightly higher, the diurnal variation of the wind speed was simulated well. On the whole, the simulation result in summer was better than that in winter.

The 10-m wind vector fields and 2-m temperature fields of the innermost WRF domain on 16th of February 2012 were shown in Figure 9. At 0200 LST, the main direction of the wind was northeast, while the wind turned to north at 0600 LST. The wind speed increased as the land surface was heated by the sun. It was noted that at 1800 LST, the wind direction within the built-up areas of Shenyang was different from that in the suburb area; the wind in the built-up areas was northwestern, while the wind in the suburb area was northern. At 2200 LST the wind turned to north in the simulation region. We can also obviously find the urban heat island effect in Shenyang at 0200 LST, 1800 LST, and 2200 LST; the 2-m temperature over the built-up area was significantly higher than that of the suburb area.

The potential temperature profile, height-time section of the potential temperature profile, and time series of boundary layer height simulated by WRF were given in Figure 6. The surface-based inversion formed at 0400 LST, and the intensity of inversion below the height of 500 m was 2.3 K per 100 m height, which was stronger than the inversion that happened in winter at the same time; and the depth of the inversion layer was also thicker than that of winter (Figures 6(a) and 6(b)). At 0800 LST, since the ground was warmed by the sunlight, a thin mixed layer was generated near ground, and the temperature inversion was detached from the ground at the height of 250 m. After 0800 LST, the low-level inversion layer disappeared gradually; at 1200 LST, a mixed layer was formed in the lower atmosphere from ground to the height of 700 m (Figure 6(b)). In the afternoon, the mixed layer extended well to 1.5 km height (Figure 6(d)), favoring the vertical dilution of pollutant. At night, the ground was cooled by the nocturnal radiation, causing the generation of the inversion stable layer near the ground. Comparing to the nocturnal stable layer of the winter case, the boundary layer height of the summer experiment was lower from 0000 LST to 0600 LST because of the calm ambient wind background (Figures 6(c) and 6(d)).

4.2.2. Simulation of CFD Experiments in Summer

Similar to the winter CFD numerical experiments, two experiments with different atmospheric boundary structures were designed at 1400 LST (unstable case) and 2200 LST (stable case).

Figure 10 illustrated the simulated concentration (ppm) and wind vector fields at different height at 1400 LST on 13th of August, when the ambient wind direction was north. It was noted that two counterclockwise vortices formed behind building-A, which was a typical double-eddy circulation reported by Hunter et al. [44] and Zhang et al. [45]. Comparing to the wind vector field at = 2.5 m, the wind vector field of 22.5 m was much simpler, and the north flow was dominant in most of the region.

At the lower level, the distributions of pollutant concentration were consistent with the flow pattern. When the point source was located at P1 or P2, the pollutants were transported and trapped in the circulations of the two vortexes behind building-A.

Comparing the concentration fields of point sources P1 and P2 at two heights, it was interesting to find that the concentration within the red rectangle in Figure 10 was higher at the height of 22.5 m than the lower level (2.5 m height). The phenomenon can be explained by the wind fields at different height. At the height of 2.5 m, since red rectangle region was surrounded by buildings, where the wind fields were detached from the ambient boundary wind and were dominant by the east wind; therefore the pollutant released from the northern P1 and P2 point sources cannot be transported to this region directly. In the contrary, at the height of 22.5 m, the red rectangle region was dominant by the north wind; therefore, the pollutant can be easily transported there.

At 2200 LST, when the south ambient wind blew, a quite different flow pattern was presented within the built-up area (Figure 11) compared to that of 1400 LST. A strong divergence zone can be found in the south of building-A at the height of 2.5 m, indicated by the red dashed rectangle in Figure 11. When the airflow blew from the south, a typical step-up notch was formed by building-A (81 m) and building-B (24 m); as the vertical wind field presented in Figure 12, the south inlet wind was blocked by the taller building (building-A) and descended along the windward wall of building-A and then diverged near the ground.

At the lower level, the wind field around the point sources (P1 and P2) was totally rearranged by the effects of the buildings around. At the height of 22.5 m, since the sparser building distribution, the wind field was more close to the ambient wind. Therefore, the wind fields of the height of 2.5 m and 22.5 m were quite different, as well as the pollutant distributions. At = 2.5 m, most of the pollutant released from P1 and P2 source was transported to south by the divergence zone formed on the south side building-A.

In these two cases, it was found that the change of the ambient wind direction can significantly affect the flow fields and pollutant distributions within the building cluster.

5. Discussion and Conclusions

This study examined urban airflow and dispersion pattern of pollutant released from hypothetical sources in Northeastern University of Shenyang, China, by using a WRF-CFD coupled model. Two WRF numerical experiments (Summer-Run and Winter-Run) were performed to study the typical features of the atmospheric boundary layer structures of different season. Then the coupled WRF-CFD was used to study how the atmospheric boundary structure, ambient wind, and source position affect the wind field and pollutant distribution within the built-up areas. It was found that the pollutant dispersion pattern and the pollution degree dispersion were complicated in the presence of real building clusters and under varying atmospheric boundary structures, ambient wind, and locations of pollution sources.

From this preliminary numerical study, it was found that the atmospheric boundary structure played a crucial role on the pollution within the building cluster, which determined the potential turbulent diffusion ability of the atmospheric surface layer. For example, within a nocturnal stable atmospheric boundary condition, although the ambient wind speed was increased, the pollution within the buildings may be more serious because of the weaker turbulent diffusion.

The change of the ambient wind direction can significantly affect the dispersion pattern of pollutants, which was a more sensitive factor than the ambient wind speed. Therefore, when a building cluster was designed, it was important to consider the frequency of the ambient wind direction there.

Under a given atmospheric state, the location of the pollution source would dramatically determine the pollution patterns within the built-up areas. When the pollutant source was set up on the leeward side of the building, the pollutant dispersion pattern was dominated by turbulent diffusion process which was affected by atmospheric boundary structures, while the convection process dominated when the pollutant source was set up in a more open area. Besides that, the distribution of pollutant concentration was driven by the flow pattern.

Finally, it should be noted that the wind fields and dispersion patterns within the building clusters were complicated, which cannot be obtained from the ambient wind. The WRF-CFD numerical evaluation can be used as a reliable method to understand the complicated flow and dispersion within the built-up areas.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the China Meteorological Administration Special Public Welfare Research Fund under Grant no. GYHY201106033 and the National Natural Science Foundation of China under Grant no. 41175004.