An inflection point-based method for estimating metrics of mangrove phenology combining climatic factors and Landsat NDVI time series

Mangroves are the most productive coastal habitat with long-term storage of abundant biomass and organic carbon (Bernardino et al. 2020). The carbon sequestration potential of mangroves is significantly influenced by their vegetation phenology. The cyclical biological processes of leaf initiation, growth, and shedding regulate the occurrence of photosynthesis, thereby controlling carbon uptake by mangroves. Though the mangroves are evergreen and retain their leaves throughout the year, the peak of the growth period is unique to species and also dependent upon their regional environment conditions. Climate change seriously threatens mangroves through temperature increase, sea-level rise, and changes in rainfall patterns (Ward et al. 2016). The study of phenological traits informs the growth and development characteristics of mangrove species, which helps to understand how mangroves will likely be affected by changing climatic factors (Pastor-Guzman et al. 2018). The search for a sustainable approach to estimate the dynamics of mangrove phenology converges towards the advances in remote sensing technology. The appropriateness of remote sensing techniques becomes even more evident in the domain of mangrove phenology due to the logistical and practical challenges that arise while conducting field surveys in marshy habitats (Wang et al. 2019).

The recent advances in optical remote sensing and the free data policy implemented in 2008 created a surge in research opportunities in environmental remote sensing using the Landsat data archive. The recent Landsat-8 and 9 missions, a collaborative program by NASA (National Aeronautics and Space Administration) and the U.S. Geological Survey (USGS), enabled remote sensing scientists to understand the challenges and prospects of using satellite data in diverse domains of the environment (Wulder et al. 2012). In the context of phenology-related research, optical satellite missions like Landsat (Nguyen et al. 2020), Moderate Resolution Imaging Spectroradiometer (MODIS) (Gong et al. 2015) and Sentinel-2 (Sun et al. 2022) showed their potential and limitations through various research outcomes.

The research on vegetative phenology using satellite data focuses on NDVI (Normalized Difference Vegetation Index), due to its capacity to track seasonal and inter-annual growth variations (Huang et al. 2021). Previous researchers noted that NDVI values exhibited the highest correlation with field-measured mangrove canopy cover compared to any other spectral index values (Tran et al. 2022). The NDVI records the mangrove phenology changes due to subtle variations in climatic factors using the spectral bands in the near-infrared and red wavelength regions (Aji et al. 2023). The challenges related to the lower data availability due to the clouds in the rainy season were tackled by various techniques applied for reconstructing the NDVI time series datasets (Colditz et al. 2008). The capacity of harmonic analysis to handle the constraints of limited optical datasets and the requirement to capture the phenology variations in shorter time windows directed the attention of researchers towards the HANTS (Harmonic Analysis of Time Series) reconstruction method in the recent past (Roerink et al. 2000; Zhou et al. 2021).

In most of the previous research, the inputs for the estimations of mangrove phenology were only the annual time series data of NDVI as a proxy measure of phenology (Ma et al. 2020). Hence the approach could not reveal the influences of climatic factors on the growth patterns of mangroves. The strong connection between the physiological aspects of mangrove growth and climatic factors, especially rainfall and temperature, demands a thorough understanding of their nature. The timing of monsoon onset and the pattern of temperature change control the leaf growth cycle of mangroves, while increased freshwater availability extends their leaf-growing seasons (Dannenberg et al. 2015; Pastor-Guzman et al. 2018). There were attempts to quantify the combined influences of rainfall and temperature through climatic diagrams like Ombrothermic diagrams, which help to derive the season’s quantitative impacts on various biological phenomena (Saleem Khan et al. 2014). However, the combined effects of rainfall and temperature were not analyzed thoroughly in the mangrove phenology-related research.

