Impact of Urbanization on the Predictions of Urban Meteorology and Air Pollutants over Four Major North American Cities

: The sensitivities of meteorological and chemical predictions to urban effects over four major North American cities are investigated using the high-resolution (2.5-km) Environment and Climate Change Canada’s air quality model with the Town Energy Balance (TEB) scheme. Comparisons between the model simulation results with and without the TEB effect show that urbanization has great impacts on surface heat ﬂuxes, vertical diffusivity, air temperature, humidity, atmospheric boundary layer height, land-lake circulation, air pollutants concentrations and Air Quality Health Index. The impacts have strong diurnal variabilities, and are very different in summer and winter. While the diurnal variations of the impacts share some similarities over each city, the magnitudes can be very different. The underlying mechanisms of the impacts are investigated. The TEB impacts on the predictions of meteorological and air pollutants over Toronto are evaluated against ground-based observations. The results show that the TEB scheme leads to a great improvement in biases and root-mean-square deviations in temperature and humidity predictions in downtown, uptown and suburban areas in the early morning and nighttime. The scheme also leads to a big improvement of predictions of NO x , PM 2.5 and ground-level ozone in the downtown, uptown and industrial areas in the early morning and nighttime.


Introduction
Rapid urbanization around the world leads to significant changes in urban settings [1,2]. The residential, industrial and commercial buildings and paved streets in urban areas significantly change the urban land use and land cover, and consequentially change the radiation balance, thermal distribution, and aerodynamics in urban areas [3][4][5]. In addition, the waste heat produced by heating or cooling buildings and by vehicles also modifies the urban surface energy balance. These changes in urban settings, along with the daily anthropogenic activities in urban areas, impose great impacts on urban meteorology and air quality [6].
Because urban meteorology is greatly affected by the surface energy balance [6,7], the heat fluxes from urban settings and daily anthropogenic activities, along with the momentum fluxes [8], can create distinct urban meteorological conditions such as meso-and micro-scale urban heat island (UHI) phenomenon [3,4], urban flooding, changes in precipitation patterns, and street canyon winds [6,[9][10][11][12]. In addition to the impact on urban meteorology, urbanization can also lead to elevated concentrations levels of gaseous pollutants and aerosols due to the release of chemical species into the atmosphere from industry activities, burning fossil fuels by vehicles [13,14] and heating buildings [15]. These directly released species can produce other pollutants such as ground-level ozone through photochemical reactions. The elevated concentration of the pollutants can perturb the energy balance near the surface by affecting the aerosol-cloud-precipitation micro-physics. In addition, the transport of species as well as the formation and transformation of aerosols in urban areas are influenced by the vertical turbulent mixing and the variational local circulation associated with the UHI effect.
Due to high urban population density and heavy dependence on service and infrastructure [1], people living in major urban centers are extremely vulnerable to hazardous weather events ranging from severe thunderstorms and snowstorms to cold and heat waves, and air pollution [16][17][18][19]. In addition to the risks and disruptions to people's daily life, these extreme weather events can also lead to huge losses in property and infrastructure [19]. In order to better prepare and respond to these events, accurate forecasts of the extreme weather events and heavy air pollution episodes in urban areas are needed. The forecast results can be used by different user groups for health, safety, comfort, security, planning and decision-making. Furthermore, the developed forecast systems, along with the archived long-term urban meteorology data sets, can be used for scenario-based urban planning projects to mitigate the UHI effect by minimizing energy consumption, and developing strategies for urban climate change mitigation and adaptation [20] The big diversity of size and shape of buildings and other compositions, as well as the highly heterogeneous spatial distribution of urban canopy elements impose a great challenge in creating a surface datum for the general descriptions of the radiative, thermal, moisture and aerodynamic impacts of the urban canopy elements on meteorological fields in all the urban areas [6]. In addition, these features also make measuring the impact of each individual urban component on urban environment at ground level difficult as other impacts can intervene in the measurement. It is also very difficult to estimate the heat fluxes and emissions of chemical species associated with anthropogenic activities due to high spatiotemporal distributions of the sources. Complicated physical, dynamic and thermal structures and turbulent motions with a wide range of spatial and temporal scales in urban boundary layer (UBL) add additional uncertainties to urban meteorology studies. Furthermore, because observing systems are not uniformly distributed in general, evaluation and validation can be difficult.
In spite of the challenges mentioned above, the growing needs for accurate predictions of urban meteorology and air quality for the comfort, safety and health of urban dwellers, along with the advent of powerful supercomputers and the advanced geographic information systems, have motivated numerous studies on the urbanization impact on the urban meteorological and chemical environment. In these studies, the urbanized high resolution numerical models have been widely employed to simulate the complicated physical, dynamical and chemical interactions between urban canopy and the atmosphere in the UBL [21][22][23][24][25][26][27][28][29].
To describe the urbanization impact on environment, different urban schemes, from simple and idealized schemes such as empirical models [30][31][32] and vegetation schemes [33] to more realistic and complicated urban canopy models such as the single-layer Town Energy Balance (TEB) scheme [21] and multi-layer urban scheme [23], have been proposed. Although parameterizations are used to describe the urban geometric properties, the urban canopy models are able to compute the contribution of each urban setting element to the urban surface energy balance based on the urban coverage types, urban canopy structure, urban fabric thermal and radiative properties, as well as the radiative interactions between each element based on physical laws [21,26]. Simulation results show that the urban canopy schemes improve the model predictions of meteorological fields and some chemical species in certain time periods of a day over urban areas [26][27][28][29]34].
Because previous studies focused on the improvement of the model prediction skills by the urban schemes in a particular urban area during a particular period, the mechanism underlying the urban impacts on meteorological and chemical fields have not been fully explored, and the common and different features of the impacts on different urban environment in different season have not been discussed. One purpose of this work is to examine the sensitivities of urban meteorological and chemical fields in different urban areas in both summer and winter to the TEB scheme, and to interpret the sensitivity results by applying analytical results about the relationships between the urban impacts and urban heat fluxes, vertical mixing, and the vertical profile of reference state [35]. Understanding the features of the urban impacts and the influence of physical and dynamical processes associated with urban impacts can not only shed physical insight into the urban impacts, but also be very important for the further improvement of the model prediction of urban meteorological and chemical fields. Another purpose is to evaluate the sensitivities to better understand the improvement in model predictions of urban meteorological and chemical fields by the TEB scheme by comparing the scores of predictions with and without the TEB scheme.
The paper is arranged as follows. The air quality model developed by the Environment and Climate Change Canada (ECCC) and the TEB scheme are introduced in Section 2. The surface heat fluxes in the urban area and the impact of the TEB scheme on predictions of urban meteorology and chemical species are examined and the associated mechanisms are discussed in Sections 3 and 4, respectively. In Section 5, the improvements in the model predictions due to the TEB scheme in the Great Toronto Area (GTA) are evaluated. Summary and discussions are given in the last section.

