Environ Monit Assess (2014) 186:4531–4542 DOI 10.1007/s10661-014-3717-6

Geo-spatial analysis of the temporal trends of kharif crop phenology metrics over India and its relationships with rainfall parameters Abhishek Chakraborty & M. V. R. Seshasai & V. K. Dadhwal

Received: 8 September 2013 / Accepted: 6 March 2014 / Published online: 29 March 2014 # Springer International Publishing Switzerland 2014

Abstract The Global Inventory Modeling and Mapping Studies bimonthly Normalized Difference Vegetation Index (NDVI) data of 8×8 km spatial resolution for the period of 1982–2006 were analyzed to detect the trends of crop phenology metrics (start of the growing season (SGS), seasonal NDVI amplitude (AMP), seasonally integrated NDVI (SiNDVI)) during kharif season (June to October) and their relationships with the amount of rainfall and the number of rainy days over Indian subcontinent. Direction and magnitude of trends were analyzed at pixel level using the Mann– Kendall test and further assessed at meteorological subdivision level using field significance test (α=0.1). Significant pre-occurrence of the SGS was observed over northern (Punjab, Haryana) and central (Marathwada, Vidarbha and Madhya Maharashtra) parts, whereas delay was found over southern (Rayalaseema, Coastal Andhra Pradesh) and eastern (Bihar, Gangetic West Bengal and Sub-Himalayan West Bengal) parts of India. North, west, and central India showed significant increasing trends of SiNDVI, corroborating the kharif food grain production performance during the time frame. Significant temporal correlation (α=0.1) between the rainfall/number of rainy days and crop phenology metrics was observed over the rainfed region of India. About 35–40 % of the study area showed significant correlation between the SGS and the A. Chakraborty (*) : M. V. R. Seshasai : V. K. Dadhwal National Remote Sensing Centre, Indian Space Research Organization, Balanagar, Hyderabad 500037, India e-mail: [email protected]

rainfall/number of rainy days during June to August. June month rainfall/number of rainy days was found to be the most sensitive to the SGS. The amount of rainfall and the number of rainy days during monsoon were found to have significant influence over the SiNDVI in 24–30 % of the study area. The crop phenology metrics had significant correlation with the number of rainy days over the larger areas than that of the rainfall amount. Keywords Phenology . GIMMS NDVI . Median trend . Rainfall . Rainy days

Introduction Crop phenology is a study of recurring pattern of crop growth and development. The changes in crop phenology have long been recognized as being the most important early indicator of the impact of climate change on ecosystem (White et al. 1997). Variations in crop phenology have broad impacts on terrestrial ecosystems and human societies by altering global carbon, water and nitrogen cycles, crop production, duration of pollination season and distribution of diseases, etc. (Penuelas and Fiella 2001). Therefore, phenological study has recently become an important focus for ecological and climatic research (Menzel et al. 2001; Cleland et al. 2007). Monitoring of vegetation phenological events can be performed by three possible ways, i.e., in situ observations, bioclimatic models, and remote sensing techniques (Schaber and Badeck 2003; Fisher et al. 2007).

4532

Field-based observations of phenological events are expensive, time consuming, subjected to uncertainty due to operators bias, and thus not feasible to extend to large geographical areas. The bioclimatic models are specific to crop species and require accurate and exhaustive vegetation and climatic data along with local calibration to extend over larger scale. The satellite remote sensing technique provides a unique perspective of the planet and allows for regular, even daily/sub-daily, monitoring of the entire global land surface. As the most frequently used satellite sensors for monitoring phenological events have relatively large “footprints” on the land surface, they gather data about the entire ecosystems or regions rather than individual species. Remote sensing of phenology can reveal broad-scale phenological trends that would be difficult, if not impossible, to detect from the ground. The satellite data, after standardization and calibration across the sensors, could provide continuous information over large temporal span. Thus, the remote sensing techniques have emerged as low-cost alternatives of phenological studies with good temporal repeatability over large and inaccessible regions. Extensive research works were conducted to find out the temporal trends of phenology metrics using remote sensing techniques. Crop phenology is usually addressed through temporal monitoring of the satellitederived Normalized Difference Vegetation Index (NDVI) time series, available since the launch of the Earth Observation Satellites (EOS). NDVI is calculated as the normalized difference between near-infrared and red bands, and quantifies the photosynthetic capacity of plant canopies, thus, indicating the amount of vegetation present in one place. Coarse resolution sensor cannot provide directly the phenological events; rather, different phenology analogues such as the start of the growing season (SGS), the end of the growing season (EGS), the length of growing season (LGS), the seasonal NDVI amplitude (AMP), the seasonally integrated NDVI (SiNDVI), etc. could be derived. The time series NDVI had been used to detect increasing trends of the LGS over North America and Eurasia (Zhou et al. 2001) and changes of vegetation activities over Inner Mongolian steppe (Lee et al. 2002), Europe (Stockli and Vidale 2004), and China (Piao et al. 2006), Central Asia (De Beur and Henebry 2004; Kariyeva and van Leeuwen 2011), and Sahel and Soudan (Heumann et al. 2007). Apart from NDVI, vegetation phenology was also been addressed by satellite-derived chlorophyll index (Dash et al. 2010).

Environ Monit Assess (2014) 186:4531–4542

