Shifts in Growing Season of Tropical Deciduous Forests as Driven by El Niño and La Niña during 2001–2016

Abstract

This study investigated the spatiotemporal dynamics of tropical deciduous forest including dry dipterocarp forest (DDF) and mixed deciduous forest (MDF) and its phenological changes in responses to El Niño and La Niña during 2001–2016. Based on time series of Normalized Difference Vegetation Index (NDVI) extracted from Moderate Resolution Imaging Spectroradiometer (MODIS), the start of growing season (SOS), the end of growing season (EOS), and length of growing season (LOS) were derived. In absence of climatic fluctuation, the SOS of DDF commonly started on 106 ± 7 DOY, delayed to 132 DOY in El Niño year (2010) and advanced to 87 DOY in La Niña year (2011). Thus, there was a delay of about 19 to 33 days in El Niño and an earlier onset of about 13 to 27 days in La Niña year. The SOS of MDF started almost same time as of DDF on the 107 ± 7 DOY during the neutral years and delayed to 127 DOY during El Niño, advanced to 92 DOY in La Niña year. The SOS of MDF was delayed by about 12 to 28 days in El Niño and was earlier about 8 to 22 days in La Niña. Corresponding to these shifts in SOS and LOS of both DDF and MDF were also induced by the El Niño–Southern Oscillation (ENSO).

1. Introduction

Tropical forests contain about 25% of the carbon in the terrestrial biosphere and account for 34% of Earth’s gross primary production [1]. Unlike temperate forests where temperatures fluctuate widely during the course of a year, their variation in tropical forest is modest. The trees of tropical forest are thus adapted to grow in a relatively narrow temperature range. Hence, the relative impact of climate warming is likely greater in the tropics than in other regions because predicted changes in temperature are large compared to normal inter-annual variations [2]. In addition, changes in precipitation patterns such as a shift toward more extreme events and extended droughts under climate change may result in loss of tropical forest and in large amount of CO₂ released to the atmosphere [3,4].

Among various types of tropical forests, tropical deciduous forest occupies about 43% of the forest area in the tropical belt with great diversity of species [5]. It provides valuable services involving biodiversity, water resources and carbon sinks. However, like other forest ecosystems, these services are being affected by climate change and variability. Cavaleri [6] for example, reported a decline in the carbon sink of tropical rainforest during El Niño because of reduced photosynthesis and increased respiration rates. Strong climatic disturbances can severely reduce forest biomass, and if the frequency and intensity of these events increases beyond historical averages, these changing disturbance regimes have the capacity to significantly reduce forest biomass, resulting in a net source of carbon to the atmosphere. The study on Atlantic tropical moist forest in a long-term experimental plot also indicated a rapid biomass decline associated with El Niño events [7]. Many others studies have indicated that strong El Niño events have negative impacts on forest ecosystems, which could result in significant increasing level of tree mortality, changing plant phenology and carbon flux [3,4,8,9,10,11,12].

Forest phenology including the start, end and duration of growing season is an important indicator of vegetation response to climate change [13]. However, the information on tropical deciduous forest response to such climate extremes is sparse. Preliminary assessment over Southeast Asia suggests that decreasing precipitation and unusual high temperature in relation to severe drought (El Niño event) results in a significant reduction the CO₂ uptake [14]. Such reduced CO₂ uptake in tropical forests may be caused by reduced photosynthetic activity, shortened growing season, or a combination of both [6,15]. Further improvements in our understanding of tropical forests responding to climatic drivers and its links to ecosystem functions is thus needed.

In this study we applied the Normalized Difference Vegetation Index (NDVI), one of the most widely used indices, to quantify the interannual variation in canopy phenology of a tropical deciduous forest in Northern Thailand and to evaluate canopy response to extreme climate events. The indices derived from remote sensing were compared against in situ observations of canopy phenology. The insight gained through such observations will help to improve our understanding of key feedback mechanisms and our ability to predict vegetation dynamics under climate changes scenarios.

2. Materials and Methods

2.1. Study Area

