Tracking Long-term Phenological Shift in Response to Climatic Parameters in Chitwan National Park, Nepal

normalized difference vegetation index (NDVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) was used to study the phenological shift in Chitwan National Park of Nepal, which is home to unique biological resources, in response to two major climatic drivers: temperature and precipitation. The four transition stages of greenness onset, maturity onset, senescence onset, and growing period were determined by fitting spatially averaged NDVI values using the phenofit package of R. It was found that the greenness and maturity onsets have been delayed over the years while the growing period has seen fluctuations due to variations in senescence onset. Precipitation was correlated positively with NDVI while temperature was negatively correlated with NDVI. Moreover, the rainfall one month earlier better explained the NDVI variability than the amount of rainfall in the same month because of the stronger correlation. Overall, this study indicates that climate variability is affecting the phenology of vegetation, and the results can help in performing suitable checks and assessments of the ecosystem in Nepal.


Introduction
The study of plant and animal life cycles in relation to the seasons is known as phenology. It is the study of the occurrence of life cycle events at the population level, with a particular emphasis on how they have been altered in response to climate change. Long-term records are often used, and activities such as greening, flowering, hatching, and leaf fall are all included. (1) Remote sensing (RS) phenology can provide information on large-scale phenological trends that would be extremely arduous to detect from the ground. Satellites allow regular monitoring of the global land surface, and the thus obtained remotely sensed information can be used in assessing phenological data to appraise critical patterns such as crop conditions, drought severity, and wildfire risk, as well as tracking invasive species, infectious diseases, and insect pests. (2) As the present world is greatly affected by the threats of global warming, climate change, and deforestation, it has become crucial to monitor environmental changes, and the study of phenology is a primary tool for keeping records of life cycle trends and the impacts of climate change on ecosystems. There has been evidence that over the last few decades, the timing of seasonal activities of many plant and animal species has altered and that these shifts are influenced by climate change. (3)(4)(5) However, in Nepal, this critical issue of phenological change has not been considered as an important research topic, despite being a country that has much of its land covered by forest and hosts critically endangered species. The studies previously performed in Nepal focused on studying either growing trends using ground data or normalized difference vegetation index (NDVI) patterns to describe climatic fluctuations. Our study provides a novel approach in integrating information derived from a remote sensor to study plant phenological patterns and analyzing the changes while incorporating weather data. The study is also focused on an area that is relatively little affected by human intervention, enabling the changes in nature to be attributed to natural phenomena.
Plant phenology has been proposed as an indicator of climatic differences and global changes by the European Environmental Agency and the Intergovernmental Panel on Climate Change (IPCC). (6) Field-based ecological studies have demonstrated that vegetation phenology tends to follow relatively well-defined temporal patterns. Some vegetation types exhibit multiple modes of growth and senescence within a single annual cycle. Therefore, RS-based methods need to be sufficiently flexible to allow for this type of variability. (7) Thus, phenological events are sensitive to environmental parameters such as temperature, rainfall, and pressure. (8) Schmid et al. (9) used Google Earth Engine (GEE) in conjunction with Landsat 5 and 8 images to study the development of NDVI in three study areas across Germany, the Biosphere Reserve Schwäbische Alb, National Park Hainich, and Biosphere Reserve Schorfheide-Chorin, from 1984 to 2016 and found that the decline in NDVI was linked to large-and small-scale changes in land use and associated with declining biodiversity. Workie and Debella (10) used time-series NDVI generated from 8-day Moderate Resolution Imaging Spectroradiometer (MODIS) data. The aggregated NDVI was used along with temperature and rainfall data for 12 ecoregions in Ethiopia in the GEE environment, which was followed by Fourier smoothing to eliminate various noises and fitting with sigmoid vegetation growth functions to identify phenophases. (10) The study, carried out over 14 years, found that the start of the growing period became earlier, and the growing period increased in length for most of the ecoregions in Ethiopia. Potter (11) analyzed the change in Landsat NDVI between 1985 and 2017 within the Santa Rosa Mountains Wilderness in the southern deserts of California at four elevation zones between 500 and 2500 m, which showed that green cover dropped notably in the below-average precipitation years of 2002,2007,2014,2015, and 2016, whereas it increased sharply in the above-average precipitation years of 1998, 2005, 2010, and 2017. This was supported by the work of Kelly and Goulden. (12) Pan et al. (13) studied the impacts of climate change on the wetland vegetation of Dunhuang Yangguan National Nature Reserve in northwest China by comparison and analysis of satellitederived NDVI and climate change factors (temperature, precipitation, and snow depth). The findings revealed that, over the last 30 years, global warming has accelerated while precipitation has remained relatively unchanged. For the arid region, snowmelt contributed more to the increase in wetland vegetation in the region than the increase in precipitation and temperature. Using MODIS nadir bidirectional reflectance distribution function (BRDF)-adjusted reflectance (NBAR), Zhang et al. (14) showed that urban areas had the earliest green-up and the latest dormancy, while croplands had the opposite pattern. Wolf et al. (15) highlighted the importance of other abiotic variables e.g., light infiltration and nutrient concentrations, which highly influenced the diversity of plants in an area. Many plant species bloomed earlier in response to decreased biodiversity, and the species with the largest ranges during peak flowering times had the most robust responses to changing biodiversity.
RS-based phenological studies involve two approaches: extracting trends in vegetation indices (VIs) or computing the phenophase time, which can then be applied to identify the patterns over space and time. (10) Smoothing is involved in extracting the trend in order to overcome noises, which can be done using i) statistical, ii) curve fitting, and iii) data transformation techniques. (16) Data can be smoothed by different methods such as the harmonic Fourier transform, (17) piecewise logistic function, (14) and polynomial curve fitting. (18) Jones et al. (19) analyzed global phenology cycles over a six-year record using satellite passive microwave RS-based vegetation optical depth retrievals derived from daily time-series brightness temperature measurements. Satellite-derived VIs and phenological metrics have been found to be comparable with the results of near-surface RS and thus enable large-scale comparison and analysis. (20) The derived indices have also been found to be in agreement with in situ data. (21) The phenology of Nepal and its drivers have been rarely studied. Although the country hosts a wide diversity of vegetation including rare species, the phenology of the vegetation has hardly been discussed and researched. Chitwan National Park (CNP) is known throughout the world for its rare flora and fauna as well as its outstanding natural features. The park aims to protect many endangered wildlife species, including the Royal Bengal tiger and the world's second-largest viable population of greater one-horned rhinoceros. Multiple forms of woodland, wetland, and grassland habitats can be found in the park's core region and buffer zone, resulting in rich habitat diversity. The park protects the delicate Churiya Hill ecosystem in the south and the lowland ecosystems of the inner Terai valley. The World Heritage Convention declared CNP as a World Heritage Site in 1984, acknowledging its unique biological resources (UNESCO/IUCN-2003). (22) Additionally, understanding how vegetation responds to climate change will help improve the management of natural resources and the development of efficient climate change adaptation strategies.
Working on large sets of satellite images in edge devices is labor-intensive, tedious, and often requires a lot of storage. Using several computers to manage a huge volume of Earth observation data over large spatial and temporal scales is one possibility, but it is also an expensive investment for a few times of usage. A cloud-based public computing platform will be the ideal solution for all of the problems mentioned above. This will save time, money, and resources by addressing the problems of infrastructure, expense, and computing time all at once. GEE is a cloud-based computing platform with a large volume of satellite data. This approach gives users free and easy access to data and makes it possible to analyze large amounts of data in seconds. (23) This study is a novel work in Nepal, which tracks the long-term phenological shift in CNP, a protected region hosting diverse ecosystems, by coupling remotely sensed datasets. In this paper, we present results obtained by deriving NDVI from satellite RS-based observations and analyzing how the peak and changes in NDVI have altered in response to changing climatic parameters (temperature and precipitation) over the years.

