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

Abstract

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 significant reduction of the pollutants concentration, namely PM_2.5, NO₂ and O₃. 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 significant 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 benefits (environmental, health and economic) determined by air-pollution control plans or other integrated policies targeting air quality, energy and climate goals.

1. 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₂), 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₃), 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₂) 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₁₀ 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₂, PM_2.5 and ozone (O₃)) 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.

2. 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₂, PM_2.5 and O₃ at 4 km resolution, which were used (outside GAINS-Italy) for the following health-impact and cost assessment.

2.2. 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. Measures adopted in the 2030 “With Additional Measures” (WAM) scenario.

Sector	Name	Description
Power plants	Coal phase–out	Gradual phasing–out of coal power plants to be completed by 2025
Residential/Services sector	Replacement of biomass systems	Renewal of old biomass heating systems with efficient and low–emission technologies
Residential/Services sector	Energy efficiency in buildings	Tighten minimum standards for building (for example, Nearly Zero Energy Buildings)
Transport	Public transport promotion	Promote public transportation to reduce private transport and renew bus fleet
Transport	Electric vehicles	Increase the spread of electric vehicles for private urban mobility
Transport	Renew fleet for freight vehicles	Promote the use of methane/liquefied natural gas (LNG)–powered heavy duty trucks. Promote the use of LNG in maritime transport
Agriculture	Incorporate fertilizers	Incorporate urea–based fertilizers
Agriculture	Ban on new waste lagoons	Ban on constructing new waste lagoons
Agriculture	Slurry	Measures to reduce spread slurry and its incorporation
Agriculture	Spreading of solid manure	Incorporation of manure distributed on the surface

, 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.

2.3. 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]:

$$Case{s}_{Polj}^{mortality}= M_0\times P_j\times \left({\mathrm{e}}^{{\beta }_{Poll}{\mathsf{\Delta }}_{Polj}}-1\right)$$

(1)

where M₀ 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₂ and O₃, respectively, 0, 20 and 70 µg/m³. NO₂ and O₃ 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. Relative risks—World Health Organization. Values are calculated for people aged more than 30 years.

Pollutant Metric	Range of Concentration	Mortality Cause	ICD-10 Codes	RR (95% CI) per 10 μg/m3	Reference(s)
PM2.5 annual mean	All	All-natural causes	A00–R99	1.066 (1.040–1.093)	[45]
		Cardiovascular diseases	I00–I99	1.10 (1.05–1.15)	[45]
		Respiratory diseases	J00–J99	1.10 (0.98–1.24)	[45]
		Lung cancer	C32–C34	1.09 (1.04–1.14)	[45,46]
NO2 annual mean	All	All-natural causes	A00–R99	1.055 (1.031–1.080)	[47]
O3 April–September mean of MDA8 *	>70 μg/m3	Respiratory diseases	J00–J99	1.014 (1.005–1.024)	[48]

* MDA8 = daily maxima of 8–hour running means. ICD-10 = International Statistical Classification of Diseases and Related Health Problems 10th Revision. RR = relative risk.

), and they were compared with values obtained by considering national RRs for PM_2.5 and NO₂ resulting from the national epidemiological study MED-HISS (Mediterranean Health Interview Survey Studies: long-term exposure to air pollution and health and surveillance) [49] (

Table 3. Relative risks—alternative exposure assessment. Values are calculated for people aged more than 35 years.

Pollutant Metric	Range of Concentration	Mortality Cause	ICD-10 Codes	RR (95% CI) per 10 μg/m3	Reference
PM2.5 annual mean	All	All–natural causes	A00–R99	1.03 (1.02–1.05)	[49]
		Cardiovascular diseases	I00–I99	1.03 (1.00–1.06)
		Respiratory diseases	J00–J99	1.04 (0.96–1.13)
		Lung cancer	C32–C34	1.17 (1.09–1.24)
NO2 annual mean	All	All–natural causes	A00–R99	1.02 (1.00–1.03)

). 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).

2.4. 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:

$$YOLL= {\displaystyle \sum }_i{E}_i\ast AYL$$

(2)

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. Cost values considered in the economic health impact assessment (price year 2010).

Outcome	Age Group	Adjusted Reference Value (EU28) (€2010)	Adjusted Value for Italy (€2010)	Reference Value (€2005)
Premature deaths (value of a year of life lost)	>30 years	62,633/144,371	88,930/204,986	57,000/133,000 [50]
Income elasticity is 0.8 for adjustment

and adjusted for Italy and for each Italian Region with the following equation:

$$VOL{Y}_j= VOL{Y}_{ref}{\left(\frac{GD{P}_j}{GD{P}_{EU}}\right)}^{\beta }$$

(3)

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].

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.

3. Results