Lampang Province is located in Northern Thailand with an area of 12,534 square kilometers (Figure 1). The Asian monsoon with two distinct seasons governs the climate in Lampang; a wet (May–October) and a dry season (November–April). The mean annual precipitation during the study period 2001–2016 is 1231 mm. The mean annual temperature is 25 °C, with a monthly minimum mean of 21–22 °C occurring in December–January, and a maximum monthly mean of 29.5 °C occurring in April [16].

Forest cover during the study period was on average 68.5% of total province area. Deciduous forest accounts for 68.3% of the overall forest coverage (Figure 1) [17]. Lampang is located on the area with the altitude ranging from 114 m to 1939 m above mean sea level. The main forest types in this study area are deciduous forests, mixed deciduous forest, and evergreen forest [18]. We selected Lampang as our study area for three reasons. Firstly, Lampang has a climate regime that is representative of northern Thailand and has experienced El Niño events during study period [19]. Secondly, there are long term field observations of leaf area index (LAI) in the Mea Mo teak plantation (2001–2012). Thirdly, the study area consists of a complex mountainous topography and the forest diversity varies according to different elevation. For example, lower elevation (300 m to 800 m) is dominated by the dry dipterocarp forest and Tectona grandis L. f (teak), and the higher elevation (800 m to 1000 m) is dominated by the mixed association of deciduous and evergreen hardwood, and the region above 1000 m is dominated by primary and evergreen with Pinus kesiya Roy. Ex Gord. (Pinaceae, pine) [20]. The complex topography, forest ecosystem and functioning make it ideal for evaluating the spatial variation in growing season of deciduous forest.

2.2. Data Set

2.2.1. MODIS Fata

This study used surface reflectance of MOD09Q1 from MODIS bands 1 (Red) and band 2 (NIR) at 250 m resolution, and MOD09A1 at 500 m resolution imagery captured in an 8-day period (2001–2016). The data were downloaded from the EROS Data Center, US Geological Survey [21,22]. Cloud cover was present in MOD09Q1 images, which limited the potential of the images for ground information extraction. Removing and replacing cloud-contaminated pixels was then performed following the method of Hoan and Tateishi [23]. Cloud removal processing for band 1 and band 2 of MODIS MOD09Q1 was divided into three main steps: cloud mask extraction, interpolation to remove cloud coverage and median smoothing. For cloud mask, combination of internal cloud algorithm flag (bit 10) and cloud shadow (bit 2) from 500 m State Flags layer of MOD09A1 (500 m) product was selected. The interpolation method for replacing cloud-contaminated areas was applied for each pixel (see detail in [23]). The NDVI was then calculated by using Red and NIR band [24].

2.2.2. Landsat Imagery

The dataset Operational Land Imager (OLI-8) images were downloaded from U.S. Geological Survey [25] by using the Global Visualization Viewer (GloVis). The images with less than 10% cloud cover were chosen to minimize the effect of cloud to interpreting forest areas. The red, green, blue, NIR, SWIR-1, and SWIR-2 bands were used in the forest classification. List of LANDSAT imagery used in this study is shown in

Table 1. Specifications of Operational Land Imager-8 (OLI) imagery used in this study.

No	Landsat Path/Row	Date Acquired
1	Landsat OLI-8 130/47	16-01-2017
2	Landsat OLI-8 130/47	16-01-2017
3	Landsat OLI-8 131/47	08-02-2017
4	Landsat OLI-8 131/47	08-02-2017

.

2.2.3. Digital Elevation Models

The Shuttle Radar Topography Mission (STRM) 30 m × 30 m resolution obtained from the U.S. Geological Survey (USGS) were used in this study for topographic correction [26]. STRM was resampled using a nearest neighborhood transformation [27,28].

2.2.4. Climate Variables

Local climate variables including daily maximum air temperature and precipitation were obtained for the period of 2001–2016 from Thai Meteorological Department. There are total of 123 meteorological stations over Thailand. However, there is only one station located within the study area (18°17′ N; 99°31′ E), and this was used as a representative of climate in Lampang. It should be noted that this station was only used to evaluate climatic trends and variation of the study area. It is not intended to characterize the climate of the whole region. The internal consistency and temporal outliers check were performed for climate data set. The missing value and outlier values on data set was checked and removed [29,30]. The daily maximum temperatures and precipitation were averaged and aggregated. Annual anomaly was calculated as the difference between the 16-year average and the individual year [31,32]. A positive anomaly value indicates that the observed value was greater than the average, while a negative anomaly indicates that the observed value was less than the average for the period 2001–2016.