Many research results also showed that phenology is strongly linked to the climatic factors (rainfall and temperature) along with land cover changes and anthropogenic activities. Zhao and Schwartz (2003) found that the changes in the onset of spring were highly correlated with changes in temperature in Wisconsin, USA, and the findings corroborated by Estrella et al. (2007) for varieties of agricultural crops. Significant correlation was observed between NDVI and eco-climatic parameters in China, and that NDVI–growing degree-days correlation was found to be stronger than NDVI–rainfall correlation (Li et al. 2002). The responses of vegetation to changes in the climatic parameters over the northern mid- and high-latitudinal zones were analyzed in many previous research studies (Gong and Shi 2003; Slayback et al. 2003). The increase in plant growth was found to be associated with the lengthening of the active growing seasons over the period of 1981–1991 (Myneni et al. 1997). The rise in temperature was found to cause an increase in NDVI over the northern mid- and highlatitudinal zones (Ichii et al. 2002). The yearly maximum NDVI was found to have high correlation with temperature at high latitudes of the northern hemisphere (Qi 1999). Mao et al. (2012) found that trends of NDVI during 1982–2009 over northeast China was spatially heterogeneous and correspond to regional climatic characteristics of different season. The NDVI–temperature correlation was found to be stronger than the NDVI– rainfall during the study period. Limited studies were found on this issue over parts of the India region. Sehgal et al. (2011) derived the trends of phenology metrics over the Indo-Gangetic Plain (IGP) using time series AVHRR NDVI data and found a significant increasing trend of the maximum seasonal NDVI values both for kharif and rabi (November to April) season. The time of occurrence of peak NDVI during kharif season preoccurred over the northern part and delayed in the eastern part of IGP. The spatial patterns of vegetation phenology and its relation with climatic controls over the eight contrasting forest type over India were studied by Prasad et al. (2007). The study identified precipitation as a major control over temperature in the tropical forest. There had been substantial changes in the agricultural crop production and cropping pattern in India after the green revolution by the introduction of high-yielding varieties, chemical fertilizers, pesticides, irrigation infrastructures, land reforms, etc. Further, the monsoon rainfall patterns have also altered drastically due to climate change and global warming (Das et al. 2013). These has

Environ Monit Assess (2014) 186:4531–4542

a cascading effect on the time of crop sowing, crop vigor, and greenness, i.e., crop phenology. This study was a complete enumeration of the changes of kharif crop phenology over India involving an analysis of long-term Global Inventory Modeling and Mapping Studies (GIMMS) NDVI data of 1982–2006. The temporal trend of crop phenology metrics was calculated at pixel level to study the spatial patterns of the trend and also assessed at aggregated level over the meteorological subdivisions. A spatial analysis of the significant correlation between the amount/distribution of rainfall and the crop phenology metrics was also performed to find out the areas having significant control of rainfall towards crop phenology.

Study area and methodology The study was carried out over 25 meteorological subdivisions of India Meteorological Department (IMD) spread across 16 states of India (Fig. 1). The selected meteorological subdivisions have large contiguous areas of agricultural land and crop phenology during kharif Fig. 1 Study area with six different calibration zones. Meteorological subdivisions of the present study: 1 Punjab, 2 Haryana Chandigarh Delhi (Haryana), 3 West Uttar Pradesh (West UP), 4 East Uttar Pradesh (East UP), 5 West Rajasthan, 6 East Rajasthan, 7 West Madhya Pradesh (West MP), 8 East Madhya Pradesh (East MP), 9 Gujarat region Daman Dadra Nagar Haveli (Gujarat), 10 Sourashtra Kutch Diu (Sourashtra), 11 Chhattisgarh, 12 Jharkhand, 13 Bihar, 14 SubHimalayan West Bengal (SHWB), 15 Gangetic West Bengal (Gangetic WB), 16 Orissa, 17 Telangana, 18 Coastal Andhra Pradesh (Coastal AP), 19 Rayalaseema, 20 Madhya Maharastra, 21 Marathwada, 22 Vidarbha, 23 North Interior Karnataka (NI Karnataka), 24 South Interior Karnataka (SI Karnataka), 25 Tamil Nadu and Pondicherry (Tamil Nadu)

4533

season could easily be indentified using coarse resolution time series satellite data. Different steps followed in the present study are provided as below. Processing of NDVI time series The AVHRR GIMMS bimonthly NDVI dataset of 8 km spatial resolution for the period of 1982 to 2006 (Tucker et al. 2005) was analyzed. The dataset had correction against residual sensor degradation, inter-sensor calibration differences, distortions caused by persistent cloud cover globally, solar zenith angle and viewing angle effects due to satellite drift, etc. The maximum value composite of 15 days was used to reduce cloud cover while maintaining temporal resolution for accurate phenological estimates. The TIMESAT software was used to generate smooth time series of NDVI as well as to estimate the vegetation phenology for the study area (Jönsson and Eklundh 2002, 2004). Agricultural practices over the vast India subcontinent differs greatly in terms of its time of sowing, cropping pattern, cropping intensity, crop management, etc. So, a single equation to fit NDVI temporal

4534

curve may not hold good across the Indian subcontinent. Thus, the study area was divided into different calibration zones, based on NDVI time series principal component analysis (TS-PCA), to select smoothening techniques and its associated parameters settings specific to the zones (Heumann et al. 2007). The potential kharif agricultural area mask was prepared from multiyear land use land cover map using Advanced Wide Field Sensor (AWiFS) data. To cope up with the difference in spatial resolution between GIMMS NDVI (8 km) and AWiFS (60 m), majority rule approach was adopted to aggregate the potential kharif agricultural area mask to coarse spatial resolution of GIMMS NDVI. The mask was applied to extract kharif agricultural area over India. The NDVI images of the first fortnight of each month of the year 1984, 1995, and 2006 (probable normal year from each decade) were selected to perform TS-PCA. The first three TS-PCA components explained 96.8 % of variance. Six zones for calibration (Fig. 1) were identified iteratively by isodata clustering technique using these first three TS-PCA components. The Savitzsky–Golay filter (Chen et al. 2004) was used with different seasonality parameters, envelope adaptation strengths, and window sizes to smoothen of temporal NDVI profile over the calibration zones. The smoothened NDVI and the original NDVI of total time series were compared with ordinary least square techniques and majority of the area show R2 >0.9, significant at α= 0.05. Thus, the smoothening technique did not cause significant change to the original NDVI values.

