High Spatial Resolution Assessment of the Effect of the Spanish National Air Pollution Control Programme on Street-Level NO 2 Concentrations in Three Neighborhoods of Madrid (Spain) Using Mesoscale and CFD Modelling

: Current European legislation aims to reduce the air pollutants emitted by European countries in the coming years. In this context, this article studies the effects on air quality of the measures considered for 2030 in the Spanish National Air Pollution Control Programme (NAPCP). Three different emission scenarios are investigated: a scenario with the emissions in 2016 and two other scenarios, one with existing measures in the current legislation (WEM2030) and another one considering the additional measures of NAPCP (WAM2030). Previous studies have addressed this issue at a national level, but this study assesses the impact at the street scale in three neighborhoods in Madrid, Spain. NO 2 concentrations are modelled at high spatial resolution by means of a methodology based on Computational Fluid Dynamic (CFD) simulations driven by mesoscale meteorological and air quality modelling. Spatial averages of annual mean NO 2 concentrations are only estimated to be below 40 µg/m 3 in all three neighborhoods for the WAM2030 emission scenarios. However, for two of the three neighborhoods, there are still zones (4–12% of the study areas) where the annual concentration is higher than 40 µg/m 3 . This highlights the importance of considering microscale simulations to assess the impacts of emission reduction measures on urban air quality.

