Long term influence of groundwater preservation policy on stubble burning and air pollution over North-West India

Stubble burning (SB) has been a major source of seasonal aerosol loading and pollution over northern India. The aftereffects of groundwater preservation act i.e., post 2010 era (2011–2020) has seen delay in crop harvesting thereby shifting the peak SB to May (Wheat SB) and to November (Paddy SB) by 8–10 and 10–12 days compared to pre-2010. Groundwater storage depletion rate of 29.2 mm yr−1 was observed over the region. Post 2010 era shows an increase of 1.4% in wheat SB and 21% in Paddy SB fires over Punjab and Haryana with 70% of PM2.5 air mass clusters (high probability > 0.8) advecting to the downwind regions leading to 23–26% increase in PM2.5 and 4–6% in aerosol loading over National Capital Region (NCR). Although the objective of water conservation policy was supposed to preserve the groundwater by delaying the paddy transplantation and sowing, on the contrary the implementation of this policy has seen groundwater storage after 2013 depleting at a rate of 29.2 mmyr−1 over these regions. Post policy implementation has led to shift and shrinking of harvest window with increased occurrences in SB fires which also increase associated particulate matter pollution over North India.

www.nature.com/scientificreports/ external advection and about 50-75% PM 2.5 concentration is contributed during the SB period from the upwind region 15 . It is estimated that 66,200 PM 2.5 attributable deaths (6.1%) are caused by stubble burning in India in 2015 16 . Introduction of post 2009 GW act, has observed a 10-day shift in timing of post-monsoon SB over NW India 17 and its influence on surface air pollution (PM 2.5 ) over Delhi is estimated ranging between 7-78% of the maximum observed PM 2.5 enhancements in Delhi to fires during post-monsoon season 18 . We investigated the impact of GW act on the changes & pattern in stubble burning, its impact on air pollution in the region and status of groundwater storage thereafter.