2.3. Methods

This study consists of three main parts: forest classification, evaluation of the difference between in situ and satellite based phenological metrics, and the influence of the El Niño–Southern Oscillation (ENSO) to phenological metrics variation. A detailed description of the data, methods and outputs used in the study is provided in Figure 2.

2.3.1. Forest Classification in Lampang

For pre-processing LANDSAT imagery, the top-of-atmosphere (TOA) reflectance and atmospheric correction (dark pixel subtraction method) were applied [33]. The cloud, cloud shadow and water pixels were identified by applying the Function of Mask (Fmask) algorithm [34,35] and were excluded from further analyses. Additionally, Landsat scenes contain jagged pixels along scene edges, which could result in incorrect reflectance values along the edge of imagery. These jagged pixels were removed using a 450 m buffer applied from the edge of the mask inward [36,37]. The Otsu threshold was then performed based on NDVI calculated from Red and NIR band of Landsat to separate the forest and non-forest areas [38]. The topographic correction with Statistical Empirical Correction method was performed on each stratified forest area to remove the topographic relief effect from Landsat imagery [28]. The forest object was then used to classify the forest type [39].

The stratified random sampling approach was performed in this study to estimate the number of samples points for each class [28,40]. The selection of a random sample of training data across all land use types ensures that samples have class proportions representative of each forest type. In total, there were 458 sample points selected for all classes (Figure 3). These 458 points included 50 points for evergreen forest (EF), 133 points for mixed deciduous forest (MDF), 50 points for dry dipterocarp forest (DDF), and 225 points for other vegetation types and non-forest class. A total of 40 of the 458 points were validated during the field survey. The rest were selected manually using a combination of the time series of MODIS NDVI during 2001–2016, the interpretation of high-resolution imagery, Google Earth images, and aerial photographs in Quantum GIS. To avoid the spatial autocorrelation a minimum distance of 2000 m was required between selected points [28]. A minimum of 50 training samples were used for each forest class [28].

The Random Forest (RF) Classifier, a widely used non-parametric machine learning classifier, was then applied for forest classification. RF is based on a tree classifier and grows many classification trees [28,39]. From the trained samples, the random Forest and raster packages were used to map forest types from the set of 06 spectral predictor variables (n = 500 classification trees). For the RF classification, 70% of sample points (320 points) were used to train the model, while the remaining 30% of sample points (138 points) were used to validate the classifier (Figure 3). The overall accuracy and Kappa statistic were considered for assessing the accuracy of RF classifier. The two major of classified deciduous forest types including MDF and DDF were used as mask to analysis their phenology in response to climate extreme event.

2.3.2. NDVI Time Series for Phenological Analysis

To examine the temporal signatures of deciduous forest classes, we extracted Normalized Difference Vegetation Index (NDVI) from MODIS to produce the profiles of vegetation dynamics during each growing season. NDVI was calculated from Red and Near infrared as described by [24]. The remaining noise in NDVI time series data was removed by using the Savitzky-Golay built within the TIMESAT software. The smoothed NDVI series was used to extract phenological metrics by using the threshold-based method [41,42,43,44]. In this study, we examined mean value of NDVI time series during 2001–2016 of two different forest types: DDF and MDF. The average of NDVI pixel locations of DDF and MDF was confirmed by checking the coordination from field survey.

2.3.3. Determination of Forest Phenological Variables from Satellite

TIMESAT software was used to extract the tropical deciduous forest phenology metrics for each year (2001–2015). An adaptive Savitsky-Golay filtering function in TIMESAT software was used to fit the curve of time series NDVI, which has been widely used for phenology monitoring [41,43,45,46]. Three phenological metrics were extracted as follows:

• Start of growing season (SOS): This is defined as the date of leaf unfolding (day of year, DOY) and this study considered SOS as a date when NDVI of the left edge has increases 20 percentage measured from the left minimum point.
• End of the season (EOS): This is defined as the dates of leaf discoloration (day of year, DOY) and leaf fall at the end of season. This study considered EOS as a date when NDVI of the right edge has decreases to 20 percentage of the right minimum level.
• Length of the season (LOS): This is the duration (number of days) from the start to the end of the season.

