Temporal Variability of MODIS Phenological Indices in the Temperate Rainforest of Northern Patagonia

Abstract

Western Patagonia harbors unique and sparsely studied terrestrial ecosystems that are threatened by land use changes and exposure to basin-scale climatic variability. We assessed the performance of two satellite vegetation indices derived from MODIS–Terra, EVI (Enhanced Vegetation Index) and NDVI (Normalized Difference Vegetation Index), over the northern and southern sectors of the Chiloé Island System (CIS) to advance our understanding of vegetation dynamics in the region. Then we examined their time-varying relationships with two climatic indices indicative of tropical and extratropical influence, the ENSO (El Niño–Southern Oscillation) and the Antarctic Oscillation (AAO) index, respectively. The 17-year time series showed that only EVI captured the seasonal pattern characteristic of temperate regions, with low (high) phenological activity during Autumn-Winter (Spring–Summer). NDVI saturated during the season of high productivity and failed to capture the seasonal cycle. Temporal patterns in productivity showed a weakened seasonal cycle during the past decade, particularly over the northern sector. We observed a non-stationary association between EVI and both climatic indices. Significant co-variation between EVI and the Niño–Southern Oscillation index in the annual band persisted from 2001 until 2008–2009; annual coherence with AAO prevailed from 2013 onwards and the 2009–2012 period was characterized by coherence between EVI and both climate indices over longer temporal scales. Our results suggest that the influence of large-scale climatic variability on local weather patterns drives phenological responses in the northern and southern regions of the CIS. The imprint of climatic variability on patterns of primary production across the CIS may be underpinned by spatial differences in the anthropogenic modification of this ecosystem, as the northern sector is strongly modified by forestry and agriculture. We highlight the need for field validation of satellite indices around areas of high biomass and high endemism, located in the southern sector of the island, in order to enhance the utility of satellite vegetation indices in the conservation and management of austral ecosystems.

1. Introduction

Anthropogenic activities have widespread and profound effects on organisms, ecosystems, the biosphere and the climate system [1,2,3,4]. Primary production is the basis of terrestrial ecosystems; trees represent 80% of the plant biomass worldwide and close to 50–60% of the annual net primary production [5]. Despite the global extent of human impacts, the effects of climate over regional scales are highly asymmetric and mainly manifested through phenological changes such as the timing of spring flowering and leaf unfolding, or differential response between early-blooming herbaceous and late-blooming woody plants (see [4]). Thus systematic observations of tree-dominated ecosystems are key to understanding long-term effects of human activities and responses of terrestrial primary production to climate variability [6].

South America is home to the largest austral temperate rainforest in the world [7], which includes extensive stands within the CIS that are dominated by Valdivian and North-Patagonian type forests. Located on the western side of northern Patagonia (41–43${}^{\circ }$S), the CIS includes the largest island of the Patagonian archipielago, which harbors numerous endemic species [7]. The CIS is recognized as a biodiversity “hotspot” [8]. However, intense forestry and agricultural activities have degraded terrestrial ecosystems, fragmenting or replacing the native temperate rainforests and peatlands of the island, particularly over the northern sector [9]. Human activities have degraded old–growth forests and driven changes in specific ecosystem services such as recreation and ecotourism potential [10,11]. The climate of the CIS and across Patagonia is controlled by the low-level westerly wind flow of the Southern Hemisphere, which produces spatial variability of synoptic-scale precipitation [12,13]. The maritime supply of moisture and the orographic enhancement of the westerlies upstream of the Andean slopes maintain high year-round annual precipitation (2000–3000 mm; maximum during austral Autumn–Winter) and a stable mean annual temperature around 10 ${}^{\circ }$C. [13]. A drying trend has been observed across Patagonia over the past half-century, partly driven by anthropogenic climate change [14]. The exposure to climatic and anthropogenic effects has motivated some conservation efforts to protect and manage the unique biodiversity of the temperate ecosystems of western Patagonia, which are still sparsely studied [15,16].