Results
Long term variation in stubble fires. Punjab (Fig. 1a,b). From satellite based daily fire counts, it is observed that average PoM SB (~ 16,000 fires yr −1 ) are 4-5 times higher than the PrM SB (~ 3600 fires yr −1 ). The non-parametric Mann-Kendall test shows a significant increasing trend with a rate of 0.21% yr −1 and 1.51% yr −1 (P < 0.05) for PrM SB and PoM SB respectively during 2002-2020. Over the past 19 years, paddy stubble fires have increased by 36% (245 fires yr −1 ) as compared to wheat stubble fires of 4% (15 fires yr −1 ) (Fig. 1a,b). Post-2010 era, witnessed PoM SB fires increased by 21% (compared to pre-2010 era) while PrM SB fires increased by 1.4% (compared to pre-2010 era). Monthly fire counts of post-2010 era, reveal that the April SB fires have decreased (0.4% yr −1 ) and fires in May have increased (0.65% yr −1 ) while decrease of stubble fires in October (2.1% yr −1 ) and a significant jump in November (3.2% yr −1 ) (Fig. 1a,b). This clearly establishes the fact that the distribution of majority stubble fires after year 2010 has shifted from April to May and substantial shift from October to November. Fire density helps in understanding the spatial distribution, variation and intensity of SB fires (Fig. 2). The fire density has been classified as High Fire Prone region (HFP) (fire density > 300), Moderate Fire Prone region (MFP) (150-300 fire density) while Low Fire Prone region (LFP) (50-150 fire density). It is observed that during PrM SB , maximum fire density is in LFP zone which has increased by 21% while during PoM SB the hotspot region under HFP has increased by ~ 51% while MFP region by 32% from pre-2010 era ( Fig. 2a-d). This is clear indication that in-spite of the enforcement of the GW act (post-2010 era) there is relatively marginal increase in fire occurrences in PrM SB while substantial jump in PoM SB burning, implying that the stubble burning practice is continuing. This increase and peak shift of burning pattern (1st week of November) in PoM SB season coincide with the festive season in India with high anthropogenic carbonaceous dust, smoke and aerosol particulate matter (PM 2.5 , PM 10 ) adversely affecting air quality and are advecting to regions of National Capital Region (NCR) and parts of Indo-Gangetic Plain.
Shift in pattern of stubble fires. During pre-2010 era, PrM SB practice over Punjab and Haryana started from 2nd week of April-2nd week of May (~ 31 days) (Fig. 3a, (Fig. 3a,b). Now, the length of PrM SB episode is from the 2nd week of April and lasts till 3rd week of May (~ 39 days) and PoM SB window is from 2nd week of October to about 4th week of November (~ 48 days). Peak stubble burning during the post-2010 era was recorded during 1-4 May & 2-6 May and 31st October-4th November & 2-8 November in Haryana and Punjab respectively. Similar PoM SB burning peak during post 2010 era was also observed elsewhere 19 . This clearly establishes that stubble burning events have   In order to understand the relationship of the SB fires to AOD and PM 2.5 , its daily average variation during 2010-2020 over NCR region was analyzed (Fig. 5). The background emissions particularly PM 2.5 in last six years have reduced due to implementation of several policies like Graded Response Action Plan (GRAP), Comprehensive Action Plan (CAP), and National Clean Air Programme (NCAP) for prevention, control and mitigation of air pollution in Delhi and NCR 20 . This also resulted in a decrease in the number of days that exceeds the Indian National Assessment Quality Standards in and around Delhi and its nearby cities 21 . Further, 1457 out of 1542 coal and oil-driven industries have shifted to greener fuel, compressed natural gas (CNG) in New Delhi 22   . This substantiates the fact that with 32% increase of SB fires only in November intensifies aerosol PM emission over NCR by 23-26% affecting regional fog, visibility and health hazards.

Discussion
Satellite based MODIS observations clearly capture heavy smoke and dust particulates emitted from SB which are dispersed & advected to downwind regions of IGP (hundreds of km) during PoM SB ( Supplementary Fig. 5b on right) due to synoptic meteorology ( Supplementary Fig. 1b). From Potential Source Contribution Function (PSCF) and Concentration Weighted Trajectory (CWT) analysis, it is also evident that moderate to high PM 2.5 concentration of 160-240 µgm −3 are advected upto central and parts of east IGP (Fig. 6) These north westerly winds transport the pollutants to NCR and get trapped owing to low boundary layer (< 500 m), low wind speed and direction, relative humidity and local meteorology which also prevents vertical mixing of pollutants 25,26 . Observation from Calipso shows the strong contribution of smoke aerosol that dominates the NCR and western IGP and present within 1 km of height ( Supplementary Fig. 6). Similar results were also observed wherein smoke plumes originating from stubble fires are mostly concentrated between surface and below 1 km altitude and depending on the steady meteorological conditions take around 3-5 days to get transported in well homogenized manner over distance of 800 km in south eastern regions (NCR, IGP) 6,27 . It is observed from PSCF and CWT  Fig. 5a on left). Clearly there is no external advection of aerosol particulate from hotspot region over IGP during NSB period ( Supplementary Fig. 5a,b on left). Noting that the area and production under rice and wheat in Punjab and Haryana did not change drastically after 2009 (production increase not more than 8%) 28   www.nature.com/scientificreports/ All the rivers flowing in Punjab and Haryana results from Himalayan snow melt, majorly contributing 58% and 50% to Punjab and Haryana river system respectively. The river system (Sutlej, Beas, Ravi, Chenab, Jhelum) in Punjab is spatially distributed with average flow rate ~ 6930 m 3 /s and the flow rate in Haryana region 30,31 (Yamuna, Ghaggar) is ~ 2950 m 3 /s. This infers that there is good recharge of groundwater in Punjab compared to Haryana and hence, Punjab has comparatively good groundwater condition than Haryana region. The results from GRACE/GLDAS data have shown that Punjab is witnessing low depletion in GWS, than in Haryana (Supplementary Fig. 2a-d). This ascertains that despite having better GW recharge and good amount of rainfall, Punjab is also facing similar chronic water condition like Haryana. Therefore, the GRACE and CHIRPS data analysis validates this fact and features the utter plight of groundwater storage condition in both Punjab and Haryana region.
We also investigated the proportional air mass clusters originating from hotspot region and influencing & arriving at receptor site (Delhi) 32 . During PrM SB around 40% of air mass clusters originating from hotspot region (north-west) are arriving at receptor site (New Delhi), 30% from oceanic region (south-west), 30% from western UP (north-east) and 10% from IGP (south-east) (Fig. 6a left). In PoM SB around 70% of air mass clusters from hotspot region are arriving at receptor site (New Delhi), 20% from IGP, 5% from oceanic (Fig. 6a right). This clearly establishes the analysis that most of the pollutants and aerosols reaching Delhi and NCR originate from the hotspot SB region. Further, potential sources of PM 2.5 observed over NCR were estimated using Potential Source Contribution Function (PSCF) and Concentration Weighted Trajectory (CWT) at a height of 500 m a.g.l as the aerosols and particulates are transported within the lower troposphere 6 . The high PSCF values (> 0.6) are considered as maximum probability sources to the receptor site and high CWT values (> 100 µgm −3 ) are the high strength sources (Fig. 6b). PSCF analysis confirms a strong probability (> 0.8) of PM 2.5 sources to be located in North-west region (Fig. 6a right) and high PM 2.5 air parcels are spread over the hotspot and nearby region (Fig. 6b). During PoM, SB about 70% of PM 2.5 trajectories are advected from hotspot region to the receptor site. We estimated the average contribution of Diwali event (considered ± 2 days from date of Diwali) to PM 2.5 concentration during the period of 2016-2020 for daily data over 6   www.nature.com/scientificreports/ concentration over Delhi 33 . This influx from hotspot region abruptly increases PM 2.5 concentration over NCR (250-300 μgm −3 ) especially Delhi (350-400μgm −3 ). The estimates also show that during 2010-2020, the PM 2.5 CWT reaching Delhi have increased at a rate of 1% (PrM SB ) and 3% (PoM SB ) per year. The analysis shows that NCR air quality is not only affected from stubble burning but also a combined effect of synoptic low and steady winds, low boundary layer, favorable wind direction and stable atmosphere. In order to improve the air quality over NCR and northern India, there needs to be long term and short-term steps to be taken. Government of India under National Action Plan for Climate Change (NAPCC) had enforced National Policy for Management of Crop residue (NPMCR) in 2014. Government is implementing ways to help farmers through use of subsidizing machinery equipment to cut stubble, promote diversified utilization of stubble for in-situ crop management, many industrial application and other uses. While the positive effect is being visible by way of slight improvement in air quality over NCR but will take time for substantial improvement. Under the short-term step, adopt rainwater harvesting & other water conservation technologies for optimal use of available water resources, renovation & use of village ponds, restore field distributaries & increase irrigation area under canal system, shifting from monoculture to diversified crop pattern and adoption to replace current rice variety with shorter duration rice varieties. This would make harvest 10-15 days (in September month) earlier than current practice period. During this time, the winds are moderate (3-4 m/s), substantial high temperature (5-6 C more than in November) and unfavorable wind direction (westerly to south-westerly) which does not transported and advect smoke to the NCR and adjoining areas in Northern India thereby reducing the number of high polluted episodes (PM 2.5 > 200 µgm −3 ) and not coinciding with the festive period 34 . This will also help in reduction in contribution of stubble burning share to PM 2.5 by 50% 34 .

