Impact of Physics Parameterizations on High-Resolution Air Quality Simulations over the Paris Region

Lei Jiang 1,2,*, Bertrand Bessagnet 1,† , Frederik Meleux 1, Frederic Tognet 1 and Florian Couvidat 1 1 INERIS, National Institute for Industrial Environment and Risks, Parc Technologique ALATA, 60550 Verneuil-en-Halatte, France; bertrand.bessagnet@citepa.org (B.B.); Frederik.MELEUX@ineris.fr (F.M.); Frederic.TOGNET@ineris.fr (F.T.); Florian.COUVIDAT@ineris.fr (F.C.) 2 Sorbonne University, UPMC-Paris 06, DE129, 75005 Paris, France * Correspondence: lei.jiang@ineris.fr † Now at Citepa: Technical Reference Center for Air Pollution and Climate Change, 42, rue de Paradis, 75010 Paris, France.


Introduction
Nowadays, around 55% of the world's population lives in urban areas, and this number is expected to increase by 68% by 2050. The European Environment Agency annual report has pointed out around 25% of the European urban population are exposed to air quality exceeding the European Union air quality standards [1], and air pollution is the leading preventable risk factors for premature death in Europe, being responsible for 400,000 deaths per year directly or indirectly. In addition, built-up surfaces that are mainly composed of artificial buildings and cement pavement which are clearly distinguished from natural surfaces favor heat accumulation [2,3]. Recent research shows that in urban regions, Urban Heat Islands (UHIs) can be observed by a difference of about 1.5-2 • C compared to rural area temperatures, and about 2.0-2.5 • C difference can be observed over densely built-up areas [4]. The UHI has a negative effect on the planetary boundary layer (PBL) height [5] and surface wind speed, which affects the transport and dispersion of pollutants. Overall, the UHI changes the local atmospheric dynamic [6], significantly influences local air pollution, reduces visibility, increases water usage, and enhances heat-related morbidity [7,8]. In modeling studies, detailed information on urban parameters is critical for the simulation of the UHI effect. Studies conducted in Taiwan have found that the urban canopy model can improve the prediction of UHI intensity, boundary layer development, land-sea breeze [9], and precipitation [5]. Air pollution prediction is crucial in order to inform a large fraction of the population and take adequate short-and long-term measures to improve air quality. The use of air quality modelling tools at high resolution is now common but remains a challenge, because buildings create local modifications in air flows and modify surface energy budgets and roughness, leading to very complex processes which make air quality simulation particularly difficult in the near-ground layers [10][11][12][13]. Chemical transport models usually require numerical weather prediction (NWP) models, which provide meteorological input data [14,15]. High-resolution urban canopy parameters databases are becoming more and more available for cities [16,17]. NWP models with different physics parameterizations have been widely used in recent decades to simulate air quality over urban regions, especially the widely used community mesoscale Advanced Research Weather Research and Forecasting (ARW-WRF) model. Studies show that physics parameterizations play a critical role both in meteorological and air quality simulation or forecast [18][19][20][21][22]. For example, the Planetary Boundary Layer (PBL) height is determined by the heat and momentum exchanges between the PBL and the surface. In meteorological models, the PBL height is calculated by the PBL scheme. Therefore, the PBL scheme is a crucial factor in simulating the formation and evolution of air pollution [23][24][25]. Studies on different PBL parameterization schemes have shown that an accurate description of the meteorological conditions within the PBL via an appropriate parameterization scheme is important for air pollution modeling [26]. Some studies have presented the influences of meteorological models on pollutant concentrations; they estimate that uncertainties of PBL heights from the different models are one of the major sources of the differences in Particulate Matter (PM) concentrations [27,28]. Kleczeket et al. [29] compared three non-TKE (Turbulent Kinetic Energy) schemes and four TKE closures schemes; they found that non-TKE schemes tend to produce higher temperatures and wind speeds than in TKE schemes, especially during nighttime. Barlage et al. [13] also had a similar result between TKE and non-TKE schemes at night and found that PBL schemes have little effect on model performances in urban locations during daytime. Some other studies indicate that most PBL schemes slightly overestimate near-ground wind speed [30]; a possible reason could be the omission of the effect of unresolved topographic features on the momentum flux [31].
For the first time, the chemistry transport model CHIMERE model has been used with a very low first layer of about 4 m to simulate the air quality. In order to understand the impact of physics parameterizations and improve high-resolution air quality simulation, various WRF configurations have been tested over the Paris region during a critical pollution episode occurring from 27th November to 4th December 2016. In this study, we examine three different PBL schemes: (1) the Mellor-Yamada-Janjic scheme (MYJ), (2) the YonSei University scheme (YSU), and (3) the Bougeault and Lacarrere scheme (Boulac); and three canopy schemes, (1) a reference scheme (Bulk) which does not consider urban canopy parameters, (2) a single urban canopy model (UCM) which simply defines the urban geometry as two-dimensional street canyons but considers the three-dimensional nature of urban morphologies, and (3) a multilayer urban canopy model considering building effect Atmosphere 2020, 11,618 3 of 23 parameterization (BEP). BEP was developed by Martilli et al. [32], and allows the direct interaction between buildings and low atmospheric layers and therefore is more realistic in reproducing urban effects [16]. This work will inform the air quality community on our ability to assess and understand the atmospheric transport of air pollutants over cities for operational forecast with air quality models.