Extraction of crop phenology metrics Three phenology metrics, i.e., SGS, AMP, and the seasonal greenness or SiNDVI of kharif season, were extracted from the NDVI time series using TIMESAT. The SGS was estimated as the time point at which the value of the fitted function had first exceeds 20 % of the distance between minimum and maximum at the rising limb. Likewise, the End of the Growing Season (EGS) was determined as the time point at which the value of the fitted function had first exceeded 20 % of the distance between minimum and maximum at the falling limb. The AMP was computed as the difference between the base and maximum NDVI values. The SiNDVI was defined as the area under the vegetation curve during the growing season and above the mean of the base value. The detailed descriptions of the

Environ Monit Assess (2014) 186:4531–4542

extraction of crop phenology matrics were discussed in Jönsson and Eklundh (2002). Calculation of temporal trends of crop phenology metrics The anomalies of the crop phenology metrics from its long-term mean were calculated and trends of these time series anomalies both at pixel and meteorological subdivision level were analyzed. The nonparametric Mann– Kendall test (Mann 1945; Kendall 1975) was used to detect presence of significant trends (α=0.1) of phenology metrics, while the Sen’s method (Sen 1968) was used to compute the magnitude of the trend. The detailed descriptions of the methodology for performing the Mann–Kendall test and estimation of Sen’s slope were provided in the Online Resources. For each pixel, and separately for the time series anomalies of SGS, AMP, and SiNDVI, the Mann– Kendall test (at α=0.1) was performed for the null hypothesis of no trend against the alternative hypothesis of increasing/decreasing trend. While analyzing phenology metrics during kharif season over a large region like India to identify trends, it was also important to ascertain the computed trend across the region using the field significance test (Livezey and Chen 1983). Field or global significance addresses the question of whether the number of independent or local tests for significance, reporting significant trends, could have occurred by chance. Field significance testing is important when applying multiple hypothesis tests for significance simultaneously to satellite-derived crop biophysical parameters (phenology) across a geographic extent, the “field” (Wilks 2006). The satellite-based NDVI data is often spatially correlated; thus, it increases the number of erroneous rejections of the local null hypothesis when performing multiple independent hypothesis tests (Wilks 2006). Field significance testing can determine whether the locations reporting rejection of the null hypothesis at the local level are statistically significant at the global level given spatial autocorrelation within the data. Further information on the field significance testing methodology is provided in the Online Resources (Das et al. 2013). Processing of gridded rainfall data The daily rainfall data of 24 years (1982–2005) generated by IMD at a grid size of 0.5° latitude × 0.5°

Environ Monit Assess (2014) 186:4531–4542

longitude were used (Rajeevan and Bhate 2009). The daily rainfall data were converted to the monthly cumulative rainfall data over Indian Summer Monsoon period (June to September). The monthly data were further cumulated over monsoon season to get total monsoon rainfall. The daily rainfall, if exceeds 2.5 mm, was considered as rainy day (Mahapatra and Mohanty 2006). The total number of rainy days in each monsoon month (June to September) and monsoon season as a whole were calculated for the period of 1982 to 2005. The time series anomalies of rainfall and number of rainy days were calculated using absolute deviation from its respective long-term average. Thus, four derivatives of rainfall were developed as the time series anomalies of monthly rainfall and rainy days, monsoon rainfall, and rainy days. Correlation between rainfall and crop phenology metrics The time series anomalies of rainfall derivatives and crop phenology metrics were brought to 10 km spatial resolution for the common period of 1982–2005. The temporal correlations between the rainfall parameters and crop phenology metrics over the different time periods were calculated using Pearson product moment. Further, spatial patterns of significant correlations at α= 0.1 were analyzed to find out the areas sensitive to rainfall.

Results and discussions Spatial patterns of temporal trends of crop phenology metrics Pixel wise temporal trends of the kharif SGS, AMP, and SiNDVI over different meteorological subdivisions were shown in Table 1. The numbers of pixels having trends (increasing/decreasing/ no trend) were presented as N in Table 1. The numbers of pixels showing significant trends (at α=0.1) both in case of increasing and decreasing were shown in brackets. Average magnitudes of the trend (the Sen’s method) separately for increasing and decreasing categories over each meteorological subdivision were presented as slope in Table 1. The meteorological subdivisions which passed the “field significance test” at α=0.1 were shown in italics.

4535

The spatial patterns of temporal trends of the kharif SGS was presented in Fig. 2a. Decreasing trend (>−0.3 days/year) of the SGS, i.e., pre-occurrence of kharif season, was observed over most of the part of north, central, and western India. Increasing trend of the SGS, i.e., delay in kharif season (>0.3 days/year), was found in eastern and southern parts of India. About 15.35 % of agricultural areas under the present study showed significant trend (increasing/decreasing) of the SGS, 10.60 % pertain decreasing and 4.75 % to increasing trend. The field significance test over the meteorological subdivisions revealed significant decreasing trends over Punjab (−0.77 days/year), Haryana (−0.57 days/year), West UP (−0.78 days/year), Marathwada (−0.75 days/year), Vidarbha (−0.61 days/ year), and Madhya Maharashtra (−0.1 days/year). Whereas, significant increasing trends were observed over Rayalaseema (0.94 days/year), Coastal Andhra Pradesh (1.24 days/year), Bihar (1.1 days/year), Gangetic West Bengal (0.94 days/tear) and SubHimalayan West Bengal (1.6 days/year). Spatial patterns of temporal trends of kharif seasonal greenness or SiNDVI showed (Fig. 2b) greening or increasing trend over large contiguous areas of north, central, western, and deccan plateau regions of India. Southern and eastern part showed decreasing or browning trend of the SiNDVI. About 25.75 % of the total agricultural area under the present study showed significant (increasing/decreasing) trends of SiNDVI, within which significant greening (increasing) and browning (decreasing) trends were found in 17.79 and 7.96 % of areas. The field significance test showed statistically significant increasing trends of SiNDVI (0.02 to 0.05 year−1) over Punjab, Haryana, West and East Uttar Pradesh, West and East Rajasthan, West Madhya Pradesh, Sourashtra, Rayalaseema, Marathwada, Vidarbha, Bihar, and Sub-Himalayan West Bengal. On the other hand, statistically significant decreasing trends of SiNDVI (−0.02 to −0.04 year−1) were found over Tamil Nadu, South Interior Karnataka, Coastal Andhra Pradesh, Madhya Maharashtra, Chhattisgarh, Gujarat, and Gangetic West Bengal. East Madhya Pradesh passed the field significant test both for increasing and decreasing categories. The northern part of East Madhya Pradesh showed significant increasing trend where as southern part gave rise significant decreasing trend of the SiNDVI. The temporal dynamics of all India kharif area and production of food grain during 1982–2006 were presented in Fig. 3. The kharif sown area remained