Study area
With a total area of 952.63 sq. km, CNP is located between 27° 34' and 27° 68' N and 83° 87' and 84° 74' E, while the buffer zone extends further from 27° 28' to 27° 70' N and 83° 83' to 84° 77' E as shown in Fig. 1. It lies in the subtropical inner Terai lowlands of south-central Nepal and spans portions of four districts, namely, Chitwan, Nawalparasi, Parsa, and Makawanpur. (22,24) Established in 1973, it was granted the status of a World Heritage Site in 1984 for its natural and cultural heritage. The park consists of diverse ecosystems including the Churia Hill forests, oxbow lakes, and the floodplains of the Rapti, Reu, and Narayani rivers. (25) The central location of CNP makes it a perfect case study for this novel work and validates the satellite RS-based phenology and the understanding of the changes in NDVI over time.

Satellite data
The satellite data for this study were acquired from the MODIS aboard the National Aeronautics and Space Administration (NASA) Terra and Aqua satellites. The satellites have a sun-synchronous orbit and a temporal resolution of 1-2 days. The NBAR product provides 500-m-resolution reflectance data of MODIS bands 1-7. These are adjusted using a bidirectional reflectance distribution function to model the values as if they were collected from a nadir view. MODIS NBAR data have the characteristics that atmospheric contamination is reduced and cloud cover is explicitly masked. (8) The seven spectral bands are explicitly designed for land surface monitoring. Table 1 shows the details of the MODIS NBAR bands along with potential applications.