GEM-MACH
The numerical model used in this work is GEM-MACH (Global Environmental Multiscale Model-Modeling Air quality and Chemistry), a multi-scale chemical weather forecast model developed by ECCC. It is composed of dynamics, physics, and on-line chemistry and transport modules [36][37][38][39], and has been employed to address a variety of tropospheric air pollution problems ranging from ground-level ozone to particulate matter (PM) to acid rain by simulating the time evolution of size-resolved PM, photochemical oxidants, and acid deposition, as well as their interactions through gaseous, aqueous, and heterogeneous reactions and physical processes. Pollutant emissions include point-source and area-sources emissions, and the details of the emission inventories and process of transforming the inventory data into model-ready emissions can be found in [40,41] In this work, the GEM-MACH with the 2.5-km grid spacing (502 × 402 grid cells) is used. The model domain in our simulation covers Southern Ontario and part of north-east of U.S., and four big cities including Toronto, New York City, Chicago and Detroit are within the domain. The model has 81 levels with the model lid at about 65 km above the ground. The lateral boundary conditions and initial conditions of our simulations are from the 10-km regional GEM-MACH model simulations. The chemical tracers are initialized only at the beginning of simulations. The meteorological variables are initialized at the beginning of each 24-h simulation cycle at 00:00 UTC. The 10-km regional model has 750 × 620 grid cells and has the same vertical model level and resolution as the 2.5-km GEM-MACH. The initial conditions of meteorological fields in the 10-km regional model are provided by the data assimilation analyses produced by the ensemble variational (EnVAR) data assimilation system at foure synoptic hours. The EnVAR system developed by ECCC produces analyses by combining model forecasts, observation and time-dependent model error correlation statistics through the variational method. The surface variables including the lake surface temperature are initialized using land surface data assimilation analyses. In total, 21 surface types defined in ISBA are used in both TEB and non-TEB simulations. Numerical simulation results over the four big urban areas are used to reveal common and different features of the urbanization impacts on meteorological and chemical fields over different urban environment. Because these urban areas are close to lakes or ocean, the simulations results are also used to examine the urban impacts on lake breeze effect. More detailed description of the GEM-MACH configuration and settings can be found in [41].
To examine the impact of the TEB scheme in summer and winter, simulations with and without the TEB scheme are carried out for July 2015 and January 2016. In the TEB and non-TEB simulations, the same initial conditions, chemical emissions and (vertical and horizontal) boundary conditions are used. In the sensitivity studies, the averaged values of meteorological and chemical variables over 9 grid cells in downtown Toronto, 14 in downtown Chicago, 12 in downtown New York, and 13 in downtown Detroit are used To examine the impacts of the TEB scheme on the meteorological and chemical environment over urban areas.

The TEB Scheme
The TEB scheme is the single-layer canopy model in which interactions between urban canopy and the atmosphere occur on the lowest model level (thermodynamic level) which is above the roof level [21]. Wind, temperature and humidity within urban canyons are uniformly distributed in vertical but vary in horizontal. In the model, one roof, generic wall and road are used. The effect of street orientation is partially taken into account through averaging variables in all directions. Unlike empirical models and vegetation schemes in which the urban settings are absent, and the impacts of urban surface are estimated either by applying the observation-based statistics or by modifying physical parameters, the TEB model includes detailed descriptions of the 3-D shape of buildings, calculations of sensible and latent heat fluxes by taking into the consideration of heat conduction within each urban component including roofs, roads and walls as well as the coverage of water and snow on roofs and roads, and radiative interactions among the urban components based on physical laws. In the model, the evolution of the temperature of the layers of urban settings and water reservoir is described by a set of prognostic equations. Unlike the meteorological variables in GEM-MACH which are initialized every 24 h, the prognostic variables in the TEB scheme are initialized only at the beginning of the model integration.
In this work, the TEB scheme is chosen as the urban model in GEM-MACH. The details of generating the geometric, thermal, and radiative parameters for the TEB scheme are discussed in [29]. Among these parameters, building fraction and urban canyon aspect ratio are two important geometric parameters in computing the urban surface heat fluxes. Figure 1 shows that the spatial distributions of the building fractions in Toronto, Chicago, New York and Detroit are inhomogeneous with big values in downtown areas in Toronto and Detroit, and in the Brooklyn area in New York City.

Surface Heat Fluxes in Urban Areas
In an urban area, there are many sources of heat flux that contribute to the energy balance and the UHI phenomenon. The heat fluxes from vehicles, air conditioning systems, and the elements of urban setting including roofs, walls, and roads are discussed in this section.

Heat Released from Vehicles
Heat produced by both on-road and off-road vehicles are one of important contributors to the urban surface heat fluxes. Several approaches have been proposed to compute the heat fluxes from traffic [42][43][44]. In this work the traffic heat flux is computed using the approach proposed in [43]. This approach is based on the statistical relationship between heat fluxes and emissions of NO x (summation of NO 2 and NO) and CO from vehicles [43]. Because the hourly surface emissions of NO x and CO from traffic have been produced for GEM-MACH [40], the traffic heat fluxes can be produced easily as the spatial allocation of heat fluxes is not needed. More importantly, because spatiotemporal distributions of heat fluxes based on this approach are consistent with those of chemical emissions, the impact of traffic heat fluxes on the variation of chemical species can be better described in GEM-MACH. Furthermore, because the hourly emissions are produced for each day of a week in each month, the traffic heat fluxes produced by this approach can describe the differences in traffic heat flux between weekdays and weekend associated with the apparent different traffic patterns.
The diurnal variations of the weekday and weekend day mean traffic heat fluxes in July 2015 in the downtown area of Toronto, Chicago, New York, and Detroit are shown in Figure 2. It can be seen from the figure that although the diurnal variations of the heat fluxes in the four cities are similar with large values between 8:00 a.m. and 17:00 p.m. (local time) and small values at other times, differences can be identified. While there are two peaks in weekday heat flux in Toronto, one at around 8:00 a.m. and another at around 4:00 p.m., corresponding to the morning and afternoon rush hours, there is only one peak in three U.S. cities at around 6:00 p.m. The difference is due to the different diurnal profiles used for allocating emission inventory to hourly emissions, and it may be associated with the different traffic patterns in these cities. In addition, the magnitudes of the traffic heat fluxes during the daytime are quite different due to different traffic volumes. While the traffic heat flux maximum in New York City can be as much as 48 w/m 2 in July, it is, however, only 25 w/m 2 in Detroit. Figure 2 also shows that while the diurnal variation patterns of traffic heat fluxes in weekdays and weekend are quite similar, the magnitudes are quite different. The traffic heat fluxes in weekdays are stronger than those in weekend in the four cities, and the daytime traffic heat fluxes are much stronger in weekdays than in weekend in Toronto and Chicago. Furthermore, there is a time lag between the weekday and weekend maximum except Toronto. Weekend heat fluxes in other three urban centers reach to maximum few hours earlier than weekday heat fluxes. These features can be also seen in the diurnal variations of traffic heat flux in winter (not shown).

Heat Flux Produced by Air Conditioning Systems
In addition to traffic heat flux, the heat flux released from air conditioning systems in summer is calculated using a simple scheme based on the energy balance within buildings. In the scheme, the energy produced is equal to the heat penetrating into buildings through walls and roofs, and the scheme is turned on when the urban canyon temperature is above 30 • C. Because the period with temperature above 30 • C is very short, and high temperature occurs only in few days in July 2015, the heat flux is small in our simulation.

Heat Fluxes Produced by the Urban Settings
The three components of the urban setting contribute to the urban surface energy balance by absorbing/releasing sensible and latent heat from/to the atmosphere. In the TEB scheme, several steps are taken in computing the averaged urban turbulent heat fluxes [21]. In the first step, the surface temperatures of roofs, walls, and roads are computed from prognostic equations. These equations contain the effects of sensible and latent heat fluxes, the interactions of the short-and long-wave radiations between the component's surface and air within the urban canyon, and the heat exchange between layers in the components and the temperature within buildings. In the second step, the computed temperatures at the surface of walls and roads and the temperature and humidity at the lowest model level are used to compute the canyon temperature and humidity based on the energy balance equations of sensible and latent energy at the top of urban canyon. Then the computed temperatures and humidity are used to compute the sensible and latent heat fluxes from each urban component in the third step. Finally these fluxes and the geometric parameters of the urban canyon are used to compute the aggregated urban turbulent heat fluxes in the last step. The computed town energy fluxes are included in the bottom boundary condition for the temperature equation in GEM-MACH.
The monthly mean diurnal variation of the traffic heat flux and the heat fluxes (summation of sensible and latent heat flux) at the surface of each urban setting component in the four downtown areas in July are presented in Figure 3. The figure shows that the diurnal variation pattern and magnitudes of the heat fluxes from each urban component in the four urban centers are quite similar in summer. Among the three components, roofs absorb the largest heat flux during the daytime, and walls are the smallest contributor. In the figure, the two peak values of the heat flux from walls correspond to the solar radiations received by the two walls in urban canyon at sunrise and sunset, respectively. Compared to the heat fluxes from the three components, traffic heat fluxes has much smaller magnitude during the daytime. Among the four urban centers, New York has the largest heat flux absorption by roofs and roads. This difference is associated with different urban settings.
In winter, however, roofs are no longer the largest contributor. Walls can absorb much more heat fluxes than roofs before the noontime, and similar heat fluxes after the noontime mainly due to low sun angle in winter. Roads absorb much less heat fluxes in winter. The heat fluxes produced by roads are almost negligible in downtown areas of Toronto and Detroit. In the downtown areas of New York City and Chicago, they are much smaller than the heat fluxes absorbed by roofs and walls during the daytime, Comparison between heat fluxes in wummer and winter suggests that the peak values of the aggregated heat fluxes in winter are just about 30% of the peak values in summer.

