The Italian National Air Pollution Control Programme: Air Quality, Health Impact and Cost Assessment

: Air pollution is the primary environmental cause of death globally. To improve air quality and reduce health impacts, the National Emission Ceilings Directive requires Member States of the European Union to provide National Air Pollution Control Programmes, including emission reduction measures aimed to achieve binding commitments for the years 2020 and 2030. Integrated assessment models are pivotal to assess the reduction of pollutants concentrations determined by measures implemented or foreseen for emission reduction. Here we discuss scenarios elaborated for year 2030 in the Italian National Air Pollution Control Programme, considering 2010 as reference year. The two scenarios, “With Measures” and “With Additional Measures”, show a signiﬁcant reduction of the pollutants concentration, namely PM 2.5 , NO 2 and O 3 . The scenarios are here also used to provide an integrated approach for calculating the effect of the program on health impacts (mortality) and related costs. Avoidable attributable cases and associated costs are here reported at both the national and regional level and provide a signiﬁcant framework to assess air-pollution reduction measures with an integrated approach. The procedure proposed may be therefore further developed and applied to assess the overall positive beneﬁts (environmental, health and economic) determined by air-pollution control plans or other integrated policies targeting air quality, energy and climate goals.


Introduction
In order to contribute to the reduction of air pollution and its negative impacts on human health and environment, in 2016, the National Emission Ceilings Directive (NECD) [1] set national emission reduction commitments for atmospheric pollutants, namely sulfur dioxide (SO 2 ), nitrogen oxides (NO X ), particulate matter with diameter of 2.5 µm or less (PM2.5), non-methane volatile organic compounds (NMVOC) and ammonia (NH 3 ), to be fulfilled by 2020 and by 2030. These reductions must be achieved through the adoption and implementation of a National Air Pollution Control Programme (NAPCP), drawn up following specific guidelines [2]. The NAPCP is required to include the quantification of the effect of Policies and Measures (PaMs) aimed at reducing emissions and improving air quality, both for already adopted PaMs (scenario "With Measures", WM) and for potential additional PaMs (scenario "With Additional Measures", WAM) required to comply with the NECD and other relevant legal obligations. Member States are also encouraged to report projected cost-benefit estimates.
As in Europe, air-pollution control plans are key tools used worldwide by administrations in charge of air-quality management, with concurrent responsibilities between central and regional levels. For example, in the USA, Canada and Brazil, the component states/provinces/territories have the main responsibility in adopting Air Quality Implementation Plans, supported by the Federal administration, if needed [3][4][5]. China's Clean Air Action, implemented since 2013, has been fostering urgent mitigation measures for air pollution in megacities [6,7]. In India, the 2019 National Clean Air Programme foresees the formulation of city-specific action plans for 102 cities in exceedance of national air-quality standards [8].
In the majority of control plans, the evaluation of the effectiveness of planned measures for emission reduction, at national and/or regional scale, is based on modeling studies, where the relative impact on pollutants' atmospheric concentration of alternative scenarios can be compared by varying the contribution of specific emission sources in agreement with the proposed intervention actions. In this respect, a worst case scenario can also be modeled, considering no new actions at all, as well as best case scenarios, where several additional measures are implemented. The comparison of future pollutant concentrations (modeled) to past reference concentrations (observed or modeled) might inform on the variations of air quality and hence on the effectiveness of the studied measures [9][10][11]. As the reduction of greenhouse gases has an impact on air quality, the integrated evaluation of climate change effects is becoming a standard practice in air-quality planning, including the NECD, and modeling assessment [12][13][14].
Among the various effects of air pollution, adverse impacts on human health are the most relevant and studied, in terms of both direct outcomes (mortality, morbidity) and indirect impacts (health and social costs). Scenario assessments of air quality are becoming largely used to derive also the health impact and related cost assessment, thus providing quantitative figures of the consequences of air-quality policies and plans [15][16][17].
At the national level, the literature on this topic is scant so far, and only a few studies are available that consider both health and economic impacts of emission reduction measures [18][19][20]. The VIIAS project (Methods for the Integrated Assessment of Environment and Health Impact of Air Pollution, www.viias.it) carried out the health-impact assessment (HIA) of nitrogen dioxide (NO 2 ) and PM 2.5 by using Italian national model scenarios for 2020 [21], while recently we modeled the effect of a subnational plan (Po Basin Agreement) on national air-quality levels and exceedances of limit values in 2020 and 2030 [22]. The effect of regional-scale air-quality plans on air quality in 2024 in the Po Valley were evaluated by Raffaelli et al. [23], while recently De Angelis et al. [24] combined a source apportionment procedure with integrated assessment modeling to evaluate the specific effect of the Po Basin Agreement on PM 10 concentrations in the Lombardy region in 2020. Significantly, Bert et al. [25] reported recently on the urgency to perform an institutionalized HIA for Italy at the national level, in order to provide a coherent framework also for reaching United Nation sustainable development goals.
In this manuscript we report, for the first time to our knowledge, the assessment carried out for the Italian NAPCP [26] on the effect of PaMs on air quality (with specific focus on NO 2 , PM 2.5 and ozone (O 3 )) and provided the related human health impacts (by considering mortality, by cause), and related monetary costs, in 2030 in Italy, considering two alternative future scenarios, WM and WAM. These selected future scenarios are part of the officially submitted national control program for the NECD. The procedure applied in this manuscript, by evaluating both the environmental effects and the associated health impacts and relative costs, may provide additional relevant information to be used by decision makers when adopting air quality plans and in evaluating integrated policies on energy, climate and air pollution.

