Spatiotemporal PAH patterns in size-fractionated particles (PM_>10-PM_0.1) from Northern Thailand biomass burning via Sentinel-2

Overviews

This study examines PM-bound PAH emissions from biomass combustion in forested and agricultural regions through satellite-based analysis. Figure 1 presents the conceptual framework, outlining the steps from data processing to emissions estimation. We applied the GSR procedure, combined with the soil-adjusted vegetation index (SAVI) and the modified normalized difference water index (MNDWI), to improve classification accuracy and reduce errors. After accurately measuring the burned area, spatial data were used to estimate PAH emissions using emission factors (EFs) specific to each of the six particle size fractions, thereby quantifying the total emissions.

Figure 1. Flowchart of the methodology. GISTDA: Geo-Informatics and Space Technology Development Agency; PAHs: polycyclic aromatic hydrocarbons.

Data collections

The GSR procedure was utilized to delineate burned areas during the dry seasons (January-May) from 2020 to 2024. We chose the B4, B8A, and B11 spectral bands of Sentinel-2 imagery to identify burned areas and reduce analysis errors. The SAVI, which is specifically designed to measure plant density in areas with sparse vegetation or very bright soil background, was also estimated based on the normalized difference vegetation index (NDVI)^[22,23], originally developed by Huete (1988)^[24], as shown in Equation (1).

The Land Development Department (LDD) of Thailand provided land use and land cover (LULC) data to support the classification of emission sources. The study focused on forested regions and agricultural lands, particularly rice, maize, and sugarcane plantations, which are major factors affecting regional air quality. Finally, due to the limited availability of local data, the EFs required to estimate PAH concentrations were derived from a comprehensive literature review.

The location of the study area

This study focuses on the upper northern part of Thailand, which includes nine provinces: Chiang Mai, Chiang Rai, Lamphun, Lampang, Nan, Phrae, Mae Hong Son, Phayao, and Uttaradit. One of the main contributing factors to air pollution in this region is open biomass burning, especially forest fires and crop residue burning, particularly from maize plantations, which is a key economic crop in the area. Maize plantations have rapidly expanded to meet the demand of the livestock feed industry in neighboring countries, resulting in forest encroachment for monoculture farming and the burning of post-harvest crop residues^[25], which is a significant source of transboundary haze pollution. Additionally, the mountainous terrain in the region facilitates the occurrence of forest fires and exacerbates air pollution due to limited air circulation in valley areas [Figure 2].

Figure 2. Study area characteristics: (A) Geographic location; (B) Topography; and (C) Land use and land cover.

Identification of burned areas

This study developed an integrated analytical workflow combining cloud-based computing, multispectral satellite imagery, and machine learning algorithms to classify burned areas across a large, diverse landscape^[26]. Google Colab was employed as the primary computational environment for data preprocessing and model execution^[27,28]. Its scalable cloud computing resources, such as graphics processing unit (GPU)/tensor processing unit (TPU) acceleration, support efficient data handling. Additionally, seamless integration with Google Drive allows the rapid implementation of analytical workflows, making it ideal for processing large remote sensing datasets^[29].

Reference burned-area information from the Geo-Informatics and Space Technology Development Agency (GISTDA), derived from 60 m Landsat-8 imagery, was used to develop the training dataset. To manage computational complexity, the study area was divided into four sections. Using the GIS Point Sampling Tool, 16,000 training samples were generated, evenly split between burned and unburned classes. The RF algorithm was selected after a comparative performance analysis with other machine learning methods in a representative sub-area, showing superior performance in assessing burned scars. The model was configured with 1,000 trees using the smileRandomForest algorithm on the Google Earth Engine (GEE) platform. For model training and validation, the dataset was randomly split into 80% training and 20% testing sets to ensure high classification reliability. The spectral bands B4, B8A, and B11 were primarily used to generate color composites for visual validation against GISTDA reference data. Furthermore, SAVI and MNDWI were implemented in Google Colab as post-classification filters rather than input variables. SAVI was applied to distinguish burned areas from high-reflectance soil or sparse vegetation, while MNDWI was used to filter out wetlands and inundated agricultural areas, particularly rice paddies, to minimize evaluation errors. This approach improved classification performance and provided a strong foundation for further analyses, including emissions estimation and environmental impact assessment [Figure 3]^[30-32].

Figure 3. Training datasets and reference burned area polygons: (A) Spatial distribution; (B) Burning and unburning points; and (C) Reference burned area polygons.

Accuracy assessment

