Non-linear response of PM2.5 to changes in NOx and NH3 emissions in the Po basin (Italy): consequences for air quality plans

Air pollution is one of the main causes of damages to human health in Europe, with an estimate of about 380 000 premature deaths per year in the EU28, as the result of exposure to fine particulate matter (PM2.5) only. In this work, we focus on one specific region in Europe, the Po basin, a region where chemical regimes are the most complex, showing important non-linear processes, especially those related to interactions between NOx and NH3. We analyse the sensitivity of PM2.5 concentration to NOx and NH3 emissions by means of a set of EMEP model simulations performed with different levels of emission reductions, from 25 % up to a total switch-off of those emissions. Both single and combined precursor reduction scenarios are applied to determine the most efficient emission reduction strategies and quantify the interactions between NOx and NH3 emission reductions. The results confirmed the peculiarity of secondary PM2.5 formation in the Po basin, characterised by contrasting chemical regimes within distances of a few (hundred) kilometres, as well as non-linear responses to emission reductions during wintertime. One of the striking results is the slight increase in the PM2.5 concentration levels when NOx emission reductions are applied in NOx-rich areas, such as the surroundings of Bergamo. The increased oxidative capacity of the atmosphere is the cause of the increase in PM2.5 induced by a reduction in NOx emission. This process could have contributed to the absence of a significant PM2.5 concentration decrease during the COVID-19 lockdowns in many European cities. It is important to account for this process when designing air quality plans, since it could well lead to transitionary increases in PM2.5 at some locations in winter as NOx emission reduction measures are gradually implemented. While PM2.5 chemical regimes, determined by the relative importance of the NOx vs. NH3 responses to emission reductions, show large variations seasonally and spatially, they are not very sensitive to moderate (up to 50 %– 60 %) emission reductions. Beyond 25 % emission reduction strength, responses of PM2.5 concentrations to NOx emission reductions become non-linear in certain areas of the Po basin mainly during wintertime.