Conclusion
Using long term satellite and ground measurements, we have investigated the trends in stubble fires, associated particulate matter pollution over NCR, downwind regions and linkage with groundwater storage. The present water preservation act was implemented to preserve the groundwater by delaying the paddy nursery and transplantation. However, this policy has led to delay in harvest & shrinking of harvest window with increased SB fire occurrences (15.8%) and shift in SB by 8-10 days (wheat) and 10-12 days (Paddy). The particulate matter pollution associated with SB fires coincides with cultural festivals in India thereby increased pollution over North India. There is no improvement in groundwater storage change after 2013 and superfluous usage for irrigation is leading to depletion of groundwater storage at a rate of 29.2 mmyr −1 over these regions and the GW act is marginalized in preventing the depletion of groundwater storage. Almost 21% increase in Paddy SB fires over a decade added additional 23-26% of annual PM 2.5 concentration over NCR. This high influx of aerosol particulate matter over NCR and downwind regions (particularly during PoM SB ) is deteriorating air quality in north India to as much of 4-5 times the permissible Indian standard. However, due to implementation of NPMCR policy & awareness programs among farmers by Government, the SB fires incidents and related pollution over NCR are on marginal decline (5%) over the last 4 years. It is envisaged that the approach and steps proposed in this study could be potential options for policy planners. This in turn will have positive impacts on air pollution with reduced mortality and morbidity rates.