under the Directive on Ambient Air Quality and Cleaner Air for Europe (Air Quality Directive, 2008/50/EC). In this context, Spain has developed a national programme (National Air Pollution Control Programme, NAPCP) which includes a set of existing (WEM scenario) and additional measures (WAM scenarios) to further reduce pollutant emissions and meet the target emissions proposed in the NEC directive. The WEM emission scenario assumes that only existing and planned policy measures are implemented, whereas the WAM scenario includes additional measures that can be applied. The WAM scenario of the Spanish NAPCP considers 50 measures split between eight sectors, including electricity generation, road and non-road transport, improvements in energy efficiency, and agriculture. The NAPCP provides emission estimates for the years 2020, 2025, and 2030. For the year 2030, the national reductions in NOx emissions are estimated to be 7% and 33% for the WEM and WAM scenario, respectively (with reductions of 7% and 48% for road transport emissions). More details of the measures can be found in [2].
Citizens living in large metropolitan areas are exposed to elevated concentrations of NO2 mainly because of traffic emissions [1]. However, estimating the exposure of urban population is a challenge due to the high spatial variability of air pollutant concentrations [3,4]. The spatial distribution of concentrations is very heterogeneous as a consequence of the complex air flows in the street due to the interaction between the atmosphere and urban obstacles. Strong gradients of pollutant concentrations have been found in field experiments [5][6][7] and in model studies [8][9][10][11][12][13][14][15][16][17]. This also means that the spatial representativeness of urban air quality monitoring stations (AQMS) is very limited [18][19][20]. Highresolution modelling is, therefore, necessary to capture the spatial variability of air pollutant concentrations and estimating population exposure appropriately. Computational fluid dynamic (CFD) models are adequate tools for this because they have a spatial resolution high enough to capture such variability and, in fact, have been successfully applied to simulate pollutant dispersion in real urban environments [11,14,16,[21][22][23][24][25]. The main limitation of these models is the high computational cost required that limits the size of the domain, and, in most of cases, makes it unfeasible to carry out unsteady simulations of long time periods (e.g., months, a year). This issue is important in a policy context because annual statistics are relevant for air quality directives. However, post-processing methodologies have been developed and successfully applied to compute high-resolution maps of long-period (e.g., yearly) average concentrations [15,17,22,[26][27][28][29][30].
In this context, the following research question is posed: How does the Spanish emission reduction programme affect yearly-averaged concentrations at the street level?
The impact of the Spanish National Air Pollution Control Programme (NAPCP) on air quality and health impacts has been previously assessed at a national level at a coarse resolution (10 × 10 km 2 ) [2,31]. However, the distribution of concentrations of pollutants such as NO2 in urban areas at microscale get smoothed out in coarse resolution treatments. The main objective of the present study is to assess the effectiveness of the Spanish NAPCP on reducing NO2 concentrations at street level in different urban areas at highspatial resolution. The main novelty for achieving this objective is the use of a novel methodology that combines CFD simulations with mesoscale simulations to obtain high spatial resolution distributions of annual average NO2 concentrations at pedestrian level in three different neighborhoods of Madrid, Spain. The main hypothesis is that microscale modelling is needed to assess the effects of the Spanish National Air Pollution Control Programme (NAPCP) on annual mean NO2 concentrations in urban areas because pollutant concentrations get smoothed out at coarse resolutions. Therefore, despite the annual mean limit value not being exceeded in terms of the spatially-averaged concentration (coarse resolution), the limit value can be exceeded in some areas at microscale within a neighborhood. This paper focuses on NO2 since the European limit value for NO2 was exceeded at several AQMSs in Madrid-including the three stations studied here-and this study aims to investigate the effectiveness of emission reduction scenarios on these NO2 exceedances. Three emission scenarios considering 2016 emissions, emissions with existing measures in the current legislation in 2030 (WEM2030), and emissions with the additional measures of the NAPCP in 2030 (WAM2030) are simulated, and the comparison between WEM2030 and 2016 results and between WAM2030 and 2016 results used to quantify the effectiveness of the measures in terms of NO2 concentration at street level.

Description of the Study Urban Areas
Three neighborhoods of Madrid, Spain, are selected around three of the Air Quality Monitoring Stations (AQMSs) that usually record the highest concentration levels within the city. These AQMSs are: Plaza Elíptica (PE), Escuelas Aguirre (EA), and Plaza del Carmen (PC) (Figure 1a) and belong to the AQMS network of the Madrid city council. The annual average of NO2 concentration recorded at these AQMSs were above 40 µg m −3 , which is the annual average limit value for EU legislation (56 µg m −3 for PE, 57 µg m −3 for EA and 46 µg m −3 for PC) [32].
The first study area is located in the South-West of Madrid city around the trafficoriented AQMS "Plaza Elíptica" (PE). It is a highly polluted microenvironment with a complex urban geometry. The center of the area is a heavily trafficked roundabout with a freeway passing under it through a tunnel ( Figure 1b). In addition, many pedestrians transit this area. More details can be found in [22].
The second study neighborhood is located in the central area of Madrid city around the "Escuelas Aguirre" (EA) traffic AQMS (Figure 1c). This area is characterized by a large green urban area (El Retiro park) and avenues and streets with intense road traffic. More details can be found in [15].
Finally, the third study area is located in the center of Madrid around the AQMS "Plaza del Carmen" (PC) representing an urban background AQMS (Figure 1d). The study area is characterized by a wide pedestrian zone, though this station is situated close to Gran Vía Avenue, which had intense traffic in 2016 (year used in this work) and recorded NO2 concentrations are influenced by these traffic emissions [25].

NO2 Concentration Modelling
Unsteady CFD simulations of the time evolution of NO2 over one year are not feasible due to the huge computational resource required. To address the issue, annual average NOx concentrations are computed using a numerical methodology (WA CFD-RANS) based on a set of steady simulations for 16 wind directions. Similar methodologies were successfully applied for the EA and PE study areas by [15,22], respectively. NO2 concentrations are computed from NOx concentrations using the NO2/NOx ratios for each hour. This approach was successfully applied by [17] using the NOx and NO2 concentrations recorded at the AQMS within the study area. In the present study, the NO2/NOx ratios used are obtained from the concentration at ground level computed by mesoscale simulations in the mesoscale cells where the microscale domains are located. CFD model and set-up and the methodology used to estimate annual average concentrations in this study are described in Sections 2.2.1 and 2.2.2, respectively.

CFD Model Description and Set-Up
The CFD simulations are based on Reynolds-averaged Navier-Stokes (RANS) equations with a k-ε turbulence closure, where k is turbulent kinetic energy and ε is its dissipation rate. The aerodynamic effects of the vegetation in the study areas are modelled by means a sink of momentum and sink and source terms in the k and ε equations in the numerical cells with vegetation [24,33,34]. Other sources of turbulence such as turbulence induced by traffic are neglected, meaning that the diffusion of pollutants is underestimated. Pollutant diffusivity is related to turbulent Schmidt number (Sc). The optimum value of Sc is between 0.2 and 1.3, depending on flow properties and geometry [35]. The decrease or increase of Sc allows increasing or decreasing the diffusivity, respectively. In these simulations, in order to minimize the error of the underestimation of the diffusion of pollutants, a low value of Schmidt number (0.3), which increases the turbulent diffusion, is selected [22].
The real urban morphologies of the three neighborhoods are simulated. For the PE area (Figure 1e), the dimensions of the numerical domain are 1300 m × 1300 and it is discretized by an irregular polyhedral mesh of 8.3 × 10 6 grid points, including hexahedral prism layer close to buildings. The spatial resolution is finer than 1 m close to buildings and emissions and the cell size progressively increases in all directions with an expansion ratio lower than 1.3 up to a resolution of 5 m outside of the study area.
The dimensions of the numerical domain of the EA area are 700 m × 800 m and it is discretized by an irregular polyhedral mesh of 3 × 10 6 numerical cells, including hexahedral prism layer close to the buildings. The spatial resolution is about 1 m close to the buildings and each street and building are resolved with at least 10 grid points.
The PC study zone is 1 km × 1 km and is discretized by an irregular polyhedral mesh with 9.3 × 10 6 grid points. The study zone is discretized with a resolution of 2 m with refinement of around 1 m close to buildings and ground.
The three domains were built taking into account the best practice guideline of COST Action 732 [36] (e.g., the heights of simulation domains are >6 times the height of the tallest buildings). The validity of the meshes is verified by means of tests of the grid-independence of the results made by [22] for the mesh used for the PE domain, by [15] for the mesh used for the EA domain, and by [25] for the mesh used for the PC domain.
Sixteen wind directions (N, NNE, and so on, every 22.5 degrees) are simulated and the lateral boundaries are defined as velocity inlet or pressure outlet depending on the forcing wind direction. For each simulation, vertical profiles of wind speed (u), turbulent kinetic energy (k), and its dissipation rate (ε) for neutral atmospheric conditions are imposed [37] (Equations (1)-(3)): where u* is the friction velocity, κ is von Karman's constant (0.4), z0 is the roughness length, and Cµ is a model constant (0.09). These profiles are widely used in CFD simulation over urban environments [10,15,23]. A transport equation for non-reactive pollutants is used to simulate the NOx dispersion. NOx (and NO2) concentrations are set to 0 at boundaries and the background concentrations are subsequently added in the methodology for estimating annual average concentration (Section 2.2.2). Only NOx traffic emissions inside the numerical domain are considered in CFD simulations and contributions of other sources outside of the domain are included through the background concentration provided by mesoscale simulations. Traffic emissions are located in the roads and distributed taking into account the mean daily traffic for each street. For the PE simulations, traffic emissions used are distributed taking into account the mean daily traffic of each street based on the emissions used in [38]. For the EA and PC simulations, the PE temporal evolution of traffic emissions is used and the traffic emission factors are modified in order to provide the same annual average NOx concentration as that recorded at the AQMS location of each area [15]. For WEM2030 and WAM2030 scenarios, NOx traffic emissions in CFD simulations are reduced considering the reductions estimated for the road transport sector (reduction of 7% and 48% for WEM2030 and WAM2030, respectively). Note that one of the measures of WAM2030 is the promotion of electric vehicles.

Methodology for Estimating Annual Average Concentration
The methodology described above can be only applied for non-reactive pollutants. For this reason, NOx is modelled, which is considered as a non-reactive pollutant. NOx is usually defined as NO + NO2 and its behavior is similar to a tracer due to the fast interconversion between NO and NO2 [39]. However, NO2 is reactive and NO2 concentrations are estimated using the NO2/NOx ratios obtained from mesoscale simulations at ground level in the mesoscale cells where the microscale domains are located. This assumption, although less accurate than using a full chemical model, allows the calculation of the annual average concentrations. In addition, NO2/NOx ratios were also used for estimating NO2 concentrations by [17,22], although in these cases the ratios were computed using concentrations recorded at AQMSs. In addition, the thermal effects are assumed to be negligible with respect to dynamical effects and the concentration for each hour depends only on emissions, wind speed, and background concentration for that hour. The methodology has three main steps (Equation (4)): (1) for each hour, a CFD simulation (NOx_CFD) is selected according to the wind direction obtained from the mesoscale meteorological simulations for that hour (WD(t)); (2) NOx concentration distributions are modified considering the ratio between reference velocity in CFD simulations (uref(CFD)) and the reference velocity obtained from the mesoscale simulations for that hour (uref(meso,t)) and (3) background concentration simulated by the mesoscale air quality model for that hour (NOxbackground(t)) is added.
Finally, the distributions of NO2 are computed as, , where NO2(meso,t) and NOx(meso,t) are NO2 and NOx concentrations at ground level obtained from mesoscale simulations corresponding to the mesoscale grid cells where the microscale domains are located. In this study, the friction velocity is adopted as reference velocity. The authors of [22] found better agreement with observations using this approach than using wind speed at a certain height as reference velocity.
The CHIMERE chemistry transport model [40] was used to estimate the background concentrations of NO and NO2. The model was applied to a domain covering an area of approximately 310 km × 330 km centred on Madrid with a spatial resolution of 0.03° × 0.03° (approx. 3 km × 3 km). This simulation was nested within simulations for two larger domains (Iberian Peninsula at 0.1° × 0.1° and Europe at 0.15° × 0.15°; described in [2]). The meteorological fields were adapted from simulations at the European Centre for Medium-Range Weather Forecasts, ECMWF, (www.ecmwf.int (accessed date: 10 January 2022)) known as the Integrated Forecasting System (IFS), for 2016, and obtained from the MARS archive at ECMWF through the access provided by AEMET for research projects. Boundary conditions for the European domain were taken from LMDZ-INCA [13] and GOCART [14] global model climatology. All the simulations used the same boundary conditions. In the air quality mesoscale modelling approach, the pollutant concentrations estimated by CHIMERE model were combined with observations at AQMS in order to minimize biases and uncertainties in the concentration estimates. A similar mesoscale modelling configuration was successfully used by [2]. An evaluation of the simulation used in the present study using the hourly NO2 concentrations recorded at the urban background AQMSs within the domain gave the following statistical metrics:

Evaluation of the Modelling Appproach
Firstly, the modelling approach using CFD simulations is evaluated for 2016. For this purpose, the time series of NO2 concentrations recorded at the three AQMSs are compared with modelled NO2 concentrations. It is noteworthy that previous evaluations of spatial distributions of NO2 concentrations computed with a similar methodology were carried out for the same numerical domains. The authors of [15] estimated the spatial distribution of time-averaged NO2 concentrations over approximately three weeks in two different time periods (January-February 2011 and November 2014) for the EA numerical domain. The time-averaged NO2 concentrations were successfully compared with data from two experimental campaigns deploying 26 and 95 passive samplers in 2011 and 2014, respectively. The authors of [38] estimated the spatial distribution of time-averaged NO2 concentrations over approximately three weeks in February 2015 for the PE numerical domain. The modelling approach showed a good agreement with the time-averaged NO2 concentrations recorded by 72 passive samplers distributed throughout the study area. Unlike previous studies, where background concentrations were derived from the closest urban background AQMS, the present work obtains background NO2 concentrations from air quality mesoscale simulations. In fact, background concentrations of emission reduction scenarios can only be derived from mesoscale simulations. Figure 2 shows the time series and scatter plots of modelled and observed NO2 concentrations at the locations of the AQMSs. The mean diurnal concentration profiles are computed for the three AQMSs in order to understand and compare the temporal distribution of observed and modeled concentrations (Figure 3). It can be seen that the average daily behavior of the observed concentrations is captured by the model. In addition, the performance of the modelling approach is quantified by means of statistical metrics such as the normalized mean-square error (NMSE), the fraction of prediction that are within a factor of two of the observation (FAC2), the correlation coefficient (R), and the mean fractional bias (MFB) (Equations (6)-(9)).  [15,22,28,33,41]. Note that MFB is only shown for the PE case because for the EA and PC areas the traffic emission factors were modified in order to give MFB = 0 for 2016 scenario.

Impacts of Emission Reduction Scenarios on NO2 Concentrations
In order to assess the impact of WEM2030 and WAM2030 scenarios projected for 2030 on annual average NO2 concentrations within the three neighborhoods, annual average NO2 concentrations for 2016 are compared with the concentrations obtained for the WEM2030 and WAM2030 scenarios under the same meteorological conditions. Future projections have not considered future temperature changes, ozone projections, or changes in the surrounding landscape around the city. Both the spatial average and the spatial distribution of annual NO2 concentrations are studied for the three emission scenarios. In addition, the extension of high-polluted zones that can affect the population exposure is determined within the study neighborhoods. The areas above the annual average limit value for NO2 (40 µg m −3 ) are estimated for each neighborhood and the ratio between this area (Aover40) and the total study area (Atotal) for each neighborhood is investigated.  Spatially-averaged NO2 concentrations for the PE study area are estimated for the three emission scenarios (Figure 5a). For 2016, NO2 is 58.5 µg m −3 , which is well above the annual average limit value established in EU legislation. The emission reductions of both scenarios succeed in reducing the spatially-averaged NO2 concentration below the limit values, although for the WEM2030 scenario, the concentration is very close to this limit (39.4 µg m −3 ). However, due to the strong spatial variability, areas with annual average NO2 concentrations above 40 µg m −3 are found within the study zone, not only for the 2016 scenario (the total area is above the limit value), but also for WEM2030 and WAM2030 scenario (Figure 5b). In particular, this area is large for the WEM2030 scenario, corresponding to 34.7% of the total area. For the WAM2030 scenario, in spite of the spatial average of NO2 concentration being well below the limit value, 4.2% of the study area has a concentration above 40 µg m −3 . This zone is located in the freeway at the exit of the tunnel. Figure 5. (a) Spatially-averaged annual NO2 concentrations for each study area for the three emission scenarios. (b) Ratio (in %) between the area above the annual average limit value for NO2 (40 µg m −3 ) (Aover40) and the total study area (Atotal) for each neighborhood and for the three emission scenarios.
Figures 4d-f show the maps of annual average NO2 concentrations in the EA study area for the three emission scenarios. In this area, strong gradients of concentration are also observed for all emission scenarios with the highest-polluted zones in the streets and avenues between the buildings and the lowest concentrations in El Retiro park in the south of the study area. Despite the emission reductions, for the WEM2030 scenario, high concentration values are found in some streets in the northern part of the domain. A more pronounced reduction of concentrations is seen for the WAM2030 scenario.
The spatially-averaged NO2 concentration over the EA study area (Figure 5a) is 58.3 µg m −3 , very similar to that of the PE study area. In this case, the emission reductions of the WEM2030 scenario do not succeed in reducing the spatially-averaged NO2 concentration (41.2 µg m −3 ) below the annual average limit value. Unlike the WEM2030 scenario, the spatially-averaged concentration is well below the limit value for the WAM2030 scenario (27.3 µg m −3 ). In this study area, the spatial variability is more pronounced and the areas above 40 µg m −3 for the WEM2030 and WAM2030 scenarios are larger than those found in the PE area. Aover40 is 41.2% of the total area for the WEM2030 scenario and the 12.3% for the WAM2030 scenario (Figure 5b). Figures 4g-i show the maps of annual average NO2 concentrations in the PC study area for the three emission scenarios. In this area, the largest concentrations are found in Gran Vía avenue and the lowest values in the pedestrianized zones. The spatially-averaged NO2 concentrations are 44.6 µg m −3 , 28.2 µg m −3 , and 19.9 µg m −3 for the 2016, WEM2030, and WAM2030 scenarios (Figure 5a). These values are lower than those in the other two neighborhoods, although for the 2016 scenario, the average concentration exceeds the limit value. Regarding the percentages of areas with concentrations above 40 µg m −3 , the emission reduction measures for the WAM2030 scenario succeed in decreasing them to 4.7% of the study area for the WEM2030 scenario and to 0% for the WAM2030 scenario (Figure 5b).

Discussion and Conclusions
The numerical methodology based on CFD simulations and the use of mesoscale modelling to obtain hourly atmospheric conditions (wind direction, reference wind velocity, and background concentration) permits the evaluation of the impacts of emission reduction measures on air quality under the same conditions and at high spatial resolution. The methodology followed in the present study focuses on NO2 since the European limit value for NO2 was exceeded at several AQMSs in Madrid-including the three stations studied here [32]-and this study aims to investigate the effectiveness of emission reduction scenarios on these NO2 exceedances. In addition, this methodology could be extended to other non-reactive pollutants such as PM10 or PM2.5 in future studies. In this study, reduction emission measures implemented at a national level are evaluated. The results estimates that the scenario with additional measures (WAM2030 scenario) projected for 2030 reduces NO2 concentration to such an extent that the spatially-average of NO2 concentration over the entire neighborhoods does not exceed the EU annual average limit value for NO2 for the three study areas. For the WEM2030 scenario, the NO2 spatially-averaged over the EA study area exceeds the annual limit value and in the case of PE it is just below this limit. This information is helpful for policymakers because not only are the reductions of spatially-averaged annual mean NO2 concentrations (similar to concentrations at coarse spatial resolution) estimated, but also the areas within the neighborhoods that have annual mean NO2 concentrations above the limit value (despite the spatially-averaged concentration not being exceeded). Therefore, policymakers can use this information to adapt and/or design new local measures for improving air quality in those areas with exceedances (e.g., limitation of the traffic or vegetation barriers for reducing the exposure of pedestrian). From another perspective, this study highlights a weakness of the EU limit values, since the directive does not specify at which spatial scale the NO2 concentrations must be estimated. This leaves several questions open, like whether the limit value can be exceeded over a small area and if so, how small this area can be. This is in agreement with [2], who used mesoscale simulations for the Iberian Peninsula, that, of the three air quality zones (areas designated under European legislation for managing air quality) that exceeded the annual mean limit value for NO2 for the 2016 scenario, two air quality zones exceeded the limit value for the WEM scenario and none for the WAM scenario. However, as shown by [4], to estimate population exposure and air quality assessment, it is important to take into account the spatial variability of NO2 concentrations within each neighborhood. Despite the annual mean limit value not being exceeded in any of the study neighborhoods in terms of the spatially-averaged NO2 concentrations for the WAM2030 scenario, there are areas with concentrations above 40 µg m −3 within two neighborhoods (4% for PE and 12% for EA). This indicates that if the population was homogenously distributed over the area, 4% of the people in the PE study area and 12% of the people in the EA study would be exposed to a concentration above the limit value. For the WEM2030 scenario, these percentages are larger with values of 35%, 41%, and 5% for the PE, EA, and PC study areas, respectively. All of these percentage depends on the dimensions of the study areas and this issue is important to bear in mind for the interpretation of the results. However, the spatially-averaged concentrations over the microscale domains are similar to the values provided for mesoscale models where the spatial resolution is similar to the size of the numerical domains. Therefore, this study shows that the annual mean limit could be exceeded in some areas within the mesoscale cells in spite of the spatially-averaged NO2 concentration being below the limit value. This issue highlights the importance of studying air quality at the street level. In future studies, the movements of people within the study areas could be considered to estimate the personal exposure. The authors of [3] used CFD simulations and microsimulations of pedestrian movements to estimate the pedestrian exposure. Combining the methodology used in the present study and similar computations of pedestrian movements, future studies could estimate the impact of strategies to improve air quality in terms of personal exposure.
A limitation of this study is the computation of emission reductions for WEM2030 and WAM2030 scenarios for CFD simulations. The same emission reduction factors computed for mesoscale simulations for traffic sector has been applied to the CFD simulations, but the traffic emissions may have a more complex behavior inside the city [42]. However, the aim is to investigate the same emission reduction scenarios at street level that were previously studied at a national level [2] and so the approach used is considered appropriate. Regarding meteorological conditions, modelled concentrations for 2030 used the meteorology for 2016 because the objective is to only evaluate the impacts due to the changes in emissions.
This methodology can be used to identify highly-polluted areas within the neighborhoods and aid the design of new local strategies focused on decreasing pollutant concentrations or reducing population exposure in those areas. In addition to the measures studied in this paper (scenarios planned at national level), other local measures are being introduced in the cities such as areas with traffic restrictions (Low Emission Zones, LEZ) [43], nature-based solutions [44], or other mitigation passive measures [45]. The methodology used in this study could also be used to assess the impacts of these measures or a combination of several local measures with the scenarios studied here on the annual average NO2 concentrations within the streets and the exceedances of the annual average limit value.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. the NAPCP and related discussion. We also acknowledge MITERD for providing data from air quality stations and information. We are also grateful for the services offered by the European Center for Medium-Range Weather Forecasts (ECMWF), including the provision of meteorological modeling data with thanks also to AEMET for managing access to this information. The authors acknowledge the useful comments from the anonymous reviewers that helped improve the manuscript.