The analysis of phenological metrics was performed for each pixel in this study. Note that one growing season cycle of tropical deciduous forest starts from current year and ends in the following year. Therefore, only 15 years of phenological metrics were extracted from NDVI time series of 2001–2016. To avoid long-term changes pixels classified as deciduous forest, only pixels for which the SOS ranged between 15 and 180 DOY for the whole study period were selected for further analysis.

2.3.4. Determination of in situ Derived Forest Phenological Metrics

The daily in situ LAI was aggregated to 8-day temporal resolution by using the mean function, with the same temporal resolution to MODIS NDVI. The Savitsky-Golay filtering was also performed with LAI data set. The phenological metrics were then extracted by the same method of satellite-based phenological metrics.

2.3.5. Validation of the Phenological Metrics

The ground observation site was located in Mae Mo District, Lampang Province (18°25′ N, 99°43′ E), in teak plantations (Tectona grandis Linn. f.) which was a major forest type in Lampang. In situ observation of leaf area index (LAI) from January 2001 to December 2012 was used to validate the time series of NDVI derived from satellite imagery [47]. To avoid positional error in MODIS data, the 3 × 3 pixels NDVI time series was extracted corresponding to the center of in situ measurement location. The daily in situ LAI was aggregated to 8-day temporal resolution by using the mean function. The phenological metrics extracted from MODIS and LAI are henceforth referred to as ${\mathrm{SOS}}_{\mathrm{NDVI}}$, ${\mathrm{LOS}}_{\mathrm{NDVI}}$ and ${\mathrm{SOS}}_{\mathrm{LAI}}$, ${\mathrm{LOS}}_{\mathrm{LAI}}$. To understand more relationship between LAI and NDVI in term of phenology, we extracted the phenological metrics by using same threshold-based of 20 percentage amplitude. Linear regression was used to examine the overall relationship between the phenology metrics from in situ measurements and satellite-derived estimates.

2.3.6. Assessing the ENSO-Related Patterns in Annual Phenological Metrics

To examine the influence of ENSO events on tropical deciduous forest vegetation dynamics, the difference of phenological metrics during El Niño and La Niña years were evaluated for the study area. The extreme climate and neutral years were determined using the following criteria: (1) any year with five consecutive Oceanic Niño Index (ONI) periods in excess of +0.5 °C (−0.5 °C) is an El Niño (La Niña) year [48], (2) comparing against the annual anomaly climate at Lampang (see detail in Section 3.1). To evaluate the influence of ENSO events on canopy phenology, the difference between El Niño and La Niña impact patterns were evaluated for the study area using a simple difference between average El Niño or La Niña year and neutral year SOS and LOS. To investigate potential differences in ENSO impacts on different forest types, these differences were grouped by forest type and the significance of difference was compared based on one standard deviation (±1SD) and analysis of variance (ANOVA) at the confidence level of 95% (p ≤ 0.05).

3. Results and Discussion

3.1. Variations in Temperature and Precipitation during 2001–2016

A clear seasonality of temperatures and precipitation is a main feature of climate characteristic in Northern Thailand, including in Lampang (Figure 4). Three clear seasons are as follows: (1) rainy or southwest monsoon season (mid-May to mid-October) with August to September as wettest period; (2) cool or northeast monsoon season (mid-October to mid-February), is the mild temperature period of the year with the coolest period in December and January; and (3) summer or pre-monsoon season (mid-February to mid-May). This is the transitional period from the northeast to southwest monsoons. March to May is the hottest period of the year with maximum temperatures often near or exceeding 40 °C. The subsequent onset of rainy season significantly reduces the temperatures from mid-May, which leads to intensive rainfall from mid-May until early October.

Based on the criterion that any year with five consecutive the Oceanic Niño Index (ONI) periods, beginning with the preceding August through October (ASO) period and ending with February through April (FMA) of the current year, El Niño or La Niña years are defined as any year excess of +0.5 °C (−0.5 °C) [48]. All other years are defined as neutral years. Following this criteria, the years 2001, 2006, 2008, 2011, 2012 are considered as La Niña, the years 2003, 2005, 2007, 2010, 2015, 2016 are El Niño year (with 2010 and 2016 defined as strong El Niño years).