Methods
Fire data collection and analysis. MODIS Thermal Anomalies/Fire locations-Collection 6 daily data of Aqua and Terra combined (MCD14DL) (distributed by LANCE FIRMS) having spatial resolution of 1 km x 1 km and a swath of 1200 km is used. Daily data from 2002-2020 for PrM SB and PoM SB was collected and classified into high confidence (> 80%), medium confidence (30-80%) and low confidence (< 30%) data. In order to capture the SB events where the fires are generally small (< 1 hectare of agricultural field) data > 30% confidence level data was used so as to capture smallest fire event. The fire density is calculated using MODIS fire data at 1 km resolution. The Climate modelling grid (CMG) fire data are gridded statistical summaries of fire pixel information intended for use in regional and global modelling and other large-scale studies and are generated at 0.25 deg. resolution for the time periods. Hence, for estimating regional fire density we used at standard 0.25 deg. square grid. For anlaysis, Pre-Monsoon (PrM SB ) (i.e., wheat stubble) is considered from 1April-31May every year while Post-monsoon (PoM SB ) (i.e., paddy stubble) is considered from 1 October-30 November every year. PrSB 20 is the sub-set of SB period where the burning is substantial (i.e., fires > 20 day −1 during PrM SB ) and PoSB 50 period (fires > 50 day −1 during PoM SB ). NSB is ± 30 days from the PrM SB or PoM SB period.
Aerosol optical depth. Daily Level-2, Collection-6 data (https:// ladsw eb. modaps. eosdis. nasa. gov/) MODIS/Aqua (MOD04_3k) and MODIS/Terra (MYD04_3k) AOD with a spatial resolution of 3 km and swath resolution approximately 2330 km wide was used to study the aerosol distribution in the region. The expected error in AOD is around ± 0.05 35 . For analysis, daily data has been collected, processed and analyzed for PrM SB , PoM SB and non-SB period during 2002-2020.
Air quality data. For analyzing the effect of SB on particulate concentration (PM 2.5 ) over NCR, the quality assured ground PM 2.5 continuous data available for more than 4 years was considered (https:// app. cpcbc cr. com/ ccr/#/ caaqm-dashb oard/ caaqm-landi ng/ caaqm-data-avail abili ty) is found over 19  Factors like temporal gaps and missing data can create discrepancies in the well datasets. Hence, it is imperative to select suitable wells for pre-processing so that requisite and accurate results can be achieved 38,39 . In the study, all 391 well data was used. GRACE derived groundwater storage and well groundwater level data were compared using Pearson coefficient. To make sure that groundwater storage depletion rate has not been affected by any other climatic factors, we investigated the TRMM (Tropical Rainfall Measuring Mission) rainfall data for these three consecutive periods (2003-2009, 2010-2013, 2014-2016). The average precipitation rate (cm) estimated over Haryana and Punjab region for the above time period is 53.5 cm, 61.2 cm and 54.1 cm respectively. The ground recharge rate for 2010-2013 may have increased but overall, over the years, there is seen a decrease in rate of replenishment due to overall irrigation practices. Climate Hazard group Infrared Precipitation with Stations (CHIPRS) is a quasi-global, reliable and contemporary dataset that is highly helpful in identifying early warning signatures like seasonal drought etc. It uses the information from satellite for representing the rain gauge data distributed sparsely in different locations, features 0.05° × 0.05° spatial resolution and measures on daily, monthly and pentad time scales. This data utilizes the TRMM (Tropical Rainfall Measuring Mission Multi-Satellite Precipitation) data of TMPA3B42v7 for measuring the rainfall estimates derived from CCD (Cold Cloud Duration) globally. It uses an approach of smart interpolation and cofunction with anomaly derived from climatology of high resolution. CHIRPS uniquely produces two products in 2 phase processes; in first phase it yields rainfall product with delay in two days and the gauge data from World Meteorological Organization's Global Telecommunication System (GTS) amalgated with the rainfall estimates (pentad) derived from CCD. While in the second phase the final product is obtained with delay in 3 weeks, and the pentad/monthly best available rainfall estimates derived from CCD is merged with the station data (derived monthly/pentad 40 . The Standardized Precipitation Index (SPI) uses CHIRPS monthly datasets for shedding light on drought conditions at various scales. It is an emerging method for estimating the drought index, and was formulated by the integrated effort of Tom Mckee, Nolan Doesken and John Kleist of Colorado Climate Center 41,42 . They distributed drought condition into several classes depending on its respective SPI value (Supplementary table 3). In this study we used SPI calculation on 6-month time scale (SPI-6) using CHIPRS (https:// www. chc. ucsb. edu/ data/ chirps), as this particular time scale is sensitive to the changes in rainfall (both seasonal and inter-seasonal). The 6-month SPI calculation includes current month precipitation and then compares with the precipitation of the exact six-month period of previous years. For instance; the SPI value at April-September month is compared with total precipitation of the same six months (April-September) period from previous years.
Air mass back trajectory and source apportionment. We modelled 5-day backward air mass trajectories at receptor site Delhi at 500 m a.g.l separately for PrM SB  www.nature.com/scientificreports/ sis archived data (ftp:// arlftp. arlhq. noaa. gov/ pub/ archi ves/ reana lysis). To investigate the sources, the air mass trajectories were computed for clusters. The PM 2.5 emitted source regions are identified using Concentration Weighted Trajectory (CWT) algorithm by incorporating actual ground measured PM 2.5 data (Source PM 2.5 data: CPCB). Further, to understand dispersion and advection of the SB emissions to influence regions, we used Potential Source Contribution Function (PSCF) which is a receptor model that incorporates meteorological information in its analysis scheme to produce a probability field that can be used to determine areas of the potential source contribution.
Statistical tests. Non-parametric Mann-Kendall test was conducted for the trend analysis with statistical significance established at p < 0.05. Pearson's coefficient was used to measure the correlation between satellite observed and ground in-situ data.