PHENOLOGY ANALYSIS OF FOREST VEGETATION TO ENVIRONMENTAL VARIABLES DURING PRE- AND POST-MONSOON SEASONS IN WESTERN HIMALAYAN REGION OF INDIA

To assess the phenological changes in Moist Deciduous Forest (MDF) of western Himalayan region of India, we carried out NDVI time series analysis from 2013 to 2015 using Landsat 8 OLI data. We used the vegetation index differencing method to calculate the change in NDVI (NDVIchange) during pre and post monsoon seasons and these changes were used to assess the phenological behaviour of MDF by taking the effect of a set of environmental variables into account. To understand the effect of environmental variables on change in phenology, we designed a linear regression analysis with sample-based NDVIchange values as the response variable and elevation aspect, and Land Surface Temperature (LST) as explanatory variables. The Landsat-8 derived phenology transition stages were validated by calculating the phenology variation from Nov 2008 to April 2009 using Landsat-7 which has the same spatial resolution as Landsat-8. The Landsat-7 derived NDVI trajectories were plotted in accordance with MODIS derived phenology stages (from Nov 2008 to April 2009) of MDF. Results indicate that the Landsat -8 derived NDVI trajectories describing the phenology variation of MDF during spring, monsoon autumn and winter seasons agreed closely with Landsat-7 and MODIS derived phenology transition from Nov 2008 to April 2009. Furthermore, statistical analysis showed statistically significant correlations (p < 0.05) amongst the environmental variables and the NDVIchange between full greenness and maximum frequency stage of Onset of Greenness (OG) activity.. The major change in NDVI was observed in medium (600 to 650 m) and maximum (650 to 750 m) elevation areas. The change in LST showed also to be highly influential. The results of this study can be used for large scale monitoring of difficult-to-reach mountainous forests, with additional implications in biodiversity assessment. By means of a sufficient amount of available cloud-free imagery, detailed phenological trends across mountainous forests could be explained.


INTRODUCTION
In humid subtropical ecosystem of Western Himalayan region of India, forest phenology patterns are potentially governed by climatic conditions (Joshi et al., 2001).Due to this, spatiotemporal patterns in phenological changes are comprehensively considered as important indicators of changes in forest biodiversity (Cleland et al., 2007).Furthermore, Western Himalayan region of India is mainly dominated by Moist Deciduous Forest (MDF) (Jeganathan et al. 2010a, b) and faces four different major seasons (summer, monsoon, autumn and winter) throughout the year.Thus, different growing season dynamics are also strongly connected with environmental variables such as topographic variations, Land Surface Temperature (LST) and rainfall.Due to this fact, detailed and precise description of phenological behaviour of MDF is essential to maintain the forest biodiversity.
The most widely used remote sensing satellite for the study of phenological pattern are Moderate Resolution Imaging Spectroradiometer (MODIS) (Roy et al., 2002) and Advanced Very High Resolution Radiometer (AVHRR) (Townshend and Tucker, 1981) which offer imagery at varying spatial resolution between 250 m and 8 Km with nearly daily temporal resolution.Due to coarser spatial resolution, these datasets are suitable only for large scale analysis of vegetation dynamics (Melaas et al., 2013;Ju and Masek, 2016).Various studies during the past two decades exploited time series NDVI data to track seasonal changes in plant, crop and forest phenology (Moulin et al., 1997;Jonsson & Eklundh, 2002;de Beurs & Henebry, 2010).However, these studies are less reliable at local scale analysis and in areas with mixed land cover where remote sensing sensors with coarser spatial resolution provide surface reflectance with mixture of land cover types, plant types and plant species (Liang et al., 2011).In contrast to this, Landsat imagery offers a possibility for a medium spatial resolution, long term datasets for local scale analysis of forest phenology.Landsat archive has been continuously providing data at medium spatial resolution of 30 m since 1982.it has given impetus to phenology based studies (Woodcock et al., 2008;Li et. al., 2013;Ke at. al., 2015).Moreover, vegetation index differencing is widely used method to study changes in vegetation phenology.This method gives prominence to differences in spectral response of different features (Townshend and Justice, 1995;Fan et al., 2015).
In this paper, we have assessed the phenological changes across MDF in western Himalayan region of India using Landsat 8 Operational Land Imager OLI-derived NDVI and Thermal Infrared Sensor (TIRS)-derived LST.The objective of this study was to examine the effect of environmental variables including elevation, aspect and LST on phenological changes of MDF by regression analysis during post and pre-monsoon seasons from Nov-2013 to Apr-2015.In short, this paper combined vegetation index differencing method to detect phenological changes between mature leaf (Nov) and leaf flush (Apr) duration and then examined the effect of topographic elements and change in LST on phenological behaviour of MDF.

