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

: 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 ﬂuctuation, 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). 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.


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 2 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 2 uptake [14]. Such reduced CO 2 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.

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)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(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 Forests 2018, 9,448 3 of 20 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.

of 20
plantation (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(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.

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)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(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].

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)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(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].

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.

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].

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.

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.

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 (

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]. 6 of 20 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.

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 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.

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.
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.

2.
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.

3.
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.

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.

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 SOS NDVI , LOS NDVI and SOS LAI , LOS 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.

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).

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. 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).

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.  1  173  345  152  324  131  303  110  282  88  260  67  239  46  218  25  197  3  175  347  154  326  133  305  112  284  90  262  69  241  51  223  30    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. 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)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(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.

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%.

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%.

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 LAI 0.013e .
, R 2 = 0.80. For this deciduous teak forest, the NDVI becomes saturated at LAI around ≥ 2 m 2 /m 2 (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.

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 LAI = 0.013e 6.53NDVI , R 2 = 0.80. For this deciduous teak forest, the NDVI becomes saturated at LAI around ≥ 2 m 2 /m 2 (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 and at same thresholds-based (20%). Overall, the SOS derived from MODIS NDVI agreed well with that derived from LAI observation data. However, relationship between and was rather weak. On average, is 18 days earlier than and is 30 days longer than , respectively. It is noted that and were delayed about 4 days and 22 days, respectively in El Niño year (2010).

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  Figure 8 shows the relationship between SOS LAI and SOS 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 LOS LAI and LOS NDVI was rather weak. On average, SOS NDVI is 18 days earlier than SOS LAI and LOS NDVI is 30 days longer than LOS LAI , respectively. It is noted that SOS LAI and SOS NDVI were delayed about 4 days and 22 days, respectively in El Niño year (2010).  Figure 8 shows the relationship between and at same thresholds-based (20%). Overall, the SOS derived from MODIS NDVI agreed well with that derived from LAI observation data. However, relationship between and was rather weak. On average, is 18 days earlier than and is 30 days longer than , respectively. It is noted that and were delayed about 4 days and 22 days, respectively in El Niño year (2010).

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

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: 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.

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. 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.

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.
12 of 20 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.

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.     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. 13 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.

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 (Figures 12 and 13). We can see that the earliest SOS was in March (DOY 80-90) which mainly occurred in southern part of Lampang. The latest SOS 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 SOS ranged between DOY 80 and DOY 140. Similar to SOS , SOS occurred earlier in the lower altitude of the southern part compared to northern parts of Lampang ( Figure 12).
For LOS, LOS mainly ranged from 280 to 330 days while LOS from 290 to 340 days. There is no obvious difference in LOS and LOS 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 SOS in neutral years during 2001-2015 occurred at DOY 106.4 ± 7, similarily to SOS at DOY 106 ± 7. On the other hand, the average of LOS in neutral years was about 307 ± 14 days, slightly shorter than the LOS with 319 ± 13 days ( Table 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 (Figures 12 and 13). We can see that the earliest SOS DDF was in March (DOY 80-90) which mainly occurred in southern part of Lampang. The latest SOS 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 SOS MDF ranged between DOY 80 and DOY 140. Similar to SOS DDF, SOS MDF occurred earlier in the lower altitude of the southern part compared to northern parts of Lampang ( Figure 12).
For LOS, LOS DDF mainly ranged from 280 to 330 days while LOS MDF from 290 to 340 days. There is no obvious difference in LOS DDF and LOS 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 SOS DDF in neutral years during 2001-2015 occurred at DOY 106.4 ± 7, similarily to SOS MDF at DOY 106 ± 7. On the other hand, the average of LOS DDF in neutral years was about 307 ± 14 days, slightly shorter than the LOS MDF with 319 ± 13 days ( Table 2). (c) (d) Figure 12. The spatial difference of SOS and SOS between El Niño (a,c) and La Niña (b,d) respectively, compared to that of the neutral year. The significant difference is based on 1SD. Figure 12. The spatial difference of SOS DDF and SOS MDF between El Niño (a,c) and La Niña (b,d) respectively, compared to that of the neutral year. The significant difference is based on 1SD. (c) (d) Figure 13. The spatial different of LOS and LOS between neutral year to El Niño (a,c), La Niña (b,d) respectively. The significant difference is based on 1SD. Figure 13. The spatial different of LOS DDF and LOS MDF between neutral year to El Niño (a,c), La Niña (b,d) respectively. The significant difference is based on 1SD. 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 SOS DDF and SOS MDF in 2004 were 122 ± 12 DOY and 120 ± 13 DOY, a delay of 16 and 13 days, respectively. The average of LOS DDF and LOS 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.

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 SOS DDF and SOS MDF (20-26 days) and shorter of LOS DDF and LOS 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 SOS DDF and SOS MDF (15-20 days) and longer of LOS DDF and LOS 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.