Climate data
Amongst the climate parameters, precipitation and temperature were used to analyze how the phenological trend has changed under the influence of these parameters. Time-series rainfall data were generated using Climate Hazards Group Infrared Precipitation with Station (CHIRPS) data for the study period. CHIRPS data incorporate 0.05° resolution satellite imagery with in situ station data to create gridded rainfall time series for trend analysis and seasonal drought monitoring. (26) Land surface temperatures were derived from the MOD11A1 Version 6 product, which provides daily per-pixel land surface temperature and emissivity with 1 km spatial resolution in a 1200 by 1200 km grid. (27) The pixel-valued temperature was reduced spatially to generate a mean for the study area. Figure 2 shows the monthly average temperature and Fig. 3 shows the average and total precipitation for the CNP region.

Vegetation index
NDVI is a graphical indicator of vegetation greenness, which is mathematically obtained as the ratio of the spectral reflectance difference between the near-infrared (NIR) and red bands to where R NIR and R Red are the spectral reflectances in the NIR and red bands, respectively. (13)   NDVI is a common index used in RS studies of vegetation. NIR bands are more sensitive than other bands to vegetation because chlorophyll reflects more NIR and green light than light with other wavelengths and absorbs more red and blue light. (2) Because all the datasets used are available from the GEE repository, the GEE platform was used for all the calculations and extraction of the data.

Phenology
The study of cyclic biological events is known as phenology. It is the study of the occurrence of life cycle events at the population level, with a particular emphasis on how they react to climate change. Long-term records are often used, and activities such as flowering, leaf fall, hatching, and annual migration are often included. (1) In today's environment, where deforestation is a major problem, climate change is a huge concern, and wildfires have become more frequent, it is critical to monitor the remaining green cover. (28) The wide and frequent coverage of satellite RS-based NDVI makes it the best indicator for studying long-term phenology. The four key transition dates in the annual cycle of vegetation phenology can be inferred from RS. (8,14) 1. Green-up: Indicated by the increase in NDVI, the period at which photosynthetic activity commences. 2. Maturity: Indicated by the peak in NDVI, the period at which the green leaf area is maximum. 3. Senescence: Indicated by the decline in NDVI, the period at which photosynthetic activity and green coverage begin to subside rapidly. 4. Dormancy: Indicated by the trough in NDVI, the date at which corporeal function nearly ceases.