Materials and Methods
2.1. Integrated Assessment Model for Impacts of Air Quality: National Integrated Model to Support International Negotiation on the Atmospheric Pollution (MINNI) The model assessment was conducted with the MINNI system (National Integrated Model to support International Negotiation on the Atmospheric Pollution), an Integrated Assessment Model (IAM) suite developed by ENEA with ARIANET (private company) and IIASA (International Institute for Applied Systems Analysis) on behalf of the Italian Ministry for Environment. The system was widely described in previous papers [27][28][29]; for the sake of brevity, hereinafter the main features of the model are reported. In MINNI, atmospheric science is linked to impacts of emission abatement measures on human health and ecosystems and related costs, through different independent and interconnected components: the AMS (Atmospheric Modeling System) and the GAINS-Italy (Greenhouse Gas and Air Pollution Interactions and Synergies Model over Italy) model. The AMS computes gas and aerosol transport, diffusion and chemical reactions in the atmosphere with the three-dimensional Eulerian model FARM (Flexible Air Quality Regional Model, [30,31]), with proper emission and meteorological data in input. The system application on Italy has a horizontal spatial resolution of 4 km and hourly temporal resolution. The GAINS-Italy model elaborates air pollutant and greenhouse gas (GHG) emission scenarios from 1990 to 2050 with a 5-year time interval, explores different cost-effective multi-pollutant emission control strategies to meet human health and/or environmental targets on air quality and calculates air quality and impacts on health and ecosystems. In GAINS-Italy, simplified mathematical relations between emission and concentrations/depositions (Atmospheric Transfer Matrices) are obtained through ad hoc AMS simulations and used for fast-response assessment of effects of emission reductions. In this study, MINNI was deployed with complete AMS annual simulations in three different scenarios (2010, 2030 WM and 2030 WAM), to obtain concentration fields of NO 2 , PM 2.5 and O 3 at 4 km resolution, which were used (outside GAINS-Italy) for the following health-impact and cost assessment.

NAPCP Scenarios Modeling
The MINNI system is extensively used in evaluating emission reduction measures, to check future compliance both to emission targets, as for the new NECD for the years 2020 and 2030, and to air-quality standards for ambient concentrations [32]. All simulations were run on the ENEA CRESCO High Performance Computing Infrastructure [33].
The preparation of an emission scenario with the GAINS-Italy model requires the definition of levels of both energy and non-energy anthropogenic activities and of a control strategy, detailed by sector, activity, abatement technology and pollutant. In elaborating the NAPCP, two emission scenarios (WM and WAM) were prepared to check the NEC emission target compliance in the years 2020 and 2030. The baseline reference WM scenario includes all policies and measures put in place before 2015 [34], while the WAM scenario is based on the last National Energy Strategy [35] that enables the achievement of the national objectives regarding energy efficiency, GHGs and renewable sources. Both the energy scenarios were elaborated with the TIMES-Italy model (The Integrated MARKAL-EFOM1 System/EFOM (Energy Flow Optimization Model) [36][37][38]). The emissions elaborated with the GAINS-Italy model were harmonized with the two versions of the official national emission inventory, which were the most recent available when the NAPCP was elaborated [39,40]. The years 2005 and 2010 were considered as base years for harmonization, being the most updated baseline inventories available when the first WM scenario was developed.
The main additional measures are presented in the following Table 1, while more details can be found in the draft version of the NAPCP [41].
The air quality scenarios were elaborated by the AMS, whose configuration adopted for the Italian NAPCP is summarized in Supplementary Materials Table S1. It is important to underline that the variations observed between the years 2010 and 2030 are linked to emission variations having kept constant the meteorological year, as required by the NECD [1]. Therefore, 3 annual simulations were conducted, each fed by the same annual meteorology (year 2010) and by an individual annual emission input (2010, 2030 WM and 2030 WAM). As indicated in Supplementary Materials Table S1, several sources of natural emissions were considered.