In this section, we present results for emission reductions, air quality improvements, health impact and costs. In Tables 8–10, reductions from 2010 to 2030 WM and from 2010 to 2030 WAM are related to absolute figures of 2010, while reductions from 2030 WM to 2030 WAM are related to absolute figures of 2030 WM, the latter in order to better show the effect of additional measures on 2030 “baseline” outcomes.

A map of names of Italian regions and cities, as cited in the text, is available at World Atlas [54].

3.1. 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. Difference between Greenhouse Gas and Air Pollution Interactions and Synergies Model over Italy (GAINS-Italy) emission estimates and the national emission inventory submissions 2017 and 2018 (INV_sub2017 and INV_sub2018, respectively), for the years 2005 and 2010, for the five National Emission Ceilings Directive (NECD) pollutants.

	Year 2005		Year 2010
Pollutant	GAINS-Italy vs. INV_sub2017	GAINS-Italy vs. INV_sub2018	GAINS-Italy vs. INV_sub2017	GAINS-Italy vs. INV_sub2018
SO2	−3.65%	−3.85%	−4.11%	−4.41%
NOX	−1.19%	−4.83%	2.30%	−0.23%
PM2.5	6.19%	1.52%	4.03%	0.75%
NMVOC	−0.78%	−8.65%	0.09%	−10.25%
NH3	1.25%	0.62%	0.45%	1.03%

, 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. National emissions by sector for the year 2010, the 2030 “With Measures” (WM) and WAM scenarios and their comparison.

SO2 Emissions
Sector	2010 (kt/yr)	2030 WM (kt/yr)	2030 WM–2010 (%)	2030 WAM (kt/yr)	2030 WAM–2010 (%)	WAM–WM (%)
Stationary Sources	169.95	90.77	−47%	69.09	−59%	−24%
Civil	7.84	5.23	−33%	5.22	−33%	0%
Road Transport	0.68	0.57	−15%	0.44	−35%	−23%
Off–Road Transport	0.40	0.43	6%	0.42	6%	0%
Maritime	29.17	10.59	−64%	3.33	−89%	−69%
Waste	0.38	0.38	0%	0.38	0%	0%
TOTAL	208.41	107.96	−48%	78.88	−62%	−27%
NOx Emissions
Sector	2010 (kt/yr)	2030 WM (kt/yr)	2030 WM–2010 (%)	2030 WAM (kt/yr)	2030 WAM–2010 (%)	WAM–WM (%)
Stationary Sources	235.40	145.58	−38%	127.92	−46%	−12%
Civil	66.41	53.87	−19%	43.28	−35%	−20%
Road Transport	485.52	151.34	−69%	126.70	−74%	−16%
Off–Road Transport	81.85	30.07	−63%	29.52	−64%	−2%
Maritime	98.43	75.14	−24%	37.98	−61%	−49%
Waste	1.20	1.20	0%	1.20	0%	0%
TOTAL	968.80	457.20	−53%	366.59	−62%	−20%
PM2.5 emissions
Sector	2010 (kt/yr)	2030 WM (kt/yr)	2030 WM–2010 (%)	2030 WAM (kt/yr)	2030 WAM–2010 (%)	WAM–WM (%)
Stationary Sources	22.82	14.84	−35%	14.03	−39%	−5%
Civil	113.52	69.17	–39%	61.00	–46%	–12%
Road Transport	27.58	6.45	−77%	5.34	−81%	−17%
Off–Road Transport	14.32	8.80	−39%	3.73	−74%	−58%
Other Sources	17.87	18.35	3%	18.24	2%	−1%
TOTAL	196.11	117.61	−40%	102.33	−48%	−13%
NMVOC emissions
Sector	2010 (kt/yr)	2030 WM (kt/yr)	2030 WM–2010 (%)	2030 WAM (kt/yr)	2030 WAM–2010 (%)	WAM–WM (%)
Stationary Sources	85.07	72.60	−15%	66.35	−22%	−9%
Civil	183.44	106.37	−42%	93.13	−49%	−12%
Solvent Use	394.03	372.98	−5%	331.32	−16%	−11%
Road Transport	190.94	86.15	−55%	52.43	−73%	−39%
Off−Road Transport	83.71	26.47	−68%	21.32	−75%	−19%
Other	63.20	58.18	−8%	45.82	−28%	−21%
TOTAL	1000.39	722.75	−28%	610.37	−39%	−16%
NH3 emissions
Sector	2010 (kt/yr)	2030 WM (kt/yr)	2030 WM–2010 (%)	2030 WAM (kt/yr)	2030 WAM –2010 (%)	WAM–WM (%)
Livestock	313.14	297.26	−5%	289.96	−7%	−2%
Fertilizer Use	50.32	63.78	27%	46.85	−7%	−27%
Transport and Other Sectors	24.97	17.11	−31%	17.69	−29%	3%
TOTAL	388.43	378.15	−3%	354.50	−9%	−6%