Cohen’s kappa coefficient was used to assess the accuracy of the burned area identified from the reference data. This coefficient is widely employed to evaluate the degree of agreement between two categorical datasets beyond that expected by chance^[32]. A kappa value of 1.0 indicates perfect agreement, whereas a value of 0.0 signifies no agreement beyond chance. The detailed results are presented in Table 1.

Assessment of the air emission inventory

Air emissions from biomass burning, specifically forest fires and crop residues, were calculated using the equation derived from Giglio et al. (2006)^[33], as shown in Equation (2). Both categories were quantified using this Equation.

where EM_i,j (ton) is the emission of pollutant (i) from area (j) and M_j (kg) is the amount of burned biomass in area (j). EF_i,j is the emission factor for pollutant (i) from area (j) (g/kg of dry matter), taken from the report by Samae et al. (2020)^[34]. M_j (kg) for forest fires and crop residues was calculated using the Equations (3) and (4).

where M_j (kg) is the amount of burned biomass in forest areas, A is the burned area (km²), B is the biomass density (kg_{dry mass}/km²), and C is the burning efficiency. The values of biomass density and burning efficiency for forest fires are shown in Table 2.

where M_j (kg) is the amount of burned biomass from crop residues, A is the burned area in crop areas (km²), and B is the biomass density in crop areas (kg_{dry mass}/km²). Distinct parameters were applied for forest and agricultural areas, as detailed in Table 2. E is the burning efficiency, and the corresponding values for crop residues are also provided in Table 2.

The carcinogenicity of PAHs

PAHs are a class of organic compounds composed of multiple fused benzene rings, primarily generated during the incomplete combustion of organic materials such as biomass, coal, petroleum products, and organic waste^[40,41]. Many PAHs have been classified as carcinogenic by the International Agency for Research on Cancer (IARC)^[42,43], as shown in Table 3. Benzo[a]pyrene (BaP) has been designated as a Group 1 carcinogen (carcinogenic to humans), whereas dibenzo[a,h]anthracene (DBahA) has been classified as Group 2A (probably carcinogenic). Compounds including benzo[a]anthracene (BaA), benzo[b]fluoranthene (BbF), benzo[k]fluoranthene (BkF), and chrysene (Chr) are categorized as Group 2B (possibly carcinogenic)^[44]. According to the World Health Organization (WHO) Air Quality Guidelines, the annual average concentration of BaP in ambient air is set at 1 ng/m³, above which an increased lifetime cancer risk has been reported^[45,46]. To quantitatively assess human exposure in this study, the toxicity equivalent concentration (TEQ) was evaluated, which estimates the carcinogenic potential of the measured PAH mixture relative to BaP. During the smoke haze season, local BaP-equivalent levels (TEQ) exceeded the WHO guideline, ranging from 0.46 to 2.06 ng/m³ in rural areas, and from 0.38 to 1.33 ng/m³ in urban areas^[47]. These results demonstrate that during intensive burning episodes, carcinogenic PAH concentrations can exceed the WHO safety limit by up to 2-fold in rural environments and 1.3-fold in urban environments, confirming a high respiratory health risk for local populations^[47]. In this study, although direct ambient concentrations were not measured, the estimated emissions, particularly from forest fires, suggest that local BaP levels during burning episodes could exceed this guideline.

BaP has been classified as a Group 1 carcinogen by the IARC^[42,43], with strong evidence supporting its carcinogenicity in humans^[48-50]. It is primarily produced during the incomplete combustion of organic materials, such as cigarette smoke, vehicle exhaust, and open biomass burning. Smaller PM was found to carry higher levels of carcinogens than larger particles. The highest PAH emissions were recorded in 2023. It was found that BaP was present at 2.08 g in PM_1.0-2.5, followed by 1.21 g in PM_0.1-1.0. The lowest concentration (0.11 g) was detected in PM smaller than 0.1 µm. However, studies suggest that BaP can significantly increase the risk of lung and other cancers, even at low concentrations^[49], particularly when accumulated in PM_2.5, where it can induce oxidative stress and DNA damage, factors associated with cardiovascular, respiratory, and other diseases^[50].

Time-series changes in burned areas during the dry seasons and identification of burned areas

The integration of the GSR procedure with SAVI and MNDWI proved essential for accurate detection. In terms of spectral response, SAVI digital numbers exceed 3,000 in characterized forest areas, while values between 0 and 2,000 indicate rice plantations during the pre-cultivation phase. Applying these thresholds yielded MNDWI values surpassing 10,000, enabling effective identification of burned areas [Figure 4]. These results align with reference datasets, confirming the accuracy of the classification method.