WRF Model Description
Several studies [33][34][35] demonstrate the importance of physics parameterizations in driving air quality modelling outputs. As the first step of this work, the WRF model (Version 3.9.1) is used to simulate the meteorological conditions. The WRF model provides several options for the parameterizations (details can be found in the WRF user manual) of physical processes, including multiple longwave and shortwave radiation schemes, surface layer schemes, land surface schemes, urban canopy schemes, and boundary layer schemes. The initial and boundary conditions are issued from the Global Forecasting System (GFS) analysis data from the National Centers for Environmental Prediction (NCEP), available at a 0.25 • × 0.25 • resolution at six-hourly time steps. The spectral nudging in WRF simulation is a way to constrain large-scale circulation and constrain the model to be more consistent with observations [36]. However, strong nudging may filter out extreme episodes, since nudging pushes the model toward a relatively smooth state [37]. Therefore, a weak spectral nudging above PBL has been tested to improve the accuracy of the downscaled fields but not filter out the simulation of extreme episodes.
The experiment set includes 11 simulations, which are conducted for a winter period from 27th November to 4th December 2016 and are designed to test the sensitivity of the surface meteorological variables and pollutant concentrations at these different physics parametrizations. Table 1 summarizes the parametrizations which are tested in the study. The UCM is a single urban canopy layer model used to consider the effects of urban geometry on the surface energy balance and wind shear for urban regions [38,39]. This model includes shadows from buildings; canyon orientation; the diurnal variation in the azimuth angle; the reflection of short and long wave radiation; wind profiles in the canopy layer; anthropogenic heating associated with energy consumption by human activities; and multi-layer heat transfer equations for roof, wall, and road surfaces. The BEP is the multilayer canopy layer scheme. This scheme considers the sub-grid wall, roof, and road effects on radiation and fluxes but does not include the energy exchange between the inside and outside of buildings.
The MYJ PBL scheme uses the 1.5-order turbulence closure model of Mellor and Yamada [40] to represent turbulence above the surface layer [41,42]. The MYJ scheme determines eddy viscosity where k is the von Karman constant set to 0.41, z is the altitude, h is the boundary layer height, and w s is the vertical scale given by the similarity formulae. The emissions of pollutants include several different gaseous and aerosol species. Those emissions can be split into four parts: anthropogenic emissions, biogenetic emissions, mineral dust emissions, and fire emissions. The anthropogenic emissions play a critical role in the CTM model, since they include all human activities, especially in megacities, and they are the only source we can reduce. For Europe region studies, anthropogenic emissions of gases and PM are generated with the Netherlands Organization for Applied Scientific Research (TNO) emission inventory [47], which is updated every six months. Biogenic emissions are estimated with the Model of Emissions and Gases and Aerosols from Nature (MEGAN). Mineral dust emissions are parameterized by Menut et al. [48]. The resolution TNO emission inventory is 0.125 • × 0.0625 • longitude-latitude or approximately 7 × 7 km, and the spatial resolutions of the emission data are recomputed from the TNO emission inventory for each domain. Driven by those input data, CHIMERE provides the hourly concentrations of tens of species at a regional scale; detailed descriptions of the model can be found in several previous studies [43,[47][48][49][50]. The aerosol module accounts for both inorganic and organic species, such as sulfate, nitrate, ammonium, primary organic matter, element carbon, secondary organic aerosol, sea salt, dust, and water. The aerosol size is assuming discrete aerosol size sections and that the particles are internally mixed. Nine diameter bins, from 10 nm to 10 µm, with 5 bins below 1 µm, are used [51].

Domains Setup and Observations Data
In order to calculate realistic meteorological input variables and pollutant concentrations, the system is configured with four nested-grid domains ( Figure 1) over Europe, France, the north of France, and Ile-de-France, with cell sizes of 45 (d01), 15 (d03), 5 (d02), and 1.67 km (d04), respectively. The finest grid domain (d04) used has 83 × 77 grid points and covers the whole Ile-de-France region, including urban, sub-urban, and rural areas. The WRF simulation is discretized vertically on 52 vertical levels from the surface to the top (50hPa), with 20 levels below 1000 m and the lowest level at 2 m, the complete name list of the WRF simulation can be found in List.S1. In order to take full advantage of the physical parameterizations, the CHIMERE simulation is performed with the same horizontal domains as the WRF; the 52 vertical levels from the WRF simulation are projected onto 33 levels in CHIMERE. The ground layer of the CHIMERE domains is at 999.5 hPa from about 4 m up to over 5000 m, with 15 levels below 1000 m. Such a low ground level could effectively reproduce the transport and diffusion of pollutants in a real urban environment.
Atmosphere 2020, 11, x FOR PEER REVIEW 5 of 25 to over 5000 m, with 15 levels below 1000 m. Such a low ground level could effectively reproduce the transport and diffusion of pollutants in a real urban environment. Figure 1. Four nested domains-from d01, the largest one, at a low resolution to d04, the finest one, at a high resolution.
Two types of observation data are used in the present study: (i) three meteorological observation stations operated by METEO FRANCE and SIRTA were used to evaluate the simulated surface meteorological variables from the model; and (ii) the hourly surface concentrations of criteria pollutants (NO2, PM2.5, and PM10) issued from AIRPARIF, the official Paris air quality network. The location of the three meteorological observation stations (Montsouris, Roissy, Sirta) and the four air quality monitoring stations can be found in Figure 2; they are urban or sub-urban background stations. The FR04143 air quality monitoring station is located in the city center. "Roissy" is a synoptic station of the Paris international Airport; "Montsouris" is an historical meteorological station located within Paris, operating since 1869; "Sirta" is a recent station located in a suburban environment.

Urban Parameters and Nudging Tests
The urban canopy model in the WRF simulation affects the modeling of several important meteorological variables for air quality issues (e.g., 10 m wind speed W10, 2 m surface temperature T2, and PBL height). It changes the vertical distribution of atmospheric pollutants by modifying the boundary layer height and mixing. Geometric and thermal parameters, including anthropogenic Two types of observation data are used in the present study: (i) three meteorological observation stations operated by METEO FRANCE and SIRTA were used to evaluate the simulated surface meteorological variables from the model; and (ii) the hourly surface concentrations of criteria pollutants (NO 2 , PM 2.5 , and PM 10 ) issued from AIRPARIF, the official Paris air quality network. The location of the three meteorological observation stations (Montsouris, Roissy, Sirta) and the four air quality monitoring stations can be found in Figure 2; they are urban or sub-urban background stations. The FR04143 air quality monitoring station is located in the city center. "Roissy" is a synoptic station of the Paris international Airport; "Montsouris" is an historical meteorological station located within Paris, operating since 1869; "Sirta" is a recent station located in a suburban environment.
Atmosphere 2020, 11, x FOR PEER REVIEW 5 of 25 to over 5000 m, with 15 levels below 1000 m. Such a low ground level could effectively reproduce the transport and diffusion of pollutants in a real urban environment. Two types of observation data are used in the present study: (i) three meteorological observation stations operated by METEO FRANCE and SIRTA were used to evaluate the simulated surface meteorological variables from the model; and (ii) the hourly surface concentrations of criteria pollutants (NO2, PM2.5, and PM10) issued from AIRPARIF, the official Paris air quality network. The location of the three meteorological observation stations (Montsouris, Roissy, Sirta) and the four air quality monitoring stations can be found in Figure 2; they are urban or sub-urban background stations. The FR04143 air quality monitoring station is located in the city center. "Roissy" is a synoptic station of the Paris international Airport; "Montsouris" is an historical meteorological station located within Paris, operating since 1869; "Sirta" is a recent station located in a suburban environment.

Urban Parameters and Nudging Tests
The urban canopy model in the WRF simulation affects the modeling of several important meteorological variables for air quality issues (e.g., 10 m wind speed W10, 2 m surface temperature T2, and PBL height). It changes the vertical distribution of atmospheric pollutants by modifying the boundary layer height and mixing. Geometric and thermal parameters, including anthropogenic

Urban Parameters and Nudging Tests
The urban canopy model in the WRF simulation affects the modeling of several important meteorological variables for air quality issues (e.g., 10 m wind speed W10, 2 m surface temperature T2, and PBL height). It changes the vertical distribution of atmospheric pollutants by modifying the boundary layer height and mixing. Geometric and thermal parameters, including anthropogenic Atmosphere 2020, 11, 618 6 of 23 heat, thermal conductivities, and heat capacity, have a significant effect on the transfer of momentum between the atmosphere and urban surface [52]. However, the representativeness of these parameters aggregated over a CTM grid cell remains highly uncertain, especially with a limited number of urban land-use types. The urban parameters (Table S1) used here are based on a previous study in Paris region [53,54] and default urban parameterization (URBPARM). The TBL tables in the WRF code are defined as URBAN-PAM1 and URBAN-PAM2, respectively. The RRTMG shortwave and longwave radiation scheme, MO similarity surface layer scheme, Noah land surface model, BEP urban canopy model, and MYJ boundary layer scheme are used for all simulations in this section. Figure 3 presents the time series of hourly-averaged T2 and W10 at Sirta station. The results show that simply modifying the values of the parameters has a strong impact on the outputs. The T2 presents a more positive correlation (0.86 than 0.81) in URBAN-PAM1 than in URBAN-PAM2, while the W10 has a more negative bias (−1.51 m than −1.42 m s −1 ) in URBAN-PAM1 than in URBAN-PAM2. A weak spectral nudging above the boundary layer for all the domains is also tested in this part. Figure 4 indicates that the spectral nudging technique is effective in improving the model skills. The performance benefits in Paris from the assimilation of meteorological data at Roissy airport in the global meteorological fields. Thus, URBAN-PAM1 and spectral nudging are used for all our sensitivity tests, which are discussed in the following sections.
Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 25 heat, thermal conductivities, and heat capacity, have a significant effect on the transfer of momentum between the atmosphere and urban surface [52]. However, the representativeness of these parameters aggregated over a CTM grid cell remains highly uncertain, especially with a limited number of urban land-use types. The urban parameters (Table S1) used here are based on a previous study in Paris region [53,54] and default urban parameterization (URBPARM). The TBL tables in the WRF code are defined as URBAN-PAM1 and URBAN-PAM2, respectively. The RRTMG shortwave and longwave radiation scheme, MO similarity surface layer scheme, Noah land surface model, BEP urban canopy model, and MYJ boundary layer scheme are used for all simulations in this section. Figure 3 presents the time series of hourly-averaged T2 and W10 at Sirta station. The results show that simply modifying the values of the parameters has a strong impact on the outputs. The T2 presents a more positive correlation (0.86 than 0.81) in URBAN-PAM1 than in URBAN-PAM2, while the W10 has a more negative bias (−1.51m than −1.42 m s −1 ) in URBAN-PAM1 than in URBAN-PAM2. A weak spectral nudging above the boundary layer for all the domains is also tested in this part. Figure 4 indicates that the spectral nudging technique is effective in improving the model skills. The performance benefits in Paris from the assimilation of meteorological data at Roissy airport in the global meteorological fields. Thus, URBAN-PAM1 and spectral nudging are used for all our sensitivity tests, which are discussed in the following sections.

Ground Meteorological Variables
The RRTMG shortwave and longwave radiation schemes, the MO similarity surface layer scheme, and the Noah land surface model are used for all the simulations in this section. The whole period (27th November to 4th December 2016) is divided into two periods: a Particulate matter Episode (PE) from 30th November to 2nd December; and other days than the PE, which are defined as regular days (RD). No rainfall occurred over both then PE and the RD. Figure 5 presents the maps of averaged T2 simulated based on the Bulk, UCM, and BEP schemes in this period. All the schemes show a higher T2 in the urban region than the rural, and the BEP scheme displays the strongest UHI effect in the urban region, which gives a positive bias of up to 1 °C compared to the other two schemes. The Bulk and BEP schemes perform similarly in rural areas and display a negative bias (1-3 °C by Bulk and 1-4 °C by BEP) compared with urban areas, but the UCM scheme gives the highest T2 in rural areas, displaying a negative bias of about 0.5-1 °C compared with urban regions. A possible reason could be the overestimation of surface heat flux (SHF) in the rural areas. Figure 6 presents the averaged SHF in this period. The Bulk and BEP schemes display obvious differences between urban and rural regions, with a negative bias (around −25 to −35 W m −2 ) in rural compared to urban areas. The UCM scheme just produces a negative bias of −10 and −15 W m −2 , respectively, in rural and urban areas, which could give an explanation of the overestimation of T2 in rural areas.

Ground Meteorological Variables
The RRTMG shortwave and longwave radiation schemes, the MO similarity surface layer scheme, and the Noah land surface model are used for all the simulations in this section. The whole period (27th November to 4th December 2016) is divided into two periods: a Particulate matter Episode (PE) from 30th November to 2nd December; and other days than the PE, which are defined as regular days (RD). No rainfall occurred over both then PE and the RD. Figure 5 presents the maps of averaged T2 simulated based on the Bulk, UCM, and BEP schemes in this period. All the schemes show a higher T2 in the urban region than the rural, and the BEP scheme displays the strongest UHI effect in the urban region, which gives a positive bias of up to 1 • C compared to the other two schemes. The Bulk and BEP schemes perform similarly in rural areas and display a negative bias (1-3 • C by Bulk and 1-4 • C by BEP) compared with urban areas, but the UCM scheme gives the highest T2 in rural areas, displaying a negative bias of about 0.5-1 • C compared with urban regions. A possible reason could be the overestimation of surface heat flux (SHF) in the rural areas. Figure 6 presents the averaged SHF in this period. The Bulk and BEP schemes display obvious differences between urban and rural regions, with a negative bias (around −25 to −35 W m −2 ) in rural compared to urban areas. The UCM scheme just produces a negative bias of −10 and −15 W m −2 , respectively, in rural and urban areas, which could give an explanation of the overestimation of T2 in rural areas.

Ground Meteorological Variables
The RRTMG shortwave and longwave radiation schemes, the MO similarity surface layer scheme, and the Noah land surface model are used for all the simulations in this section. The whole period (27th November to 4th December 2016) is divided into two periods: a Particulate matter Episode (PE) from 30th November to 2nd December; and other days than the PE, which are defined as regular days (RD). No rainfall occurred over both then PE and the RD. Figure 5 presents the maps of averaged T2 simulated based on the Bulk, UCM, and BEP schemes in this period. All the schemes show a higher T2 in the urban region than the rural, and the BEP scheme displays the strongest UHI effect in the urban region, which gives a positive bias of up to 1 °C compared to the other two schemes. The Bulk and BEP schemes perform similarly in rural areas and display a negative bias (1-3 °C by Bulk and 1-4 °C by BEP) compared with urban areas, but the UCM scheme gives the highest T2 in rural areas, displaying a negative bias of about 0.5-1 °C compared with urban regions. A possible reason could be the overestimation of surface heat flux (SHF) in the rural areas. Figure 6 presents the averaged SHF in this period. The Bulk and BEP schemes display obvious differences between urban and rural regions, with a negative bias (around −25 to −35 W m −2 ) in rural compared to urban areas. The UCM scheme just produces a negative bias of −10 and −15 W m −2 , respectively, in rural and urban areas, which could give an explanation of the overestimation of T2 in rural areas.  Results at the three meteorological stations are compared with the observations to evaluate the canopy scheme's performance. Figure 7 presents the time series of hourly-averaged T2 during the whole period. The observed T2 varies between −1.4 and 9.3 °C, −3.9 and 7.8 °C, and −2.1 to 8.4 °C, with an average of 4.0, 2.6, and 4.9 °C at Montsouris, Roissy, and Sirta stations, respectively. Compared with the observations, all three schemes tended to underestimate the T2 during the period 16:00 pm-10:00 am UTC on the 27th and 28th for all stations, while higher positive values were observed around midnight of the 1st December for Montsouris and Roissy stations, especially at Roissy, with a positive bias from 1.8 to 4.9 °C during the PE. The BEP scheme simulates relatively higher T2 at PE for all the stations because the BEP scheme simulates a stronger UHI than the other two schemes. All the schemes have the best performance at the Sirta site. Compared with the observations from three stations, the differences in the modelled T2 are minimal between three stations during the PE. However, there are significant differences in the observed T2 between the stations. For example, the minimum observed T2 are −1.4, −3.9, and −0.3 °C during the midnight of 1st December at Montsouris, Roissy, and Sirta stations, respectively. A possible reason could be the use of a single urban category in our WRF simulations, but in the real environment Montsouris and Sirta stations are closer to commercial and residential areas and Roissy is a synoptic station close to the airport, which could be assimilated as a rural environment with a weaker UHI effect, especially during nighttime. Results at the three meteorological stations are compared with the observations to evaluate the canopy scheme's performance. Figure 7 presents the time series of hourly-averaged T2 during the whole period. The observed T2 varies between −1.4 and 9.3 • C, −3.9 and 7.8 • C, and −2.1 to 8.4 • C, with an average of 4.0, 2.6, and 4.9 • C at Montsouris, Roissy, and Sirta stations, respectively. Compared with the observations, all three schemes tended to underestimate the T2 during the period 16:00 p.m.-10:00 a.m. UTC on the 27th and 28th for all stations, while higher positive values were observed around midnight of the 1st December for Montsouris and Roissy stations, especially at Roissy, with a positive bias from 1.8 to 4.9 • C during the PE. The BEP scheme simulates relatively higher T2 at PE for all the stations because the BEP scheme simulates a stronger UHI than the other two schemes. All the schemes have the best performance at the Sirta site. Compared with the observations from three stations, the differences in the modelled T2 are minimal between three stations during the PE. However, there are significant differences in the observed T2 between the stations. For example, the minimum observed T2 are −1.4, −3.9, and −0.3 • C during the midnight of 1st December at Montsouris, Roissy, and Sirta stations, respectively. A possible reason could be the use of a single urban category in our WRF simulations, but in the real environment Montsouris and Sirta stations are closer to commercial and residential areas and Roissy is a synoptic station close to the airport, which could be assimilated as a rural environment with a weaker UHI effect, especially during nighttime.  Figure 8 shows the maps of the averaged W10 simulated for the Bulk, UCM, and BEP schemes in this period. All the schemes give a higher surface wind speed in rural areas than in urban regions; the Bulk scheme shows the weakest gradient from urban to rural areas, with a negative bias for 0.51 m s −1 ; UCM gives a relative higher negative bias than Bulk from 1 to 2 m s −1 ; and the BEP scheme had an  Figure 8 shows the maps of the averaged W10 simulated for the Bulk, UCM, and BEP schemes in this period. All the schemes give a higher surface wind speed in rural areas than in urban regions; the Bulk scheme shows the weakest gradient from urban to rural areas, with a negative bias for 0.51 m s −1 ; UCM gives a relative higher negative bias than Bulk from 1 to 2 m s −1 ; and the BEP scheme had an obviously strongest impact from cities, displaying a negative bias around 3.5 m s −1 . The comparison between the observations and models for W10 are displayed in Figure 9. The observed W10 varied between 0.6 m and 7.6 m s −1 , 0.1 m and 9.4 m s −1 , and 0.2 m and 5.5 m s −1 , with an average of 2.3 m, 1.5 m, and 2.8 m s −1 in three meteorological observation stations (Montsouris, Roissy, Sirta), respectively. Both the observed and modelled W10 display much lower values during the PE than the RD. In contrast with the T2, the W10 does not display remarkable differences between stations. The BEP scheme tends to produce the lowest W10 throughout the whole period for all the stations. Bulk and UCM performed well for both the PE and RD. A possible reason could be due to the roughness length, which is not directly dependent on urban morphology for the UCM scheme, while the momentum sinks and drag force estimated by the BEP scheme depend on urban morphology in this study. Low wind speed leads to weak natural circulation and the accumulation of heat in the lower urban boundary layer, which in turn enhances the warming of the urban area [55]. This would explain why BEP has a relative higher T2 simulation than the other schemes during the PE. This sensitivity test indicates that more realistic urban morphologies are needed to improve the model performance, because the parameters used in the urban canopy scheme (e.g., albedo, heat capacity, and roughness length) strongly influence the surface energy balance and heat absorption.

Air Quality Modeling
In this study, NO2, PM2.5, and PM10 concentrations are used to evaluate the impact of physics parameterizations on air quality simulations. NO2, which comes mainly from automobile exhaust

Air Quality Modeling
In this study, NO 2 , PM 2.5 , and PM 10 concentrations are used to evaluate the impact of physics parameterizations on air quality simulations. NO 2 , which comes mainly from automobile exhaust emissions in cities, is selected to be an indicator of local pollution. Comparing the evolution and spatial patterns of NO 2 versus the PM 2.5 and PM 10 gives an indication of the local contribution during this episode. Figure 10 presents the differences in pollutants distribution between each canopy scheme. The results show that, compared to the Bulk scheme, the pollutant concentrations from the BEP scheme were lower in the city center, but displayed an opposite behavior compared to the UCM scheme. The UCM scheme displayed higher pollutant concentrations in the urban area and lower pollutant concentrations in rural areas. Figure 11 presents the time series of NO 2 mass concentrations for all the experiments compared to the observations for the four stations (the time series of PM 2.5 and PM 10 mass concentrations can be found in Figure S1). The average biases of the four stations is given in Table 2, and the time series of the biases are displayed in Figure S2. The results show that the BEP scheme gives the best performance for almost the whole period, especially during the PE. All the schemes can effectively predict air pollution during the RD, but they usually overestimate pollutant concentrations for all the stations during the PE; all the schemes had their worst performance in NO 2 prediction, which means that the canopy schemes do not simulate well the emission of local air pollution. Figure 12 gives the vertical diffusion coefficient at daytime and nighttime in the RD and PE. The Bulk scheme has the highest vertical diffusion rate among the three schemes during the RD. Figure 13 presents the vertical profile of the NO 2 and PM 2.5 mass concentrations under 200 m; both NO 2 and PM 2.5 show the lowest vertical mass concentrations at low layers by the Bulk scheme and the highest vertical mass concentrations by the BEP scheme during the RD, which agrees with the vertical profile of Kz in CHIMERE. At 06:00 UTC on 30th November, the Kz displays very low values at ground layers for all three schemes, but rises rapidly in the BEP scheme from the third layer. Figure 13c also presents high pollutant concentrations at the ground layer for all three schemes but decreases rapidly in the BEP scheme from the third layer. Looking at the NO 2 concentrations, there is clearly an added value of using a low first model layer in the model to enhance the concentrations close to emission sources related to ground sources such as road traffic. In the models, it is more common to use a first layer of about 20 to 50 m [56,57]. Figure S3 presents the hourly timeseries of NO 2 as an average of all the air quality monitoring stations between the observations and simulations, with first layers of about 4, 18, and 40 m (based on the best physics parameterizations in this study). An extremely low vertical diffusion rate leads to high concentrations of pollutants, which is one reason for the strong overestimations at the surface layers at this time. In general, the BEP scheme has the best performance among the three schemes, but all the schemes strongly overestimate the concentrations at the beginning of the PE, especially the UCM scheme. A low near-ground wind speed usually indicates low diffusion rates, which means the BEP scheme involves higher pollutant concentrations than for the other two schemes, but our study shows the pollutants concentrations is not found to be positively correlated with the W10. A possible explanation could be the underestimation of PBL height, which leads to a low diffusion rate at this period; a detail discussion about the impact of the PBL height is in the next section.
concentrations at the beginning of the PE, especially the UCM scheme. A low near-ground wind speed usually indicates low diffusion rates, which means the BEP scheme involves higher pollutant concentrations than for the other two schemes, but our study shows the pollutants concentrations is not found to be positively correlated with the W10. A possible explanation could be the underestimation of PBL height, which leads to a low diffusion rate at this period; a detail discussion about the impact of the PBL height is in the next section.

PBL Height
WRF version 3.9.1 provides 14 options for the PBL scheme, following the approach of LeMoine et al. [58,59]. Briefly, they can be categorized as TKE schemes and non-TKE schemes. Two TKE-based schemes and one non-TKE-based PBL scheme are selected for these PBL sensitivity simulations. The TKE schemes chosen in this assessment are the MYJ and Boulac schemes, and the non-TKE scheme is the YSU scheme. Unified PBL sensitivities are not possible in this section with the canopy schemes from the previous section, because the BEP urban canopy scheme requires the simultaneous use of a TKE scheme and the Boulac PBL scheme is designed for coupling with the multi-urban canopy model. Regarding the configuration of the canopy models in Section 3.2, each canopy scheme will be coupled with two PBL schemes, giving the YSU-Bulk and MYJ-Bulk schemes, the YSU-UCM and MYJ-UCM schemes, and the Boulac-BEP and MYJ-BEP schemes. All the other physics options are the same. For the observed PBL height, we have the complete dataset for Roissy station and the last four days data for Sirta station. Figure 14 displays a comparison between the observed Lidar data and the 15 min averaged modeled output of PBL height. Generally, all the PBL schemes systematically underestimate PBL height for the whole time series, especially during the PE and the beginning of RD (from 27th November 00:00 UTC). This explains the strong overestimations of pollutants during the PE in Section 3.2 and is consistent with the positive biases of pollutant concentrations at the beginning of the RD in Figure 10. If we consider the simulation of PBL during the PE, the average observed PBL height was around 110 m and is only relatively well reproduced by the Boulac-BEP, scheme even though it still underestimates the PBL height by around 40 m on average. The other PBL schemes reach the minimum PBL height (set up to 20 m in the model) during the whole PE and strongly underestimate PBL height, especially at night and early in the morning. In general, the Boulac-BEP scheme gives the best performance, and it is the only option which can effectively estimate the evolution of the PBL height during the PE in this sensitivity test. For the other PBL schemes, the agreement between the simulations and observations is poor. MYJ displays the poorest correlation compared to the observed soundings, regardless of other parameter settings. represent Roissy and Sirta station, respectively; 1, 2, and 3 represent each PBL scheme coupled with each urban canopy scheme.

Ground Meteorological Variables
A comparison of the performances of various PBL schemes in simulating the surface meteorological variables (T2 and W10) is presented in this section. Figures 15 and 16 shows the time series of T2 and W10 based on several PBL schemes. The YSU-Bulk and Boulac-BEP have a slight positive bias compared to MYJ, but the YSU-UCM scheme had a cooling effect compared to the MYJ-UCM scheme during the PE. In the MYJ-BEP and Boulac-BEP schemes, it can be observed that the models always underestimate the wind speed compared to the observations in this period. The Boulac-BEP scheme produces a lower W10 than MYJ-BEP throughout the whole period, especially during RD. In contrast, W10 does not show significant differences between the MYJ and YSU schemes coupled with the other two Bulk and UCM canopy models. Table 3 summarizes the performances of the models, including the mean bias (MB) and mean Pearson correlation coefficient (R) of meteorological variables. These schemes display a high correlation which approximately reproduces the measured trend. The three experiments with the MYJ boundary layer scheme display a negative bias in W10 in Roissy and a negative bias in T2 in both Roissy and Sirta stations. A less negative bias is noticed with the simulations using the Boulac-BEP scheme compared with the others for the whole period, but it still produces a positive bias during nighttime in the PE. The reason could be that the Boulac PBL scheme takes the turbulent exchange affected by complex terrain into account in the represent Roissy and Sirta station, respectively; 1, 2, and 3 represent each PBL scheme coupled with each urban canopy scheme.

Ground Meteorological Variables
A comparison of the performances of various PBL schemes in simulating the surface meteorological variables (T2 and W10) is presented in this section. Figures 15 and 16 shows the time series of T2 and W10 based on several PBL schemes. The YSU-Bulk and Boulac-BEP have a slight positive bias compared to MYJ, but the YSU-UCM scheme had a cooling effect compared to the MYJ-UCM scheme during the PE. In the MYJ-BEP and Boulac-BEP schemes, it can be observed that the models always underestimate the wind speed compared to the observations in this period. The Boulac-BEP scheme produces a lower W10 than MYJ-BEP throughout the whole period, especially during RD. In contrast, W10 does not show significant differences between the MYJ and YSU schemes coupled with the other two Bulk and UCM canopy models. Table 3 summarizes the performances of the models, including the mean bias (MB) and mean Pearson correlation coefficient (R) of meteorological variables. These schemes display a high correlation which approximately reproduces the measured trend. The three experiments with the MYJ boundary layer scheme display a negative bias in W10 in Roissy and a negative bias in T2 in both Roissy and Sirta stations. A less negative bias is noticed with the simulations using the Boulac-BEP scheme compared with the others for the whole period, but it still produces a positive bias during nighttime in the PE. The reason could be that the Boulac PBL scheme takes the turbulent exchange affected by complex terrain into account in the calculation of the turbulent diffusion coefficient, leading to a smaller bias in T2 prediction [60]. The Boulac-BEP scheme slightly improves the results with respect to the MYJ-BEP scheme for the W10, but still massively underestimates by 1.53 m s −1 at Roissy station. The correlation indicates the PBL schemes coupled with the urban canopy schemes (MYJ-UCM, YSU-UCM, MYJ-BEP, and Boulac-BEP) are more accurate than those with no urban canopy scheme (MYJ-Bulk and YSU-Bulk). In contrast to the PBL height results, the model performances for T2 and W10 based on various PBL schemes are mixed. Therefore, it is not straightforward to agree on a definitive best choice for the PBL scheme.     Figure 17; the PM2.5 and PM10 can be found in Figures S4 and  S5). In general, all the schemes have a better performance for the RD than the PE period, and almost all the schemes strongly overestimate pollutant concentrations at the beginning of the PE, except the Boulac-BEP scheme. The bias between the modelled and observed pollutant concentrations can be The impacts of the PBL schemes on air quality are analyzed in this section. The time series of the model predictions of the NO 2 , PM 2.5 , PM 10 concentrations for all the schemes are compared with the observations (the NO 2 can be found in Figure 17; the PM 2.5 and PM 10 can be found in Figures S4 and  S5). In general, all the schemes have a better performance for the RD than the PE period, and almost all the schemes strongly overestimate pollutant concentrations at the beginning of the PE, except the Boulac-BEP scheme. The bias between the modelled and observed pollutant concentrations can be found in Figures S6-S8. All the schemes give a worse performance for NO 2 and PM 2.5 prediction than for PM 10 . The PBL height reaches the minimum value at 18:00 UTC on 29th November, with the YSU-Bulk scheme corresponding to an underestimation of 500 m; such a difference creates a large bias in air pollutant concentrations. The same phenomenon occurs for the MYJ-Bulk, MYJ-UCM, YSU-UCM, and MYJ-BEP schemes at 00:00 UTC on 30th November. However, the overestimation of pollutant concentrations does not occur at this time but at 06:00 UTC on 1st December, because this episode is largely influenced by local sources (road traffic and residential heating). On 1st December 2016 at 06:00 UTC (a weekday), the morning rush hour leads to a lag effect in the bias, explaining why all the schemes have their worst performance for the NO 2 and PM 2.5 predictions. The averaged statistics for the four air quality monitoring stations are summarized in Table 4. Figure 18 gives the scatterplot between the observed and modelled NO 2 ; the scatterplot for the PM 2.5 and PM 10 can be found in Figure S9. Overall, the MYJ scheme overestimates the mass concentrations of all pollutants, whatever the selected urban canopy model. The Boulac-BEP scheme has a significantly better performance than other schemes for all the evaluated variables. The YSU-Bulk scheme has the lowest mean bias (−1.3 µg m −3 ) in PM 10 prediction, but its mean linear correlation coefficients of NO 2 , PM 2.5 , and PM 10 were 0.52, 0.40, and 0.53, respectively. The mean biases for the NO 2 and PM 2.5 prediction were 17.6 and 46.0 µg m −3 , respectively. The low correlation for all the analyzed variables, positive biases for NO 2 and PM 2.5 prediction, and negative bias for PM 10 prediction indicate that this scheme under and over-estimates the concentrations simultaneously, which means that the YSU-Bulk scheme cannot be performed effectively for the prediction of short-term PE in urban regions. The Boulac-BEP scheme was better for all the analyzed variables, both in terms of bias and correlation statistics; its mean biases for NO 2 , PM 2.5 , and PM 10 prediction were −5.1, 1.2, and −8.6 µg m −3 , respectively. This indicates that both the urban canopy and PBL schemes play an important role in the PE forecast for urban regions.