41 184 (6) 510 (40) 899 (58) 414 (48) 743 (36) 895 (36) 343 (10) 212 (37) 213 (30) 457 (36) 266 (15) 783 (215) 133 (50) 429 (118) 466 (78) 659 (74) 261 (116) 283 (97) 251 (26) 223 (11) 290 (15) 299 (55) 362 (79) 237 (52)

0.113 581 (275) 0.252 499 (262) 0.412 761 (323) 0.379 1,629 (653) 0.626 722 (234) 0.363 1,252 (561) 0.387 1,488 (1,043) 0.435 415 (188) 0.651 163 (54) 0.619 371 (163) 0.645 392 (126) 0.611 264 (54) 1.083 957 (547) 1.605 168 (96) 0.938 413 (115) 0.793 308 (87) 0.598 543 (151) 1.244 166 (72) 0.936 314 (103) 0.577 252 (25) 0.393 231 (52) 0.352 565 (250) 0.666 323 (127) 0.796 252 (64) 1.093 180 (56)

Slope N ( )

0.0032 635 (350) 0.021 768 (548) 0.0039 433 (261) 0.027 491 (183) 0.0028 653 (255) 0.019 955 (355) 0.0026 1,571 (856) 0.021 1,154 (109) 0.0029 600 (242) 0.029 711 (92) 0.0044 797 (410) 0.038 1,094 (111) 0.0066 1,381 (726) 0.038 939 (87) 0.0031 416 (143) 0.026 675 (72) 0.0032 129 (51) 0.040 193 (45) 0.0041 397 (184) 0.046 273 (31) 0.0022 244 (35) 0.019 424 (48) 0.0019 260 (37) 0.018 367 (47) 0.0044 724 (359) 0.026 365 (46) 0.0045 105 (34) 0.021 55 (7) 0.0024 221 (29) 0.011 283 (35) 0.0022 306 (41) 0.018 456 (68) 0.0020 534 (113) 0.022 454 (49) 0.0036 137 (38) 0.024 161 (25) 0.0028 296 (92) 0.027 131 (23) 0.0016 409 (79) 0.024 678 (163) 0.0019 393 (107) 0.028 534 (104) 0.0027 683 (407) 0.045 625 (144) 0.0026 310 (94) 0.024 385 (77) 0.0020 205 (50) 0.019 242 (49) 0.0023 86 (19) 0.019 220 (35)

N()

N() 230 (44) 196 (37) 745 (295) 439 (41) 453 (41) 609 (68) 365 (116) 638 (275) 278 (123) 166 (20) 549 (191) 407 (106) 232 (37) 21 (4) 339 (100) 660 (179) 609 (145) 355 (141) 113 (16) 767 (274) 526 (164) 352 (133) 348 (116) 391 (163) 418 (157)

−0.766 −0.567 −0.782 −0.456 −0.507 −0.445 −0.511 −0.638 −0.813 −0.517 −0.659 −0.706 −0.616 −0.498 −0.556 −0.946 −0.631 −0.948 −0.821 −0.999 −0.748 −0.613 −0.819 −0.721 −0.904

AMP Slope

Slope N ( )

N()

Slope

SGS

AMP

SGS

SiNDVI

Decreasing trend

Increasing trend

−0.0017 −0.0020 −0.0028 −0.0013 −0.0018 −0.0016 −0.0025 −0.0029 −0.0033 −0.0021 −0.0028 −0.0023 −0.0017 −0.0015 −0.0023 −0.0026 −0.0021 −0.0031 −0.0019 −0.0028 −0.0025 −0.0027 −0.0024 −0.0027 −0.0026

Slope 143 (16) 95 (16) 528 (120) 252 (33) 168 (14) 420 (67) 361 (67) 576 (153) 272 (145) 106 (12) 678 (212) 403 (85) 286 (49) 40 (7) 398 (122) 614 (138) 587 (109) 316 (112) 89 (12) 569 (181) 347 (54) 232 (57) 336 (93) 436 (175) 413 (193)

N()

SiNDVI

−0.009 −0.010 −0.015 −0.012 −0.013 −0.015 −0.023 −0.027 −0.040 −0.021 −0.034 −0.025 −0.021 −0.016 −0.018 −0.025 −0.022 −0.033 −0.021 −0.028 −0.025 −0.029 −0.025 −0.028 −0.033

Slope 2 20 46 27 577 33 27 63 519 280 161 209 179 97 103 376 129 368 262 498 83 19 413 335 872

N

0 0 5 12 527 9 8 28 483 229 101 171 129 96 63 330 90 269 249 408 83 17 426 296 731

N

33 167 330 257 934 653 119 89 523 263 120 179 308 140 196 378 121 337 291 449 100 19 451 298 830

N

SGS AMP SiNDVI

No trend

811 695 1,511 2,080 1,702 1,870 1,861 1,081 924 766 1,042 842 1,318 285 815 1,298 1,242 790 676 1,427 840 934 1,097 939 1,329

N

Total pixels

N number of pixels having increasing/decreasing/no trend; pixels showing significant trend at α=0.1 are presented in parentheses; average magnitude of trend (increasing or decreasing) estimated by the Sen’s method are presented as slope (days/year for SGS; year−1 for AMP and SiNDVI); meteorological subdivisions passing the “field significant test” at α=0.1 are marked as italics