Health Impact Assessment
Attributable cases for each health endpoint and for each pollutant were calculated in the open-source statistical software R [42], according to the following function [43]: where M 0 represents the baseline mortality, P j is the exposed population resident in the grid cell j th , β Poll is the concentration-response coefficient for the specific pollutant of interest and ∆ Polj is the variation of pollutant concentration above a selected threshold, specific for each pollutant, at the grid cell j th . Here we selected the following counterfactual values for PM 2.5 , NO 2 and O 3 , respectively, 0, 20 and 70 µg/m 3 . NO 2 and O 3 values were selected according to the limits for health protection suggested by World Health Organization (WHO) when performing a health impact assessment on air pollutants. For PM 2.5 , the counterfactual value was selected by considering that no safe concentration for this pollutant can be assumed from the latest epidemiological analyses [44]. β Poll is calculated from the relative risk (RR) specific for each pollutant and each cause of death. Here we considered the International Statistical Classification of Diseases and Related Health Problems 10th Revision (ICD-10) codes for selection of relevant causes of death. The attributable cases were calculated by considering the RRs proposed by WHO for long-term effects, overall mortality and cause specific of air pollution (References [45][46][47][48] and Table 2), and they were compared with values obtained by considering national RRs for PM 2.5 and NO 2 resulting from the national epidemiological study MED-HISS (Mediterranean Health Interview Survey Studies: longterm exposure to air pollution and health and surveillance) [49] (Table 3). In this study, we considered the baseline mortality for all the population of age above a threshold (30 years in References [45][46][47][48], and 35 years in Reference [49]), using values specific for each Italian region. These data were calculated starting from the Italian National Institute of Statistics (ISTAT) public data, available for 2010, which is the reference year we selected in our simulation. For future scenarios, attributable impacts were calculated by considering a static population, i.e., constant overall distribution among the different class ages and constant baseline mortality. Moreover, relative risks were considered constant in the future. The differences obtained are therefore related to the expected air-pollutant-concentrations reductions and do not consider possible variation in population structure (i.e., increasing ageing of the population).

Cost Assessment
In the present study, we monetized the air pollution economic impacts firstly estimating the Years of Life Lost (YOLLs) [50] and then attributing an average value of a life year (VOLY), selecting costs widely used in previous EU-wide impact assessments [51].
Briefly, YOLLs associated to the attributable cases were calculated with the following equation: where E i are the number of attributable cases per pollutant associated to the i th grid cells belonging to the Italian Regions and AYL the average years lost, set at 10.3 years for all studied pollutants, according to Reference [50].
The YOLLs were valued considering the costs reported in Table 4 and adjusted for Italy and for each Italian Region with the following equation: where the VOLY ref derives from Holland (2014), is the income elasticity of 0.8 (as suggested in Reference [51]) and GDP j and GDP EU (Gross Domestic Product) per capita refers to Italy/Italian Regions and to Europe for the year 2010, respectively [52]. Income elasticity is 0.8 for adjustment The costs associated to attributable cases for all natural causes premature mortality induced by exposure to air pollutants were calculated according to Reference [53], by multiplying the regional YOLL for the regional average adjusted VOLY.

Results
In this section, we present results for emission reductions, air quality improvements, health impact and costs. In Tables 8-10 A map of names of Italian regions and cities, as cited in the text, is available at World Atlas [54].

