Detection of delay in post-monsoon agricultural burning across Punjab, India: potential drivers and consequences for air quality

Since the Green Revolution in the mid-1960s, a widespread transition to a rice–wheat rotation in the Indian state of Punjab has led to steady increases in crop yield and production. After harvest of the summer monsoon rice crop, the burning of excess crop residue in Punjab from October to November allows for rapid preparation of fields for sowing of the winter wheat crop. Here we use daily satellite remote sensing data to show that the timing of peak post-monsoon fire activity in Punjab and regional aerosol optical depth (AOD) has shifted later by approximately two weeks in Punjab from 2003 to 2016. This shift is consistent with delays of 11–15 d in the timing of maximum greenness of the monsoon crop and smaller delays of 4–6 d in the timing of minimum greenness during the monsoon-to-winter crop transition period. The resulting compression of the harvest-to-sowing period coincides with a 42% increase in total burning and 55% increase in regional AOD. Potential drivers of these trends include agricultural intensification and a recent groundwater policy that delays sowing of the monsoon crop. The delay and amplification of burning into the late post-monsoon season suggest greater air quality degradation and public health consequences across the densely populated Indo-Gangetic Plain.


Introduction
Rapid increases in mechanized harvesting in the Indo-Gangetic Plain (IGP) since the mid-1980s, together with steady increases in crop production, have led many farmers to burn the abundant residue left behind by this practice (Badarinath et al 2006). Such burning is a quick, cheap, and efficient method to ready the fields for the next crop. However, the smoke from post-monsoon crop residue burning, primarily during October to November, amplifies severe haze events in the region , such as that observed in early November 2016 . Of particular concern is the observed increase in aerosol loading associated with an upward trend in post-monsoon burned area and with a shift toward a later peak in post-monsoon fires in northwestern India (Thumaty et al 2015, Jethva et al 2018, 2019, Liu et al 2019, Balwinder-Singh et al 2019. Here we use daily satellite remote sensing data to quantify the temporal shift toward later burning in the state of Punjab, the 'breadbasket' of India. Such a shift would have implications for air quality degradation, since peak burning is more likely to coincide with meteorological conditions that are favorable in amplifying persistent haze, such as weak winds, low mixing layer heights, temperature inversions, and high relative humidity (Ojha et al 2020).
Agricultural intensification of rice and wheat in India has led to over two-fold and three-fold increases, respectively, in crop yield since the Green Revolution in the mid-1960s. In the western IGP, the predominant rice-wheat rotation is highly productive . Punjab, an agricultural state in northwestern India, contributes more than one-fifth of rice and one-third of wheat to the central grain pool in India, and thus generates large amounts of crop residue annually. Since the mid-to-late 1980s, farmers have increasingly used mechanized harvesting methods in preference to sickle-based manual harvesting in order to reduce labor costs and save time (Badarinath et al 2006. The use of combine harvesters, however, leaves behind an abundance of scattered and root-bound residue that is difficult to remove and thus often burned post-harvest to prepare for timely sowing of the next crop . The burning allows for quick disposal of crop residues and shortens the harvest-to-sowing transition from the kharif (monsoon crop) to rabi (winter crop) season. A quicker transition between crops also allows for earlier sowing of wheat during post-monsoon to avoid springtime heat (Lobell et al 2013).
However, the burning of post-monsoon rice residue can severely degrade air quality downwind of the agricultural fires over the IGP (Badarinath et al 2006, Jethva et al 2018, Sarkar et al 2018. In particular, smoke from rice residue burning in October and November may account for more than half the fine particulate matter (PM 2.5 ) concentrations in the Delhi National Capital Region , which already experiences intense urban pollution from local and other regional sources (Amann et al 2017). A temporal shift in fire activity to later in the year could exacerbate air quality degradation since late autumn-to-winter meteorology in the IGP favors smog formation due to weak winds, frequent temperature inversions, and a shallow boundary layer (Choudhury et al 2007, Saraf et al 2010. Observations from the moderate resolution imaging spectroradiometer (MODIS), aboard NASA's Terra and Aqua satellites, have been extensively used to investigate fire activity, crop yields, production, and phenology, and land use change detection. However, MODIS multi-day composites (8 d, 16 d) typically analyzed are insufficient to capture and resolve rapid changes in crop phenology (Zhao et al 2009).
Here we use daily active fire and surface reflectance data from MODIS to investigate trends in agricultural activity in Punjab. Specifically, we quantify the delays in post-monsoon agricultural fire activity and determine whether the seasonal cycle of monsoon to post-monsoon vegetation greenness reveals similar delays. We conclude with a discussion of the potential drivers of these interannual changes and an analysis of the consequences for regional air quality.

Study region
The IGP is home to over 700 million people (appendix S1.4, which is available online at https:// stacks.iop.org/ERL/16/014014/mmedia), many of whom rely on agricultural productivity of the densely cropped belt of northern India and parts of Pakistan, Nepal, and Bangladesh for livelihood and food security. Relative to other double-cropped states in northern India, such as Haryana, Uttar Pradesh, and Bihar, Punjab has the highest rice-wheat productivity  and is spatially more homogenous in terms of fire intensity (figure 1(a)), rice-wheat yields, and topography (Azzari et al 2017). Here we focus on Punjab during the post-monsoon rice residue burning season (defined here as September 20 to November 30), when fields are prepared for winter wheat sowing. To a lesser degree, we examine the pre-monsoon wheat residue burning season (April 1 to May 31), when fields are prepared for monsoon rice sowing (figure 1(b)).

Active fires and vegetation indices
For analysis of fire activity, we sum daily 1 km maximum fire radiative power (FRP), a proxy for fire intensity, derived from MODIS/Terra and Aqua (MOD14A1/MYD14A1, Collection 6). We also compare FRP with MODIS-derived fire counts and burned area and MODIS-based fire emissions from the Global Fire Emissions Database, version 4 with small fires (GFEDv4s) (table S1). For analysis of vegetation greenness, we use daily 500 m MOD-IS/Terra surface reflectance (MOD09GA, Collection 6) to derive two vegetation indices, the normalized difference vegetation index (NDVI) and normalized burn ratio (NBR): where ρ i is the surface reflectance of MODIS band i. The wavelength range of the bands is as follows: 620-670 nm for band 1 (red), 841-876 nm for band 2 (near infrared), and 2105-2155 nm for band 7 (shortwave infrared). These active fire and surface reflectance datasets are described in more detail in appendix S1.

Statistical analysis
We estimate linear trends with residuals bootstrapping. Unlike the linear regression t-test, which assumes that the residuals are normally distributed, bootstrapping preserves and resamples from the sample residuals distribution. To obtain a sample distribution, 1000 iterations are performed in which residuals are randomly sampled with replacement for each iteration and the dependent variable y re-fit using linear regression.

Characterizing the temporal progression of agricultural fires
We characterize the progression of the pre-monsoon and post-monsoon burning seasons, defined in section 2.1, of each year in order to assess interannual temporal trends. First, let X = {x t | x 1 , x 2 , x 3 . . . x n } denote the daily time series of a fire metric, such as FRP, during a given burning season lasting n days. We define the pre-monsoon (April 1-May 31) and post-monsoon (September 20-November 30) burning seasons with broad windows in order to capture all possible seasonal fire activity. To derive the peak burning date, k peak , we fit Gaussian density curves to X, thus smoothing potential noise in FRP due to inconsistencies in observing area caused by cloud and haze cover: where g(t) is the Gaussian function to be optimized, t is days of the burning season expressed as 1 to n total days, µ is the mean of t, σ is the standard deviation of t, and γ is an arbitrary scaling parameter. k peak is the day t that corresponds to the maximum value of g (t). We then use the optim function from the R stats package to minimize non-linear least squares of g(t) and X and to estimate the µ, σ, and γ parameters that yield the optimal Gaussian fit. As first guesses of the three parameters for the optim function, we set σ to 7, γ to 1, and µ to our estimate of the midpoint of the burning season, k midpoint,w . We estimate k midpoint,w using a simple weighted average approach. We weight t by x t and take the average to obtain k midpoint,w : Following Zhang et al (2014), we also estimate the start, midpoint, and end of the pre-monsoon and post-monsoon burning seasons for each year. We define a sequence of partial sums, Y = x t . We normalize y k by y n , or the sum of FRP during the entire burning season. We then approximate k β , or the first day when normalized Y has surpassed breakpoint β: As in Zhang et al (2014), we define arbitrary breakpoints, β = 0.1, 0.5, and 0.9, to represent the start (k β=0.1 or k start ), midpoint (k β=0.5 or k midpoint ), and end (k β=0.9 or k end ), respectively, of the burning season.

Tracking crop phenology with NDVI and NBR
NDVI is widely used to characterize the cycling in vegetation growth, land cover change, and crop productivity (Justice et al 1985, Yengoh et al 2015. NBR, while typically used in burned area and burn severity classification (Key and Benson 2006), is analogous to NDVI, which relies on the visible red reflectance instead of the shortwave infrared (SWIR) reflectance. A major advantage of NBR is that compared to visible wavelengths, SWIR wavelengths can better discriminate between vegetation and bare soil (Asner andLobell 2000, Chen et al 2005) and are less susceptible to atmospheric interference from smoke aerosols and thin clouds (Avery and Berlin 1992, Eva and Lambin 1998, Roy et al 1999. Here we use NBR as a complement to NDVI to track crop phenology with variations in vegetation greenness. We estimate the timing of crop maturation, or maximum greenness, during the monsoon growing season with both the daily median NDVI and NBR time series. Assuming that the seasonal progression in the crop cycle is similar across years, the timing of peak greenness in the growing season diagnoses the timing of the overall growing season. To estimate the timing of the maximum monsoon greenness with the noisy daily time series, we apply weighted cubic splines smoothing with bootstrapping on time steps within a defined window that straddles the day of monsoon peak greenness. Cubic splines smoothing stitches together piecewise third-order polynomial interpolation between 'knots,' or selected experimental points, and has been used extensively for crop phenology applications (Jain et al 2013, Mondal et al 2014. We apply weights to the NDVI and NBR time series using the daily fraction of 'usable' pixels, or those uncontaminated by clouds or thick haze (hereafter referred to as usable fraction) in the study area. This weighting follows from our greater confidence in daily median NDVI and NBR on clearer days versus cloudier and/or hazier days. Prior to bootstrapping, we make initial guesses of the four local maxima and minima: monsoon and winter peak greenness and pre-monsoon and post-monsoon trough greenness. We use these initial guesses to center a window of 300 d. Using a smoothing parameter of 0.75, we smooth the vegetation index time series with weighted cubic splines within the defined window and estimate the bootstrapped mean timing of maximum NDVI or NBR for each year. We repeat this process to estimate the earliest date when fields are ready to sow the winter crop, or trough greenness, during the post-monsoon transition period.

Regional aerosol optical depth exceedances
To quantify enhancements in regional air quality degradation during the post-monsoon burning season, we use MODIS/Terra Deep Blue retrievals of aerosol optical depth (AOD) over Punjab, Haryana, Delhi, and western Uttar Pradesh (i.e. encompassing the aerosol source and downwind transport regions of the IGP) (Levy et al 2013; appendix S1.3). Previous studies have shown that during the post-monsoon burning season, AOD co-varies with FRP, visibility, and PM, thus indicating the level of surface air quality degradation across north India (Vadrevu et al 2011. In order to minimize the contribution of background AOD, we analyze regionally averaged AOD 'exceedances'-that is, the daily spatial mean of AOD increments above the mean AOD +1σ for each pixel and season across Punjab, Haryana, Delhi, and western Uttar Pradesh. We analyze these daily mean AOD exceedances within the k start (FRP) and k end (FRP) window to isolate the effect of agricultural burning. To estimate the timing of peak AOD exceedances, or k peak (AOD), we apply Gaussian density curve optimization to values within this window expanded by two weeks. Such expansion ensures that the optimization is not thrown off by high AOD days isolated at the beginning or end of the season.

Trends in seasonal agricultural fire activity
The bimodal distribution of peak agricultural fire activity in both pre-monsoon and post-monsoon periods is limited to northwestern India, primarily in Punjab, as well as northern Haryana (figures 1 and S1). Generally, 80% of post-monsoon fires in Punjab are set within an approximate three-week window (23 ± 3 d) from mid-October to early November. We  (table S3). In contrast, we generally find no consistent statistically significant delays in the pre-monsoon burning season in Punjab (table S2).
Spatially, the post-monsoon temporal shift is larger in magnitude in districts in western Punjab than in eastern Punjab ( figure S3). Moreover, the 14year trends in total fire intensity for each 3 d block within this window signal a shift in the peak burning period, with decreasing FRP in mid-to-late October and increasing FRP in early November (figure 2). The magnitude of peak post-monsoon fire activity, indicated by the maximum 3 d block sums of FRP, has more than doubled over the 14-year period, an increase that may be partly attributed to some homogenization in the timing of burning across districts.

Trends in vegetation greenness from monsoon to post-monsoon
We also examine whether vegetation greenness in Punjab show similar shifts during the monsoon growing season and post-monsoon harvest-to-sowing transition period. Whereas the timing of minimum NBR and NDVI occurs after near-completion of post-monsoon burning in mid-to-late November, the temporal maximum of these vegetation indices occurs near the end of the monsoon around late August or early September ( figure 1(b)), indicating crop maturation. In Punjab, the timing of maximum NDVI and NBR shows an overall delay of 11-15 d, with a large, abrupt shift of 7-9 d around 2008-09 (figures 3(a) and (b)). Concurrently, there is an evident increasing trend in maximum monsoon NBR (0.06 decade −1 , 95% CI: [0.04, 0.08]) and NDVI (0.07 decade −1 , 95% CI: [0.05, 0.09]), consistent with steady increases in annual total kharif rice production in Punjab of 0.13 Tg yr −1 (95% CI: [0.09, 0.17]) (figures 3(b) and S5, table S5). Such increases in peak NBR and NDVI also suggest greater quantities of crop residue, which may lead to amplified fire intensity and emissions. In contrast to the shift in maximum NBR and NDVI, we find a smaller delay of 4-6 d in the timing of the minimum values of these indices during post-monsoon (figures 3(c) and (d), table S5), indicating that the shift in the monsoon growing season is greater than the corresponding shift in the timing of the earliest date when fields are ready for winter wheat sowing. In addition, we find that the duration from the start of the burning season to trough post-monsoon greenness has decreased by 0.71 d yr −1 (95% CI: [−1.02, −0.39]), providing evidence for a shortened harvest-to-sowing period (figure S6). Taken together, our results suggest that the temporal shifts in post-monsoon burning are likely associated with later sowing and harvesting of the monsoon crop.

The utility of NBR as a vegetation index
We have so far considered NBR and NDVI as complementary vegetation indices. Here we further demonstrate the utility of NBR for tracking crop phenology, particularly in resolving the troughs of the crop cycle. The weaker detrended correlations (r = 0.23 ± 0.39) between the two vegetation indices during transition months between the kharif and rabi seasons (May, June, October, and November) compared to other months (r = 0.88 ± 0.12) support the notion that NDVI more poorly resolves and tends to 'flatten' the troughs of the double-crop cycle curve ( figure S7). Moreover, the monthly distributions of detrended r(NDVI, NBR) values closely follow variations in greenness in the double-crop cycle, with greater correlation during seasons of crop growth. This pattern of correlation suggests that the performance of NDVI depends on the level of greenness in-field and that NDVI values at or near-minimum greenness should be interpreted Here exceedances are defined as the spatially averaged AOD increments above the mean AOD +1σ for each season and pixel. The black line shows the timing of peak daily AOD exceedances. The red lines represent the timing of the start, peak, and end of the post-monsoon burning season, based on daily FRP (same as in figure 2). The text inset shows the linear trend in the k peak (AOD), statistically significant at the 95% confidence level. Example of thick haze over the western IGP on November 6, 2016, as observed by MODIS/Terra, shown as (b) true color (NASA/Worldview; https://worldview.earthdata.nasa.gov/) and (c) Deep Blue AOD at 550 nm from the MOD04_L2 product. with caution. However, without further investigation, the noisiness and saturation of NDVI should not be generalized to all other vegetation indices that rely partly on visible bands, including the enhanced vegetation index (EVI) and soil-adjusted vegetation index (SAVI).

Trends in post-monsoon regional aerosol optical depth
To quantify the consequences of the delays in postmonsoon agricultural fire activity for regional air quality, we assess AOD exceedances during the main burning period bounded by k start (FRP) and k end (FRP). Within this window, post-monsoon AOD exceedances have increased by 55% from 2003 to 2016, likely associated with the reported upward trend in fire intensity (figure 4). Similar to the magnitude of the delay in k peak (FRP), the timing of the peak in AOD, k peak (AOD), has shifted by 0.8 d yr −1 (95% CI: [0.47, 1.13]), or ∼11 d during the 14-year period. The delay and increase in post-monsoon agricultural fire activity appear to drive the coherent shifting pattern in heavy aerosol loading episodes (higher AOD exceedances), notably observed in early November after 2008, despite the variability in AOD impacted by meteorology and other pollution sources, such as fireworks during the Diwali festival. Diwali lasts several days, and its timing is highly variable from year to year (October-November), following the lunar calendar. On average, the absolute difference between the timing of Diwali and the peak post-monsoon burning date is 7 d, with a range of 1-20 d ( figure S4). While Diwali fireworks, concentrated over dense population centers, can severely exacerbate urban pollution for a few days, crop residue burning contributes to widespread regional pollution for weeks (Mukherjee et al 2020, Ojha et al 2020.

Implications of delays in post-monsoon fire activity
We find that the peak fire intensity of the postmonsoon burning season in Punjab has shifted later in time by over two weeks from 2003 to 2016, with a 42% increase in overall fire intensity. This delay is gradual, likely influenced by steady increases in crop production and mechanization, which yield higher amounts of excess crop residue. We hypothesize that a shortened harvest-to-sowing turnaround time after kharif rice harvests has amplified this increase by making it difficult for farmers to prepare fields for timely sowing of rabi wheat. The optimal time to sow wheat in Punjab is late October to early November (Balwinder-Singh et al 2016, Liu et al 2019, yet co-occurring post-monsoon fires indicate that fields are often not ready at this time, particularly in recent years. Since fire is a quick and cheap method to remove the leftover residue generated by combine harvesters, farmers may have even greater incentive to burn crop residue, especially if harvests are delayed past the optimal date to sow wheat. Consistent with this hypothesis, we find that high fire intensity days preferentially occur during the latter half of the fire season, when the optimal window for sowing is shrinking. As postmonsoon fires increase in response to mechanization and pressures to sow on time, the burning season gradually trends later, further compressing the harvest-to-sowing window and increasing fire intensity rates. As a result, winter wheat sow dates across the region will likely homogenize, collapsing around a small optimal window to mitigate crop losses from increasing temperatures from February to March (Lobell et al 2012).
Additionally, we estimate a 55% increase in regional AOD exceedances and ∼11 d delay in the timing of peak AOD within the post-monsoon burning period from 2003 to 2016. Delays in the post-monsoon burning season also suggest that high fire activity periods may increasingly coincide with late-autumn/winter meteorological conditions that favor severe fog/smog and haze events across the IGP (Dey 2018). Dense fog formation peaks in winter (December to January) over the IGP (Ghude et al 2017, Dey 2018, Gautam and Singh 2018, but in recent years there appears to be an increasing tendency in dense fog episodes observed earlier in November, coinciding with the buildup of intense smoke associated with crop residue burning activity (figure S8). Further analysis of the long-term trends in meteorological conditions that are conducive to fog, such as western disturbances (Gautam et al 2007), is prerequisite to assessing possible feedbacks between meteorology and smoke from late October to November. Aside from increasing exposure to high regional PM concentrations both locally and in urban centers downwind, crop residue burning depletes soil moisture and decreases roadside visibility (Badarinath et al 2006, Sidhu et al 2015, Sinha et al 2015. In spite of bans, such burning continues to persist and gain traction (Tallis et al 2017). New technology that simultaneously reuses crop residue as mulch cover and incorporates seeds into the bare soil has been tested as an alternative to slash-and-burn methods of managing crop residue (Sidhu et al 2015, Tallis et al 2017, Shyamsundar et al 2019.

Potential drivers of delays in the rice-wheat rotation
Delays in the post-monsoon burning season are consistent with such shifts in the timing of monsoon peak greenness (11-15 d) and post-monsoon trough greenness (4-6 d), though the latter is of lesser magnitude. Unlike the steady shifts seen in post-monsoon burning, an abrupt delay of roughly one week occurring around 2008-09 dominates the overall delay in the timing of monsoon peak greenness, with relatively little change thereafter. Abrupt delays of similar magnitude are also apparent in the timing of the start of the post-monsoon burning season. Here we consider whether policy changes implemented around this time may have contributed toward these abrupt shifts. In 2009, in order to counteract severe groundwater depletion driven by low monsoon rainfall and widespread agricultural intensification, the Government of Punjab enacted the 'Preservation of Sub-Soil Water Act' (ordinance in 2008), which prohibits sowing rice nurseries before May 10 and transplanting the resulting rice seedlings to flooded paddies before June 10 (Ramanathan et al 2005, Singh 2009, Tripathi et al 2016, Asoka et al 2017. The Act delays the onset of water-intensive agricultural practices that would otherwise coincide with warm temperatures and high pre-monsoon evapotranspiration rates, which lead to excessive usage of the groundwater supply from tube wells and other reservoirs (Humphreys et al 2010).
Another policy that could be related to the shift is the 2008 all-India implementation of the Mahatma Gandhi National Rural Employment Guarantee Act (MGNREGA), a measure that provides a social security net to rural workers (Reddy et al 2014) and may have decreased the seasonal migration of workers to Punjab and led to labor shortages there (Singh 2009). Such shortages may have delayed the sowing of rice and incentivized use of combine harvesters, which may in turn explain the increase in crop residue burning. However, the already widespread transition to mechanized harvesting in Punjab, with diminishing dependence on manual labor, suggests that MGNREGA may have had a smaller impact on the timing of harvest and burning. Finally, variations in the timing of monsoon onset and withdrawal may be partly responsible for the interannual variability in these observed shifts, such as the early monsoon onset and rice maturation in 2013, but do not appear to drive the overall one-week delay in peak monsoon greenness from the 2003-2007 to 2008-2016 time periods (figure S9). It is important to note that here we do not establish direct causality with the groundwater policy, MGNREGA, or monsoon rainfall variability, but suggest a relationship that needs to be further explored in the field. Figure S10 summarizes the potential drivers and implications of the delay in and amplification of post-monsoon fire activity associated with double-crop cycle. Although this study focuses on Punjab, similar temporal shifts may have also occurred in adjoining states, such as Haryana, that experienced both severe groundwater depletion and increases in rice production and postmonsoon fire activity (Asoka et

Conclusion
In summary, we show robust, statistically significant temporal shifts of over two weeks in the timing of peak fire activity during the post-monsoon burning period in Punjab over a 14-year period from 2003 to 2016, and smaller delays of 11-15 d in monsoon peak greenness and 4-6 d in post-monsoon trough greenness. We estimate the timing of peak FRP and regional AOD exceedances by optimizing the Gaussian mean and the start, midpoint, and end of the burning season by using the partial sums of FRP. We further demonstrate the viability and applicability of using daily MODIS surface reflectance to characterize crop cycles and the utility of NBR as a useful complement to NDVI for quantifying these vegetation changes. We hypothesize that while the gradual delays in the post-monsoon burning season are likely linked to agricultural intensification and increasing mechanization, the abrupt delay of one week around 2008-09 seen in the monsoon crop growing season appears to coincide with the state-wide groundwater policy. The unintended consequences of these temporal shifts in the doublecrop cycle may be severe. First, a shortened harvestto-sowing period may further encourage farmers to burn crop residues in order to sow winter wheat on time. Second, the timing of peak crop residue burning may increasingly coincide with winter meteorology that favors severe smog events downwind across the IGP, where we diagnose a 55% increase in AOD exceedances, defined as the increment of AOD above the mean +1σ, over 2003-2016. For Punjab, alternative technology that combines the co-benefits of incorporating wheat seeds with rice residue and eliminating the need to burn residue, as well as switching to less water-intensive and stubble-producing crops, may alleviate the double bind of having to conserve groundwater while reducing public health exposure to smoke from post-monsoon fires.

Data Availability
All satellite-derived data used in this study are publicly available. MODIS-derived datasets can be accessed through NASA Earthdata (https:// search.earthdata.nasa.gov/) and Google Earth Engine (Gorelick et al 2017) (https://earthengine.google. com/). The Global Fire Emissions Database, version 4s, (GFEDv4s) and MODIS and VIIRS active fire geolocations are available from GFED (www.globalfiredata.org/), University of Maryland (http://fuoco.geog.umd.edu/), and NASA Fire Information for Resource Management System (FIRMS) (https://firms.modaps.eosdis.nasa. gov/). The analyzed data that support the findings of this study are available upon reasonable request from the authors.
The data that support the findings of this study are available upon reasonable request from the authors.