Phenology is defined defined as the study of the timing of periodic biological (e.g., seasonal cycles in primary production) and non-biological events as climate change and increasing human activities [4,17]. For example studies have documented the effects of recent climate change on the phenology, e.g, advanced the spring and delayed the arrival of winter with consequences for ecological processes, agriculture, forestry, human health, and the global economy [18]. Recent satellite-based phenological studies have highlighted the relevance of studying the relationship between vegetation indices and productivity in different terrestrial ecosystems [3,19,20]. For this purpose, transformations of multiple spectral components (e.g., vegetation indices) have been developed to study the contribution of vegetation properties using remote sensing. Two of these indices, NDVI (Normalized Difference Vegetation Index) and EVI (Enhanced Vegetation Index) (see Methods for details) provide useful radiometric and biophysical information for land surface characterization using different spectral saturation effects [21]. These satellite-based vegetation indices have allowed the study of phenological variation over large spatial and temporal scales [22]. In this way, satellite-based vegetation indices may provide a powerful starting point to understand basic ecosystem processes and their variability in poorly studied regions of high conservation value such as the CIS. Our study objetives were (i) compares the performance of EVI and NDVI in capturing phenological variability across the complex mosaic of vegetation types in the productive CIS terrestrial ecosystem; and (ii) examines their dynamic relationship with large-scale climate indices.

2. Methods

2.1. Satellite and Climate Data

We used data from MODIS (Moderate Resolution Imaging Spectroradiometer) aboard the Terra (EOS AM) satellite to evaluate phenological patterns between 2000 and 2016, which allowed us to examine patterns of interannual variability. This sensor has 36 spectral bands in Visible (400–600 nm), Near-Infrared (NIR1: 841–875 nm; NIR2: 1230–1250 nm) and Shortwave-Infrared (SWIR1: 1628–1652 nm, SWIR2: 2105–2155 nm) with high temporal (1–2 days) and spatial (250, 500, 1000 m) resolution [23,24]. NDVI and EVI are computed from atmospherically corrected surface reflectance that has been masked for water, clouds, heavy aerosols and cloud shadows [24]. We selected these vegetation indices following their ability and sensitivity to monitor community composition, quantifying the relationship between the MODIS vegetation indices and LAI (Leaf Area Index). Our assessment includes unique ecosystems such as peatlands and bogs, which are relevant players in the global carbon cycle [25], and stands of endemic, broad-leaved evergreen rainforest species (e.g., Eucryphia cordifolia) [26].

To evaluate the spatio-temporal phenological variability in the Chiloé Island System (Figure 1), we used 16-day composites from MODIS (MOD13Q1 version 6) produced by the Earth Observing System (EOS) program. NDVI is the continuity index from the existing National Oceanic and Atmospheric Administration Advanced Very High Resolution Radiometer (NOAA-AVHRR), while the EVI is closer to a functional index because of its better performance in the Near-Infrared (NIR) part of the electromagnetic spectrum and greater sensitivity over regions of high biomass [24].

The Normalized Difference Vegetation Index (Equation (1)) is a normalized ratio of the NIR and red MODIS bands:

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$$

(1)

where ${\rho }_{NIR}$ and ${\rho }_{RED}$ are the surface bidirectional reflectance factors for their respective bands [21].

The three-band ratio EVI (Equation (2)) is defined as:

$$EVI=2.5\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{1+{\rho }_{NIR}+(6\times {\rho }_{RED}-7.5\times {\rho }_{BLUE})}$$

(2)

where $\left(\rho \right)$ are atmospherically corrected surface reflectances in $\left(NIR\right)$, $\left(RED\right)$ and $\left(BLUE\right)$ bands derived from MODIS, and the coefficient 1 accounts for canopy background scattering and the blue and red coefficients, 6 and 7.5, minimize residual aerosol variations [27].