Emission Reductions
The first step in the elaboration of an emission scenario is the harmonization process with the official emission inventory for the base reference years. The results of the harmonization are presented in Table 5, showing that the inventory is well reproduced by the GAINS-Italy model, with reasonable differences in total emissions both for the years 2005 and 2010. The emission scenarios were elaborated both for 2020 and 2030, in order to check compliance with NECD emission targets at the required years. Due to the COVID-19 pandemic, we consider the current 2020 scenario not reliable. The lockdown measures adopted in Italy during 2020 are too recent, variable in space (administrative regions) and unstable in time to allow an estimate of the final impact on emission reductions and air-quality variation, and we therefore do not report results for this modeled year. Though the economic crisis due to the pandemic is still ongoing, first observed signals during the lockdown exit process indicate a "rebound effect" with recovery of activity levels (road traffic, industrial production and energy consumption [55]); therefore, we assume that 2030 emission scenarios, whose results are summarized in Table 6, can still be considered valid.  SO 2 emissions show a decrease from 2010 to 2030 of 48% and 62% in the WM and WAM scenarios, respectively, mainly driven by stationary sources, in particular by power plants, and the maritime sector. The same sectors mostly contribute to the WAM-WM variation. Compared to the year 2010, a large reduction is also foreseen for NO X emissions, with a decrease of 53% and 62% in the WM and WAM scenarios, respectively, driven by the road transport sector. The additional reduction in the WAM scenario is due to maritime and road transport sectors. The PM 2.5 emission variation between WM and WAM (13%) is lower than for the previous pollutants. This reduction is driven by the civil sector, which remains the main contributor also in the WAM scenario, followed by the off-road transport. The NMVOC variation mainly derives by the solvent use sector, due to a revision of the emission factor and not to an additional measure [40]. Ammonia is the pollutant with the lowest reductions, both between 2010 and 2030 and between WM and WAM. The emission reduction at year 2030 mainly derives from the incorporation of urea-based fertilizers in the agriculture sector, followed by reductions in emission from livestock activities.
Emission reduction goals fixed by the NECD are achieved for SO 2 in the WM scenario and for all pollutants in the WAM scenario (Table 7), but with a narrow margin.

