Impact of post-monsoon crop residue burning on PM_2.5 over northern India: optimizing emissions using a high-density in situ surface observation network

The impact of post-monsoon crop residue burning (CRB) on surface PM_2.5 concentrations over the Punjab–Haryana–Delhi (PHD) region in northern India was investigated using a regional meteorology–chemistry model, NHM(WRF)-Chem, and a high-density in situ surface observation network comprising Compact and Useful PM_2.5 Instrument with Gas Sensors (CUPI-G) stations. We optimized CRB emissions from 1 to 15 November 2022 using the model and surface PM_2.5 observational data. The CUPI-G data from Punjab were found to be crucial for CRB emission optimization, as the CRB emissions in northern India in October and November are predominantly originating from Punjab, accounting for 80 % of the CRB emissions. The new emission inventory is referred to as OFEv1.0, with 12 h time resolution, in daytime (05:30–17:30 IST) and nighttime (17:30–05:30 IST). The total emissions in OFEv1.0, such as PM_2.5, CO, organic carbon, and black carbon, were consistent with previous studies. OFEv1.0 substantially increased emissions relative to those calculated from satellite fire observation data (prior emissions). We showed that the prior PM_2.5 emissions were underestimated by approximately 8.6 times in the period 1–15 November 2022 and sometimes obscured completely due to clouds or thick smoke/haze on 8 and 10 November 2022. Large differences in optimized daytime and nighttime emissions indicated the importance of diurnal variations. Daytime emissions were larger than nighttime emissions on some days but not on others, indicating that diurnal variation shape may have differed each day. The mean contribution of CRB to surface PM_2.5 over PHD was 30 %–34 %, which increased to 50 %–56 % during plume events that transported pollutants from Punjab to Haryana and to Delhi. Due to insufficient performance of the meteorological model simulation on 8 and 9 November 2022 and the lack of measurement sites on the southern side of Punjab, emission optimization was not successful in the case of increased PM_2.5 concentrations observed in Haryana on these days.

1 Introduction

Delhi is a megacity with severe air pollution, and substantial efforts have been undertaken to understand the reasons underlying the rise in surface concentrations of air pollutants, source apportionments, impacts on human health, and mitigation policies and their effects (Rizwan et al., 2013; Guttikunda and Goel, 2013; Ghude et al., 2016; De Vito et al., 2018; Singh et al., 2019; Yadav et al., 2022; Lan et al., 2022; Guttikunda et al., 2023). Source apportionment studies indicate that vehicle exhaust, anthropogenic dust, biomass burning, and industry contribute approximately equally (10 %–30 %) to surface PM_2.5 concentrations in Delhi, although the dominant sector differs depending on the study (Yadav et al., 2022; Guttikunda et al., 2023). Delhi's air quality worsens during the post-monsoon to winter period, which is associated with (1) weaker wind speeds than other times of the year and increased emissions from space heaters (Guttikunda and Gurjar, 2012; Chowdhury et al., 2017, 2019; Guttikunda et al., 2023), (2) use of fireworks during Diwali festivities (Singh et al., 2019), and (3) crop residue burning (CRB) upwind of Delhi (Cusworth et al., 2018; Beig et al., 2020; Takigawa et al., 2020; Liu et al., 2020; Singh et al., 2023; Hayashida, 2023; Mangaraj et al., 2025). Delhi's air quality is influenced by CRB of the kharif crop, which is grown in the monsoon season (June–September) and harvested in the post-monsoon season (October–November) in Punjab and Haryana states under the prevailing northwesterly winds.

The post-monsoon CRB emissions have increased since implementation of the groundwater conservation policy in 2009, which delays planting of kharif crops and thus harvesting, leaving farmers with insufficient time to remove residues before planting rabi crops (grown in winter), resulting in farmers burning their stubble (Balwinder-Singh et al., 2019; Hayashida, 2023; Mukherjee et al., 2023). However, a consensus has yet to be reached on the impact of post-monsoon CRB emissions on regional air quality. For example, Cusworth et al. (2018) reported that CRB contributed 7.0 %–78 % of PM_2.5 primary components in Delhi depending on the year and selected emission inventory. Wiedinmyer et al. (2023) summarized commonly used emission inventories of open biomass burning such as the Fire Inventory from National Center for Atmospheric Research (FINNv2.5; Wiedinmyer et al., 2023), the Global Fire Emissions Database (GFED4.0s; van der Werf et al., 2017), Fire Energetics and Emissions Research (FEER; Ichoku and Ellison, 2014), the Global Fire Assimilation System (GFASv1.2; Kaiser et al., 2012), and the Quick Fire Emissions Dataset (QFED v2.5; Darmenov and da Silva, 2015). They reported that there were substantial variations in the estimation of open biomass burning emissions among them over the south and southeast Asian region, and the relative magnitudes also varied substantially among species (see Fig. 4 of Wiedinmyer et al., 2023). They concluded that determining the cause of different fire emissions in the region is a target for their future research.

Accuracy in emission inventories is crucial for the better prediction of source apportionment using 3D chemical transport models. Constructing emission inventories can take a bottom-up or top-down approach. Bottom-up inventories are based on the amount of fuel (activity data) multiplied by the emissions of chemical compounds per unit mass of fuel (emission factor). Active fire data or burned area products observed by polar-orbiting satellites are commonly employed as the activity data for bottom-up inventories of open biomass burning (e.g., Kaiser et al., 2012; Giglio et al., 2013; van der Werf et al., 2017; Beig et al., 2020; Liu et al., 2020; Singh et al., 2020; Wiedinmyer et al., 2023). Although substantial efforts have been undertaken to improve inventory data, large differences exist among emission inventories (e.g., Cusworth et al., 2018; Wiedinmyer et al., 2023). In the case of post-monsoon CRB in northern India, two main issues need to be resolved concerning emission inventories: (1) fire counts observed by satellites can be underestimated due to the presence of clouds or thick smoke/haze, and (2) determining diurnal variations is impossible because polar-orbiting satellites travel twice during the day and the night, usually around noon and midnight (Takigawa et al., 2020; Liu et al., 2020). In fact, Takigawa et al. (2020) demonstrated that 3D dispersion simulations employing emission data based on MODerate resolution Imaging Spectroradiometer (MODIS) fire radiative power (FRP) did not reproduce high PM_2.5 concentration episodes observed in Delhi in late October and early November 2019. To improve post-monsoon CRB emission estimations in northern India, Singh et al. (2020), Beig et al. (2020), and Liu et al. (2020) considered small fires utilizing Visible Infrared Imaging Radiometer Suite (VIIRS) data, with a finer spatial resolution (375 m) than MODIS (1 km) data, which resulted in a 109 % increase in FRP compared to using MODIS data alone in 2017 (Liu et al., 2020). Beig et al. (2020) combined geostationary satellite data with lower resolution (4 km) but higher temporal variation (15–30 min) to overcome the diurnal variation issue. Liu et al. (2020) tried to overcome both issues by (1) assuming a Gaussian distribution for spatial and temporal variations of CRB emissions and (2) employing a household survey of > 2000 farmers.

Another way to improve emission estimations is through a top-down approach, which minimizes the cost function between simulations and observations by adjusting for emission fluxes (e.g., Maki et al., 2010, 2011; Yumimoto and Takemura, 2015). In many cases, satellite data are employed to constrain emission fluxes due to their large spatial coverage (Elguindi et al., 2020). Surface observational data, which include more chemical compounds and have higher temporal variations than satellite data, can be utilized for limited regions (i.e., only over land because surface observations are scarce over oceans) if the spatial coverage is sufficiently high. For example, Henze et al. (2009) optimized inorganic PM_2.5 precursor emissions over the United States using Interagency Monitoring of PROtected Visual Environment (IMPROVE) monitoring datasets. However, to date, no studies have developed top-down inventories of post-monsoon CRB emissions in northern India.

India hosts a nationwide network of surface air quality monitoring stations, named Continuous Ambient Air Quality Monitoring (CAAQM) stations, and data are provided by the Central Pollution Control Board (CPCB). Because CAAQM stations are mostly located in large cities, surface PM_2.5 concentration data in CRB source areas have not been available. To address this issue, Singh et al. (2023) distributed low-cost sensors for air quality monitoring, known as Compact and Useful PM_2.5 Instrument with Gas Sensors (CUPI-G), over rural and farmland areas in the Punjab–Haryana–Delhi (PHD) region, under the Aakash project (Hayashida, 2023). Together with meteorological analysis datasets, these authors successfully identified two transport events of air pollutants from Punjab to Haryana and to Delhi in 2022.

Based on the work by Singh et al. (2023), the current study aimed to quantify the impact of CRB emissions on surface PM_2.5 concentrations (i.e., the so-called source–receptor relationship of PM_2.5) over the PHD region using a 3D regional meteorology–chemistry model (Kajino et al., 2019). Additionally, we applied an emission optimization technique to develop a top-down CRB emission inventory by resolving underestimations due to clouds and thick smoke/haze and considering diurnal variations (12 h resolution, differentiating daytime and nighttime) using CUPI-G station data (Singh et al., 2023). To the best of our knowledge, this is the first study to apply a top-down approach to constrain post-monsoon CRB emissions in northern India.