Punjab Haryana West UP East UP West Rajasthan East Rajasthan West MP East MP Gujrat Sourashtra Chhattisgarh Jharkhand Bihar SHWB Gangetic WB Orissa Telangana Coastal AP Rayalaseema Madhya Maharastra Marathwada Vidarbha NI Karnataka SI Karnataka Tamil Nadu

Meteorological subdivision

Table 1 Pixel wise temporal trends of the start of growing season (SGS), the seasonal NDVI amplitude (AMP), and the seasonal integrated NDVI (SiNDVI) during kharif season using GIMMS NDVI data for the period of 1982–2006

4536 Environ Monit Assess (2014) 186:4531–4542

Environ Monit Assess (2014) 186:4531–4542

4537

Fig. 2 Spatial patterns of median trends (Sen’s slope) of a start of the growing season (SGS); b seasonally integrated NDVI (SiNDVI) and; c seasonal NDVI amplitude of kharif using GIMMS NDVI time series of 1982–2006

increasing trends over north, east, and central regions of India, whereas the southern part showed decreasing trends. It was found that 30.76 % of the total agricultural area under the study was having significant trend of the AMP, wherein 20.16 % of it pertaining to increasing and 10.60 % to decreasing trend of the AMP. The field significance test at α=0.1 found statistically significant increasing trends of the AMP (0.002 to 0.007 year−1) over Punjab, Haryana, West and East Rajasthan, East Uttar Pradesh, Bihar, Sub-Himalayan West Bengal, West Madhya Pradesh, Sourastra, North Interior Karnataka, and Rayalaseema. On the other hand, significant decreasing trends of the AMP (−0.002 to

almost same during the period of 1982–2006. The sown area oscillated around 80 Mha with two dips during two major drought incidences of 1987–1988 and 2002–2003. But there was a steady increase in total production of the food grain crop from 79.87 Mtons (1981–1982) to 109.87 Mtons (2005– 2006). It had occurred because of better plant growth and yield of the kharif crop leading to large-scale greening trend over India. Thus, the crop statistics also confirmed the greening trends over India as found from satellite observations. Spatial patterns of temporal trends of kharif season NDVI amplitude or crop vigor (Fig. 2c) showed

120

100 Production (Mtons)

110

80

Soybean area of Madhya Pradesh (%)

100 60 90 40 80 20

70

0

2005-06

2003-04

2001-02

1999-00

1997-98

1995-96

1993-94

1991-92

1989-90

1987-88

1985-86

1983-84

60

% of kharif agricultural area

Area (Mha)

1981-82

Area (M ha) / Production (M tons)

Fig. 3 Area, production of kharif food grains of India, and percentage of kharif agricultural area under soybean cultivation in Madhya Pradesh State during 1982–2006. (Source: Directorate of Economics and Statistics, Department of Agriculture and Cooperation, India; The Soybean Processor Association)

4538

−0.003 year−1) were observed over Tamil Nadu, South Interior Karnataka, Coastal Andhra Pradesh, Marathwada, Madhya Maharastra, Orissa, Gujarat, and Jharkhand. The field significance test was found to be passed both for increasing and decreasing trend in case of Telangana, Vidarbha, Chhattisgarh, Gangetic West Bengal, East Madhya Pradesh, and West Uttar Pradesh indicating that part of these meteorological subdivisions possess both significant increasing and decreasing trends of kharif seasonal NDVI amplitude. The significant increasing trends of crop vigor (AMP) also leads to greening trend and the statistics of kharif food grain production (Fig. 3) supported this observation. It is important to note here that the greening trend could occur because of better crop growing environment such as irrigation/rainfall, fertilization, crop protection, etc. In addition to this, introduction of a new crop with more luxurious growth may also cause change in the trend of AMP. Very high increasing trend in the AMP (0.0066 year−1) was observed over West Madhya Pradesh. Majority of the areas of Madhya Pradesh consist of vertisol or black soil. The general practices over there was to keep the vertisol fallow during the rainy season (kharif) and take the crop on residual soil moisture in post rainy season (rabi) season because of the difficulties in soil preparation before monsoon for timely sowing of kharif crop. The soybean crop (Glycine max) was introduced successfully during kharif season particularly in western part of Madhya Pradesh in early 1980s. Thus, kharif fallow and other competing crop like Bengal gram (Arachis hypogaea) with low canopy cover were replaced by soybean. Figure 3 showed the temporal dynamics of percentage of soybean crop to total kharif agriculture area of Madhya Pradesh State. The soybean crop consisted of less than 10 % of the total kharif agricultural area in Madhya Pradesh before 1984. It had steadily increased to more than 50 % by 1997 and reached a plateau afterwards. The soybean crop has very high canopy cover with luxurious growth and it has contributed in substantial increase in vegetation fractions over the area during kharif season leading to high increasing trend of crop vigor (AMP) over the West Madhya Pradesh. Spatial patterns of correlation between crop phenology and rainfall The changes in crop phenology may be attributed to different anthropogenic activities such as introduction of

Environ Monit Assess (2014) 186:4531–4542

irrigation, change in crop type/variety, crop management, etc. in addition to variation in climatic parameters. In temperate zones, biological activities and crop growth are mostly controlled by temperature (Menzel 2002), while at the lower latitude, rainfall and evapotranspiration are the main drivers (Linderholm et al. 2005). India being in the subtropical region, monsoon rainfall has major control over the agricultural production. About 2/3 of the Indian agricultural land comes under rainfed condition. Crop growth and development in these areas is highly dependent on the amount and temporal distribution of rainfall. Thus, crop phenology during the kharif season is mainly dependent on the monsoon rainfall pattern.