Study area
The study area is located in the Western Himalayan region of Doon valley, Uttarakhand, India (Latitude 29 0 55' to 30 0 30' N and Longitude 77 0 35' to 78 0 24'E) (Figure 1).The Doon valley is surrounded by hills and has varied range of MDF, mainly dominated by Shorea robusta (Sal tree) along with Syzygium spp., Terminalia spp., Ehretia spp., and Litsea spp as codominant species.The average maximum and minimum temperatures of Doon valley are 27.65 0 c and 13.8 0 c, respectively, with average maxima in June (40 0 c) and average minima in January (1.80 0 c) (Mandal and Joshi, 2015).This region faces four different seasons annually: Summer seasons is from March to May and monsoon season starts from June onwards and lasts up to end of September.An average annual rainfall of 2025.43 mm has been recorded, with most of the annual rainfall occurred in monsoon season.Furthermore, post monsoon season covers October and November, while winter season duration is from December to February.Cloud free scenes of Landsat-8 were selected approximately for similar acquisition dates of Landsat-7 scenes (table 1).The LST was calculated from Landsat-8 TIRS bands using ATCOR 3. In addition, topographic variables such as elevation and aspect were generated using Cartosat-I version 3.1 Digital Elevation Model (DEM) (Ritter, 1987) provided by Bhuvan portal of Indian Space research Organization (ISRO).The spatial resolution of DEM was 30m and it was orthorectified and preprocessed by the data provider.