Figure 4. DN values of SAVI and MNDWI: (A) SAVI values greater than 3,000; (B) SAVI values between 0 and 2,000; (C) MNDWI values greater than 1,000; (D) MNDWI values exceeding 1,000, differentiated by color-coded thematic classes. SAVI: Soil-adjusted vegetation index; MNDWI: modified normalized difference water index.

Analysis of the spatiotemporal dynamics reveals that the extent of burned areas peaked in 2023 (74,908.50 km²), followed by 2024 (74,606.43 km²), 2020 (72,347.33 km²), and 2022 (72,192.36 km²). The minimum extent was recorded in 2021 (56,066.36 km²). This decline was likely due to the combined effects of the COVID-19 pandemic and La Niña [Table 4]. During this period, regional meteorological data indicated significantly higher-than-average precipitation and humidity levels consistent with La Niña conditions, which naturally suppressed open burning activities. Simultaneously, socioeconomic restrictions and reduced agricultural labor mobility during the pandemic served as extrinsic factors that further limited the incidence of fire ignitions. These observations align with recent statistical analyses in the region, which report strong correlations among increased rainfall indices, economic slowdowns, and reduced aerosol loading during 2021. Overall, variability in burned area patterns was predominantly governed by meteorological anomalies, specifically the El Niño-Southern Oscillation (ENSO). The impact of the El Niño phase was pronounced during the 2023-2024 period, characterized by elevated temperatures and precipitation deficits. However, the persistence of extensive burning in 2024, despite the onset of La Niña conditions (typically associated with higher rainfall), suggests that meteorological drivers were not the sole determinants. Anthropogenic factors, particularly socioeconomic recovery following the COVID-19 pandemic, played a critical modulating role. Economic deceleration during the pandemic dampened demand for agricultural products, thereby suppressing open residue burning. Further quantitative details are provided in Table 4.

Moreover, the time series of burned area, as illustrated in Figure 5, indicates that February and March experienced the largest burned areas, influenced by meteorological conditions, particularly high temperatures and low relative humidity. In addition, terrain characteristics, such as surrounding mountains, were found to increase burning activity^[51]. These factors also hinder fire suppression efforts, allowing fires to spread rapidly. The combination of steep terrain and unfavorable meteorological conditions was found to increase combustion efficiency, leading to the release of large amounts of pollutants.

Figure 5. Monthly changes in burned area: (A) 2020; (B) 2021; (C) 2022; (D) 2023; and (E) 2024.

The spatial distribution of burned area

Figures 5-7 display the spatial distribution of burned areas in forests, as well as rice, maize, and sugarcane plantations. From these figures, the characteristics and distribution patterns of burned areas were identified. Throughout the study period, a consistent temporal variation in burned area distribution was observed. To effectively highlight this spatial pattern, data from 2023, when the highest total burned area was recorded, were selected. The resulting map indicates that geographical characteristics strongly influence the distribution of burned areas. As shown in Figure 5, forest areas exhibit recurring concentrations of burned zones across all provinces. This trend is likely related to the large proportion of forest area in northern Thailand, which accounts for approximately 63.66% of the total land area^[51].

Figure 6. Burned areas in forest regions in 2023: (A) January, (B) February, (C) March, and (D) April.

Figure 7. Burned areas in rice plantations in 2023: (A) January, (B) February, (C) March, and (D) April.

In contrast, the extent of burned areas associated with crop residues is lower than that observed in forests, particularly in rice plantations, as illustrated in Figure 6. This discrepancy is mainly due to differences in planting and harvesting schedules among agricultural systems. The dry season generally does not coincide with the rice harvesting period^[52-54], as most farmers in northern Thailand plant rice during the rainy season. Specifically, rice is typically planted between May and October and harvested between November and December. Consequently, rice straw burning usually occurs in November and December. However, a rise in burned area is observed in April, which is associated with land preparation for off-season rice cultivation in certain regions.

In maize plantation areas, burned areas are widely distributed across the region, in contrast to other agricultural regions, as shown in Figure 7. This widespread distribution can largely be attributed to the agronomic advantages of maize, particularly its ability to thrive without irrigation^[55]. As a result, this characteristic increases the likelihood of forest encroachment. In contrast, burning in sugarcane plantation areas is relatively limited, primarily because sugarcane farming is less prevalent in northern Thailand. This is partly due to the lack of sugarcane processing facilities in the region. Consequently, fires in sugarcane-growing areas are recorded at very low levels, with crop-specific economic and infrastructural factors strongly influencing agricultural burning across different farming systems.

