Missing emissions from post-monsoon agricultural fires in northwestern India: regional limitations of MODIS burned area and active fire products

A rising source of outdoor emissions in northwestern India is crop residue burning, occurring after the monsoon (kharif) and winter (rabi) crop harvests. In particular, post-monsoon rice residue burning, which occurs annually from October to November and is linked to increasing mechanization, coincides with meteorological conditions that enhance short-term air quality degradation. Here we examine the Global Fire Emissions Database (GFED), whose bottom-up emissions are based on the 500-m burned area product, MCD64A1, derived from Moderate Resolution Imaging Spectroradiometer (MODIS) observations. Using a household survey from 2016, we find that MCD64A1 tends to underestimate burned area in many surveyed villages, leading to poor representation of small, scattered fires and consequent spatial biases in model results. To more accurately allocate such small fires and resolve sub-village heterogeneity, we use an experimental hybrid MODIS-Landsat method (ModL2T) to map burned area at 30-m spatial resolution, which results in 44 ± 21% higher burned area than MCD64A1 and up to 105 ± 52% increase in dry matter emissions over GFEDv4s. In our validation and assessments, we find that ModL2T performs better relative to MCD64A1 in terms of bias and omission error, but may introduce commission error due to conflation of burning with harvest and still underestimate burned area due to Landsat’s coarse temporal resolution (every 16 days). We conclude that while MODIS and Landsat provide more than two decades worth of observations, their spatio-temporal resolution is too coarse to overcome several region-specific challenges: small median landholding size (1–3 ha), quick harvest-to-sowing turnover period, prevalence of partial burning, and increasing haziness. To further constrain agricultural fire emissions in northwestern India and improve model estimates of associated public health impacts, integration of finer resolution imagery, as well as better understanding of the spatial patterns in burn rates, burn practices, and fuel loading, is requisite.


