Assessment of Air Quality Impacts from Wood-Fired Pizza Ovens: Comparison with Major Urban Emission Sources

1. Introduction

Urban air quality remains a critical environmental and public health issue, particularly in densely populated metropolitan areas where a complex mixture of emission sources contributes to atmospheric pollution [1,2]. While major ones such as traffic and residential heating have been extensively investigated [3], comparatively less attention has been devoted to small-scale commercial activities that may collectively contribute locally to air pollution. Among these, wood-fired ovens used in pizzerias represent a potentially relevant source in urban centers [4], where they are both highly concentrated and largely unregulated in terms of emissions.

Previous experimental studies have demonstrated that wood-fired ovens contribute to the emissions of several major air pollutants, especially fine particulate matter (PM_2.5) [5], polycyclic aromatic hydrocarbons (PAHs) [4], and black carbon [6]. In addition, these appliances emit well-known tracers of biomass combustion, such as benzo(a)pyrene (BaP) [4] and levoglucosan (LG) [7], which are widely used to identify and apportion biomass combustion sources in atmospheric studies [8,9]. The emission factors outlined in Bergomi et al. [4] have been applied in local emission inventories, which estimate contributions in urban areas of around 5% for PM₁₀ and more than 30% for selected PAHs [10]. However, the integration of such emission factors into urban-scale air quality assessments remains limited, mainly due to the lack of high-resolution spatial data on the distribution of these activities.

Recent advances in data science and geospatial analysis provide new opportunities to overcome these limitations by enabling the identification and mapping of sources, such as wood-fired ovens used in pizzerias, which are typically difficult to locate due to the absence of dedicated registries or official databases. In this context, the present study leverages data-driven and artificial intelligence techniques to reconstruct a geolocation dataset of facilities employing wood-fired ovens across the Municipality of Milan. This information is then combined with experimentally derived emission factors for wood-fired ovens [4] to perform a detailed assessment of their impact on urban air quality, including the dispersion of biomass-burning tracers, which may improve source attribution in urban environments.

Specifically, a Lagrangian puff dispersion modelling approach was adopted using the CALPUFF modelling system, driven by meteorological data obtained from the WRF (Weather Research and Forecasting) model. While CALPUFF is commonly used for source-specific and multi-source air quality assessments, applications involving a very large number of distributed emitters often become computationally demanding [11]. To overcome this limitation, the present work employs a hybrid modelling framework in which a reference concentration field generated by CALPUFF is linearly combined according to the location and emission strengths of each source. This approach exploits the linearity of dispersion processes for primary pollutants and enables the efficient estimation of cumulative concentration fields from many emitters while preserving the atmospheric transport and dispersion patterns reproduced by the original model.

Similar source-superposition principles and reduced-complexity approaches have been exploited in previous studies to support urban air-quality assessment and scenario analysis while reducing computational requirements [11,12]. However, previous reduced-complexity applications generally relied on urban background or city-increment estimates derived from measurements at air quality monitoring stations. While effective for improving the representation of urban concentration gradients, such approaches do not explicitly account for the combined effects of the spatiotemporal variability of emissions and meteorological conditions on pollutant accumulation and transport. In contrast, the framework proposed in this study performs continuous simulations throughout the entire year and dynamically represents the contribution of multiple source categories without relying on monitoring-derived background concentrations.

Moreover, the new EU Ambient Air Quality Directive (EU 2024/2881) [13] significantly strengthens the role of air quality modelling and points towards the use of high-resolution urban assessments, including microscale applications in areas characterized by pollution hotspots and strong spatial concentration gradients [13]. The Directive places modelling at the core of air quality assessment, requiring Member States to combine monitoring, indicative measurements and modelling tools in a more integrated way. It also emphasizes the need for assessments that can adequately capture population exposure, particularly in locations where pollutant concentrations may vary significantly over short distances.

The framework developed in this study enabled the simulation of pollutant dispersion at a high spatial resolution, with a specific focus on the pollutants of greatest concern in the investigated area: PM₁₀, PM_2.5, and NO_x. The modelling approach also provided a basis for interpreting the spatial variability of biomass-burning tracers. This aspect is particularly relevant in urban contexts where the contribution of residential biomass heating is relatively low compared with rural areas, yet measurable concentrations of both LG and BaP are still observed [14,15]. The model outputs were then compared with measured air quality data from the regional monitoring network of ARPA Lombardia, enabling an assessment of the relative contribution of wood-fired ovens to observed pollutant concentrations.

The main aim of this study is therefore to quantify the potential impact of wood-fired pizzerias on urban air quality in Milan, combining innovative data acquisition methods, experimental emission characterization, and advanced atmospheric dispersion modelling. By addressing a largely overlooked emission source, this work contributes to a more comprehensive understanding of urban pollution dynamics and may support future policies and regulatory strategies targeting small-scale combustion sources.

2. Materials and Methods

The area of interest of this study is the Municipality of Milan (Figure 1).

As one of the largest metropolitan areas in Italy, Milan faces significant air quality challenges associated with dense urbanization, intensive mobility, and high energy consumption.

The INEMAR (INventario EMissioni ARia) emission inventory is the official regional system used in Lombardy to estimate and monitor atmospheric emissions from both anthropogenic and natural sources. It provides a detailed, spatially resolved database of air pollutant emissions, supporting air quality assessment, environmental planning, and policy development. The INEMAR inventory therefore provides a crucial scientific basis for developing evidence-based strategies aimed at reducing emissions, improving public health, and achieving sustainability and decarbonization goals [10].

For the city of Milan, the inventory quantifies emissions from a wide range of sources, including road transport, heating systems, industrial activities, energy production, waste management, agriculture, and commercial activities. Among the latter, pizzerias represent a specific emission source due to the combustion and cooking processes occurring in pizza ovens, particularly wood-fired ovens. Indeed, according to the latest release of INEMAR, pizzerias with wood-fired ovens contribute to the total emissions of several air pollutants (Figure 2).