The model, observation data, and the emission optimization method are described in Sect. 2.1, 2.2, and 2.3, respectively. The model validation using the observation data is presented in Sect. 3.1. The results of emission optimization are intercompared with each other, and the best estimate is referred to as Optimized Fire Emission v1.0 (OFEv1.0) in Sect. 3.2. The OFEv1.0 data are compared against other bottom-up inventories in Sect. 3.3. The reconstructed movements of polluted air masses are presented in Sect. 3.4, and the contributions of CRB to surface PM_2.5 concentrations are quantified in Sect. 3.5. Uncertainties in OFEv1.0 associated with parameters used in the optimization method and with the planetary boundary layer (PBL) simulations are discussed in Sect. 3.6. Concluding remarks and future issues are summarized in Sect. 4.

2 Methods

2.1 Numerical models and simulation settings

The Japan Meteorological Agency (JMA)'s regional meteorology–chemistry model, named Non-Hydrostatic Model (NHM)-Chem v1.0 (Kajino et al., 2019, 2021), was utilized herein. NHM-Chem v1.0 is a chemical transport model (CTM) coupled offline with NHM, the previous version of JMA's numerical weather prediction (NWP) model (Saito et al., 2006, 2007). The CTM part of offline-coupled NHM-Chem can be employed with other meteorological models (Nakata et al., 2021; Kajino et al., 2022; Sato et al., 2023a, b). In this study, the CTM part offline-coupled with the Weather Research and Forecast model (WRF v4.1.5; Skamarock et al., 2019), named NHM(WRF)-Chem, was employed in this study because NHM(WRF)-Chem exhibited a slightly superior performance compared to NHM-Chem when evaluated against observed time series of surface PM_2.5 in northern India.

Figure 1 presents the model domains utilized in the simulation. The parent domain (domain 1; D01) covered all of India, with a horizontal grid resolution of Δx = 30 km to resolve the transport phenomena associated with synoptic scale circulations. The nested domain (domain 2; D02) covered the northwestern part of India, with Δx = 6 km to resolve local transport phenomena such as mountain–valley circulation and mesoscale cloud processes. There were 32 vertical levels from the surface to 50 hPa for WRF and 40 vertical levels from the surface to 18 km for CTM using terrain-following coordinates. Two-way nesting was applied for WRF, with the initial and boundary conditions provided by the National Center for Environmental Prediction (NCEP) final operational global analysis data (available from https://rda.ucar.edu/datasets/ds083.2, last accessed: 7 May 2024) (1° × 1°, 6 hourly). The NCEP final analysis data were also utilized for grid nudging over D01. The climatological chemistry data provided by the Meteorological Research Institute Chemistry-Climate Model version 2 (MRI-CCM2; Deushi and Shibata, 2011) (TL159, Δx ∼ 1.125°, monthly) and Model of Aerosol Species in the Global Atmosphere mark-2 (MASINGAR mk-2; Tanaka et al., 2003; Tanaka and Ogi, 2017; Yumimoto et al., 2017) (TL159, Δx ∼ 1.125°, monthly) were employed for the initial and boundary concentrations of gaseous and aerosol species over D01, respectively. The boundary concentrations of gaseous and aerosol species over D02 were provided by the hourly simulation results of D01. The simulation period was from 10 October 2022 (00:00 UTC) to 15 November 2022 (00:00 UTC), with a 5 d spin-up period, resulting in an analysis period of 1 month from 15 October 2022 (00:00 UTC) to 15 November 2022 (00:00 UTC). The time interval of WRF output and CTM input/output was 1 h.

Figure 1 Parent model domain (a, Δx = 30 km) and nested model domain (b, Δx = 6 km) with terrestrial elevations. Areas below 0 m indicate the water body grids in the model. Symbols in panel (b) indicate observation sites used in the study, which were provided by the Compact and Useful PM_2.5 Instrument with Gas Sensors (CUPI-G) observation network (colored circles) and the Central Pollution Control Board (CPCB) in India (white squares). Colors in the CUPI-G sites indicate (blue) “source”, (green) “intermediate”, and (orange) “Delhi National Capital Region (NCR)” as categorized by Singh et al. (2023). The names of states in India, such as Punjab, Haryana, Uttar Pradesh, Rajasthan, and Delhi, are depicted in panel (b).

Among physics modules in WRF, we employed Morrison's double-moment cloud microphysics scheme (Morrison et al., 2009), the RRTMG-K scheme for shortwave and longwave radiation (Baek, 2017), the Mellor–Yamada–Janjic scheme for PBL turbulence (Janjic, 1994), and the Unified Noah Land–Surface Model (Chen and Dudhia, 2001) for both domains in WRF. The Grell–Freitas ensemble scheme (Grell and Freitas, 2014) was utilized for sub-grid-scale cumulus parameterization only over D01. For the aerosol and chemistry modules in the CTM of NHM-Chem, we employed the same approach as Kajino et al. (2019, 2021), adopting the 5-category nonequilibrium method for aerosol representation. Within the 5-category method, aerosols are categorized into five categories or modes, including the Aitken mode (ATK), soot-free accumulation mode (ACM), internal mixture of soot aggregates (AGR), mineral and anthropogenic dust (DU), and sea salt particles (SS), and changes in the size distribution and chemical components in each category due to emissions, secondary production, new particle formation, advection, turbulent diffusion, and dry and wet deposition processes are solved dynamically using a triple-moment modal dynamics approach (Kajino, 2011).

As demonstrated in Fig. 2, we considered anthropogenic and natural emissions for the CTM simulation. We employed the Regional Emission inventory in Asia (REAS) version 3.2.1 (monthly, Δx = 0.25°, base year = 1950–2015) (Kurokawa and Ohara, 2020) as the anthropogenic emissions for D01 and D02. We utilized 2015 values for the simulation because 2015 was the latest dataset available. The near-real-time open biomass burning emissions provided by the Global Fire Assimilation System (GFAS; daily, Δx = 0.1°, 2003 to present) version 1.2 (Kaiser et al., 2012) were utilized for the simulation of D02. Note that fire emissions were not considered in the D01 simulation because we assumed that fire emissions from outside D02 had little influence on surface PM_2.5 concentrations over the targeted region during the simulation period. For natural sources, we employed the scheme of Han et al. (2004) for mineral dust deflation, the scheme of Clarke et al. (2006) for sea salt emissions, and the Model of Emissions of Gases and Aerosols from Nature (MEGAN2; Guenther et al., 2006) for biogenic emissions. No volcanic SO₂ emissions were considered since there were no active volcanoes in the region. The hourly and vertical profiles of Li et al. (2017) were applied for the anthropogenic emissions of REASv3.2.1. Natural emissions, such as mineral dust, sea salt, and biogenic compounds, which vary in time, were derived from simulated hourly surface meteorological variables. No diurnal variations are considered in the open biomass burning emissions (GFAS), and their vertical distribution was assumed to be uniform up to 1 km above ground level, based on the data described in Tang et al. (2022). The emissions of NO_x, SO₂, NH₃, nonmethane volatile organic compounds (NMVOCs), black carbon (BC), primary organic carbon (OC), primary PM_2.5, and primary PM₁₀ were utilized in the simulation, as described by Kajino et al. (2021). However, only carbon monoxide (CO) data are depicted in Fig. 2 to illustrate the spatial variations of anthropogenic and open biomass burning emissions in the study area.

Figure 2 (a) Anthropogenic CO emissions ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) of November 2015, provided by REASv3.2.1 over D01. (b) Same as (a) but over D02. (c) Same as (b) but the simulation period-mean open biomass burning CO, provided by the Global Fire Assimilation System (GFAS). The black arrows in the center and right panels indicate Delhi. The pink areas in the right panel indicate the five regions for area tags, North Punjab (NP), West Punjab (WP), Central Punjab (CP), East Punjab (EP), and South Punjab and North Haryana (SP&NH).

2.2 Observational data

2.2.1 Monitoring data provided by CPCB