Phenological trends over the years
The phenophases were identified for six distinct periods:  Table 2 shows the vital periods in the plant cycle including the greenness onset, maturity onset, and senescence onset. The table shows how the phenological stages changed over time. Both of the earlier phenophases, the greenness onset and maturity onset, have been delayed while the growing period has fluctuated over the years. Compared with the year 2001, the greenness onset was delayed by 19 days, and maturity onset was delayed by 24 days in 2019. Moreover, the growing period has significantly shortened from 196 to 180 days. Among the years, the greenness onset showed the least fluctuation as it usually started toward the end of April. Although the maturity onset was in June for all periods, senescence was significantly altered, resulting in the varied growing period. Table 3 shows the peak NDVI date and the date of peaking for different periods. The peak value of NDVI has fluctuated very little, but the dates for the peaking ranged from early July to mid-September. The date was much later for the period 2019-2020, and from Fig. 4, we can observe that it also lacked a distinctive phenological pattern, unlike other years. Table 4 shows the correlation between NDVI and the climatic parameters. Precipitation shows a positive correlation with NDVI while temperature shows a negative correlation. This   The precipitation period is a month earlier than the NDVI period.

Effect on NDVI due to climatic parameters
indicates that as the vegetation starts greening up, the temperature decreases and vice versa. Moreover, the average precipitation one month earlier seemed to affect NDVI more than that of the same month because it showed a stronger correlation. Early precipitation would therefore result in the early greening of plants, and extended precipitation would result in an extended growing period. The temperature was weakly correlated with NDVI, so the effect of temperature on NDVI is not as strong as the effect of precipitation, which shows a stronger correlation.

Discussion and Conclusion
The findings presented in this paper were obtained by deriving NDVI from satellite RS observations and examining how it has changed over time in response to the changing climatic factors of temperature and precipitation. We used a cloud-based computing platform (GEE) to acquire data for CNP and applied fitting algorithms using the phenofit package of R. The spatial mean NDVI of the park was derived to determine phenological cycles for six different periods. The changes in NDVI were assessed using the correlation with the drivers of climate change, i.e., temperature and precipitation. Precipitation of the same month showed moderate correlation with NDVI (average 0.416), while the precipitation of the previous month had a correlation coefficient of 0.60, signifying a stronger correlation and thus more effect on NDVI. The temperature was weakly correlated with NDVI, implying that vegetation is less affected by temperature. Greenness onset showed the least variation over the years, usually starting toward the end of April. While maturation usually began in June, senescence was greatly altered, resulting in a varying growth period. Despite Landsat providing better spatial resolution, MODIS data were used due to the better temporal resolution. Note that the ground phenology of plants is not studied by the officials of CNP, which makes it impossible to compare any satellite-derived data with ground truth data. We are hopeful that the insights from this paper will help conservationists become more aware of the need to track these changes.
We compared our work with the research conducted by Zhang et al., (7) who monitored the vegetation phenology in the northeastern United States. Greenness onset was similar, occurring around April in both works, while maturity was reached at different periods after greenness onset. The differences and delays could be explained by differences in the elevation, types of forest cover, and climatic conditions of the study areas. Workie and Debella (10) studied the phenological shift behavior of Ethiopia's different ecoregions for the years 2002 and 2015. For most of Ethiopia's ecoregions, the growing period started earlier and the growing period was longer in 2015. The work was performed for two specific years, making it impossible to conclude whether the changes were continuous or just an annual fluctuation, in contrast to our work, which showed the trends over the years. However, the findings of the climatic influence on phenological shifts were very similar, with both works indicating a strong positive correlation between rainfall and NDVI, and thus the change in phenophases was confirmed to be influenced by the shift in the rainfall trend to earlier in the year Because the temperature was weakly correlated with NDVI, Workie and Debella also explained the behavior by implying that greening up in vegetation is less affected by temperature drops, given the importance of photosynthetic behavior to plants. This may imply that the minimum temperature is sufficient to catalyze compound responses, while the light energy related to the temperature is sufficient to meet the energy prerequisites for photosynthetic action. Rather than actuating NDVI, any increase in temperature promotes evapotranspiration, making vegetation lose its leaf dampness content and gradually become light green or wither. Pan et al. (13) presented the impacts of climate change on the wetland vegetation of Dunhuang Yangguan National Nature Reserve (YNNR) in northwest China by comparing two approaches: trend and correlation analyses and time series analysis. The wetland vegetated area (WVA) and the spatial mean NDVI (mNDVI) of the entire wetland vegetation exhibited an increasing trend. Trend and correlation analyses of the annual maximum values were found to be more useful than time series analysis. The study showed that in dry areas, the most significant requirement of wetland vegetation is water accessibility in soils, which is identified with surface water confinement and the release of groundwater. Thus, their work illustrates the limitation of our work, in which we only focused on the major climatic parameters affecting the phenological cycle. Thus, to determine if the cycle has been altered over the long term, other parameters should also be taken into consideration.
We encountered some problems while conducting the study. The first major problem was obtaining cloud-free data for the monsoon season. Because the data volume was large, it was not possible to use individual datasets. The phenofit package of R provided a solution to the problem by fitting a phenological curve, which ignores extreme values. The curve not only helped to remove these outliers but also assisted in defining a clear transition cycle. Another problem was obtaining climatic data over a longer period. Because station data were not available to the public, it was necessary to obtain the climatic data from the accessible remotely sensed CHIRPS data and land surface emissivity. CNP comprises different ecosystems including wetlands and grasslands and thus has different phenological cycles. Because this study considered the spatial mean NDVI of the whole park, our results do not represent those of an individual ecosystem. Thus, it is recommended that the discrete ecosystems are evaluated, as well as how their phenomena have changed under the influence of climatic parameters. No comparisons were made in this work between ground observations and the RS-based results because of the coarse resolution of MODIS data and the field-based methods being species-specific. Also, different species of plants and their numbers present in an area will affect the cycle. Thus, considering biotic components such as biodiversity would yield additional findings. The phenological patterns are characterized by the species of plants, which in turn are influenced by altitude, weather phenomena, and natural availability. Incorporating altitude variations will yield some valuable insights into which belts of Nepal (plains, hills, or mountains) are most severely affected by the changes in climatic parameters, resulting in a difference in growth patterns.
Finally, our research has implications for ecological studies because the findings provide a wealth of opportunities to understand ecosystems when coupled with species distributions and will help show how CNP has changed over the years. Our approach is also capable of identifying phenologic behavior characterized by multiple growth and senescence periods within a single year. By taking this as an initial step, in further research on CNP, the area of study can be expanded to the hilly and mountainous regions in the east and west. Installation of suitably equipped weather stations and high-resolution data capture will be very effective in understanding specific vegetation responses to the changing climate in the Himalayan country of Nepal.

