Impacts of vegetation onset time on the net primary productivity in a mountainous island in Pacific Asia

Vegetation phenology reflects the response of a terrestrial ecosystem to climate change. In this study, we attempt to evaluate the El Niño/La Niña-Southern Oscillation (ENSO)-associated temporal dynamics of the vegetation onset and its influence on the net primary productivity (NPP) in a subtropical island (Taiwan) of Pacific Asia. We utilized a decade-long (2001–2010) time series of photosynthetically active vegetation cover (PV) data, which were derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance data, to delineate the vegetation phenology. These data served as inputs for the phenological analysis toolbox TIMESAT. The results indicated that the delayed vegetation onset time was directly influenced by a dry spring (February and March) in which less than 40 mm of rainfall was received. This seasonal drought impeded vegetation growth in the subsequent growing season, most likely due to delayed impacts of moisture stress related to the preceding ENSO events. The significant correlations obtained between the annual MODIS NPP and both the vegetation onset time and the length of the growing season may imply that the accumulated rainfall in the spring season governs the annual NPP. The model simulations revealed that the frequency and intensity of the ENSO-related spring droughts might increase, which would result in cascading effects on the ecosystem metabolism.


Introduction
Net primary productivity (NPP) is the net carbon (C) dioxide retained in vegetation from the atmosphere. NPP quantifies the production of new plant material, including new biomass, hydrocarbon emissions, root exudates and transfers to mycorrhizae (Chapin et al 2002). NPP is a temporally accumulative variable and is often used to measure the C net uptake in a calendar year. Moreover, NPP varies across time and space, and its interannual dynamics are Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. strongly influenced by climatic variations (Potter et al 2011), vegetation, land use (Cao et al 2004) and disturbances (Huang and Asner 2010). This measurement is often used as a salient indicator to monitor the metabolism of an ecosystem (Crabtree et al 2009).
Remote sensing has been frequently utilized to monitor vegetation phenology and ecosystem productivity over large regions (Running et al 2004, Jeong et al 2011. Several lines of evidence suggest that temperature and precipitation are two major climatic factors that govern NPP spatiotemporal dynamics (Nemani et al 2003, Zhang et al 2004. The interannual variability of vegetation greenness may also be subject to large-scale irregular, however crucial, climatic anomalies due to El Niño/La Niña-Southern Oscillation (ENSO) cycles (Wharton et al 2009). Studies have shown that the earlier start and the increasing length of the growing season are positively correlated with the spatial patterns of temperature during spring and autumn in mid-and high-latitude regions (Chen et al 2005, Jeong et al 2011. During ENSO warm events (El Niño), warmer and greener conditions dominate in the spring over North America, Far East Asia and, to a certain extent, central Europe (Buermann et al 2003). In contrast, the extended droughts associated with El Niño at 5-to 10-year frequencies in Amazonia have formed a more spatially homogeneous greenness surface. However, there is more spatial variability during La Niña due to higher rainfall (Poveda and Salazar 2004). These events can cause decadal-scale impacts on the regional C cycle and a persistent change in the forest canopy structure (Saatchi et al 2013). However, the influences of ENSO-related events on tropical and subtropical islands, which are the major biodiversity hotspots (Kier et al 2009), have not been widely studied (Nagai et al 2007). To date, the ecosystem responses to regular and irregular climate variabilities remain largely unknown, which constitutes a gap in the understanding of their adaptation and feedback under changing climate conditions.
The ENSO-related events that have occurred over the past half-century have strongly influenced the climate of the northwestern Pacific (Chen et al 2008). The preceding cold season El Niño (November-February) is often associated with a subsequently wetter spring (February-March) with more heavy rainfall events and a greater amount of plum rain (May-June) in Pacific Asia and vice versa for La Niña conditions . In this region, ENSO-related events can affect the water availability in spring and subsequently impact the ecosystem productivity in the summer growing season because the vegetation growth and development are sensitive to the preseason timing and amount of precipitation (Tao et al 2004, Naylor et al 2007, Woodward et al 2008. Therefore, we utilized a decade-long (2001-2010) record of remotely sensed phenological metrics (the vegetation onset and offset times and length of the growing season), NPP, temperature, and precipitation data to investigate their impacts on ecosystem productivity. The vegetated areas of Taiwan were selected as the study region due to their spatial extent, which covers tropical and subtropical regions and exhibits a wide topographical (bioclimatic) gradient with diverse vegetation types. This area provides a valuable opportunity to investigate the vegetation responses to regular and irregular climatic variations.

Study area
Taiwan, which lies across the Tropic of Cancer, is a 36 000-km 2 tropical/subtropical mountainous island near the continent of Eurasia ( figure 1(a)). The elevation ranges from 0 to 3952 m a.s.l. The mean (± standard deviation, sd) annual air temperature and annual precipitation are 22 ± 0.25 • C and 2500 ± 560 mm, respectively (Shih 1991). Taiwan experiences prevailing northeast monsoons in winter and southwest monsoons in summer (Yen and Chen 2000). Summer (May-October), which is the primary growing season for the region, receives 75% of the annual precipitation with warm temperatures up to 29 • C. Spring is relatively dry over southwestern Taiwan to the leeside of the Central Mountain Range due to the prevailing northeast monsoons.
Five mountain ranges run from north to south and bisect the island; there are more than 200 peaks over 3000 m a.s.l. in Taiwan. The dominant vegetation types along the elevation gradient are farmland on the coastal plains and tablelands, evergreen broadleaf forests at low and mid-elevations (200-2000 m a.s.l.), and mixed and conifer forest at mid-and high elevations (1100 m a.s.l. and above) ( figure 1(b)). Approximately 60% (21 600 km 2 ) of the land is occupied by forests consisting of natural (73% of the total forested areas), plantation (20%), and bamboo (7%) forests. The developed lands (farmlands (29%, 10 400 km 2 ) and urbanized areas (6.1%, 2200 km 2 )) are generally situated below 800 m a.s.l. (Chang et al 2011).

Climate data
Monthly continuous temperature (n = 136) and precipitation (n = 390) data during the period 1981-2010 were obtained from the Data Bank for Atmospheric Research of Taiwan (http://dbar.ttfri.narl.org.tw/) ( figure 1(b)). The monthly temperature anomalies (difference from the long-term mean) were calculated from interpolated gridded surface for the study region (Chang et al 2011). The monthly precipitation (MP) surfaces were generated using the ordinary kriging function with a spherical model (ArcGIS v. 9.3, Environmental Systems Research Institute, Inc., CA, USA). Kriging is a geostatistical method for the estimation of the optimal unbiased environmental variables at an unsampled location and can be formulated as: where Z(MP u ) is a kriged estimate of the MP at the unknown location u, λ ui is the weight for the observed locations i, which sum to unity to ensure that no bias will minimize the estimation variance, and Z(MP i ) is the observed MP at location i (Zhang and Srinivasan 2009). Studies have shown that this approach can be an effective tool for the prediction of climate-related spatial uncertainties (Diodato 2005, Vienneau et al 2009. The spatial resolution of these MP data was set to 1 km × 1 km. A dry spring (February and March) was defined as an area receiving less than 40 mm following Chen et al (2009). To delineate the interannual precipitation variability in Taiwan during the study period, the standardized precipitation index (SPI, Hayes et al 1999) was generated using MP during the period 1981-2010. The SPI was defined as the transformation of the precipitation time series into a standardized normal distribution (z-distribution): where MP i , MP, and MP σ are the precipitation for month i and the mean and the standard deviation of MP over a specific time period, respectively (Llyod-Hughes and Saunders 2002). Negative and positive SPI values indicate dry and wet conditions, respectively. To properly identify the drought conditions in the region, we used the SPI with a three-month time scale (SPI3) as recommended by Chen et al (2009). The southern oscillation represents the fluctuation of the measured atmospheric pressure at mean sea level between Darwin in Northern Australia and Papeete in Tahiti. The Southern Oscillation Index (SOI) is commonly used to assess ENSO, with negative (<−1.0) and positive (>1.0) values indicating warm (El Niño) and cold (La Niña) events, respectively (Jin et al 2005). Ropelewski and Halpert (1987) found that the dynamics of the SOI are associated with precipitation anomalies in tropical and subtropical regions. The monthly SOI for the period 1981-2010 was directly downloaded from the NOAA Climate Prediction Center (www.cpc.ncep.noaa. gov/) to investigate its relationship with the MP in Taiwan.

MODIS data
We applied spectral mixture analysis on the MODIS (the Moderate Resolution Imaging Spectroradiometer) eight-day L3 500-m surface reflectance time-series (2001-2010) product (MOD09A1) to delineate the vegetation phenology. Two tiles (H28V06 and H29V06) were required to cover the entire study area. The method has been commonly used to decompose image pixels into main surface components, e.g., photosynthetically active vegetation (PV), non-photosynthetically active vegetation (NPV), and soil and rock outcrop (SRO) (Adams et al 1993). This method assumes that the reflectance value of a pixel is a linear combination of the pure-component cover fractions (i.e., the endmembers). An automated probability-based spectral mixture analysis approach with three spectral endmember bundles (i.e., PV, NPV, and SRO) was used to obtain the proportions of sub-pixel cover fractions for each pixel (for details, see Asner and Heidebrecht 2002). A set of spectral libraries (n PV = 580, n NPV = 267 and n SRO = 256) was required for the spectral unmixing, which was collected using a field spectroradiometer (Analytical Spectral Devices, Inc., Colorado, USA) or a spaceborne hyperspectral imager (Hyperion, http://glovis.usgs.gov/). These measurements were convoluted to match the MODIS spectral profiles (see Huang et al 2013 for technical details). Monthly PV images (n = 120) were selected from 460 images (3-4 images/month) between 2001 and 2010; the criterion was based on a maximum value composite procedure (Huete et al 2002). We used PV instead of NDVI, or other commonly used vegetation indices, e.g., the Enhanced Vegetation Index (EVI), because it was more intuitive and very sensitive (r = 0.70-0.92, Chang et al 2013) to the monthly mean temperature and precipitation for different land-cover types.
All of the images were resampled to a resolution of 1 km for the following phenological analysis. Annual MODIS NPP (MOD17A3) 1-km products from the Terra platform for the period 2001-2010 were acquired from the NASA Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/lpdaac/get data/data pool) (Running et al 2004). The quality of the product was satisfactory and comparable with field observations across different biomes (

Phenological analysis
The phenological metrics (the onset and offset times and the length of the growing season) were derived from 10-year MODIS PV time series using TIMESAT v. 3.0 , which is a software package that finds the best mathematical fitting function to characterize vegetation temporal dynamics (figure 2). The Savitzky-Golay smoothing filter was selected because it exhibited superior performance in the preservation of the vegetation temporal dynamics and a minimization of the atmospheric effects . Using PV data and several specific land-cover settings, we can obtain ten full phenological cycles (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and derive the major phenological metrics of Taiwan. The start (onset)/end (offset) of the growing season (presented as the Julian date [JD]) was defined by referring to a fitted function as the point in time setting to 10% of the distance between the left/right minimum and the maximum, respectively. The approach has been implemented in previous studies (Brady et al 2007, Eklundh et al 2009, Kariyeva and van Leeuwen 2011. The length of the growing season was determined by the time from the start to the end of the growing season (figure 2, Jönsson and Eklundh 2004). The relationships among the vegetation onset time, the vegetation offset time, the length of the growing season, ENSO-related spring rainfall, and the interannual dynamics of MODIS NPP were investigated.

Climate variability and PV time series
Compared with the average of the recent three decades , the mean annual temperature of Taiwan was  3(b)). However, the SOI had no significant correlations (p > 0.05) with monthly temperature anomaly nor SPI3. The pattern of preceding La Niña occurrences in the cold season (November-February) during 2000-2001, 2001-2002, 2007-2008, and 2008-2009 and a drier spring (February-March) was observed in 2001, 2002, 2008, and 2009. In contrast, the preceding El Niño events in the cold season during 2004-2005, 2006-2007, and 2009-2010 were followed by a wetter spring in 2005, 2007, and 2010 (figure 3). In the period 1981-2010, the most apparent correlations were between the SOI from the preceding September to the following March and monthly rainfall in February-March (r < −0.50, p < 0.05, figure 4). We also observed significant correlations between the November-February SOI and the subsequent July rainfall and between September-October SOI and October rainfall of the same year ( figure 4).

Spatial patterns of phenological metrics
The dates of onset in Taiwan ranged from approximately JD 80 (late March) to 160 (early June); the earlier and later dates appeared in the south plain areas and the mid-altitude regions (∼1500 m a.s.l.), respectively (figure 5). The subsequent offset dates were from JD 240 (late August) to 340 (early December). The earlier offset dates were distributed in the western coastal plain (primarily farmland); the later ones were mainly observed in the southwestern and southeastern regions (forests, figure 5). For the entire study region, the reduced spring rainfall led to a later vegetation onset (r = −0.60, p < 0.05, figure 7(a)) and a shorter growing season (r = 0.65, p < 0.05, figure 7(b)). The areas receiving sufficient spring rainfall with an earlier onset ( figure 7(a)) generally results in a longer growing season (r = −0.77, p < 0.01, figures 5 and 7(d)), and vice versa. However, there were no significant correlations (p > 0.05) between offset time and spring rainfall and the length of growing season. The average length of the growing season at mid-altitudes (210 days) was approximately 40% longer than that at the western coastal plain and at high altitudes (<150 days, figure 5).
We further investigated the spatial pattern of vegetation onset and its relationship to spring rainfall during the period 2001-2010 (figure 5). The areas of delayed onset (late May, ≥JD 140) were similar (the agreement was approximately 51% and up to 79%) to those areas that received spring rainfall of less than 40 mm in central and southwestern Taiwan (dry zones). The delayed onset areas occupied 33% of Taiwan in 2001 and extended to a maximum of 38% in 200238% in . From 200338% in to 2006, the percentage of delayed onset regions decreased substantially from 30% to 13% and then increased from 25-28% in 2007-2009 to 32% in 2010 (figure 5).
Along the elevation gradient, the NPP generally increased from the coastal plain (approximately 5-450 m a.s.l., <4.0 Mg C ha −1 yr −1 ) to mid-elevations (300-2200 m a.s.l., 8.0-10.0 Mg C ha −1 yr −1 ), decreasing at high elevations (>2700 m a.s.l., <3.5 Mg C ha −1 yr −1 , figures 1(b) and 5). The maximum annual NPP was observed in southern Taiwan (12.0 Mg C ha −1 yr −1 ). The variation (sd) in the C net uptake was higher in southwestern Taiwan (1.5 Mg C ha −1 yr −1 ) than northeastern Taiwan (0.5 Mg C ha −1 yr −1 ). For the entire island, significant (p < 0.05) positive and negative correlations between MODIS NPP and spring rainfall and the mean onset time were found, respectively (figures 7(c) and (e)); but there was no significant correlation (r = −0.10, p = 0.61) between annual NPP and the mean offset time. As expected, a pronounced positive correlation (r = 0.88, p < 0.01, figure 7(f)) was also observed between the annual mean NPP and the mean length of the growing season. However, a weak relationship (r = 0.35, p = 0.28) between annual NPP and precipitation in the summer growing season was found.

Climate fluctuations and PV
The global circulation has substantial effects on regional climate conditions. In this study, although there were no significant relationships between the SOI, temperature anomalies, and the SPI3 (dryness), the preceding dominant La Niña (El Niño) events in the cold season induced drier (wetter) springs in Taiwan (figure 3, also see . The time-lag effects of ENSO events in the preceding cold season significantly altered the following spring rainfall over Taiwan and affected the vegetation growth (figures 4 and 5). The abundance of PV in the study region reflects the complex interaction of human intervention and climate. Generally, at low elevations, the physical environment is very suitable for plant growth with sufficient precipitation and warm temperatures throughout the year. However, this flat region yielded the lowest annual mean PV (figure 5), primarily due to urbanization and agricultural practices. The precipitation increased along the elevation gradient up to ∼2500 m a.s.l., likely due to the orographic effect, which also reflects the abundance of PV. At higher elevations (>2700 m a.s.l.), the low temperatures downregulate plant growth; the alpine ecosystems are dominated by grasslands or mountain dwarf bamboo forests with relatively low PV (figures 1(b) and 5). However, the interannual variability of the dominant vegetation types along the elevation gradient is quite synchronous (figure 6) and possibly associated with seasonality and the availability of water, i.e., the lowest annual mean PV occurred in 2002, which resulted in long-lasting droughts (1520 mm yr −1 ), compared with the wetter years during the period 2005-2008 (2800-3500 mm yr −1 ) (figures 3, 5 and 6). Our findings illuminate the sensitivity of tropical/subtropical ecosystems to climatic fluctuations and the time-lag of accumulated rainfall (Prasad et al 2007, Meir andWoodward 2010).

ENSO-related spring rainfall
Previous studies indicated that ENSO-effected climate variables significantly regulated interannual and spatial variability in vegetation phenology and NPP (Nemani et al 2003, Van der Werf et al 2006. In temperate northeast Asia (42 • N), the years with an earlier spring onset time related to the ENSO-induced warmer temperatures (March-May) in stronger El Niño years revealed that air temperature, rather than precipitation, affected vegetation phenology (Park et al 2012). In the tropical Amazon forest (7 • N-14 • S), severe drought during El Niño periods led to large C emissions by elevating forest flammability, mortality and inhibiting plant growth (Asner et al 2000, Nepstad et al 2004. The interannual variations of vegetation activities over tropical rainforest regions (5 • N-10 • S) were also dominated by the availability of water associated with ENSO events (Nagai et al 2007). As a whole, the interannual variations in ENSO-related rainfall are crucial for ecosystem development across tropical and subtropical regions (Myneni et al 1996, Erasmi et al 2009. Potts (2003) found that the tree mortality rates in tropical rainforests in Southeast Asia were higher (7.6%) than the background mortality (2.4%) after the severe drought in 1997-1998 associated with El Niño events. Although no large-scale tree dieback events were reported, water shortages can still be a limiting factor for plant growth in this tropical/subtropical mountainous island, particularly in the central and south regions during the dry season.
We found that the spatial dynamics of a delayed onset time (i.e., after JD 140) corresponded to areas that received spring (February-March) rainfall of less than 40 mm (less than 53.3% of the normal rainfall). The timing and amount of ENSO-related spring rainfall, in conjunction with warmer temperatures, affected the bioclimatic regimes and subsequently delayed plant growth (

Interannual variability of NPP
Ecosystem productivity is primarily regulated by temperature, precipitation and solar radiation (Nemani et al 2003, Running et al 2004). In this study, the interannual variability in the NPP was approximately 28%; the lowest productivity occurred in 2002 (23.3 Tg C) and the highest in 2010 (29.7 Tg C) (figure 5), which may be beyond the normal climatic fluctuations (Running et al 2004). This could be the result of intertwining El Niño and La Niña events during the observation period. The global-scale simulation showed that the terrestrial NPP is negatively correlated with El Niño events in the tropics, which suffered drier and warmer climatic conditions (Cao et al 2004, Mohamed et al 2004. This result was confirmed by observational studies conducted in South and East Asia (Saigusa et al 2008, Bala et al 2013. Moreover, preceding El Niño (La Niña) events typically resulted in a subsequent wetter (drier) spring in Taiwan (i.e., time-lag effect) (figures 3 and 4, also see ). These biophysical memory effects caused by the ENSO-related rainfall patterns (figure 4) would amplify the vegetation dynamics and affected the vegetation phenology and productivity in the following year (figures 5 and 7). It is surprising that no strong relationship was found between NPP and the primary growing season precipitation in this tropical/subtropical mountainous island. This region receives more than 1500 mm of rainfall, even in a dry year ; this amount should be sufficient to sustain most of the biome types (Knapp and Smith 2001). In addition, approximately 40% of summer precipitation is from typhoon-associated short, high-intensity rain events (Chen and Fan 2003), which provide little contribution to plant growth in this mountainous island. A majority of the water is exported to the ocean in a very short time period (Huang et al 2012). This illuminates the pivotal role of spring rainfall in net C uptake in this bioclimatic region.
Drought-induced damages are pervasive and not only limited to forests in water-limited arid and semi-arid regions (Allen et al 2010, Misson et al 2011. In tropical forests, prolonged ENSO-related droughts may permanently alter the forest structures and lead to substantial effects on the NPP , Nepstad et al 2004, Van der Werf et al 2006, Brando et al 2008. More frequent and intense El Niño events may occur in the foreseeable future in many humid tropical and subtropical regions (Tsonis et al 2005), which will likely result in recurrent droughts that may impede the vegetation recovery (Saatchi et al 2013). This result would increase the interannual fluctuations in the vegetation phenology and consequently the interannual variability of the NPP and C budget.

Conclusions
In this study, we derived vegetation phenological metrics (the onset and offset times of the summer growing season and the length of summer growing season) for Taiwan using monthly composites of PV derived from MODIS surface reflectance data. The results indicated that the delayed vegetation onset time was directly influenced by the ENSO-related dry spring (February-March) (i.e., rainfall amount of less than 40 mm). Our remote sensing analysis indicated that a variation in the vegetation onset time might have collateral effects on the length of the summer growing season and annual NPP. Future climate predictions suggest that the frequency and intensity of extreme ENSO-related meteorological events in tropical and subtropical regions might be amplified, which could result in cascading effects on the terrestrial C budget in the near future. Woodward Running S W 2005 Improvements of the MODIS terrestrial gross and net primary production global data set Remote Sens. Environ. 95 164-76 Zhou L, Kaufmann R K, Tian Y, Myneni R B and Tucker C J 2003 Relation between interannual variations in satellite measures of northern forest greenness and climate between 1982-1999 J. Geophys. Res. 108 4004