Establish relationship between multi-temporal NDVI and phenology activity of MDF:
The figure 2 summarizes the methodology of this study.Mean, standard deviation, minimum and maximum values of NDVI were calculated for 6 scenes of Landsat-7 and similarly for both time lines of the Landsat-8 datasets.Furthermore, mean values of NDVI trajectories were plotted from Nov to Apr for Landsat-7 (figure 3) and Landsat-8 (figure 4).Two time lines for Landsat-8 data (figure 4) were selected to establish the relationship between phenological activity of MDF and mean values of NDVI.Since Landsat-7 data was used as a reference purpose, comparison was done between Landsat-7 and Landsat-8-based NDVI trajectories.( Where Y represents the response variable, is the intercept, , , , are regression coefficient of the respective explanatory, and ε is random error term.The coefficients were estimated by OLS method.Each global model was screened by an exhaustive model selection procedure via repeated model fitting of all numerous possible sub-models as described by Barton (2015).The model selection procedure followed Latifi et al. (2016).Second-order Akaike Information Criterion (AICc) was calculated for all NDVI change between Nov (A) to Apr (F) phenology stages to identify the sub-model with highest AICc.AICc takes the general form of AIC=2k-2 ln (L), where k is the number of parameters in the model and L denotes the maximization of a likelihood function.The AICc is an extension of AIC in the presence of small sample sizes or when number of fitted parameters is moderate to large (Hurvich and Tsai, 1989;Johnson and Omland, 2004).'MuMin' package in R was applied to calculate the subsets, AICc and select submodels.

Phenology variation of MDF
This section describes the results for NDVI trajectories along with phenology stages of MDF for Landsat-7 and Landsat-8 data separately.Further sections will cover statistical analysis between NDVI change and environmental variables.The curve (figure 3) for Landat-7 derived NDVI closely agreed with MODIS-derived NDVI curve previously published by Upadhyay et al., (2013) for Maximum NDVI value achieved in November 2008 was 0.861 as well as 0.714 for April 2009.These values were closely in agreement with finding of Upadhyay et al., (2013).Since monsoon season ends in September, MDF shows full greenness (mature leaf) in the month of November which leads to higher NDVI values.In contrast to this, April belongs to pre-monsoon season and is the starting point of leaf flush activity which is defined as Onset of Greenness (OG) activity, in which minimum NDVI value was obtained (table 1).However, higher (39.7 mm) rainfall was observed for mid-February 2009 as compared to December 2008 (0.0 mm) and January 2009 (5.2 mm) (table 1), due to which there is a sudden increase in NDVI value during mid-February 2009.Table 1 shows that Doon valley occasionally faces winter rainfall.Therefore similar rise in NDVI values was also observed by Upadhyay et al., (2013).However, they used coarser spatial resolution datasets to understand the phenology of MDF.In contrary, our study used medium spatial resolution dataset (30 m) and additionally estimated the effect of environmental variables causing change in phenological activities of MDF using a regression approach.
Based upon the available cloud free data, similar type of Landsat-8 derived NDVI trajectories were drawn from Nov 2013 to Apr 2014 as well as between Nov 2014 and Apr 2015.
Here, we only show the results from Nov 2013 to Apr 2014, since the final selected regression model included the same variables as the other studied time span (table 2) as well as the rainfall data was only available up to year 2014.Results indicated that trend in the variation of NDVI values is approximately similar for both timelines except for some dates (B and C in figure 4), in which cloud free scenes for similar dates were unavailable.
Table 1: Selected Landsat-7 and Landsat-8 datasets along with rainfall and phenology activity of corresponding months

Statistics between Nov 2013 and Apr 2014
Correlation between the response variable and the environmental variables was shown by Pearson correlation coefficient.A linear and statistically significant correlation was observed between NDVI change and environmental variables.Elevation and ∆TM were significantly correlated with NDVI change (p < 0.05).Moreover, elevation and ∆TM also showed to be positively correlated, whereas the aspect was negatively correlated with NDVI change .In addition, interaction effect was observed between elevation and ∆TM in final selected model (table 2), which suggests their relevance as compared to the interaction between aspect and ∆TM in linear regression model (equation 2).The figure 5 (a) also represents the interaction among elevation and ∆TM for NDVI change .Among the ∆TM range, higher association between elevation and NDVI change was observed within medium and minimum ∆TM categories.

Date
However, the standard error of estimation (i.e.confidence bands) and slope of the regression plot in minimum ∆TM category were affected by small number of observations.In case of univariate relationship between NDVI change and aspect, standard error was substantially larger in northeast direction.It was also observed that maximum variations in NDVI occurred in those areas which are at medium (600 to 650 m) and maximum (650 to 750 m) elevation ranges.

CONCLUSION
This study used a vegetation index differencing method to develop relationships between changing phenology of MDF and a set of commonly-studied environmental variables.We have found that elevation and ∆TM are significantly correlated with NDVI change (p<0.05).The Landsat 8 OLI-based phenological trajectories showed unique phenological characteristics of MDF.Two specific phenological phases of full greenness (mature leaf) and maximum frequency stage of OG activity (leaf flush) of MDF were accurately identified from the temporal profiles of Landsat-7 and Landsat-8 derived NDVI.
Phenology is a major characteristic to be explored in order to sustain forest biodiversity within mountainous forest regions.
Though our results are confined to a reserved forest range, the adopted approach may be applied to other forest areas as well for comparative analyses.Furthermore, the outcomes of this model-assisted analysis provide a cost effective solution to monitor ecosystem dynamics in tedious mountain forest regions.
Coarser spatial resolution data derived from MODIS, AVHRR, and SPOT sensors were utilized for high temporal observations at larger spatial scales.However, significant variations in phenology occur at medium and high spatial resolution which may not be addressed by such sensors.Thus Landsat-7 and Landsat-8 (OLI and TIRS) datasets are potentially appropriate to address fine scale changes in phenology.The existence of other sources of medium resolution optical data such as those from Sentinel 2 warrants a global coverage and temporal continuity of such data for studying vegetation phenology in remote mountainous areas.

Figure 1 :
Figure 1: Location of study area along with sample points (N=100) within Thano forest site.2.2 Methodology 2.2.1 Acquisition and preprocessing of remote sensing datastes : Since major phenological changes for MDF takes place from November to April (Jeganathan et al., 2010b) and a recent study of Upadhyay et al., (2013) has also shown the phenological stages of MDF for Doon valley from November 2008 to April 2009.Thus, for this study we have acquired nearby dates Landsat-7 derived surface reflectance NDVI from November 2008 to April 2009.Landsat-7 data was used for reference purpose only.Temporal time series of Landsat-8 derived surface reflectance NDVI were acquired from United States Geological Survey (USGS).We acquired 7 scenes of processed NDVI from Nov-2013 to Apr-2014 and 5 scenes from Nov-2014 to Apr-2015.Cloud free scenes of Landsat-8 were selected approximately for similar acquisition dates of Landsat-7 scenes (table1).The LST was calculated from Landsat-8 TIRS bands using ATCOR 3. In addition, topographic variables such as elevation and aspect were generated using Cartosat-I version 3.1 Digital Elevation Model (DEM)(Ritter, 1987) provided by Bhuvan portal of Indian Space research Organization (ISRO).The spatial resolution of DEM was 30m and it was orthorectified and preprocessed by the data provider.
Figure 2: Flow chart of methodology 2.2.3 Modelling relationship between NDVI change and environmental variables: The effect of environmental variables (tpography and LST) on changes in phenology between Nov and Apr was modelled as a function of the spectral feature (NDVI change ) derived from vegetation index differencing method (equation 2).A global Ordinary Least Squares (OLS) model was initially set by taking NDVI change as a dependent variable.Independent variables included elevation, aspect and change in LST (∆TM).

Figure 3 :
Figure 3: Landsat-7 derived NDVI trajectories along with phenological activities from Nov 2008 to Apr 2009 OG: Onset of Greenness (leaf flush durtation) and ES: End of Senescence (leaf fall duration)

Figure 4 :
Figure 4: Landsat-8 derived NDVI trajectories along with phenological activities for Nov 2013 to Apr 2014 and Nov 2014 to Apr 2015 durations.
Figure 5: Graphical representation of outputs between Nov 2013 to Apr 2014.(a) ∆TM map representing variation in Max_∆TM (maximum change in LST), Med_∆TM (medium change in LST) and Min_∆TM (minimum change in LST) categories (b) Elevation map representing variation in topography (c) Cross sectional plots of the NDVI change model with an intraction between elevation and ∆TM with emphasis on elevation (d) Univariate relationship of NDVI change with elevation, aspect and ∆TM.Furthermore, estimated univariate association of NDVI change with elevation, aspect and ∆TM are shown in figure 5 (d).The figure5(a) also represents the interaction among elevation and ∆TM for NDVI change .Among the ∆TM range, higher association between elevation and NDVI change was observed within medium and minimum ∆TM categories.