Start of the growing season and rainfall The growing season starts with the crop sowing operation. The SGS is influenced by the soil wetness and its persistence. Thus, rainfall amount and its distribution play a key role towards the SGS. Indian summer monsoon starts in the month of June. Rainfall activities during June and the subsequent months decide the crop sowing operation particularly in the rainfed situation. Any deficiency of rainfall during the early monsoon season delays the crop sowing activities. So negative anomalies of rainfall and number of rainy days during early part of the monsoon season would cause positive anomalies, i.e., delay in the SGS and vice versa. It was observed that substantial area (26 %) under the present study was found to have significant negative correlation between the anomalies of June month rainfall and the SGS (Fig. 4 and Table 2). Large contiguous areas over East and West Rajasthan, Gujarat, Sourashtra, Haryana, East and West Uttar Pradesh, East and West Madhya Pradesh, Vidarbha, Marathwada, and Telangana showed significant negative relationships. July and August rainfall were significantly correlated with the SGS in only 6 and 5.1 % of the study area, respectively. Cumulative rainfall of June–July and June–August were found to be correlated the SGS in 13 and 7.8 % of the study area, respectively. It was found that the areas with significant correlation with the SGS and rainfall over different time scale were not mutually exclusive. So, unions of these areas (35.3 % of the study area) had significant correlation between the SGS and rainfall of June or July or August or June–July or June–August.

Environ Monit Assess (2014) 186:4531–4542

SGS and June rainfall

SiNDVI and June/June-July /June-August/June-September rainfall

SGS and June/July /August /June –July /June -August rainfall

SiNDVI and June /June-July /June-August/June-September rainy days

4539

SGS and June rainy days

SGS and June/July /August /June –July /June -August rainy days

AMP and June/June-July /June-August/June-September rainfall

AMP and June/June-July /June-August/June-September rainy days

Fig. 4 Spatial patterns of significant correlation (α=0.1) between time series anomalies of the crop phenology metrics and rainfall/number of rainy days

The anomalies of the number of rainy days during June were found to have significant negative correlation with anomalies of the SGS over the 30.1 % of the study area (Fig. 4 and Table 2). The spatial patterns were found to be similar to that of the SGS and rainfall amount. The anomalies of the number of rainy days of July and August were found to have significant negative correlation with anomalies of the SGS over the 9.8 and 5.2 % of the area under the study, respectively. The correlations between the anomalies of cumulative number of rainy days during June–July and June–August with the anomalies of the SGS were also calculated; 21.5 and 13.5 % area under the study showed significant relationship, respectively. Almost 41.6 % areas of the present study showed significant negative correlation between the SGS and the number of rainy days during June or July or August or June–July or June–August. The bio-window of the SGS during kharif is from third week of June to second week of August with a peak in mid-July. June month rainfall/number of rainy days was found to be crucial as it showed significant correlation to the SGS over large contiguous areas in northwest, northern, and central regions of India. Rainfall

during the subsequent months might have importance for crop establishment and development but had less control over the SGS. This also confirmed the lag relationship between the rainfall and the SGS. The eastern region, west and east coast of India being high rainfall zones (>250 mm) did not show any significant relation. The anomalies of rainfall did not have significant effect on the SGS anomalies over these high rainfall zones. Seasonal greenness and rainfall The SiNDVI represents the seasonal greenness and seasonal vegetation activities. Any deficiency of rainfall has negative effect on the SiNDVI and vice versa provided there is no other source of water supply. The SiNDVI was calculated as integrated area under the vegetation curve above the mean of the base value (“Spatial patterns of correlation between crop phenology and rainfall” section) and thus SiNDVI was a cumulative property of crop phenology. So, the relationships between the SiNDVI and rainfall/number of rainy days cumulated over the monsoon period (June–September) were evaluated.

4540

Environ Monit Assess (2014) 186:4531–4542

Table 2 Percentage of the study area showed significant (α=0.1) correlation between the crop phenology metrics and rainfall parameters

any significant relation because of high average monsoon rainfall over there.

Crop phenology metrics

Time period

Rainfall Rainy days

Seasonal crop vigour and rainfall

Star of the growing season

June

26.0

30.1

July

6.0

9.8

August June-July June-August

Seasonal NDVI amplitude (crop vigor)

5.1

5.2

13.0

21.5

7.8

13.5

June/July/August/June– 35.3 July/June–August June 5.4

41.6

June–July

9.4

10.8

June–August

9.6

13.4

June–September

9.8

14.6

18.2

22.2

9.5

10.2

June/June–July/June– August/June– September Seasonally integrated June NDVI (seasonal June–July greenness) June–August

5.4

12.3

16.7

11.6

16.6

June–September

13.1

20.4

June/June–July/June– August/June– September

24.0

29.4

The rainfall/number of rainy days of June had significant correlation with the SiNDVI over 10 % areas under the study (Table 2). The areas with significant relationship had increased as cumulative rainfall/number of rainy days over June–July, June–August, and June– September were regressed with the SiNDVI. It was observed that 13.1 and 20.4 % of the study area showed significant correlation between the SiNDVI and rainfall/ number of rainy days during June–September, respectively. It was also seen that 24 and 29.4 % of the study area had significant correlation between the SiNDVI and the rainfall/number of rainy days during June or June–July or June–August or June–September, respectively (Fig. 4 and Table 2). The areas having significant relation were found to be concentrated over West and East Rajasthan, Gujarat, Sourashtra, West Madhya Pradesh, Haryana, East and West Uttar Pradesh, Vidarbha, Marathwada, Rayalaseema, and Telangana. These areas grossly covered the arid and semiarid regions of India with little to no irrigation support. Considerable portion of rainfed region in the eastern part of the country did not show