Climate data analysis at Lampang station shows that the maximum temperature anomaly (+1.04 °C) was observed during the El Nino year (2010), which also had a pronounced negative precipitation anomaly (−133.9 mm). Precipitation anomalies during the period of study were mostly negative, particularly in El Niño events, whereas the precipitation anomalies values were highest in 2006 and 2011 (which indicates the wettest conditions) during the study period (Figure 5). Based on these results we identified 2010 and 2011 as strong El Nino and La Nina years, respectively, for use as extreme events in the impact analysis. The remaining years (2001–2009 and 2012–2015, including the weak El Niño and La Niña) with less variation of anomaly temperature and anomaly precipitation were considered as neutral years in this study. Two phenological metrics during neutral and extreme event years were extracted for each pixel, then the difference between them was investigated.

3.2. Forest Classification

RF classification for various different forest types from OLI-8 is presented in Figure 6. The classification of land use in Lampang 2017 was conducted with five land use and land cover classes: EF, MDF, DDF, other vegetation and non-forest. The total areas of classified deciduous forest were 878,831 ha (70.4% of total areas), where MDF and DDF occupied 45.4% and 24.9%, respectively. Areas of other land use types including evergreen forest, other vegetation and non-forest were 370,428 ha (29.7% of total areas). The machine learning RF classifier can be implemented effectively for forest classification in this study site. The overall accuracy of the classified map was 87.8%.

3.3. Variations in NDVI during 2001–2016

3.3.1. Relationship between Satellite-Based NDVI and Observed LAI of Teak Forest Plantation

In this study, satellite-based NDVI was compared against the observed LAI in a local teak forest plantation. NDVI showed a strong exponential relationship with the observed LAI (Figure 7). During the dry season, LAI decreased to the minimum value of about zero, while minimum NDVI remained around 0.3–0.4. However, NDVI seemed to be more sensitive than LAI during the onset of growing season, as NDVI increased more rapidly and reached the maximum value earlier than LAI (Figure 7a). The relationship between LAI and NDVI can be expressed as $\mathrm{LAI}=0.013{\mathrm{e}}^{6.53\mathrm{NDVI}}$, R² = 0.80. For this deciduous teak forest, the NDVI becomes saturated at LAI around ≥ 2 m²/m² (Figure 7b). Potithep et al. [49] also found an exponential relationship between NDVI and in situ LAI in a deciduous forest and NDVI was saturated at values over 0.8. Thus, it is confirmed that the pre-processing NDVI data set has a significant relationship to LAI and its capability for capturing the pattern of seasonality tropical deciduous forest in this study.

Figure 8 shows the relationship between $\mathrm{SO}{\mathrm{S}}_{\mathrm{LAI}}$ and $\mathrm{SO}{\mathrm{S}}_{\mathrm{NDVI}}$ at same thresholds-based (20%). Overall, the SOS derived from MODIS NDVI agreed well with that derived from LAI observation data. However, relationship between $\mathrm{LO}{\mathrm{S}}_{\mathrm{LAI}}$ and $\mathrm{LO}{\mathrm{S}}_{\mathrm{NDVI}}$was rather weak. On average, $\mathrm{SO}{\mathrm{S}}_{\mathrm{NDVI}}$ is 18 days earlier than $\mathrm{SO}{\mathrm{S}}_{\mathrm{LAI}}$ and $\mathrm{LO}{\mathrm{S}}_{\mathrm{NDVI}}$ is 30 days longer than $\mathrm{LO}{\mathrm{S}}_{\mathrm{LAI}}$, respectively. It is noted that $\mathrm{SO}{\mathrm{S}}_{\mathrm{LAI}}$ and $\mathrm{SO}{\mathrm{S}}_{\mathrm{NDVI}}$ were delayed about 4 days and 22 days, respectively in El Niño year (2010).

3.3.2. Temporal Variations of NDVI in Lampang Province during 2001–2016