Consistent with previous findings^[56,57], many studies in Thailand continue to rely on dNBR for burned area detection^[58-60]. However, the popular dNBR method depends heavily on the selection of appropriate pre- and post-fire imagery, limiting its applicability for monthly or high-frequency monitoring^[61,62]. In contrast, machine learning approaches such as the GSR procedure can process large, continuous datasets more effectively and provide scalable, high-resolution results.

Accuracy assessment of the activity data

The accuracy assessment of the burned area classification is illustrated in Figure 8, which compares the classification results with the reference burned area data and the GSR procedure. The visual comparison shows strong correspondence between the detected burned areas and the reference datasets in both full-scene and close-up views [Figure 9].

Figure 8. Burned areas in maize plantations in 2023: (A) January, (B) February, (C) March, and (D) April.

Figure 9. Map of the burned area as of March 2024 in Mae Hong Son: (A) Overall view; (B) Close-up of the area outlined by the black box in (A); and (C) Detailed view of the same area, including reference data and Sentinel-2 imagery.

In addition to visual interpretation, the quantitative accuracy assessment confirms strong performance. The average kappa coefficient values were 0.82 in January, 0.80 in February and March, 0.85 in April, and 0.88 in May, indicating consistently high classification reliability throughout the fire season. The accuracy assessment is further summarized in Table 5, which presents the confusion matrix and corresponding performance metrics.

The PM-bound PAHS emissions

The estimation of PM-bound PAH emissions, presented in Tables 6 and 7, was derived by integrating high-resolution burned area data generated in this study with specific EFs and combustion parameters (e.g., fuel density and combustion completeness) sourced from established literature^[34-39]. An assessment of PM-bound PAH emissions from forest fires showed that Chr was the most prevalent PAH across all size-fractionated PM, particularly in PM_1.0-2.5 and PM_2.5-10. These particle size ranges are critical because they can penetrate deeply into the human respiratory system, increasing the risk of adverse health effects^[63,64]. PM₁₀ exhibited a different trend compared with smaller particles, with naphthalene (Nap) identified as the dominant compound. This result indicates that PAH compounds attach differently to particles of varying sizes. It is important to note that PM-bound PAH emissions during the dry seasons from 2020 to 2024 showed similar trends. Consequently, the results from 2023 were used to represent the overall findings, as this year recorded the highest emissions. Detailed information on PM-bound PAH emissions from forest fires during the 2023 dry season is shown in Table 6. The results are categorized into high-molecular-weight (HMW) and low-molecular-weight (LMW) PAHs.

From Table 6, a significant presence of Chr was observed in both the PM_1.0-2.5 and PM_2.5-10 fractions. This may be attributed to condensation processes and the absorption behavior of these compounds during combustion. This dominance is primarily driven by volatility distribution and gas-particle partitioning mechanisms. As a HMW semi-volatile PAH, Chr undergoes rapid condensation from the gas phase onto the surfaces of pre-existing particles as the combustion plume cools. This process preferentially occurs on accumulation-mode particles, which provide a large surface area for adsorption and facilitate the adsorption of HMW PAHs such as Chr. In particular, PM_2.5 has a relatively large surface area, which enhances its ability to adsorb semi-volatile organic compounds^[65]. PM₁₀ was found to contain compounds with HMW and specific chemical properties. In crop residues [Table 7], Chr was also identified as the dominant PAH, although it was released in lower amounts than in forest fires. The highest PAH emissions were recorded in the PM_1.0-2.5 fraction (0.9829 g) and the PM_0.5-1.0 fraction (0.6699 g). These results indicate that crop residues, particularly rice straw, maize stalks, and sugarcane leaves, released Chr within the PM_1.0-2.5 and PM_0.5-1.0 fractions. These PM size ranges can remain airborne for extended periods and disperse over large areas, thereby increasing their potential impact on air quality and public health^[66]. However, for crop residues, the emission characteristics were consistent with those of forest fires, but the quantities were significantly lower, as shown in Table 7.

From Table 7, Chr exhibited the highest emissions across most size-fractionated PM, except for PM₁₀, where the maximum value was only 0.0111 g. In this fraction, Act was the most abundant compound (0.2629 g). This finding contrasts with forest fire emissions, where PM₁₀ typically shows the highest Nap concentration. In crop residues, however, Act was the dominant compound in this size fraction. This difference is mainly attributed to variations in biomass fuel types and combustion conditions^[67]. In forest areas, biomass materials such as wood debris and dry leaves, when burned under variable temperatures, commonly produce PAHs. Although Nap is classified as an LMW PAH, it has been found to associate with larger particles, such as PM₁₀, under specific combustion conditions^[68]. Conversely, in crop residues, combustion of different biomass types leads to higher emissions of Chr, including in larger particle sizes. In addition, BaP concentrations were influenced by landscape characteristics, with areas of dense forest cover presenting higher pollution levels^[69,70], as shown in Figure 10.