We used two climate indices that influence the delivery of moisture to southern South America [14,28]: ENSO (El Niño–Southern Oscillation) and AAO (Antarctic Oscillation Index). For ENSO we used monthly values of the Niño 3.4 index obtained from the NOAA Working Group on Surface Pressure (www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Data/nino34.long.anom.data). Niño 3.4 values are the average monthly SST anomalies across the Equatorial Pacific, from around the dateline to the coast of Peru [14,29]. We obtained monthly data of AAO from the NOAA Climate Prediction Center (http://www.cpc.noaa.gov). The AAO is calculated by projecting the monthly mean 700-hPa height anomalies onto the leading EOF of monthly mean height anomalies south of 20${}^{\circ }$S for the 1979–2000 period [30]. The monthly climatic indices (ENSO and AAO) were linearly interpolated to match the observational time step of the vegetation indices from MODIS for joint analyses.

2.2. Statistical Analyses

We used a wavelet analysis to assess the statistical significance and quantify the principal rhythmic components in transient temporal dynamics of primary productivity derived from EVI and NDVI time series. Wavelet is a family of functions derived from a “mother wavelet” that provides information on the distribution of variance between frequency and different time locations (see [31]). In this study we used a Continuous Morlet wavelet that has been widely used in the study of ecological signals compared to other decomposition schemes [31,32]. This time–frequency Morlet wavelet is defined as:

$$\psi \left(t\right)={\pi }^{-1/4}exp\left(-i{\omega }_{0}t\right)exp(-t^2/2)$$

(3)

The Morlet wavelet $\psi \left(t\right)$ is the result of a sinusoidal complex $exp(-2i2{\omega }_{0}t)$ by a Gaussian envelope $exp(-t^2/2)$ where ${\omega }_0$ is the central angular frequency of the wavelet. The term ${\pi }^{-1/4}$ is a normalization factor to ensure unit variance [31]. To quantify statistical relationships between two temporal series, wavelet coherence can be computed. This analysis provides local information about where two non-stationary signals ($x\left(t\right)$ and $y\left(t\right)$) are correlated at a particular period (or frequency) [31,33]. The wavelet coherency is defined as the cross-spectrum normalized by the spectrum of each signal:

$$R_{x,y}(a,t)=\frac{<W_{x,y}(a,t)>}{<W_{x,x}{(a,t)}^{1/2}><W_{y,y}{(a,t)}^{1/2}>}$$

(4)

where $<>$ denotes a smoothing operator in time and scale and $W_{x,y}(a,t)$ is the cross-wavelet transform. $R_{x,y}(a,t)$ = 1 when there is a perfect linear relationship at a particular time and frequency between the two signals [33]. To explore the temporal variability in phenological patterns derived from both EVI and NDVI vegetational indices from satellite data, we separated the CIS into two sub–regions: north and south of 42.4${}^{\circ }$S , and averaged all valid pixels across each (see Figure 1).

3. Results

Time series of indices of vegetation showed seasonal and interannual variability over the study period. Across the northern sector NDVI showed a similar temporal trend, with minimum (NDVI${}_{min}$ = 0.33) and maximum (NDVI${}_{max}$ = 0.84) values in June 2001 and May 2016, respectively (Figure 2a; red line). Over the same area, EVI fluctuated between EVI${}_{min}$ = 0.22 in June 2001 and EVI${}_{max}$ = 0.49 during January 2001 (Figure 2b; red line). For the southern region, the minimum (NDVI${}_{min}$ = 0.41) and maximum (NDVI${}_{max}$ = 0.83) NDVI values were observed in June 2001 and April 2013, respectively (Figure 2a; blue line), while the minimum EVI value (EVI${}_{min}$ = 0.25) was observed during July 2002 and EVI${}_{max}$ = 0.54 was found during December 2012 (Figure 2b; blue line). The structure of raw EVI and NDVI time series suggested a similar temporal pattern for both indices, but differences emerged when we examined their annual cycles. The annual cycle of EVI indicated a pattern consistent with seasonal variation, with high values during austral Spring (October–December) and Summer (January–March) and minimal values during austral Autumn (April–June). The annual cycle of NDVI, however, showed values closer to saturation and greater variability throughout the year (Figure 2a).

The local power spectrum of EVI indicated the presence of a significant annual mode in the northern area that decayed in spectral intensity towards the end of the study period (Figure 3a, left panel). The annual component of temporal variability was significant for NDVI only for the 2001–2006 period (Figure 3b, left panel). The global wavelet power spectrum revealed differential patterns in the main rhythmic components between vegetation indices for the northern region. The EVI time series was dominated by the annual component throughout the 17-year study period (Figure 3a, right panel) while NDVI showed two significant modes, annual and bi-annual, which did not persist over time (Figure 3b, right panel). For the southern region, the local power spectrum for EVI showed an intense significant annual band that persisted throughout the study period (Figure 3c, left panel), while for NDVI the significant 1-year periodicity was observed only during the 2001–2006 period (Figure 3d, left panel). For both indices, the annual cycle was the only significant periodic component, as indicated by the global power spectrum (Figure 3c,d right panels). Our results indicated that EVI clearly outperforms NDVI as an indicator of phenological changes across our study region. Hence, we focused our analyses to assess the influence of climatic variability on phenological variation in the northern and southern sector of the CIS using only the EVI time series.

The temporal co-evolution of EVI and climate indices showed complementary patterns across the two study regions of the CIS. For the northern sector, wavelet coherence between EVI and ENSO showed transient significant coherence in the annual cycle between 2003–2006, 2012 and 2015, while a broad period of coherence for the 2–2.5 year period was observed for the 2004–2011 period (Figure 4a). The association between EVI and AAO showed significant coherence for two periods and frequencies: a 2-year period for the 2002–2005 and 2008–2012 periods and for the annual cycle in 2012–2015 (Figure 4b). Strong and significant coherence between EVI and ENSO was observed in southern region for a 2-year cycle in 2002–2003 and 2006–2012 (Figure 4c). Significant coherence was observed in the annual cycles of 2003, 2004–2006 and 2012–2013. Significant coherence between EVI and AAO was restricted to the 2-yr band during 2002–2005 and 2010–2012, and for the annual cycle to the 2013–2015 period (Figure 4d).

4. Discussion

Our analyses of two satellite-derived vegetation indices, NDVI and EVI, showed clear differences in their ability to capture the year-round photosynthetic activity across the northern and southern sectors of the CIS temperate ecosystem. It is known that NDVI is susceptible to sources of errors and uncertainty [21] and EVI has been widely adopted by the MODIS Land Discipline Group (https://modis-land.gsfc.nasa.gov/) to study the properties and dynamics of terrestrial photosynthetic communities [35]. EVI is less sensitive to atmospheric and soil background contamination through differences in vegetation absorption of the red and maximum reflectance in the NIR. Several other sources of bias are associated with satellite estimates of ecosystem functioning, fundamentally through their interaction with different sources of noise. For example, noise can be driven by contamination of the remote sensing signal through its interaction with ozone, water vapor, aerosols and other atmospheric constituents [36]. Also, the complex topography of some parts of the region [35] and the presence of persistent cloudiness and strong winds from the cyclonic depressions associated with the westerlies [13] can also affect the performance of vegetation indices derived from satellites [37]. Despite these challenges, suitable corrections to this type of contamination and other operational errors allow the generation of ecologically relevant measurements of ecosystem functioning at different scales and should be incorporated in future studies of the region [19].

Our results also showed the difference in the dynamics captured by NDVI/EVI, which stem from their algorithmic design. EVI was designed to minimize the effects of the atmosphere and canopy background, and consequently it does not become saturated as NDVI does in areas with large amounts of chlorophyll [21,36]. The design of the algorithm could explain the apparent lack of significant annual oscillations in the wavelet analysis of NDVI for both sub-regions. Differences between EVI and NDVI responses could also originate in the saturation of LAI values because of the particular structure of native vegetation, while EVI can remain sensitive to the canopy architecture of austral broadleafed vegetation [21,38]. Similar remote sensing studies have been carried out in the northern hemisphere at high latitudes. For example, neither NDVI nor EVI were able to correctly capture the onset and termination of the growing season in evergreen needleleaf forests across central north America, which was improved incorporating land surface temperature [39]. Also, a study of NDVI in evergreen boreal forests achieved an improved phenological signal through atmospheric corrections [40] or direct improvements of the algorithm [41]. Further research is needed to extend the use of EVI to assess the phenological contribution of different vegetation components in the temperate rainforest of the Chiloé Island System, and its potential connection and/or synchrony with the productivity of the adjacent marine ecosystem—recent studies suggest that, in this region, there is a modulation of sea surface chlorophyll by the continental input of freshwater and its terrigenous materials [42,43]. The effects of climate on phenological shifts have largely been observed over larger spatial and temporal scales than our study [1].

Although the “phenology-climate connection” operates at different scales, i.e., local-to-global scale, it is intrinsically related to climatic variability [44]. According to different models and observations, primary productivity in high latitude terrestrial ecosystems has experienced an increase during the last 30 years [45]. Moreover, an analysis of NDVI derived from Advanced Very High Resolution Radiometer (AVHRR) data [46] showed different trends in photosynthetic intensity globally and in response to global change, highlighting profound impacts on ecosystem fluxes, including matter and energy [47]. The interplay between ENSO and AAO describes the large-scale distribution of atmospheric mass from the tropics to the extratropics around the Southern Hemisphere [30,48], and thus the supply of maritime moisture to the CIS [13]. During our study period and over the past four decades ENSO and AAO were negatively correlated, a pattern that was briefly interrupted during the 2015 ENSO event, which bears the fingerprint of anthropogenic climate change [14,48]. Positive (negative) phases of AAO (ENSO) are associated with a curtailed supply of moisture, while during the opposite conditions, wetter conditions prevail [13,14,49]. For example, recent studies on the frequency of fires across South America showed a clear influence of the AAO, with positive-AAO conditions associated with drier, more fire-prone conditions across Patagonia [49,50]. The near-neutral ENSO conditions during most of the MODIS mission were disturbed by a few extreme climatic events; the short-lived 2010–2011 La Niña and the strong 2015 El Niño event, together with a mild La Niña in 2008 that was quickly followed by a mild El Niño event in 2009. Chiloé Island System could be used as sentinels phenoregions to monitor the impacts of climate change on these austral broadleaf forests [26,51].

The recent pattern of climatic variability was reflected in the pattern of co-oscillation with vegetation properties across the CIS (Figure 2c) [14]. In both northern and southern sectors there was clear coherence between the annual cycles of EVI and ENSO during the neutral 2003–2006 period, while their coherence with the SAM over the annual cycle was only observed during the period prior to the 2015 ENSO event. Coherence of EVI and both climatic indices was observed over multiple and longer time scales over the 2004–2012 period, particularly for the southern sector. The time-varying relationships between climate and EVI suggest that a phenology-climate association is more evident between EVI and ENSO, particularly in the northern sector, suggesting an important role of tropical atmospheric dynamics on the delivery of moisture to the CIS [12,13]. It is worth noting that the recent coupling of positive phases of the AAO and ENSO generated extreme conditions through high solar irradiation and reduced precipitation and runoff, that cascaded into large harmful algal blooms in the Inner Sea of Chiloé and elsewhere around western Patagonia [14,52]. These disruptions, however, were not evident in the variation of EVI patterns, and were likely ameliorated over the CIS by synoptic-scale processes [14]. Finally, the clear attenuation of the annual cycle in EVI, particularly in the northern sector, deserves further attention from future studies on the impact of a changing climate on the primary productivity of western Patagonia [53]. The temporal extension of our observation does not allow us to identify long-term trends and/or oscillations (e.g., decadal cycles). However, our results indicate that the phenological activity of the CIS clearly corresponds to a temperate pattern, with high (low) activity during the Spring-Summer seasons (Winter months) [54], with clear evidence of interannual coherence with climatic forcing.

5. Conclusions

The EVI showed better performance in capturing phenological patterns compared to NDVI in a highly productive temperate region in northern Patagonia. Quantification of relevant ecosystem processes, such as the temporal variation in primary productivity, provides us with important information for field-based process studies. Our results highlight that despite the maritime weather pattern that dominates the CIS, a marked phenological response is evident in the satellite-based vegetation indices. Our analysis highlighted the connection between the remotely-sensed vegetation index and large-scale climate variability through multiscale coherence with ENSO and AAO. Intense anthropogenic activity across the CIS can modify the observed trends, mainly due to changes in land use following agricultural activities, logging, reforestation and urbanization. Hence we suggest that our results should be taken with caution, as field data are required to validate our observations. Moreover, special emphasis should be paid to validation efforts around forested areas that combine high biomass with high endemism, and in other ecosystems that constitute large reservoirs of carbon, such as peatlands. Establishing the presence of spatial and temporal structure, in association with multiple land uses, will remain an important challenge for future studies.