*Article* **Impact of Land Use on Atmospheric Particulate Matter Concentrations: A Case Study of the Beijing–Tianjin–Hebei Region, China**

**Haoran Zhai 1,2 , Jiaqi Yao 2,3,\* , Guanghui Wang 2,4 and Xinming Tang 2,3**

	- Xuzhou 221116, China

**Abstract:** The increasing frequency of human activities has accelerated changes in land use types and consequently affected the atmospheric environment. In this manuscript, we analyze the relationships between the particulate matter concentration and land use changes in the Beijing–Tianjin–Hebei (BTH) region, China, from 2015 to 2018. The experimental results indicate that (1) an improved sine function model can suitably fit the periodic changes in the particulate matter concentration, with the average R<sup>2</sup> value increasing to 0.65 from the traditional model value of 0.49, while each model coefficient effectively estimates the change characteristics of each stage. (2) Among all land use types, the particulate matter concentrations in construction land and farmland are high, with a large annual difference between high and low values. The concentration decreases slowly in spring and summer but increases rapidly in autumn and winter. The concentrations in forestland and grassland are the lowest; the difference between high and low values is small for these land use types, and the concentration fluctuation pattern is relatively uniform. Natural sources greatly influence the concentration fluctuations, among which frequent dusty weather conditions in spring impose a greater influence on forestland and grassland than on the other land use types. (3) The landscape pattern of land use exerts a significant influence on the particulate matter concentration. Generally, the lower the aggregation degree of patches is, the higher the fragmentation degree is, the more complex the shape is, the higher the landscape abundance is, and the lower the particulate matter concentration is. The higher the construction land concentration is, the more easily emission sources can be aggregated to increase the particulate matter concentration. However, when forestland areas are suitably connected, this land use type can play a notable role in inhibiting particulate matter concentration aggravation. This conclusion is of great relevance to urban land use planning and sustainable development.

**Keywords:** land use and land cover; PM2.5; PM10; spatiotemporal characteristics; air pollution; remote sensing; China

#### **1. Introduction**

Land not only provides an important basis of human survival but is also an essential resource for human development [1]. Land use and land cover are two different but closely related concepts. Land use refers to the dynamic and purposeful utilization of natural land resources by humans, while land cover constitutes the synthesis of natural and humanmade buildings covering the Earth's surface. Both of these concepts are important land system attributes [2,3]. All human activities are inseparable from land. These activities

**Citation:** Zhai, H.; Yao, J.; Wang, G.; Tang, X. Impact of Land Use on Atmospheric Particulate Matter Concentrations: A Case Study of the Beijing–Tianjin–Hebei Region, China. *Atmosphere* **2022**, *13*, 391. https:// doi.org/10.3390/atmos13030391

Academic Editors: Elena Hristova, Manousos Ioannis Manousakas, Anikó Angyal, Maria Gini and Rajasekhar Balasubramanian

Received: 9 December 2021 Accepted: 21 February 2022 Published: 26 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

have affected the distribution patterns and change mechanisms of the original land system and exerted far-reaching influences on resources, the environment, geographical processes, and biodiversity [4,5]. The urban environment is the result of the combined influence of industrial production and urban construction. Land use types have increasingly changed from natural types, such as farmland, forestland, and grassland, to artificial types, which is one of the important characteristics of urbanization. The increasing frequency of human activities has accelerated land use type changes, altered the original urban landform and regional energy balance, damaged the natural ecological functions of the surface system, and driven air flow exchange variations between surface particles and the atmosphere, affecting the generation, diffusion, dilution, and collection of particles and thus affecting the atmospheric environment [6,7]. The compositional structure, content, and properties of the atmosphere can change due to changes in land use/land cover, and the chemical properties and processes of the atmosphere can also be affected, which can impact the gas-generation mechanism, thus influencing the circulation of atmospheric materials within a given region. In the long term, these effects may alter the local climate and even affect the global climate [8]. Landscape patterns describe the spatial distribution of landscape patches of different shapes and sizes, which is an important manifestation of landscape heterogeneity and can reflect the effects of various ecological processes. The study of landscape patterns can elucidate the inherent spatial distribution of seemingly disorderly patches and quantitatively describe their changes over time [9]. Analyzing changes in land use/land cover types and the evolution of landscape patterns can provide insights into the extent, directions, and characteristics of changes in particulate matter [10]. Environmental protection, environmental improvement, and the optimization of land use patterns in development to improve urban air quality are all major issues related to the effects of land use/land cover on atmospheric particulate matter.

In recent years, the contribution of human activity-related emissions to particulate matter has gradually increased. Vehicle exhaust, road dust, industrial emissions, domestic emissions, fuel combustion-related emissions and building dust continuously drive particulate matter emissions into the atmosphere, and these activities mainly occur in construction land areas. Buildings in cities also have a significant impact on particulate matter dispersion [11,12]. Many research results revealed that the higher the construction land proportion is, the higher the air pollution degree. Xu et al. [8] found that the construction land and road area proportion in the Changsha–Zhuzhou–Xiangtan urban agglomeration was significantly positively correlated with the concentrations of NO<sup>2</sup> and particles with aerodynamic diameters of 2.5 µm or smaller (PM2.5), and the influence on NO<sup>2</sup> was greater than that on PM2.5. In contrast, the influence on particles with an aerodynamic diameter of 10 µm or smaller (PM10) was unstable, which may have been because vehicle exhaust emissions comprised the main NO<sup>2</sup> source and were closely related to construction land and road areas. Peng et al. [13] determined that forestland in Chengdu exerted an obvious influence on NO<sup>2</sup> concentration variations; urban land was highly correlated with the highest total suspended particulate matter concentration; and the spatial distributions of NO<sup>2</sup> and total suspended particulate matter concentrations were highly consistent with those of urban land, industrial and mining land, and transportation land areas. Font et al. [14] reported that PM<sup>10</sup> greatly increased during and after road expansion in London, but the impact on PM2.5 was relatively limited, indicating that construction activities contributed more notably to PM10. Zhai et al. [15] observed a highly positive correlation between the NO<sup>2</sup> column concentration in the troposphere and the coverage rates of impervious surfaces in various cities and municipal districts, and the spatial distribution patterns in these two hot-spot areas were consistent. The mean NO<sup>2</sup> column concentration in the troposphere above open, compact, intensive, and highly intensive municipal districts increased sequentially, and the standardized concentration index of the NO<sup>2</sup> column in the troposphere above areas exhibiting impervious surface expansion exhibited an upward trend. Tang et al. [10] found that hazy days in Beijing showed obvious positive correlations with construction land, residential areas, industrial and mining land, and transportation

land. There was a significant positive correlation between the building area and PM1.0 concentration in the 0.5- and 1-km buffer zones. Wei et al. [16] determined that an obvious negative correlation existed between urban land use and pollution-free weather; however, positive correlations with other weather conditions were found, and these correlations gradually increased. Mo et al. [17] reported that the concentrations of PM2.5, PM10, and other particulate matter types in Beijing were positively correlated with the resident population density, regional gross domestic product (GDP), and other factors that significantly reflected the urbanization level.

Forestland and grassland, as typical vegetation types, are generally considered to reduce particulate matter. First, human activities in forestland and grassland are limited; few anthropogenic emission sources occur; and natural emissions are the main source of particulate matter. Second, vegetation cover can effectively separate the surface and atmosphere, thus preventing floating dust near the surface from rising into the atmosphere. Third, vegetation with a complex canopy structure can reduce the wind speed and prevent particles from entering local areas. In addition, particles fall onto vegetation surfaces due to dry sedimentation and are adsorbed or captured by leaves, stems, and other organs, which can play notable roles in dust retention. Xu et al. [8] found that forestland, green land, and cultivated land areas were negatively correlated with NO<sup>2</sup> and PM2.5 concentrations. Wei et al. [16] reported that forestland was positively correlated with nonpolluted and slightly polluted weather conditions but negatively correlated with highly polluted weather conditions. Shi et al. [11] determined that a decrease in cultivated land and forestland areas and an increase in construction land areas contributed to a decrease in the number of foggy days but aggravated the haze occurrence degree and frequency. Mo et al. [17] observed that the trend of the concentrations of PM2.5, PM10, and other particulate matter types in Beijing was the opposite of the development trend of forest coverage, exhibiting a significant negative correlation. Tang et al. [10] found that an obvious negative correlation existed between ecological land and cultivated land areas and the number of haze days in Beijing, and an obvious negative correlation between the green space area proportion and PM1.0 concentration in 0.5- and 1-km buffer zones was identified. Sun [18] reported that the PM<sup>10</sup> concentrations associated with various land use types in the Pearl River Delta region followed the ascending order of forestland, grassland, farmland, construction land, and desert wasteland. Generally, the influence of farmland on particulate matter is bidirectional. During the crop growing season, vegetation cover causes farmland to exert a certain dust-retention effect, but the corresponding efficiency may be lower than that of forestland and grassland [19,20]. After crop harvesting, dust emitted during straw burning represents an important source of particulate matter [21,22]. Moreover, in winter, without vegetation coverage, a large amount of dust becomes exposed on the surface and is easily lifted into the atmosphere [19,20]. The influence of water on particles is complicated. On the one hand, when particles fall on the water surface due to sedimentation, these particles are not resuspended. On the other hand, the presence of water causes the temperature to decrease and the humidity to increase, which may affect the accumulation and diffusion of particulate matter [20,23]. Due to the lack of vegetation coverage, a large amount of floating dust is exposed near the surface of unused land and is easily lifted into the atmosphere under the action of wind. Generally, floating dust readily becomes the main source of particles, and the corresponding contribution to particles with a large aerodynamic equivalent diameter is relatively high.

The temporal and spatial variability of particulate matter concentrations is a complex system that is driven by a combination of influencing factors. The effects of meteorological factors, human activities, and topographic conditions on particulate matter have been studied by several scholars [24–27]. In this manuscript, we focus on the intrinsic patterns of land use, and it is essential to explore its impact on the atmospheric environment to guide land use planning scientifically. In this paper, land use distribution and changes are analyzed using land use classification data in the Beijing–Tianjin–Hebei (BTH) region in 2015–2018. The differences in different land use types on the temporal changes in atmospheric particulate matter concentrations are explored with the help of an improved sine function model. The correlation between landscape patterns and atmospheric particulate matter is analyzed using landscape ecology and statistical methods. The results of this paper are important for the prevention and control of atmospheric pollution and the sustainable development of resources and the environment and can also provide a scientific reference for the optimal allocation of land resources.

#### **2. Methods**

#### *2.1. Optimized Particle Concentration Fitting Model*

Systematic mathematical models can be used to fit concentration change patterns. Comparison of the parameters in the fitted results among different regions is a powerful tool for analyzing the changes in pollutant concentrations. The concentration exhibits a significant periodicity in annual cycles, with a peak and a trough in each cycle. To analyze the characteristics of the variation in particulate matter concentration over time quantitatively, many scholars fitted periodic quantitative mathematical function models, among which the sinusoidal function model is the most widely adopted [28–32]. The traditional sine model is based on the least-squares method and comprises linear and sine functions. The general form of the sine function model is shown in Figure 1a. The corresponding equation is defined as follows:

