Impacts of Urban Layouts and Open Space on Urban Ventilation Evaluated by Concentration Decay Method

Qun Wang 1, Mats Sandberg 2, Yuanyuan Lin 1, Shi Yin 3 ID and Jian Hang 1,* 1 School of Atmospheric Sciences, Sun Yat-Sen University, Guangzhou 510275, China; wangqun_sysu@163.com (Q.W.); linyy8@mail2.sysu.edu.cn (Y.L.) 2 Laboratory of Ventilation and Air Quality, University of Gävle, SE-80176 Gävle, Sweden; Mats.Sandberg@hig.se 3 Department of Mechanical Engineering, The University of Hong Kong, PokFuLam Road, Hong Kong, China; shiyin103@gmail.com * Correspondence: hangj3@mail.sysu.edu.cn; Tel.: +86-20-8411-0375


Introduction
The increase in number of vehicles in cities and the ongoing urbanization worldwide are causing more concerns about urban air pollution.The urban canopy layer (UCL) is defined as the outdoor air volume below the rooftops of buildings [1].Raising UCL ventilation capacity by the surrounding atmosphere with relatively clean air has been regarded as one effective approach to diluting environmental pollutants [2][3][4][5][6][7][8] and improving the urban thermal environment in the hot summer [9,10] as well as reducing human exposure to outdoor pollutants [11][12][13].
In recent years, various ventilation indices have been applied for UCL ventilation assessment, such as pollutant retention time and purging flow rate [2,7,23], age of air and ventilation efficiency [3][4][5][6][7]23], in-canopy velocity and exchange velocity [22,25], net escape velocity [31], etc.Similar to indoor ventilation [32][33][34][35][36], the key point of applying these ventilation indices is based on the UCL ventilation processes as below: The surrounding external air is relatively clean and can be transported into urban areas to aid pollutant dilution with physical processes of horizontal dilution, vertical transport, recirculation of contaminants, and turbulent diffusion.In this sense, the tracer gas technique originated from indoor ventilation sciences has been used to predict outdoor ventilation by analyzing the final steady-state pollutant concentrations or the transient concentration decay history, which illustrates the processes of pollutant dilution and ventilation.For example, the concept of purging flow rate (PFR) was adopted to assess the net flow rate of flushing the urban domain [2] or the entire urban canopy layer [23] induced by the convection and turbulent diffusion.Moreover, the urban age of air [3][4][5][6][7] represents how long the external air can reach a place after it enters the UCL space.Hang et al. [3] first adopted the homogeneous emission method into CFD simulations [32] to calculate the urban age of air for quantifying the characteristics of wind supplying external air into UCL space for pollutant dilution.
In particular, the air change rate per hour is one of the most widely used indoor ventilation concepts [32][33][34][35][36], calculated by ACH = 3600Q T /Vol, where Q T is the total volumetric flow rate and Vol is the room volume.Later it is adopted to quantify the volumetric air exchange rate of the entire UCL volume per hour, including two-dimensional street canyon ventilation assessment [37,38] and three-dimensional UCL ventilation modeling [39][40][41][42][43][44].Previous researchers usually calculated outdoor ACH indexes based on the volumetric flow rate obtained by integrating the mean value of the component of the mean velocity normal to the UCL boundaries and the effective flow rate due to turbulence based on integrating half the standard deviation (rms-value) of the velocity fluctuations on street roofs [37][38][39][40][41][42][43][44].Specially, vertical turbulent exchange across a street roof has been verified to significantly influence UCL ventilation and pollutant removal in urban-like geometries [37][38][39][40] because street roofs are open at the top with a large area.However, the mean velocity often exhibits recirculation and the velocity fluctuations across the UCL open roof are bidirectional and therefore air or pollutant can return to the given UCL space several times.Thus, these two ACH indexes do not represent the entire UCL air volume that is really exchanged ACH times in an hour by external air, i.e., they are not the actual air change rate per hour in UCL models calculated by the actual flow rate flushing or leaving this space and never returning again.There are efficiency problems in UCL ventilation assessment [3,5].
The purpose of this paper is to predict the actual or net air change rate per hour in UCL models and how it is related to urban morphologies and atmospheric conditions.The net or actual air change rate is defined as the exchange of air flushing the space and never returning to the UCL space again.It depends on several time scales, i.e., the overall turn over time (air volume (Vol) divided by the total volumetric flow rate (Q T )), purging time for air moving from entrances to exits, time constant for local vortices and turbulent mixing.The concentration decay method has been widely used to predict the actual ACH and age of air in indoor environments in which the concentration temporal profile accords with the exponential decay law and the concentration decay history shows the pollutant dilution and ventilation processes in room ventilation [32][33][34][35][36]45,46].However, to date, this method has not been introduced to quantify the actual or net ACH index and age of air in UCL models.
Therefore, as a novelty, this paper aims to verify whether the temporal decay profile of concentration in UCL models accords with the exponential decay law and explore the effectiveness of applying concentration decay method for outdoor ventilation assessment.As a start, the typical medium-dense urban canopy layers (H/W = 1, λf = λ p = 0.25) are first studied.The effects of urban size, ambient wind directions, overall urban forms, and open space arrangements are evaluated under neutral atmospheric conditions.