9310
P.  to changes in NO x and NH 3 emissions fraction is often dominating the total concentration of particulate matter in urban areas as shown by e.g. Beekmann et al. (2015) for the Greater Paris region and by De Meij et al. (2006 or Larsen et al. (2012) in northern Italy, hence the importance of understanding the complex chemical processes that lead to its formation. In particular, it is key to identify the precursors involved in these reactions in order to target the right sectors of activity in air quality plans to effectively reduce pollution levels. According to the EDGAR estimates for 2015 for Italy, about 90 % of the NH 3 is directly emitted in the atmosphere by the agriculture sector, while SO x precursors are predominantly released by the energy production and use (industrial) sectors (EDGAR, 2020). For NO x , emissions are spread among various sectors, with transport (50 %), industry (40 %) and agriculture (4 %) being the main ones. The gaseous precursors of secondary organic aerosols (SOA) include a vast range of NMVOCs among which are biogenic terpenes and anthropogenic aromatics. The main sources of aromatics in Italy were in 2012 transport (58 %), use of fuels and solvents (32 %), and domestic heating (15 %) (EDGAR, 2020).
Regarding SIA, early works used box models with thermodynamic schemes to address the sensitivity of ammonium nitrate and sulfate concentrations to gaseous NH 3 , NO x and SO 2 emissions (Watson et al., 1994;Blanchard and Hidy, 2003;Pozzer et al., 2017;Guo et al., 2018;Nenes et al., 2020). These models were later on integrated into chemical transport models (CTMs), in particular to address the benefit of additional NH 3 emission reductions in addition to already ongoing SO 2 and NO x emission reductions. For North America, Makar et al. (2009) simulated with a regional CTM that a 30 % reduction of ammonia emissions would lead to about 1 µg m −3 reduction in PM 2.5 . For Europe, Bessagnet et al. (2014) simulated the effects of a 30 % NH 3 emission reduction in addition to those foreseen by the Gothenburg Protocol for 2030 and found that the G ratio defined as the ratio between free ammonia and total nitrate (Ansari and Pandis, 1998) was a good predictor for the efficiency of NH 3 reductions on SIA concentrations. These sensitivities to emission reductions are often governed by complex chemical mechanisms. A well-known phenomenon is the release of free ammonia as a result of decreased SO 2 emissions and sulfate formation, which allows for the formation of additional particulate nitrate, as described for example by Blanchard and Hidy (2003) and Shah et al. (2018) for wintertime PM 2.5 over the eastern United States. For eastern China, Fu et al. (2017) and Lachatre et al. (2019) showed both from modelling and satellite observations that such processes lead to strongly increased ammonia tropospheric columns. Finally, several works compared CTM simulations to specific observations. For instance, Pay et al. (2012) showed that the G ratio was generally underestimated over Europe, inducing that the Caliope model they used could probably overestimate the efficiency of NH 3 emission reductions. Petetin et al. (2016) came to a similar conclusion comparing CHIMERE CTM simulations to observations in the Paris region. They ascribed this underestimation to missing NH 3 emissions especially during warmer periods.
The formation of SOA results from even more complex reactions involving photo-chemical oxidation (as for SIA), nitration, fragmentation and oligomerisation of gaseous precursors or secondary products (Kroll and Seinfeld, 2008;Shrivastava et al., 2017). Models generally use simplified parameterisations to calculate the SOA formation yield from various classes of parent VOCs (Tsigaridis et al., 2014), and comparison with measurements often shows that SOA sources are still missing in models (Huang et al., 2020;Tsimpidi, 2016).
In this work, we focus on one specific region in Europe, the Po basin. In a companion paper (Clappier et al., 2021) that analyses PM secondary formation chemical regimes across Europe, the Po basin is clearly identified as a peculiar area where the chemical regime distributions are the most complex, showing non-linear processes (Thunis et al., 2013(Thunis et al., , 2015Carnevale, 2020;, especially those related to interactions between NO x and NH 3 . The Po basin is also one of the pollution hot spots in Europe where the number of days above the limit values prescribed by the European Ambient Air Quality Directives (AAQD) for PM 10 is yet largely exceeded (EEA, 2020). This situation results from the high emission density in this region and also from the geographical setting of the area, in the border of the Alps and Apennines mountain ranges that lead to very weak winds in the area, favouring the accumulation of atmospheric pollutants.
We focus the present analysis on the NH 3 -NO x chemical processes and describe their spatial and seasonal variability, which could help to design more effective mitigation strategies. We start by describing the modelling set-up and detail the series of simulations required to perform our analysis. Section 3 provides a brief overview of the modelled basecase concentrations. In Sect. 4, we analyse the sensitivity of PM 2.5 concentrations to NH 3 and NO x emissions. Section 5 provides an analysis of the non-linearity in PM 2.5 response to these emissions, while in Sect. 6 we discuss the implications of these results in terms of mitigation measures and design of air quality plans. Conclusions are finally proposed.

Modelling set-up
The modelling study is performed with the EMEP air quality model, version rv4_17 . The emission input consists of gridded annual national emissions (SO 2 , NO, NO 2 , NH 3 , NMVOC, CO and primary PM 2.5 ) at 0.1 × 0.1 • resolution, based on data reported every year by parties to the Convention on Long-range Transboundary Air Pollution (CLRTAP). These emissions are provided for 10 anthropogenic source sectors classified by SNAP (Selected Nomenclature for Air Pollution) codes (EMEP, 2003). Meteorological input data are based on forecasts from the Integrated Forecasting System (IFS), a global operational forecasting model from the European Centre for Medium-Range Weather Forecasts (ECMWF). Meteorological fields are retrieved at a 0.1 × 0.1 • longitude-latitude resolution and are interpolated to the 50 × 50 km 2 polar stereographic grid projection (EMEP, 2011).
The gas-phase chemistry is based on the evolution of the so-called "EMEP scheme" as described in  and references therein. The chemical scheme couples the sulfur and nitrogen chemistry to the photochemistry using about 140 reactions between 70 species (Andersson-Sköld and Simpson, 1999;Simpson et al., 2012). In the EMEP Status Report 1/2004 (Fagerli et al., 2004) the reactions that cover acidification, eutrophication and ammonium chemistry are described. The aqueous-phase chemistry describes the formation of sulfate in clouds via SO 2 oxidation by ozone and H 2 O 2 and catalysed by metal ions. An important pathway of particulate nitrate formation is through the hydrolysis of N 2 O 5 on wet aerosol surfaces that converts NO x into HNO 3 . More information on the chemical equations is given in Simpson et al. (2012), Sect. 7.
For inorganic aerosols, EMEP uses the MARS equilibrium module to calculate the partitioning between gas and finemode aerosol phase in the system of SO 2− 4 , HNO 3 , NO − 3 , NH 3 and NH + 4 (Binkowski and Shankar, 1995). Aerosol water is calculated to account for particle water within the PM 2.5 mass, which depends on the mass of soluble PM fraction and on the type of salt mixture in particles. Sea salt (sodium chloride) and dust components are not accounted for by MARS, which might lead to PM underestimations close to coastal sites and where the dust contribution is important. More information on the gas and aerosol partitioning is given in Simpson et al. (2012), Sect. 7.6.
Regarding secondary organic aerosols (SOA), the Em-Chem09soa scheme is used, which is a simplified version of the so-called volatility basis set (VBS) approach (Robinson et al., 2007;Donahue et al., 2009). The VBS mechanism is discussed in detail in Bergström et al. (2012). The main difference between the VBS schemes and EmChem09soa is that all primary organic aerosol (POA) emissions are treated as non-volatile in EmChem09soa. This is done to keep the emission totals of both PM and VOC components the same as in the official emission inventories. The semi-volatile biogenic and anthropogenic SOA species are assumed to further oxidise (also known as ageing process) in the atmosphere by OH reactions. This will lead to a reduction in volatility for the SOA and thereby increased partitioning to the particle phase.
More information on SOA is given in Simpson et al. (2012), Sect. 7.7.
The modelling domain covers the entire Po basin ( Fig. 1) with a 0.1 by 0.1 • resolution (polar stereographic projection centred at 60 • N) and includes 20 vertical levels. The initial and background concentrations for ozone are based on Logan (1998) climatology, as described in Simpson et al. (2003). For the other species, background/initial conditions are set within the model using functions based on observations (Simpson et al., 2003;Fagerli et al., 2004). The simulations cover the entire meteorological year 2015. We will not discuss the validation of the base-case simulation, as this is available in other publications (e.g. Simpson et al., 2012) and regular status reports by EMEP (https://emep.int/ mscw/mscw_publications.html, last access: 1 June 2021).
In this work, we simulated a series of 24 scenarios in which NO x and NH 3 emissions were reduced independently or simultaneously by 25 %, 50 %, 75 % and 100 % from the base-case reference levels. Emission reductions were applied over the entire Po basin domain for a complete meteorological year (2015).

Spatial and temporal focus
Results are generally presented in terms of maps, but three locations within the domain were selected for a more detailed analysis. The locations are Bergamo (Be) in the northern part of the domain, Mantua (Ma; also known as Mantova) in the central eastern part of the Po basin and Bologna (Bo) in its southern part. As described in the following sections, we will see that these locations show very different behaviours in terms of response to emission changes.
We also aggregate results into two seasons: winter and summer, which cover the period from November to February and from April to September, respectively. These two seasons are characteristic of different chemical regimes as illustrated in the following sections. The process to define the temporal bounds of these two seasons is discussed in Appendix A. The two remaining months (March and October) represent transition periods and are not considered in our analysis.
As we only analyse processes involving inorganic gasphase precursors, our focus is on secondary inorganic PM 2.5 , although most of the results are expressed in terms of total PM 2.5 concentrations. The impact of NO x emission reductions on SOA concentration is only briefly discussed in Sect. 5.

Indicators
To describe the interactions between NH 3 and NO x emissions, we use the relationship proposed by Stein and Alpert (1993) and Thunis and Clappier (2014). This relation expresses the change in concentration resulting from a reduction of both precursors NO x and NH 3 simultaneously, as the sum of two single concentration changes and an interaction term, as follows: where C stands for the PM 2.5 concentration change (reference minus scenario) for a percentage α emission reduction (thus the term C/α is defined as positive for a concentration reduction, consecutive to an emission reduction) andĈ for the interaction term. We then scale each of these terms by the emission reduction (α) to generate potential impacts (P ) (Thunis and Clappier, 2014).
This division by the factor α is a means to extrapolate virtually the impact resulting from any percentage emission reduction to 100 %. Potential impacts facilitate the comparison of concentration changes obtained for different emission reduction levels. Indeed, equal potential impacts imply a linear relationship between emission reductions and concentration changes. For example, P α = P β ⇒ C α = α β C β for α and β, two emission reduction levels.
The overall potential impact is therefore the sum of two single potential impacts and one interaction term.
A relation between the potential impacts of combined emission reductions at two levels of intensity is obtained by writing Eq. (2) for two reduction levels α and β and subtracting the two equations. This leads to the following relation: or equivalently via Eq. (2), While the two first terms on the right-hand side of Eq. (4) represent the single potential impacts, the remaining righthand side terms quantify the magnitude of non-linearities.
NH 3 are the single non-linearities associated with NO x and NH 3 emissions, respectively, between levels α and β; andP β−α NO x NH 3 represents the incremental change in the NO x -NH 3 interaction between levels α and β.
Information about non-linearity is important to design air quality plans as it informs on the robustness of a given response, i.e. whether or not this response remains valid over a certain range and type of emission reductions. Because air quality models often provide responses for a limited set of scenarios that are then used as a basis to interpolate/extrapolate the responses to other emission reduction levels, robustness shall always be carefully assessed.
In the next section, we present the baseline results in terms of spatial and temporal variations.

Baseline concentrations of PM 2.5 and gaseous inorganic precursors
Before analysing the impact of emission changes on concentrations, it is worth having a look at the baseline concentration fields. In Fig. 1, the yearly averaged PM 2.5 concentration fields show a widespread pollution plume covering most of the area, with peak values extending in its central part. The maximum modelled yearly values reach 29 µgm −3 , which represents an average between maximum winter values (maximum of 59 µgm −3 ) and minimum summer values (17 µgm −3 ). The seasonal fields of PM 2.5 clearly show that high yearly average values mostly result from the winter season contributions when more stable atmospheric conditions lead to stagnant conditions, favouring the accumulation of particulate matter in the area (Pernigotti et al., 2014;Raffaelli et al., 2020). The increased emissions from the residential sector (heating, especially wood burning) foster this process (Ricciardelli et al., 2017;Hakimzadeh et al., 2020). Wintertime low temperatures also favour the partitioning of semi-volatile components (e.g. ammonium nitrate) towards the particulate phase. Overall, the relative contribution of secondary inorganic particles (SIA) ranges between 40 % and 50 %, regardless of the season, and is quite homogeneously distributed spatially over the entire area. Strategies targeting SIA have therefore the potential to abate about half of the total PM 2.5 concentration.
As mentioned earlier, the secondary inorganic fraction of PM 2.5 results from complex atmospheric processes that involve gaseous precursors (mainly SO 2 , NO x and NH 3 ), which can be summarised by the two following chemical pathways: where (g) means gas phase, (aq) means aqueous phase, and (p) means particulate matter, and the character →→ symbolises a chemical pathway that summarises a set of underlying reactions.
The second of these pathways (Reaction R2) is generally slower than the first one (Reaction R1), with the NO 2 oxidation specific time constant being typically some hours to a day and that of SO 2 being 1 to several days (Seinfeld and Pandis, 2006).
The spatial fields for the seasonal average concentrations of these precursors (Fig. 2) reflect their emission spatial patterns, resulting in a NO 2 -rich area that comprises Milan plus its northern districts (and to a lesser extent, Turin), while NH 3 is more abundant in the central part of the Po basin, east of Milan, where intensive agriculture practices take place. Finally, high SO 2 concentrations are collocated with the NO 2rich areas nearby Milan but with an additional zone around the harbour city of Genoa (also known as Genova; along the south coast), reflecting the more important SO 2 emissions from the shipping sector there. However, SO 2 concentrations are about 1 order of magnitude below those of NO 2 . Seasonal variations are well marked for NO 2 and SO 2 concentrations -not so much in terms of minimum and maximum values but rather in terms of spread with an extended high-concentration zone during wintertime. In contrast, NH 3 concentrations remain very similar in summer and winter, both in terms of values and spatial distribution.
In next sections, we simulate a series of emission reduction scenarios to analyse the response of PM 2.5 concentration to single and combined reductions of NH 3 and NO x .
4 Analysis of the SIA formation chemical regimes

Seasonal trends
From a strategic point of view, it is important to know whether NO x (mostly emitted by the transport, industry and residential sectors) or NH 3 (mostly emitted by agriculture) needs to be reduced in priority in order to reach effective results on particulate pollution mitigation. In Clappier et al. (2021), we have observed a great heterogeneity in SIA formation chemical regimes across the Po basin, with different regimes being present in limited geographical areas. Here we intend to look at these regimes in more detail.
To analyse in detail the chemical regimes in the Po basin, we compare the two single potential impacts P 25 % NO x , P 25 % NH 3 obtained for moderate emission reductions of 25 %. Figure 3 provides a spatial overview of the difference between these two potential impacts (P 25 % NO x − P 25 % NH 3 ). This indicator tells whether reductions of NO x or NH 3 will lead to the greatest PM 2.5 concentration abatement, i.e. if the regime is rather NO x -or NH 3 -sensitive, with positive and negative values, respectively. During summertime (Fig. 3 -left), the entire area is under weak NO x -sensitive conditions with a maximum intensity in its central part, between Bergamo and Mantua.
During wintertime (Fig. 3 -right), the situation is contrasted with a wide and intense NH 3 -sensitive area that appears around and south-eastwards of Bergamo. This area includes big cities like Milan. Other (not as marked) NH 3sensitive regime zones appear near coastal areas. Most strongly NO x -sensitive areas are located in the eastern parts of the domain, north of Bologna and Venice. NH 3 -sensitive regimes are generally collocated with the NO 2 -and SO 2rich areas (Fig. 2), whereas NO x -sensitive regimes coincide with NH 3 -rich areas. The cases of the three selected cities (Bergamo, Mantua and Bologna) representative of the NH 3sensitive, NO x -sensitive and neutral regimes, respectively, are further analysed below.
The chemical regimes deducted from the results of emission reduction scenarios can be compared with the maps of the G ratio (Fig. 4), defined by Ansari and Pandis (1998) as the ratio between free ammonia (NH 3 and NH + 4 ) and total nitrate (HNO 3 + NO − 3 ) after neutralisation of H 2 SO 4 . Values of the G ratio below 1 indicate a NH 3 -limited chemical regime, while values above 1 characterise a HNO 3 -limited chemical regime.
During summer, the G-ratio values well above unity indicate a HNO 3 limited chemical regime across the Po basin. This corresponds to a NO x -sensitive chemical regime in this region (Fig. 3). Moreover, the location of the G-ratio maximum between Bergamo and Mantua spatially coincides with the most pronounced NO x -sensitive regime and to a maximum of NH 3 concentrations of about 20 µg m −3 (Fig. 2). Indeed, NO x emission reductions lead to HNO 3 concentration reductions, which is the limiting factor in NH 4 NO 3 formation according to the G ratio. During winter, the G ratio still shows large values in the region south-east of Bergamo, but, contrary to summer, the chemical regime is clearly NH 3 -sensitive (Fig. 3). More generally, G-ratio values remain above unity over the whole Po basin, while both NO xand NH 3 -sensitive chemical regimes prevail in different areas. Therefore, the G ratio, related to the abundance of total free ammonia and total nitrate, provides information which differs from that obtained by determining the distribution of the NH 3 -and NO x -emission-sensitive chemical regimes. These differences illustrate the impossibility to directly use the G ratio for air quality management, an interesting result in itself. We will further discuss this interesting behaviour later when addressing non-linearity in Sect. 4.3.

Impact of the emission reduction strength
In this section, we repeat the analysis of Sect. 4.1 for yearly average concentrations but looking at the step changes of regimes as we progressively reduce emissions from the basecase situation. (Fig. 5). Chemical regimes are well in place for a 25 % level reduction (top left) and are only slightly perturbed from 25 % to 50 % with a reinforcement of the NO x  limited regime (top right). Despite this slight change in intensity, the regimes therefore keep the same spatial patterns. From 50 % onward, chemical regimes tend to attenuate and reverse themselves from 75 % to 100 % (bottom right).
In other words, locations that are NH 3 -sensitive for the first steps emission reductions will become NO x -sensitive for the last steps emission reductions, and vice versa.

A summarised overview: the PM 2.5 isopleths
Like isopleth plots that show the variations in the O 3 concentrations as a function of NO x and VOC concentrations (Dodge, 1977), similar plots can be created for PM 2.5 concentrations as a function of NO x and NH 3 emissions. Simulation results have indeed often been presented as 2D isopleths of PM 2.5 or nitrate as a function of precursor emissions, which allows showing in a comprehensive manner their sensitivity and also in a qualitative manner non-linear effects (for example Watson et al., 1994, over the USA;and Xing et al., 2018, over the Beijing-Tianjin-Hebei region in China). Figure 6 shows the PM 2.5 isopleths obtained through an interpolation among the 25 simulation concentration values (these 25 simulations correspond to the white square symbols in each isopleth diagram) at the three locations Be, Ma and Bo previously defined. The x and y axes represent the strengths of the NH 3 and NO x emission sources, respectively. With this type of graphical representation, it is possible to visualise the response of PM 2.5 to a NO x emission change by moving vertically, the reaction to a NH 3 emission change by moving horizontally or the reaction to a combined NO x -NH 3 emission change by moving diagonally. The larger the number of isopleths we cross on the path (high gradient), the larger the expected impact from an emission reduction will be. A simple theoretical model to generate and interpret these isopleths is proposed in Appendix B.
From the analysis of the isopleths, we note the following points. -In general, the isopleths show a regular pattern with a progressive decrease in PM 2.5 concentration when either the NO x or NH 3 emissions are reduced, with the only exception being Bergamo during wintertime, where NO x reductions up to 70 % lead to a small increase in PM 2.5 (Fig. 6 -top left), whatever the reduction in NH 3 emission is. We discuss this particular feature later in this section.
-The diagram areas can be divided into two zones separated by a ridge (dashed line in Fig. 6). Above the ridge line, PM 2.5 is more sensitive to NH 3 , while below the ridge line, it is more sensitive to NO x emission reductions. The orientation of the ridge (tending to vertical or horizontal) informs on the type of chemical regime (NO x -or NH 3 -sensitive, respectively).
-The efficiency of a NO x vs. NH 3 emission reduction varies across locations. We can compare the efficiency of a given reduction by looking at the horizontal (for NO x ) and vertical (for NH 3 ) gradients. To support this comparison, we included in each diagram dashed oblique lines that connect similar PM 2.5 concentration values for single NO x and NH 3 emission reductions. The more vertical these lines are, the larger the NH 3 abatement impact compared to the NO x abatement impact is. Conversely, the more horizontal they are, the larger the NO x abatement impact compared to the NH 3 abatement impact is. For moderate emission reductions (up to 50 %, top-right corner), different behaviours are observed: while in Bergamo PM 2.5 is more sensitive to NH 3 reductions, it is more sensitive to NO x reductions in Mantua, and it is equally sensitive to both precursors in Bologna. This corresponds to the spatial patterns of NO x -and NH 3 -sensitive regimes depicted in Fig. 3.
-Winter-and summertime isopleths show completely different patterns in Bergamo, whereas at the two other locations, they remain similar.
-At Mantua where moderate NO x reductions (e.g. 50 %) are the most efficient among the three sites, NH 3 emission reductions are more efficient than NO x emission reductions for larger additional reductions (going for example from 75 % to 100 %). At Bergamo, NH 3 reductions are the most effective for moderate reductions, whereas NO x reductions become more effective for larger reductions, as seen by the isopleths spacing. This confirms the finding of reversed chemical regimes for larger additional emission reductions detailed in the previous section and illustrated in Fig. 5.
The special pattern of Bergamo's PM 2.5 isopleths during wintertime needs some additional discussion. The increase in the inorganic fraction of PM 2.5 as a response to NO x reductions during wintertime has already been noted by several authors (e.g. Le et al., 2020;Sheng et al., 2018). It has been related to an increase in the oxidising capacity of the atmosphere and in particular to increased ozone levels. This is due to the prevailing titration of O 3 by NO in wintertime high-NO x conditions and in the absence of photochemical ozone production due to reduced solar radiation (Kleinman et al., 1991).
The impact of NO x emission reductions on the concentration of various pollutants in Bergamo during wintertime is illustrated in Fig. 7. As expected, a 50 % reduction in NO x emissions leads to a decrease in NO 2 concentration (from 47 to 28 µg m −3 , i.e. a factor of 1.7). In contrast, O 3 concentration increases from 8 to 16 µg m −3 , roughly a factor of 2. These compensating changes result in a small increase in NO 3 radical production (Reaction R4), which is the initial step of the major pathway of wintertime HNO 3 and particulate nitrate formation (Kenagy et al., 2018).
In this pathway, the NO 3 radical formation is followed by combination with NO 2 to form N 2 O 5 , a reversible process, and heterogeneous HNO 3 formation on wet particle surfaces. NO 3 (g) + NO 2 (g) ↔ N 2 O 5 (g) (R5) Figure 6. PM 2.5 isopleths during winter (a, b, c) and summer (d, e, f) at the three locations of interest (see maps). PM 2.5 concentrations are expressed in µg m −3 as a function of the intensity of the NO x (y axis) and NH 3 (x axis) emissions, respectively. The overlaid dashed oblique lines on each diagram connect similar PM 2.5 concentration values for single NO x and NH 3 reductions. The more vertical these lines are, the larger the NH 3 abatement impact compared to the NO x abatement impact is; the more horizontal they are, the larger the NO x abatement impact compared to the NH 3 abatement impact is. The dashed line delineates the ridge between the NH 3 -and NO x -sensitive areas.
The NO 3 radical has three major rapid sinks: reaction with NO, photolysis, and reaction with NMVOCs, especially terpenes.
Reactions (R3) to (R10) induce additional dependence of HNO 3 formation on NO x species on top of Reaction (R1), but which partly cancel out, as they are both involved in formation and sink processes.
SOA is formed through a series of chemical reactions of gaseous precursors -mainly volatile, intermediate volatility or semi-volatile organic compounds (VOCs) with the oxidants O 3 , OH and nitrate radical (NO 3 ) (Li et al., 2011). Putting all the arguments together, it follows that wintertime ammonium nitrate formation over Bergamo is most probably controlled by NO 3 radical formation (Reaction R4). The fact that this behaviour is observed in Bergamo and not in Mantua or Bologna is due to the much larger NO 2 levels in the Bergamo-Milan area (above 50 µg m −3 during winter, Fig. 2). Such large NO 2 levels are also simulated locally over the Turin area and also lead to a slightly NH 3 -sensitive regime there despite a G ratio well above unity. Beyond 50 % NO x reduction, NH 4 NO 3 formation decreases because NO 2 decreases more rapidly than ozone increases up to its maximum (at 75 % NO x emission reduction; see Fig. 7).
The negative response of NH 4 NO 3 to NO x emission reductions in Be during wintertime explains the apparent discrepancies with the analysis of G ratio, which indicates that NH 4 NO 3 is strongly HNO 3 limited. Simply, the HNO 3limited chemical regime cannot be extrapolated to sensitivity to NO x emissions in the case of the above shown negative response. Total nitrate is less abundant than free ammonia (defined as NH 3 (g) + NH + 4 (p) − 2SO 2− 4 (p)), but NO x emission reductions up to about 50 % do not reduce its concentration, and NH 3 emission reductions are thus more efficient. In this respect, the G ratio cannot provide information about negative responses.
At Bergamo during winter, the increase in PM 2.5 (+1.8 µg m −3 ) arising from a 50 % reduction in NO x emission also results from an increase in sulfate (+0.3 µg m −3 ) and in SOA (+0.6 µgm −3 ) concentrations. Both sulfate and SOA concentrations are closely related to O 3 concentrations (Fig. 7). The sulfate increase is comparable in magnitude to the nitrate increase, even if sulfate levels are much smaller than nitrate ones. Figure 7 actually shows a strikingly similar response of sulfate, SOA and ozone to NO x emission reductions (given that SO 2 and NMVOC emissions are held constant). Indeed the prevailing wintertime aqueous production of H 2 SO 4 requires oxidants and in particular ozone (Le et al., 2020;Sheng et al., 2018). In addition, the formation of SOA in both the gas and particulate phases also requires oxidants (Vahedpour et al., 2011;Huang et al., 2020;Feng et al., 2016Feng et al., , 2019Li et al., 2011;Tsimpidi et al., 2010). Pinder et al. (2008) also note an oxidant limitation for SIA formation over the eastern United states for the 2000 to 2020 period, but in their simulations, it mainly affects sulfate that increases as a result of NO x emission reductions while nitrate decreases. This is due to a more important sulfate-to-nitrate ratio in the eastern United States than over the Po basin. Fu et al. (2020) derive from combined measurements and modelling that wintertime nitrate during haze events in the North China Plain (NCP) is nearly insensitive to 30 % NO x emission reductions, because increased ozone levels increase the NO x to HNO 3 conversion efficiency. Following these authors, this conversion also involves the homogeneous HNO 3 formation via the NO 2 + OH. This reaction also could play a role in the Po basin, in addition to the heterogeneous pathway. Also Leung et al. (2020) simulate that wintertime nitrate abatement in the NCP is buffered with respect to emission reductions by increased oxidant build-up but also by sulfate to nitrate conversion by liberating free NH 3 through sulfate concentration reduction, which can then enhance nitrate formation. Womack et al. (2019) find an oxidant limitation of nitrate formation over wintertime Utah (USA) and show that nitrate concentration diminishes when reducing VOC emissions. Clappier et al. (2021) highlighted the specificities of the Po basin area within Europe. They showed that non-linearities are present in this region. One peculiarity of the Po basin is a marked difference between the chemical regimes encountered within a confined area, which have implications on the linearity of PM 2.5 responses to emission changes. In this section, we analyse in more details these non-linearities.

Analysis of non-linearities
Information on the NO x -NH 3 interaction term at the reduction level α = 25 %,P 25 % NO x NH 3 [= P 25 % NO x NH 3 − (P 25 % NO x + P 25 % NH 3 )] (first non-linear term in Eq. 4), is provided in Fig. 8. At α = 25 %, the interaction term is negative (or null) across the entire modelling domain (most data points are below the 1:1 line) regardless of the chemical regimes and averages to approximately −10 %, as indicated by the linear fit slope (= 0.9). In relative terms, this interaction term is also almost constant regardless of the season (not shown).
This negativity can be explained by the fact that a reduction of only NO x implies a reduction of both NO − 3 and NH + 4 , and the same happens when reducing only NH 3 ; therefore a simultaneous reduction of both precursors is lower than the sum of the two. Single impacts would therefore lead to an overestimation (of about 10 %) in PM 2.5 reduction if added up to extrapolate linearly the impact of combined 25 % NO x and NH 3 emission reductions on yearly averaged PM 2.5 concentrations. This result is expected for what concerns particulate NH 4 NO 3 , as a consequence of the gas-particle equilibrium described in Reaction (R1), although non-linear rela-tionships between NO x emissions and HNO 3 concentrations also play a role. Qualitatively, this negative interaction is also highlighted by the hyperbolic shapes of the PM 2.5 isopleths determined for three different sites of the domain (Fig. 6). As discussed in Appendix B, linearity would result in isopleths parallel to the descending diagonal lines.
When emission reductions increase from 25 % to 50 %, three additional non-linear terms are generated (three last terms in Eq. 4). Figures 9 and 10 provide an overview of these non-linear terms during wintertime and summertime, respectively. The top-left panel of each figure represents their sum, i.e. the total non-linearity generated between 25 % and 50 % emission reduction. The right column shows the nonlinearities associated with NO x and NH 3 , while the bottomleft panel reports the non-linear interaction between the two precursors.
Overall, non-linearities are more important during wintertime than during summertime. This is true both in absolute and relative (not shown) terms. Non-linearities tend to be the largest in between areas that are characterised by wellmarked NH 3 -or NO x -sensitive regimes (indicated by the blue and red drawn contours). This can be explained by the fact that when one of the two components (NO x or NH 3 ) is in large excess (compared to the other one), reductions of this compound have then generally only little impact, implying that both single and combined reductions only involve one compound and are therefore similar.
During wintertime, the overall non-linearity ( Fig. 9 top left) is largely dominated by the single NO x -related nonlinearity (P 50 %-25 % NO x , Fig. 9 top right) -a singularity in Europe as the Po basin is the only area where this occurs to this extent (Clappier et al., 2021). In the region of Bergamo, the NO x non-linearity remains weak despite the peculiar PM 2.5 responses to NO x emission reductions (i.e. an increase in PM 2.5 concentrations for NO x emission reduction up to 50 %, Fig. 7). In this NH 3 -sensitive region, this behaviour can be explained by the strong oxidant limitation of HNO 3 formation outlined above (Fig. 7). It is worth mentioning that, although atypical, this behaviour is quasi-linear with PM 2.5 responses that remain proportional to the emission reduction strength (up to 50 %) but with a negative slope. Similarly to NO x , NH 3 single non-linearities (Fig. 9, bottom right) are positive but weaker, indicating slightly larger potential impacts for emission reductions in the range 25 %-50 % than in the 0 %-25 % range. Finally, the NO x -NH 3 non-linear interaction terms (Fig. 9, bottom left) are mostly negative, indicating that the interaction term for 50 % emission reductions is more negative than the corresponding term for 25 %, pointing out to a strengthening of the NO x -NH 3 non-linearity when more intense emission reductions are considered.
In relative terms, the overall wintertime non-linearity terms increase by about 30 % when emission reductions increase from 25 % and 50 % (Fig. 11 top, blue points). Note that these non-linearity terms reach larger values in some places of the modelling domain as highlighted by the data Figure 9. Wintertime maps of the overall non-linearity term (a) and of its components expressed as potential impacts between 25 % and 50 %: the single NO x non-linearity term (a, b), the single NH 3 non-linearity term (d) and the NO x -NH 3 interaction term (c). The three locations of interest, Bergamo, Mantua and Bologna, are indicated by their first two letters, while other cities are indicated with their first letter for convenience (Venice, Milan, Turin and Genoa). The hand-drawn contours roughly indicate the central position of the NH 3 (blue) and NO x (red) sensitive regime areas (see Fig. 3 -right). Figure 10. Same as Fig. 9 but for summertime. point dispersion. In a previous work, Thunis et al. (2015) quantified the non-linearity of model responses to emission reductions in three areas in Europe, among which is the Po Basin. One of their conclusions was that non-linearities remain relatively low for yearly averaged responses. Although the results presented here show important non-linearities, these occur mainly during wintertime and are limited to specific areas. It is also worth noting that these non-linear responses (for moderate emission reductions up to 50 %) only occur in the Po Basin (Clappier et al., 2021).
During summertime, the magnitude of non-linear terms is smaller than during wintertime. The overall non-linearity term (top left) is dominated by the NH 3 -NO x non-linear interaction term (P 50 %-25 % NO x NH 3 , Fig. 10, bottom left). In summer, the PM 2.5 concentration is NO x -sensitive in almost all the domain (Fig. 3), and emission reductions do not lead to shifts in the chemical regimes. In this case, the interaction terms become more important. The non-linear interaction terms are negative everywhere in the Po basin, implying again that the interaction terms at 50 % are more nega- Figure 11. Changes in the overall non-linearity terms from 25 % to 50 % (a), from 50 % to 75 % (b) and from 75 % to 100 % (c) in NO x and NH 3 . The overall non-linearity is visualised as the distance from the 1 : 1 diagonal, i.e. the difference between the overall potential impact at two levels of emission reduction, x and y axis. Each point represents one land grid cell within the domain for wintertime (blue) and summertime (red). The "fit" parameter indicates the slope of the regression line, while R2 and RMS provide information on the coefficient of determination and the root mean square error, respectively. tive than at 25 % emission reduction (reinforcement of the non-linearities when more intense emission reductions are considered). Most negative values appear in the western and northern part of domain. Figure 11 shows the increase in the overall non-linear terms for emission reduction steps starting from different points, from 25 % to 50 %, from 50 % to 75 % and finally from 75 % to 100 %. Regardless of the emission reduction step, summertime non-linearities remain small all over the domain, with the regression slope close to 0.9 and very limited data point dispersions (i.e. low RMSE). Wintertime non-linearities further increase significantly from 50 % to 75 % reduction levels (regression fit parameter close to 1.21) but tend to stabilise between 75 % and 100 % reduction levels (regression fit parameter close to 1.05). It is interesting to note that potential impacts in winter increase for all segments (all winter points are above the 1 : 1 line), indicating that the same percentage reduction (25 %) gains progressively more impact when more intense reductions are considered.

Discussion
When designing air quality plans, it is important to identify the key precursors on which to act in priority to hit a specific air quality target but also to understand the consequences of these choices for various seasons (temporal variations), locations (spatial variability), emission reduction levels (strength) and strategies (combined or single emission reductions). As the information to make this decision is generally incomplete, assessing the robustness of the available model responses is essential. From the results presented here, a few key points appear.
The seasonal and spatial variabilities in the response of PM 2.5 to the reduction of NO x and NH 3 emissions are extremely large, with different and sometimes opposite responses to emission changes. Yearly averages do not represent the appropriate time window to evaluate the impact of such emission reductions, and a focus on wintertime (November to February) seems to be the right option, especially because concentrations are larger during this period of the year.
The responses of PM 2.5 to emission reduction plans that cover the whole area (i.e. uniform emission reductions are applied everywhere in the domain) vary from location to location: opposite responses occur within a few hundred kilometres for some reduction levels. In the region of Bergamo, the PM 2.5 response to NO x emission reductions can be negative, meaning an increase in PM 2.5 when reducing the NO x emissions. It is important to combine NO x and NH 3 emission reductions in winter or to go for stronger emission reductions to make sure these unwanted effects are limited.
Despite quite important non-linearities, PM 2.5 responses to emission reductions are not chaotic. Indeed, regardless of the emission reduction level, the non-linear terms related to NH 3 emission reduction and to NO x -NH 3 interactions are quite uniform spatially. This is not the case of NO x emission reduction, for which care must be taken to ensure that the detailed response of PM 2.5 is captured.
Although they are location specific, PM 2.5 isopleths represent an interesting tool to assess the impact of different NO x and NH 3 emission reductions on PM 2.5 concentration. They indeed allow visualising in one single diagram the impact of any type of reductions on concentrations in a single grid cell (or set of grid cells). We must however remember that these isopleths derive from uniform emission reduction over the whole domain. Comparing sites where PM 2.5 responses to the same "domain-wide" policy are different, it appears challenging to define a single domain-wide policy efficiently reducing PM at all locations.

Conclusions
In this work, we analysed the sensitivity of PM 2.5 to NO x and NH 3 emissions by means of a set of EMEP simulations performed with different levels of emission reductions, from 25 % up to a total switch-off of those emissions. Both single and combined precursor reduction scenarios were applied to determine the most efficient emission reduction strategies and quantify the interactions between NO x and NH 3 emission reductions. Our results confirm the peculiarity of secondary inorganic PM 2.5 formation in the Po basin suggested by Clappier et al. (2021), characterised by contrasting chemical regimes within distances of a few (hundred) kilometres, as well as non-linear responses to emission reductions during wintertime. One of the striking results is the increase in the PM 2.5 concentration levels when NO x emission reductions are applied in NO x -rich areas, such as the surroundings of Bergamo. The isopleths in the graphs showing PM 2.5 , nitrate, sulfate, SOA and O 3 concentrations as a function of NH 3 and NO x emissions (Fig. 7) indicate that the increased oxidative capacity of the atmosphere is the cause of the increase in PM 2.5 induced by a reduction in NO x emission of up to 50 %. This process can have contributed to the absence of significant PM 2.5 concentration decrease during the COVID-19 lockdowns in many European cities (EEA, 2020;Putaud et al., 2021). It is important to account for this process when designing air quality plans, since it could well lead to transitionary increases in PM 2.5 at some locations in winter as NO x emission reduction measures are gradually implemented. At this type of location, mitigation measures that would be optimal in the long term might not be efficient in the short term.
Joint analyses of PM sensitivity to emissions and the G ratio can give a clue if a NH 3 -sensitive chemical regime is due to either a lack of NH 3 or to a non-linear and negative response of HNO 3 concentration to NO x emissions. In this latter case, the chemical regime is NH 3 -sensitive in terms of NH 3 and NO x emission reductions, but it can be HNO 3 limited in terms of the G ratio, as observed for the Bergamo-Milan region. Inversely, a G ratio greater than 1, indicating HNO 3 limitation to particulate nitrate formation, does not necessarily indicate a larger sensitivity to NO x than to NH 3 emissions. Thus, our results show the impossibility to directly use the G ratio for air quality management, an interesting result in itself. While PM 2.5 chemical regimes (determined by the relative importance of the NO x vs. NH 3 responses to emission reductions) show large variations seasonally and spatially, they are not very sensitive to moderate (up to 50 %) emission reductions. Beyond 25 % emission reduction strength, responses of PM 2.5 concentrations to NO x emission reductions become non-linear in certain areas of the Po basin mainly during wintertime.
Since sulfate concentrations are low in the Po basin, the impact of SO 2 emission reductions was not evaluated here. However, the simulations performed in Clappier et al. (2021) indicate that air quality improvement plans addressing SO 2 emissions may still lead to additional PM 2.5 decreases. Further works should also test if NMVOC emission would further affect the concentration of oxidants and subsequently of nitrate (and sulfate) during winter. This depends on the fraction of ozone formed photochemically in the Po basin, compared to the one transported from outside by advection or entrainment.
Finally, it would be important to compare the results obtained in this work from the EMEPrv4_17 model with similar results obtained from other models. With its complex setting, the Po basin represents a good test case for such intercomparisons.

9324
P.  to changes in NO x and NH 3 emissions Appendix B: A theoretical example for the isopleths To facilitate the interpretation of the isopleths diagrams, we use a simple theoretical example that mimics the complex reactions process schematised by Reactions (R1) and (R2) above. Our simplified process is described by the following relation: C c (x, y) = min E A [x] , E B y , where C c is the concentration of a compound "c" that is given by the minimum between two emitted species A and B. The concentration depends on the strengths of these emissions, specified by the parameters x and y. For each emission strength (x or y), the two emission species are proportional to their full-scale value (100 %): E A (x) = xE A [100] and E B (y) = yE B [100], respectively. If we choose E A (100) E B (100), we create a B-sensitive environment (Fig. B1 -left column) and inversely (Fig. 13 -right column). If we select mixed situations, representing for example an average of days, during which we alternate between A-and B-sensitive regimes, we obtain the two bottom isopleth diagrams that represent cases where a larger number of A-sensitive (right) or B-sensitive (left) events are recorded. Although extremely simple, these diagrams illustrate properties that are observed on the real test cases.
Let us take the example of an A-sensitive regime. Similar observations can be made in the case of a B-sensitive regime. We note the following points.
-The diagram area can roughly be divided into two zones separated by a ridge: a top-left triangle where the sensitivity to emission reductions of species "B" dominates and a bottom-right triangle where the sensitivity to A dominates. The slope of the ridge (larger or less than 1) informs on the type of regime.
-In the case of a single A-sensitive day (top right) with E B (100) = 2E A (100), the concentration of compound "c" remains unchanged for emission reductions of B up to 50 %, while its concentration react in a linear way to emission changes in A from 0 % to 100 %. Between the base case and a reduction level of 50 %, the A gradient is therefore larger than the B gradient. This implies that the B gradient is larger than the A gradient between the 50 % and 100 % reduction levels because we know that for a 100 % reduction of A or B, the concentration must be zero. In our simple example, the gradient of B is zero from 0 % to 50 % but is twice as large as A between 50 % and 100 %.
-While the combination of several events (e.g. days) characterised by different regimes leads to smoother isopleths (bottom), the same characteristics can be noted. In particular, the inclination (tending to the horizontal or vertical) provides information on the type of chemical regime. Code availability. The Open Source EMEP model and all necessary input data to run the model are available through this link: https://github.com/metno/emep-ctm (Norwegian Meteorological Institute, 2021).
Data availability. Data are in the process of being transferred to a public repository. In the meantime they are available upon request to the authors.
Author contributions. PT and AC designed the methodology, analysed the results and drafted a first version of the paper. MB and JPP contributed to the analysis and interpretation of the results. ADM and JM performed the simulations. CC contributed to the data treatment and visualization aspects. All co-authors critically reviewed the paper.
Competing interests. The authors declare that they have no conflict of interest.
Review statement. This paper was edited by Astrid Kiendler-Scharr and reviewed by two anonymous referees.