About the Authors
Aman KC received his B.E. degree in geomatics from Tribhuvan University, Nepal, in 2019. Since December 2020, he has been working as a surveyor in Kanchan Rural Municipality in Nepal. His responsibilities include planning field surveys and assisting in the collection of data with the aid of land cover maps by employing emerging technologies and methods for managing conflicted lands. His research interests are in environmental remote sensing, environmental and land cover changes, terrestrial ecosystems, and phenological studies. (amankc929@gmail.com) Tri Dev Acharya is currently a postdoctoral researcher at the Institute of Transportation Studies (ITS), University of California, Davis. He is involved in the UC Davis Hydrogen Study Group, which is studying future hydrogen system design, scaling-up, and optimization within California. He received his B.E. degree in geomatics from Kathmandu University, Nepal, in 2011 and a combined M.S. and Ph.D. degree from the Department of Civil Engineering, Kangwon National University, Korea, in 2018. He worked as a postdoctoral researcher in Korea and China, focusing on geospatial data preparation, modeling, and simulation of land cover, surface water, disasters, and transportation. (tdacharya@ucdavis.edu) Nimisha Wagle received her B.E. degree in geomatics from Kathmandu University, Nepal, in 2017. She is a survey officer at the Survey Department, Government of Nepal. Her interests are related to geospatial data preparation, machine/deep learning, and mapping for agriculture, land cover, and surface water. (nimisha.wagle@nepal.gov.np)