In Lampang the rainy season starts around March to April, the same time as the deciduous forest starts to bud. A gradual increase in NDVI is also observed during this period (65–100 DOY). Rainy season ends around October to November, corresponding to the decrease of NDVI values during this period (Figure 9). The timing of the minimum NDVI value varied among different forest types. The minimum NDVI value of DDF occurred in March/April and it is lower than the minimum NDVI value of MDF in general. On the other hand, the maximum NDVI value of DDF occurred in June/July and it is also lower than the maximum NDVI of MDF. The seasonal pattern for growing season each year cycle can be described as follows:

NDVI was generally at its lowest (0.3–0.4) during the hottest period of the year (February to April), which is comparable to the field observations of LAI. NDVI then increased rapidly for 4–5 weeks, overlapping with the period of decreasing in temperature and increasing in precipitation resulting in leaf expansion. NDVI reached their maximum peak (0.8–0.9) in July–August, corresponding to maximum and saturation LAI value. NDVI then gradually decreased to below 0.5 and LAI below 2.0 in November, overlapping with the period with lower rainfall and cooler temperature.

3.4. Variations of Phenological Metrics and the Effects of ENSO

3.4.1. Temporal Variations

The temporal variations in tropical deciduous forest phenology metrics were examined for both the dry dipterocarp (DDF) and the mixed deciduous forests (MDF) pixels which were confirmed from field survey locations. The results are shown in Figure 10.

The average start of growing season for dry dipterocarp (${\mathrm{SOS}}_{\mathrm{DDF}}$) and start growing season for mixed deciduous forests (${\mathrm{SOS}}_{\mathrm{MDF}}$) during 2001–2015 were at DOY 113 ± 16, and 113 ± 15, respectively. On the other hand, length growing season of dry dipterocarp (${\mathrm{LOS}}_{\mathrm{DDF}}$) was about 302 ± 18 days, slightly shorter than the ${\mathrm{LOS}}_{\mathrm{MDF}}$ of 309 ± 17 days. During the El Niño (Year 2010), ${\mathrm{SOS}}_{\mathrm{DDF}}$ and ${\mathrm{SOS}}_{\mathrm{MDF}}$was delayed about 26 days and 23 days, respectively and ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$was 19 days and 30 days shorter than the neutral year, respectively. On the other hand, during La Niña (Year 2011) ${\mathrm{SOS}}_{\mathrm{DDF}}$and ${\mathrm{SOS}}_{\mathrm{MDF}}$was 23 and 18 days advanced, whereas ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$was 24 and 13 days, respectively, longer than in the neutral years. An example of the difference in phenological metrics between El Niño (Year 2010) and La Niña (Year 2011) to average of phenology metrics during 2001–2015 is illustrated in Figure 11. It is obvious that ENSO significantly affects the start of the season while the effects of ENSO on the end of season was not clear as that of the start of the season.

3.4.2. Spatial Variations

In order to investigate more details of phenological metrics for two different deciduous forest types (DDF and MDF), we created the mask for different forest types, then phenological metrics for each forest types were extracted and analyzed separately (Figure 12 and Figure 13). We can see that the earliest ${\mathrm{SOS}}_{\mathrm{DDF}}$was in March (DOY 80–90) which mainly occurred in southern part of Lampang. The latest ${\mathrm{SOS}}_{\mathrm{DDF}}$mainly appeared in May (DOY 130–140). This mainly occurred in the northern part of Lampang which may be caused by differences in precipitation and/or altitude (Figure 12). On the other hand, the average ${\mathrm{SOS}}_{\mathrm{MDF}}$ranged between DOY 80 and DOY 140. Similar to ${\mathrm{SOS}}_{\mathrm{DDF}}$_, ${\mathrm{SOS}}_{\mathrm{MDF}}$occurred earlier in the lower altitude of the southern part compared to northern parts of Lampang (Figure 12).