Turbulence Models for Urban Airflow Modeling
As outlined by the review papers, CFD simulations have been widely applied in urban flow/dispersion modeling outdoor in the last two decades [47][48][49].Large eddy simulations (LES) are known to perform better in predicting turbulence than the Reynolds-Averaged Navier-Stokes (RANS) approaches.However, there are still some challenges involved with widely applying LES because of the strongly increased computational requirements, the development of advanced sub-grid scale models, and the difficulty in specifying appropriate time-dependent inlet and wall boundary conditions [16,17,20,26,38,39].Therefore, in quantitative work one has to adopt time-averaged turbulence models (i.e., RANS approaches).Actually, in spite of their deficiencies in predicting turbulence, steady RANS models have become the most popular CFD approaches for urban flow modeling, and the standard k-ε model is one of the most widely adopted in the literature [2,3,[5][6][7][9][10][11][21][22][23], with successful validation by wind tunnel measurements and good performance in predicting mean flows.This paper selects the standard k-ε model and the CFD validation case is conducted to evaluate the reliability of CFD methodologies by wind tunnel data.

Model Description in the CFD Validation Case
As shown in Figure 1a, Brown et al. [50] measured turbulent airflows in a seven-row and 11-column cubic building array with a parallel approaching wind.Building width (B), building height (H), and street width (W) are the same (B = H = W = 15 cm, H/W = 1, building packing densities λ p = λ f = 0.25, Lx = 13H, Ly = 21H).The scale ratio to full-scale models is 1:200.x, y, and z are the stream-wise, span-wise (lateral), and vertical directions.x/H = 0 represents the location of windward street opening.y/H = 0 is the vertical symmetric plane of the middle column.Point Vi represents the center point of the secondary street No. i.The measured vertical profiles of time-averaged (or mean) velocity components (stream-wise velocity u(z), vertical velocity w(z)), and turbulence kinetic energy κ(z) at Points Vi are used to evaluate CFD simulations.
In the CFD validation case, the full-scale seven-row building array is numerically investigated (B = H = W = 30 m, λ p = λ f = 0.25, H/W = 1; see Figure 1b).According to the literature [23,[51][52][53], the lateral urban boundaries hardly affect airflows in the middle column as the wind tunnel model is sufficiently wide in the lateral (y) direction (Ly = 21H).Thus, it is effective to only consider half of this middle column (shaded area in Figure 1a) to reduce the computational requirement.Figure 1b shows model geometry, grid arrangements, CFD domain, and boundary conditions.In particular, the distances of UCL boundaries to the domain inlet, domain outlet, and domain top are 6.7H, 40.3H, and 9.0H, respectively.Zero normal gradient condition is used at the domain top, domain outlet, and two lateral domain symmetry boundaries.At the domain inlet, the measured power-law profile of time-averaged (or mean) velocity U 0 (z) in the upstream free flow is adopted (see Equation (1a)) (Brown et al. [51]).Moreover, the vertical profiles of turbulent kinetic energy k(z) and its dissipation rate (ε) at the domain inlet are calculated by Equations (1b) and (1c) [23,[51][52][53]: where U ref (=3.0 ms −1 ) is the reference velocity at the building height (H = 30 m), C µ is a constant (=0.09), u * is the friction velocity (=0.24 ms −1 ), and κ v is von Karman's constant (=0.41).In the CFD validation case, the full-scale seven-row building array is numerically investigated (B = H = W = 30 m, λp = λf = 0.25, H/W = 1; see Figure 1b).According to the literature [23,[51][52][53], the lateral urban boundaries hardly affect airflows in the middle column as the wind tunnel model is sufficiently wide in the lateral (y) direction (Ly = 21H).Thus, it is effective to only consider half of this middle column (shaded area in Figure 1a) to reduce the computational requirement.Figure 1b shows model geometry, grid arrangements, CFD domain, and boundary conditions.In particular, the distances of UCL boundaries to the domain inlet, domain outlet, and domain top are 6.7H, 40.3H, and 9.0H, respectively.Zero normal gradient condition is used at the domain top, domain outlet, and two lateral domain symmetry boundaries.At the domain inlet, the measured power-law For this CFD validation case, hexahedral cells of 531,657 and the minimum cell sizes of 0.5 m (H/60) are used (the medium grid).Following the CFD guidelines [54,55], four hexahedral cells exist below the pedestrian level (z = 0 m to 2 m).The grid expansion ratio from wall surfaces to the surrounding is 1.15.With such grid arrangements, the normalized distance from wall surfaces y + (y + = yu τ /ν) ranges from 60 to 1000 at most regions of wall surfaces within UCL space.
No slip boundary condition with standard wall function [56] is used at wall surfaces.The practice guidelines for setting the upstream and downstream ground are followed to reproduce a horizontally homogeneous atmospheric boundary layer surrounding urban areas [46,57]: The roughness height k S and the roughness constant C S are correlated by the aerodynamic roughness length z 0 : Note that the distance between the center point P of the ground adjacent cell and the ground surface is Y P = 0.25 m.Fluent 6.3 does not allow k S to be larger than Y P [57].Thus, according to van Hooff and Blocken [46], a user-defined function is used to set the roughness constant C S = 4 because Fluent 6.3 does not allow it to be greater than 1.Then as C S = 4 and z 0 = 0.1 m, k S = 0.245 m is smaller than Y P = 0.25 m.
ANSYS FLUENT with the finite volume method is used to predict the steady isothermal urban airflows [56].The SIMPLE scheme is adopted to couple the pressure to the velocity field.The first-order upwind scheme is not appropriate for all transported quantities since the spatial gradients of the quantities tend to become diffusive due to a large viscosity.Thus, as recommended by the literature [54,55], all transport equations are discretized by the second-order upwind scheme for better numerical accuracy.The under-relaxation factors for pressure term, momentum term, k, and ε terms are 0.3, 0.7, 0.5, and 0.5, respectively.The solutions do not stop until all residuals become constant.The fine grid with a minimum size of 0.2 m (H/150) is also adopted to perform a grid independence study.The above CFD settings on domain sizes, grid arrangements and boundary conditions fulfill the major requirements given by CFD guidelines [54,55].In particular, steady CFD simulations are first carried out about 4000 iterations with the first-order upwind scheme, then continued with the second-order upwind scheme until all residuals became constant.The residuals reached the following minimum values or less: 10 −4 for the continuity equation, 0.5 × 10 −5 for the velocity components and k, 0.5 × 10 −5 and 0.5 × 10 −4 for pollutant concentration and ε.After solving the steady-state airflow field, the uniform initial concentration C(0) is defined in the entire urban canopy layer, then the unsteady concentration decay history is calculated and recorded for ventilation assessment as the solved flow field remains constant.

Model Description in All CFD Test Cases
As displayed in Table 1 and Figure 2a-d, two groups of medium-dense UCL models (H = B = W = 30 m, λ p = λ f = 0.25) are investigated in CFD simulations.The row and column numbers are referred to the numbers of buildings along the stream-wise direction (x, main streets) and span-wise direction (y, secondary streets).Wind directions are represented by the angles (θ = 0 • to 90 • ) between the approaching wind and the main streets.In Group I (Table 1 and Figure At the domain outlets and domain top, the zero normal gradient boundaries are used. In Group II (Table 1 and Figure 2d), for UCL models with five rows and five columns, open space arrangements are included.In total 24 test cases of Case Here Oij represents only one building of position i-j (2-1, 2-2, 2-3, 2-4, 3-3, or 3-4) (see Figure 2d) is removed to attain an open space effect.Computational domain and boundary conditions are similar to Group I.
For all the test cases in Table 1, the grid arrangements are similar to the medium grid of the CFD validation case.The total number of hexahedral cells is from 531,657 to 3,360,096.The above CFD settings on domain sizes, grid arrangements, and boundary conditions fulfill the major requirements recommended by the best CFD guidelines.[23,[37][38][39][40][41][42] To analyze the flow field, the velocity components and wind speed are all normalized by the reference velocity in the upstream free flow at the same height, i.e., Equation (1a) at the domain inlet U 0 (z) = U re f × (z/H) 0.16 .For example, for velocity at z = 15 m = 0.5H, the velocity is normalized by dividing it by U 0 (z = 15 m) = 2.685 ms −1 .

Volumetric Flow Rates and the Corresponding ACH
To quantify the ventilation capacity, the mean volumetric flow rates are calculated by integrating the normal air velocity across all UCL boundaries (Equation (4a)); moreover, the effective flow rate through street roofs due to turbulent exchange is defined by integrating the fluctuation velocity across open street roofs (Equation (4b)) [23,[37][38][39][40][41][42]: where, in Equation (4a), → V is the velocity vector, consisting of three time-averaged (or mean) velocity components (

and
→ n and A are the normal direction and area, respectively, of UCL boundaries.In Equation (4b), A roo f is the area of street roofs, is the fluctuation velocity based on the approximation of isotropic turbulence in all k-ε turbulent models where u , v , w are the stream-wise, span-wise, vertical velocity fluctuations It is worth noting that the values of +0.5 and −0.5 in Equation (4b) represent the fact that turbulent fluctuations induce the same upward and downward air fluxes across street roofs [23,[37][38][39][40][41][42].As pollution sources exist within UCL space, the literature confirmed that turbulent fluctuations across street roofs significantly contribute to pollutant removal since the upward pollutant flux through street roofs usually exceeds the downward pollutant flux because the pollutant concentration below street roofs is higher than that above them [21,31].Due to the much greater total area of street roofs, the effective flow rates by turbulent exchange across street roofs Q roo f (turb) are also important for the urban canopy layer ventilation compared to mean flows across UCL boundaries (Equation (4a)).Thus, this paper mainly adopts Q roo f (turb) to assess the effect of turbulent fluctuations on UCL ventilation [23,[37][38][39][40][41][42][43][44].

Actual or Net ACH and Urban Age of Air by the Concentration Decay Method
In urban ventilation sciences, the flow rates across a control surface are usually defined and recorded by integrating the normal velocity through it.Generally the complicated circulating flow is generated in an urban area, hence there are fluid particles coming back into UCL space across UCL boundaries.Moreover turbulent fluctuations induce upward and downward fluxes across the open street roofs, representing a large amount of fluid particles leave and re-enter UCL space.Thus, similar to the indoor environment [35,36], fluid particles flowing across UCL boundaries can be divided into two groups, one is for air leaving UCL space for never returning again and the other for air returning or revisiting UCL space (it has been in the UCL space before).The first group is much more significant to UCL ventilation and pollutant dilution.Therefore, only a fraction of the flow rates defined in Equations (4a) and (4b) contributes to flushing UCL space by external air or diluting pollutants within it.Furthermore, the ACH indexes defined in Equation ( 6) do not represent the entire UCL air volume that is really exchanged ACH times in an hour by external air.
To attain and assess the actual and net ACH induced by mean flows (ACH T ) and turbulent exchange (ACH turb ) across UCL boundaries, the concentration decay method is introduced into CFD simulations of UCL ventilation modeling.Similarly with indoor ACH by concentration decay method [32,45,46], if the predicted temporal decay profile of spatial mean concentration in the entire UCL volume accords with the exponential decay law, this decay rate is correlated to air change rate per hour in UCL models.
After the steady-state flow field is solved, a uniform initial time-averaged concentration of tracer gas (CO, carbon monoxide) is defined in the entire UCL air volume at time of t = 0 s (C(0) = 1.225 × 10 −7 kg/m 3 , see Figure 2e).Then the transient concentration decay process C(t) is numerically simulated and recorded while the flow field keeps steady.
The unsteady governing equation of time-averaged concentration C(t) is: where u j is the time-averaged (or mean) velocity components (u, v, w), and D m and D t are the molecular and turbulent diffusivity of pollutants or tracer gas.Here D t = ν t /Sc t , ν t is the kinematic eddy viscosity and Sc t is the turbulent Schmidt number (Sc t = 0.7) [3][4][5][6][7][21][22][23].
The decay duration of all test cases is sufficiently long (400 s).The default time step is dt = 1 s.In the CFD validation case, time steps of dt = 0.5 s and 2 s are also tested.The second-order upwind scheme is used for Equation (7) with the under-relaxation factors of 1.0.For each time step, the residual of Equation ( 7) reaches a value below 10 −10 .For Equation (7), the inflow pollutant concentration at the domain inlet is set as zero and zero normal flux condition is used at all wall surfaces; moreover, zero normal gradient condition is applied at the domain top, domain outlet, and domain lateral boundaries.
The net and actual ACH index is calculated by the decay rate of spatial mean concentration in the entire UCL space (<C(t)>/C(0)) (Equation ( 8)) [32][33][34]: where C(0) = 1.225 × 10 −7 kg/m 3 is a uniform initial concentration, Vol is the entire volume of UCL space, <C(t)> is the spatial mean concentration in the entire UCL space, and 3600 denotes one hour or 3600 s.Moreover, the concentration decay method can be used to predict urban age of air (τ p s) at a given point (see Equation ( 9)) [32,36]: Obviously, the physical meanings of Equations ( 8) and ( 9) are that the larger decay rate (K) of <C(t)>/C(0) represents greater ACH and the pollutant dilution rate of the entire UCL volume; moreover, a larger K at a point denotes that it requires a shorter time for external air to reach this point (i.e., the age of air is smaller).

Evaluation of CFD Simulations Using Wind Tunnel Data
To validate CFD simulations, Figure 3 compares wind tunnel data and CFD results in the CFD validation case by using the fine and medium grid, including time-averaged stream-wise velocity u(z) and turbulent kinetic energy k(z) at Point V2 and Point V5.Compared to wind tunnel data, the standard k-ε model with present grid arrangements can provide results of mean flows u(z) in good agreement with wind tunnel data, but does a little worse at predicting k(z) and can only predict the shape of vertical profile well.Such findings are similar with the same CFD validation studies the literature [23,[51][52][53].In addition, CFD results with the fine grid are only slightly different from those with the present medium grid.Given the results from the validation tests, we conclude that the application of the standard k-ε model with present grid arrangement is acceptable for the purposes of our research.There are relatively large discrepancies of stream-wise velocity between the simulated and measured results at the higher area.The possible reason is the ratio of vertical grid size to the stream-wise grid size is relatively large at much higher levels above the building roofs (i.e., at the higher area).However, the simulation results below and near the building roofs are performed in good quality, which is more important for UCL ventilation assessment.
literature [23,[51][52][53].In addition, CFD results with the fine grid are only slightly different from those with the present medium grid.Given the results from the validation tests, we conclude that the application of the standard k-ε model with present grid arrangement is acceptable for the purposes of our research.There are relatively large discrepancies of stream-wise velocity between the simulated and measured results at the higher area.The possible reason is the ratio of vertical grid size to the stream-wise grid size is relatively large at much higher levels above the building roofs (i.e., at the higher area).However, the simulation results below and near the building roofs are performed in good quality, which is more important for UCL ventilation assessment.
To further quantitatively evaluate the CFD models with the present grid arrangement and standard k-ε model, several statistical performance metrics are calculated, including the mean value, the standard deviation (named St dev.here), the fraction of predictions (here it is the CFD result) within a factor of two of observations (wind tunnel experiment data) named FAC2, the normalized mean error (NMSE), the fraction bias (FB), and the correlation coefficient (R).Results of ( ) at point V2 and V5 and ( ) at point V5 are shown in Table 2.According to the recommended reference criteria [52], for ( ) prediction, a good performing simulation model is supposed to meet the statistical metrics standards as below: FAC2 ≥ 0.5; NMSE ≤ 1.5; −0.3 ≤ FB ≤ 0.3.All the results satisfy the recommended criteria except the FB value of ( ) at point V5, which exceeds the upper limit, while its FAC2 just reaches the threshold.Furthermore, ( ) has a quite low R between wind tunnel data and CFD results.Results suggest a poorer quality of numerical prediction of ( ) than ( ).As for ( ), although it is overestimated as the FR is positive, the value of R is particularly high (0.93) at both V2 and V5, which implies a credible prediction of ( ).Assessing models' acceptance requires considering all relevant performance metrics rather than one specific index, so the CFD models are considered to have quite a satisfactory performance on the whole.Finally, with this CFD setup, Figure 3c shows the effects of time steps on the decay history of spatial mean concentration (<C(t)>/C(0)).The decay rates for three time steps (dt = 0.5, 1, 2 s) are almost the same, confirming that the present time step (dt = 1 s) is good enough for ACH prediction.

Flow and Concentration Decay in an Example Case [5-5, 0°] (θ = 0°)
As an example to analyze the concentration decay history related to UCL ventilation, Figure 4 shows 3D streamline, normalized velocity (    V u v w ), normalized lateral velocity ( v ) in the plane of z = 0.05H, and normalized concentration (C(t)/C(0)) at a time of 10 s to 400 s in z = 0.5H in Case [5-5, 0°].Here the velocity and lateral velocity are normalized by the freestream velocity at the same height of z = 0.05H (see Equation ( 1)), similar to the literature [17,18].Then Figure 5 displays the concentration history (C(t)/C(0)) and the age of air (τp) at various points in z = 0.5H in Case [5-5, 0°].In the main streets (Figure 4a), the flow is channeled toward downstream regions.In the secondary streets 3D vortices and helical flows exist in building wake regions.Across the lateral UCL boundary (at y = To further quantitatively evaluate the CFD models with the present grid arrangement and standard k-ε model, several statistical performance metrics are calculated, including the mean value, the standard deviation (named St dev.here), the fraction of predictions (here it is the CFD result) within a factor of two of observations (wind tunnel experiment data) named FAC2, the normalized mean error (NMSE), the fraction bias (FB), and the correlation coefficient (R).Results of u(z) at point V2 and V5 and k(z) at point V5 are shown in Table 2.According to the recommended reference criteria [52], for u(z) prediction, a good performing simulation model is supposed to meet the statistical metrics standards as below: FAC2 ≥ 0.5; NMSE ≤ 1.5; −0.3 ≤ FB ≤ 0.3.All the results satisfy the recommended criteria except the FB value of k(z) at point V5, which exceeds the upper limit, while its FAC2 just reaches the threshold.Furthermore, k(z) has a quite low R between wind tunnel data and CFD results.Results suggest a poorer quality of numerical prediction of k(z) than u(z).As for u(z), although it is overestimated as the FR is positive, the value of R is particularly high (0.93) at both V2 and V5, which implies a credible prediction of u(z).Assessing models' acceptance requires considering all relevant performance metrics rather than one specific index, so the CFD models are considered to have quite a satisfactory performance on the whole.Finally, with this CFD setup, Figure 3c shows the effects of time steps on the decay history of spatial mean concentration (<C(t)>/C(0)).The decay rates for three time steps (dt = 0.5, 1, 2 s) are almost the same, confirming that the present time step (dt = 1 s) is good enough for ACH prediction.

Flow and Concentration Decay in an Example
As an example to analyze the concentration decay history related to UCL ventilation, Figure 4 shows 3D streamline, normalized velocity (V = u 2 + v 2 + w 2 ), normalized lateral velocity (v) in the plane of z = 0.05H, and normalized concentration (C(t)/C(0)) at a time of 10 s to 400 s in z = 0.5H in Case [5-5, 0 • ].Here the velocity and lateral velocity are normalized by the freestream velocity at the same height of z = 0.05H (see Equation ( 1)), similar to the literature [17,18].Then Figure 5 displays the concentration history (C(t)/C(0)) and the age of air (τ p ) at various points in z = 0.5H in Case [5-5, 0 • ].In the main streets (Figure 4a), the flow is channeled toward downstream regions.In the secondary streets 3D vortices and helical flows exist in building wake regions.Across the lateral UCL boundary (at y = 135 m), there are lateral helical airflows leaving or re-entering UCL volume.Points with bigger concentration decay rates (K) experience better ventilation, a greater dilution rate, and smaller age of air (τ p ). Figures 4b and 5a,b show that the ventilation in upstream regions is better than in downstream regions (Figure 5a), and that in the main streets it is better than in the secondary streets (Figure 5a); those (Point S1b-S4b) near the lateral UCL boundary are better than those (Point S1o-S4o) in urban center regions (Figure 5b). Figure 5b displays that K at Point S4b is near to Point S3b, and that at Point S4b near the UCL lateral boundaries it is much bigger than at Point S4o because there is a strong helical inflow across lateral boundaries, bringing in external air to help with pollutant dilution at Point S4b (Figure 4a).Similarly, in Figure 5c, the τ p at Points M1b to M4b in the main streets is much smaller than Points S1a to S4a and Points S1o to S4o in the secondary streets.Moreover, air at Points S1b to S4b near lateral UCL boundaries is younger than Points S1o to S4o far from lateral UCL boundaries.More importantly, the τ p at Points S4b, S4a, and S4o (89.5 to 249.4 s) verifies that the lateral inflow across UCL lateral boundaries significantly improves the ventilation at Point S4b.
There are qualitative differences between the decay curves in Figure 5a,b.Analyzing the behavior of the decay curves in Figure 5a, for example, after some time the decay rates at Points S1a and M1b are similar to each other.This is a manifestation that there is a coupling interaction between the two locations (more explanation can be found by comparing Figure 11.2 in [32]).However, in Figure 5b all curves exhibit different decay rates, which shows that is not feedback (no backflow against the wind direction) that connects the different locations (see also Chapter 11 in [32]).
inflow across lateral boundaries, bringing in external air to help with pollutant dilution at Point S4b (Figure 4a).Similarly, in Figure 5c, the τp at Points M1b to M4b in the main streets is much smaller than Points S1a to S4a and Points S1o to S4o in the secondary streets.Moreover, air at Points S1b to S4b near lateral UCL boundaries is younger than Points S1o to S4o far from lateral UCL boundaries.More importantly, the τp at Points S4b, S4a, and S4o (89.5 to 249.4 s) verifies that the lateral inflow across UCL lateral boundaries significantly improves the ventilation at Point S4b.There are qualitative differences between the decay curves in Figure 5a,b.Analyzing the behavior of the decay curves in Figure 5a, for example, after some time the decay rates at Points S1a and M1b are similar to each other.This is a manifestation that there is a coupling interaction between the two locations (more explanation can be found by comparing Figure 11.2 in [32]).However, in Figure 5b all curves exhibit different decay rates, which shows that is not feedback (no backflow against the wind direction) that connects the different locations (see also Chapter 11 in [32]).

Effect of Overall Urban Form and Ambient Wind Direction on UCL Ventilation
As an example, Figure 6 shows 3D streamline, normalized velocity, C(t)/C(0), and τp at various points in z = 0.5H in Case [5-5, 45°].Here the velocity is normalized by the freestream velocity at the same height of z = 0.5H.Obviously wind brings clean air into an urban area downstream for pollutant dilution.There are recirculation regions with small wind speed (Figure 6a).Thus, the concentration decay rate (the age of air) is relatively small (large) in downstream regions and recirculation regions (Figure 6b,c).

Effect of Overall Urban Form and Ambient Wind Direction on UCL Ventilation
As an example, Figure 6 shows 3D streamline, normalized velocity, C(t)/C(0), and τ p at various points in z = 0.5H in Case [5-5, 45 • ].Here the velocity is normalized by the freestream velocity at the same height of z = 0.5H.Obviously wind brings clean air into an urban area downstream for pollutant dilution.There are recirculation regions with small wind speed (Figure 6a).Thus, the concentration decay rate (the age of air) is relatively small (large) in downstream regions and recirculation regions (Figure 6b,c).
Then Figure 7a,b further displays the normalized velocity in Case [7-7, θ • ] (square overall urban form) and Case [10-5, θ • ] (rectangular overall urban form) with similar total UCL air volume.Obviously ambient wind directions significantly influence the flow pattern for both overall urban forms (Figure 7a,b).As θ = 0 • , there are lateral flows across lateral UCL boundaries.With oblique winds, the flows enter UCL across two sides in upstream regions and leave across the other two toward downstream.The velocity is relatively small in recirculation regions and the presence of street crossings produces considerable momentum and scalar exchange between neighbor streets.In particular, Figure 7a   There are qualitative differences between the decay curves in Figure 5a,b.Analyzing the behavior of the decay curves in Figure 5a, for example, after some time the decay rates at Points S1a and M1b are similar to each other.This is a manifestation that there is a coupling interaction between the two locations (more explanation can be found by comparing Figure 11.2 in [32]).However, in Figure 5b all curves exhibit different decay rates, which shows that is not feedback (no backflow against the wind direction) that connects the different locations (see also Chapter 11 in [32]).

Effect of Overall Urban Form and Ambient Wind Direction on UCL Ventilation
As an example, Figure 6 shows 3D streamline, normalized velocity, C(t)/C(0), and τp at various points in z = 0.5H in Case [5-5, 45°].Here the velocity is normalized by the freestream velocity at the same height of z = 0.5H.Obviously wind brings clean air into an urban area downstream for pollutant dilution.There are recirculation regions with small wind speed (Figure 6a).Thus, the concentration decay rate (the age of air) is relatively small (large) in downstream regions and recirculation regions (Figure 6b,c).Obviously ambient wind directions significantly influence the flow pattern for both overall urban forms (Figure 7a,b).As θ = 0°, there are lateral flows across lateral UCL boundaries.With oblique winds, the flows enter UCL across two sides in upstream regions and leave across the other two toward downstream.The velocity is relatively small in recirculation regions and the presence of street crossings produces considerable momentum and scalar exchange between neighbor streets.In particular, Figure 7a confirms that, for square overall urban form, UCL models with θ = 45° experience more recirculation regions and quicker wind reduction than those with θ = 0°.It is consistent with the literature that building arrays with θ = 45° produce greater flow resistances and worse UCL ventilation performance than θ = 0° [2,7,23,58,59].Such characteristics with θ = 45° produce adverse effects on its ventilation performance in contrast to θ = 0°.

Effect of Open Space Arrangements and Ambient Wind Direction
Open space arrangements have been regarded as one possible way to improve UCL ventilation.As shown in Table 1 (Group II), this subsection investigates 24 test cases with six kinds of open space arrangements under four wind directions (θ = 0 • , 15 • , 30 • , 45 • ).Note that Oij represents the building of position i-j is removed for better ventilation.Figure 8a displays the 3D streamline in two example cases; Figure 8b further summarizes the net ACH by the concentration decay method in all 24 test cases.

Effect of Open Space Arrangements and Ambient Wind Direction
Open space arrangements have been regarded as one possible way to improve UCL ventilation.As shown in Table 1 (Group II), this subsection investigates 24 test cases with six kinds of open space arrangements under four wind directions (θ = 0°, 15°, 30°, 45°).Note that Oij represents the building of position i-j is removed for better ventilation.Figure 8a displays the 3D streamline in two example cases; Figure 8b further summarizes the net ACH by the concentration decay method in all 24 test cases.] in Figure 8c).This percentage at some points can be from 12% to 65%, showing that τ p significantly decreases in the regions near these points.Obviously, the effects of open space on flow pattern satisfy the variations of τ p (Figure 8c,d) and overall ACH (Figure 8b).

Conclusions
As a novelty, this paper confirms the temporal decay profile of concentration in UCL models accords with the exponential decay law, and it is effective to introduce the concentration decay method into CFD simulations to predict the net air change rate (ACH) flushing UCL space and never returning.Street-scale (~100 m), medium-dense (λ f = λ p = 0.25, H/W = 1), urban-like geometries are studied under neutral atmospheric conditions.The standard k-ε model is first successfully evaluated by wind tunnel data.Then the flow pattern, the concentration decay rate and age of air (τ p ) at some points in UCL space, the net ACH are analyzed.
With a parallel approaching wind (θ = 0 • ), UCL models with larger urban sizes (Lx = Ly = 390 m) experience smaller ACH (~16.8-18.6 h −1 ) than smaller UCL models (Lx = Ly = 270 m, ACH~21.1-30.2h −1 ).The urban age of air τ p at points near upstream/lateral UCL boundaries is smaller (or the air is younger) than far from them.For the square overall urban form, the parallel wind attains greater net ACH than non-parallel wind (θ = 15 • , 30 • , 45 • ).For the rectangular overall urban form (Lx = 570 m, Ly = 270 m), ACH is greater than the square overall urban form (Lx = Ly = 390 m) under most wind directions (~21.1-28.7 h −1 as θ = 30 • to 90 • ), and the ventilation is the best as the approaching wind is parallel to the shorter urban size (~28.Similar to the purging flow rate [2,23], ACH calculated by the concentration decay approach has been proven effective to evaluate the effects of urban morphologies on the overall UCL ventilation capacity induced by mean flows and turbulent diffusions, which seems to be a better ventilation index than ACH calculated by the volumetric flow rates integrating the normal mean velocity or fluctuation velocity across UCL boundaries.Turbulent diffusions across open street roofs have been proven to significantly contribute to UCL ventilation, but further investigations are still required to analyze the net ventilation efficiency by mean flows and turbulence.
V j velocity vector Vol control volume x, y, z stream-wise, span-wise, vertical directions u, v, w stream-wise, lateral, vertical velocity components u , v , w stream-wise, lateral, vertical velocity fluctuations

Figure 1 .
Figure 1.(a) Model description of wind tunnel test; (b) setups in CFD validation case.

Figure 1 .
Figure 1.(a) Model description of wind tunnel test; (b) setups in CFD validation case.
2a), test cases without open space are named as Case [row number-column number, wind direction θ • ].For cases with a parallel approaching wind (θ = 0 • , Figure 2b), three cases are included (i.e., Case [5-5, 0 • ], Case [7-7, 0 • ], Case [10-5, 0 • ]) with only half of computational domain simulated.The distances from UCL boundaries to the domain top, domain outlet, domain inlet, and the domain lateral boundary are 9H, 40.3H, 6.7H, and 10H, respectively.Zero normal gradient boundary condition is used at the domain outlet, domain top, and the lateral domain boundary.At the domain inlet, Equation (1) is used to provide the boundary condition.For test cases with oblique wind (θ = 15 • , 30 • , 45 • , 60 • , 75 • , 90 • , Figure 2c), the full CFD domains are used.There are two domain inlets and two domain outlets.The distances from UCL boundaries to the domain top, domain outlets, domain inlets are 9H, 41H and 6.7H.At two domain inlets, the vertical profiles of time-averaged (or mean) velocity components in Equations (3a)-(3c) and turbulent quantities defined in Equations (1b) and (1c) are used to define boundary conditions:

Figure 3 .
Figure 3. CFD results and experimental data in the CFD validation case: (a,b) Vertical profiles of u(z) and k(z) at Point V2 and/or Point V5.(c) <C(t)>/C(0) with various time steps.
confirms that, for square overall urban form, UCL models with θ = 45 • experience more recirculation regions and quicker wind reduction than those with θ = 0 • .It is consistent with the literature that building arrays with θ = 45 • produce greater flow resistances and worse UCL ventilation performance than θ = 0 • [2,7,23,58,59].Such characteristics with θ = 45 • produce adverse effects on its ventilation performance in contrast to θ = 0 • .Then Figure 7c,d displays the net ACH by concentration decay method, ACH T and ACH turb in cases of Group I. Since Case [5-5, θ • ] and Case [7-7, θ • ] are of square overall urban form, the flows with θ = 0 • , 15 • , 30 • , 45 • are the same with those with θ = 90 • , 75 • , 60 • , 45 • .Obviously for Case [5-5, θ • ] (Lx = Ly = 270 m) and Case [7-7, θ • ] (Lx = Ly = 390 m), θ = 0 • attains smaller ACH T and ACH turb but bigger net ACH than θ = 15 • , 30 • , 45 • .Therefore, the parallel wind (θ = 0 • ) experiences the best overall UCL ventilation with the highest net ventilation efficiency.It can be explained that UCL models with non-parallel approaching wind (for example θ = 45 • in Figure 7a) experience more recirculation regions and weaker wind due to the greater blockage induced by buildings, which can reduce the flow rates flushing through UCL space that never return or decrease the actual or net air change rate of the entire UCL space.Then for Case [10-5, θ • ] with a rectangular overall urban form (Lx = 570 m, Ly = 270 m), as θ varies from 0 • to 90 • , ACH T rises from 16.1 h −1 to 29.0 h −1 , ACH turb first increases from θ = 0 • to θ = 45 • then decreases to θ = 90 • , and the net ACH rises from 14.3 h −1 to 28.7 h −1 .Thus θ = 90 • obtains the best overall UCL ventilation in Case [10-5, θ • ], in which the approaching wind is parallel to the shorter urban size (Ly = 270 m).Finally, the net ACH by the concentration decay method are always much smaller than the sum of ACH turb and ACH T .This verifies that UCL ventilation efficiency is limited.As discussed in Section 2.4.2, the recirculation flows in UCL space and turbulent fluctuations across street roofs induce a significant fraction of fluid particles to return or revisit UCL space across UCL boundaries after they leave it.However, UCL ventilation mainly depends on the flow rates flushing UCL space and leaving it to never return.Thus only a fraction of ACH turb and ACH T contributes to UCL ventilation and the actual air change rate is limited due to the ventilation efficiency problems.Atmosphere 2017, 8, 169 14 of 24 (c) Age of air at these points at z = 0.5H in Case [5-5, 0°]
7 h −1 as θ = 90 • ).With open space arrangements, the dilution capacity near/surrounding open space and in its downstream regions is enhanced; moreover, most open space arrangements are more effective at improving UCL ventilation under oblique wind directions (θ = 15 • , 30 • , 45 • ) than θ = 0 • , and the best ventilation improvement by open space appears as θ = 45 • .
at Points Vi are used to evaluate CFD simulations.

Table 2 .
Statistical performance metrics in CFD validation case.

Table 2 .
Statistical performance metrics in CFD validation case.
Atmosphere 2017, 8, 169 20 of 25 has smaller urban size and requires shorter time for the approaching wind to flow through and will be exchanged more times by the external air within one hour.Rectangular urban form (Case [10-5, θ°]) attains greater ACH (~21.1-28.7 h −1 ) than the square urban form (Case [7-7, θ°]~16.8-18.6 h −1 ) for most ambient wind directions except θ = 0° and 15° (~14.3 and 16.9 h −1 ).It can be explained by two example cases: for Case [10-5, 0°] there are 10 rows of buildings for the approaching wind to flush, much longer than Case [10-5, 90°], in which there are only five rows of buildings for wind to flow through.