Although the contribution to the total emissions of Lombardy can be less relevant, the large number of pizzerias operating within the urban areas may result in a non-negligible contribution to local air pollutant concentrations, especially in densely populated districts.

2.1. Emission Factors from Wood-Fired Pizza Ovens

The air pollutant mass fluxes used to quantify emissions from wood-fired pizza ovens (Table 1) were derived from emission factors reported in [4]. These values were selected as reference emission factors because they currently constitute the only experimentally based quantification of emissions from wood-fired pizza ovens available in the scientific literature.

The dilution tunnel method was used to sample emissions from the wood-fired ovens. NO_x mass fluxes were calculated from measurements performed on hot flue gases, whereas particulate matter mass fluxes were derived from diluted samples collected from diluted flue gases. This approach accounts for the condensation of semi-volatile compounds, which contribute significantly to total particulate emissions of wood-fired ovens and should therefore be considered in dispersion assessments.

Emissions of the chemical tracers, benzo[a]pyrene (BaP) and levoglucosan (LG), were estimated as a fraction of the total PM₁₀ emissions, as these compounds are not explicitly included among the species available for dispersion simulations in the CALPUFF model. Consequently, only their particulate-bound fractions were considered in modelling simulations. These fractions were derived from contributions stemming from the original experimental study, which reported that BaP and LG account for 0.0111% and 3.78% of particulate matter emissions, respectively [7].

2.2. Number and Location of Wood-Fired Pizza Ovens in Milan

To estimate emissions generated by wood-fired pizza ovens, the first step was to construct an inventory of the emission sources within the study area (Metropolitan City of Milan). Since no official registry exists identifying food service activities equipped with wood-fired ovens, a dedicated data-driven methodology was developed to identify these establishments based on official business records and online data sources accessed either openly or under commercial licensing agreements. The approach combines administrative data, online business directories (i.e., structured listings of commercial activities such as OpenStreetMap [16] and The Data Appeal [17]), user-generated content (i.e., unstructured contributions such as reviews, photos, and other media shared by users on Google Places [18]), and artificial intelligence techniques within a multi-stage identification pipeline aimed at constructing a reliable inventory of potential emission sources. Access to The Data Appeal and the Google Places API was made possible through contractual licensing agreements between ARIA and the respective data providers.

The starting point of the analysis was the PARIX business register [19], a national administrative database containing official records of registered economic activities, including their legal status, economic activitiy classification according to the Italian statistical coding system (ATECO, “Attività Economiche”), and basic identifying information such as business name and registered address. From this source, 22,342 active food service establishments belonging to the selected ATECO categories (i.e., food service businesses with and without table service, codes 56.11.11 and 56.11.12 respectively) were extracted as the reference population. At this stage, businesses were identified exclusively through their registered name and address.

The resulting dataset was subsequently enriched through the integration of activity-type information (such as restaurant type, cuisine type, etc.) available from multiple complementary data sources, namely OpenStreetMap (OSM), The Data Appeal (TDA), and Google Places. Since no common identifier linking records across the different data sources was available, records from the complementary sources were reconciled with the PARIX businesses through a progressive record-matching procedure based on address and business name normalization, and both exact and threshold-based fuzzy string-matching algorithms. This process enabled the association of each PARIX establishment with the corresponding information available from the complementary sources, whenever a reliable match was found.

A preliminary candidate selection was then performed using structured semantic information available from these external sources, including business categories, Google Places types, and OSM tags related to restaurants, pizzerias, bakeries, and other activities potentially associated with wood-fired cooking. This filtering stage reduced the initial data population from 22,342 establishments to 9,993 candidate businesses (44.7%), while 12,349 establishments (55.3%) were filtered out due to the absence of keywords indicative of wood-fired cooking in their primary metadata. The aim of this stage was to exclude businesses clearly unrelated to the target activity, thereby substantially reducing the number of establishments subjected to the successive computationally intensive analyses.

For each record in the candidate set, up to five images and ten user reviews were retrieved through the Google Places API under the commercial agreement. User reviews were available for 8,486 establishments (84.9% of the candidate set), while 8,381 establishments (83.9%) had at least one image available for visual analysis. Images were analyzed using a CLIP-based vision-language model to identify visual evidence of traditional wood-fired ovens and using Optical Character Recognition (OCR) to detect textual references appearing in signs, menus, or other images. In parallel, user reviews were analysed using Natural Language Processing techniques. In addition to conventional keyword matching, a multilingual Natural Language Inference (NLI) model was fine-tuned to determine whether the semantic content of each review provided evidence supporting the hypothesis that the restaurant used a wood-fired oven. This approach improves robustness compared with simple lexical matching by enabling the recognition of semantically equivalent expressions in multiple languages and reducing ambiguities associated with context-dependent mentions.

Finally, the outputs of the different analytical modules were combined through a conservative logical aggregation strategy, whereby positive evidence generated by any independent source was considered sufficient to classify an establishment as equipped with a wood-fired oven. The proposed methodology identified 6,339 establishments as potentially equipped with a wood-fired oven, corresponding to 63.4% of the candidate establishments analyzed through the multimodal pipeline. Among the identified establishments, 6,319 (63.2%) were detected through visual analysis, 167 (1.7%) through OCR-based textual evidence, and 126 (1.3%) through semantic analysis of user reviews. Note that some establishments were detected by more than one method; the counts above therefore sum to more than the total of 6,339 unique identified establishments in the metropolitan area. As a final step, only establishments falling within the administrative boundaries of the Municipality of Milan were retained, resulting in a dataset comprising of 3,743 wood-fired ovens. The resulting inventory constitutes the basis for atmospheric dispersion modelling presented in the following sections.

2.3. Additional Major Emission Sources in the Area

In Milan’s urban environment, the INEMAR emission inventory estimates that more than 80% of NO_x emissions originate from road traffic and non-industrial combustion. These sources, together with population-density-related activities such as tobacco smoking and fireworks, also account for a large share of primary PM₁₀ emissions [10].