, can still be considered valid.

SO₂ 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₂ in the WM scenario and for all pollutants in the WAM scenario (

Table 7. Verification of fulfilment of emission reduction targets set by the NECD for the year 2030.

	2030 Emission Reductions Respect to the Base Year 2005
Pollutant	NECD Target	2030_WM	2030_WAM
SO2	−71%	−73%	−80%
NOX	−65%	−63%	−70%
PM2.5	−40%	−33%	−42%
NMVOC	−46%	−43%	−50%
NH3	−16%	−11%	−17%

), but with a narrow margin.

3.2. Air-Quality Improvement

The 2010 annual average NO₂ 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₃ (MDA8) (Supplementary Materials Figure S2) show only minor spatially defined hotspots. This is due to atmospheric chemistry, generating O₃ 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₂, NOx and NMVOCs. These two latter species are precursors also of O₃, and their combined decrease is reflected in lower O₃ 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.

3.3. Health Impact

3.3.1. Mortality

Table 8. National-level mortality (number of deaths).

Pollutant Metric	Health Outcome	2010 (min–max)	2030 WM (min–max)	2030 WM–2010 (%)	2030 WAM (min–max)	2030 WAM– 2010 (%)	WAM–WM (%)
PM2.5 annual mean	Mortality all–natural causes	58,867 (35,379–83,670)	37,335 (22,608–52,656)	−37%	34,666 (21,013–48,840)	−41%	−7%
	Mortality cardiovascular diseases	46,960 (22,936–72,106)	26,817 (13,392–40,277)	−43%	22,847 (11,427–34,264)	−51%	−15%
	Mortality respiratory diseases	7396 (0–18,948)	4223 (0–10,173)	−43%	3598 (0–8632)	−51%	−15%
	Mortality lung cancer	4040 (1753–6436)	2337 (1037–3642)	−42%	2168 (964–3372)	−46%	−7%
NO2 annual mean	Mortality all–natural causes	11,769 (6566–17,301)	1727 (972–2513)	−85%	793 (449–1149)	−93%	−53%
O3 April–September mean of MDA8	Mortality respiratory diseases	2692 (945–4702)	1851 (654–3211)	−31%	1725 (610–2990)	−36%	−7%

and

Table 9. Sensitivity of mortality results to alternative RR.