Heat Flux Differences between TEB and Non-TEB Simulations
The UHI phenomenon is closely associated with the change of surface heat fluxes by urbanization. To assess the impact of TEB on surface heat fluxes, the monthly mean diurnal variations of heat fluxes from the TEB and non-TEB simulations are computed, and their differences in the four urban centers in July 2015 and January 2016 are presented in Figure 4. It can be seen from the figure that the heat flux differences in the four urban centers in the same month have a similar diurnal variation pattern with a negative value during the daytime peaking around the noontime, and a positive value during the early morning at nighttime in both winter and summer. The differences in summer are similar to the differences between the suburban and rural area shown in [45]. Among the four urban centers, New York City has the largest differences in both summer and winter, reflecting different settings in urban center of New York and in other urban centers.  The comparison between the summer and winter results suggests that during the early morning and at nighttime, the heat fluxes enhanced by urban components are quite similar with about 25 W/m 2 enhancement. However, during the daytime, the reduction of heat flux by urban components are very different in both magnitude and variation pattern. In summer, the magnitude of the reduction is about 150 W/m 2 which is much larger than 10 W/m 2 reduction in winter, and the reduction lasts about 19 h (from 8:00 a.m. to 18:00 p.m.), much longer than 3 h (from 10:00 a.m. to 1:00 p.m.) in winter. The differences are due to the differences in the strength and time period of solar radiation as well as the canyon temperature in the two seasons. It should be pointed out that because the diagnostic equations in the TEB scheme is re-initialized only once, the differences in Figure 5  The heat flux enhancements in the early morning and at nighttime are responsible for The UHI phenomenon, and they are closely associated with the heat storage of urban fabric. The magnitude of heat storage of urban fabric is determined by many factors among which the thermal admittance (Y) is an important one. Y is defined as the square root of the product of thermal conductivity and heat capacity, and is used to measure the surface's ability to absorb/release heat from/to the atmosphere. Other factors affecting the urban heat storage include density, structure, and spatial distribution of buildings, building aspect ratio and canyon aspect ratio, interior temperature in buildings, eddy diffusion in canyons, and moisture availability [7,46] which are taken into account in computing the energy in urban canyon in the TEB scheme. The heat flux enhancements shown in Figure 5 suggest that urban canyons absorb more heat from the atmosphere during the daytime and release the stored heat back into the atmosphere at the nighttime and in the early morning. As the motions of the atmosphere in the lower atmospheric boundary layer (ABL) are driven by surface heat fluxes, the urban effects would have big impacts on urban meteorology and pollutant transport. These impacts are discussed in the following section.

Impacts of the TEB Scheme on Urban Meteorology
Urbanization affects the urban meteorology in different ways. Large friction velocity due to buildings in urban areas enhances downward momentum flux near the surface and affects the urban wind fields. The heat fluxes from urban settings and anthropogenic activities can have big impacts on urban temperature, vertical mixing, static stability of the atmosphere, the ABL height and humidity. Through the coupling between thermodynamic and dynamic processes, the change in horizontal temperature gradient between urban areas and adjacent rural areas or lakes and oceans due to the urban heat island (UHI) effect, can lead to the change of urban meso-scale circulations. Furthermore, the changes of local and non-local vertical diffusivity by the changes in surface urban heat flux, the static stability of the atmosphere and friction velocity can introduce additional impacts on urban temperature, wind and humidity.
The impact of urbanization on the urban meteorology can be measured in different ways. Due to the lack of long term and well covered observations in pre-urbanization period at the same geographic location, the meteorological field differences between the urban area and adjacent rural areas are often used to assess the urbanization impacts. The disadvantage of this approach is that the urban-rural differences may contain the contribution of the effects of different incoming radiative forcing due to different cloud coverage, and different meteorological fields associated with geographic locations and water bodies as well as the downwind effects of the urban air itself. With the numerical model, the impact can be examined by comparing the diurnal variations of monthly mean differences between the TEB and non-TEB simulations at the same geographic locations. Because the comparisons are made at the same location, the extraneous effects involved in the first approach can be avoided, and the magnitude of the urbanization impact can be better assessed. This approach is used in our work to assess the impact of the TEB scheme over the downtown areas of Toronto, New York City, Detroit and Chicago.

Impact on Vertical Diffusivity
Vertical turbulent mixing plays an important role in the dynamic process in the ABL. It can be affected by the changes of urban heat flux and momentum flux. Therefore, understanding the impact on turbulent mixing can help better understand influence of urbanization on dynamic and thermodynamic processes as well as the transport of air pollutants in urban areas.
The diffusion coefficient K is used to describe the strength of the Reynolds stress ρw c produced by unresolved turbulent mixing, where c is the unresolved scalar C. It affects the evolution of C through the diffusion term ∂(K∂C/∂)/∂z. In GEM-MACH, K is derived from the turbulent kinetic energy (TKE)-a prognostic variable in the TKE equation. In the scheme the surface heat flux affects TKE through the ABL convective velocity (w 3 * ) and Obuhkov length used in the lower boundary condition as well as the Richardson number in the statistic stability function [47]. The friction velocity (surface momemtum) affects KTE through the boundary condition at the bottom.
In order to examine the impact of the TEB scheme on the vertical diffusivity, monthly mean diurnal variations of K differences between the TEB and non-TEB simulations are computed. The differences at the lowest model level (24 m above the ground) over the four urban centers in July and January are shown in Figure 5. The figure shows that in July, the diurnal variation patterns and magnitudes of the differences over the four urban centers are similar with a negative value (between 7:00 a.m. to 4:00 p.m.) during the daytime peaking at 8:00 a.m. and a positive value in the early morning and at nighttime. Among the four urban centers, New York City has the largest differences during the daytime due to big surface heat flux reduction ( Figure 4). In January, the diurnal variation patterns over the four centers are similar. However, the magnitudes of the differences over the three U.S. urban centers are very different from that over Toronto. The former ones are just half of the latter in the early morning and at nighttime, but much larger magnitudes during the daytime. In addition, the negative differences over the three U.S. urban centers last much longer.
The comparison of the diffusion coefficient differences between summer and winter suggests that the reductions of K by urbanization in summer and winter are different in both magnitude and variation pattern. In summer, the magnitude of the reduction over the four urban centers during the daytime is about 2 m 2 /s which is larger than 1 m 2 /s in winter over the three U.S. urban centers, and much larger than 0.2 m 2 /s over Toronto. Furthermore, the maximum reduction occurs around 8:00 a.m. in summer which is much earlier than 12:00 a.m. in winter. In addition to the differences during the daytime, K enhancement in the early morning and at nighttime over Toronto is twice as much as that in summer. All these differences are due to the very different diurnal variation pattern and magnitude of the heat flux difference between summer and winter shown in Figure 4.
In GEM-MACH, K at the lowest momentum level (about half level above the lowest model level) is used to compute the time tendency of temperature and chemical species at the lowest model level. Therefore the K differences at this model level need to be examined. The monthly mean diurnal variation of this difference over the four urban centers is presented in Figure 6. The comparison between Figures 5 and 6 shows that the differences at the two levels are very different. While the winter differences have negative values during the daytime at the lowest model level, they are positive over 24 h at the lowest momentum level. In summer, the differences have negative values over the four urban centers at the lowest model level, they have big negative values only over new York City and extremely small negative values over Toronto. It will be shown that these features can help understand the sensitivity of concentration of chemical species to urbanization.
The impact of urbanization on vertical diffusivity is not limited to the levels near the surface. The variation of the impact with height is discussed in Section S1 in the Supplementary Information (SI).