$$y = \mathbf{A}\_0 \mathbf{x} + \mathbf{B}\_0 + \mathbf{C}\_0 \cdot \sin\left(\frac{2\pi}{\mathbf{D}\_0} (\mathbf{x} + \mathbf{E}\_0)\right) \tag{1}$$

where *y* is the particulate matter concentration, and the first two terms on the right side of the equal sign constitute a linear function, i.e., Ao*x* + Bo, representing the interannual variation trend of the particulate matter concentration. The third term is a sine function, i.e., <sup>C</sup>o· sin 2π Do (*x* + Eo) , which is applied to describe the periodic variation characteristics of the particulate matter concentration on a monthly scale. Scholars evaluated and applied this model in monthly analyses of mean particulate matter concentrations at the country, region, and city spatial scales and confirmed that it can describe the general temporal variation patterns of particulate matter concentrations. However, due to its relatively simple function shape, the model is limited in its ability to provide detailed information on the variation in particulate matter. To capture the actual trends in changes in particulate matter concentrations more closely, describe temporal changes in monthly mean concentrations more accurately, and determine the inherent mechanisms, this paper improves the traditional sinusoidal function model in the following aspects [33]:


(4) The particulate matter concentration variation during a given cycle is often not symmetrical or uniform, and the concentration increase and decrease rates in this case are not equal. In particular, the time from one wave peak to the next wave trough is not the same as the time from one wave trough to the next wave peak. According to this feature, an oscillation term with the same period and a π/2-phase difference is added to *x* in the sine function; *x* is replaced with *x* + H<sup>t</sup> · sin( 2π Dt (*x* + Et)); and the left/right shifting of the peak or trough is controlled by the magnitude of the offset factor H<sup>t</sup> .

**Figure 1.** Schematic diagram of the traditional (**a**) and improved (**b**) sine function models.

After the above 4 improvements, the general form of the improved sine function model is schematically shown in Figure 1b. The improved function is reorganized, and an improved sine function model equation is obtained as follows:

$$y = \left(\mathbf{A}\_{\mathrm{l}}\mathbf{x}^{2} + \mathbf{B}\_{\mathrm{l}}\mathbf{x} + \mathbf{C}\_{\mathrm{l}}\right) \cdot \left(\mathbf{1} + \mathbf{D}\_{\mathrm{l}} \left(\frac{1}{2} (\mathbf{1} + \sin\left(\frac{2\pi}{\mathbf{E}\_{\mathrm{l}}} \left(\mathbf{x} + \mathbf{H}\_{\mathrm{l}} \cdot \sin\left(\frac{2\pi}{\mathbf{E}\_{\mathrm{l}}} (\mathbf{x} + \mathbf{G}\_{\mathrm{l}})\right) + \mathbf{G}\_{\mathrm{l}}\right)\right)\right)\right)^{\mathrm{F}} \tag{2}$$

where *y* is the particulate matter concentration, and the terms on the right side of the equal sign comprise a quadratic function and the improved sine function that captures uneven fluctuations. There are a total of 8 parameters from A<sup>t</sup> to H<sup>t</sup> in the formula, and the meaning and description of each parameter are shown in Table 1.





#### *2.2. Landscape Pattern Index*

In this paper, the landscape pattern changes reflected by common landscape pattern indices at the landscape and patch levels were evaluated. It should be noted that certain indices may yield different meanings at different levels and should be examined separately [34,35].

#### 2.2.1. Landscape Level

At the landscape level, the calculation methods and specific meanings of the patch density (PD), largest patch index (LPI), edge density (ED), contagion index (CONTAG), splitting index (SPLIT), Shannon's diversity index (SHDI), Simpson's diversity index (SIDI), aggregation index (AI), and eight other landscape pattern indices are explained below.

PD is equal to the number of patches in the landscape divided by the total landscape area (m<sup>2</sup> ), which is then multiplied by a scaling factor of 1,000,000 (conversion into 100 hectares), with a value range of (0, +∞). Moreover, PD is expressed as the number of patches per 100 hectares. When each unit involves an individual patch, the PD value is the largest, which indicates that based on the number of patches per unit area, landscapes of different sizes can be compared. PD can be calculated as follows:

$$\text{PD} = \frac{\text{N}}{\text{A}} \cdot (1,000,000) \tag{3}$$

where N is the number of patches in the landscape, and A is the total landscape area (m<sup>2</sup> ).

LPI is expressed as a percentage and is equal to the percentage coverage of the largest patch in the area, with a value range of (0, 100). When the largest patch in the landscape decreases in coverage, the LPI gradually approaches 0. When the whole landscape comprises only one patch, i.e., the largest patch accounts for 100% of the landscape area, LPI = 100. This index is a simple way to measure dominance. The LPI equation is

$$\text{LPI} = \frac{\max(\mathbf{a\_{ij}})}{\mathbf{A}} \cdot (100) \tag{4}$$

where aij is the area of patch ij (m<sup>2</sup> ), and A is the total landscape area (m<sup>2</sup> ).

ED is equal to the sum of the lengths (m) of all edge segments related to the corresponding patch types divided by the total landscape area (m<sup>2</sup> ), which is then multiplied by 10,000 (conversion into hectares), with a value range of (0, +∞) in meters per hectare. When no graded edges occur in the landscape, i.e., when the whole landscape comprises only one patch, ED = 0. This index can be adopted to compare landscapes of different sizes. ED can be calculated as

$$\text{ED} = \frac{\sum\_{\mathbf{k}=1}^{\text{m}} \mathbf{e}\_{\mathbf{ik}}}{\mathbf{A}} \cdot (10,000) \tag{5}$$

where eik is the total length (m) of the edge segments in the landscape of patch type i, and A is the total landscape area (m<sup>2</sup> ).

CONTAG represents the highest possible contagion degree observed for a given number of patch types, with a value range of (0, 100) expressed as a percentage. When the patch types are decomposed and dispersed to the greatest extent, CONTAG approaches 0. CONTAG = 100 when all patch types exhibit the maximum aggregation degree. CONTAG reflects the aggregation or expansion degrees of different patches in the landscape. A high value indicates a good relationship between certain major patch types, while a low value indicates that the landscape encompasses a multielement intensive pattern with a high degree of landscape fragmentation. The CONTAG calculation equation is

$$\text{CONTAG} = \left[ 1 + \frac{\sum\_{\text{i=1}}^{\text{m}} \sum\_{\text{k=1}}^{\text{m}} \left( \mathbf{P}\_{\text{i}} \cdot \frac{\mathbf{g}\_{\text{ik}}}{\sum\_{\text{k=1}}^{\text{m}} \mathbf{g}\_{\text{ik}}} \right) \cdot \left( \ln \left( \mathbf{P}\_{\text{i}} \cdot \frac{\mathbf{g}\_{\text{ik}}}{\sum\_{\text{k=1}}^{\text{m}} \mathbf{g}\_{\text{ik}}} \right) \right)}{2 \ln(\text{m})} \right] \cdot (100) \tag{6}$$

where P<sup>i</sup> is the proportion of the landscape occupied by patch type I; gik is the number of connections between pixels of patch types i and k based on the double-counting method; and m is the number of patch types in the landscape.

SPLIT is equal to the square of the total landscape area (m<sup>2</sup> ) divided by the sum of the squared patch areas (m<sup>2</sup> ), and its value range is [1, n<sup>2</sup> ], where n is the total number of pixels in the area (dimensionless). When the landscape comprises a single patch, SPLIT = 1. SPLIT reaches its maximum value when the landscape is subdivided to the greatest extent. SPLIT is based on the regional distribution of accumulated patches and can characterize the degrees of separation between various patch types. SPLIT can be calculated as follows:

$$\text{SPLIT} = \frac{\text{A}^2}{\sum\_{1=1}^n \sum\_{\text{j=1}}^n \text{a}\_{\text{ij}}^2} \tag{7}$$

where aij is the area of patch ij (m<sup>2</sup> ), and A is the total landscape area (m<sup>2</sup> ).

SHDI denotes the proportional abundance of each patch type among all negative patch types multiplied by the sum of the proportions, and its value range is [0, +∞). When the landscape contains only one patch type, SHDI = 0. As the proportional distribution of areas among various patch types becomes increasingly uniform, SHDI increases. SHDI can reflect landscape heterogeneity and is particularly sensitive to an unbalanced distribution of various patch types in the landscape. If the land use diversity and fragmentation degree are high, the information content in the uncertainty and the SHDI value will increase. The SHDI equation is

$$\text{SHDI} = -\sum\_{\mathbf{i}=1}^{\text{m}} (\mathbf{P}\_{\mathbf{i}} \cdot \ln \mathbf{P}\_{\mathbf{i}}) \tag{8}$$

where P<sup>i</sup> is the proportion of the landscape occupied by patch type i.

SIDI is equal to 1 minus the sum of the proportional abundance of the square of each patch type among all patch types, with a unitless range of [0, 1]. When the landscape contains only one patch type, SIDI = 0. With an increasing number of different patch types (i.e., increasing patch abundance), SIDI approaches 1. SIDI is another popular diversity index borrowed from community ecology that is insensitive to the presence of rare types and provides a more intuitive explanation than does the Shannon index. SIDI can be calculated as follows:

$$\text{SIDI} = 1 - \sum\_{i=1}^{m} \mathbf{P}\_i^2 \tag{9}$$

where P<sup>i</sup> is the proportion of the landscape occupied by patch type i.

AI denotes the number of similar adjacencies of a certain type divided by the maximum possible number of similar adjacencies of this type, which is then multiplied by the landscape proportion of this type. Subsequently, the total value for all types is multiplied by 100 (conversion into a percentage), and AI has a value range of [0, 100] expressed as a percentage. For AI = 0, the patch types are decomposed to the greatest extent. With increasing landscape aggregation, AI increases. When the landscape comprises a single patch type, AI equals 100. AI represents the frequency of different patch types occurring side by side on the map, which can reflect the patch aggregation degree. The AI equation is

$$\text{AI} = \left[ \sum\_{i=1}^{\text{m}} \left( \frac{\text{g}\_{\text{ii}}}{\max(\text{g}\_{\text{ii}})} \right) \mathbf{P}\_{\text{i}} \right] \cdot (100) \tag{10}$$

where gii is the number of similar adjacencies between pixels of type i based on the singlecounting method; max(gii) is the maximum value of gii; and P<sup>i</sup> is the landscape proportion of type i.

#### 2.2.2. Class Level

At the class level, five landscape pattern indices, including the percentage of the landscape (PLAND), PD, LPI, largest shape index (LSI), and SPLIT, were selected, and their calculation methods and specific meanings are explained below.

PLAND is the percentage of the area of the corresponding patch type in the total landscape area, with a range of (0, 100) expressed as a percentage. When the number of corresponding patch types in the landscape decreases, PLAND approaches 0. When the whole landscape comprises a single patch type, PLAND = 100. PLAND quantifies the proportional abundance of each patch type in the landscape, which is an important measure of the landscape composition in many ecological applications. The PLAND calculation equation is

$$\text{PLAND} = \text{P}\_{\text{i}} = \frac{\sum\_{\text{j=1}}^{\text{n}} \text{a}\_{\text{i}}}{\text{A}} \cdot (100) \tag{11}$$

where P<sup>i</sup> is the landscape proportion occupied by patch type i; aij is the area of patch ij (m<sup>2</sup> ); and A is the total landscape area (m<sup>2</sup> ).

PD is equal to the number of patches of the corresponding patch type divided by the total landscape area (m<sup>2</sup> ) and multiplied by 1,000,000 (converted to 100 hectares), and its value range is (0, +∞) in units of the patch number per 100 hectares. When each cell is an individual patch, PD is the largest, and the final cell size determines the maximum number of patches per unit area. PD represents the number of patches per unit area, which is helpful for comparing landscapes of different sizes. The formula of PD is

$$\text{PD} = \frac{\mathbf{n}\_{\text{i}}}{\text{A}} \cdot (1,000,000) \tag{12}$$

where n<sup>i</sup> is the number of patches i of the corresponding type, and A is the total landscape area (m<sup>2</sup> ).

LPI denotes the percentage coverage of the largest patch of the corresponding type in the total area, with a value range of (0, 100) expressed as a percentage. When the largest patch area of the corresponding type decreases, the LPI gradually approaches 0. When the whole landscape comprises individual patches of corresponding patch types, LPI = 100. This is a simple method to measure the dominant position. LPI can be calculated as

$$\text{LPI} = \frac{\max\_{\mathbf{j}=1}^{\text{max}} (\mathbf{a}\_{\mathbf{j}})}{\mathbf{A}} \cdot (100) \tag{13}$$

where aij is the area of patch ij (m<sup>2</sup> ), and A is the total landscape area (m<sup>2</sup> ).

LSI is equal to 0.25 multiplied by the sum of the whole landscape boundary and the lengths (m) of all edge segments of the corresponding patch type in the landscape boundary, which is then divided by the square root (m<sup>2</sup> ) of the total landscape area. The LSI value range is [1, +∞] and is unitless. When the landscape contains a single square patch of the corresponding type, LSI = 1. With increasing edge lengths of the corresponding patch type in the landscape, LSI increases infinitely, which provides a standardized measure of the total edges or ED. The LSI calculation equation is

$$\text{LSI} = \frac{0.25 \cdot \sum\_{\mathbf{k}=1}^{m} \mathbf{e}\_{\mathbf{ik}}^{\*}}{\sqrt{\mathbf{A}}} \tag{14}$$

where e ∗ ik is the total length (m) of lateral edges between patch types i and k, and A is the total landscape area (m<sup>2</sup> ).

SPLIT is determined as the square of the total landscape area (m<sup>2</sup> ) divided by the sum of the squares of all patch areas (m<sup>2</sup> ), and its value range is [1, n<sup>2</sup> ], where n is the total number of pixels in the area (dimensionless). When the landscape contains a single patch type, SPLIT = 1. As the area of the focal patch type gradually decreases or the focal patch type is subdivided into smaller patches, SPLIT increases. SPLIT is based on the regional distribution of accumulated patches and can characterize the degrees of separation between various patch types. The SPLIT equation is

$$\text{SPLIT} = \frac{\text{A}^2}{\sum\_{\text{j=1}}^{\text{n}} \text{a}\_{\text{ij}}^2} \tag{15}$$

where aij is the area of patch ij (m<sup>2</sup> ), and A is the total landscape area (m<sup>2</sup> ).

#### **3. Results and Discussion**

#### *3.1. Overview of the Study Area and Data*

The BTH region is located between 113◦270~119◦510 E and 36◦050~42◦400 N (Figure 2). The terrain is inclined, exhibiting high elevations in the northwest and low elevations in the southeast, with complex and heterogeneous landforms, including plateaus, mountains, hills, basins, and plains. From northwest to southeast, the region is roughly divided into the Bashang Plateau, Yanshan-Taihang Mountains, and Haihe Plain. The BTH region experiences a temperate of a semihumid continental monsoon climate, with distinct winters and summers, and more than 67% of precipitation is concentrated in the summer. The BTH region covers an area of approximately 217,200 square kilometers, which accounts for 2.3% of the total land area in China. This area hosts a dynamic economy, the highest degree of

economic openness, the strongest innovation ability, and the highest foreign population absorption level in China. However, many particulate matter indicators seriously exceed relevant standards due to the dense population and notable industrial agglomeration, and smog and other air pollution phenomena frequently occur in the BTH region due to its rapid economic development. Air pollution seriously threatens the ecological environment and human health and affects the residential quality of life. The pollution conditions in the BTH region are highly representative of those in China and worldwide.

**Figure 2.** Geographical location of the study area.

#### 3.1.1. Land Use Data

The land use/land cover data analyzed in this experiment included the C6 version of MCD12Q1, which is a land cover product of the Moderate Resolution Imaging Spectroradiometer (MODIS, https://modis.gsfc.nasa.gov/, accessed on 1 December 2021). The temporal resolution of the data was 1 year; the spatial resolution was 500 m; and a data mosaic of four scenes could cover the entire BTH region. Since 2001, relevant data were updated at annual intervals and stored in the HDF4 file format. The projection method is sinusoidal projection, which includes 8-day comprehensive MODIS observation results obtained from Terra and Aqua satellite MODIS sensors throughout the year. Combined with the characteristics of the surface reflectivity and surface temperature, the data were analyzed with a supervised decision tree classification algorithm based on a high-quality land cover training sample database. The land cover classification schemes for the dataset include the International Geosphere-Biosphere Programme (IGBP) scheme, University of Maryland (College Park, MD, USA), leaf area index (LAI)/fraction of the absorbed photosynthetically active radiation (FPAR), net primary production (NPP), plant functional type (PFT), and Food and Agriculture Organization-Land Cover Classification System (FAO-LCCS, https://www.fao.org/home/en/, accessed on 1 December 2021) (Rome, Italy). The PFT classification scheme, including water bodies, evergreen coniferous forestland, evergreen broad-leaved forestland, deciduous broad-leaved forestland, shrub/forestland, grassland, cereal crops, broad-leaved crops, built-up areas, snow and ice, bare land, and

other types, was adopted in this experiment [36,37]. Many experiments show that the MCD12Q1 dataset has a high classification accuracy. For more information on the accuracy of the dataset products, please refer to Sulla-Menashe et al. [36,37]. Before conducting the experiment, we first reclassified the data and organized them into six types (Table 2).

**Table 2.** Land use types and descriptions after reclassification.


#### 3.1.2. Air Quality Monitoring Data

The air quality data obtained in this experiment included real-time monitoring data from monitoring stations. The China National Environmental Monitoring Center (Beijing, China) releases automatic hourly air quality monitoring results from all monitoring stations and cities in China to the public through the National Urban Air Quality Real-time Publishing Platform (http://106.37.208.233:20035/) accessed on 1 December 2021. In this experiment, 215 stations distributed in the BTH region and surrounding cities were selected (as shown in Figure 3). From 0:00 on 1 December 2014, to 23:00 on 31 December 2018, 35,808 h of hourly PM2.5 and PM<sup>10</sup> concentration data were acquired. First, the data were calculated and integrated. For the small amount of missing data, hourly average values on other days in the corresponding months were used to obtain PM2.5 and PM<sup>10</sup> values at each point. Then, the ordinary kriging interpolation method was applied to obtain the spatial distribution of particulate matter concentrations in each time interval within the study area, with a spatial resolution of 500 m. Moreover, the PM2.5/PM<sup>10</sup> ratio at each spatial position and each region was calculated. Finally, the spatial mean values of particulate matter concentrations in the whole study area in each region and for each land use type at different time scales were obtained by the statistical partitioning method. The ordinary kriging method dynamically determines the values of variables according to an optimization criterion function in the interpolation process to ensure that the interpolation function is in the best state; this process takes into account both the positional relationship between the observed and estimated points and the relative positional relationship among the observed points and is more effective when the number of points is relatively sparse.

#### 3.1.3. Accuracy Verification of the Optimization Model

To verify the accuracy of the improved model, both the traditional and improved models were adopted to fit the monthly mean values of PM2.5 and PM<sup>10</sup> concentrations in the whole study area and in cities. In addition, the goodness-of-fit R<sup>2</sup> value was compared between these two models, and the results are listed in Table 3. Among the PM2.5 and PM<sup>10</sup> fitting results, the average R<sup>2</sup> values of the traditional model are 0.51 and 0.49, respectively, and the average R<sup>2</sup> values of the improved model are 0.65 and 0.57, respectively, which indicate increases of approximately 0.14 and 0.08, respectively. The goodness of fit of the improved model is considerably higher than that of the traditional model, which verifies that the improved model achieves a higher fitting accuracy, can extract more information from the data, and is more suitable for quantitative analyses of the variation characteristics of particulate matter concentrations. In addition, compared to the traditional sine function, the improved function better quantifies the general trends, amplitudes, deformation degrees, deviation degrees, and other parameters. In addition to the particulate matter concentration, this method can be extended to other elements conforming to the periodic variation pattern involving a single peak and a single valley, such as the vegetation

coverage, monthly mean temperature, and monthly mean precipitation. The corresponding patterns of variation can be examined based on the obtained fitting results.

**Figure 3.** Spatial distribution of the air quality monitoring sites in the study area and surrounding areas.

**Table 3.** Comparison of R<sup>2</sup> between the traditional and improved models.


*3.2. Influence of Land Use on the Atmospheric Particulate Matter Concentration*

3.2.1. Spatiotemporal Characteristics of Land Use

The spatial distribution of land use in the study area for each year is shown in Figure 4a–d. The eastern and southern parts of the study area are on the Haihe Plain, and the construction land is concentrated in the built-up areas of the cities and the surrounding villages, while the agricultural land is mostly between the cities and villages. The northeastern and southwestern parts of the study area are in the Yanshan Mountain Range and Taihang

Mountain Range, and forestland is mostly distributed in areas with higher elevations and higher slopes, while grassland is distributed in areas with lower elevations and gentler slopes. The northwestern part of the study area is in the Bashang Plateau, and grassland is the most dominant land use type. The locations where land use changes occurred are shown in Figure 4e. Forestland, as the largest land use type with the largest new area, occurred mainly in the mountainous areas. The increase in grassland occurred mostly on the gentle slopes of the mountains. New farmland is mainly present in the northwestern Bashang Plateau, the eastern side of the Taihang Mountains, and the eastern coastal areas. The new construction land mainly occurs on the periphery of the built-up areas of the cities and is the spatial manifestation of urban expansion.

**Figure 4.** Spatial distributions of land use in 2015 (**a**), 2016 (**b**), 2017 (**c**), and 2018 (**d**) and a schematic representation of the land use changes that occurred (**e**).

The area of each land use type in descending order is farmland, grassland, forestland, construction land, water bodies, and unused land, among which the areas of farmland and grassland are decreasing; the areas of forestland and construction land are increasing at relatively fast and slow rates, respectively; and the area of water bodies and unused land is small and has remained relatively unchanged. The land use conversion matrix is shown in Table 4. Approximately 4.5% of the study area underwent land use change in mid-2015–2018, with grassland and forestland being the two relatively active types. A total of 8.8% of grassland was transformed, mainly to forestland and farmland. A total of 4.8% of forestland was transformed, mainly to grassland. A total of 2.6% of farmland was transformed into grassland, forestland, and construction land. Construction land is the most stable of all types and is generally not converted to other types. Smaller areas of water bodies and unused land underwent negligible conversion.


**Table 4.** Land use transfer matrix from 2015 to 2018.

3.2.2. Comparison of Atmospheric Particulate Matter Concentrations between Various Land Use Types

The spatial distributions of each land use type in the study area from 2015 to 2018 were extracted, and the extracted distributions were superimposed onto the interpolation data of the monthly mean values of the particulate matter concentration in the corresponding year and month to determine the monthly mean particulate matter concentrations corresponding to each land use type. Moreover, the average values for the same month in each year were calculated. The results are shown in Figure 5. The results indicate obvious differences in particulate matter concentrations among the various land use types. Among the land use types, the particulate matter concentrations corresponding to construction land and farmland were higher than those corresponding to the other types. Among the various land use types, construction land encompasses areas with the highest human activity intensity and motor vehicle exhaust, road dust, industrial emissions, domestic emissions, fuel combustion-related emissions, and construction dust, all of which are closely related to particulate matter emissions. In addition, the building density and height in urban built-up areas are high. The resultant blocking and friction effects weaken regional winds, greatly affecting particulate matter diffusion [11,12]. The influence of farmland on particulate matter exhibited seasonal differences. During the crop growing season, farmland is mostly covered with vegetation, which inhibits the emission of surface particulate matter to a certain extent, while leaves generate a certain dust-retention effect. However, due to the relatively low heights and comparatively neat arrangements of crops, this dust-retention effect may be lower than those of forestland and grassland [19,20]. After harvest, the dust emitted during straw burning represents an important source of particulate matter [21,22]. Especially during the cold winter season, farmland lacking vegetation coverage is mostly bare, and a large amount of dust becomes exposed and is easily emitted into the atmosphere under windy weather conditions, resulting in a notable increase in the wintertime particulate matter concentration [19,20,38]. The particulate matter concentrations corresponding to forestland and grassland were significantly lower than those corresponding to the other land use types, and the intensity of particulate matter removal was positively related to the vegetation growth degree, which was considerably higher in summer than in winter. Vegetation stems, leaves, and other organs remove particles suspended in air via retention, attachment, and adhesion mechanisms [39]. First, the plant canopy can block airflow and locally reduce the wind speed, which can result in the deposition of atmospheric particles onto the surfaces of leaves. Moreover, the turbulence between branches and leaves is very high, which makes it easier for particles to collide with and contact branches and leaves, thus increasing the settling rate [40]. Second, the surface structures of plant leaves and bark are very rough, and certain groove-like tissues, cilia, waxy layers, etc. can intercept particles that become embedded into their rough surfaces, and the deposition effect remains relatively consistent [41]. Third, plant leaves often secrete sticky substances, which can cause particle adhesion to their surfaces. This dust retention effect remains the most stable and can withstand rain erosion. Vegetation can also create an environment conducive to particle sedimentation. Tree canopies generate shade effects, and leaves generate transpiration

effects, which can achieve cooling effects. Moreover, vegetation can increase the relative air humidity and adjust the local microclimate, thus altering the viscosity and quality of particles, shortening the particle suspension time and accelerating the sedimentation process. Furthermore, vegetation effectively inhibits chemical reactions of particles and reduces the generation of secondary particles [42]. In addition, the PM<sup>10</sup> concentrations in forestland and grassland areas increased the most in spring, which is mainly attributed to frequent dusty weather conditions. The effect of water on the particulate matter concentration is complicated. On the one hand, the particulate matter emission levels of water bodies were lower than those of the other types, and water simultaneously exerted a certain adsorption effect on suspended particulate matter. On the other hand, water evaporation can increase humidity, which can affect the secondary formation, accumulation, and diffusion of particulate matter [19,43–45]. The vegetation coverage in unused land areas is low, and bare soil or sand on the surface can easily be entrained by wind, resulting in particulate pollution. However, because the unused land area in the study region is small and mostly distributed in the coastal areas and northwestern Bashang Plateau, unused land areas are mostly surrounded by water bodies and grasslands. Therefore, the particulate matter concentration mostly varied between the levels corresponding to water bodies and grasslands.

**Figure 5.** Monthly mean values of the PM2.5 and PM<sup>10</sup> concentrations corresponding to each land use type.

3.2.3. Variation Patterns of the Atmospheric Particulate Matter Concentration for Different Land Use Types

The analysis of the temporal variation in the particulate matter concentration corresponding to each land use type in the study area revealed that the particulate matter concentration corresponding to each land use type was similar to that in the whole study area, and an obvious periodic variation pattern occurred. Therefore, the improved sine function model was employed to fit the particulate matter concentration corresponding to each land use type, and the inherent relationship between the land use type and particulate matter concentration was explored by comparing the fitting parameters between various land use types. The results are shown in Figures 6 and 7.

**Figure 6.** Fitting results of the monthly PM2.5 concentration for each land use type: (**a**) water; (**b**) forestland; (**c**) grassland; (**d**) farmland; (**e**) construction land; and (**f**) unused land.

**Figure 7.** Fitting results of the monthly PM<sup>10</sup> concentration for each land use type: (**a**) water; (**b**) forestland; (**c**) grassland; (**d**) farmland; (**e**) construction land; and (**f**) unused land.

Fitting parameters were separately obtained for the whole study area and each city, and the results are summarized in Table 5.


First, the goodness-of-fit R<sup>2</sup> value was evaluated. The average R<sup>2</sup> values of the PM2.5 and PM<sup>10</sup> concentrations corresponding to the six land use types are 0.70 and 0.53, respectively, verifying that the improved model can better reflect the change pattern of the monthly mean particulate matter concentration for each land use type. The fitting R<sup>2</sup> value of the PM2.5 concentration for all land use types is above 0.6, and the fitting effect is good. The fitting R<sup>2</sup> value of the PM<sup>10</sup> concentration exceeds 0.6 for construction land; is above 0.5 for water bodies, farmland, and unused land; and is below 0.5 for forestland and grasslands. Generally, the fitting R<sup>2</sup> value of the PM2.5 concentration for all land use types is higher than that of the PM<sup>10</sup> concentration, which once again demonstrates that the periodic fluctuation patterns of the PM2.5 concentration were stronger than those of the PM<sup>10</sup> concentration. Among all types, construction land exhibits the highest R values, which is attributable to the close relationship between the seasonal patterns of human activities. Especially after the heating period in winter, the particulate matter concentration corresponding to construction land sharply increased, resulting in the largest differences between high and low values for all land use types and the most notable periodicity. The fitting R<sup>2</sup> values for forestland and grassland are relatively low. On the one hand, periodic fluctuations were slightly limited in forestland and grassland areas due to the presence of fewer particulate matter emission sources and fewer anthropogenic emissions. On the other hand, particulate matter concentrations corresponding to forestland and grassland were significantly lower than those corresponding to the other land use types, and random factors imposed a great influence, producing certain randomness and consequently reducing the periodic fluctuation patterns, thus yielding slightly poor fitting effects.

The quadratic coefficient A<sup>t</sup> and linear coefficient B<sup>t</sup> reflect the overall change trends of particulate matter concentrations in the whole study area and in cities. The PM2.5 concentrations corresponding to all land use types exhibited downward trends, with those of water bodies exhibiting an accelerated downward trend; those in forestland, grassland, and unused land all exhibiting gradual downward trends; and those in both construction land and farmland first increasing and then decreasing. Additionally, the axis of symmetry occurred in 2015, i.e., an accelerated downward trend was exhibited throughout most of the research period. In terms of the overall change trends of PM<sup>10</sup> concentrations for various land use types, the PM<sup>10</sup> concentrations corresponding to water bodies, farmland, and construction land decreased at an accelerated rate; the concentration corresponding to unused land decreased at a low rate; the concentrations in forestland and grassland areas first decreased and then increased; and the axis of symmetry occurred in 2017. However, the average PM<sup>10</sup> concentration was the lowest among all types, with a limited decline and increase; thus, the PM<sup>10</sup> concentrations basically remained at low levels. Overall, although the overall average concentrations in the study area exhibited downward trends, the change rates varied between different land use types. The particulate matter

concentrations corresponding to construction land and farmland were relatively high but simultaneously decreased at higher rates, and the decrease rates were accelerated. The particulate matter concentrations in forestland and grassland areas were low; the decrease rates of the PM2.5 concentration were decelerated; and the PM<sup>10</sup> concentrations slightly increased after decreasing to a constant level. The particulate matter concentrations in areas containing water bodies and unused land were lower than those in farmland and construction land areas but higher than those in forestland and grassland areas, with the particulate matter concentration corresponding to water bodies rapidly decreasing, while the particulate matter concentration corresponding to grassland decreased gradually.

The amplitude D<sup>t</sup> is the ratio of the focal concentration fluctuation amplitude to that of the overall concentration. Construction land and farmland exhibited significantly higher values than did the other land use types because these two types are the most affected by human activities. Therefore, these land use types exhibited large differences between winter and summer. Moreover, the amplitude of the PM2.5 concentration for all land use types was obviously larger than that of the PM<sup>10</sup> concentration, which also indicates that the PM2.5 concentration is more closely related to human sources.

The deformation factor F<sup>t</sup> reflects the durations of periods with high and low concentrations. The larger the value is, the shorter the durations of high and low concentrations are. The deformation factor F<sup>t</sup> of the PM2.5 concentration for each type was larger than 1, which suggests that the duration of low-value periods was longer than that of highvalue periods, with the F<sup>t</sup> values of construction land and farmland being significantly higher than those of the other land use types. This finding was attributable to the maximal intensity of human activities in construction land and farmland areas producing sharp increases in particulate matter concentrations in winter. This pattern was quite different from the low-value periods in summer. The relatively high values further shortened the duration of high-value periods. However, the F<sup>t</sup> values of the PM<sup>10</sup> concentration for water bodies, forestland, grassland, and unused land were smaller than 1, which indicates that the duration of high-value periods was longer than that of low-value periods for these land use types. This pattern occurred because the frequent dusty weather conditions in spring prolonged the duration of high-PM10-concentration periods, and high values continued to occur in winter and spring, which imposed a particularly notable impact on water bodies, forestland, grassland, and unused land. Farmland and construction land areas are more affected by human activities, and the F<sup>t</sup> values were larger than 1, i.e., the duration of low-value periods was increased.

The phase G<sup>t</sup> reflects the time when a peak value occurs during a period. The larger the value, the earlier the peak and valley values are observed. There was little difference in the phases of PM2.5 concentrations among different land use types, which indicates that the peak and valley values for each land use type occurred at similar times. Regarding the PM<sup>10</sup> concentration model, the phases G<sup>t</sup> of forestland and grassland were obviously lower than those of the other land use types, which was attributed to the low vegetation coverage in spring. In addition, frequent dust events slightly increased the PM<sup>10</sup> concentration for all land use types, while the spring season exhibited relatively low vegetation coverage, which exerted a notable impact on vegetation types such as forestland and grassland. Therefore, the PM<sup>10</sup> concentrations in forestland and grassland areas were higher than the notable peak values observed in March and April, while the peak values for the other land use types often occurred in December or January, resulting in higher values of the phase G<sup>t</sup> in forestland and grassland areas.

The offset factor H<sup>t</sup> controls the slight shift in the peaks and valleys from left to right. A positive value indicates a rightward shift in the valleys; a negative value indicates a leftward shift in the valleys; and the larger the absolute value is, the greater the shift. The shift factor H<sup>t</sup> of the PM2.5 concentration for each land use type was larger than 0, which indicates that for all land use types, each valley was far from the previous peak and that the next peak was closer. Notably, the concentration slowly decreased and subsequently increased rapidly. The PM<sup>10</sup> concentration indicates that the H<sup>t</sup> values in forestland and

grassland areas were smaller than or close to 0, i.e., the peak values occurred later. In addition, the concentrations decreased faster and increased slower, which is similar to the reason why the G<sup>t</sup> values of the PM<sup>10</sup> concentrations of forestland and grassland were smaller than those of the other land use types. Among all types, the H<sup>t</sup> value of farmland was the largest because farmland reduced particulate matter concentrations during the crop growth period. Hence, the particulate matter concentration exhibited a steady and gradual downward trend. However, after crops were harvested, surface dust became exposed due to the absence of vegetation coverage and was emitted into the air, and the particulate matter concentration rapidly increased.

#### *3.3. Influence of Landscape Pattern on the Atmospheric Particulate Matter Concentration*

According to the division of county-level administrative regions, the study area was divided into 200 county-level units, and the average PM2.5 and PM<sup>10</sup> concentrations throughout the whole year and during the four seasons of 2018 in each district and county were separately determined. The landscape pattern indices of each district and county at the landscape and class levels were obtained in Fragstats 4.2 software (University of Massachusetts Amherst, Amherst, MA, USA), and the influence of the landscape pattern on the atmospheric particulate matter concentration was evaluated by correlation analysis. Considering that water bodies and unused land accounted for only 1.3% and 0.2%, respectively, of the total area of the study region and imposed limited regulation effects on the large-scale particulate matter concentration, only the four most important land use types, namely, forestland, grassland, farmland, and construction land, were examined at the class level.

#### 3.3.1. Atmospheric Particulate Matter Concentration Distribution

The average PM2.5 and PM<sup>10</sup> concentrations in all districts and counties throughout the whole year and in the four seasons were statistically analyzed. The results are shown in Figure 8. The results reveal that the overall particulate matter concentration exhibited a spatial distribution trend similar to that in the whole region, in which the annual average PM2.5 concentration exhibited a trend of high values in the southeast and low values in the northwest. High values were mainly concentrated in the Shijiazhuang–Baoding and Handan–Xingtai areas, while the concentrations in Zhangjiakou, Chengde, and other cities in the north were low. During each season, the spatial distribution patterns of the seasonal average concentrations in winter and spring in all districts and counties were similar to those of the annual average concentrations, while in spring and summer, excluding the Shijiazhuang–Baoding and Handan–Xingtai areas, relatively high concentrations were observed in Beijing, Tianjin, Langfang, and other cities. The annual average trend of the PM<sup>10</sup> concentration was similar to that of the PM2.5 concentration, exhibiting the same distribution trend of high values in the southeast and low values in the northwest. The difference is that Tangshan city on the eastern coast also contained relatively high-concentration areas, with the lowest concentrations mainly concentrated in Zhangjiakou city and Chengde city in the north. The distribution pattern during each season was similar to that of the annual average concentration. Moreover, the frequent dusty weather conditions in spring caused a slight increase in the concentration in Northwest China, and simultaneously, the concentration in Tangshan increased, while relatively high-concentration areas moved slightly northwards.

**Figure 8.** Spatial distributions of the PM2.5 and PM<sup>10</sup> concentrations at the county level.

#### 3.3.2. Effects of Land Use Types on Atmospheric Particulate Matter Concentrations

The landscape pattern indices at the landscape level are shown in Figure 9, and the correlation between the particle concentration and indices at the landscape level is shown in Table 6. At the landscape level, each selected index exhibited a certain correlation with the particulate matter concentration, but the correlation degrees were lower than those at the class level, indicating that the distribution pattern of the whole landscape influenced the particulate matter concentration, but the effect degree was weaker than that observed in specific land cover types. Among the indices, LPI, AI, and CONTAG mostly exhibited positive correlations with the particulate matter concentration, while PD, ED, SPLIT, SHDI, and SIDI mostly showed negative correlations with the particulate matter concentration. The remaining indices exhibited weak correlations. The lower the aggregation degree is, the higher the fragmentation degree, the more complex the shape, and the richer the land use types are, the more closely the landscape resembles a natural landscape. Conversely, the opposite characteristics indicate an artificial landscape. Compared to artificial landscapes, natural landscapes contain fewer emission sources, and natural landscapes can simultaneously reduce and limit particulate matter concentrations. In a comparison of the seasons, the correlation of each index in winter and summer was generally higher than that in spring and autumn, which may be attributed to the high particulate matter concentration and notable spatial heterogeneity in winter. In contrast, in summer, the vegetation coverage was the highest; the reduction effect of vegetation on the particulate matter concentration was the strongest; and the landscape pattern distribution exerted a notable influence.

**Figure 9.** Spatial distributions of the landscape pattern indices at the landscape level.


\* indicates a significant correlation at the 0.05 level, and \*\* indicates a significant correlation at the 0.01 level.

The landscape pattern indices of forestland, grassland, farmland, and construction land at the class level are shown in Figure 10, and the correlation between the particle concentration and indices at the class level is shown in Table 7. PLAND, PD, LPI, and LSI all exhibited negative correlations with the particulate matter concentration, with LSI showing the highest correlation. In contrast, SPLIT exhibited no significant correlation, indicating that the larger the area, density, dominance, and shape complexity of forestland patches were, the lower the particulate matter concentration was, and the complexity of the patch shape imposed the greatest influence on the particulate matter concentration. If the forestland area within a given region is large, this land use type can become the dominant type in the region. However, if the forestland area is concentrated and connected as one unit, this area can reduce and inhibit particles. Among all seasons, the correlation in summer was slightly higher than that in the other seasons, which was related to the high vegetation coverage in forestland areas in summer and the strongest particulate matter

inhibition effect of lush trees. Moreover, the influence of forestland areas on the PM<sup>10</sup> concentration exceeded that on the PM2.5 concentration.

**Figure 10.** Spatial distributions of the landscape pattern indices at the class level.


**Table 7.** Correlations between the particle concentration and indices at the class level.

\* indicates a significant correlation at the 0.05 level, and \*\* indicates a significant correlation at the 0.01 level.

PLAND and LPI reflect the proportion, dominance degree, and patch connectivity of grassland patches, respectively, and these indices decreased spatially from northwest to southeast. High values were mainly concentrated on the Bashang Plateau and Taihang Mountains, indicating that these areas contained high grassland proportions and high dominance degrees and good connectivity of grassland patches. PD reflects the density of grassland patches, and high PD values were mostly concentrated in the Yanshan Mountains and eastern coastal areas, indicating that although the grassland patches in these areas are small, many patches are present, and their distribution is fragmented. LSI reflects the shape of grassland patches, and high LSI values largely occurred in the Yanshan Mountains and Taihang Mountains, where grassland and forestland are staggered, with complex shapes, dense patches, and small patch sizes. High SPLIT values primarily occurred in the Haihe Plain in the southeast, where the grassland coverage is low, and the distribution is scattered.

Among the various indices for farmland patches, PLAND and LPI exhibited strong positive correlations with the particulate matter concentration, while PD and LSI showed negative correlations. SPLIT exhibited no significant correlation, which indicates that the higher the farmland proportion in the region, the higher the connectivity and the more complete the patches, the higher the particulate matter concentration, and corresponding regions mainly occurred in the Haihe Plain. However, the particulate matter concentration was low in areas with high densities, scattered distributions, and irregular shapes of farmland patches. Among all seasons, PLAND, PD, and LPI exhibited the highest correlations in winter or spring, and there was less vegetation coverage in farmland patches in winter and spring. In areas where farmland dominated, bare surfaces became the main emission particulate matter source, and the resultant correlation was high. LSI showed the highest negative correlation in summer, which may have occurred because the more complex the shapes of farmland patches are, the closer the landscape is to a natural landscape. These areas were dominated by mountains, and farmlands were mostly distributed among many forestlands and grasslands. In summer, the farmland coverage intensity increased. Moreover, the forestland and grassland in the surrounding areas played a notable role in blocking and settling particles, thus reducing the particulate matter concentration.

PLAND and LPI were mostly positively correlated with the PM2.5 concentration but not significantly correlated with the PM<sup>10</sup> concentration. Notably, the PM2.5 concentration was higher in areas with high proportions of construction land, once again verifying that human activities comprise the main PM2.5 source, while PM<sup>10</sup> is more related to natural factors. PD and LSI were positively correlated with the particulate matter concentration,

indicating that the denser the construction land patches and the more complex the shape are, the higher the particulate matter concentration is. A dense distribution of construction land patches can lead to the concentration of emission sources. In contrast, a complex shape indicates that the city is highly developed, and the human activity intensity is high, which produces higher particulate emissions and does not facilitate particulate matter diffusion. SPLIT was negatively correlated with the particulate matter concentration, which indicates that in areas where the distribution of construction land patches is scattered, construction land is mostly separated by other land use types, and it is difficult to establish centralized emission sources. Moreover, the other land use types limited particulate matter concentrations, and the particulate matter concentrations were therefore low. Among the seasons, PLAND, LPI, and SPLIT exhibited the highest correlations with the PM2.5 concentration in spring, which demonstrates that the influence of construction land on the PM2.5 concentration was the strongest in spring. In the other seasons, the PM2.5 concentration was more notably influenced by other factors. PD and LSI showed the highest correlations with the particulate matter concentration in winter, which may be related to the overall high particulate matter concentration and highest spatial heterogeneity in winter.

#### **4. Implications and Limitations**

As urbanization continues to accelerate, the extent to which human activities contribute to particulate matter emissions is increasing. A large amount of harmful substances are emitted into the atmosphere, which poses a great danger to production and living activities and human health. The problem of atmospheric pollution is receiving more and more attention from relevant departments and individuals. China, especially the eastern part of the country, is one of the more serious regions in the world. Many cities have concentrations of atmospheric pollutants that far exceed the guideline values provided by the World Health Organization (Geneva, Switzerland). At the same time, with the rapid expansion of built-up areas, many cities are gradually approaching saturation in terms of space. Therefore, to control air pollution, a rational planning and layout of the limited space available is necessary. In this experiment, an attempt was made to establish the relationship between land use and particulate matter concentrations using geography, landscape ecology, and statistics. Firstly, the differences in the patterns of change in particulate matter concentrations on each land use type were analyzed, and then the relationship between landscape patterns and particulate matter concentrations was compared for each type. These analyses can provide some theoretical support for the optimization of land use structures. This can be useful for urban planning in the study area and can also provide a reference for other regions and cities. The results obtained from the experiments may be deficient due to limitations in the data sources or experimental methods. The spatial distribution of particulate matter concentrations is obtained by spatial interpolation of the Ordinary Kriging method based on the measured data at the stations. Due to the limited number of stations and the uneven spatial distribution, the interpolation results may have some errors. The experiment covers the period from 2015 to 2018, which is a limited time span. It should also be noted that the spatial and temporal characteristics of particulate matter concentrations are a complex system driven by a combination of influencing factors. Land use affects air quality both through anthropogenic emissions and agricultural activities that alter pollutant emission potential, and through interactions with the atmosphere, such as turbulence stimulation and wind speed reduction. In this manuscript, we have not explored the influence of meteorological factors in depth. All of the above will be focused on in subsequent experiments.

#### **5. Conclusions**

In this manuscript, we collated land use and atmospheric particulate matter concentration data pertaining to Beijing, Tianjin, and Hebei in China and analyzed the corresponding temporal and spatial changes and correlations. The conclusions are as follows:


**Author Contributions:** J.Y. is the corresponding author. H.Z. contributed to improving the methodology and wrote the manuscript. X.T. helped edit and improve the manuscript. G.W. contributed to testing the methodology. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was jointly supported by grants from the Application Demonstration System of GaoFen Remote Sensing Mapping of China (No. 42-Y30B04-9001-19/21); the Chinese Postdoctoral Science Foundation (No. 2021M693782).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Mahmoud Fathy ElSharkawy \* and Osama Ahmed Ibrahim**

Department of Environmental Health, College of Public Health, Imam Abdul Rahman Bin Faisal University, Dammam 34212, Saudi Arabia; olabib@iau.edu.sa

**\*** Correspondence: msharkawy@iau.edu.sa

**Abstract:** The emission of cooking fumes becomes a serious concern due to the fast development of the restaurant business because it harms the health of restaurant workers and customers and damages the outdoor air quality. This study was conducted to evaluate the impact of restaurant emissions on ambient air quality. Twenty restaurants with four different types of food cooking were selected in Dammam City, which represents a densely populated urban city in Saudi Arabia. Levels of five air pollutants were simultaneously measured in the restaurants' chimneys and in the surrounding ambient air. The highest mean levels of CO (64.8 ± 44.3 ppm), CO<sup>2</sup> (916.7 ± 463.4 ppm), VOCs (105.1 ± 61.3 ppm), NO<sup>2</sup> (4.2 ± 2.4 ppm), and SO<sup>2</sup> (8.0 ± 7.4 ppm) were recorded in chimneys of the grilling restaurants. Similarly, the highest levels of all pollutants were recorded in the areas adjacent to the grilling restaurants rather than other types.

**Keywords:** restaurants; chimney emission; combustion efficiency; ambient air pollution; emission standards

#### **1. Introduction**

Worldwide, the cooking process consumes huge amounts of energy, especially in developing countries [1,2]. Several types of fuels are usually used for cooking including natural gas, charcoal, wood, kerosene, liquefied petroleum gas, electricity, biogas, and biomass [3,4]. Consequently, large amounts of harmful air pollutants and greenhouse gases are emitted daily during the cooking processes [5,6]. Restaurants represent the most important site for cooking where many local and foreign people tend to spend a lot of their time. The emission of cooking fumes has become more serious due to the fast development of the restaurant business [7]. The emitted pollutants from restaurants not only harm the health of restaurant workers and customers but also represent a great contributor to the outdoor air pollution levels [8–12]. Complaints against cooking fume/odor emissions from restaurants have been increasing and recorded in some areas of the world [13].

For several years, great concern has been given to the cooking fumes, particularly in the highly crowded cities where restaurants are usually located in densely populated areas that are very close to residential and other sensitive buildings [13]. Emission of pollutants from restaurants results from heating and cooking operations where several types of food are cooked and different types of fuels are used [14]. The amounts and composition of pollutants emitted from those sources depend greatly on the cooking materials, cooking styles, and even cooking fuel [15]. For example, charcoal is used extensively for barbecuing in most restaurants in the world because it has high heating value, is cheap compared to other types of fuels, can be easily stored, and gives a unique flavor and texture to the food [16–18]. Charcoal contains various types of organic and inorganic compounds such as hydrocarbons, sulfur, water, and oxygen along and numerous trace elements [19–21]. Therefore, the combustion of charcoal creates a considerable amount of airborne toxic elements both in the solid and gaseous states. The coal-tars and soot (fine black particulate matter) have been documented as human carcinogens since the

**Citation:** ElSharkawy, M.F.; Ibrahim, O.A. Impact of the Restaurant Chimney Emissions on the Outdoor Air Quality. *Atmosphere* **2022**, *13*, 261. https://doi.org/10.3390/ atmos13020261

Academic Editors: Begoña Artíñano, Manousos Ioannis Manousakas, Anikó Angyal, Maria Gini and Elena Hristova

Received: 17 December 2021 Accepted: 31 January 2022 Published: 3 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

late 1700s [22,23]. Several previous studies revealed that the combustion of charcoal is considered a potential source of volatile organic compounds (VOCs) and polycyclic aromatic hydrocarbon (PAHs) [24,25] that have several adverse health effects including carcinogenicity in addition to its contribution to the formation of photochemical groundlevel ozone [26,27]. Sulfur—as a component of fuels that occurs primarily in coal, petrol, kerosene, and diesel—can produce sulfur dioxide gas (SO2) when combusted during the cooking process or any high-temperature combustion. The presence of SO<sup>2</sup> in the air leads to irritations of the mucous membranes and the eyes, as well as chronic bronchitis [28]. Combustion of charcoal is also considered a source of carbon monoxide (CO) in the air [29] that is considered a chemical asphyxiant for humans [30]. The other cooking activities in restaurants, such as charbroiling, frying, and baking, are also considered sources of the same or different air pollutants that make significant contributions to both indoor and outdoor air pollution [31,32].

Unfortunately, a lot of the exhaust outlets of ventilation ducting systems in restaurants are always not located at favorable locations, in particular in densely populated urban regions near the sensitive receptors, such as residential premises, schools, or clinics [13]. Effective control measures should be taken from the formal governmental agencies with the cooperation of owners and operators of the restaurants to ensure that no visible cooking fumes nor objectionable odor would be emitted causing any harm or forms of pollution. In this regard, appropriate high-performance air pollution control equipment must be installed at the kitchen ventilation system of the food premises for treating cooking fume emissions before being discharged to the outdoor environment as well as it is considered a cost-effective way to reduce indoor air pollution and the related health problems [33,34]. As a general guideline, the control equipment of the restaurant must be installed directly above the stoves and cooking appliances and properly connected with the exhaust ducts to prevent cooking fume from leaking through possible cracks. Moreover, the ducts must be connected with exhaust fans of adequate capacity [35]. The range hood is one of the most common types of ventilation [36]. It has the advantage of providing constant ventilation for the smoke to escape [37]. The more effective type of cooker-hood is the one that extracts the contaminated air from the cooking zone and ejects it to the ambient environment [38]. Chimneys, that must be extended to above the roof of the restaurant, are more effective because they largely prevent the smoke from entering the kitchen or any other internal site of the restaurant [39,40]. Additionally, the chimney plays an active role in the performance of the stove and in reducing emissions by influencing the overall air-to-fuel ratio and subsequently the production of CO and/or particulate matter (PM) [41]. Additionally, chimneys keep flue gas separated from ambient conditions, providing a longer residence time of the gas within a heated environment [42].

Numerous previous studies have been conducted concerning the impact of cooking emissions on the indoor environment of restaurants in developed and developing countries [43–45]. Despite their importance, data on the impact of cooking emissions on the direct surrounding environment and air pollution levels are still very scarce. This study was conducted to fill this gap by studying the impact of the restaurant emissions on the outdoor ambient air pollution in a densely populated urban city representing a developing country. It was conducted to quantify emissions of different pollutants from the chimneys of various restaurants and simultaneously levels of the same pollutants in the ambient air. It was aiming also to guide the owners and operators of restaurants, food businesses, and corresponding governmental agencies in helping them understand and apply the best practical control measures to minimize these emissions, thereby preventing air pollution problems.

#### **2. Materials and Methods**

#### *2.1. Study Area and Period*

This study was conducted in Dammam City in the Eastern province of the Kingdom of Saudi Arabia (KSA). Generally, KSA is characterized by the presence of a wide variety and large number of restaurants in all cities, particularly those with high dense populations or visitors such as Mecca, Riyadh, and Dammam because of the large numbers of immigrants, foreign workers, and people visiting the country to perform Hajj and Omrah. Dammam is considered one of the most important cities in the Kingdom. It is the capital, the major seaport, and one of the most populated cities in the Eastern Province of KSA. It is also a major administrative center for the Saudi oil industry (Aramco) and about 40% of the industrial activity of KSA is located in it. There is an increase in the migration of people to Dammam city for obtaining jobs and studying because of the presence of a large number of industries, universities, and different governmental centers. Due to the rapid population growth in the city, there is an increasing demand for food and, consequently, the number of restaurants is also increasing. All restaurants in Dammam are located close to the residential premises, schools, hospitals, and other sensitive receptors.

A variety of cooking methods are used in restaurants, but it differs from one to another according to the type of food that characterizes each restaurant. For example, the main cooking methods of some restaurants include stir frying, simmering, steaming, roasting, smoking, and stewing, while in other restaurants the main cooking methods are grilling, broiling, and deep frying. Twenty restaurants were selected in Dammam for this study representing four different types of food cooking; grilling (such as chicken or meat grilled on charcoal), frying (such as fried chicken), cooking (such as cooked rice and vegetables), and baking (such as pizza and pastry). Five restaurants from each type were selected and the twenty restaurants were contributed to the same criteria and specifications, except the type of food cooking. All selected restaurants were located on the ground floor of a residential building of 3–4 floors with a chimney extended to above the building roof. They are installed in densely populated areas adjacent to moderate traffic activity streets (about 500 cars/h) and far from any other air pollution sources such as industrial activity or any other restaurants. All restaurants have nearly the same size and number of cooking appliances. The objective of this selection was to remove any factor that could affect the results of a comparison between the emission of pollutants from the restaurants' chimneys and the ambient levels of the same pollutant. For confidence, confirmation, and comparison, a building with the same characteristics was selected in an area far from any restaurants but with the same traffic activity. This area was considered a "control area".

Climatically, Dammam has a hot desert climate. The winter temperatures range from mild to warm, while the summer temperatures are extremely hot, usually exceeding 40 ◦C (104 ◦F) for about six months. Rainfall in Dammam is generally sparse and usually occurs in small amounts in December. Heavy thunderstorms are not uncommon in winter. For this reason, this study was conducted during the six warm and hot months (April–August) of the year 2019.

#### *2.2. Measurement of Chimney Emissions*

The cooking fumes of all selected restaurants were extracted through an exhaust hood and then discharged into the atmosphere near the surrounding neighbors. Owners of restaurants did not want to make a hole in the chimney. Therefore, the sampling probe was placed directly near the outlet center of the outdoor chimney and paralleled with the direction of the chimney. Generally, for measuring the quality of the combustion of the cooking tools, the probe of a flue gas analyzer (electronic sensor) was applied to the cooking tool chimney where levels of gaseous pollutants can be measured directly on site. The advantage of the electronic sensors is that they can be used for a long time and their usage lies in real-time measurement. The air pollution content is determined using sensors where there are whole ranges of measuring principles that can be employed such as flame or photo-ionization detection for organic species, chemiluminescence for oxides of nitrogen, non-dispersive infrared for carbon monoxide, Fourier transform infrared for sulfur dioxide, etc. A digital readout indicates the measured value at the spot. These devices need to be calibrated before each monitoring session. The calibration occurs through a test gas of known pollutant concentration [42].

During our study, the Lancom 4 Portable Flue Gas Analyzer was used for measuring levels of combustion gases and cooking vapors in the exhaust chimneys of all restaurants. This analyzer meets the US EPA CTM 034 reference method, and it has data acquisition and analysis software. The analyzer is composed of a monitor, a probe hose of 3 m length, and a probe pipe of 0.3 m length. Its monitor can read up to 17 measurement parameters and it is a useful tool to observe trends. Moreover, it is free from any bias that can be caused by substances in the waste gas. The measurement specifications of this analyzer are illustrated in Table 1. The Quality Management System of Land Instruments International is approved to BS EN ISO 9001 for the design, manufacture, and on-site servicing. For quality assurance and quality control (QA/QC), the analyzer was recently calibrated by the manufacturer themselves with a certificate of conformity and calibration No. 21572853, and calibration before and after measurement using standards that are traceable to certified reference materials was conducted. Carbon dioxide (CO2), carbon monoxide (CO), nitrogen dioxide (NO2), and sulfur dioxide (SO2) were representing the combustion gases, while the volatile organic compounds (VOCs) were representing cooking vapors. Concentrations of these five pollutants were measured in parts per million (ppm). The exhaust measurements were carried out at least twice during the peak cooking period. A real-life photograph of the sampling setup while sampling was being undertaken is shown in Figure 1.

**Table 1.** Measurement specifications of the Lancom 4 Portable Flue Gas Analyzer.


\* Calibration per ASTM D-6522 or LAND factory procedure.

**Figure 1.** A real-life photograph of the sampling setup (authors personal contribution). TVOCs 0 to 10,000 ppm 0.1 ppm <25 sec <2% per mo **Figure 1.** A real-life photograph of the sampling setup (authors personal contribution).

**Table 2.** Measurement specifications of the Gray Wolf's DirectSense Gas Detector.

**Detection** 

SO2 0–20.0 ppm 0.2 ppm <25 sec <2% per mo NO2 0–20.0 ppm 0.1 ppm <20 sec <2% per mo NO 0–200 ppm 1 ppm <20 sec <2% per mo CO 0–500 ppm 1 ppm <35 sec <2% per mo CO2 0 to 10,000 ppm ±3% rdg ±50 ppm <25 sec <2% per mo

**T90 Response**  **Sensor Drift** 

**Parameter Range Limit of** 

Simultaneously with the chimney measurements, levels of the same pollutants in

distance was based on the actual presence of inhabitants' rooms on the residential building roof at 10 m from the chimney in some of the selected restaurants' buildings. The above five air pollutants were directly measured by the Gray Wolf's DirectSense® (Shelton, CT, USA) mobile PC-based products, AdvancedSense™ (Shelton, CT, USA) meters, and Wolf Pack™ (Shelton, CT, USA) area monitor. This monitor is composed of multigas detectors, and it is equipped with a wireless radio frequency modem that allows the unit to communicate and transmit readings and other information on a real-time basis with a remotely located base controller. Reliably measure key specific pollutants (VOCs, CO, O3, NO2, NH3, HCHO, etc.; choose from 25+ gas sensors), as well as particulate, ventilation rates (CO2 and airflow), differential pressure (DP), and more. High-performance, fast-response instrumentation for consistent use over portable, long-term, and continuous testing applications. In stand-alone operation, it is a rugged, weather-resistant, portable monitor that can run over 24 h on either rechargeable lithium-ion or alkaline batteries. The probe dimensions are 2 in. (5 cm) diameter × 12.5 in. (30 cm) length. Concentrations of the five pollutants were also measured in ppm. The measurement specifications of this gas detector are illustrated in Table 2. For quality assurance and quality control (QA/QC), the detector was recently calibrated by the manufacturer with a certificate of conformity and calibration No. 03-1291. Similarly, for each measuring point in the ambient around the chimney, at least three measurements for each pollutant were under-

*2.3. Measurement of Outdoor Ambient Air Pollution* 

taken for 2 h.

For each one of the four types of restaurants (grilling, frying, cooking, and baking), the measurements were conducted for inside all chimneys of each restaurant that were installed on the restaurant's building roof and simultaneously from the ambient air outside the stack. Some restaurants have only one chimney while others have two chimneys. For each chimney, at least three measurements for each pollutant were performed for 2 h. Each one of the 20 selected restaurants was visited twice, and the monitoring process was performed during the evening period because this period represents the rush hours and maximum activity of cooking and food preparation inside each type of restaurant.

#### *2.3. Measurement of Outdoor Ambient Air Pollution*

Simultaneously with the chimney measurements, levels of the same pollutants in the outdoor air were measured at 10–20 m from the chimney in the downwind direction to study the effect of emitted pollutants in the ambient air levels. The selection of this distance was based on the actual presence of inhabitants' rooms on the residential building roof at 10 m from the chimney in some of the selected restaurants' buildings. The above five air pollutants were directly measured by the Gray Wolf's DirectSense® (Shelton, CT, USA) mobile PC-based products, AdvancedSense™ (Shelton, CT, USA) meters, and Wolf Pack™ (Shelton, CT, USA) area monitor. This monitor is composed of multi-gas detectors, and it is equipped with a wireless radio frequency modem that allows the unit to communicate and transmit readings and other information on a real-time basis with a remotely located base controller. Reliably measure key specific pollutants (VOCs, CO, O3, NO2, NH3, HCHO, etc.; choose from 25+ gas sensors), as well as particulate, ventilation rates (CO<sup>2</sup> and airflow), differential pressure (DP), and more. High-performance, fast-response instrumentation for consistent use over portable, long-term, and continuous testing applications. In stand-alone operation, it is a rugged, weather-resistant, portable monitor that can run over 24 h on either rechargeable lithium-ion or alkaline batteries. The probe dimensions are 2 in. (5 cm) diameter × 12.5 in. (30 cm) length. Concentrations of the five pollutants were also measured in ppm. The measurement specifications of this gas detector are illustrated in Table 2. For quality assurance and quality control (QA/QC), the detector was recently calibrated by the manufacturer with a certificate of conformity and calibration No. 03-1291. Similarly, for each measuring point in the ambient around the chimney, at least three measurements for each pollutant were undertaken for 2 h.


**Table 2.** Measurement specifications of the Gray Wolf's DirectSense Gas Detector.

#### *2.4. Measuring of Meteorological Factors*

So far as the dispersion of pollutants from a chimney is concerned, temporal wind distribution is the most important factor for the concentration buildup of air pollutants in the surrounding air basin. The most important meteorological condition in this study was the prevailing wind at the time of measurement in the study area. Before conducting any measurements, the prevailing wind direction was recorded by the Kestrel 4500 electronic weather station (Kestrelmeters, Boothwyn, PA, USA). This tool calculates crosswind and headwind/tailwind regarding a user-set target heading and stores the information along with all the other environmental readings in its 1400 data point memory. The smoke exhaust is mixed fast with ambient air, in which high temperature and relative humidity (RH) is not a major issue in the measurements. However, wind speed, temperature, and RH were also measured by the same tool.

#### **3. Results and Discussion**

#### *3.1. Chimney Emissions*

Figure 2 represents the mean levels of the measured air pollutants in the chimney exhaust of the four methods of food cooking (grilling, frying, cooking, and baking). The highest levels of all pollutants were emitted from the grilling chimneys followed by frying and baking while the lowest levels were emitted from the cooking ones. Inside the grilling chimneys, the highest mean levels ± standard deviation (SD) of CO, CO2, VOCs, NO2, and SO<sup>2</sup> were (64.8 ± 44.3 ppm), (916.7 ± 463.4 ppm), (105.1 ± 61.3 ppm), (4.2 ± 2.4 ppm), and (8.0 ± 7.4 ppm), respectively, while in the cooking chimneys the lowest mean levels were (8.3 ± 4.4 ppm), (555.2 ± 108.7 ppm), (17.7 ± 7.1 ppm), (1.4 ± 0.6 ppm), and (1.2 ± 1.0 ppm), respectively. In most restaurants, food is prepared under high temperatures when grilled or fried whereas most fire-based cooking is based on the combustion of various fuel types (e.g., coal, natural gas, liquefied petroleum gas (LPG), and electrical energy) [3,46]. The most important chemical processes during the high-temperature treatment of food are the degradation of sugars, pyrolysis of proteins and amino acids, and the degradation of fats [44]. Several previous studies reported that the burning of charcoal is the major source of emission of air pollutants and offensive odorants in the atmosphere [16,20,47]. For example, a recent study was conducted to quantify and characterize the gaseous emissions from charcoal combustion in a brick barbecue grill revealed that emissions of CO, CO2, NOx, acid gases, NH3, and VOCs from the combustion of charcoal were higher than those of the other fuels and appliances [16]. Another study which was conducted in Portugal to assess levels of VOCs in the exhaust stacks on the roofs of a university canteen and a charcoal-grilled chicken restaurant concluded that the cooking fumes of the barbecued chicken contribute to emissions of VOCs higher than those of the university canteen [16]. Although the frying pan is different from the charcoal-burner, it can be used to heat food for high temperatures, and consequently, excess air pollutants, such as PM, CO, and VOCs, are released in the atmosphere. Numerous previous studies revealed that emission of pollutants from frying food on a hot steel pan and broiling food on steel bars above a charcoal burner was always higher than those of any other methods of cooking [12,48,49]. The results of my study are quite like most of these studies. *Atmosphere* **2022**, *13*, 261 7 of 14

Combustion efficiency is influenced by the fuel quality and the combustion chamber characteristics of a stove. A simple determination of the combustion efficiency can

value of 0.1 or lower for this ratio is a good indication of the combustion efficiency [28]. In the present study, the CO/CO2 for grilling, frying, cooking, and baking was 0.071, 0.022, 0.015, and 0.013, respectively, which reflects the good quality of combustion cham-

Applying the independent t-test for comparing statistically between means of pollutants, indicated that there is a significant difference (*p* < 0.05) between mean levels of all pollutants, except CO2, emitted from grilling and the other three methods of cooking. As for CO2, there is a significant difference (*p* < 0.05) between grilling and both frying and cooking, while there is no significant difference (*p* > 0.05) between grilling and baking. The presence of the absence of statistical differences for the other three methods differs from one pollutant to another as shown in Table 3. This means that the emission of air pollutants from the grilling process is much higher than those of the other cooking methods, and it reflects the great contribution of the grilling process in emitting air pol-

**Frying Cooking Baking** 

Grilling 0.000 \* 0.000 \* 0.005 \* Frying 0.023 \* 0.146 Cooking 0.943

Grilling 0.027 \* 0.008 \* 0.252 Frying 0.193 0.760 Cooking 0.236

VOCs Grilling 0.002 \* 0.000 \* 0.003 \*

lutants from restaurants with comparing to other types of food cooking.

**Table 3.** Independent T-test for mean levels of pollutants in chimney exhaust.

**Food Preparation Type** 

**Figure 2.** Mean levels of air pollutants in different restaurant chimneys. Frying 0.000 \* 0.004 \* **Figure 2.** Mean levels of air pollutants in different restaurant chimneys.

bers in all selected restaurants.

**Pollutant** 

CO

CO2

Combustion efficiency is influenced by the fuel quality and the combustion chamber characteristics of a stove. A simple determination of the combustion efficiency can be conducted with the calculation of the CO/CO<sup>2</sup> ratio. It is known that the mass of pollutants can be related to the mass of burnt fuel or the ratio between CO and CO2. The value of 0.1 or lower for this ratio is a good indication of the combustion efficiency [28]. In the present study, the CO/CO<sup>2</sup> for grilling, frying, cooking, and baking was 0.071, 0.022, 0.015, and 0.013, respectively, which reflects the good quality of combustion chambers in all selected restaurants.

Applying the independent t-test for comparing statistically between means of pollutants, indicated that there is a significant difference (*p* < 0.05) between mean levels of all pollutants, except CO2, emitted from grilling and the other three methods of cooking. As for CO2, there is a significant difference (*p* < 0.05) between grilling and both frying and cooking, while there is no significant difference (*p* > 0.05) between grilling and baking. The presence of the absence of statistical differences for the other three methods differs from one pollutant to another as shown in Table 3. This means that the emission of air pollutants from the grilling process is much higher than those of the other cooking methods, and it reflects the great contribution of the grilling process in emitting air pollutants from restaurants with comparing to other types of food cooking.


**Table 3.** Independent *t*-test for mean levels of pollutants in chimney exhaust.

\* The mean difference is significant at the 0.05 level.

Unfortunately, there are no emission standards for restaurants' chimneys. The emission standards promulgated in the U.S. by the EPA, Europe, and some countries of Asia are standards intended to control air pollution from several industries [50–53]. From these standards, we selected the nearest industries to restaurants such as coal-fired power plants and municipal waste combustors (MWCs) to compare the results of my study. The purpose of this comparison is a trial to set a range of safe limits for protecting people's health against restaurant fumes. Table 4 indicates the results of our study compared with the selected emission standards for only three pollutants: CO, NO2, and SO2. I did not find any standards for VOCs and CO<sup>2</sup> for the same industries. All mean levels of my study were much lower than the selected standards.


**Table 4.** Mean levels of pollutants in chimneys compared to emission standards.

#### *3.2. Outdoor Air Quality*

Figure 3 represents mean levels of air pollutants in the outdoor air at the surrounding areas of the selected chimneys, in addition to the control area. Similarly, with the chimneys results, the highest levels of CO (5.4 ± 1.4 ppm), CO<sup>2</sup> (427.1 ± 86.8 ppm), VOCs (0.31 ± 0.23 ppm), NO<sup>2</sup> (0.044 ± 0.029), and SO<sup>2</sup> (0.18 ± 0.07 ppm) were emitted from the grilling chimneys followed by frying and baking while the lowest levels were emitted from the cooking chimneys. It is shown that Figure 2 has completely the same trend as Figure 1, which indicates the direct effect of restaurant chimney exhaust in the adjacent outdoor air quality levels. It can be confirmed by the lowest levels of all pollutants that were recorded in the control area as shown in Figure 2. This means that any negative or positive change in the combustion efficiency or the internal cooking process of any restaurant will be accompanied by the same change in the outer atmosphere. No doubt, this conclusion will help the decision-makers and regulators to effectively inspect the cooking emissions from restaurants. *Atmosphere* **2022**, *13*, 261 9 of 14

**Figure 3.** Mean levels of air pollutants in the surrounding areas outdoor of restaurants. **Figure 3.** Mean levels of air pollutants in the surrounding areas outdoor of restaurants.

Statistically, the one-way ANOVA test was used to compare means of pollutants at the outdoor air for the four types of chimneys and the control area as shown in Table 5. Similarly, and surprisingly there is nearly the same statistical significance of chimneys results. Significant differences (*p* < 0.05) between mean levels of all pollutants, except CO2 and SO2, emitted from grilling were found with the other three methods of cooking. It confirms again the role of the grilling process in polluting the indoor and outdoor air. On the other hand, except for grilling restaurants, there is no significant difference (*p* > 0.05) between the control area and areas of frying, cooking, and baking restaurants. This can be explained by the considerable emission of the studied pollutants from the traffic activity. Statistically, the one-way ANOVA test was used to compare means of pollutants at the outdoor air for the four types of chimneys and the control area as shown in Table 5. Similarly, and surprisingly there is nearly the same statistical significance of chimneys results. Significant differences (*p* < 0.05) between mean levels of all pollutants, except CO<sup>2</sup> and SO2, emitted from grilling were found with the other three methods of cooking. It confirms again the role of the grilling process in polluting the indoor and outdoor air. On the other hand, except for grilling restaurants, there is no significant difference (*p* > 0.05) between the control area and areas of frying, cooking, and baking restaurants. This can be explained by the considerable emission of the studied pollutants from the traffic activity.

Mean concentrations of CO, NO2, SO2, and VOCs of the outdoor air during this study were compared with their Air Quality Guidelines (AQG) as adopted by the Saudi Environmental Law [54] and the WHO guidelines [55] and presented in Table 6. Levels Mean concentrations of CO, NO2, SO2, and VOCs of the outdoor air during this study were compared with their Air Quality Guidelines (AQG) as adopted by the Saudi

> of CO and NO2 were lower in their AQGs. Levels of VOCs were higher than their AQGs in areas of grilling and frying restaurants, while levels of SO2 exceeded their AQGs in

**Pollutant Food Preparation Type Frying Cooking Baking Control** 

Grilling 0.000 \* 0.000 \* 0.000 \* 0.000 \* Frying 0.965 0.976 0.934 Cooking 0.999 0.961 Baking .968

Grilling 0.009 0.002 \* 0.132 0.000 \* Frying 0.569 0.870 0.044 \* Cooking 0.573 0.118 Baking 0.082

Grilling 0.000 \* 0.000 \* 0.000 \* 0.000 \* Frying 0.033 \* 0.124 0.085 Cooking 0.998 0.934

sulfur in both charcoal and diesel fuel which is usually used in baking stoves.

**Table 5.** One-way ANOVA test for mean levels of pollutants in the outdoor ambient air.

CO

CO2

VOC

Environmental Law [54] and the WHO guidelines [55] and presented in Table 6. Levels of CO and NO<sup>2</sup> were lower in their AQGs. Levels of VOCs were higher than their AQGs in areas of grilling and frying restaurants, while levels of SO<sup>2</sup> exceeded their AQGs in areas of grilling and baking restaurants. This can be easily clarified by the presence of sulfur in both charcoal and diesel fuel which is usually used in baking stoves.


**Table 5.** One-way ANOVA test for mean levels of pollutants in the outdoor ambient air.

\* The mean difference is significant at the 0.05 level.



However, exact quantification of the contribution of restaurants' emissions to outdoor air is very scarce. Few previous studies were conducted to study the impact of cookstove smoke on ambient air quality. For example, a field study was conducted in four randomly selected households in two rural locations of southern Nepal during April 2017. This study revealed that 66% of particulate matter is less than 2.5 microns (PM2.5) and 80% of the black carbon emissions from biomass cookstoves directly escape into ambient air [56]. Another study was also conducted in rural Nepal revealed that a range of 6–58% of the particulate matter emitted from the open design cookstoves is liberated to the outdoor atmosphere [57].

Cooking emissions are produced from the stove used for cooking and the emissions produced by cooking the food itself. Characteristics of both stove and the food being cooked influence cooking emissions type and concentration levels. Emissions from the stove can vary significantly depending on the fuel source [58]. For example, gas burners produce higher particle concentrations, formaldehyde (HCHO), CO, and NO<sup>2</sup> when compared to electric stoves [59,60]. Solid fuel combustion in cookstoves emits a complex mixture of particulate and gaseous species, some of these pollutants contribute to levels of commonly regulated pollutants in the ambient environment [61]. Many biomass fuels and coal also contain low concentrations of chlorine that lead to low levels of emissions of dioxins and furans [62]. Besides stove and fuel source characteristics, the type of food, method of cooking, and cooking temperature can also impact the type and intensity of the cooking emission. For example, high-heat cooking activities such as broiling and frying can produce acrolein, polycyclic aromatic hydrocarbons (PAHs), and particulates, while it has also been demonstrated that the process of charbroiling and the practice of cooking fatty foods (such as high-fat hamburgers) yield higher particle emission concentrations compared to lower-heat cooking and low-fat foods [63,64].

Access to clean cooking fuels and technologies is essential for maintaining human health and achieving environmental sustainability, particularly in developing counties. A recent study has been conducted to for the first time the environmental sustainability of household cooking, focusing on remote communities in developing countries in the Southeast Asia-Pacific (SEAP) region and considering both life cycle and local impacts. To guide rural development policies, the impacts of the following cooking fuels were considered: liquefied petroleum gas, kerosene, wood, charcoal, crop residues, biogas, and electricity. Results of the study revealed that biogas from manure is environmentally the most sustainable cooking fuel, while fuelwood is the best option for climate change, with relatively low other impacts, apart from freshwater eutrophication. Cooking using electricity is the worst option since it is typically generated from diesel in off-grid communities. LPG and kerosene have higher resource depletion and land use impacts compared to biomass fuels derived from waste. Solid biomass fuels (fuelwood, charcoal, and crop residues) have high freshwater eutrophication, terrestrial ecotoxicity, and human toxicity. In addition, direct emissions from their combustion cause significant local health and environmental impacts [65].

Many intervention strategies can be used to effectively mitigate the emissions of pollutants from restaurants and protect people's health against the restaurant fumes both inside the restaurant and in its chimneys before discharging their contaminants to the ambient air. For example, separate exhaust systems must be provided to those cooking operations giving rise to oily fume and strong odor emissions and treat the emissions with separate control equipment such as venturi and activated carbon. Control equipment must be installed directly above the stoves and properly connected with the exhaust ducts to prevent cooking fume from leaking through possible cracks. For exhaust outlets near the sensitive receptors, the air pollution problem would still exist even after the application of advanced control technologies. To avoid air nuisance likely caused to the air-sensitive receivers, the owners and operators of the restaurants and food business should refrain from choosing these sites for their business. Suitable siting or positioning of the outlet of the exhaust system is of paramount importance to avoid causing or contributing to air pollution. The exhaust outlet of the restaurant chimney must be installed as high as possible for upward discharge.

#### **4. Conclusions**

The wide and fast spread of restaurants in all urban areas of the world cannot be dispensed or neglected, particularly in densely populated areas. The emission of pollutants from the restaurant chimneys has a considerable and direct effect on the outdoor ambient air, particularly the grilling process that emits pollutants at a much higher rate than those of the other food cooking methods used in restaurants. Any negative or positive change in the combustion efficiency or the internal cooking process of restaurants will be accompanied by the same change in the outer atmosphere. Fortunately, the combustion chambers and processes in all selected restaurants for this study were working efficiently, and most of the emitted pollutants were lower than their standards. The result of this study is expected to help the decision-makers and regulators to effectively inspect the emissions of pollutants from restaurants for protecting people's health against restaurant fumes and helping the restaurants' owners to take the correct actions for reducing levels of air pollution both inside the restaurant and in its chimneys before discharging their contaminants to the outer atmosphere. Furthermore, more studies must be conducted to separately study the effect of each type of fuel that is used in restaurants on the outdoor air quality.

**Author Contributions:** Methodology, M.F.E. and O.A.I.; formal analysis, M.F.E. and O.A.I.; data curation, M.F.E.; writing—original draft preparation, M.F.E.; writing—review and editing, M.F.E.; supervision, M.F.E.; project administration, M.F.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Ethical Approval:** The approval to conduct this study was obtained from the Municipality of the Eastern Province, Saudi Arabia. Each restaurant owner was provided with information about the study and the purpose of the study before conducting any sampling step. All the aspects of the subjects were kept confidential and used only for the study purpose.

#### **References**