Pollutant Metric	Health Outcome	2010 (min–max)	2030 WM (min–max)	2030 WM–2010 (%)	2030 WAM (min–max)	2030 WAM–2010 (%)	WAM–WM (%)
PM2.5 annual mean	Mortality all-natural causes	26,448 (17,574−44,367)	16,950 (11,296–28,269)	−36%	15,760 (10,508–26,264)	−40%	−7%
	Mortality cardiovascular diseases	13,633 (0–27,653)	8031 (0–16,074)	−41%	6857 (0–13,711)	−50%	−15%
	Mortality respiratory diseases	2876 (0–9751)	1687 (0–5495)	−41%	1439 (0–4677)	−50%	−15%
	Mortality lung cancer	7927 (4040–11,568)	4427 (2337–6266)	−44%	4095 (2168–5780	−48%	−7%
NO2 annual mean	Mortality all-natural causes	4216 (0–6351)	627 (0–941)	−85%	299 (0–447)	−93%	−52%

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₂ (all natural causes) and O₃ (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 WM scenario leads to significant reductions for all these pollutants: 37% (equal to 21,532 avoided cases), 85% (10,042) and 31% (841), for PM_2.5, NO₂ and O₃, 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₂ and 967 for O₃ considering the 2010 reference scenario. The relative reduction (close to 100%) of airborne NO₂ 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₂ 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 Figure 1, Figure 2 and Figure 3, for PM_2.5, NO_x and O₃ 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₂ and O₃, mortalities occur in a small fraction of the territory (Figure 2 and Figure 3). For NO₂, in 2010, events are located at the main Italian metropolitan areas, i.e., Naples, Rome and Milan, at specific urban areas in the Po Valley and along highway connections between Milan, Venice and Bologna. In 2030, events occur only in correspondence of the metropolitan area of Milan, which is one of the areas mostly interested in the pollutant reductions between 2010 and 2030 and between WM and WAM. Few events related to O₃ are recorded at the metropolitan areas of Milan, Venice, Rome and Naples.

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.

3.3.2. Mortality Rate

Mortality rates are reported as the average of non-zero values obtained on grid cells of the computational domain (

Table 10. National-level mortality rate over 10,000 inhabitants; national average (both sex); and variations between 2010, 2030 WM and 2030 WAM.

Pollutant Metric	Health Outcome	2010	2030 WM	2030 WM–2010 (%)	2030 WAM	2030 WAM–2010 (%)	WAM–WM (%)
PM2.5 annual mean	Mortality all-natural causes	7.25	4.73	−35%	4.43	−39%	−6%
	Mortality cardiovascular diseases	4.53	3.16	−30%	2.96	−35%	−7%
	Mortality respiratory diseases	0.71	0.50	−30%	0.47	−35%	−6%
	Mortality lung cancer	0.43	0.30	−30%	0.28	−34%	−7%
NO2 annual mean	Mortality all-natural causes	3.34	2.47	−26%	2.05	−39%	−17%
O3 April–September mean of MDA8	Mortality respiratory diseases	0.40	0.26	−36%	0.24	−40%	−7%

). Like for the previous section, Table 10 shows the absolute and relative mortality rate variations between 2010 and 2030 and between WM and WAM.

Considerations similar to those on absolute mortality can be made for the normalized figures, with the highest impact deriving from PM_2.5 (about 7.25 deaths each 10,000 inhabitants), followed by NO₂ (3.34) and O₃ (0.40). Relative reductions between 2010 and 2030 and between WM and WAM are similar to those of Table 8, except for NO₂. For this pollutant, in 2030 WM (WAM), the rate is reduced by 26% (39%) with respect to 2010, compared to 85% (93%) of absolute mortality.

The maps of mortality rates (PM_2.5 all-natural causes–PM_2.5 concentrations and variations in the three simulations are reported in Supplementary Materials Figure S2) for the 2010 scenario and for the expected variations in 2030 WM vs. 2010 and 2030 WAM vs. 2010 scenarios are shown in 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₂ (Supplementary Materials Figure S3), and O₃ (the maps are not reported since the values of standardized mortality are always below 1) and are not here furthermore discussed.

3.4. 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. 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.

Regions	Benefits 2030 WM–2010(mln €)	Additional Benefits WAM–WM (mln €)	2030 WM–2010 % of GDP2010	WAM–WM % of GDP2010
Abruzzo	350 (205–506)	44 (26–64)	1.14%	0.14%
Basilicata	83 (49–119)	7 (4–10)	0.74%	0.06%
Calabria	242 (142–349)	22 (13–32)	0.73%	0.07%
Campania	2122 (1227–3079)	157 (92–227)	2.06%	0.15%
Emilia–Romagna	2617 (1510–3814)	311 (183–449)	1.90%	0.23%
Friuli-Venezia Giulia	510 (297–741)	60 (35–86)	1.46%	0.17%
Lazio	4163 (2388–6056)	310 (180–448)	2.22%	0.17%
Liguria	638 (371–924)	73 (43–105)	1.38%	0.16%
Lombardia	11,925 (6773–17,575)	1697 (973–2472)	3.41%	0.49%
Marche	465 (272–674)	53 (31–77)	1.18%	0.14%
Molise	59 (35–85)	6 (3–8)	0.88%	0.09%
Piemonte	2382 (1377–3461)	364 (212–524)	1.91%	0.29%
Puglia	667 (389–964)	50 (29–72)	0.95%	0.07%
Sardegna	250 (147–361)	16 (10–24)	0.76%	0.05%
Sicilia	663 (386–960)	36 (20–53)	0.75%	0.04%
Toscana	1210 (705–1753)	149 (88–215)	1.15%	0.14%
Trentino-Alto Adige	187 (110–269)	21 (13–30)	0.98%	0.11%
Umbria	250 (146–362)	32 (19–46)	1.13%	0.15%
Valle D’Aosta	21 (12–30)	3 (2–4)	0.44%	0.06%
Veneto	2944 (1705–4293)	354 (208–512)	2.06%	0.25%
ITALY	29,691 (17,088–43,320)	3415 (1981–4947)	1.84%	0.21%

). 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₂ (16%) and O₃ (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₂-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₃ 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₂ (31% and 25%) and O₃ (3% in both scenarios). The PM_2.5 and NO₂ contributions to total benefits of Lombardia and Lazio are almost the same, indicating a large decrease in the attributable deaths to NO₂, as already underlined in Section 3.3.1.

4. 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₂ and O₃, respectively, in 2010. With respect to the present study, where the number of events is substantially the same for NO₂ and O₃, 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³ in VIIAS and 0 µg/m³ 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₁₀ 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₁₀ vs. PM_2.5), counterfactual concentrations used (20 vs. 0 µg/m³), 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 (Figure 1, Figure 2 and Figure 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₂. The absolute number of avoided events from PM_2.5 is more than the double of that from NO₂, 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₂, thus allowing for higher potential in reducing mortality for particulate exposure. This, in turn, is reflected in the high relative reduction of NO₂ mortality: The planned policies reach already almost the maximum feasible effect considering the threshold level of 20 µg/m³ 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₃, for NO₂, 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₂ 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₃, 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 country-specific 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.

5. 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.