Figure 10. Spatial distribution of BaP from (A) PM_<0.1, (B) PM_0.1-1.0, and (C) PM_1.0-2.5 in 2023 following a forest fire. BaP: Benzo[a]pyrene; PM: particulate matter.

Crop residues were found to release lower levels of carcinogenic PAHs compared with forest fires. Moreover, BaP recorded a maximum emission of only 0.0822 g in the PM_1.0-2.5 fraction, compared with 2.08 g emitted by forest fires. While Chr remained the most prevalent PAH within the HMW group resulting from the combustion of crop residues, its concentration in PM_1.0-2.5 was only 0.9829 g, which was 180 times lower than that emitted by forest fires. Similarly, BbF, which showed emissions of 48.09 g from forest fires in the PM_1.0-2.5 fraction, was detected at only 0.5966 g from crop residues. The analysis indicated that forest fires were the primary source of carcinogenic PAHs, especially in PM_2.5, which can penetrate deeply into the human respiratory system and accumulate within the body. Consequently, the management and control of forest fires were identified as critical for reducing public health risks during the dry season. The spatial distribution of these emissions is presented in Figure 11.

Figure 11. Spatial distribution of BaP from (A) PM_<0.1, (B) PM_0.1-1.0, and (C) PM_1.0-2.5 in 2023 from crop residue burning. BaP: Benzo[a]pyrene; PM: particulate matter.

The analysis of PAH emissions revealed that Chr was the dominant PAH in most particle size fractions, particularly in the PM_1.0-2.5 and PM_2.5-10 ranges, whereas Nap was the dominant compound in the > 10 μm fraction. This pattern aligns with the combustion characteristics of mixed biomass in the region but differs from those reported in other contexts^[71-73]. These differences can be attributed to variations in biomass type, combustion temperature, moisture content, and chemical composition, all of which influence PAH profiles^[68,74]. Forest fires released higher amounts of PAHs across almost all compounds and particle size fractions. This is due to the larger biomass fuel load and the higher content of woody material and dry leaves, which are rich sources of HMW PAHs^[47]. Chr was identified as the dominant PAH from crop residue burning, such as rice straw and maize stalks, but was released in lower amounts due to the smaller fuel mass and lower carbon density of the biomass^[70]. In terms of health impacts, the high levels of Chr and BaP in the PM_1.0-2.5 fraction were identified as a significant concern. Their persistence in the atmosphere increases the potential for long-range transport, contributing to both local and transboundary air pollution^[75,76]. Emissions and accumulation of PAHs were strongly influenced by landform and meteorological conditions. In northern Thailand, mountainous terrain and valley structures were shown to restrict pollutant dispersion^[77]. During the dry season, low relative humidity and high temperatures were found to enhance combustion efficiency, thereby increasing pollutant emissions^[78].

To quantitatively assess human exposure, we calculated the TEQ of PAHs using BaP as a reference. During the biomass burning period, the estimated BaP-equivalent concentrations ranged from 0.86 to 2.76 ng/m³ across different size fractions. These values exceeded the WHO air quality guideline for BaP (1 ng/m³, annual average) by approximately 2-3 times. This indicates that during intense burning episodes, populations in the study area are exposed to carcinogenic PAH levels that significantly exceed international health-based thresholds, highlighting substantial inhalation risk. The results indicate that PAHs attach to particles varying from PM_>10 to PM_0.1. PM_>10 particles are larger and more likely to remain in the upper respiratory tract. In contrast, PM_2.5 and PM_0.1 particles are smaller and can get deeper into the lungs. Some PAHs, such as BaP and DahA, can bind to PM_0.1 and enter the bloodstream, where they can induce carcinogenic and mutagenic effects. Therefore, long-term exposure to PAHs attached to PM_2.5 and PM_0.1 may cause their accumulation in the lungs and bloodstream, leading to oxidative stress, inflammation, and DNA damage. These repeated exposures significantly elevate the risk of non-communicable diseases (such as lung cancer, cardiovascular disease, and COPD), particularly among vulnerable populations, including children, pregnant women, older adults, and individuals with pre-existing health conditions.