The AMP, computed as the difference between the base and maximum seasonal NDVI value, represents the seasonal crop vigor. The seasonal crop vigor may change due to change in crop variety/type, increase in crop fraction within the satellite footprint, and alteration of crop growth environment. Large proportion of agricultural lands showed significant increasing trend of the AMP (“Spatial patterns of temporal trends of crop phenology metrics” section) indicating better crop growth and vigor. It was observed that the maximum peak of the NDVI temporal profile occurred in the month of September in majority of the region. Thus, cumulative rainfall/ number of rainy days during June, June–July, June–August, and June–September were correlated with the AMP. It was observed that only 5.4 % of the area under the study showed significant correlation between the anomalies of June month rainfall/number of rainy days and the AMP (Table 2). The significant pixels were highly scattered and did not show any clustering. The areas having significant correlation had increased to 9.4 and 9.6 % of the study area in case of relationship between anomalies of the AMP and anomalies of the rainfall during June–July and June–August, respectively. Whereas, only 9.8 % of the area showed positive correlation between the anomalies of monsoon (June– September) rainfall and the AMP. On the other hand, 5.4 and 10.8 % of the area under the present study showed significant positive correlation between the anomalies of the AMP and number of rainy days during June and June–July, respectively. Further, 13.4 and 14.6 % of the area were found to have significant positive correlation between anomalies of AMP and rainy days during June–August and June–September (monsoon), respectively. Almost 18.2 and 22.2 % of the study areas were found to have significant correlation between the AMP and rainfall/number of rainy days of either June or June–July or June–August or June– September, respectively (Fig. 4 and Table 2). Cumulative rainfall and number of rainy days of June–August were found to have control on the AMP over the maximum areas under the study. There was marginal increase in area having significant correlation

Environ Monit Assess (2014) 186:4531–4542

between the AMP and the monsoon rainfall/number of rainy days. North-western and central Indian dry track along with eastern India covering West Bengal, Bihar, Jharkhand, and Orissa showed significant relationship.

Conclusions The present study utilized the time series GIMMS NDVI to find out the spatial patterns of temporal trends of kharif season crop phenology metrics over India. The significant influences of the amount of rainfall and the number of rainy days over the crop phenology metrics at different time lags were also analyzed. The satellite observations of 25 years detected significant shift in different crop phenology metrics. Significant preoccurrence of the SGS (0.1 to 0.7 days/year) was found in large contiguous areas of Punjab, Haryana, West Uttar Pradesh, Marathwada, Vidarbha, and Madhya Maharashtra. Delay in the SGS (0.9–1.6 days/year) was found in Rayalaseema, Coastal Andhra Pradesh, Bihar, Gangetic West Bengal, and Sub-Himalayan West Bengal. Two kinds of greening trends over Indian region during kharif season were observed. Increasing SiNDVI along with increase in the AMP was observed over Punjab, Haryana, West and East Uttar Pradesh, West and East Rajasthan, West and East Madhya Pradesh, Bihar, Sub-Himalayan West Bengal, Sourashtra and Kutch, and Rayalaseema. Whereas, increasing SiNDVI along with decrease in the AMP was found over Marathawada and Vidarbha implying increase of the length of the growing period. Significant browning trends during kharif season were observed in most of the south and eastern part of India covering Tamil Nadu, South Interior Karnataka, Coastal Andhra Pradesh, Madhya Maharashtra, Gujarat, Chhattisgarh, Jharkhand, and Gangetic West Bengal. Large contiguous areas of north, west, and central India were found to have significant correlation between the time series anomalies of rainfall/number of rainy days and crop phenology metrics. June month rainfall and number of rainy days were found to be more sensitive to the SGS with maximum area (25–30 %) under its influence. Almost 35–40 % of the study area showed significant correlation between the SGS and the rainfall/ number of rainy days during June to August month. The rainfall and the number of rainy days during monsoon period significantly controlled the SiNDVI over the 24–

4541

30 % of the area under the present study. Rainfall was found to have significant influence to the crop phenology metrics over the rainfed regions of India covering the arid and semiarid tracks. High rainfall region of India covering eastern India and east and west coast did not show any significant correlation between the rainfall parameters and crop phenology metrics. The temporal distribution of rainfall in terms of the number of rainy days was found to have significant relationship to crop phenology over the larger areas than that of the amount of rainfall. The present study demonstrated the changing patterns of Indian agriculture in terms of its crop calendar, cropping pattern, crop type, and net sown area along with its sensitivity towards Indian Summer Monsoon. The scope needs to be extended to quantify separately the causative factors (climatic and anthropogenic) towards the changes in crop phenology metrics. Acknowledgments IMD is duly acknowledged for providing gridded rainfall data. The TIMESAT software is available online in http://www.nateko.lu.se/personal/Lags.Eklundh/TIMESAT/ timesat. Tucker et al. 2004 is duly acknowledged for providing GIMMS NDVI dataset. Thanks to the anonymous reviewers and the editor for their valuable suggestions.