Most research contributed towards estimating vegetative phenology focused on a fixed percentage-based thresholding approach of phenological descriptors like NDVI and paused several challenges regarding its reliability (Jönsson & Eklundh 2004; Dash et al. 2010). The omission of climatic influences for the estimation of phenological metrics and the fixed thresholding approach were the two major gaps that led to the formulation of the present problem. The framed research questions addressed the scope of quantifying the combined effect of rainfall and temperature as a statistical descriptor. The study also attempted to include the climatic factors together with NDVI for the phenological metrics estimations. The further part of the study focused on developing an algorithm for estimating phenological metrics using the curves generated by the quantitative descriptors of climatic factors and mangrove growth. The present study utilized the facility of the cloud-enabled Google Earth Engine platform (GEE, https://earthengine.google.com) for generating the time-series datasets of NDVI using cloud-free images of Landsat-8 of Pichavaram regions (Li et al. 2019).

The study has two specific objectives, and the first objective was the application of the HANTS reconstruction algorithm to generate the time-series data sets of NDVI using Landsat-8 data of the Pichavaram region from 2013 to 2020. The study’s second objective focuses on devising a method of phenology estimation combining the statistical descriptor of climate and NDVI using the common inflection points and verifying the outcomes with MODIS GPP values of Pichavaram mangroves.

As per the Köppen climate classification, Pichavaram is classified as a sub-humid region distinguished by its hot summers (where temperatures exceed ⁠). Rainfall in the region is primarily caused by the northeast monsoon season, leading to an average yearly rainfall total of 1,310 mm. The climate patterns in the area delineate four distinct seasons: post-monsoon (January–March), summer (April–June), pre-monsoon (July–September), and monsoon (October–December). The present study adopts this categorization of the seasons for executing the objectives (Selvam 2003).

The Landsat-8 satellite images with a spatial resolution of 30 m over the Pichavaram region, accessed through the GEE platform, were the prime datasets utilized for the study. The selection of satellite data acquisition dates (Table 1) spans from 2013 to 2020. The research also guaranteed the inclusion of a minimum of nine cloud-free images for reconstructing NDVI values per year, encompassing at least one image from each season, except for the monsoon period of 2014. The NDVI images for all the chosen dates for HANTS analysis were generated using red (Band-4) and near-infrared (Band-5) bands of Landsat-8. The hourly ERA5 reanalysis data of spatial resolution of 30 km was the source of rainfall and temperature datasets for the selected span of the study. The MODIS operational GPP product MOD17A2H of 500 m resolution was the data source for verifying the performance of the approach developed in the study. The primary objective of this study was to estimate phenological metrics, which inherently have a temporal nature and hence mitigate the limitations associated with coarse spatial resolution temperature and rainfall data from ERA5, as well as the MODIS GPP product.

In this study, we considered the phenological year as the period between the start of the previous year’s pre-monsoon season (1 July) and the end of the following year’s summer season (30 June). For each such phenological year, the four phenological metrics were estimated: Start of Season (SoS), Peak of Season (PoS), End of Season (EoS), and Length of Season (LoS). The present study develops a new algorithm (Algorithm 1) for obtaining phenological metrics of SoS, PoS and EoS by utilizing the inputs of the 8-day average values of the z-score sum of temperature and rainfall and the 8-day average of the reconstructed NDVI values. Initially, HANTS was employed to reconstruct the daily NDVI values for the respective phenological year. Subsequently, these reconstructed daily NDVI values were utilized to compute the 8-day average NDVI values. The simplicity and reduced computational demands of the HANTS algorithm make it an effective tool for reconstructing NDVI time series. It excels especially when handling shorter gaps in data. By utilizing Fourier analysis, HANTS breaks down NDVI datasets into sinusoidal elements, filtering out anomalies and addressing data gaps or noise. Furthermore, the resulting reconstructed time-series data, with its amplitude and phase details, provides valuable insights into phenological characteristics such as annual growth cycle counts and other factors crucial for understanding mangrove growth. The reason for considering the ‘8-day average’ for the above data is its synchronization in duration with the verification product of MODIS 8-day GPP. This could reduce the error induced by daily datasets. The inputs of the algorithm are the plots of reconstructed NDVI (8-day average) and the z-score sum of temperature and rainfall (8-day average) for the particular phenological year. The algorithm necessitates computing both first-order and second-order derivatives of the inputs. The first-order derivative represents the rate of change of a function at a specific point. In this algorithm, we calculated the first-order derivatives of both the z-score sum and NDVI values to identify regions where the rate of change aligns with the indicators of phenological metrics. However, a second derivative reveals the curvature of a function at a given point. Specifically, it indicates whether the function is convex upwards or downwards and conveys the rate of change at which the functional change occurs. Estimating the first-order derivatives of z-score sum and NDVI identifies the potential regions of occurrence of phenological metrics by applying the necessary conditions for each dataset. Further, second-order derivatives facilitate the narrowing down of the regions of metrics to a single day based on sufficient conditions. The algorithm selects the earliest day that meets the criteria in the case of multiple occurrences.

Algorithm 1: Algorithm for mangrove phenology estimation

The methodology adopted ensures the convergence towards the global maxima or minima and estimates phenological events with a maximum error of 8 days. The outcome of the developed algorithm was verified with an independent data source like the MODIS GPP product. The present algorithm eliminates the error ambiguities due to a single data source, such as NDVI with fixed percentage thresholding criteria, which was the approach followed in past studies (Prihantono et al. 2022). On the other hand, Algorithm 1 is a method which recovers daily NDVI values by HANTS first and then computes 8-day average NDVI values so as to improve the data flow. This averaging also follows the MODIS 8-day GPP verification product that helps in enhancing synchronization among data sets and reducing any possible errors that might appear when using daily data sets. Using 8-day moving sums of z-scores for temperature and rainfall and NDVI, the algorithm’s first-order derivatives were computed from which phenological characteristics may exist from necessary conditions. Further, given second-order derivatives, the regions are refined to single-day intervals, thus making certain the convergence to global maximum or minimum places with a maximum error being 8 days. It is noteworthy that, the algorithm takes only the earliest date among the dates consistent with the conditions as a reference when multiple occurrences arise, thereby contributing to the continuous and correct estimation of phenological events. This methodology is also being used to improve on the past schemes that were purely NDVI based with fixed percentage thresholding criteria. However, further vetting is being done to address errors using independent data such as the MODIS GPP product which also offers another reality check.

The present section includes the details of the performance of the HANTS algorithm for the reconstruction of the NDVI time series. Further, the details of climatic quantifiers for understanding the combined effect of rainfall and temperature were presented through Ombrothermic diagrams and statistical quantifiers. The end of the section includes the performance of the phenological metrics algorithm and analysis of the temporal behavior of the derived metrics.

The details of the control parameter selected for HANTS are listed in Table 2.

Table 2

Details of control parameters selected for HANTS implementation

Parameters .	Details .	Value .
Base Period	The period for which the signal to be reconstructed	365 days
Number of frequencies	Total number of harmonics (excluding zero frequency) considered for reconstructing the NDVI data	4
Hi or Lo suppression flag	Flag to indicate the direction of outliers. The low NDVI were eliminated in the current study	Low
Invalid data rejection threshold	The NDVI values below a certain threshold are considered as invalid
Degree of over-determinedness	The minimum number of additional data points provided	2
Fit error tolerance	Permissible fit error between actual values and reconstructed values	0.05

Parameters .	Details .	Value .
Base Period	The period for which the signal to be reconstructed	365 days
Number of frequencies	Total number of harmonics (excluding zero frequency) considered for reconstructing the NDVI data	4
Hi or Lo suppression flag	Flag to indicate the direction of outliers. The low NDVI were eliminated in the current study	Low
Invalid data rejection threshold	The NDVI values below a certain threshold are considered as invalid
Degree of over-determinedness	The minimum number of additional data points provided	2
Fit error tolerance	Permissible fit error between actual values and reconstructed values	0.05

View Large

The study identified the optimal number of frequencies equal to four for HANTS execution and found it appropriate for natural vegetation growth cycle-related studies like mangroves as per previous research (Dash et al. 2010). The HANTS reconstructed 8-day average NDVI values from 2013 to 2020 were considered for relating to climatic factors and phenological metrics estimations.

The analysis of the relationship between the area under the curve of NDVI and climatic factors (Figure 4(d)–4(f)) revealed that the variation in the percentage of wet days (outcome of combined analysis of rainfall and temperature using Ombrothermic diagram) is more related to the phenological variations compared to the individual relationships of total rainfall and mean temperature. The analysis of temporal behavior of the mangrove growth pattern identifies 2017 as the phenological year of the least growth period with less rainfall and, subsequently lower percentage of wet days (Figure 4(b,c)). The phenological year 2017 (July 01, 2016, to June 30, 2017), which comprises the 2016 monsoon period, witnessed a drought throughout south India, and hence the pattern analyzed was in accordance with the actual situations. Also, the phenological year 2016 (1 July 2015 to 30 June 2016), comprising the 2015 monsoon period with maximum rainfall, exhibited peak values of mangrove growth (area under curve (NDVI)) and the percentage of wet days. The relationship between the percentage of wet days and the area under the NDVI curve emphasizes the importance of the combined effect of temperature and rainfall. It hence necessitates further investigations in a shorter temporal window (daily) than the Ombrothermic diagram (monthly) to estimate phenological metrics.

The present study utilized Landsat-8-derived NDVI as a proxy parameter of mangrove phenology and attempted to evaluate the simultaneous influences of rainfall and temperature on the dynamics of NDVI. The HANTS reconstruction algorithm with a four-harmonic frequency was found to be appropriate for the seamless data generation of NDVI throughout the selected period of 2013–2020. The Ombrothermic diagram-derived percentage of wet days revealed the significant correlation (⁠⁠) with mangrove phenology, which was insignificant with the single climatic factor. As the Ombrothermic diagram was for monthly datasets, another quantifier called the sum of the z-score based on statistical aspects of daily climatic datasets were estimated. The new approach for estimating the metrics of phenology considering the combined effect of rainfall and temperature together with NDVI was found to be well mapped with the actual scenario. The synchronization of our results with MODIS GPP values verified the reliability of the new method. The combined utilization of ‘8-day average NDVI’ and ‘8-day average z-score sum’ reduce the chance of error due to single datasets. The examination of mangrove phenology metrics from 2014 to 2020 showed a consistent pattern in the timing of SoS, PoS, and EoS across various seasons throughout these years with an average LoS of 230 days. It also highlighted a clear correlation between monsoon rainfall patterns and mangrove phenology, particularly evident in the relationship between PoS and EoS. The results revealed the influences of drought in 2016 through the metrics of 2017 with minimum LoS and least GPP values compared to other years. The PoS recorded the optimal temperature of ⁠, which is required for mangrove growth, and the descriptors of phenology extracted through the developed methodology were in tune with the climatic pattern. Hence the combined approach of integrating climatic factors and NDVI proved to be effective for mangrove phenology estimations for other regions. The study’s lacunae included the utilization of only 7 years of satellite datasets for evaluating phenological metrics patterns, as well as the coarse resolution of MODIS GPP used for verifying occurrences of phenological metrics. These factors impeded the present study’s ability to accurately estimate the precise delay period between phenological events and their correlation with climate factors. The availability of a longer time-series dataset could enhance the assessment of phenological trends and offer deeper insights into climate change impacts. However, the present study focuses predominantly on the robustness of the developed methodology for analyzing mangrove phenology, with minimal limitations arising from the absence of an extended time series dataset in this specific context. The future prospectus includes the replacement of NDVI with GPP values for a longer period which may reveal finer details of mangrove phenology. Phenology determined by GPP offers a better grasp of shifts in carbon sequestration, surface energy, and water balances compared to phenology estimated based on NDVI values. Moreover, integrating flux tower information could provide a more accurate approach to verifying the current results. The mangrove species-wise phenology metrics estimation using finer resolution remote sensing images may enhance the knowledge of species’ response to various climatic factors and their favorable conditions of growth.

The authors express their gratitude to the European Union’s Copernicus Climate Change Service for granting them access to ERA5 datasets. Additionally, they would like to extend their appreciation to the team behind the USGS/NASA Landsat Program for their efforts in the Landsat mission and for providing free data access.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.