Air-Quality Improvement
The 2010 annual average NO 2 map (Supplementary Materials Figure S2) shows the highest concentrations in large urban areas Milan, Turin, Rome, Naples and medium cities within and outside the Po Valley, where the combined effect of emissions from domestic heating, urban and highway road traffic is largest. By 2030, a widespread reduction in concentrations is observed, driven by a major renewal of the vehicle fleet due to the introduction of the Euro 6 Phase 2 standard. The WAM scenario leads to a visible reduction of the concentration at the hotspot in the center of Milan metropolitan area, possibly due to the increase of the share of electric vehicles for private urban mobility.
Reference scenario maps for annual average of PM 2.5 and April-September average of daily maxima of 8 h running means for O 3 (MDA8) (Supplementary Materials Figure S2) show only minor spatially defined hotspots. This is due to atmospheric chemistry, generating O 3 and secondary PM 2.5 also far from local emission sources, and transport at regional scale (tens-hundreds kilometers). The WM and WAM scenarios introduce extensive and significant reductions of PM 2.5 , due to reduced emissions of both primary PM (planned renewal of biomass heating systems) and precursor species of secondary PM, mainly SO 2 , NOx and NMVOCs. These two latter species are precursors also of O 3 , and their combined decrease is reflected in lower O 3 concentrations in future scenarios.
The historical model (2010 baseline scenario) was not corrected with assimilation of corresponding observed values; therefore, concentration results should be analyzed mainly in terms of their evolution and variation between baseline and 2030 alternative scenarios, as further discussed in Section 4. Tables 8 and 9 show the aggregated results obtained at the national-level for mortality in 2010 and in the two 2030 scenarios. Mortality data are reported by considering average, min and max RR, as reported in Table 8, to account for the possible lowest and highest potential attributable cases. Table 8 shows the variations in absolute and relative mortality between 2010 and 2030 and between 2030WM and 2030WAM to better clarify the possible expected effects of these two alternative scenarios. In 2010, attributable cases for PM 2.5 (all natural causes), NO 2 (all natural causes) and O 3 (respiratory disease) are 58,867 (lower and higher value equal to 35,379 and 83,670 cases), 11,769 (min 6566-max 17301) and 2692 (min 945-max 4702), respectively. Considering cause specific deaths, calculated only for PM 2.5 , the highest share of additional burden of death is due to cardiovascular diseases with 46,960 cases (about 72% of total cases), while respiratory diseases and lung cancer account for 7396 (11%) and 4040 (7%), respectively.
The WAM scenario leads to significant additional mitigations for all pollutants, with about 24201 avoided cases for all-cause PM 2.5 mortality, 10976 for NO 2 and 967 for O 3 considering the 2010 reference scenario. The relative reduction (close to 100%) of airborne NO 2 concentration expected in 2030 determines the highest significant reduction among the analyzed pollutants (−93% of attributable cases). Noteworthy, the additional measures of the WAM scenario could lead to significant reduction in mortality for all the pollutants here analyzed. The highest reduction is to be expected on NO 2 attributable deaths (−53% WAM vs. WM scenario), while the highest number of avoidable cases is related to additional measures implemented to reduce PM 2.5 (3970 avoided cases for cardiovascular mortality accounting for an additional 15% reduction for this specific cause of death). Results at regional level are not here discussed but are reported in Supplementary Materials Table S2.
The maps of attributable deaths and variations among the modeled scenarios, for each pollutant, are shown in Figures 1-3, for PM 2.5 , NO x and O 3 respectively. PM 2.5 attributable deaths are spread all over the nation, with hotspots on large and medium urban areas. The reduction of PM 2.5 attributable deaths in 2030 is significant especially in the Po Valley and the urban areas of Florence, Rome and Naples. The WAM scenario underlines additional effects on the Milan metropolitan area. For NO 2 and O 3 , mortalities occur in a small fraction of the territory (Figures 2 and 3). For NO 2 Table 9 underlines the importance of also using nationally relevant RRs for a HIA procedure. Compared to international RRs, national ones reduce the expected impacts of air pollutants; however, the relative burden of PM-induced lung cancer events is increased. The differences here reported may also be related to the differences in the population age considered in international and national RRs evaluation procedures. While international RRs are derived from epidemiological data of all population over 30 years old, Italian values are derived from a population of age over 35 years old and from the municipality of Rome.

Mortality Rate
Mortality rates are reported as the average of non-zero values obtained on grid cells of the computational domain (Table 10). Like for the previous section, Table 10 shows the absolute and relative mortality rate variations between 2010 and 2030 and between WM and WAM.   Figure 4. Rates of mortality are more homogenously spread over the national territory if compared to the attributable cases (Figure 1 (2010) vs. Figure 4 (2010)). Mortality rates follow the spatial pattern of PM concentration (Supplementary Materials Figure S1), with a large portion of the national territory that is characterized by more than five attributable cases per 10.000 inhabitants and with the Po Valley and other municipalities with rates above 12 events per 10.000 inhabitants. The expected reduction in mortality rate in 2030 scenarios is higher for the Po Valley area but significant all over Italy. Similar considerations can be made for the mortality rate of NO 2 (Supplementary Materials Figure S3), and O 3 (the maps are not reported since the values of standardized mortality are always below 1) and are not here furthermore discussed.

Costs
Air pollution externalities, here calculated as monetary cost of YOLL, are reported as avoided costs by region, comparing the 2030-WM scenario with the 2010 reference scenario and the WAM scenario with WM (Table 11). The total benefits associated to the reduction in attributable cases for the pollutants here considered in the 2030 WM scenario, compared to the base year 2010, reach 29,691 million € for Italy, with clear differences among regions and pollutants, with an additional total benefits of 3415 million € in the WAM scenario. In the 2030-WM scenario, the benefits due to the avoided deaths respect to the year 2010 have an incidence on the national GDP2010 of 1.84%, with an extra incidence of 0.21% in the WAM scenario.
It is evident from the Supplementary Materials Table S4 that, in the year 2010, PM 2.5 causes the majority of total damage costs, with differences among regions. On a national level, PM 2.5 contributes to 80% of total damage costs, followed by NO 2 (16%) and O 3 (4%). These average values considerably differ among regions. The lowest contribution of PM 2.5 , for example, is observed in Lombardia (68%) and Lazio (70%), due to a higher share of NO 2 -related costs, 30% and 26%, respectively. In most regions, almost the whole contribution to total costs (around 95%) can be ascribable to PM 2.5 . The O 3 contribution spans from 3% to 7%. The picture of the 2030 scenarios retraces what was described for year 2010, as shown in Supplementary Materials Table S5. On average, PM 2.5 contributes to 66% and 72% of total benefits of the WM and WAM, respectively, followed by NO 2 (31% and 25%) and O 3 (3% in both scenarios). The PM 2.5 and NO 2 contributions to total benefits of Lombardia and Lazio are almost the same, indicating a large decrease in the attributable deaths to NO 2 , as already underlined in Section 3.3.1. Table 11. Total and per pollutant benefit-costs (mln €) of the 2030-WM scenario, compared to the base year 2010, and of the 2030-WAM scenario, compared to the 2030 WM; incidence of the benefits as percentage of the regional and national GDP2010.