Impact On Temperature
The impact of the TEB scheme on urban temperature can be measured at different height including the level within the canopy layer, at top of canopy layer and within the urban boundary layer (UBL0. The temperature differences at these levels are used to measure the canopy, surface and UBL UHI effects [4]. In the following we will examine the temperature differences at the lowest model level and at at 1.5 m. The monthly mean diurnal variations of the temperature differences between TEB and non-TEB simulations at the lowest model level are shown in Figure 7. In the figure, about 1 • C to 1.5 • C temperature enhancement by the TEB scheme appears in the early morning until 8:00 a.m. and at nighttime over the four urban centers in both summer and winter. Corresponding to the reduction of surface heat flux, the temperature difference drops quickly from 8:00 a.m. until 10:00 a.m. in summer and 12:00 a.m. in winter, and then gradually decreases. The difference remains positive over the four urban centers in winter and also in summer except over New York City where the temperature difference is negative between 8:00 a.m. and 4:00 p.m. The positive temperature differences during the daytime suggest that urbanization can lead to the UHI effect not only in the early morning and at nighttime, but also during the daytime when the urban surface heat fluxes are reduced by the TEB scheme ( Figure 4). The comparison of the summer and winter cases suggests that in the early morning and at nighttime the UHI effect is stronger in summer than in winter, but weaker during the daytime except over the urban center of Detroit.
Due to the use of the same initial and lateral and upper boundary conditions, absence of direct feedback of chemistry to solar radiation in GEM-MACH and negligible humidity differences above 400 m, the incoming solar radiation at the lowest model level is the same in the TEB and non-TEB simulations. Thus urbanization affects urban temperature by changing the surface heat flux, vertical mixing and wind. Because the temperature difference between the two adjacent model grids (2.5-km) is less than 0.2 • C, and the change of wind speed by the TEB scheme is about 0.4 (m/s) in average (see Section 4.3.1), the change of the temperature advection over one hour due to the wind change is very small (0.08 • C) The contributions of the changes in surface heat flux and diffusivity are discussed in the following part.
In order to understand the mechanism underlying the impact of the change of surface heat fluxes and vertical mixing on urban temperature, Ren and Stroud [35] employed a theoretical model to show explicitly how the changes in urban surface flux and vertical diffusivity as well as the reference state affect urban temperature. In a semi-infinite domain, the analytical solution with a time-varying K = K 0 η(t), where T is a dimensionless variable, to the simple model shows that the temperature response ∆T to the change of surface heat fluxes and vertical mixing can be described by the following equation where T = t τ η(t )dt . ∆K and ∆J are the differences of diffusion coefficient and surface heat fluxes between the TEB and non-TEB simulations shown in Figures 5 and 7, respectively. Equation (2) shows that the Green function (G) decreases as height increases, and is small/big when diffusivity is weak/strong in general. The work of Ren and Stroud shows that ∂G/∂z is negative below 80 m. The temperature differences shown in Figure 7 are produced by the combined effects of ∆K and ∆J. Ren and Craig examine the contribution of each effect based on Equation (1) with the typical values of ∆K, ∂T noTEB /∂z and K. The results show that because ∆K is positive in the early morning and at nighttime, and the vertical gradient of the reference temperature is also positive due to temperature reversion, the enhanced ∆K by the TEB scheme leads to warm effect. However, the temperature enhancement by ∆K is about 0.2 • C due to weak temperature difference between the lowest model level and the level above. The temperature increments produced by positive ∆J during the period have similar magnitude to that shown in Figure 7. Therefore, ∆J is the dominant contributor to the UHI phenomenon. The results also show that the impact of ∆K on ∆T is also small during the daytime.
The small impact can be also understood by examining the temperature differences between the noontime and 16:00 p.m. During this period, ∆K at the lowest momentum level is positive except over New York City and increases quickly in both summer and winter. According to Equation (1), the positive ∆K would lead to temperature decrease because the gradients of T noTEB and G are negative. However the temperature differences during this period increase rather than decrease as ∆K increases. This suggests that ∆J is a dominant contributor to the variation of ∆T during the daytime.
According to Equation (1) the impact of ∆J on urban temperature is modulated by G. This modulation effect can be applied to explain the different temperature response in the early morning and during the daytime in summer. Equation (1) shows that the drop of ∆T during the daytime from its value in the early morning is a consequence of reduction of ∆J. Although the magnitude of the surface heat flux reduction during the daytime in summer is about seven times larger than the magnitude of heat flux enhancement during the early morning and nighttime (Figure 4), the corresponding temperature drop shown in Figure 7 has the same magnitude as that of temperature enhancement. Because the diffusivity during the daytime is about 20 times larger than the diffusivity in the early morning, the strong diffusivity would damp the impact of ∆J through G according to Equation (1), and this explains the disproportionate temperature response to the surface heat flux changes during the daytime.
The modulation effect on the impact of ∆J on ∆T through G can be also applied to explain the different temperature responses to the surface heat flux changes in summer and winter. Although Figure 5 shows that in the early morning, the surface heat flux enhancements in the downtown areas of Toronto and Detroit, are stronger in winter, the corresponding temperature enhancements shown in Figure 7, however, are about 0.5 • C lower. According to Equation (1), the weaker temperature responses are due to a stronger diffusivity in winter (not shown). During the daytime, the surface heat flux reductions by the TEB scheme are about seven times stronger in summer than those in winter. Furthermore, the corresponding temperature drops in winter and summer are similar over the downtown areas of Toronto, Chicago and Detroit, and have about 1 • C difference over New York City. The strong damping effect by G associated with strong diffusivity in summer is responsible for the weaker temperature responses.
Equation (1) shows that the temperature response to ∆J is the convolution of G 0 and ∆J. This suggests that the positive temperature response in the early morning can impact the temperature response during the daytime. This impact and the strong damping effect by diffusivity explain the positive ∆T in winter over the four urban centers and in summer over the three urban centers during the daytime when ∆J is negative.
The fact that the temperature at a given time is the convolution of G 0 and ∆J can be also employed to understand the different evolution pattern of temperature ( Figure 7) and heat flux differences ( Figure 5). While Figure 4 shows that corresponding to the diurnal variation of radiation, the evolution patters of heat flux differences are nearly symmetric around noontime, such symmetry, however, does not exist in the diurnal variation of temperature differences in summer. The variation of the temperature is much slower in the afternoon than in the morning. The convolution of G 0 and ∆J destroys such symmetry even though the evolution of diffusion coefficient is symmetric around noontime.
In GEM, the lowest model level over downtown Toronto is at around 24 m, and the temperature at the screen-level (1.5 m) is calculated diagnostically. The impact of the TEB scheme on the temperature at this level has similar diurnal variation pattern as that at the lowest model level but with larger magnitudes (2.5 • C).
The impact of urbanization on temperature is not limited to the lowest model level, and the variation of the impact with height is discussed in Section S2.

Impact on Humidity
In the TEB scheme, the impact of roofs and roads on urban hydrology is considered in the evolution equation of water-reservoir [21], and the turbulent fluxes of heat and humidity are calculated in a similar way. Due to the small impact of the change of wind, the urbanization affects the urban humidity through surface heat flux and the vertical mixing. Because the surface of urban components in the TEB scheme is impervious and the liquid precipitation intercepted by roads and roofs goes into underground sewer systems, the surface humidity fluxes in the TEB simulation should be much weaker than that in the non-TEB simulation in urban areas.
In order to examine the impact of the TEB scheme on humidity, the diurnal variations of the monthly mean urban relative humidity difference over the four urban centers are computed, and the differences at the lowest model level are presented in Figure 8. It can be seen from the figure that urbanization leads to the reductions of humidity in both summer and winter, and the reduction is larger in winter than in summer. The variation of the differences in summer is smooth except over New York City. In winter, the TEB scheme leads to continuous drop of relative humidity except over Toronto.
The impact on humidity can lead to not only the change of latent heat flux, but also changes in the uptake of water onto hydroscopic particles, as well as chemical reactions such as O1D Specific humidity is another variable measuring the humidity level in the atmosphere. Comparison between the TEB and non-TEB simulations shows that the TEB scheme also leads to the reduction of specific humidity over the four urban centers. In winter, the reduction at the lowest model level is about 2% of the specific humidity in the non-TEB simulation between 8:00 p.m. to 8:00 a.m. over the four urban centers, and 14% over New York , 7% over other centers between 10:00 a.m. to 7:00 p.m. In summer, the reduction is also about 2% of the specific humidity in the non-TEB simulation between 8:00 p.m. to 8:00 a.m. over the four centers, and is about 3% over New York, 6% over Toronto, and 10% over Chicago and Detroit. The TEB impacts on relative humidity and specific humidity decay rapidly with height, and they become negligible above 400 m.

Impact on the ABL Height
Because the top of the ABL separates the boundary layer and free troposphere, the ABL height describes the vertical extent of mixing process generated within the ABL, and thus imposes constraints on the exchange of momentum, energy and chemical species between the ABL and free atmosphere [48]. Studies show that the ABL height has a strong negative correlation with air pollution in big cities [49,50]. In GEM-MACH, the ABL height is involved in the TKE scheme through w * in the surface boundary condition in the unstable case [47]. It is also utilized in two places: in the calculation of the plume-rise height for anthropogenic major point source emissions and forest fire emissions, and in assigning the minimum of diffusion coefficient within the ABL.
In numerical models, the ABL height is not a prognostic variable. Different methods have been used to calculate the ABL height, and each method produces quite different results [51]. Among these methods, potential temperature (PT)-based method and the Richardson number (R i )-based methods are two popular methods. Because the PT-based method is used to compute the ABL height in GEM-MACH, the differences in vertical temperature profiles between the TEB and non-TEB simulations suggest that the TEB model would have a great impact on the ABL height.
The monthly mean diurnal variations of the ABL height differences between the TEB and non-TEB simulations are shown in Figure 9. It can be seen from the figure that the TEB scheme has similar impact on the ABL height in summer over the four urban centers. Large positive difference drops rapidly from 2:00 a.m. and becomes negative around 8:00 a.m. It then starts to increase until 4:00 p.m. and has very small changes between 4:00 p.m. and 8:00 p.m. Rapid increase occurs between 8:00 p.m. and 2:00 a.m. The diurnal variations of the differences in winter are different from those in summer. The positive differences in the early morning are much smaller and variations are much slower than those in summer. Rapid drops occur at around 8:00 a.m. and continue until noontime. The drops are much deeper over New York City and Detroit than over Toronto and Chicago. While the differences over New York City and Detroit start to increase quickly after noontime and reach maximum value around 7:00 p.m., they remain around zero over Toronto and Chicago until 4:00 p.m. The enhancement of the ABL height in the late afternoon by the TEB scheme is much larger in winter than in summer.

Toronto
To show the sensitivity of the ABL height calculated using the R i -based method, ABL height differences based on this method are computed offline, and their diurnal variations are shown in Figure 10. In the figure, the differences are zero at 8:00 p.m. due to re-initialization. The figure shows that in summer, the differences are quite similar between 4:00 p.m. to 8:00 a.m. over the four urban areas. During this period, the TEB scheme leads to the ABL height enhancement. The enhancements increase rapidly from 4:00 p.m. until 10:00 p.m., then drop quickly between 10:00 p.m. and 11:00 p.m. They continue to drop after 11:00 p.m. but with a much slower rate until 7:00 a.m. and then increase slightly until 8:00 a.m. Between 8:00 a.m. and 4:00 p.m. the difference over New York City is different from those over the other three urban areas. During this period, while the TEB scheme leads to as much as 100 m reduction over New York City, but small enhancement in other three urban areas (except a small reduction around 9:00 a.m. over Toronto). In winter, the TEB scheme leads to the enhancement of the ABL height over the whole day over the four urban centers. The differences over the four areas are very similar with a large enhancement in the afternoon. The evolution patterns between 1:00 a.m. and 6:00 p.m. are similar to those in summer except over New York City. The comparison between Figures 9 and 10 suggests that the sensitivities of the ABL height based on the two approaches are quite different. In summer, while the large PT-based ABL height enhancements appear in the early morning, and large reductions appear in the morning over the four urban areas, the large R i -based enhancements appear at the nighttime, and large reductions appear only over New York City. In winter, while the large PT-based ABL height enhancements appear between 8:00 p.m. to 10:00 p.m., and reductions appear in the early afternoon over the four urban areas, the large R i -based enhancements appear in the late afternoon, and there is no reduction.

Impact on Urban Wind Fields
Urbanization affects wind in the ABL by exerting a dragging force by buildings, by changing the horizontal pressure gradients and vertical diffusivity. In the following, the impact of the TEB scheme on wind will be examined by comparing wind speed from simulations with and without the TEB schemes. As the four urban centers are near big water bodies, the impact on lake breezes will also be examined.

Differences in Wind Speed
The monthly mean diurnal variations of wind speed differences at the lowest model level over the four urban centers are shown in Figure 11. In summer, The evolution patterns and the magnitudes of the differences over the four centers are different. Among the four centers, Toronto has the largest magnitude (0.4 m/s) of the difference in the early morning. In winter, the evolution pattern of the differences over the four urban centers are similar between 8:00 a.m. and 2:00 p.m. with a maximum value around noontime. While the differences have positive and negative values in summer, they are positive most of the time.

Differences in Lake Breezes
Lake breezes are driven by the heat contrast between land and adjacent water. As the urbanization changes the energy balance over land, it modulates the strength and direction of lake breezes, and consequentially affects the transport of pollutants. Because the four urban centers are close to big water bodies, the numerical simulation results with and without TEB can be employed to investigate the TEB effect on lake breezes.
The impact of the TEB scheme on the land-to-lake circulation can be understood by examining the height cross sections of the monthly mean wind differences between TEB and non-TEB simulations at four hours in July. The simulation results show that at the lower levels, the TEB effect weakens the city-to-lake breezes in the early morning and at nighttime due to the surface heat flux enhancement, and lake-to-city breezes during the daytime due to the surface heat flux reduction in the downtown area. Opposite effects can be identified few levels above. The weakened city-to-lake circulation in the early morning and at nighttime should have an impact on the transport of urban pollutants between urban areas and water bodies, and the implications of the change of the local circulation for the transport of urban pollutants in Toronto have been discussed in [41].

Impacts on Urban Pollutants and AQHI
Air pollutants have adverse effects on human health and cause damages on trees, vegetables and buildings. Their impacts are measured by the Air Quality Health index (AQHI) developed by the Environment and Climate Change Canada (https://www.canada.ca/en/environment-climatechange/services/air-quality-health-index/, accessed on 10 February 2020). Due to high population density and high concentration of air pollutants in the urban areas, the impacts of urbanization on urban chemical environment have been widely studied. In the following part the impacts of the TEB scheme on the evolution of some major urban pollutants and AQHI are examined by comparing the TEB and non-TEB simulation results.
Several factors including boundary conditions, short-wave solar radiation and meteorological conditions can affect the regional model simulations of major air pollutants. In the TEB and non-TEB simulations, the same chemical and meteorological lateral and upper boundary conditions are used. Although the same surface emissions are used in TEB and non-TEB simulations, the dry deposition velocities can be different over urban areas due to different surface momentum flux. Therefore, surface chemical boundary conditions can be different in the two simulations. Due to the lack of direct feedback of chemistry to solar radiation and extremely small differences in humidity above 400 m, the incoming solar radiation should be roughly the same in the two simulations. However, the reflected solar radiations in the two simulations over urban areas can be different due to different albedo. Thus changes in meteorological fields, chemical boundary conditions at bottom, chemical reactions, and reflected solar radiations can lead to differences in the concentration of air pollutants between the two simulations. These differences are examined in this section.

Impact on Carbon Monoxide
Among the major air pollutants, carbon monoxide (CO) can be taken as a passive tracer. The concentration differences between the TEB and non-TEB simulations can be attributed to the change of meteorological fields due to the TEB effects. The TEB impacts on the CO mixing ratio can be seen from Figure 12 which shows the weekday mean diurnal variations of the mixing ratio differences between the TEB and non-TEB simulations. The figure shows that in summer, the diurnal variation of CO mixing ratio differences over the four urban centers are quite similar. The TEB scheme leads to 50 µg/kg mixing ratio decrease in the early morning and 100 µg/kg (about 28%) decrease in the late afternoon. Between 8:00 a.m. and noontime, the TEB scheme has a small impact. In winter, the diurnal variation patterns of the differences are also similar over the four urban centers with two drops of the differences: one in the early morning between 0:00 a.m. and 4:00 a.m. and another at around noontime. However, the magnitudes of the two drops are very different over different urban centers. Over New York City, the first drop is more than 400 µg/kg (25%), four times larger than that over Detroit, and the second drop is about 600 µg/kg (50%), which is five times larger that over Detroit.
The figure also shows that the TEB has very different impacts in summer and winter. In summer, the differences over the four urban centers vary slowly in the early morning, but vary rapidly in winter, especially over Toronto and New York City. During the daytime, big difference appears in the late afternoon in summer, but appears around 1:00 p.m. in winter. Overall, the TEB scheme leads to much larger mixing ratio decrease in winter time than in summer time. The changes of both wind and vertical mixing by the TEB scheme can contribute to the mixing ratio difference shown in Figure 12. Although their contributions cannot be separated, the contribution of the changes in wind can be assessed based on the magnitude of ∆U · ∇C, where ∇C is the horizontal gradient of mixing ratio. Because the magnitude of this term is less than 30 µg/kg/h in the early morning over Toronto and New York City, it cannot account for 100 µg/kg/h drops. Therefore the contribution associated with the change of vertical mixing is dominant.
Because the impact of diffusivity on tracer concentration can also be described by the 1-D diffusion model, Equation (1) can be applied to interpret the numerical results shown in Figure 12. Because the same boundary and initial conditions are used in the TEB and non-TEB simulations, the differences in mixing ratios between the two simulations in the 1-D model can be written as Unlike the temperature gradient which is positive in the early morning, the gradient of the CO concentration is always negative. Thus according to Equation (4), positive/negative ∆K at the lowest momentum level would lead to the decrease/increase of mixing ratio (because the gradient of G is negative). Comparison between Figures 6 and 12 suggests that it is indeed the case. Positive ∆K in winter corresponds to negative ∆C over the four urban areas over 23 h. In summer, positive ∆K also corresponds to negative ∆C. Furthermore, the negative ∆K between 8:00 a.m. and 12:00 a.m. in summer over New York City and Toronto corresponds to (small) positive ∆C. According to Equation (4), the impact of change of vertical diffusivity on tracers is modulated by the vertical gradients of mixing ratio of the non-TEB simulation and G. For a given ∆K, large/small magnitude of the vertical gradients of C nonTEB and G would enhance/weaken the impact. To examine the modulation effects, CO mixing ratio difference between the lowest model level and one level above in the non-TEB simulation (∆C CO,2−1 ) is used to represent the vertical gradient of mixing ratio. Because the TEB scheme has a larger impact on mixing ratio over urban centers of New York City and Toronto, the diurnal variations of (∆C CO,2−1 ) over the two urban centers are presented in Figure 13 along with the mixing ratio differences between the TEB and non-TEB simulations (∆C CO ) at the lowest model level. It can be seen from the figure that except in the early morning in summer, the diurnal variation patterns of ∆C CO and ∆C CO,2−1 are quite similar. The correlation between ∆C CO,2−1 and ∆C CO is more robust in the afternoon than in the morning due to larger ∆K, and the correlation in winter is more robust than in summer. The modulation of the vertical gradient of G can be identified by comparing the variation of δC CO during the daytime in summer. According to Equation (4) the negative ∆K shown in Figure 6 over New York City and Toronto between 8:00 a.m. and 12:00 a.m. in summer would enhance mixing ratio. Although the magnitude of ∆K over New York City is much larger than that in Toronto, the enhancements over the two urban centers are quite similar (Figure 12). Because the values of ∆C CO,2−1 over the two centers are similar, the similar mixing ratio enhancements should be due to a weaker vertical gradient of G. This weaker gradient is caused by a stronger diffusivity over New York City which is much stronger than the diffusivity over Toronto. The numerical results show that the TEB scheme has similar impacts on PM 2.5 .

Impact on NO 2 , NO x and O 3 (g)
Because NO 2 and O 3 (g) are involved in photochemical reactions, they cannot be taken as a passive tracers. Thus both ∆K and the change of chemical reactions due to the change of mixing ratio by ∆K can affect the mixing ratio of the two species. Their impacts can be seen in Figure 14 which shows the diurnal variations of the weekday mean NO 2 mixing ratio differences (∆C NO 2 ) in summer and winter. The comparison between this figure and Figure 12 suggests that the evolution patterns of the mixing ratio differences of CO and NO 2 are quite similar over the four urban centers, suggesting that ∆K plays a major role in affecting the mixing ratio. Numerical results also show that in summer the variation of mixing ratio of NO 2 correlates closely with the vertical gradient except in the early morning (between 0:00 a.m. to 2:00 a.m.). In winter, strong positive correlation can be identified except between 1:00 p.m. and 6:00 p.m. over Toronto and between noontime and 4:00 p.m. over New York City during which the correlations are negative (not shown). O 3 (g) is produced by the photochemical reactions of nitrogen oxides (NO 2 ) and volatile organic components (VOCs) during the daytime, and is destructed by reacting with NO to form NO 2 . Therefore the diurnal variation of O 3 (g) difference (∆C O 3 (g) ) should be opposite to the variation of the NO x (summation of NO and NO 2 ) difference as shown in Figure 15. In summer, the variations of the O 3 (g) differences correlates negatively with the variations of the NO x difference. However, the magnitudes of the two differences are different. In the early morning, the magnitude of NO x can be 10 µg/kg more than that of ∆C O 3 (g) , and about 5 µg/kg at nighttime. In winter, however, the magnitudes of ∆C O 3 (g) are much smaller than those of ∆C NO x except in the afternoon over urban center of Detroit, and ∆C O 3 (g) has no response to the strong variation of ∆C NO x in the morning and at nighttime. Numerical results show that the modulations of ∆C O 3 (g),2−1 and ∆C O 3 (g),2−1 to ∆C O 3 (g) and ∆C O 3 (g) are similar to the modulation of ∆NO 2,2−1 to ∆NO 2 .

Impact on AQHI
AQHI is a scale designed by ECCC and Health Canada to measure the health impact of short-term exposure to air pollution. The concentrations of three major pollutants including ground-level ozone, NO 2 and PM 2.5 are used to calculate AQHI through the following formula Because the value of AQHI increases exponentially with the increase of the concentrations of the three pollutants, high value of AQHI means great health risk associated with the air quality. Because AQHI is a nonlinear function of the concentrations of the three pollutants, the monthly mean difference of AQHI is not the linear combination of the monthly mean differences of the three pollutants.
To show the impact of the TEB scheme on AQHI, monthly mean hourly AQHIs in July and January are computed, and the weekday mean diurnal variations of AQHI and its components (AQ NO 2 , AQ O 3 and AQ AF ) are presented in Figures 16 and 17 for July and January, respectively. It can be seen from the two figures that while the TEB effect leads to the reduction of AQHI in summer between late afternoon and early morning in summer, it leads to the reduction in the whole day in winter. The impacts of the TEB scheme on AQHI are stronger in summer than in winter except over New York City where the reduction of AQHI can be as big as 1. The impacts are similar over the four urban centers in summer. and over Toronto, Chicago and Detroit in winter. It can been seen from Figures 16 and 17 that in summer, the enhanced AQHI by O 3 (g) enhancement in the early morning and at nighttime is offset by the reduced AQHI by the reductions of PM 2.5 and NO 2 . In winter, the contribution of PM 2.5 is extremely small and the reduction of AQHI is entirely due to the contribution of the reduction of NO 2 . In both summer and winter, the reduction of NO 2 by the TEB scheme plays a dominant role in reducing AQHI.

Toronto
The numerical results also show that the reduction of AQHI in summer is much larger in weekend than in weekdays between late afternoon and early morning, but it is smaller in winter during this period (not shown).

Evaluation of the Urbanization Impacts on Meteorological and Chemical Predictions
The purpose of implementing the TEB scheme in GEM-MACH is to improve the accuracy of numerical prediction of urban meteorological and chemical fields. While the results in preceding sections show robust responses of both meteorological and chemical fields to the TEB scheme in urban areas, and the features of the responses can be interpreted by analytical results, these impacts still need to be evaluated against observations to justify the improvements by the TEB scheme. In this section, the ground based observations are used to evaluate the model simulations of temperature, humidity, wind and major pollutants in the GTA in July 2015.

Meteorological Observations
Observations used for the evaluation are from Mesonet Automated Land Stations Observations, a high-resolution monitoring system across Southern Ontario developed by the ECCC by adding 40 compact stations, 10 ATMOS stations, and three standard automated MSC Auto8 stations to the existing observation networks during the Pan-, and Para-pan American Games period in 2015.
The minute-by-minute observations of temperature, wind and humidity in July are averaged at the beginning of each hour with a 5-min time window. Based on their locations, the averaged observations are then grouped into six categories including downtown, uptown, water front, island, suburban and rural areas in the GTA. The grouped data are employed to evaluate both TEB and non-TEB model simulations, and to assess the TEB impact on urban meteorological fields.

Observations of Chemical Species
Data from the U.S. EPA AirNow program (www.AirNow.gov, accessed on 12 March 2019) contains real-time air quality observations from over 2000 monitoring stations and forecasts for more than 300 cities across the United States, Canada, and Mexico. The AirNow data in the GTA is grouped into five categories including downtown, uptown, lakeshore, suburban and industrial areas within the GTA.

Statistics of P-O
To evaluate the temperature at 1.5 m from TEB and non-TEB simulations, the P-O statistics of temperature and relative humidity including bias, root-mean-square deviation (RMSD) and correlation coefficient are computed for the whole month of July 2015 for each group of observations and presented in Figure 18. The figure shows that both TEB and non-TEB simulations lead to negative temperature biases except over rural area where the bias of the TEB simulation is slightly positive suggesting that both simulations under-predict the urban temperature. The TEB scheme leads to significant improvements in temperature biases in downtown, uptown and suburban areas. The temperature biases of the TEB simulation in these areas are just about 20% of those of the non-TEB simulation. The TEB scheme also leads to the improvement of RMSD of temperature in all the areas except the rural area. The RMSD in the downtown area is about 50% of the RMSD of the non-TEB simulation. However only small improvement in correlation coefficient by the TEB scheme can be found in the downtown area, and there is virtually no impact in other areas.
Significant improvement of relative humidity bias by the TEB simulation can be seen in downtown, uptown and water-front areas, and small improvement can be found over island and suburban area. However, the TEB scheme produces slightly worse bias over the rural area. Except in the rural area, the non-TEB simulation produces positive biases suggesting that relative humidity is over-predicated in general. The TEB-simulation however, under-predicts the relative humidity in downtown, suburban and rural areas. RMSD is significantly improved by the TEB simulation in the downtown area and noticeable improvement can be found in uptown and island areas. Improvement in correlation coefficient of temperature and relative humidity by the TEB scheme is found only over the downtown area.
The statistics of P-O of wind speed shows that the TEB scheme has virtually no impact on the biases, RMSDs and correlation coefficients over the all areas. The results show that while both simulations produce a smaller bias (−1.5 m/s) and RMSD (2 m/s) over downtown area, they produce big negative biases (−4-−5 m/s) and RMSD (5 m/s) over other areas. They also produce smaller correlation coefficient (0.4) over downtown area than other areas (not shown).

Diurnal Variation of P-O Statistics
The results in Section 4 show that there is a strong diurnal variations of sensitivities of both meteorological and chemical fields to the TEB scheme, and large impacts are found in the early morning and at nighttime. These time-dependent sensitivities can be evaluated by examining the P-O statistics at each forecast hour. It should be pointed out that due to different averaging process, the mean values of the scores at each forecast hour are not identical to the scores shown in preceding subsection.
The monthly mean diurnal variations of P-O bias and RMSD of temperature over the six areas are shown in Figure 19. It can be seen from the figure that the TEB scheme leads to the persistent improvement in bias in all areas except the rural area between 8:00 p.m. to 8:00 a.m. Significant improvement can be seen in downtown and uptown areas and noticeable improvement in other areas. During this period, temperature biases of TEB and non-TEB simulations are negative and have very small values over the rural area.
The diurnal variations of RMSD of temperature of TEB and non-TEB simulations shown in Figure 19 suggest that the TEB scheme leads to significant improvement on RMSD in nearly the whole day in the downtown area, big improvement between 8:00 p.m. and 8:00 a.m. and noticeable improvement at other times in uptown, and noticeable improvement only between 8:00 p.m. and 8:00 a.m. over in areas except the rural area. The temperature correlation coefficients is improved by the TEB scheme only over the downtown area (not shown). For relative humidity, significant improvements in bias by the TEB scheme can be found in downtown and uptown areas between 8:00 p.m. to 6:00 a.m. (Figure 20). Noticeable improvements during this period can be also found in island and waterfront areas. Outside of this period, the TEB scheme produces larger biases in the downtown area, but smaller biases in uptown and island areas. Biases in the suburban area show some improvement only between 8:00 p.m. and 1:00 a.m. Figure 20 also shows that the TEB scheme leads to significant improvement of RMSD of relative humidity in the downtown area and big improvement in uptown and island areas over 24 h. Small improvements can be found in suburban and waterfront areas. While there is significant improvement of correlation coefficient over 24 h in the downtown area by the TEB scheme, no improvement is found in other areas (not shown).  Figure 19 but for relative humidity.

Evaluation of Predictions of Air Pollutants
In order to evaluate the model predictions of chemical species, observations from AirNow in July 2015 are used to evaluate both TEB and non-TEB simulation results of O 3 , NO 2 , NO x and PM 2.5 over downtown, uptown, suburban, lakeshore and industrial areas.

P-O Statistics
To show the overall evaluation results, the P-O biases, RMSDs and correlation coefficients of NO 2 , NO x and O 3 are computed over the whole month of July, 2015 for each categorized area. The results are presented in Figure 21. For NO 2 , Figure 21 shows both TEB and non-TEB simulations produce positive bias except over uptown area suggesting that both simulations over-predicate the NO 2 mixing ratio. The non-TEB simulation produces more than 4 ppb biases over the downtown and industrial areas, and the inclusion of the TEB scheme improves the bias of NO 2 over all the areas, particularly over downtown and uptown areas where bias is reduce by 50% and 90%, respectively. The figure also shows that the non-TEB simulation produces more than 4 ppb RMSDs over the downtown, industrial and lakeshore areas, and the TEB scheme leads to a 40% improvement of RMSD over the downtown area and noticeable improvements over other areas. However it does not lead to any improvement on correlation coefficients. Big positive biases of NO x produced by both TEB and non-TEB simulations over downtown and industrial areas can be identified in the Figure 21. However, the biases over other areas are very small. This is different from the biases of NO 2 which have big values over lakeshore and suburban areas. The figure shows that the TEB scheme can greatly reduce the biases over all the areas except the lakeshore area. It leads to a 30% reduction of RMSDs over downtown and 24% over industrial areas but it has virtually no impact on correlation coefficients.
The variation of O 3 (g) is opposite to the variation of NO 2 . Therefore biases of O 3 (g) and NO 2 should have different signs. Statistics shown in Figure 21 suggests that it is indeed the case. The figure shows that both TEB and non-TEB simulations under-predict O 3 (g) mixing ratio over all the areas with big negative biases over downtown, uptown and industrial areas. Although the TEB scheme improves the negative biases over downtown and uptown areas, and produce small improvement over other areas, the improvements are smaller than those of NO 2 . The TEB scheme also leads to a reduction of RMSD over the downtown area and small reduction over uptown and industrial areas. However, it has virtually no impact on correlation coefficients. P-O statistics of PM 2.5 is discussed in Section S3.

Diurnal Variation of P-O Statistics
The impact of the TEB scheme on the prediction of chemical tracers can be better understood by examining the diurnal variation of P-O statistics. The monthly mean diurnal variation of P-O bias of NO 2 (left panel in Figure 22) shows that NO 2 mixing ratios are over predicated in both TEB and non-TEB simulations from 8:00 p.m. to 10:00 a.m. over all the areas except the uptown area where NO 2 mixing ratio is over-predicated in the non-TEB simulation between 8:00 p.m. to 10:00 a.m., and between 2:00 a.m. to 10:00 a.m. in TEB simulation. The TEB scheme leads to small under-predictions between 10:00 a.m. to 6:00 p.m. over the downtown area, and 10:00 a.m. to 11:00 p.m. over the uptown area. A big improvement (>30%) in bias by the TEB scheme can be found over downtown and industrial areas between 8:00 p.m. to 8:00 a.m., and over the uptown area between 2:00 a.m. to 8:00 a.m. Noticeable reduction of bias can be identified between 8:00 p.m. and 8:00 a.m. over the lakeshore area. The figure shows Similar improvement in RMSD by the TEB scheme. However, the scheme leads to no improvement in correlation coefficient. The diurnal variations of P-O biases of NO x of both TEB and non-TEB simulations are similar to those of NO 2 over all the areas except over the lakeshore area where biases are positive from 1:00 a.m. to 10:00 a.m., and are large and negative around the noontime (not shown).

NO2 (P-O) bias
Unlike the positive biases of NO 2 and NO x , Figure 23 shows that the bias of O 3 (g) is negative most of the time. Very small positive biases can be found over uptown and lakeshore areas between noontime and 4:00 p.m. The improvement on biases and RMSDs of O 3 (g) by the TEB scheme can be found between 5:00 p.m. to 6:00 a.m. over the downtown, uptown areas, and between 8:00 p.m. to 6:00 a.m. over the industrial area. Small improvement over lakeshore and suburban appears only between 6:00 p.m. to 10:00 p.m., much shorter than that in other areas. The RMSD is improved by the TEB scheme between 8:00 p.m. to 8:00 a.m. over the downtown and uptown areas, and between 0:00 a.m. to 8:00 a.m. over the industrial area. The TEB scheme leads no improvement in correlation coefficient. The diurnal variations of PM 2.5 biases and RMSD are discussed in Section S3.

Summary and Discussion
The Urban surface and anthropogenic activities have large impacts on weather and air quality. The impacts are investigated in this work by using ECCC's high-resolution air quality forecast model incorporating with the TEB scheme. The sensitivities of the evolution of meteorological and chemical fields to the urban effects over four major North American cities are examined by comparing the TEB and non-TEB simulation results at the same geographic location. The improvement of the predictions due to the inclusion of the urban scheme is evaluated by comparing the model simulation results against the ground based observations. The sensitivities of surface heat fluxes to the TEB scheme is examined by comparing the monthly mean diurnal variations of the fluxes from the TEB and non-TEB simulations on the lowest model level. The comparisons show that over the four urban centers the surface heat fluxes are enhanced in the early morning and at nighttime, but weakened during the daytime by the TEB scheme in both summer and winter. While the magnitudes of the enhancement are similar in summer and winter, the magnitudes of the reduction are much larger in summer than in winter. The changes in heat fluxes lead to the change in the vertical turbulent diffusivity. At the lowest model level, vertical mixing is weakened during the daytime but is enhanced almost throughout the whole day. Above the lowest model level, however, it is enhanced almost in the whole day except during the daytime over New York City in summer. In the lower half of the ABL, the TEB scheme leads to the increases of diffusivity with increasing height.
Corresponding to the diurnal variation of the surface heat flux differences, temperatures at the lowest model level are enhanced by 1 • C and 2 • C at 1.5 m in the early morning and at nighttime in summer and winter over the four urban centers. During the daytime, temperatures over the four urban centers in winter, and over Chicago and Detroit in summer are also enhanced by the TEB scheme, but the temperatures over Toronto and New York City are weakened. These features are found to be associated with the modulation effect of the vertical diffusivity. The results show that the impact on temperature is not limited to the lowest model level, but can reach as high as 500 m at 6:00 p.m. local time. Studies based on the 1-D model suggest that the change of urban surface heat flux play a major role in the UHI phenomenon and the change of urban vertical diffusivity has very small effect on urban temperature.
Because the vertical gradient of temperature is used in the definition of the ABL height, the change of the temperature profile leads to a big impact on the ABL height. The results show that the TEB scheme leads to the increase of the ABL height in the early morning and at nighttime but decrease during the daytime over the four urban centers. The magnitude of the reduction is much larger in summer than in winter. It is found that the impact on the ABL height is sensitive to the definition of the ABL height.
The simulation results show that while the TEB scheme leads to the reduction of relative humidity and specific humidity in both summer and winter, it has a small impact on wind speed. the UHI effect weakens the city-to-lake circulation in the early morning and at nighttime, and consequentially affects the transport of chemical species between the urban centers and adjacent water body.
The TEB scheme is found to have a big impact on the mixing ratio of urban chemical species. In summer it leads to the reduction of mixing ratio of major pollutants except ground-level ozone in the early morning and at nighttime, but very small impact during the daytime. In winter, the TEB scheme leads to the reduction of mixing ratio in the whole day with big reductions in the morning and early afternoon. The sensitivity of ground-level ozone to the TEB scheme is opposite to that of NO 2 . Although the enhanced ground-level ozone leads to increase the value of urban AQHI, the impacts of the TEB scheme on other major pollutants lead to the a bigger reduction of AQHI, and thus improve urban air quantity in the early morning and at nighttime in summer, and in the early morning and early afternoon in winter. Because the same surface emissions are used in the TEB and non-TEB simulations, the mixing ration differences produced by the TEB scheme can be attributed to the change of vertical diffusivity. The analytical results based on the 1-D model and 3-D model simulation results shows that the impact of the change of diffusivity is modulated by the vertical gradients of GF and the non-TEB mixing ratio. Large mixing ratio reductions are due to large vertical gradient of mixing ratio.
Our results show that although the impacts of the TEB scheme on the meteorological and chemical fields over different urban centers share some similarities, the magnitudes of the impacts are different due to different urban settings and the state of non-TEB simulation. Among the four urban centers, New York City has the largest surface heat flux reduction during the daytime in summer. The big reduction leads to not only the strongest cooling effects but also the largest diffusivity reduction among the four urban centers during the daytime. Large vertical gradients of chemical species around the noontime over New York City in winter lead to the largest reductions of the concentration of species.
The sensitivities of the model predictions of meteorological and chemical fields to the TEB scheme over the GTA are evaluated by comparing the prediction-minus-observation (P-O) statistics of TEB and non-TEB simulations. The P-O statistics show that the TEB scheme greatly improves the predictions of both temperature and humidity. It leads to the significant reduction of cold biases of temperature predictions in downtown, uptown and suburban areas as well as the significant reduction of the RMSD of temperature predictions in the downtown area, and big reduction in the uptown area between 8:00 PM and 8:00 AM. During the same period, the biases of humidity predictions over downtown, uptown and waterfront areas are also greatly improved, along with the significant improvement in RMSDs over the downtown area. Small improvement of the correlation coefficient is found only in the downtown area. The TEB scheme does not lead to any improvement in wind speed.
The predictions of some major urban pollutants over the GTA are also greatly improved by the TEB scheme. The improvement of biases and RMSDs of NO 2 prediction can be found over all area between 8:00 PM. and 8:00 AM. The improvement of biases over downtown, uptown and industrial areas are significant. Similar improvement in biases and RMSDs of the PM 2.5 can be identified. Biases and RMSDs of the ground-level ozone predictions are also improved over the all area, particularly over downtown and uptown areas during the same period. All these evaluation results show that the TEB scheme can indeed improve the predictions of urban meteorological and chemical fields.
Both analytical and numerical results show that the vertical diffusivity has a great modulation effect on the responses of temperature and chemical species to the TEB scheme. The vertical mixing in current version of GEM is calculated based on the turbulent kinetic energy scheme which is different from the non-local scheme [52]. Because the improvement of urban air quality is very sensitive to the change of vertical mixing, it would be interesting to see in future work if a different mixing scheme can further improve the urban predictions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/9/969/ s1, Figure S1: Vertical profiles of the monthly mean diffusion coefficient differences between the TEB and non-TEB simulations over July at 0 6, 12, 18 h (local time) in the downtown areas of Toronto (a), Chicago (b), New York City (c), and Detroit (d). The heights of the first and the 14th level over the Toronto urban center are about 24 m and 2766 m, respectively, Figure S2: Same as Figure S1 but for temperature., Figure S3: Monthly mean and hourly mean PM 2.5 P-O biases, root-mean-square deviation, and correlation coefficients of simulations of TEB (orange) and non-TEB (blue), Figure S4