Conclusions
This study evaluates the impact of three urban canopy parametrization schemes (Bulk, UCM, and BEP) and three boundary layer parametrization schemes coupled with urban canopy schemes

Conclusions
This study evaluates the impact of three urban canopy parametrization schemes (Bulk, UCM, and BEP) and three boundary layer parametrization schemes coupled with urban canopy schemes (MYJ-Bulk, YSU-Bulk, MYJ-UCM, YSU-UCM, MYJ-BEP, and Boulac-BEP) from the Weather Research and Forecasting (WRF Version 3. 9.1) mesoscale model on the simulations of meteorological variables and predicted pollutant concentrations issued from the CHIMERE chemical transport model. This study aims at understanding those impacts on the short-term pollution episode in the Paris region, France.
For the first time, chemistry transport models have been used with a first layer at 4 m to better account for physical processes within the suburban layer in relation to ground emissions such as road traffic. Clearly, the sharp gradient of concentrations shows the added value of using such a low level

Conclusions
This study evaluates the impact of three urban canopy parametrization schemes (Bulk, UCM, and BEP) and three boundary layer parametrization schemes coupled with urban canopy schemes (MYJ-Bulk, YSU-Bulk, MYJ-UCM, YSU-UCM, MYJ-BEP, and Boulac-BEP) from the Weather Research and Forecasting (WRF Version 3. 9.1) mesoscale model on the simulations of meteorological variables and predicted pollutant concentrations issued from the CHIMERE chemical transport model. This study aims at understanding those impacts on the short-term pollution episode in the Paris region, France.
For the first time, chemistry transport models have been used with a first layer at 4 m to better account for physical processes within the suburban layer in relation to ground emissions such as road traffic. Clearly, the sharp gradient of concentrations shows the added value of using such a low level in urban areas. This study shows the high sensitivity of air quality simulation to meteorological model outputs during stable conditions favorable to the occurrence of PM episodes.
Three urban canopy models could reproduce the T2 effectively in the daytime; the Bulk and UCM model reproduced the W10 well, but the BEP model underestimates W10 for the whole episode, indicating an important impact of buildings on the diffusion rate for this scheme. The UCM scheme overestimates the surface heat index in rural areas, leading to a slight negative bias for T2 in rural versus urban areas. T2 displays ambiguous results, while all schemes reproduced results well at Sirta station but overestimated during nighttime at Montsouris and Roissy during a pollution episode. This could be due to the use of a single urban category in this study, which leads to a discrepancy that causes more spatial variability in urban parameters. NO 2 , PM 2.5 , and PM 10 were selected to evaluate the effect of physics parameterizations on air quality simulation. Large differences were found in the simulation of pollutant concentrations, with an overall overestimation during the pollution episode.
The six PBL sensitivity tests display an underestimation of the PBL height when compared with the observations issued from the Lidar. All the sensitivity tests, except Boulac-BEP, could not effectively simulate the PBL height during the pollution episode. The two PBL schemes coupled with the BEP scheme still underestimated the W10 for the whole period, but the Boulac-BEP scheme improved the simulations compared with the MYJ-BEP scheme. The Boulac-BEP scheme had a significantly better performance than the other schemes for the simulation of pollutant concentrations, indicating that both the canopy parametrization and PBL schemes have a positive effect on air quality simulation in urban regions, especially during pollution episodes. An advanced urban category could improve air quality prediction in urban regions.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/6/618/s1: Table S1. The urban parameters for the urban canopy model; Figure S1. Time series of PM 2.5 , PM 10 Figure S4. Time series between observed and modelled PM 2.5 : a, b and c represent each PBL schemes coupled with each urban canopy schemes respectively; 1, 2, 3, and 4 represent Station FR04002, FR04034, FR04143, FR04156, respectively; Figure S5. Time series between observed and modelled PM 10 : a, b and c represent each PBL schemes coupled with each urban canopy schemes respectively; 1, 2, 3, and 4 represent Station FR04002, FR04034, FR04143, FR04156, respectively; Figure S6. Time series of bias between observed and modelled NO 2 : a, b and c represent each PBL schemes coupled with each urban canopy schemes respectively; 1, 2, 3, and 4 represent Station FR04002, FR04034, FR04143, FR04156, respectively; Figure S7. Time series of bias between observed and modelled PM 2.5 : a, b and c represent each PBL schemes coupled with each urban canopy schemes respectively; 1, 2, 3, and 4 represent Station FR04002, FR04034, FR04143, FR04156, respectively; Figure S8. Time series of bias between observed and modelled PM 10 : a, b and c represent each PBL schemes coupled with each urban canopy schemes respectively; 1, 2, 3, and 4 represent Station FR04002, FR04034, FR04143, FR04156, respectively; Figure S9. Liner fit between observed and modelled PM 2.5 and PM 10 ; List.S1. Name list of WRF model.