Discussion
In the present study, we reported an integrated procedure to evaluate the expected air pollutant emissions and concentrations deriving from the first National Air Pollution Control Programme and to quantify the possible attributable adverse health cases (here we focused on mortality, but morbidity can be evaluated, as well) with the relative quantification of the expected societal costs (expressed as YOLL related costs). This is, to our knowledge, the first application of this kind of procedure on Italy.
Integrated assessment models are a potent tool for assessing the potential effects of new policies and regulations aiming at reducing air pollution, energy and climate and their relative impact on health and the environment. Nonetheless, it is relevant to recall that modeled air-pollutant concentrations may suffer from biases introduced by, among others, model-specific approximations, limits in the emission sources characterization and incorrect parametrization of meteorological processes. A way to overcome this biased estimation would be to apply a procedure to correct the modeled concentration with the data from monitoring stations, as done in some studies on long-term projection of future air quality [56][57][58][59], taking into account climate change. However, other studies [60,61] focus on the relative change of concentrations induced by emissions mitigation between the modeled future and reference years, rather than on the absolute values of pollutants concentration. This is also the perspective of the NECD, requiring explicitly to model the effects of emission policies, therefore keeping fixed other parameters influencing model concentrations [2]. In the present manuscript, a bias-correction procedure was not applied, as this would have altered the estimates of variations in concentrations.
The absolute mortality numbers here reported may be underestimated, but their variations among the three scenarios, as well as the geographical distributions, are reliable and give useful indications of the effects of the NAPCP on human health. The reliability of the procedure we applied is confirmed by the number of attributable cases for "all causes natural mortality" due to PM 2.5 exposure in 2010, resulting in fair agreement with values previously reported by the European Environment Agency (EEA) for 2011 (64,544, min 42,650-max 84,475) ( [62]).
At the Italian level, the VIIAS study [21,63] has found 21,524, 11,993 and 1858 premature deaths attributable to PM2.5, NO 2 and O 3 , respectively, in 2010. With respect to the present study, where the number of events is substantially the same for NO 2 and O 3 , in VIIAS, the lower number of events from particulate matter is partly due to lower model concentrations and partly to a different counterfactual value (10 µg/m 3 in VIIAS and 0 µg/m 3 in our study) selected for attributable cases evaluation. The counterfactual level selected here for fine PM take into consideration the latest epidemiological results [44] that show that even at concentrations close to zero the health risk (HR) associated to long-term exposure is not negligible. More in general a counterfactual value equal to zero, accounting for all the attributable deaths to the air pollutant, is to be considered as the worst case scenario for the expected impacts, indicating the maximum impacts and related costs to be mitigated and avoided. At the regional level, Carugno et al. (2017) [64] evaluated the number of attributable death for PM 10 short-term exposure in Lombardy main cities: Although the absolute numbers reported in the study are not comparable with those we reported for differences in the study rationale and procedures (such as pollutant considered (PM 10 vs. PM 2.5 ), counterfactual concentrations used (20 vs. 0 µg/m 3 ), different population (main cities vs. all region) and air-quality data origin (measured vs. modeled)), the rates of attributable deaths reported are in agreement with those we report here. In particular, the Milan urban area is characterized by the highest rate (1.53% for 100.000 inhabitants as average 2007-2010 in Reference [64] and between 15 and 20% for 10.000 inhabitants here) when compared to the rest of the regional territory (1.04 in Reference [64] vs. 12.5-15 here).
Variations of attributable deaths (Figures 1-3) clearly show that the baseline evolution of emission impacts substantially on the number of expected deaths, mainly due to the combined reduction of PM 2.5 and NO 2 . The absolute number of avoided events from PM 2.5 is more than the double of that from NO 2 , despite a lower value of the relative reduction and a similar value of the RR. This is related to the higher exposure of population to PM 2.5 with respect to NO 2 , thus allowing for higher potential in reducing mortality for particulate exposure. This, in turn, is reflected in the high relative reduction of NO 2 mortality: The planned policies reach already almost the maximum feasible effect considering the threshold level of 20 µg/m 3 accepted for this pollutant.
While reductions in mortality rate from 2010 to 2030 are similar to those of absolute mortality for PM 2.5 and O 3 , for NO 2 , the reduction in mortality rate is equal to 26% in the WM scenario and 39% in the WAM scenario, compared to 85% and 93% for WM and WAM absolute mortality reduction, respectively. This means that NO 2 concentrations and their reductions are more related to the density of population and therefore with the source of emission of this pollutant, while for PM 2.5 and O 3 , also less populated areas are expected to be impacted for the more regional distribution of these pollutants, also in relation to their secondary formation processes.
The expected reductions in mortality are reflected in the avoidable societal costs, here calculated as YOLL costs. At the national level, the implementation of the WAM scenario could determine an economical benefit equal to 2% of Italian 2010 GDP, with a major share due to the WM measures (1.8%) and a relative lower contribution to the additional measures (0.2%). At the regional level, the highest share of the economic benefits is expected in regions of the Po Valley, namely Piemonte, Lombardia, Veneto and Emilia Romagna (Table 11). Nonetheless, significant benefits are expected in other regions, such as Campania, Lazio and Toscana (characterized by significant urban areas) and for Abruzzo, Friuli-Venezia Giulia, Liguria, Marche and Umbria, where the benefits, although low as absolute value, may be higher than 1% of the relative regional GDP. Therefore, as to be expected, measures implemented at national level could have beneficial economical outcomes to different regions all over the national territory. The estimate of economic effects of the health outcome is affected by high uncertainty due to the different variables, assumptions and parametrization considered [53,65]. Moreover, there are no available national countryspecific estimates of VOLY, and European estimates were considered in this study in a simplified way, with respect to more sophisticated approaches available in the literature. Nonetheless, the values here discussed are in agreement with recently reported values from, for example, the CE Delft study [53]. The values reported in Reference [53] cover 54 Italian cities (the most densely inhabited), covering a population of around 15 million people, and reach a total damage costs of 20820 mln€ in the year 2015. These differences are possibly related to the different monetary value for each year lost, to population size, GDP or different ways to estimate pollution concentrations (measured vs. modeled); however, both studies underline the large impact of poor air quality on human welfare. Interestingly, a recent study [66] reported the benefits, in terms also of avoided mortality in relation to implementation of different energy technologies at 2030, showing the importance of considering the avoided societal and environmental costs in planning future strategies.

Conclusions
In this study, we reported, for the first time, an integrated assessment of energy, climate and air-pollution policies, analyzing their impact on air-pollutant concentrations and related consequences on human health, with the associated economic evaluation of attributable cases at the regional level. This procedure, although with some limitations due the missing of a bias-correction procedure to the modeled concentration, of country-specific VOLY and of the morbidity evaluation, may be a relevant reference for future studies. The evaluation of externalities associated to pollutants' effects should not be neglected anymore, as it is becoming a key factor in the decision-making process. Our study, therefore, provides-again, for the first time-a preliminary evaluation of the economic benefit of air pollution reduction measures of the first National Air Pollution Control Programme elaborated for the NEC Directive. The avoided socioeconomic costs may be an additional driver for selecting the most appropriate control measures integrating energy, climate and air-pollution policies to reduce air pollutants and save lives.  Table S1: Main characteristics of the modelling setup for the Italian NAPCP simulations,  Data Availability Statement: A publicly available dataset on mortality in Italy was analyzed in this study. This data can be found here: http://dati.istat.it/Index.aspx?DataSetCode=DCIS_ MORTALITA1. Pollutant concentrations data produced by MINNI are available upon request.