For LOS, ${\mathrm{LOS}}_{\mathrm{DDF}}$mainly ranged from 280 to 330 days while ${\mathrm{LOS}}_{\mathrm{MDF}}$from 290 to 340 days. There is no obvious difference in ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$between northern and southern parts as those found in SOS (Figure 13). There is a slight difference of phenology (SOS and LOS) between DDF and MDF, and their responses to ENSO year is analyzed below. The average ${\mathrm{SOS}}_{\mathrm{DDF}}$in neutral years during 2001–2015 occurred at DOY 106.4 ± 7, similarily to ${\mathrm{SOS}}_{\mathrm{MDF}}$at DOY 106 ± 7. On the other hand, the average of ${\mathrm{LOS}}_{\mathrm{DDF}}$in neutral years was about 307 ± 14 days, slightly shorter than the ${\mathrm{LOS}}_{\mathrm{MDF}}$with 319 ± 13 days (

Table 2. Results of phenological metrics difference during extreme climate events across two forest types in Lampang, Northern of Thailand. The percentage showing the pixels which significant different in temporal between average neutral years and ENSO year.

	Neutral Year (Mean ± 1SD)	El Niño Year (2010)	La Niña Year (2011)	The Area Size Delayed by El Niño (2010) (±Ha, 1SD)	The Areasize Advanced by La Niña (2011) (±Ha, 1SD)
${\mathrm{SOS}}_{\mathrm{DDF}}$ (DOY)	106.4 ± 7.1	131.8	86.5	168,337.5 (75.3%)	147,718.8 (66.0%)
${\mathrm{SOS}}_{\mathrm{MDF}}$(DOY)	106.9 ± 6.9	126.9	92.4	210,662.5 (65.8%)	153,762.5 (48.0%)
${\mathrm{LOS}}_{\mathrm{DDF}}$ (Day)	306.9 ± 13.9	279.0	321.9	214,806.3 (96.0%)	221,818.8 (99.2%)
${\mathrm{LOS}}_{\mathrm{MDF}}$ (Day)	319.0 ± 12.5	292.7	327.3	307,262.5 (96.0%)	307,931.3 (96.2%)

).

Overall, the growing season (${\mathrm{SOS}}_{\mathrm{DDF}}$) started on the 106 ± 7 DOY during the neutral years, delayed to 132 DOY during El Niño (Year 2010) and advanced to 87 DOY in La Niña year (Year 2011). It was delayed about 19 to 33 days in El Niño and was earlier about 13 to 27 days in La Niña (1SD). ${\mathrm{SOS}}_{\mathrm{MDF}}$started almost same time as of ${\mathrm{SOS}}_{\mathrm{DDF}}$ on the 107 ± 7 DOY during the neutral years but was delayed to 127 DOY and advanced to 92 DOY during El Niño and La Niña years, respectively. The ${\mathrm{SOS}}_{\mathrm{MDF}}$ was delayed by about 12 to 28 days in El Niño and was earlier by about 8 to 22 days in La Niña (1SD).

The LOS of these tropical deciduous forests, as indicated by the presence of live canopy, was significantly shorter in that El Niño year compared to neutral years in this study. The results from ANOVA test confirmed that the difference in SOS or LOS between ENSO and neutral years is statistically significant (p ≤ 0.05). During the El Niño (Year 2010), ${\mathrm{SOS}}_{\mathrm{DDF}}$ and ${\mathrm{SOS}}_{\mathrm{MDF}}$was delayed about 35 days and 29 days, respectively and ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$was 31 days and 29 days shorter than the neutral years, respectively. On the other hand, during La Niña (Year 2011) ${\mathrm{SOS}}_{\mathrm{DDF}}$and ${\mathrm{SOS}}_{\mathrm{MDF}}$was 24 and 23 days advanced, whereas ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$was 15 and 9 days, respectively, longer than in the neutral years (p ≤ 0.05). However, we acknowledge that the large sample sizes may result in significant test even if the difference is small and biologically insignificant.

The analysis presented here showed a significant impact of extreme climate events (El Niño and La Niña years) on the timing of leaf flush and the length of growing season in two tropical deciduous forest types in northern Thailand. Moreover, the degree and nature of impact was different between the forest types. In this study, during the El Niño year, the anomaly of maximum temperature significantly increased and summer monsoon rainfall significantly decreased (Figure 5). Suepa et al. [43] also indicated that at there is significant reduction of vegetation growth during EL Niño year. The mechanisms controlling leaf phenology in dry tropical forest are related to water limitation as opposed to light limitation in tropical rainforests [9]. A previous study at the same location [47] found that the interannual variation of the timing of leaf flush responded to increasing soil moisture during March to May. Recent studies of climate impacts on leaf phenology in deciduous species within Asian monsoon and neotropical climate regions observed that differences in rehydration processes following the dry season associated with differences in micro-sites water availability were key in determining the primary factors controlling leaf flush and its response to rainfall [50,51]. In our study, higher maximum temperature and lower precipitation during summer (March, April) likely result in higher water stress during the dry season, which leads to a significant delay in start of growing season during El Niño. In contrast, the higher precipitation and lower maximum temperature in dry season leads to advanced SOS during La Niña year.