CPCB provides near-real-time surface monitoring data of air pollutants and meteorological variables from CAAQM stations. There are currently 542 stations across India (https://airquality.cpcb.gov.in/ccr/#/caaqm-dashboard-all/caaqm-landing, last accessed: 7 May 2024) providing 15 min averaged data of air pollutants such as CO, SO₂, NO, NO₂, NO_x, O₃, PM_2.5, PM₁₀, benzene, toluene, xylene, ethyl benzene, m-xylene, p-xylene, methane (CH₄), NH₃, HCHO, and Hg, as well as meteorological variables such as temperature, wind speed, wind direction, relative humidity, pressure, solar radiation, and rainfall. The technical documentation can be found at https://erc.mp.gov.in/Documents/doc/Guidelines/CAAQMS_Specs_new.pdf (last accessed: 7 May 2024), and raw data are available at https://app.cpcbccr.com/ccr/#/caaqm-dashboard/caaqm-landing/ (last accessed: 7 May 2024). Former 1 h averaged data based on UTC were derived from 15 min averaged data to compare with the simulation data. We obtained the observational data available in 2022 from 40, 30, and 8 stations in Delhi, Haryana, and Punjab states, respectively, as shown in Fig. S1 and Table S1 in the Supplement.

2.2.2 High-density in situ surface observation network using CUPI-G stations

The CUPI-G field campaign was described by Singh et al. (2023), and details of the low-cost PM_2.5 sensors were presented by Nakayama et al. (2018); however, some key features are included here. A CUPI-G station consists of a low-cost PM_2.5 sensor developed by Panasonic Co., Ltd. (Osaka, Japan) (Nakayama et al., 2018) and other low-cost sensors for gases provided by AMETEK Inc. (Berwyn, PA, USA), such as the CO-B4 Carbon Monoxide Sensor, NO-B4 Nitric Oxide Sensor, NO2-A43F Nitrogen Dioxide Sensor, and the OX-A431 Oxidizing Gas Sensor (https://www.alphasens.com, last accessed: 7 May 2024). The 15 min averaged data were generated based on the 2 min raw data after a quality check, and then the former 1 h averaged data based on UTC were derived from 15 min averaged data to compare with the simulation data. There were 29 stations available in 2022. As listed in Table S1, Singh et al. (2023) categorized the CUPI-G stations into “source”, “intermediate”, and “Delhi National Capital Region (NCR)” regions, which were basically based on state boundaries, namely, Punjab, Haryana, and Delhi NCR, and some CUPI-G stations from other states. The site names, locations, and location types can be found in the supplementary material of Singh et al. (2023), and the same category definitions were adopted herein. All 14 stations in Punjab were categorized as the source region. Nine stations were categorized in the intermediate region, including eight stations from Haryana and one station from Uttar Pradesh (Muzaffarnagar). Delhi NCR included six stations, including one station in New Delhi (Jawaharlal Nehru University), two stations in Uttar Pradesh (Meerut and Aurangabad), and three stations in Haryana (two in Gurugram and one in Faridabad). Notably, the categorizations differed for CPCB and CUPI-G stations because CPCB stations were categorized on a state basis.

2.2.3 Remote sensing data, MODIS, TROPOMI, and AERONET

The Level 3 daily global 1° × 1° aerosol optical depth (AOD) at a wavelength of 550 nm and cloud fraction (CF) data of MODIS aboard the National Aeronautics and Space Administration (NASA) Terra satellite (observation time: 10:30 local time, descending mode) were employed in this study. The AOD_550_Dark_Target_Deep_Blue_Combined_Mean and the Cloud_Fraction_Mean variables of Collection 6.1 were used for the AOD and CF data, respectively. The data description is available at https://atmosphere-imager.gsfc.nasa.gov/sites/default/files/ModAtmo/L3_ATBD_C6_C61_2019_02_20.pdf (last accessed: 7 May 2024).

The 1-orbit Level 2 (5.5 km × 3.5 km) TROPOspheric Monitoring Instrument (TROPOMI) ultraviolet aerosol index (UVAI) using the 340/380 nm wavelength pair data aboard the European Space Agency (ESA) Sentinel-5 Precursor (S5P) satellite was utilized in this study. The data description is available at https://sentinel.esa.int/documents/247904/2474726/Sentinel-5P-Level-2-Product-User-Manual-Aerosol- Index-product (last accessed: 7 May 2024). Even though MODIS AOD is not retrieved in the presence of clouds, UVAI data can be available.

The ground-based measurement data of AOD from the Aerosol Robotic Network (AERONET) were employed in model validation as an independent dataset from the data utilized for emission optimization, i.e., CPCB and CUPI-G. AERONET consisted of 434 stations worldwide (two in D02) in 2022, imposing standardization of the instrument (Cimel's sun photometer, CE318 series), calibration, and data processing. We employed Level 2.0 (quality-assured) Version 3 AOD at a wavelength of 500 nm for validation of the simulated AOD at the same wavelength by assuming the Maxwell Garnett approximation for non-light-absorbing and light-absorbing internal mixtures, such as BC and mineral dust. The data and their description are found at https://aeronet.gsfc.nasa.gov/ (last accessed: 7 May 2024).

2.2.4 IAGOS

The In-service Aircraft for a Global Observing System (IAGOS) data were used to evaluate the simulated vertical profiles of atmospheric constituents. IAGOS is a European Research Infrastructure for global observations of atmospheric composition from commercial aircraft (https://www.iagos.org/, last accessed: 18 February 2025). During the analysis period, six flight profiles are available departing from and landing at Indira Gandi International (IGI) airport at Delhi, as shown in Fig. S1. CO and O₃ data were used for model evaluation in this study.

2.3 Emission optimization using tagged simulation

To optimize CRB emission fluxes, we used a simple cost function basically following the concept of the equation commonly used for Bayesian synthesis inversion (e.g., Eq. 1 in Baker et al., 2006), which is expressed as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}f& =\sum _{n=\mathrm{1}}^N{\displaystyle \frac{{\left(O_n-\left(S_{\mathrm{0},n}+{\sum }_m{x}_m\left(S_{m,n}-S_{\mathrm{0},n}\right)\right)\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_n^{\mathrm{2}}}}\\ & +\sum _{m=\mathrm{1}}^M{\displaystyle \frac{{\left(x_m-x_{\mathrm{0}}\right)}^{\mathrm{2}}}{{\left(u_m-l_m\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

where O and S are the observed and simulated PM_2.5 surface concentrations, respectively; N and M are the numbers of observational data and sensitivity tests (or tags), respectively; S_0,n is a simulation result without CRB emissions (i.e., anthropogenic and natural emissions only); S_m,n is a tagged simulation including CRB emissions of the right panel of Fig. 2 (period-mean GFAS emission) for certain times and areas; σ_n is uncertainty of the observational data; x_m is a variable to optimize the cost function with the upper (u_m) and lower (l_m) limits; and x₀ is the initial condition of x_m. Similar to Eq. (1) in Baker et al. (2006), the first and second terms on the right-hand side of the equation represent the deviations between simulations and observations and between optimized and a priori simulations, respectively. Optimization reduced both sets of deviations simultaneously. We employed the limited-memory Bryoyden–Fletcher–Goldfarb–Shanno (L-BFGS)-B algorithm (L-BFGS-B; Byrd et al., 1995; Zhu et al., 1997) to minimize the cost function. L-BFGS-B is an extension of the quasi-Newton algorithm L-BFGS (Nocedal, 1980; Liu and Nocedal, 1989), which is utilized to minimize a nonlinear function f(x) subject to l ≤ x ≤ u using the derivatives g_m = ∇f(x_m) as a key driver of the algorithm to identify the direction of steepest descent (Zhu et al., 1997). This optimizer was commonly used in previous studies of atmospheric chemistry inverse modeling (e.g., Zheng et al., 2018).

The list of sensitivity simulations performed herein is summarized in Table 1. We conducted whole-period simulations without CRB emissions (“No_CRB”) (i.e., without GFAS emissions), with GFAS emissions (“GFAS”), and using several emissions optimized by Eq. (1) (“Optimized”). Among the “Optimized” emission simulations, as described in Sect. 3.2, the best estimate was referred to as Optimized Fire Emission v1.0 (OFEv1.0). The whole-period simulations began on 15 October 2022 (00:00 UTC), using the same initial conditions simulated by the “No_CRB” test during 10–15 October 2022 (00:00 UTC), that is, the 5 d spin-up period. The tagged simulations and subsequent emission optimization were separately applied for the two plume periods defined by Singh et al. (2023), including the Plume 1 period in 2–5 November 2022 (00:00 IST; India standard time) and Plume 2 period in 8–12 November 2022 (00:00 IST). The simulation periods were 1–6 November 2022 (00:00 UTC) and 6–14 November 2022 (00:00 UTC) for Plumes 1 and 2, respectively. As shown in Fig. 2, we set five areal tags over the CRB source region, including North Punjab (NP), West Punjab (WP), Central Punjab (CP), East Punjab (EP), and South Punjab and North Haryana (SP&NH). Based on GFAS, the emission contribution from the five regions to the emissions over D02 was approximately 80 % from 15 October to 15 November 2022 (00:00 UTC). In addition to the five areal tags, 7 and 11 temporal tags were set for emission optimization for the Plume 1 and Plume 2 periods, respectively, as shown in Table 1. On each day, the temporal tags were divided into a.m. (00:00–12:00 UTC or 05:30–17:30 IST) and p.m. (12:00–24:00 UTC or 17:30–05:30 IST) to elucidate the relative abundance of daytime and nighttime fire ignition.

Table 1 Sensitivity simulations.

^a The best estimation in this study: Optimized Fire Emission v1.0 (OFEv1.0).
^b The plume periods defined by Singh et al. (2023).
^c 00:00–12:00 UTC (05:30–17:30 IST).
^d 12:00–24:00 UTC (17:30–05:30 IST).
n/a: not applicable.

In total, there were 36 and 56 tags (areal tags multiplied by temporal tags plus one; GFAS tag) for the Plume 1 and Plume 2 periods, respectively, and solving Eq. (1) in one step was found to be unsuccessful. In addition, there was no prior information regarding the upper and lower limits of x_m. Therefore, we applied multi-step optimization with smaller limit values as follows: Step 1 involved solving Eq. (1) using temporal tags only with (u_m, l_m) = (2.0, 0.5) by summing all areal tagged simulations; Step 2 involved solving Eq. (1) using areal tags only with (u_m, l_m) = (2.0, 0.5) by summing all temporal tagged simulations multiplied by the optimized x_m value of the result of Step 1; and Step 3 was the same as Step 1, but it involved summing all areal tagged simulations multiplied by the optimized x_m value of the result of Step 2. Let x_m,i be the value of x_m obtained in the ith step. This multi-step optimization was repeated until 0.009 < x_m,i < 1.001 was obtained for all m at the kth step, so that the final x_m was obtained as

$$\begin{array}{}\text{(2)}& x_m=\prod _{i=\mathrm{1}}^k{x}_{m,i}.\end{array}$$

The optimization results of Eq. (1) should be substantially altered by the selection of uncertainty of observation data σ_n. There are several ways to estimate σ_n using constant values (Maki et al., 2010) for CO₂ and constant rates (Maki et al., 2011) for p.m. In the same manner as Maki et al. (2011), we derived the uncertainty rate as 10 % from the standard deviations of 3 h mean values, so that σ_n = 0.1O_n. However, as shown later in Sect. 3.6, sensitivity tests of emission optimization using the constant uncertainty rates ranging from 1 % to 500 % and (u_m, l_m) = (2.0, 0.5), (10, 0.1), and (100, 0.01) underpredicted the emission fluxes. Therefore, in this study, constant σ_n values are selected as 20 µg m⁻³ for all data, which is 10 % of the average of observed PM_2.5 from all stations (200 µg m⁻³). Maki et al. (2010) selected constant values in time but different values in space. In this study, we used the same value for all data by assuming that the PM_2.5 values do not vary substantially in space over the region and the period. The sensitivity tests of emission optimization using the constant σ_n values ranging from 1 to 500 µg m⁻³ are also conducted as shown in Sect. 3.6. We set x₀ = 1 for all cases.

Observation site selection is essential for better estimations of CRB emissions. Herein, we only selected observation sites where data were continuously available and reliable for the whole simulation period (i.e., the data did not include sudden gaps and exhibited zero drift). Moreover, we only selected sites where the simulated values agreed well with the observed values during the period not affected by CRB, namely, 15–28 October 2022. The NHM(WRF)-Chem simulation can only predict spatiotemporal mean concentrations from steady-state emission sources, whereas observational data are affected by the local environment and sporadic emissions. If the spatial and temporal representativeness of an observation site was small or the site was largely affected by local-scale disturbances, deviations between the simulation and observations would be large. Using all observational data, optimization would reduce these deviations due to differences in spatiotemporal representativeness between the simulation and observations. In other words, if there are big gaps in the simulated and observed data during the non CRB period, our optimization would be misled to fill the gaps by altering CRB emissions. All CPCB and CUPI-G sites are listed in Table S1, as well as those selected in this study.

The optimized x_m value should vary depending on the chemical compound. Ideally, the same compound should be utilized for optimization and emission changes. For example, the optimized x_m value obtained from observed and tagged simulations for CO should be applied to optimize CO emissions. However, quality assurance and control (QA/QC) had only been performed for CUPI-G PM_2.5 data at the time of the study (Singh et al., 2023), and not all variables were available in CPCB data. Therefore, we applied the same x_m values obtained from the observed and tagged simulations of PM_2.5 for the optimized fire emissions of all primary precursor species, including NO_x, SO₂, NH₃, NMVOCs, BC, OC, PM_2.5, and PM₁₀. Thus, we assumed that the relative magnitude of the emission factors for all primary emission species was consistently estimated in GFAS emissions, and chemical reactions, transport, and deposition processes were consistently solved by NHM(WRF)-Chem; therefore, the deviations between GFAS simulations and observations only originated from discrepancies in spatiotemporal variations of satellite-derived fire detection over the region and difference in emission factors, which is common for all chemical compounds. Although this assumption was unlikely, it was the best available option.

3 Results and discussion

3.1 Evaluation of meteorological and chemical simulations with GFAS

Figure 3 presents the time series of station-mean observed and simulated PM_2.5 concentrations by “GFAS” for CPCB and CUPI-G stations over Punjab (source), Haryana (intermediate), and Delhi (Delhi NCR). The observational data differed for CPCB and CUPI-G stations on 24 October 2022, which was the date of Diwali, the festival of lights. Because CPCB stations are situated in urban locations, remarkable enhancements in PM_2.5 were observed at Punjab and Haryana CPCB stations due to festive use of firecrackers and candles, whereas no enhancements were observed at CUPI-G source and intermediate stations situated in farmland areas. Because the model did not consider sporadic emissions due to events such as Diwali, enhancements were not simulated on 24 October 2022 in these regions. Thus, optimization using CPCB data that included the Diwali day would cause misleading CRB emissions. Another remarkable feature was the large difference in observed PM_2.5 between the CPCB Punjab and CUPI-G source regions. The observed and simulated PM_2.5 for CPCB Punjab did not differ greatly, except on 8 and 9 November, whereas large underestimations were found in the simulation for the CUPI-G source region from 29 October to 10 November 2022, which was the intensive CRB period of that year. A similar feature was observed in CPCB Haryana and CUPI-G intermediate regions, where observation values at CUPI-G stations over the CRB period were larger than those observed in CPCB Haryana and those simulated by “GFAS”. These features demonstrated the successful allocation of CUPI-G stations over source and intermediate regions, thereby avoiding large anthropogenic emission sources and capturing the influence of CRB emissions.

Figure 3 Temporal variations (in UTC) of selected station-mean (black) observed and (red) “GFAS”-simulated PM_2.5 concentrations (µg m⁻³) for (a) CPCB Punjab, (b) CUPI-G source, (c) CPCB Haryana, (d) CUPI-G intermediate, (e) CPCB Delhi, and (f) CUPI-G Delhi NCR from 20 October to 15 November 2022 (00:00 UTC).

The time series differences between selected and all stations for PM_2.5 and CO are presented in Figs. S2 and S3, respectively. The importance of station selection was evident, especially for CO in Delhi (Fig. S3). The “GFAS” simulation seemed to underestimate emissions due to the underestimation of CRB emissions only for CPCB selected stations over Delhi (Fig. S3e), whereas the simulation tended to underestimate emissions almost every night for all CPCB stations over Delhi (Fig. S3f). This was likely because CPCB stations were located within the urban canopy, whereas the simulation predicted air concentrations above it. The inclusion of all CPCB station data in emission optimization would prevent the optimization of only the effect of CRB emissions. The same trends were observed in Punjab and Haryana. However, the PM_2.5 differences between selected and all stations (Fig. S2) were not overly remarkable compared to those of CO, which was likely because CO is generated by primary sources, whereas PM_2.5 consists of both primary and secondary sources. Concentrations within and above the canopy for primary species such as CO may differ substantially because emission contributions dominate air concentrations within the canopy. However, for secondary species, the contribution of transport from other regions became relatively larger, thereby reducing concentration differences within and above the canopy. The observed PM_2.5 was slightly greater in all stations than in selected stations, especially for the CUPI-G intermediate region (Fig. S2i, j); however, the differences were minor compared to those of CO.

Figure 4 presents the observed and simulated time series of selected station-mean meteorological variables (the selection of stations was based on PM_2.5 data, not CO data). Since no height information was available for the wind measurements of CPCB, height adjustments were not applied in the figure. The simulated 10 m wind data were compared with the observed data, which seemed to be < 10 m a.g.l. or within the urban canopy, as previously discussed, because the simulation yielded substantially larger values than the observations. Nevertheless, good agreement was achieved between the simulation and observations in terms of temporal variations for wind speed, except on 8 and 9 November 2022, which was likely due to clouds and precipitation, especially over Delhi. As wind direction is affected by the local environment and mesoscale meteorology, the derivation of station-mean wind direction was not overly meaningful. Nevertheless, the simulated wind direction matched well with observations over Haryana, where northwesterly winds prevailed over the Plume 1 (2–4 November 2022) and Plume 2 (8–12 November 2022) periods. In Delhi, the simulated winds were northwesterly during the two plume periods, whereas those in observations were southwesterly. This was likely due to many of the selected stations being located within the urban canopy, causing the simulated wind direction to deviate from observed values. This difference in wind direction in Delhi may have caused errors in emission optimization, but the effect could not be quantitatively derived. The observed wind directions in Punjab during the two plume periods were almost northwesterly. Similarly, the simulated values were almost northwesterly, although they sometimes deviated from observations, especially in the early stage of Plume 2 (8–9 November 2022).

Figure 4 Same as Fig. 3 but for selected CPCB station-mean (a, c, e) wind speed (m s⁻¹), (b, d, f) wind direction (degree), (g, i, k) temperature (°C), and (h, j, l) relative humidity (%, left axis) over Punjab, Haryana, and Delhi, respectively. The simulated (blue lines) and observed (sky blue bars) surface precipitation amounts (mm h⁻¹) are also depicted in the panels on the right axis (h, j, l).

Figures S4–S7 present the time series of spatial distributions of GFAS CO emissions used in the simulation, MODIS AOD and CF, and TROPOMI UVAI for the entire CRB period (27 October to 15 November 2022). The simulated surface concentrations of anthropogenic PM_2.5 and OFEv1.0 CRB PM_2.5 and wind vectors are also presented in the figures. Notably, the data in the figures were obtained around noon, but the timing differed slightly: approximately 10:30 IST for MODIS, 12:00–15:00 IST for TROPOMI, and 11:30 IST for the simulation. Additionally, MODIS and TROPOMI represent column values, whereas the simulation data are in the surface air (bottom layer of the model grids).

During the Plume 1 period (1–5 November; Figs. S5), remarkable signs of CRB emissions were observed over Punjab due to substantial underestimations of simulated PM_2.5 in the CUPI-G source region data, as shown in Fig. 3b. There were higher GFAS emissions on 2 and 4 November 2022, corresponding to large values of MODIS AOD (> 2). MODIS retrieval was unsuccessful on 3 and 5 November 2022, likely due to the presence of clouds; thus GFAS emissions were small, whereas TROPOMI UVAI indicated a substantially large aerosol burden over the PHD region (> 3). Southeasterly winds prevailed on 1 and 5 November 2022 (Fig. S5), preventing transportation of surface PM_2.5 from CRB to Haryana and Delhi (Fig. 3). However, from 2 to 4 November 2022, northwesterly winds carried air pollutants from Punjab to Delhi, which caused high surface PM_2.5 concentrations of > 500 µg m⁻³. In fact, the Plume 1 period could be divided into two events, including Plume 1A on 2–3 November 2022 and Plume 1B on 3–4 November 2022. The Plume 1A event (station mean-values = approximately 600–700 µg m⁻³) was larger than the Plume 1B event (approximately 300–400 µg m⁻³) (Fig. 3d–f), likely because the former plume directly transported pollutants from Punjab to Delhi, while the latter plume was a blowback of pollutants previously carried downwind (southeast direction) of Delhi by the former plume (Figs. 4d, f, and S5).

Similarly, the Plume 2 period (8–12 November) could be divided into two events, including Plume 2A on 8–9 November 2022 and Plume 2B on 10–12 November 2022. During the Plume 2 period (Figs. S6 and S7), remarkable signs of CRB emissions on 8–9 November 2022 were observed in Punjab, as the CUPI-G (farmland) mean PM_2.5 concentrations were > 700 µg m⁻³ (Fig. 3b), and even the CPCB (urban) mean values became much higher (300 µg m⁻³) (Fig. 3a) than the “GFAS” simulation (100 µg m⁻³) (Fig. 3a, b). However, MODIS AOD data were not available during this period, and even TROPOMI UVAI did not detect a high concentration event in Punjab on 8 November 2022, likely due to thick clouds associated with rainfall over the PHD region (Fig. 4h, j, l). CRB cannot be conducted in the rain; however, CRB might have occurred because rain was not observed in Punjab during this period (Fig. 4h). The surface PM_2.5 concentrations in Punjab were high on 8–9 November 2022, whereas those in Delhi NCR were high on 10–12 November 2022 (Fig. 3b, f). However, unlike during the Plume 1 period, there were no remarkable enhancements of MODIS AOD and TROPOMI UVAI over Punjab on 8–9 November 2022 (Fig. S6) and Delhi NCR on 10–12 November 2022 (Figs. S6 and S7) during the Plume 2 period. In the latter period, the transport of air pollutants from Punjab to Delhi was observed in the simulation (Fig. S7). By contrast, CRB plumes were not transported to the Haryana region in the former period (Fig. S6), whereas an enhancement was observed in surface PM_2.5 at CUPI-G intermediate stations (Fig. 3d). As shown in Fig. 4a–f, discrepancies between simulated and observed wind fields were enhanced in the former period (8–9 November 2022), especially for Delhi NCR (Fig. 4e, f). The meteorological simulation might have failed to reproduce the air flows associated with rainfall from thick convective clouds, hindering successful optimization of CRB emissions during this period.

3.2 Optimization of CRB emissions

Various combinations of data selection, application, and evaluation could be proposed for CRB emission optimization, including optimization using CUPI-G data and validation using CPCB data. Another approach could be emission optimization utilizing data from Punjab and validation using data from Delhi. We needed emission data of all chemical components for the simulation of PM_2.5; however, not all data were available from observations, as described in Sect. 2.3. We attempted various combinations of data selection and found that employing more data yielded better results; therefore, we decided to use PM_2.5 data from the CPCB and CUPI-G stations over all regions as an optimization starting point.

Table 2 summarizes the data selection combinations employed in optimization and the obtained x_m values from Eq. (1). “CPCB + CUPI-G” indicates optimization using PM_2.5 data from all CPCB and CUPI-G selected stations. Optimization was performed for each plume period, so that different x_m values were obtained for the two plume periods. Because the simulated wind fields deviated from observations in the Plume 2 period (8–9 November 2022), x_m for the GFAS tag during the Plume 1 period (6.87 or 6.99) was applied for the whole simulation period (15 October–15 November 2022; 00:00 UTC). For the areal tags, only time-averaged data are presented in Table 2; however, x_m values were also obtained for different temporal tags, which were applied for each period. Optimized x_m values for all temporal and areal tags are presented in Table S2.

Table 2 Combination of observational data selection utilized for emission optimization and obtained x_m values.

^a Optimization using all CPCB and CUPI-G stations.
^b Optimization using selected CUPI-G Punjab stations.
^c Optimization using only CUPI-G stations no. 8 (Thikriwala) and no. 9 (Beauscape Farm) in Punjab state.

The time series of observed and simulated all selected CPCB and CUPI-G station-average PM_2.5 concentrations are presented in Fig. 5. Statistical scores between the simulated and observed surface concentrations are listed in Table 3, such as simulation to observation median ratio (Sim. $/$ Obs.), root mean square error (RMSE), and correlation coefficient (R). The “GFAS” simulation substantially underestimated surface PM_2.5 concentrations over the entire PHD region (Fig. 5a–c; Table 3); however the simulation using CRB emissions optimized by “CPCB + CUPI-G” data (Fig. 5d–f) was substantially improved (Table 3), especially for the Plume 1 (2–4 November 2022) and Plume 2B (10–12 November 2022) periods, demonstrated by increased emissions in the EP and SP&NH regions. GFAS emissions are generally smaller than CRB fire emissions in India (Cusworth et al., 2018; Wiedinmyer et al., 2023). Because the optimized x_m values were much larger than unity, the GFAS emissions may have been underestimated in our case as well.

Figure 5 Time series of all selected CPCB and CUPI-G station-average (black) observed and simulated (red) total and (blue) No_CRB PM_2.5 using (a–c) GFAS emission, (d–f) optimized fire emission (CPCB + CUPI-G), (g–i) optimized fire emission (OFE v1.0; average of CPCB + CUPI-G and CUPI-G Punjab sites), and (j–l) optimized fire emission (average of CPCB + CUPI-G and CUPI-G station numbers 8 and 9 in Punjab) over (left to right) source (CPCB Punjab + CUPI-G Source), intermediate (CPCB Haryana + CUPI-G intermediate), and Delhi NCR (CPCB Delhi + CUPI-G Delhi NCR) regions from 1 to 15 November 2022 (00:00 UTC).

Table 3 Statistical scores comparing simulated and observed PM_2.5 over the PHD region using several CRB emission cases.

^a Simulation to observation median ratio.
^b Root mean square error (RMSE).
^c Correlation coefficient (R).
^d Optimized Fire Emission ver1.0 (OFEv1.0).

Even though CRB emissions optimized by “CPCB + CUPI-G” data substantially improved the simulation of PM_2.5, the simulation continued to greatly underestimate the observed peaks in Plume 2A (8–9 November) over the entire PHD region. Therefore, additional optimizations were conducted using CUPI-G Punjab data only. The “CPCB + CUPI-G” optimization was insensitive to emissions from the NP, WP, and CP regions (x_m < 10⁻⁵), whereas the additional optimization “CUPI-G Punjab” enhanced emissions from CP. Because the additional optimization continued to underestimate surface concentrations, especially on 8 and 9 November 2022, an additional optimization was conducted using only data from CUPI-G station nos. 8 and 9 (“CUPI-G P89”), where substantial enhancement of PM_2.5 was observed on these dates. Subsequently, we constructed two additional optimized CRB emission inventories in the same manner as the “CPCB + CUPI-G” case, which were merged with “CPCB + CUPI-G” emissions by taking the larger x_m values among CPCB + CUPI-G and the additional cases, regarded as “CPCB + CUPI-G plus CUPI-G Punjab” and “CPCB + CUPI-G plus CUPI-G P89.” Note that optimization of “CUPI-G P89” was only conducted for the Plume 2 period and the same x_m values with “CUPI-G Punjab” were used for the Plume 1 period.

The two additional simulation results obtained using “CPCB + CUPI-G plus CUPI-G Punjab” and “CPCB + CUPI-G plus CUPI-G P89” emissions are shown in the lower half of Fig. 5 (Fig. 5g–l). Simulated PM_2.5 in the source region was improved due to enhanced emissions over the CP region on 8–9 November 2022 (Fig. 5g, l). However, no increases in the simulation were observed over the intermediate region and Delhi NCR in the same period (Fig. 5h, i, k, l), likely because the simulated wind fields deviated from those observed during this period, as shown in Fig. 4. Consequently, the reasons for the enhanced PM_2.5 concentrations observed over the intermediate region and Delhi NCR on 8–9 November 2022 (Plume 2A) were not identified in the current study. As shown in Table 3, substantial improvements were achieved in the simulations using optimized emissions compared to the original simulation using GFAS. Only small differences were observed among the three simulations using optimized emissions, but “CPCB + CUPI-G plus CUPI-G Punjab” was slightly better than the other simulations over the source region for Sim./Obs, RMSE, and R, which is referred to as OFEv1.0. Comparison against AERONET AOD at the Lahore and Amity University Gurgaon stations, as shown in Fig. S8, supported the same conclusion: the optimized emissions were better than GFAS, with no great differences observed among the optimized simulations.

The simulated source contributions are illustrated in Fig. 6. The anthropogenic contribution was calculated using the brute-force method, derived from the difference between the control run and the 20 % emission reduction run multiplied by 5 to reduce the effect of nonlinearity in chemical reactions (Kajino et al., 2013): the brute-force method with a zeroed-out simulation (no anthropogenic NO_x and NMVOCs) affects oxidant concentrations substantially and thus chemical production rates of secondary PM_2.5, which deteriorates the calculation of source contribution estimations. The brute-force method is easy to implement, and its validity has been confirmed via comparisons with other sophisticated methods, such as the decoupled direct method (Hakami et al., 2004; Choi et al., 2014; Kelly et al., 2015). The brute-force method with a 20 % reduction run was employed to obtain the anthropogenic PM_2.5, whereas the zeroed-out brute-force method was applied to obtain PM_2.5 from CRB emissions. The miscellaneous contribution was the remainder after subtracting all contributions, namely, natural emissions, upper and lateral boundary effects, and numerical errors.

Figure 6 Same as Fig. 5 but for simulated contributions of different emission sources: (gray) anthropogenic emissions; CRB emissions from (red) North Punjab, (pink) West Punjab, (yellow) Central Punjab, (green) East Punjab, (sky blue) South Punjab and North Haryana; (blue) fire from other areas; and (white) miscellaneous, which contains natural emissions such as biogenic, mineral dust, and sea salt, influences from lateral and upper boundary conditions, and numerical errors.

During the first half of November 2022, as shown in Fig. 6, a small but certain contribution of CRB emissions was continuously observed in the “GFAS” simulation in the source region (10 %–20 %, Fig. 6a), whereas the simulated PM_2.5 was only affected during the two plume events in the intermediate region and Delhi NCR by 10 %–20 % (Fig. 6b, c). After emission optimization using “CPCB + CUPI-G”, approximately 50 % of surface PM_2.5 concentrations accounted for CRB emissions in the entire PHD region during the Plume 1 period (2–3 November 2022). In the intermediate region and Delhi NCR (Fig. 6e, f), contributions from SP&NH were substantially increased up to approximately 50 % during the Plume 1A event (2–3 November 2022), whereas contributions from EP were also increased up to approximately 60 % during the Plume 1B event (3–4 November 2022). In the source region, CRB emissions from CP and EP contributed almost equally to simulated PM_2.5. Even though the transport of pollutants by northwesterly winds was clearly observed during the Plume 1 period, the largest contributions to PM_2.5 in the source (Punjab) and receptor (Haryana and Delhi NCR) regions differed from one another. After additional inclusion of CUPI-G Punjab stations into the emission optimization (Fig. 6g–l), the contributions from CP were enhanced approximately from 20 % to 50 % during the Plume 2A event (8–9 November 2022) only in the source region, and the total CRB emission contributions were > 80 %. Although the observed concentrations in the receptor region became higher during the Plume 2A event, the simulation did not present any increases (Fig. 6h, i, k, l), likely due to the large discrepancy with simulated wind fields, as previously mentioned. During the Plume 2B event (10–12 November 2022), in the presence of prevailing northwesterly winds in the simulation, PM_2.5 in the receptor regions was substantially affected by CRB emissions from CP on 10–11 November 2022 (approximately 50 % of total PM_2.5) and WP, EP, and SP&NH on 11–12 November 2022 (Fig. 6h, i, k, l). The CRB emission contributions were > 80 % on 11–12 November 2022 in receptor regions (Fig. 6h, i, k, l). Similar to the Plume 1 period, the largest contributions to PM_2.5 in the source (Punjab) and receptor (Haryana and Delhi NCR) regions differed from one another during the Plume 2B event, with PM_2.5 greatly affected by emissions from NP and EP in the source region and emissions from WP, EP, and SP&NH greatly affecting the receptor region.

Figure 7 illustrates the time series of D02 total 12 h CRB CO emissions for GFAS and OFEv1.0. The OFEv1.0-to-GFAS CO emission ratio was approximately 7 (Table 2), and the two datasets are separately depicted on the left and right axes. Although the x_m values had almost no limit because optimization was repeatedly applied until the values converged, the time variation of OFEv1.0 was continuous without any sudden spikes or gaps. This indicated that the number of observational data available was sufficient for the number of tags. The optimization was successful for the recovery (or boosting) of emission fluxes, especially on 8 and 10 November 2022, when GFAS emissions may have been underestimated due to clouds or thick smoke/haze (see Fig. S6). Large differences were sometimes observed between optimized emissions in daytime (or a.m. in UTC) and nighttime (or p.m. in UTC), such as on 8 November 2022. Because each CRB event lasted only a few hours, the 12 h time resolution may not have been sufficient to reproduce CRB activity. Nevertheless, the study findings indicate the importance of diurnal variations in CRB emission estimations for air quality simulations. In addition, diurnal variation shape may not be the same on each day. For example, daytime (a.m.) emissions were larger on 2, 9, and 10 November 2022, while nighttime (p.m.) emissions were larger on 3, 4, and 8 November 2022. Judging from Fig. S5, GFAS emissions may have been underestimated due to the presence of clouds on 3 November 2022; however, OFEv1.0 did not demonstrate a substantial increase in emissions on this day (Fig. 7). This was likely because the observation stations were insensitive to CRB emissions that occurred on 3 November 2022, or not many CRB events occurred on that day.

Figure 7 Time series of D02 total 12 h CRB CO emissions (Gg CO 12 h⁻¹) of (sky blue, left axis) GFAS and (orange, right axis) OFEv1.0 from 1 to 15 November 2022 (in UTC). Time resolution of GFAS is 1 d, whereas OFEv1.0 is a 0.5 d. Daytime (day; 05:30–17:30 IST; or a.m., 00:00–12:00 UTC) and nighttime (ngt; 17:30–05:30 IST; or p.m., 12:00–24:00 UTC) are separately marked as open and closed circles, respectively.

3.3 Intercomparison of OFEv1.0 and other studies

Table 4 compares the emissions of OFEv1.0 and previous studies developed for post-monsoon CRB emissions in northern India on a total amount basis. The total GFASv1.2 emissions from 1 to 15 November 2022 (00:00 UTC) for CO, PM_2.5, OC, and BC were 111, 10.0, 5.04, and 0.517 Gg, whereas those of OFEv1.0 were 963, 86.1, 43.5, and 4.46 Gg, respectively. The OFEv1.0 emissions were 8.6 times larger than those of GFASv1.2 (same ratio for all chemical compounds, as previously explained). The daytime and nighttime emissions of OFEv1.0 did not differ considerably; for example, nighttime CO emissions (495 Gg) were 5.8 % larger than daytime CO emissions (468 Gg). The OFEv1.0-to-GFASv1.2 ratio for the non-tagged period (6.99 in Table 2) can be regarded as the general underestimation of GFASv1.2 on clear sky days. The difference between the 14 d total OFEv1.0-to-GFASv1.2 ratio (8.6 times) and the general underestimation ratio ($\mathrm{8.6}/\mathrm{6.99}$ = 23 %) can be regarded as the boosting of emissions that could not be detected due to clouds or thick smoke/haze. Liu et al. (2020) estimated that post-monsoon CRB emissions were boosted by 142 % considering cloud/haze gap fill, which was substantially larger than our estimates.

Table 4 Total emission amounts (Gg) of post-monsoon CRB emissions from northern India estimated in this study and previous studies.

^aOptimization was not performed from 15 October to 1 November 2022, but the general underestimation value of GFASv1.2 (6.99) was applied to boost emissions during the period for comparison with other emission datasets.
^b Starting from 00:00 UTC until 00:00 UTC.

Emission optimization was performed for 1–15 November 2022 (00:00 UTC); however, applying the same general underestimation ratio (6.99) to boost the emissions from 15 October to 1 November 2022 (00:00 UTC), yielded the monthly total emissions of OEFv1.0 (referred to as OFEv1.0^*) (Table 4). The OFEv1.0^* monthly total optimized emissions from 15 October to 15 November 2022 (00:00 UTC) were 1460, 130, 65.8, and 6.73 Gg for CO, PM_2.5, OC, and BC, respectively. Based on GFASv1.2, the monthly values accounted for > 85 % and approximately 90 % of the whole post-monsoon emission period (from 1 October to 1 December 2022; 00:00 UTC) over the D02 and CRB source regions (areal tagged regions), respectively. The monthly total emissions of OFEv1.0^* were consistent with those reported in previous studies, such as by Liu et al. (2020), at 65 ± 18 and 5.6 ± 1.6 Gg for OC and BC, respectively, and Beig et al. (2020), at 141.65 Gg for PM_2.5, but did not align with CO estimated by Liu et al. (2020), at 791 ± 225 Gg. Our emissions may have been overestimated because the same optimized factor obtained using PM_2.5 data was applied for all species. Certainly, optimization using observed and simulated CO data would be better in CO emission estimations. In fact, differences between observed and simulated CO during the Plume 1 and 2 events (Fig. S3a, c, e) are much smaller than those of PM_2.5 (Fig. 3a, c, e). Therefore, CO of GFASv1.2 may not be underestimated to the extent of PM_2.5.

As reported by van der Werf et al. (2017), Liu et al. (2020), and Wiedinmyer et al. (2023), additional use of VIIRS data with finer resolution (Δx = 375 m) to capture smaller scale fires, such as CRB, substantially increased emissions over South Asia compared to using only MODIS data (Δx = 1 km), such as GFASv1.2 (Kaiser et al., 2012). The values of the VIIRS- and MODIS-based inventories such as the GFED4.0s (van der Werf et al., 2017) and FINNv2.5 (Wiedinmyer et al., 2023) are also listed in Table 4. Both GFEDv4.0s and FINNv2.5 values are larger than GFASv1.2, but the magnitudes are very different. GFED4.0s values are smaller than Liu et al. (2020), Beig et al. (2020), and this study by more than half, whereas FINNv2.5 values are larger than those studies by more than twice. Substantial efforts have been made to improve the CRB emission estimates, but no consensus has been reached yet.

MODIS-based GFASv1.2 was used in emission optimization herein; however, fire emissions should also have been underestimated in the presence of clouds in VIIRS- and MODIS-based inventories such as GFED4.0s and FINNv2.5. Therefore, the emission optimization was necessary in any case, and the conclusions reached may have been the same but with different x_m values, even if GFED4.0s or FINNv2.5 data were used for the emission optimization instead of GFASv1.2. However, since the current emission optimization was based on period-mean values of GFASv1.2 emissions (Fig. 2), optimization may be improved by using GFED4.0s or FINNv2.5, which may better represent the spatial distributions of CRB than GFASv1.2 data.

3.4 Reconstructed air mass movement of anthropogenic and CRB PM_2.5 in November 2022

The time series of simulated surface PM_2.5 and satellite observation maps during the Plume 1 and Plume 2 periods are illustrated in Figs. 8 and 9, respectively. The major features were already described in detail and are summarized here. The general features of MODIS AOD and TROPOMI UVAI were similar, with higher resolution and more available data in the presence of clouds for TROPOMI UVAI. Although direct comparisons between the simulated and TROPOMI observed UVAI were not performed in this study, the plume shapes were similar during the Plume 1 period (Fig. 8). Cyclonic wind fields carried pollutants from the PHD region toward Pakistan on 1 November 2022; the air mass stagnated around the Punjab region on 2 November 2022; northwesterly winds prevailed over the PHD, transporting pollutants from Punjab to Haryana and to Delhi NCR on 3 November 2022 (Plume 1A); pollutants once carried further downwind of Delhi blew back on 4 November 2022 (Plume 1B); and cyclonic winds again carried pollutants toward Pakistan on 5 November 2022. These patterns were also observed by AERONET AOD at the Lahore station (Fig. S8a, c, e).

Figure 8 Spatial distribution of (top to bottom) simulated surface concentrations of anthropogenic PM_2.5, simulated surface concentrations of PM_2.5 from Optimized Fire Emission v1.0 (OFEv1.0), simulated surface concentrations of PM_2.5 from GFAS emission, TROPOMI UV Aerosol Index, and MODIS AOD from 1 to 5 November 2022, which includes the Plume 1 period (2–4 November 2022). The simulation time was 06:00 UTC (11:30 IST), the TROPOMI observation time was 06:30–09:30 UTC (12:00–15:00 IST), and the MODIS time was approximately 10:30 IST.

Figure 9 Same as in Fig. 8 but from 8 to 12 November 2022, which includes the Plume 2 period (8–12 November 2022).

During the Plume 2 period (Fig. 9), the general features of MODIS AOD and TROPOMI UVAI were consistent, although much more data were missing compared to the Plume 1 period (Fig. 8). Even though UVAI data were missing during this period, peaks were observed in surface PM_2.5 concentrations (Figs. 3 and 5), whereas no peaks were observed over Delhi NCR by the satellites. UVAI was >  3 and AOD was > 2 over Delhi during the Plume 1 period (Fig. 8), and these values were much lower during the Plume 2 period (Fig. 9). During the Plume 2A event, no simulated influences of CRB were observed over Haryana and Delhi, whereas anthropogenic emissions were high over these regions; however, almost no satellite data were available on 8 November 2022. Likely owing to convective activity, the wind patterns were too complex to be resolved by the low-resolution simulation (Δx = 6 km), which resulted in the simulated wind fields deviating from those observed. The air was almost stagnant from 8 to 10 November 2022. The shapes of plumes between the simulation and TROPOMI were similar on 9 November 2022; however, pollutants from TROPOMI seemed to be more widely distributed than in the simulation. Certainly, the difference in altitude between TROPOMI (column) and the simulation (surface) might have affected horizontal distributions. During the Plume 2B event, northwesterly winds prevailed over the PHD region, and pollutants were transported from Punjab to Haryana and to Delhi.

3.5 Contributions of CRB to surface PM_2.5 concentrations

Based on the simulation using OFEv1.0, the contributions of CRB to surface PM_2.5 concentrations and the source–receptor relationship among Punjab, Haryana, and Delhi NCR during the two plumes periods in November 2022 are summarized in Fig. 10. The plume duration periods are defined as follows: Plume 1A (2 November 2022 (12:00 UTC) to 3 November 2022 (11:00 UTC)), Plume 1B (3 November 2022 (12:00 UTC) to 4 November 2022 (23:00 UTC)), Plume 2A (8 November 2022 (00:00 UTC) to 9 November 2022 (23:00 UTC)), and Plume 2B (10 November 2022 (00:00 UTC) to 12 November 2022 (23:00 UTC)); thus, the durations were 24, 36, 48, and 72 h, respectively. During the Plume 1 period, the observed surface PM_2.5 concentrations were > 50 % higher than the half-month average values (1–15 November 2022; 00:00 UTC) over the receptor regions (Haryana and Delhi NCR). The mean concentration values and simulated contributions of CRB emissions from the source regions (Punjab and North Haryana) were higher during the Plume 1A event (274 µg m⁻³ and 50 % for Haryana and 419 µg m⁻³ and 50 % for Delhi NCR, respectively) than those during the Plume 1B event (219 µg m⁻³ and 40 % for Haryana and 324 µg m⁻³ and 50 % for Delhi NCR, respectively) because the Plume 1A event was the result of direct transport, while the Plume 1B event was the blowback. The mean observed concentration over Punjab and simulated contribution of CRB emission became highest during the Plume 2A event (312 µg m⁻³ and 47 %, respectively) due to the stagnated horizontal wind field predicted by the simulation. Although small increases were observed in surface PM_2.5 concentrations over receptor regions (Figs. 3c–f, 5b–c, e–f, h–i, and k–l), the optimized emissions could not reproduce the observed high concentration levels, which was likely due to low performance of the meteorological simulation of air flows associated with rainfall and convection. Thus, the source–receptor relationship among the three regions could not be identified for the Plume 2A event in this study. During the Plume 2B event, PM_2.5 from CRB emissions was directly transported from the source region to the receptor region, and the contributions became highest in Haryana (56 %) and Delhi NCR (55 %), whereas the observed PM_2.5 concentrations were not very large compared to those in the Plume 1 period or the half-month averages. This result may have been due to the larger mean wind speeds during the Plume 2B event (11–12 November 2022, in Fig. 9) than during the Plume 1A event (3 November 2022, in Fig. 8) or the amount of CRB emissions in Plume 2B being smaller than that in Plume 1A.

Figure 10 Schematic illustration of the contributions of CRB to surface PM_2.5 concentrations and the source–receptor relationship among Punjab, Haryana, and Delhi NCR during the (a) Plume 1 and (b) Plume 2 periods in November 2022. The period and duration of each plume (1A: 2 November 2022, 12:00 UTC to 3 November 2022, 11:00 UTC, 24 h; 1B: 3 November 2022, 12:00 UTC to 4 November 2022, 23:00 UTC, 36 h; 2A: 8 November 2022, 00:00 UTC to 9 November 2022, 23:00 UTC, 48 h; and 2C: 10 November 2022, 00:00 UTC to 12 November 2022, 23:00 UTC, 72 h) are noted in the figure. Blue circles, orange circles, and arrows indicate roughly simplified CRB emission flux, high concentration plumes, and wind direction, respectively, predicted by the simulation using OFEv1.0. The observed mean PM_2.5 concentrations are indicated in black or red, and gray values with brackets indicate the observed mean PM_2.5 for the first half of November 2022 (1–15 November 2022, 00:00 UTC). The values in red indicate when values exceed the periodical mean values by 50 %. The contributions of CRB from areas 1 to 5, namely, NP, WP, CP, EP, and EP&NH, are reported (in %) for each region and plume. The contributions of fire emissions from outside the PHD region are not included. The period-mean contributions of CRB from areas 1 to 5 were 34 %, 33 %, and 30 % over Punjab, Haryana, and Delhi NCR, respectively.

3.6 Uncertainties in our top-down emission optimization approach

Uncertainties in our emission optimization approach are discussed in this section, especially for the assumption of uncertainties in observation data and a priori simulation and uncertainties due to the predictability of simulated PBL.

As described in Sect. 3.2, we used the observation uncertainty σ_n = 20 µg m⁻³ (10 % of all station mean PM_2.5) and the upper and lower limits of the scaling factor x_m as (u_m, l_m) = (2.0, 0.5). The estimated total primary PM_2.5 emission from 1 to 15 November 2022 was 86.1 Gg. The sensitivity tests were conducted using σ_n ranging from 1 to 500 µg m⁻³ and (u_m, l_m) = (2.0, 0.5), (10, 0.1), and (100, 0.01), which resulted in the estimated PM_2.5 emission of 88.1 ± 3.01 Gg. On the other hand, when we used the constant observation uncertainty rate of 10 % together with (u_m, l_m) = (2.0, 0.5), the estimated total primary PM_2.5 emission was 31.6 Gg. The sensitivity tests using observation uncertainty rates ranging from 1 to 500 % and (u_m, l_m) = (2.0, 0.5), (10, 0.1), and (100, 0.01) yielded the estimated PM_2.5 emission as 30.6 ± 1.75 Gg. The simulation using the optimized emission derived using the constant uncertainty rates was not conducted, but judging from Fig. 5 and Table 3 the simulation would be substantially underpredicted because the optimized emissions using constant uncertainty rates (30.6 ± 1.75 Gg) were approximately 3 times smaller than those using constant uncertainty values (88.1 ± 3.01 Gg).

Emission flux estimates in the top-down approach are highly sensitive to the simulation of vertical profile of species, which in turn are highly sensitive to the simulation of boundary layer dynamics. In Fig. S9, vertical profiles of simulated CO and O₃ with and without OFEv1.0 were compared with the IAGOS aircraft observation over the IGI airport. On 4 November, when the CRB plume arrived at Delhi (Plume 1B), significant enhancements of pollutants due to CRB were observed up to approximately 2 km above the ground level (a.g.l.), with the highest concentrations at approximately 1 km a.g.l. for both CO and O₃. The simulated and observed profiles were consistent with each other on the day. On 7 November there were no enhancements of simulated CO and O₃ due to CRB. The simulation agreed well with the observations at 20:00 IST but was underestimated at 22:00 IST. There seemed to be some enhancements in the observed CO and O₃, which might be a sign of CRB contributions. However, due probably to insufficient performance of meteorological simulation during Plume 2A, the simulation did not reproduce higher observed CO and O₃ concentrations at around 1 km a.g.l. On 11 November, during Plume 2B, the simulated plume height (approximately 1.5 km a.g.l.) agreed well with observations, but the simulated vertical profiles of CO and O₃ below 1.5 km a.g.l. were not consistent with those observed. It is difficult to quantify the uncertainty due to predictability of PBL simulations only using surface air quality measurements and a single meteorological simulation. At the right moment, it indicated that our PBL simulation was not always wrong. In future, it will be necessary to use different meteorological simulations with different PBL schemes to quantify variations in optimized emission fluxes.

4 Conclusions and future remarks

The impact of post-monsoon CRB on surface PM_2.5 concentrations over the PHD region in northern India was investigated using a regional meteorology–chemistry model, NHM(WRF)-Chem, a high-density in situ surface observation network comprising CUPI-G stations, and the emission optimization technique. Emission optimization was applied for the Plume 1 (2–4 November 2022) and Plume 2 (8–12 November 2022) periods identified by Singh et al. (2023) using CUPI-G and meteorological analysis data.

In the source region (Punjab state), almost no enhancements were observed in surface PM_2.5 concentrations in CPCB stations mostly situated in big cities, whereas substantial increases associated with CRB were observed in CUPI-G data from stations in rural and farmland areas. Employing the CUPI-G data from Punjab enabled us to obtain optimized CRB emissions (Optimized Fire Emission v1.0; OFEv1.0) from 1 to 15 November 2022, which substantially improved the PM_2.5 simulation over the PHD region compared to using GFAS emission data. Diurnal variations of satellite-derived fire emissions have been previously unavailable. However, unlike forest fires, each CRB event lasts only a few hours; hence, information of diurnal variation may be crucial. Also, some farmers ignite fires after sunset (Liu et al., 2020). Thus, emission optimization was performed at 12 h resolution on a daily basis in the daytime (05:30–17:30 IST, or a.m., 00:00–12:00 UTC) and nighttime (17:30–05:30 IST, or p.m., 12:00–24:00 UTC).

The major findings of the study are summarized as follows:

• The total CO and PM_2.5 emissions of OFEv1.0 over the PHD region from 1 to 15 November 2022 were 963 and 86.1 Gg, respectively, which was 8.6 times larger than the original GFAS emissions. OFEv1.0 boosted CRB emissions that were substantially underestimated due to clouds or thick smoke/haze on 8 and 10 November 2022. The total emissions of OFEv1.0 were consistent with other relevant inventories for PM_2.5, CO, OC, and BC. Optimized daytime and nighttime emissions differed greatly, indicating that consideration of diurnal variations is crucial in emission estimations. Daytime emissions were larger than nighttime emissions on some days but not others, indicating that diurnal variation shape may differ for each day.

• Using OFEv1.0 and NHM(WRF)-Chem, the half-month (1–15 November 2022) mean contributions of CRB to the surface PM_2.5 concentrations over Punjab, Haryana, and Delhi NCR were 34 %, 33 %, and 30 %, respectively. As seen from Figs. 8 and 9, the emission increments by OFEv1.0 were observed over the Punjab region mainly, suggesting the origin of the CRB sources being concentrated in the state, while affecting the neighboring two states similarly. The role of CRB in Haryana was lower compared to the PM_2.5 enhancement in the region.

• The Plume 1 period was divided into two events, including Plumes 1A (2–3 November 2022) and 1B (3–4 November 2022). During the Plume 1A event, northwesterly winds prevailed over the PHD region and station-mean concentrations became larger in Haryana (274 µg m⁻³) and Delhi NCR (491 µg m⁻³), and the CRB contributions were both 50 %. The Plume 1B event was smaller than the Plume 1A event because it was a blowback of pollutants that once carried further downwind of Delhi. The station-mean concentrations of Haryana and Delhi NCR were 219 and 324 µg m⁻³, respectively, and the CRB contributions were 40 % and 43 %, respectively.

• Similarly, the Plume 2 period was divided into two events, including Plumes 2A (8–9 November 2022) and 2B (10–12 November 2022). The Plume 2B event was similar to the Plume 1A event (direct transport of pollutants from Punjab to Delhi due to northwesterly winds), with similar CRB contributions (56 % in Haryana and 55 % in Delhi); however, the PM_2.5 concentration was not as high as in Plume 1A. During the Plume 2A event, the air mass stagnated around the Punjab region, resulting in a high CRB contribution in Punjab (47 %). Although the observed PM_2.5 in Haryana was increased in Plume 2A, the simulation did not reproduce the increased emissions, likely due to low performance of the meteorological simulation or inability of the CUPI-G observation network to capture the emission source because the southern Punjab region of intense CRB was not covered by the 2022 observation sites.

Future issues are itemized as follows:

• The results of this study were obtained using a single model. Multi-model analysis is indispensable for better predictions and quantification of prediction uncertainties using different boundary conditions, meteorological models with different physical schemes, chemical transport models with different chemical schemes, different emission inventories, and different optimization techniques.

• Post-monsoon CRB emissions have a negligible impact on surface PM_2.5 on an annual basis and thus may not cause substantial long-term health effects (Guttikunda et al., 2023; Ghude et al., 2016). However, short-term exposure of vulnerable populations to high aerosol concentrations may lead to the need for emergency medical care and death (Krishna et al., 2021). Our model simulation using optimized emissions indicated a 50 % contribution to short-term exposure on a mass basis. While the 50 % value is specific to mass, corresponding information on a toxicity basis has not been determined. Toxicity per unit of aerosol mass may vary substantially depending on the chemical compound, size, mixing state (e.g., Das et al., 2020, 2021, 2023; Ching and Kajino, 2018), and emission source (e.g., Fushimi et al., 2021; Kajino et al., 2024). Further studies are warranted to determine the associations between chemical compounds and toxicity in order to perform accurate health impact studies.

• Similar issues, namely, “satellite-based observations are underestimated by clouds or thick smoke/haze” and “ground-based observations are mostly in cities and not near the fire emission source regions”, may also exist in other parts of the world. Our methodology, top-down estimation of emissions using distributed low-cost sensors, can be applied to other cases in the world.