References Chen, J., Jönsson, P., Tamura, M., Gu, Z. H., Matsushita, B., & Eklundh, L. (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sensing of Environment, 91(4), 332–344. Cleland, E. E., Chuine, I., & Menzel, A. (2007). Shifting plant phenology in response to global change. Trends in Ecology and Evolution, 22(7), 357–365. Das, P. K., Chakraborty, A., & Seshasai, M. V. R. (2013). Spatial analysis of temporal trend of rainfall and rainy days during the Indian Summer Monsoon season using daily gridded (0.5° × 0.5°) rainfall data for the period of 1971–2005. Meteorological Applications. doi:10.1002/met.1361. Dash, J., Jeganathan, C., & Atkinson, P. M. (2010). The use of MERIS terrestrial chlorophyll index to study spatio-temporal variation in vegetation phenology over India. Remote Sensing of Environment, 114(7), 1388–1402. De Beur, K. M., & Henebry, G. M. (2004). Land surface phenology, climatic variation, and institutional change: analyzing agricultural land cover change in Kazakhstan. Remote Sensing of Environment, 89(4), 497–509. Estrella, N., Sparks, T. H., & Menzel, A. (2007). Trends and temperature response in the phenology of crops in Germany. Global Change Biology, 13(8), 1737–1747. Fisher, J. I., Richarson, A. D., & Mustard, J. F. (2007). Phenology model from surface meteorology does not capture satellite

4542 based greenup estimation. Global Change Biology, 13(7), 707–721. Gong, D. Y., & Shi, P. J. (2003). Northern hemispheric NDVI variation associated with large-scale climate indices in spring. International Journal of Remote Sensing, 24(12), 2559–2566. Heumann, B. W., Seaquist, J. W., Eklundh, L., & Jonsson, J. (2007). AVHRR derived phenological change in the Sahel and Soudan, Africa, 1982–2005. Remote Sensing of Environment, 108(4), 385–392. Ichii, K., Kawabata, A., & Yamaguchi, Y. (2002). Global correlation analysis for NDVI and climatic variables and NDVI trends, 1982–1990. International Journal of Remote Sensing, 23(18), 3873–3878. Jönsson, P., & Eklundh, L. (2002). Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Transactions on Geoscience and Remote Sensing, 40(8), 1824–1832. Jönsson, P., & Eklundh, L. (2004). TIMESAT—a program for analyzing timeseries of satellite sensor data. Computers & Geosciences, 30(8), 833–845. Kariyeva, J., & van Leeuwen, W. J. D. (2011). Environmental drivers of NDVI based vegetation phenology in Central Asia. Remote Sensing, 3(2), 203–246. Kendall, M. G. (1975). Rank correlation methods (4th ed.). London: Charles Griffin. Lee, R., Yu, F., Price, K. P., Ellis, J., & Shi, P. (2002). Evaluating vegetation phenological pattern in Inner Mongolia using NDVI time series analysis. International Journal of Remote Sensing, 23(12), 2505–2512. Li, B., Tao, S., & Dawson, R. W. (2002). Relations between AVHRR NDVI and ecoclimatic parameters in China. International Journal of Remote Sensing, 23(5), 989–999. Linderholm, H. W., Walther, A., Chen, D. (2005). Growing season trends in the greater Baltic area. University of Goteborg, C69 Rapport, Goteborg. ISSN 1400-383X. Livezey, R. E., & Chen, W. Y. (1983). Statistical field significance and its determination by Monte Carlo technique. Monthly Weather Review, 111, 46–59. Mahapatra, M., & Mohanty, U. C. (2006). Spatio-temporal variability of summer monsoon rainfall over Orissa in relation to low pressure system. Journal of Earth System Science, 115(2), 203–218. Mann, H. B. (1945). Non-parametric tests against trend. Econometrica, 13, 245–259. Mao, D., Wang, Z., Luo, L., & Ren, C. (2012). Integrating AVHRR and MODIS data to monitor NDVI changes and their relationships with climatic parameters in Northeast China. International Journal of Applied Earth Observation and Geoinformation, 18, 528–536. Menzel, A. (2002). Phenology, its importance to the global change community. Editorial comments. Climatic Change, 54(4), 379–385. Menzel, A., Estrella, N., & Fabian, P. (2001). Spatial and temporal variability of the phenological seasons in Germany from 1951 to 1996. Global Change Biology, 7(6), 657–666.

Environ Monit Assess (2014) 186:4531–4542 Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in northern high latitude from 1981 to 1991. Nature, 386(17), 698–702. Penuelas, J., & Fiella, I. (2001). Response to a warming world. Science, 294(5543), 793–794. Piao, S. L., Fang, J. Y., Zhou, L. M., Ciais, P., & Zhu, B. (2006). Variation in satellite derived phenology on China’s temperate vegetation. Global Change Biology, 12(4), 672–685. Prasad, V. K., Badarinath, K. V. S., & Anuradha, E. (2007). Spatial pattern of vegetation phenology metrics and related climatic controls of eight contrasting forest types of India—analysis from remote sensing datasets. Theoretical and Applied Climatology, 89(1–2), 95–107. Qi, Y. (1999). The effect of climate change on vegetation at high latitude of the northern hemisphere, a functional analysis. Acta Ecologia Sinica, 19(4), 474–477. Rajeevan, M., & Bhate, J. (2009). A high resolution daily gridded rainfall dataset (1971–2005) for mesoscale meteorological studies. Current Science, 96(4), 558–562. Schaber, J., & Badeck, F. W. (2003). Physiology based phenology models for forest tree species in Germany. International Journal of Biometeorology, 47(4), 193–201. Sehgal, V. K., Jain, S., Aggarwal, P. K., & Jha, S. (2011). Deriving crop phenology metrics and their trends using time series NOAA-AVHRR NDVI data. Journal of Indian Society of Remote Sensing, 39(3), 373–381. Sen, P. K. (1968). Estimation of regression coefficients based on Kendall’s tau. Journal of the American Statistical Association, 63, 1379–1389. Slayback, D. A., Pinzon, J., Los, S. O., & Tucker, C. J. (2003). Northern hemisphere photosynthetic trends 1982-1999. Global Change Biology, 9(1), 1–15. Stockli, R., & Vidale, P. L. (2004). European plant phenology and climate as seen in a 20 year AVHRR land surface parameters dataset. International Journal of Remote Sensing, 25(17), 3303–3330. Tucker, C. J., Pinzon, J. E., Brown, M. E., Slayback, D. A., Pak, E. W., & Mahoney, R. (2005). An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. International Journal of Remote Sensing, 26(20), 4485–4498. White, M. A., Thornton, P. E., & Running, S. W. (1997). A continental phenology model for monitoring vegetation response to interannual climatic variability. Global Biogeochemical Cycles, 11(3), 217–234. Wilks, D. S. (2006). On “Field Significance” and the false discovery rate. Journal of Climate, 45, 1181–1189. Zhao, T. T., & Schwartz, M. D. (2003). Examining the onset of spring in Wisconsin. Climate Research, 24(1), 59– 70. Zhou, L. M., Tucker, C. J., Kaufmann, R. K., Slayback, D., Shabanov, N. V., & Myneni, R. (2001). Variation in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999. Journal of Geophysical Research – Atmospheres, 106(D17), 20069–20083.