Introduction
India is embracing agricultural mechanization to increase crop productivity and decrease labor costs in order to fulfill food security demands for its rapidly growing population (Mehta et al 2014). Agriculture in India is currently only mechanized on 40%-45% of cropland, below that of the United States, Russia, western Europe, China, and Brazil (57%-95%) (Bai 2014, Mehta et al 2014. India's projected population surge from 1.3 billion in 2015 to 1.7 billion by 2050 demands sustainable increases in crop productivity, intensity, and yield, which in turn bolsters the rise of agricultural mechanization (United Nations 2015). Traditionally, farmers manually collect crop residue to feed livestock. However, as India mechanizes, farmers are using combine harvesters, which leave behind root-bound and scattered crop residues that are labor intensive to remove (Vadrevu et al 2011, Kumar et al 2015. Gupta (2012) estimates that rice residues in 90% of area harvested by combine harvesters are burned in Punjab, from which emissions can severely degrade regional air quality seasonally (Gupta 2012, Kumar et al 2015, Jethva et al 2018. However, the air quality impacts from agricultural fires remain highly uncertain due to differences in global fire emissions inventories that are coupled with atmospheric transport models . Here we assess the challenges of using satellite observations to map burned area and active fires in order understand where current emissions estimates are most underestimated and uncertain.
We focus on the post-monsoon burning season in northwestern India. Previous work using satellite fire detections and HYSPLIT atmospheric back trajectories suggests that pre-monsoon (April-May) wheat residue burning is of less concern to the Delhi National Capital Region's air quality than post-monsoon (October-November) rice residue burning due to different atmospheric transport patterns, higher ventilation from high boundary layer conditions, and less overall fire intensity (Jethva et al 2018. Smoke plumes from post-monsoon crop residue burning, primarily originating from agricultural states Punjab and Haryana, are transported across the densely-populated Indo-Gangetic Plain (IGP) (figure 1). In general, carbonaceous particles in smoke can be transported hundreds of kilometers in the atmosphere (Sharma et al 2010, Kaskaoutis et al 2014. Aside from air quality degradation and public health impacts, crop residue burning also inhibits the productivity of the next cropping season by reducing soil quality (Gupta et al 2004). However, the short timeframe to clear fields of rice residue and sow winter wheat is a key limiting factor, thus leading to increased combine harvester use and subsequent burning (Gupta 2012, Jain et al 2014. Thus, despite restrictions on agricultural burning, farmers continue to burn crop residue due to the lack of viable, well-incentivized, and cost-effective alternatives. In this study, our primary goal is to identify and assess the region-specific challenges of estimating fire activity using MODIS-derived burned area and active fire products, which are used as input in emissions inventories, for northwestern India. As an intermediate step, we develop a hybrid MODIS-Landsat method to experimentally downscale post-monsoon agricultural burned area to 30-m from the 500-m spatial resolution of the current MODIS burned area product, MCD64A1. Next, we validate burned area estimates using household survey data and make further assessments using 375-m Visible Infrared Imaging Radiometer Suite (VIIRS) active fire detections and MODIS aerosol optical depth (AOD). We evaluate the MODIS gridded active fire product using fine-resolution (<5 m) imagery. Finally, we discuss crop residue burning practices in northwestern India in the context of policy changes and increasing mechanization and land fragmentation.

Overview
The study area consists of two neighboring agricultural states in northwestern India, Haryana and Punjab (figures 1, S1 is available online at stacks.iop.org/ERC/1/011007/mmedia). Punjab and Haryana are situated at the heart of India's 'breadbasket,' where most farmers predominantly follow a rice-wheat rotation (Sekar and Pal 2012, Naresh et al 2013, Sidhu et al 2015. Table S1 summarizes the satellite-derived surface reflectance, fire, and land cover datasets, primarily from MODIS and Landsat, used in this study.

Burned area
Previous studies on high-resolution agricultural burned area estimation in northwestern India span 1-2 years of study (Yadav et al 2014a, 2014b, PRSC 2015. Here we use Google Earth Engine (Gorelick et al 2017) to expand the study time period to 14 years, from 2003-2016, and estimate the total extent of post-monsoon agricultural burned area at 30-m spatial resolution, improving on 'baseline' 500-m MODIS MCD64A1 burned area (Giglio et al 2018) with better spatial allocation of small fires. Our hybrid MODIS-Landsat method is a simplified version of the MCD64A1 global burn mapping algorithm and GFEDv4s approach of using active fires to boost burned area from small fires (Giglio et al 2009, Randerson et al 2012. MODIS 1-km active fire locations represent endmembers of larger clusters of small fires from which we can obtain the spectral signature and apply to Landsat at higher resolution. ModL2T is described in more detail in appendix S2. Figure S3 describes the workflow for the ModL2T algorithm, which can be summarized as follows: (1) pre-process individual scenes; (2) composite cloud-free scenes in pre-fire and post-fire collections; (3) define thresholds based on the quantile intersection of normalized burn ratio (NBR), a metric used extensively in burn scar mapping, in burned and unburned agricultural areas; (4) separately derive MODIS and Landsat burned area using NBR thresholds; and (5) merge Landsat and MODIS classifications and apply an agricultural mask.
We independently validate burned area by using a 2016 household survey on farm management practices across the IGP. The survey asks participants whether crop residue is burned before planting wheat. Because the survey responses inherently distinguish between burned and unburned fields, this validation addresses the conflation of burning with harvest. We use 1112 responses from farmers in 30 Punjab and 32 Haryana villages, spanning eight districts. Because the GPS coordinates associated with each response are not located in-field, we cannot match responses to individual fields. Thus, we group responses by village name and match mean GPS coordinates with an accuracy of <10 m to village shapefiles. On average, 18±5 households were surveyed per village. We normalize the % households that burned crop residue by approximate operated landholding area. We do not account for partial burns and assume a field is entirely burned if a farmer affirms crop residue burning. For comparison, we estimate the % BA ModL2T and BA MCD64A1 of total village cultivated area based on 30-m GlobeLand30 and 500-m MODIS MCD12Q1 land cover, respectively. Due to these normalized approximations spurred by data limitations, the two metrics of % burning per village are not directly comparable on a 1:1 basis. We further assess BA MCD64A1 and BA ModL2T with simple checks using: (1) higher resolution active fire locations from VIIRS (pixel-level), (2) previous burned area estimates (district-level) and (3) satellite AOD (region-level). These assessments and their caveats are described in detail in appendix S3.4. We next estimate the maximum relative increase in agricultural dry matter (DM) emissions from the Global Fire Emissions Database (GFED), version 4 s with small fires (Randerson et al 2012, van der Werf et al 2017, by using ModL2T and MCD64A1 C6 burned area; DM emissions can be converted to other chemical species (e.g. CO 2 , CO, OC, BC) using emissions factors . While DM emissions are often calculated using fuel loading and combustion completeness estimates in addition to burned area in bottom-up approaches, we exploit the highly linear GFEDv4s DM/BA slope for each 0.25°× 0.25°grid cell to directly scale BA to DM (appendix S4).

Active fires
We use Google Earth's sparse collection of fine-resolution (<5 m) historical imagery (DigitalGlobe and CNES/ Airbus) to validate the MODIS MOD/MYD14A1 gridded active fire products. Using all publicly available DigitalGlobe and CNES/Airbus imagery, we estimate omission error using >500 identified ignition hotspots, spanning >400 1-km pixels, for 34 different days over Oct-Nov, 2010 we can pinpoint these active fires by tracing smoke plumes back to individual fields. We also categorize each ignition as a complete or partial burn to assess variations in satellite detection of fires related to the method of burning. We define complete burns as burn scars that extend across entire fields in which both intact and loose residues are burned and partial burns as circular or ring-shaped burn scars usually located in the center of fields, where loose residues are stacked. We can then assign a date to the in-progress fires based on the scene acquisition date, adjusted to local time, and determine whether MODIS indeed detected these fires on the same day.

Landholdings and mechanization
We consider ancillary data on landholdings and combine harvester use to assess trends in land fragmentation and mechanization. The Indian Department of Agriculture, Cooperation, and Farmers Welfare conducts the agricultural census and provides two online quinquennial databases: agricultural Census and Input Survey. The Agricultural Census database, which is based on census and input sample survey data, includes detailed data on landholdings in India from 1995-96 to 2010-11 (http://agcensus.dacnet.nic.in/); the Input Survey database contains information on agricultural implements and machinery, including combine harvesters, from 1996-97 to 2011-12 (http://inputsurvey.dacnet.nic.in/). In addition, the 2016 household survey asks participants about rice harvesting methods. Response choices include: fully mechanical (e.g. combine harvester), partially mechanical (e.g. thresher), and manually. We exclude 140 responses from farmers who never harvested rice. Because one main source of uncertainty is the MODIS burned area and/or active fire estimates used as input in these inventories, here we assess the region-specific limitations of these satellite-derived fire products.
We independently validate burned area with household survey data from 2016. We compare post-monsoon village-level survey crop residue burning rates, normalized by total landholding area, with BA MCD64A1 and BA ModL2T expressed as a fraction of cropland area. The village-level fraction of surveyed households that burn crop residue is moderately correlated with fractional BA ModL2T (r=0.67, p<0.05) (figures 3(c), 4(a)).
ModL2T underestimates burn rates for villages with high fractional burn rates (0.9-1), which may be partly due to partial burning and uncertainties in agricultural area mapped by GlobeLand30 ( figure 4(b)). BA MCD64A1 achieves a weaker correlation of r=0.6 (p<0.05) with higher normalized mean bias and severe underestimates in burned area for many villages, skewing its distribution toward low fractional burn rates (figure S8).
We assess omission and maximum commission errors based on the co-location of VIIRS active fire detections with BA MCD64A1 and BA ModL2T , from 2012-2016. With a higher spatial resolution (375 m) than MODIS/Terra and Aqua (1 km), VIIRS more consistently detects smaller and cooler fires (figure S9). We find that BA ModL2T yields a lower omission error (1%-8%) than BA MCD64A1 (33%-46%) (table S4). The maximum commission error is much higher for BA ModL2T (43%-55%) than BA MCD64A1 (10%-19%), but may reflect undetected active fires outside VIIRS overpasses or those obscured by thick haze or clouds. In particular, BA MCD64A1 is often unable to detect active fire hotspots in regions with prevalent partial burning, such as in  figure S9(c)). In addition, VIIRS detected 41% of grid cells burned consecutively from 2012-2016, while MODIS detected only 14% of grid cells by this criterion.
Finally, we assess detrended interannual variations in BA ModL2T and mean post-monsoon AOD from the MODIS merged Dark Target and Deep Blue product. Similar to daily FRP-AOD relationship quantified in Liu et al (2018), we find that regional BA ModL2T is weakly positively correlated with mean regional AOD (r=0.46, p=0.1), but not statistically significant (figure S12(a)). Comparatively, BA MCD64A1 is unexpectedly anticorrelated with AOD (r=−0.54, p<0.05) ( figure S12(b)). Similar correlations are observed using groundbased AOD measurements from the Aerosol Robotic Network (AERONET; figures S12(a), (d)).
3.3. Validation of active fires with fine-resolution imagery: two burning practices Two main crop residue burning practices are observed in Punjab and Haryana: complete and partial burns (Gupta 2012, Kumar et al 2015. Although farmers employ a mixture of the two practices, mapped active ignitions from available fine-resolution imagery show that complete burns are widespread in Punjab and northern Haryana, while partial burns are more pervasive in central and southeast Haryana ( figure 3(b)). Complete burns induce dark scarring over entire fields such that adjoining fields burned in this way within days of each other are starkly contrasted against the surrounding unburned landscape (figures 5(a), (b)). Partial burns leave small, circular or ring-shaped scarring in the center of fields; only ∼1/9 of the field area is scarred (figures 5(c), (d)). We find that the MODIS active fire product poorly matches in-progress fires identified from available fine-resolution imagery. Same-day omission error is 95%, with all co-locations from complete burns (table 2). Same-season omission error decreases to 75%, suggesting active fires within the same 1-km pixel were detected on other days.

Trends in landholding size, combine harvesters, and agricultural burning
The median landholding size in Haryana (1-2 ha) is smaller than that in Punjab (2-3 ha) (figure S16). After some consolidation of small landholdings from 1995-96 to 2000-01, landholdings became increasingly fragmented from 2000-01 to 2010-11. From 1996-97 to 2011-12, the number of combine harvesters increased over 20-fold from 14 664 to 297 132 in Haryana and almost 3-fold from 93 191 to 256 162 in Punjab. Based on the 2016 household survey, 72% of farmers using a combine harvester to harvest rice subsequently burned the crop residue in preparation for sowing wheat in Punjab and Haryana, compared to manual harvesting (8%) (table 3). Mechanization is less strongly linked to burning in Haryana, where only 32% of farmers using combine harvesters burned rice residue, compared to 87% in Punjab. Overall, of those who burned rice residue, 98% had used fully or partially mechanical methods of harvesting.
Overall, BA ModL2T increased by 966±179 km 2 yr −1 (p<0.05), or 82% in total, from 2003-2016. While increased Landsat scene availability (figure S2) may account for the some of the upward trend in BA ModL2T , the upward trend in BA MCD64A1 , which has no dependency on Landsat, is higher at 974±85 km 2 yr −1 (p<0.05), or 142% in total. Over the same 14-year time period, mean Oct-Nov satellite AOD increased by 39% overall, or 0.017±0.003 yr −1 (p<0.05); increased aerosol loading during the post-monsoon is also apparent from ground-based column AOD measurements from the Aerosol Robotic Network (AERONET) site at Lahore (in the neighboring Pakistan province of Punjab) (figure S13). In addition, we find an average step increase of 54%-65% in BA ModL2T and BA MCD64A1 from the 2003-07 to 2008-16 time period. Building on Thumaty et al (2015), we find that the timing of peak post-monsoon burning, indicated by 5-day block mean active fire counts, has shifted from mid-October to early-November during the 14-year period (figure S6). In 2009, the Punjab and Haryana governments officially implemented the 'Preservation of Sub-soil Water Act, 2009' (Ordinance in 2008 to counteract groundwater depletion by delaying rice transplanting to after June 10 and 15, respectively (Singh 2009, Tripathi et al 2016. In effect, this policy forces the rice harvest season to extend to mid-November (Singh 2009, PRSC 2015, which may explain the abrupt increase in burned area around 2008. Further analysis is needed to robustly quantify these temporal shifts in post-monsoon burning and substantiate the link to the groundwater policy.

MCD64A1 and ModL2T burned area: validation, assessments, and uncertainties
In northwestern India, increasing rates of post-monsoon agricultural burning enhance downwind air quality degradation and are linked to more widespread use of mechanized harvesting methods. Emissions estimates for agricultural fires in northwestern India are poorly constrained, on average ranging from 12-119 Tg OC+BC over the post-monsoon burning period among five inventories. In this study, we target the MODIS-derived burned area estimates used as input in GFEDv4s. MCD64A1, a primary input in GFEDv4s, is known to perform poorly in various agricultural regions (Giglio 2015, Hall et al 2016, Fornacca et al 2017, Lasko et al 2017, Zhu et al 2017. We combine MODIS and Landsat imagery to experimentally improve the spatial allocation of postmonsoon agricultural burned area in northwestern India for 14 years from 2003-2016. Use of Landsat imagery has been primarily limited by: (1) its low temporal resolution (16 days) and (2) storage and computing power. To minimize these limitations, we implement a hybrid MODIS-Landsat approach in GEE to rapidly process large collections of MODIS and Landsat imagery and expand the spatio-temporal range of study. Our simplified methodology is subject to several limitations, such as inconsistent Landsat availability, averaging of NBR across Landsat platforms, region-averaged NBR thresholds, and assumption that timing of the crop cycle is relatively homogeneous. Nevertheless, we find that incorporating Landsat imagery can improve the spatial allocation of small fires in northwestern India, which is important for modeling studies in which small fire emissions in close proximity to population centers can significantly impact local air quality estimates.
In comparison to MCD64A1, the ModL2T algorithm estimates on average 44±21% higher burned area in Haryana and Punjab during post-monsoon, from 2003-2016. ModL2T allocates burned area for partial burns in Haryana that are largely unaccounted for in MCD64A1. Validation of burned area with household survey data in 2016 suggests that the ModL2T algorithm can estimate burned area with increased accuracy (r=0.67, NMB=−25.7%), compared to MCD64A1 (r=0.6, NMB=−28.6%). In additional assessments, we find that BA ModL2T improves on BA MCD64A1 in terms of omission error, comparison with previous estimates of burned area, and relationship with satellite AOD, but may introduce commission errors (appendix S3.4).  Table 2. Validation of MODIS MOD/MYD14A1 active fire products using geolocations of ignitions identified from fine-resolution imagery. The spatial distribution of the ignitions is shown in figure 3(b).

Ignitions
MxD14A1, co-location with ignition pixels  4.2. Limitations of burned area and active fire algorithms in northwestern India We identify several key limitations that make the spatio-temporal resolution of MODIS, Landsat, VIIRS, and INSAT-3D insufficient for detecting active fires and accurately mapping cropland burned area in northwestern India: (1) prevalence of partial burning, (2) small landholding sizes and increasing fragmentation, (3) short duration of fires, (4) possible commission error from conflation of burning with harvest, (5) quick harvest-tosowing period lead to missed fires, and (6) increasing haziness that limits satellite observing area. Based on the two dominant types of burning practices in Punjab and Haryana, partial burning, which is prevalent in central and southern Haryana, may be more difficult to detect due to sub-landholding size fires and likely lower thermal energy. This difficulty is compounded by small median landholding sizes in Haryana (1-2 ha) and Punjab (2-3 ha). The inability of MODIS to readily detect partial burns and its tendency to homogenize over clusters of fields means that GFED4s grid cells with mostly partial burning are likely to contain a small sample of small fires, or none. This implies that the potential of the GFEDv4s small fires boost is limited in these areas in particular, and that the spatial allocation of these small fires is also not well-represented.
The GFEDv4s small fire boost relies on active fire hotspots and dNBR-based ratios from 16-day MODIS surface reflectance composites (Randerson et al 2012). This methodology assumes a linear correlation of burn severity with burned area. However, unlike wildfires, whose burn severity and burned area extent can vary greatly, cropland fires are generally controlled in burn rate, time, and area, thus limiting the upper bound of burn severity and burned area extent per fire. For cropland fires, dNBR has been used more as a threshold for burned area classification rather than a proxy for burn severity (e.g. McCarty et al 2008, 2009, Oliva and Schroeder 2015, Zhu et al 2017. However, the decline in NBR at the end of the growing season is influenced by both harvest and burning (Hall et al 2016). Clearly attributing decreases in NBR to burning remains challenging due to noise in the daily NBR timeseries. Further, the limited harvest-to-sowing turnaround period during post-monsoon means that burning may immediately follow harvest (Kumar et al 2015); we find that burn scars can disappear as soon as within several days. The low temporal availability of Landsat further increases its susceptibility to low pixel availability from haze and clouds, possibly leading to large mismatches in the satellite acquisition date between neighboring scenes. We conclude that both Landsat and MODIS surface reflectance products (8-day and 16-day) are fundamentally too temporally coarse to accurately classify burned area.
Further, post-monsoon BA ModL2T and BA MCD64A1 are not correlated with total active fire counts (r=0.05-0.1), when detrended, which may reflect differences in the active fire and burned area algorithms. Unlike burned area, active fires are derived from thermal anomalies and thus not susceptible to conflation of burning with harvest. However, slant satellite viewing geometry may dilute the signal of small, short-lasting fires, whose detection is already hindered by coarse sensor spatio-temporal resolution. In India, agricultural fires typically last no more than half an hour (Thumaty et al 2015). We find an omission error of >90% by the MODIS active fires product. VIIRS, at a higher 375-m spatial resolution, detected active fires in ∼20% more 0.02º grid cells than MODIS. Even so, VIIRS would not be able detect fires obscured by haze and clouds and those outside of its overpass time. Li et al (2018) showed that even slight differences in VIIRS and MODIS/Aqua overpasses of ∼15 min can lead to large discrepancies in active fire detections over Punjab. In addition, cloud cover and increasing haziness, indicated by AOD, can limit retrieved scenes that are usable and block active fires from satellite detection . The short return time (30 min) of INSAT-3D, a geostationary satellite, makes it ideal for capturing short-lasting agricultural fires, but its coarse 4-km spatial resolution makes detection of such fires unviable (appendix S3.4, figure S14). Additionally, analysis on the impact of smoke aerosols on public health and climate, such as on radiative forcing and aerosol-cloud interactions, will also need to be refined with fine-resolution sensors.

Future directions for improving agricultural fire emissions
The recent proliferation of finer resolution satellites, such as S-NPP (375 m and 750 m, daily, post-2012), Sentinel-2 (10-20 m, every 5 days, post-2015) and Planet (<5 m, daily, post-2016), offers added potential for active fire and burn scar detection (Drusch et al 2012, Strauss 2017. For more recent years of study, these imagery can help to better constrain the spatial and temporal variability of agricultural fire emissions in northwestern India. In particular, partial burns are difficult to detect at Landsat resolution, but potentially discernable with Sentinel or Planet imagery. Additionally, the present inability of moderate-resolution sensors to detect partial burns also raises the question of how end-users of fire emissions inventories should account for these missing emissions. More detailed on-the-ground knowledge of the amount of crop residues generated, as well as burn rates and practices, is needed to inform inventories retroactively. Differences in burn scar area from complete and partial burns also imply that separate fuel loading estimates are needed. Additional uncertainty in post-monsoon smoke OC+BC emissions, which differ by an order of magnitude among five widely-used inventories, signals a need to evaluate not only the satellite fire products used as input, but also differences in statistical boosts applied, emissions factors, and fuel consumption estimates. Due to region-specific limitations of MODIS burned area and active fire products, it is likely that atmospheric models using current global fire emissions inventories considerably underestimate smoke exposure and public health impacts from agricultural fires. Future collaborations to collect extensive ground truth data and incorporate accurate estimates of emissions can provide more robust input for policy decisions.