The response of tropical deciduous forest was also found in the other El Niño years (e.g., in 2004, Figure 5) where anomaly temperature was higher (+0.37 °C) and precipitation was much lower (−373.3 mm) than neutral years. The average of ${\mathrm{SOS}}_{\mathrm{DDF}}$and ${\mathrm{SOS}}_{\mathrm{MDF}}$in 2004 were 122 ± 12 DOY and 120 ± 13 DOY, a delay of 16 and 13 days, respectively. The average of ${\mathrm{LOS}}_{\mathrm{DDF}}$ and ${\mathrm{LOS}}_{\mathrm{MDF}}$were 277 ± 18 days and 291 ± 18 days, 30 and 28 days respectively shorter compared to the neutral years. This result confirms the effect of El Niño to shifting of tropical deciduous forest phenology metrics. From the analysis we did not find a significant difference in SOS and LOS between DDF and MDF in El Niño year. However, LOS of DDF was significantly longer (15 days) than that LOS of MDF (8 days) in La Niña years (Table 2).

In other years, La Niña (e.g., 2006) where anomaly temperature was close to zero and precipitation was much higher than the average, the variations of tropical deciduous forest phenology metrics were small. It could be explained by very small variation of maximum temperature anomaly and precipitation anomaly in March and April of year 2006. Differences in phenological response to ENSO events between the two forest types were likely related to variations in soil moisture and elevation. The deciduous dipterocarp forests and the teak plantation area are generally located in lower elevation areas ranging from 300 m to 800 m, while mixed deciduous forest extends mostly from 800 m to 1000 m [20]. At the lower elevation, deciduous dipterocarp forest canopy is more open and soils tend to be drier with thinner organic layers and coarser texture relative to the mixed deciduous forest [20]. This is consistent with lower NDVI values during the dry season in the DDF forest type relative to the MDF (Figure 9). It is noted that this study does not consider other related factors that may affect phenology as solar radiation, species diversity, soil layer and soil moisture, which should be addressed in follow-up studies. Further investigation could reveal how the micro-climate, forest species, and soil moisture differ along elevation gradients. Long-term field measurements at different topographic locations would help to enhance our understanding of micro-climatic impact on this forest ecosystem.

4. Conclusions

It is important to gain a better understanding of the phenological responses of tropical deciduous forests to extreme climate events as it will help inform the development of adaptation strategies to reduce the risk declining forest health associated with future climate change. In this study, various data types including meteorological observation, remote sensing, forest inventory maps, digital elevation models, and in situ LAI observation were applied to evaluate the response of tropical deciduous forest to extreme climate anomalies including ENSO events.

The results indicate that phenological metrics of tropical deciduous forest (including dry dipterocarp forest and mixed deciduous forest) varied in response to precipitation and temperature anomalies associated with ENSO events. During the El Niño, significant delay of ${\mathrm{SOS}}_{\mathrm{DDF}}$and ${\mathrm{SOS}}_{\mathrm{MDF}}$ (20–26 days) and shorter of ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$ (26–28 days) were found in most of areas whereas the larger difference was found in DDF. During the La Niña, significant advancement of ${\mathrm{SOS}}_{\mathrm{DDF}}$and ${\mathrm{SOS}}_{\mathrm{MDF}}$ (15–20 days) and longer of ${\mathrm{LOS}}_{\mathrm{DDF}}$and ${\mathrm{LOS}}_{\mathrm{MDF}}$ (8–15 days) were found in most of areas whereas the larger was difference found in DDF. Since frequency and intensity of extreme climate events including ENSO are predicted to increase in the future, tropical deciduous forests may become increasingly vulnerable. Furthermore, climatic impacts on the distribution and productivity of this large region could have significant feedback on global climate systems.