To validate the results obtained for wood-fired ovens, road traffic and heating systems were also considered in the modelling study. The emission estimates were taken from the regional emission inventory (Table 2).

Specifically, total road-transport emissions (Macrosector Level 07), excluding road dust resuspension (Activity Level 07), were used for road traffic. Instead, total emissions from non-industrial combustion (Macrosector 02) were considered for heating systems, excluding plants from agriculture, forestry and aquaculture facilities (Sector Level 03) and pizzerias with wood-fired ovens (Sector Level 01, Activity Level 07) to avoid double counting.

The spatial distribution in the urban areas was calculated based on the public register of thermal plants in Lombardy [20] and on highly detailed information on traffic intensity in the city of Milan made available by the Lombardy Region Administration.

2.4. Air Pollutant Dispersion Simulation Model

Atmospheric dispersion simulations were performed using the CALPUFF modelling system, driven by meteorological fields generated by the WRF (Weather Research and Forecasting) model for the year 2019. The meteorological domain consisted of a 40×36 cell grid with a horizontal resolution of 500 meters, covering the Municipality of Milan. The south-western corner of the domain was located at coordinates (502000, 5025000) in the UTM coordinate system (WGS84, Zone 32N). The same domain was adopted as the computational domain for the CALPUFF simulations.

To characterize the dispersion patterns associated with wood-fired pizza ovens, an initial reference simulation was performed by representing the emission source as a point source located near the center of the modelling domain with unit emissions (1 g s⁻¹). This configuration was adopted to minimize potential boundary effects and to obtain representative first-guess concentration fields. To validate the modelling approach used for wood-fired pizza ovens, analogous simulations were performed for heating systems, represented by a point-source boiler, and for road traffic, represented as an area source (50m × 50m).

Source-specific temporal emission profiles were implemented in the CALPUFF model (Figure 3). For wood-fired ovens, hourly emission variations were based on activity patterns associated with meal preparation periods. Heating systems emissions were modulated daily using heating degree days (HDD) as a proxy for heating demand, calculated based on a previously validated method [21].

The modulation factors $f_i$were expressed through the i-th daily contribution of HDD to the annual HDD total (Equation 1), while they were set to zero during the non-heating season, from 15 April to 15 October:

$$f_i = {\displaystyle \frac{21-\overline{T_i}}{{\sum }_{d=1}^{365}\left(21-\overline{T_d}\right)}},$$

(1)

where $\overline{T_i}$ is the average temperature of the i-th day, and $\overline{T_d}$ is the average temperature of each day of 2019. Road traffic emissions were adjusted using detailed temporal profiles derived from traffic intensity data.

Post-processing of model results was performed using CALPOST and a linear combination workflow to calculate the cumulative contribution of all wood-fired ovens, heating boilers, and 50m × 50m road traffic areas identified within the Municipality of Milan. The methodology varied according to source type, with distinct approaches adopted for road traffic and heating systems versus wood-fired ovens. Both approaches were based on the assumption that meteorological conditions are sufficiently homogenous across the metropolitan area, allowing the concentration field generated by a representative source located around the center of the domain to be spatially translated while preserving the validity of the associated meteorological conditions. The contribution of each source was therefore calculated independently and subsequently aggregated to obtain the overall concentration field. This linear combination strategy reduced computation times by 99% while allowing the inclusion of all identified sources without limitations on the number of simulated sources within each run.

The dispersion model starts from the CALMET meteorological domain previously defined. A nesting factor of 2 was applied, resulting in a grid composed of cells $c_{ij}$ measuring 250m × 250m. In the approach adopted for road traffic and heating systems, for each grid cell, a source-specific weighting factor was calculated and used to spatially distribute emissions across the modelling domain.

For road traffic, the weight associated with each grid cell $w_{ij}^{traffic}$ is proportional to the annual number of vehicles crossing the area:

$$w_{ij}^{traffic}={\displaystyle \frac{V_{ij}}{\sum _{c_{ij}\in M}{V}_{ij}}} ,$$

(2)

where, $V_{ij}$ is the annual traffic volume associated with each grid cell, and $M$ is the set of grid cells located within the municipality of Milan.

For heating systems, the weight associated with each grid cell $w_{ij}^{heating}$ is proportional to the thermal power of all boilers located within that cell:

$$w_{ij}^{heating}={\displaystyle \frac{{\displaystyle {\sum }_{k=1}^{n_{ij}}}{P}_{k,ij}}{{\sum }_{c_{ij}\in M}{\displaystyle {\sum }_{k=1}^{n_{ij}}}{P}_{k,ij}}} ,$$

(3)

where, $P_{k,ij}$ is the thermal power of boiler $k$ located in grid cell $c_{ij}$, and $n_{ij}$ is the number of boilers within grid cell $c_{ij}$.

Only emission sources located within the municipality of Milan were considered. Consequently, non-zero weights were assigned exclusively to cells belonging to the municipality:

$$w_{ij}^{\left(s\right)}=\left\{\begin{array}{c}{w}_{ij}^{\left(s\right)}, c_{ij}\in M\\ 0, c_{ij} \notin M\end{array}\right.$$

(4)

where, $s$ denotes the considered source category: traffic or heating.

These weights were then used to calculate the emission intensity within each grid cell. The share of total emission values allocated to each grid cell ($E_{p,ij}^{\left(s\right)}$) was calculated using:

$$E_{p,ij}^{\left(s\right)}=E_p^{\left(s\right)}{w}_{ij}^{\left(s\right)}$$

(5)

where $E_p^{\left(s\right)}$ is the total annual emission of pollutant $p$ from source category $s$. Annual emissions of the air pollutants investigated in the municipality of Milan were taken from the INEMAR inventory, as described in Section 2.3.

A different approach was adopted for wood-fired ovens when calculating the spatial allocation of total emissions. Although Equation 5 was also applied to this source category, the experimentally derived mass fluxes (Section 2.1) were used to estimate the emissions ($E_p^{\left(s\right)}$), which were then weighted by the total number of wood-fired ovens in each cell ($w_{ij}^{\left(s\right)}$), identified using the data-driven methodology described in Section 2.2.

Simulation results from all three sources where then post-processed with CALPOST to generate a first-guess concentration field, $C_{unit}^{\left(s\right)}(x,y)$. For each grid cell, this unit concentration field associated with road traffic and heating systems was spatially translated so that the reference source location coincided with the center of the destination cell, and subsequently scaled according to the corresponding allocated emission, $E_{p,ij}^{\left(s\right)}$. A similar procedure was applied to wood-fired ovens; however, in this case, the translations were based on the actual locations of the ovens rather than on the grid cell centers.

The contribution of each grid cell to the concentration at location $(x,y)$ was therefore:

$$C_{p,ij}^{\left(s\right)}\left(x,y\right)=E_{p,ij}^{\left(s\right)}{C}_{unit}^{\left(s\right)}\left(x-x_{ij},y-y_{ij}\right)$$

(6)

where, $\left(x_{ij}, y_{ij}\right)$ represent the coordinates of the center of the grid cell (road traffic and heating systems) or the location of the wood-fired pizza ovens.

The final concentration field for each pollutant and source category was obtained through the superposition of the contributions of all cells:

$$C_p^{\left(s\right)}(x,y)=\sum _{c_{ij}\in M}{C}_{p,ij}^{\left(s\right)}(x,y)$$

(7)

or, equivalently,

$$C_p^{\left(s\right)}\left(x,y\right)=E_p^{\left(s\right)}\sum _{c_{ij}\in M}{w}_{ij}^{\left(s\right)}{C}_{unit}^{\left(s\right)}\left(x-x_{ij},y-y_{ij}\right).$$

(8)

Simulations were carried out for PM₁₀, PM_2.5, and NO_x and the resulting concentration fields represent the annual average values attributable to wood-fired pizza ovens, traffic and residential heating throughout the Municipality of Milan.

In addition, concentration values were extracted at the location of the air-quality monitoring stations operated by ARPA Lombardia within the study area (Table 3, Figure 4). The monitoring network in Milan comprises five air-quality stations representative of different exposure settings, including urban traffic (UT) and urban background (UB) environments. The station typology and the associated pollutants/tracers investigated in the analysis are also indicated in Table 3.

These results were subsequently used to evaluate the contribution of wood-fired pizza ovens, heating systems, and road traffic to ambient pollutant concentrations and to compare model estimates with observed air quality measurements.

2.5. Modelling Approach

Figure 5 summarizes the overall methodological framework adopted in this study to assess the air quality impact of wood-fired pizza ovens. The approach integrates information from multiple sources, including geospatial data, experimentally measured emission factors, atmospheric dispersion modelling, and air quality measurements.

The workflow starts with the identification and geolocation of pizzerias and restaurants equipped with a wood-fired oven through a data-driven procedure. The resulting geospatial dataset is then combined with experimentally determined emission factors for the investigated air pollutants and tracers, allowing the calculation of source-specific emission fluxes for each identified facility. Meteorological fields generated with the WRF model were used to produce the three-dimensional meteorological input required by the CALPUFF dispersion model. Reference dispersion simulations were performed for unitary emission sources located within the modelling domain, post-processed using CALPOST to derive concentration fields that were subsequently spatially translated according to the actual location of each source and aggregated to estimate their cumulative impact. The resulting concentrations were then compared with observations from the regional air-quality monitoring network.

2.6. Measurements from the ARPA Lombardia Air Quality Network (Methodology and Analysis of PM Composition)

The air quality monitoring network operated by ARPA Lombardia consists of 83 fixed monitoring stations distributed across urban, suburban, rural background, and industrial environments, providing long-term observations of atmospheric pollutants and ensuring representative spatial coverage of the Lombardy region. Within this framework, the present study focuses on the monitoring stations located within the Municipality of Milan (Figure 4).

Mass concentrations of PM₁₀ and PM_2.5 are continuously measured using beta attenuation monitors [22], while nitrogen oxides (NO, NO₂, and NO_x) are continuously monitored according to the European reference method EN 14211, based on the chemiluminescence measurement principle [23]. At the MI-P_UB, MI-S_UT, and MI-M_UT sites, PM₁₀ samples are collected using reference gravimetric samplers and beta attenuation monitors, aimed at determining the chemical composition of atmospheric particulate matter.

In the present study, particular attention was also devoted to the analysis of specific tracers such as BaP and LG. The former is routinely monitored as a marker of incomplete combustion processes and is recognized as one of the most relevant carcinogenic components associated with atmospheric particulate matter [24,25]. Polycyclic aromatic hydrocarbons (PAHs), including BaP, were determined according to the European reference method EN 15549, with a sampling frequency of one sample every three days, as prescribed by current European air-quality legislation, [26]. Instead, LG was determined daily by ion chromatography coupled with pulsed amperometric detection (IC–PAD), providing sensitive and selective quantification of this compound in particulate matter samples.

The minimization of measurement errors and the proper estimation of measurement uncertainty require the reduction of both random and systematic sources of uncertainty. This objective is achieved through the implementation of comprehensive Quality Assurance and Quality Control (QA/QC) procedures, in accordance with the requirements of Directive 2008/50/EC and the Italian Ministerial Decree of 30 March 2017 (DM 62/2017) [27,28]. Quality Control activities include routine operational and procedural checks performed by ARPA Lombardia to verify instrument performance, calibration stability, and data validity, ensuring compliance with the data-quality objectives established by European legislation. Quality Assurance activities comprise regular participation in proficiency-testing schemes, inter-laboratory comparison exercises, blind tests, and performance audits aimed at assessing measurement accuracy, comparability, traceability, and long-term consistency. The implementation of these QA/QC procedures ensures the reliability, traceability, robustness, and scientific integrity of the air-quality measurements and chemical-speciation data used in the present study.

2.7. Source Apportionment of Particulate Matter

Among the objectives of the PrepAIR project [29], the activation of a measurement network for the chemical characterization of PM₁₀ was envisaged. This network was created based on existing monitoring stations and consists of four urban background sites - Turin, Milan, Vicenza and Bologna - and one rural background site - Schivenoglia - to which, during the project, the suburban background site of Cavallermaggiore (CN) was added. By means of specific laboratory analyses of the samples collected daily at the 6 sites, a database was created that covers a total of six years (from 1 April 2018 to 31 March 2024). A source apportionment analysis was then performed, carried out with Positive Matrix Factorization (PMF) [30].

The source apportionment analysis allowed the identification, albeit with margins of uncertainty, of the main sources and/or processes of particulate formation. In summary, 7 factors common to all receiving sites were resolved: two factors of secondary origin, Secondary Background and Winter Secondary, as well as Biomass Burning, Traffic, Salt, Soil and Anthropogenic Mix. On average over the Po basin, the two factors Soil and Salt together explain more than 15% of the reconstructed PM₁₀ mass. Soil has a contribution of more than 10%, while salt has a contribution of about 5%. The primary factors of anthropogenic origin, Traffic and Biomass Combustion have a similar contribution, around 10%. The secondary component represents about 50% of the total mass of PM₁₀ and is divided equally into the two factors that characterize it: Winter Secondary and Secondary Background. In particular, the secondary background is enriched with a carbonaceous component which probably indicates that the ammonium sulphate may have formed very far from the receptor, thus enriching itself during an aging process. As also pointed out in the analysis of the chemical composition, the contributions of the different sources are uniform over the basin. It should be noted that there is an accumulation of ammonium nitrate during the winter months, which leads the winter secondary to account for over 30% of the PM₁₀ mass in January and February.

PM₁₀ sampling has been active at the Milan site since 2013, so a long series of chemical composition is available, which combines the measurements made during the project with those collected previously. In the last decade, there has been a decline in the total mass of PM₁₀ and in the concentration of elemental carbon (EC); there is therefore a reduction in the contribution of the factor associated with traffic tailpipe emissions. The fraction of winter secondary (ammonium nitrate in particular) is instead rather stable.

3. Results

3.1. Validation of Pizza Ovens Identification Method

To evaluate the reliability of the identification methodology, a validation dataset was constructed by randomly selecting 1,713 establishments (approximately 10% of the sample). The presence or absence of a wood-fired oven was then determined through manual review, using evidence such as venue photographs, customer reviews, and business websites. This manual assessment provided the ground truth against which the automatic classification was evaluated.

A total of 276 businesses were confirmed to operate a wood-fired oven. The classifications generated by the employed methodology were then compared against the manually established ground truth labels to evaluate its ability to correctly identify potential emission sources.

Performance was assessed using Recall, Precision, and F1-score, computed against the manually established ground truth. Recall $\left(R\right)$ was defined as:

$$R={\displaystyle \frac{TP}{TP+FN}} ,$$

(9)

where $TP$ denotes true positives and $FN$ denotes false negatives, representing the proportion of establishments with a wood-fired ovens that were correctly identified.

Instead, precision $\left(P\right)$ was defined as:

$$P={\displaystyle \frac{TP}{TP+FP}} ,$$

(10)

where $FP$ denotes false positives, representing the proportion of identified establishments that were confirmed to have a wood-fired oven.

Finally, the F1-score was defined as the harmonic mean of $R$ and $P$, providing a summary measure of their trade-off:

$$F1={\displaystyle \frac{2\left(P\times R\right)}{P+R}} ,$$

(11)

Emphasis was placed on Recall, as the primary objective of our approach was to minimize the number of emission sources omitted from the inventory. In the context of atmospheric emission assessment, missing a true wood-fired oven has a greater impact than temporarily including a false positive.

The proposed methodology achieved a Recall of 77%, indicating that more than three quarters of the establishments equipped with a wood-fired oven were correctly identified (Table 4). The corresponding Precision was 20%, while the overall F1-score reached 31%, reflecting the expected trade-off resulting from a methodology intentionally optimized for coverage rather than strict selectivity.

An additional analysis was conducted to investigate the performance across different business categories. The highest results were obtained for establishments explicitly classified as pizzerias, achieving a Recall of 87% and a Precision of 50% (Table 5). Conversely, more generic business categories (e.g., generic restaurant) generated a substantially larger number of false positives, confirming that the lack of specificity of the available business metadata is the primary limiting factor for classification precision.

The validation also confirmed the effectiveness of the structured semantic filtering adopted during the candidate selection stage. Business tags collected from PARIX, OpenStreetMap, The Data Appeal, and Google Places were consistent with the manually assigned labels in 85.7% of the analyzed establishments. This result demonstrates that integrating structured metadata from multiple heterogeneous sources provides a robust basis for constructing the candidate inventory before applying computationally intensive multimodal analyses. The remaining mismatches primarily concern generic or ambiguous business categories, indicating that the specificity of source classifications remains a primary constraint on pipeline precision.

3.2. Validation of the Modelling Approach

The reliability of the parallel processing method was evaluated by comparing its results against those obtained from a conventional CALPUFF simulation in which all identified wood-fired pizza ovens were explicitly included as individual emission sources. Both approaches were applied to the same geospatial database of facilities equipped with wood-fired ovens and were driven by identical emission rates. Annual average concentration fields obtained with the two methods were then compared on a cell-by-cell basis over the entire modelling domain (Figure 5).

The comparison showed an excellent agreement between the two approaches, confirming the validity of the proposed methodology. A strong linear relationship was observed between the concentrations estimated by the conventional CALPUFF simulations and those derived from the parallel processing approach (r = 0.986). Furthermore, the regression line (y = 1.00x − 0.04) is nearly coincident with the 1:1 reference line, confirming the near-linear additivity of source contributions and the absence of significant systematic bias. Although a small dispersion is observed at higher concentration levels, the overall results demonstrate that the proposed parallel processing method can accurately reproduce the concentration fields generated by traditional multi-source CALPUFF simulations in the investigated domain, while significantly reducing computational requirements.

It should be noted that the quality of the agreement depends on the location selected for the test source prior to the linear combination. In this study, the test source was placed at the location providing the best agreement with the reference multi-source CALPUFF simulation, which was considered the most representative of the meteorological conditions across the modelling domain. Consequently, the methodology was adopted for the assessment of the cumulative impact of all wood-fired pizza ovens, heating systems, and road traffic within the study area.

3.3. Spatial Distribution of Major Air Pollutants

Annual average air pollutant concentrations obtained from the dispersion simulations were analyzed to assess the impact of wood-fired pizza ovens and other major urban emission sources on air quality within the study area. Figure 6 presents the spatial distribution of annual average PM₁₀, PM_2.5, and NO_x concentrations associated with the investigated sources. For heating systems, average concentrations were calculated considering only simulated values in the winter period during which these sources are active.

The spatial distribution of the modelled pollutant concentrations exhibits distinct patterns depending on the emission source. For wood-fired pizza ovens, the highest contributions to all the pollutants investigated are concentrated in the central districts of the city, with progressively lower concentrations towards the outskirts. This pattern reflects the spatial distribution of restaurants and pizzerias equipped with wood-fired ovens, which are predominantly located in the city center where commercial activities, population density and recreational services are more concentrated.

A different spatial pattern is observed for road traffic emissions. The concentration fields of all investigated pollutants closely follow the structure of Milan’s road network, with elevated values along the major arterial roads and ring-road system and lower contributions in less busy areas. The strong correspondence between modelled concentration gradients and the actual road layout indicates that the dispersion modelling framework adequately reproduces the spatial influence of traffic-related emissions.

Instead, emissions from heating systems show a distribution that is broadly consistent with both residential density and the spatial concentration of commercial and institutional activities, in line with the pattern observed for wood-fired ovens. Indeed, higher contributions were observed in the most densely populated and built-up areas of the municipality, where residential buildings coexist with offices, schools, retail establishments, and other service-sector activities. Conversely, lower concentrations characterize peripheral zones with lower population density and fewer commercial and institutional facilities. This pattern is evident for both particulate matter and nitrogen oxides and reflects the spatial distribution of heating demand and energy use across the city.

Although exhibiting distinct spatial patterns, the model results indicate that road traffic is the dominant contributor to particulate matter concentrations (0.075 – 10.42 µg m⁻³ PM₁₀), followed by wood-fired pizza ovens (0.085 – 3.34 µg m⁻³ PM₁₀) and, to a lesser extent, heating systems (0.0054 – 0.61 µg m⁻³ PM₁₀). These findings highlight the non-negligible role of emissions from restaurants and pizzerias as a source of primary particulate matter in the urban environment. In contrast, the contribution of wood-fired ovens to NO_x concentrations is relatively limited (0.0044 – 1.75 µg m⁻³). Road traffic remains the principal source of NO_x across the municipality (1.05 – 146.23 µg m⁻³), while heating systems also provide a significant contribution during the winter season (0.14 – 16.25 µg m⁻³). Overall, the results suggest that particulate matter concentrations arise from a combination of sources with varying relative importance, whereas NO_x levels are driven primarily by traffic emissions, with heating systems representing an additional relevant source, particularly during periods of increased demand.

3.4. Impact Assessment on Air Quality

To assess the impact of wood-fired pizza ovens relative to other major urban emission sources, simulated pollutant concentrations were extracted at selected receptor locations corresponding to the ARPA air-quality monitoring stations within the Municipality of Milan. Daily average concentrations derived from the dispersion model were then compared with the corresponding measurements collected at each monitoring site, allowing for an evaluation of the modelling framework's ability to reproduce the observed temporal variability of pollutant concentrations. Figure 7 shows the results obtained for NO_x at the MI-P_UB monitoring station (data from other stations can be found in the SM: Figure S1-S4).

The results indicate that the model can reproduce almost all the observed NO_x concentrations at the monitoring stations, accounting for between 95.16% and 123.06% of the measured values. Model performance was further supported by a moderate positive correlation between simulated and observed concentrations (r = 0.50 - 0.62), indicating a good representation of the temporal variability of NO_x levels. Road traffic was identified as the dominant contributor throughout the year, accounting for the largest share of NO_x concentrations at all monitoring sites (89.75 – 92.29%). During the winter season, heating systems also exerted a significant influence on ambient NO_x levels (12.37 – 27.67%), whereas the contribution of wood-fired pizza ovens remained negligible in comparison with the other emission sources (0.45 – 2.09%).

Different trends were observed for particulate matter, with no significant differences between PM₁₀ and PM_2.5. Unlike NO_x, the relative contribution of wood-fired pizza ovens to total simulated particulate matter concentrations was substantial at all monitoring stations (16.97 – 48.27%), indicating an influence that extends across the entire urban area. At the MI-P_UB urban background station, road traffic remained the dominant source, although wood-fired ovens also accounted for a significant share of particulate matter concentrations (Figure 8, Figure S5). A similar pattern was observed at the MI-M_UT station, where the contribution from road traffic was particularly pronounced, while the relative impact of pizzerias was lower, consistent with the nature of the station and its location outside of the city center (Figure S6).

By contrast, stations located closer to the city center, such as MI-S_UT and MI-V_UT, where the density of restaurants and pizzerias is greatest, exhibited a markedly higher contribution from wood-fired pizza ovens to total simulated particulate matter concentrations (Figure 9, Figure S7 and Figure S8). These results are consistent with the spatial patterns described in the previous section, which showed a stronger influence of wood-fired oven emissions in the central districts of Milan and a gradual decline towards the urban periphery.

This effect was particularly evident during the summer months, when heating systems are inactive and traffic-related emissions are generally reduced. Under these conditions, wood-fired pizza ovens emerged as one of the main contributors to particulate matter, in some locations exceeding the contribution of road traffic. Figure 9(d) highlights the significant contribution of wood-fired ovens to PM₁₀ concentrations at the MI-V_UT site. During the summer months, this source emerges as the dominant contributor to PM₁₀ levels, accounting for up to 80% of the total concentration. Unlike heating emissions, which are strongly seasonal, and traffic emissions, which exhibit substantial temporal variability, restaurant and pizzeria activities operate throughout the year and therefore provide a relatively constant source of primary particulate matter emissions.

Overall, the findings indicate that the contribution of wood-fired pizza ovens to urban particulate matter is strongly influenced by both local source density and seasonal conditions. In areas characterized by a high concentration of food-service activities, their impact can become comparable to, or even exceed, that of more traditionally recognized urban sources.

Despite the relatively important contribution of wood-fired pizza ovens within the simulated source categories, their contribution to the total observed particulate matter concentrations remains limited. Indeed, the model accounts for only 11.98 – 17.78% of measured PM₁₀ concentrations and 16.92 – 25.69% of the measured PM_2.5 concentrations. This is because the simulations include only primary particulate matter emissions from three of the major urban sources, whereas a substantial fraction of ambient particulate matter originates from secondary aerosol formation processes and from regional and long-range transport, which are not represented in the source-specific simulations [31]. This interpretation is supported by data from the Copernicus Atmosphere Monitoring Service (CAMS), which indicates that only 32.64% of PM₁₀ concentrations in Milan originate from local emissions within the city, whereas more than 65% are associated with contributions from outside the municipal boundaries [32].

Further support is provided by Positive Matrix Factorization (PMF) analyses performed on measurements collected at the MI-P_UB monitoring station. These analyses indicate that biomass burning and road traffic together account, on average, for 26.82% of the total PM₁₀ concentration, a fraction that could also include part of secondary particulate matter formed from precursors emitted by these sources. Therefore, the contribution estimated by the modelling framework is consistent with the range of impacts identified through receptor-based source apportionment techniques.

Overall, the agreement between the model results and PMF estimates, together with the good positive correlation observed between simulated and PMF-derived concentrations (Figure 10), suggests that the modelling framework provides a realistic representation of the primary particulate matter contribution from the investigated urban sources and is consistent with both observational evidence and previous source-apportionment studies.

3.5. Wood-Fired Pizza Ovens Impact Assessment on Selected Tracers

A similar analysis was also performed to assess the contribution of wood-fired pizza ovens to particulate-bound levoglucosan (LG), a tracer of biomass burning, and benzo(a)pyrene (BaP), a marker of incomplete combustion. Emission fluxes for these species were estimated as fractions of PM₁₀ emissions; consequently, their dispersion maps exhibit the same spatial pattern as that shown in Figure 6, with proportionally lower absolute concentrations. Hence, as was the case for particulate emissions, the contribution of wood-fired ovens to both tracers is highest in the city center and progressively decreases towards the outskirts.

The comparison between simulated and observed LG concentrations at the urban background monitoring station of MI-P_UB indicates a good level of agreement (r = 0.480), demonstrating that the model can reproduce a substantial fraction of the observed spatial and temporal variability. Instead, only a moderate level of agreement (r = 0.343) was observed between simulated and observed values at the MI-S_UT monitoring station. However, in both cases, the model captures the pronounced seasonal pattern, with higher concentrations during winter and lower values during summer, reflecting once again the influence of meteorological conditions on pollutant accumulation and dispersion (Figure 11, Figure S9).

The simulated wintertime contribution of wood-fired pizza ovens was lower in absolute terms, reflecting the influence of other biomass-burning sources, particularly residential heating, which is employed during the cold season. Nevertheless, wood-fired ovens still accounted for a substantial share of LG concentrations, contributing on average 19.70% of the total levels at MI-P_UB and 34.28% at MI-S_UT. This is in line with previous observations which showed a greater concentration of restaurants and pizzerias around the MI-S_UT monitoring station. These devices therefore represent a non-negligible additional source of biomass-burning emissions in urban environments, where the use of woody biomass for residential heating is generally less common than in rural areas.

In contrast, measured LG concentrations during the summer months frequently fall below the limit of detection, whereas the modelled contribution from wood-fired ovens remains non-negligible because pizzerias and other establishments operating these appliances remain active throughout the entire year. In some cases, the discrepancy can be explained by the fact that the simulated concentrations are below the experimental detection limit and therefore cannot be reliably quantified in particulate matter samples. In other cases, the difference between modelled and measured concentrations may be related to the atmospheric degradation of LG during the warm season. Although emitted continuously by wood-fired ovens, LG is susceptible to oxidation by hydroxyl radicals, and its degradation is enhanced under conditions of elevated solar radiation and intense photochemical activity [33]. Consequently, its atmospheric lifetime is reduced during summer, leading to lower ambient concentrations than would be expected based solely on emissions and dispersion processes.

Instead, the comparison between modelled and measured BaP concentrations shows a weaker level of agreement at both monitoring sites: r = 0.223 and r = 0.261 at the MI-P_UB and MI-S_UT monitoring sites, respectively (Table S1 and Table S2). This suggests that additional combustion sources likely play a significant role in determining the spatial and temporal variability of ambient BaP concentrations across the study area, including residential biomass burning, vehicular traffic, and industrial activities [34]. Nevertheless, the model was able to reproduce the main seasonal pattern, which broadly mirrored that observed for levoglucosan, with higher concentrations during the winter months and lower levels during summer.

Despite the weaker correlation with experimental observations, the simulations indicate a substantial contribution of wood-fired ovens to atmospheric BaP concentrations, accounting on average for 34.95% and 37.47% of the measured levels at MI-P_UB and MI-S_UT, respectively. These results highlight that wood-fired ovens can represent a significant source of combustion-related pollutants in urban environments and may contribute appreciably to the overall BaP burden. At the same time, the similar contribution estimated at the two monitoring sites, despite their different urban settings, supports the hypothesis that wood-fired pizza ovens are not the primary drivers of the spatial variability of BaP concentrations. Unlike levoglucosan, whose spatial distribution appears to be more closely associated with the density of wood-fired ovens, BaP is influenced by a broader range of combustion sources, which likely exert a stronger control on the observed concentration patterns across the city.

4. Discussion

This study represents urban air pollution as the combined effect of three major emission sources: wood-fired pizza ovens, residential and commercial heating systems, and road traffic. By coupling emission inventory and experimental data with the developed dispersion modelling framework, spatially resolved concentration profiles of NO_x and PM₁₀ were reconstructed across a representative west-to-east transect of the city of Milan (Figure 12). The resulting patterns are consistent with those expected in a densely urbanized environment [35], with road traffic emerging as the dominant contributor to NO_x levels and displaying marked peaks in correspondence with the A50 and A51 ring roads. In contrast, PM₁₀ concentrations reveal a more diversified source contribution, where traffic remains important but emissions from wood-fired pizzerias become comparatively more relevant within the urban core.

In urban environments, ambient pollutant concentrations are generally composed of a regional or urban background component and local concentration increments generated by specific emission sources [12,35]. Figure 12 highlights the capability of the modelling framework developed in this study to quantify these source-specific urban increments and distinguish the contributions of individual emission sources. The high-resolution concentration profiles reveal street-level concentration increments and distinct spatial patterns across the city, with the influence of heating systems and wood-fired pizza ovens increasing towards the urban core, while traffic-related concentration peaks were approximately symmetric with respect to the city center and coincide with major ring roads. Furthermore, the comparison between model results and experimental measurements suggests that the background concentration plays a larger role for PM₁₀ and PM_2.5 than for NO_x. This finding is consistent with previous PMF and CAMS analyses, which showed that the largest fraction of atmospheric particulate matter originates from the regional background rather than from local urban sources.

In previous studies, the city increment was derived from pollutant concentrations measured at one or more air quality monitoring stations on the previous day, thereby accounting for the temporal accumulation of pollutants in urban areas adopting the approach of SHERPA-City project [12,36,37]. Piccoli et al. [38] showed that the CAMx-LPiG hybrid model improves the representation of traffic-related NO₂ concentration gradients in Milan and reduces simulation bias relative to conventional CAMx modelling, enabling more accurate urban-scale air quality assessments [38]. In this study, however, the proposed approach makes it possible to estimate the effects of the spatiotemporal variability of emissions from traffic, residential heating, and pizzerias in relation to atmospheric accumulation and transport conditions. Unlike approaches based on monitoring data, neither measurements from air quality monitoring stations nor background concentration data are used, since the simulations are performed continuously for every day of the simulated year.

5. Conclusions

This work presented a data-driven methodology to support the construction of an inventory of establishments equipped with wood-fired ovens, addressing the lack of official information on this type of diffuse emission source. By integrating administrative business records with user generated digital information and multimodal artificial intelligence techniques, the proposed approach enabled the identification of establishments to be included in the emission inventory used for the subsequent air quality assessment.

The validation results showed that the methodology provides a reliable basis for inventory construction, achieving a Recall of 77% while deliberately prioritising the identification of potential emission sources over the reduction of false positives. This design choice is consistent with the objectives of environmental applications, where omitting actual emission sources may have a greater impact than including false positives. The analysis also confirmed that the lack of specificity in available business metadata represents the primary limiting factor for classification precision.

Although developed for the identification of wood-fired pizza ovens, the proposed methodology is inherently reusable. The same workflow can be adapted to other commercial activities for which official inventories are unavailable, integrating heterogeneous administrative records, digital content, and artificial intelligence to identify potential emission sources. Future applications will explore the transferability of this framework to vehicle body repair shops, where identifying facilities equipped with paint booths and curing ovens is a prerequisite for estimating VOC emissions. More generally, the proposed approach provides a scalable methodology to support public administrations in the development and maintenance of emission source inventories for different categories of diffuse anthropogenic emissions.

The local impact of pizzerias on urban air quality has been estimated to be in the range of 2 - 3 µg m⁻³ of PM₁₀, based on the annual average concentration. This impact is not negligible and confirms that modelling systems need to be complemented by new data-driven approaches for source identification. The comparison of the simulation results also confirmed observations from monitoring stations, PMF and findings from the CAMS modelling system, which indicate that most PM₁₀ concentrations are attributable to physicochemical transport processes associated with emission sources located outside the urban area.

Moreover, the integrated air quality assessment system developed in this paper is aligned with the objectives of the new EU Ambient Air Quality Directive (EU 2024/2881). The system combines air quality monitoring data, indicative measurements, and advanced modelling techniques within a unified framework, placing modelling at the core of the assessment process. It provides high-resolution urban-scale air quality maps and supports microscale analyses in areas affected by pollution hotspots and steep spatial concentration gradients. Furthermore, the system is specifically designed to capture population exposure with a high degree of spatial detail, enabling the identification of local variations in pollutant concentrations over short distances and supporting more accurate exposure assessments in urban environments.

The proposed modelling approach reduced computational time by approximately 99% while enabling the inclusion of all identified emission sources without limitations on the number of sources simulated in each run. The system captures the spatiotemporal variability of emissions from traffic, residential heating, and pizzerias, accounting for atmospheric accumulation and transport processes. Moreover, unlike approaches based on monitoring or externally simulated data, it operates independently of air quality measurements and background concentration data, providing continuous simulations for every day of the year.

Considering Equation 8, the proposed approach could enable the development of a modeling system capable of attributing concentration contributions to each type of emission source, as well as to its spatial location within the computational domain. This capability is particularly valuable for inverse modeling applications, which aim to reconstruct emissions through the comparison of simulated data with observations from air quality monitoring stations and satellite measurements.