Vertical Differences in the Long-Term Trends and Breakpoints of NDVI and Climate Factors in Taiwan

: This study explored the long-term trends and breakpoints of vegetation, rainfall, and temperature in Taiwan from overall and regional perspectives in terms of vertical differences from 1982 to 2012. With time-series Advanced Very-High-Resolution Radiometer (AVHRR) normalized difference vegetation index (NDVI) data and Taiwan Climate Change Estimate and Information Platform (TCCIP) gridded monthly climatic data, their vertical dynamics were investigated by employing the Breaks for Additive Seasonal and Trend (BFAST) algorithm, Pearson’s correlation analysis, and the Durbin–Watson test. The vertical differences in NDVI values presented three breakpoints and a consistent trend from positive (1982 to 1989) to negative at varied rates, and then gradually increased after 2000. In addition, a positive rainfall trend was discovered. Average and maximum temperature had similar increasing trends, while minimum temperature showed variations, especially at higher altitudes. In terms of regional variations, the vegetation growth was stable in the north but worse in the central region. Higher elevations revealed larger variations in the NDVI and temperature datasets. NDVI, along with average and minimum temperature, showed their largest changes earlier in higher altitude areas. Speciﬁcally, the increasing minimum temperature direction was more prominent in the mid-to-high-altitude areas in the eastern and central regions. Seasonal variations were observed for each region. The difference between the dry and wet seasons is becoming larger, with the smallest difference in the northern region and the largest difference in the southern region. Taiwan’s NDVI and climatic factors have a signiﬁcant negative correlation ( p < 0.05), but the maximum and minimum temperatures have signiﬁcant positive effects at low altitudes below 500 m. The northern and central regions reveal similar responses, while the south and east display different feedbacks. The results illuminate climate change evidence from assessment of the long-term dynamics of vegetation and climatic factors, providing valuable references for establishing correspondent climate-adaptive strategies in Taiwan.


Introduction
Vegetation is an important hub linking the atmosphere, hydrology, and soil; thus, continuous monitoring of vegetation dynamics is imperative in the context of global climate change [1][2][3][4][5][6][7]. As vegetation is sensitive to climate change, and mountainous regions are experiencing more pronounced climate change impacts, many studies have investigated the relationship between vegetation and climatic factors from a vertical perspective [8][9][10][11][12][13]. Located in Southeast Asia, Taiwan is a mountainous subtropical island, with over 60% of the land covered by mountains. The mountains of Taiwan play a critical role in its ecoenvironmental functions of ecosystem conservation, hydrological regulation, and water resource balance for both mountain communities and lowland residents [14,15]. Therefore, many studies have attempted to investigate Taiwan's vegetation dynamics and responses studied vegetation drought in Somalia and the Amazon rainforest, respectively. Chen et al. (2014) [19] applied the BFAST algorithm to detect changes in the Poyang Lake wetland ecosystem in the subtropical monsoon region. Watts and Laffan (2014) [44] evaluated the use of BFAST for detecting fires and floods in semi-arid climates; they found that the breakpoints detected by BFAST corresponded to known flood events, with an accuracy of between 68% and 79%.
From a vertical perspective, this study explores the long-term trends and breakpoints of vegetation and climatic factors-particularly rainfall and temperature-in Taiwan. The vertical perspective is defined as the elevation gradients in the vertical dimension of the terrain. Using AVHRR GIMMS3g NDVI data and gridded temperature and rainfall data produced by the Taiwan Climate Change Estimate Information Project (TCCIP), the vertical differences in the impacts of climatic factors on vegetation-related temporal and spatial scales in Taiwan were studied. The purposes of this study are as follows: (1) to understand the trends and breakpoints in the overall NDVI, rainfall, and temperature in Taiwan; (2) to explore the trends and breakpoints in NDVI, rainfall, and temperature at varying elevations; (3) to investigate the regional characteristics of NDVI, rainfall, and temperature; (4) to study the seasonal variations in NDVI, rainfall, and temperature; and (5) to investigate the relationship between NDVI and climatic factors on a regional scale at varying elevations.

Study Area
Taiwan is located at the junction of the subtropical Eurasian and the Philippine Sea plates. With a total area of 36,000 km 2 , and an elevation range of nearly 4000 m (Figure 1), Taiwan has five mountain ranges running north-south, with more than 200 peaks over 3000 m a.s.l. Additionally, the weather in Taiwan is influenced by the northeast and southwest monsoon systems, bringing abundant rainfall in the rainy season and frequent tropical cyclones in summer. Under such special topography and geological conditions, Taiwan is characterized by rich biological resources and ecosystems ranging from tropical to temperate zones, making the study area the only region with a lush forest ecosystem among countries that pass through the Tropic of Cancer. For the purposes of this study, Taiwan was divided into four regions. Each region's counties and cities were consistent with those from the Taiwan Climate Change Estimate Information Project (TCCIP) ( Table 1).

Research Flow
The research flow of this study is illustrated in Figure 2. First, after collection of data on climatic factors (average temperature, minimum temperature, maximum temperature, and rainfall) and NDVI, all time-series data were divided by elevations at 500 m intervals, based on the DTM (digital terrain model) data derived from the Taiwan Ministry of the Interior across the study period 1982-2012. Later, Pearson's correlation analysis, the Durbin-Watson test, and BFAST analysis were applied to those time-series data to investigate the correlation, autocorrelation, trends, and breakpoints on multiple spatial and temporal scales. Eventually, the seasonal and spatial variations were assessed in order to better understand vegetation responses to climate factors in Taiwan.

The Normalized Difference Vegetation Index (NDVI)
The time-series vegetation data from 1982 to 2012 were derived from the AVHRR GIMMS3g NDVI dataset-an 8 km spatial resolution satellite dataset distributed via th Global Land Cover Facility (http://www.landcover.org/, accessed on 1 October 2017) Monthly NDVI values were obtained via the maximum value composite (MVC) method in order to reduce the interference of clouds and the atmosphere [45]. NDVI has been widely utilized as an important indicator of vegetation conditions based on the red and near-infrared spectra [46,47]. To date, the AVHRR GIMMS3g NDVI dataset is the longes global NDVI time series, and has been carefully calibrated and widely used in large-scal

The Normalized Difference Vegetation Index (NDVI)
The time-series vegetation data from 1982 to 2012 were derived from the AVHRR GIMMS3g NDVI dataset-an 8 km spatial resolution satellite dataset distributed via the Global Land Cover Facility (http://www.landcover.org/, accessed on 1 October 2017).
Monthly NDVI values were obtained via the maximum value composite (MVC) method Remote Sens. 2021, 13, 4707 5 of 29 in order to reduce the interference of clouds and the atmosphere [45]. NDVI has been widely utilized as an important indicator of vegetation conditions based on the red and near-infrared spectra [46,47]. To date, the AVHRR GIMMS3g NDVI dataset is the longest global NDVI time series, and has been carefully calibrated and widely used in largescale phenological studies worldwide [39,46,[48][49][50][51]. The NDVI can be derived using Equation (1) [52]: where RED and NIR represent surface reflectance averaged over ranges of wavelengths in the spectrum's visible and near-infrared regions, respectively. The corresponding time series of monthly temperature and rainfall data at 5 km spatial resolution were obtained from the TCCIP's gridded dataset (https://tccip.ncdr.nat.gov.tw/, accessed on 1 October 2017). The Ministry of the Interior provided the DTM data via an open data platform (https://data.gov.tw/datasets/, accessed on 15 October 2017). The DTM data were used to produce the 500 m interval mask for extracting values from time-series NDVI and climatic data.

Break for Additive Seasonal and Trend (BFAST)
The Break for Additive Seasonal and Trend (BFAST) method was developed by Verbesselt et al. in 2010 [34], and has been successfully applied in many ecosystem-related studies on topics such as such as drought [38], hydrology [53], and climatology [54].
Specifically, the BFAST method decomposes time-series data into the trend, seasonal, and remainder components using Equation (2): where Y t represents the observed data at time t, which signifies the monthly NDVI value, rainfall, and temperature data in this study; T t is the trend component (Trend) and S t is the seasonal component (Seasonal), while e t is the remainder component (Remainder) to capture the remaining variation in the data beyond that in the seasonal and trend components. A piecewise linear relationship is assumed for T t , with sector-specific slopes and intercepts on n + 1 different sectors. Accordingly, there are m breakpoints t 1 * , . . . , t n * , defining t 0 * = 0, so that: for t * j−1 < t ≤ t * j , and where j = 1, . . . , n. The intercept α j and the slope of consecutive linear models β j can be used to determine the magnitude and direction of the abrupt change, and the slope of the gradual change, between the detected breakpoints. The magnitude of an abrupt change at a breakpoint is derived from the difference between T t at t * j−1 and t * j , so that: The magnitude of an abrupt where β j−1 represents the slope before the gradual change, whereas β j represents the slope after the gradual change.
The seasonal component across breakpoints can vary, but is rigid between breakpoints. Let the seasonal breakpoint be given by t 1 # , . . . , t p # , and define t 0 # = 0; then, for t # j−1 < t ≤ t # j , we assume that: where s is the number of observations per year referred to the period of seasonality, and γ i,j represents the effect of season i. Thus, the summation of the seasonal component S t across s successive times is exactly zero for t # j−1 < t ≤ t # j . The seasonal term can be rewritten as: where d i,j = 1 when t is in season i, and 0 otherwise. Consequently, if t is in season 0, The parameter d i,j is often called a seasonal dummy variable, used to account for the seasons in a regression model. Following a designed iteration procedure, the detection of breakpoints-including their number and their position in the time series-is based on the Bayesian information criterion (BIC) and ordinary least squares residual-based moving sum (MOSUM) method, respectively [34,42]. In this study, the BFAST method was implemented using the R environment's bfast package (version 3.5.2, http://www.R-project.org, accessed on 1 January 2018).

Pearson's Product-Moment Correlation Coefficient and Serial Autocorrelation
Pearson's product-moment correlation coefficient is a global-adapted research method developed by Pearson in 1896 [55]. The main concept of the correlation is to evaluate the linear relationship between two variables. Based on the examination of the correlation coefficient (r), the correlation between two variables-and the correlation's direction-can be revealed. Furthermore, the Durbin-Watson test can test the common issue of temporal serial autocorrelation in time-series data [56]. In this study, Pearson's product-moment correlation coefficient and the Durbin-Watson test were implemented in the R environment (Version 3.5.2, http://www.R-project.org, accessed on 1 January 2018).

Temporal Variation of NDVI, Rainfall, and Temperature in Taiwan
In this study, the BFAST method was used to analyze the long-term trends of Taiwan's overall NDVI, rainfall, and temperature over the 30 years from 1982 to 2012, on a monthly basis. A total of 372 observations were derived. For the BFAST method, a minimal section size, h, is employed to distinguish between potentially detected breaks in the trend model given as a portion relative to the sample size. On the other hand, the h parameter determines the sensitivity of BFAST. The value of h also determines the maximal number of breaks to be calculated. The number of samples needs to be at least 30 in order to be considered an efficient sample size based on general statistical rules. Therefore, considering the minimum number of observations in each section divided by the length of the time series, a series of h values from 1/2 to 1/12 were tested to determine the optimal h value for BFAST. Furthermore, the number of iterations was set to 500, as suggested by Verbesselt et al. (2010) [34].
The BFAST results are shown in Tables 2 and 3, along with Figure 3. Only NDVI and minimum temperature (Min-Temp) had breakpoints detected. When h was set to 1/5-1/12 (except 1/6), the temporal characteristics of NDVI had breakpoints. When h was 1/5, the breakpoints occurred in April 1988, August 1994, and October 2000. The slopes changed by 0.0084, −0.0205, −0.0216, and 0.0004, which revealed a positive increasing trend, two negative decreasing trends, and a relatively small positive trend, respectively. A similar trend was found when h was 1/7-1/12; the breakpoints occurred in May 1989, August 1994, and October 2000, with the corresponding intercepts 0.0065, −0.0357, −0.0221, and 0.00002. For the minimum temperature (Min-Temp), the breakpoints occurred in November 1997, with an increasing trend before (0.0175 for the slope value) and a decreasing trend after (−0.0187 for the slope value). Considering the overall data characteristics, including NDVI and climatic factors, the h value was set to 1/7 afterwards for further analysis.
The NDVI trend is closely linked to the sequence of Taiwan's development. As a developing region, Taiwan has experienced population growth with subsequent urbanization, a decoupling of the demand for agricultural land from population growth, and a transition from forest shrinking to forest expansion. Based on a study conducted by Chen et al. (2019) [57], substantial land cover changes were discovered from the historical reconstruction of Taiwan's land cover maps from 1904 to 2015. Between 1982 and 1994, urban expansion consumed agricultural and forested land in the vicinity of larger cities, responsible for the NDVI breakpoints observed in 1989 and 1994. Additionally, according to the Taiwan Forest Management Plan implemented in 1997, the annual logging volume shall not exceed 200,000 cubic meters. Therefore, the forest area of Taiwan has been under control since then, thus reflecting the upward trend in NDVI observed since 2000.
There were no breakpoints detected for the rainfall (Rain), average temperature (Avg-Temp), or maximum temperature (Max-Temp). The slopes were 0.0124 and −0.0068 for rainfall and average temperature, respectively, indicating an increasing rainfall pattern and decreasing average temperature. Although the maximum temperature had no breakpoint, it gradually decreased with a slope of −0.0170.
The increasing rainfall pattern may be associated with the positive trend found in extreme rainfall reported by Henny et al. (2021) [58], based on a 56-year dataset across Taiwan in different seasons and geographic regions; however, the slightly decreasing average temperature and maximum temperature differed from the findings of other studies [59,60]; these discrepancies may result from the different datasets and study periods.

Temporal Variations of NDVI, Rainfall, and Temperature in Taiwan at Varied Elevations
Data on NDVI and climatic factors (average temperature, minimum temperature, maximum temperature, and rainfall) were divided by elevation at 500 m intervals, based on the DTM data from the study period 1982-2012. After applying BFAST, the results shown in Figure 4 displayed a consistent pattern of breakpoints across varied elevations.
A total of three breakpoints were detected from 0 to 3500 m, and two breakpoints from 3500 to 4000 m. For elevations of 0 to 3500 m, the first breakpoint occurred in 1989, the second breakpoint was in 1994, and the third breakpoint was in 2000. The NDVI slope from 1982 to 1989 was positive in terms of trends, reflecting a positive plant growth state; from 1989 to 2000, the slope was negative, reflecting a negative plant growth state. After 2000, only the lowest and second-highest altitudes (0 to 500 m and 3000 to 3500 m) had a positive trend, and the remaining elevations from 500 to 3000 m showed a negative trend. At 3500-4000 m, the two breakpoints were the same as at other altitudes, and the trends were also similar. A total of three breakpoints were detected from 0 to 3500 m, and two breakpoints from 3500 to 4000 m. For elevations of 0 to 3500 m, the first breakpoint occurred in 1989, the second breakpoint was in 1994, and the third breakpoint was in 2000. The NDVI slope from 1982 to 1989 was positive in terms of trends, reflecting a positive plant growth state; from 1989 to 2000, the slope was negative, reflecting a negative plant growth state. After 2000, only the lowest and second-highest altitudes (0 to 500 m and 3000 to 3500 m) had a The NDVI trend also reflects the effects of land use/cover change with associated public administrative policies, such as the Veteran Agricultural Reclamation Policy (VARP) (1955−1988), the Farmland Release Policy (FRP) (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004), and the National Land Restoration Policy (NLRP) (2005−2015) [61]. Based on a study of three major Veterans Affairs Council (VAC) farms located in central Taiwan, with the implementation of the VARP, 21.15% of forestland was transformed into farmland; with the FRP policy, 7.36% of public farmland and 0.72% of private farmland changed into forestland. With the NLRP policy, 27.33% of public farmland and 3.53% of private farmland changed into forestland [61]. According to the high elevation of these three VAC farms (3000 to 3500 m), the positive trend in NDVI corresponded well with the findings of the present study.
Additionally, the positive trend in NDVI found at 0-500 m is consistent with the findings of Abdulmana et al. (2021) [59], who studied NDVI based on satellite images across Taiwan from 2000 to 2020; the authors reported a significant mean increase in decadal trends of NDVI (95% confidence interval: 0.013-0.025 unit/decade). As the effects of land use/cover would interact with the mesoscale meteorological simulation, understanding the exchange mechanism of heat and momentum plays an important role in predictive models [62,63] that can be applied to future climate change mitigation. Figure 5 shows the BFAST results of climatic factors, including rainfall (Rain), average temperature (Avg-Temp), maximum temperature (Max-Temp), and minimum temperature (Min-Temp). There were no breakpoints detected in rainfall, but the rainfall across all of Taiwan was found to be gradually increasing at all altitudes. Furthermore, based on a study of three catchments in northern Taiwan, no significant trend was found from 1980 to 2017 in annual precipitation records derived from five weather stations [64], which is consistent with the absence of breakpoints observed in the present study. Moreover, the reduction in light precipitation (<4 mm h −1 ) and an increase in heavy precipitation (>10 mm h −1 ) derived from the 1961-2005 meteorological data for Taiwan may contribute to the overall increase in rainfall at all altitudes [65].
In terms of average temperature, in contrast to elevations lower than 500 m, the elevations of 500-1000 m and 3500-4000 m gradually decreased. The remaining patterns between 1000 and 3500 m were similar, with one breakpoint occurring in 1988. The trend before 1988 was positive, and the trend after that was negative. The maximum temperature (Max-Temp) for the elevations of 0-500 m gradually increased; the other altitudes showed a cooling trend. The observed increased maximum temperature for the 0-500 m elevations coincided with the nighttime warming trend found in major urban areas [65]. The temperature changes in megacities have been studied worldwide with respect to the effects of urban heat islands and associated health and environmental impacts [60,[66][67][68][69].
There was a breakpoint between 500 and 2500 m in 1986, and the cooling trend after that was relatively gentle. Most of the minimum temperature (Min-Temp) values increased, especially in areas above 2000 m, and the areas below 2000 m only began a downward trend after 1997. There were three breakpoints at the higher altitudes of 3000-3500 m, all of which showed rising trends in 1985, 1988, and 2001. The observed pattern at higher altitudes provides a great insight into the changing climate, because higher altitude areas typically present natural status without many human disturbances. Additionally, montane regions harbor more than half of the world's biodiversity [70]. Therefore, many scholars have explored the effects of climate change on montane vegetation in Taiwan, with valuable results, including the changing phenology [71], the impacts of seasonal drought [72], and the typhoon disturbance mediation of forest structure [73,74].

Regional Variations of NDVI, Rainfall, and Temperature in Taiwan at Varied Elevations
The regional variations of NDVI, rainfall, and temperature in Taiwan at varied elevations were investigated, including the occurrence and the largest changes in the time and direction of breakpoints. Detailed results of regional variations are provided in the supplementary materials ( Figures S1-S5).   Figure 6 illustrates the occurrence of breakpoints examined via BFAST analysis for the NDVI and temperature datasets. There were no breakpoints observed for the rainfall dataset. For NDVI (Figure 6a), the northern region had only one breakpoint at higher elevation (2000-3000 m), while other regions revealed two or more breakpoints over the study periods. The central and eastern regions showed three breakpoints at elevations between 500 and 2000 m, and above 1500 m. The southern region presented two consistent breakpoints across all altitudes.

Occurrence of Breakpoints
The occurrence of NDVI breakpoints may be related to topographic characteristics and complex weather systems such as monsoons, meiyu periods, typhoons, and associated natural disaster events [75,76]. As Taiwan is located on the Pacific Rim of the seismic belt, its unique geological characteristics also make Taiwan vulnerable to natural disasters. For example, based on historical records, an average of 3-4 typhoons make landfall in Taiwan annually, and the associated natural disasters such as flooding, landslides, and debris flows have normally resulted in land deformation at various scales across Taiwan [77][78][79]. Thus, the NDVI breakpoints with the downward trends during the study periods were highly plausible, reflecting island-wide land deformation.
The average temperature (Figure 6b) showed no breakpoints for most of the northern and central regions. One and two breakpoints were observed in the eastern and southern regions, respectively. For the maximum temperature (Figure 6c), areas below 500 m in the northern and central regions, and below 1000 m in the southern and eastern regions, showed one breakpoint. One breakpoint of the minimum temperature was detected for most of the northern regions (Figure 6d). One to two breakpoints were discovered above 1500 m in the central region, and above 500 m in the south and east. Four breakpoints were observed at the highest elevation in the central region. Additionally, the northern region displayed a stable situation, with fewer occurrences, while the central region showed higher variations. Higher elevations showed larger variations in the NDVI and temperature datasets.
The temperature has long been the dominant control factor in global vegetation trends, especially in the higher latitude and altitude areas [6,7,12]. A recent study conducted by Wang and Chang (2021) [80] performed a detailed assessment of the relationships between vegetation, climatic factors, and non-climatic factors at multiple timescales across various land cover types in Taiwan; based on their findings, temperature and precipitation play important roles in vegetation growth; however, the effects vary across seasons. The conclusions of Wang and Chang (2021) [80] correspond to the statement presented by Piao et al. (2014) [81], who recognized that the relationship between vegetation and climate factors might change over time due to the changes in climate and other environmental factors [81]. In general, the eastern region showed its largest change earlier than other regions. In Figure 7b, the observed pattern of average temperature was similar to that of NDVI. The largest change occurred in the late 1980s at an elevation above 1500 m for the central, southern, and eastern regions. For the area above 2000 m in the north, and between 500 and 2500 m in the south, the largest changes occurred in the 1990s. For the maximum temperature (Figure 7c), most eastern regions showed their largest change earlier in 1986, while the other regions showed theirs in the 1990s. Mid-altitude (1000-3000 m) areas showed their largest maximum temperature change in the early 1990s, while lower altitude areas showed theirs in the late 1990s. Lower altitude areas in the southern, central, and eastern regions displayed their largest minimum temperature change in the 1980s (Figure 7d). The largest breakpoint occurred in the majority of the northern regions around the late 1990s. In short, NDVI, average temperature, and minimum temperature showed their largest changes earlier at higher altitude areas, reflecting the sensitivity and diversity of montane regions. Based on a global study of elevational diversity patterns conducted by Guo et al. (2013) [82], 443 montane regions around the Northern and Southern Hemispheres have presented both physical and physiological constraints in various ways. Thus, it is necessary to better quantify and interpret the patterns in order to understand the underlying causal mechanisms. The occurrence of NDVI breakpoints may be related to topographic characteristics and complex weather systems such as monsoons, meiyu periods, typhoons, and associated natural disaster events [75,76]. As Taiwan is located on the Pacific Rim of the seismic belt, its unique geological characteristics also make Taiwan vulnerable to natural disasters. For example, based on historical records, an average of 3-4 typhoons make landfall in Taiwan annually, and the associated natural disasters such as flooding, landslides, and debris flows have normally resulted in land deformation at various scales across Taiwan [77][78][79]. Thus, the NDVI breakpoints with the downward trends during the study periods were highly plausible, reflecting island-wide land deformation.

Largest Changes in the Time and Direction of Breakpoints
The average temperature (Figure 6b) showed no breakpoints for most of the northern and central regions. One and two breakpoints were observed in the eastern and southern   Figure 8 displays the largest change in direction of the NDVI and temperature datasets, which should be interpreted carefully, considering that the largest change in direction was a snapshot of the study period and, thus, could not reveal the whole trajectory of the time series. Nevertheless, based on our results (Figure 8a), a dominant browning of NDVI (direction < 0, shown in red) was observed in the southern and central regions, which is consistent with the long-term landslide monitoring records [83] and the increasing extreme rainfall trend found over the southwestern mountain slopes during the meiyu season [58].
In the eastern and northern regions, NDVI browning was observed at higher altitudes. A dominant decreasing direction of average and maximum temperature was found in most of the eastern and southern regions (Figure 8b,c). An increasing trend (direction > 0, shown in green) was found in the maximum and minimum temperature datasets at higher altitudes. Specifically, the increasing minimum temperature direction was more prominent in the mid-to-high-altitude areas in the eastern and central regions.

ENSO Events
A large-scale climate index, the El Niño-Southern Oscillation (ENSO) is an oceanatmosphere interaction phenomenon that moderates climatic factors with great variations worldwide. Many previous studies have indicated that rainfall variations in northeast Taiwan are influenced by the ENSO [15,84,85]. A recent study further illustrates that ENSO forcing, including monsoon systems, exhibits a coherent association with March and October-November rainfall patterns in Taiwan [86]. Since rainfall and temperature are two critical climatic factors for plants' biophysical conditions, many scholars have devoted research efforts to discerning their relationships. Chang et al. (2013) [17] revealed that ENSO-related rainfall patterns in Taiwan carried out biophysical memory effects on vegetation's time-series patterns. Based on a Holocene stalagmite record, Li et al.  [87] found that during the La Niña years, southern Taiwan experienced decreased rainfall in April-June and increased rainfall in July-September, which were associated with weaker meiyu rainfall and stronger typhoons, respectively. For central Taiwan,   [88] discovered that frequent drought-related wildfires, which impact vegetation in central Taiwan, were associated with the dry and cool weather during the El Niño years.  Figure 8 displays the largest change in direction of the NDVI and temperature datasets, which should be interpreted carefully, considering that the largest change in direction was a snapshot of the study period and, thus, could not reveal the whole trajectory of the time series. Nevertheless, based on our results (Figure 8a), a dominant browning of NDVI (direction < 0, shown in red) was observed in the southern and central regions, which is consistent with the long-term landslide monitoring records [83] and the increasing extreme rainfall trend found over the southwestern mountain slopes during the meiyu season [58].
In the eastern and northern regions, NDVI browning was observed at higher altitudes. A dominant decreasing direction of average and maximum temperature was found in most of the eastern and southern regions (Figure 8b,c). An increasing trend (direction > 0, shown in green) was found in the maximum and minimum temperature datasets at higher altitudes. Specifically, the increasing minimum temperature direction was more prominent in the mid-to-high-altitude areas in the eastern and central regions.

ENSO Events
A large-scale climate index, the El Niño-Southern Oscillation (ENSO) is an oceanatmosphere interaction phenomenon that moderates climatic factors with great variations worldwide. Many previous studies have indicated that rainfall variations in northeast Taiwan are influenced by the ENSO [15,84,85]. A recent study further illustrates that ENSO forcing, including monsoon systems, exhibits a coherent association with March and October-November rainfall patterns in Taiwan [86]. Since rainfall and temperature are two The months of ENSO events determined by the NOAA (National Oceanic and Atmospheric Administration) climate prediction center are marked in orange and green for the El Niño and La Niña events, respectively, in Figures 4 and 5 [89]. Between 1982 and 2012, the La Niña events occurred 10 times, while the El Niño events occurred 9 times; therefore, the total number of months of La Niña events accounted for 30% of the study period from 1982 to 2012 (372 months), while the total number of months of El Niño events accounted for 26% (Table 4). Additionally, between 2000 and 2012, La Niña events occurred five times, while El Niño events occurred four times; however, a larger percentage difference between ENSO events was discovered after 2000. The La Niña and El Niño events accounted for 28% and 20% of the months in 2000-2012 (156 months), respectively. The larger percentage difference implies the increasing frequency of La Niña events, with associated impacts such as weaker meiyu rainfall and stronger typhoons in Taiwan.

Regional and Seasonal Variations of NDVI, Rainfall, and Temperature in Taiwan at Varied Elevations
Four divisions (north, central, south, and east) and four seasons (spring, summer, fall, and winter) were adapted from the TCCIP. Regional and seasonal variations in NDVI, rainfall, and temperature at varied elevations were studied. The seasonal and regional linear trends of NDVI are shown in Figure 9. Most regions except the east showed higher NDVI values in fall/winter and lower values in spring and summer; however, at altitudes below 500 m-especially in the central and southern regions-the discrepancy between seasons was not as clear, associated with a relatively higher level of human occupation. Thus, the area with relatively higher human occupation contains cities and human-controlled vegetation systems such as croplands, agricultural fields, public parks, etc., which are typically maintained at all times. Thus, the variations in NDVI between seasons are due to human influences. In the northern region, an upward trend in NDVI in spring was found at altitudes below 2000 m.
Additionally, NDVI's upward trend was found in summer and winter, while there was a downward trend in autumn. For the central, southern, and eastern regions, there was an upward trend in all seasons. Apart from the regional deviations, it is noticeable that the discrepancies between seasons seemed to gradually shrink, which implies a decreasing strength of seasonality islandwide.
As shown in Figure 10, summer rainfall increased in most regions-especially in the southern and central regions-which is consistent with the findings of increasing extreme rainfall over much of Taiwan during typhoon season [58]. For the areas of the northern region below 1000 m altitude, a negative trend of spring rainfall was observed, requiring further investigation of possible associations with La Niña events and the PDO (Pacific decadal oscillation) [84,90]. The spring and winter rainfall patterns were negative in the higher altitude (1000 to 3500 m) areas, while summer and fall patterns were positive. As an increasing trend of extreme rainfall has been noted in spring and winter by Henny et al. (2021) [58], there is a need to pay attention to the variation between elevations.
Similar rainfall patterns were found for the central and southern regions; again, there were negative trends in spring and winter and positive trends in summer and fall. The positive trends in summer and fall rainfall were in accordance with the consistent increase in extreme rainfall discovered in the typhoon season [58]. Based on another long-term study by Shiu et al. (2009) [65], a reduction in light precipitation (<4 mm h −1 ) and an increase in heavy precipitation (>10 mm h −1 ) may also contribute to the observed rainfall seasonal pattern.
For the eastern region, only areas below 1500 m had a negative spring rainfall trend; however, summer, fall, and winter all showed positive trends. In addition, the highest winter rainfall was observed in the northern region, and the highest summer rainfall was found in the central/southern regions, which experienced northeast and southwest monsoons during winter and summer.
The average temperature reached its peak in summer, and then descended through fall, spring, and winter ( Figure 11). Both the northern and central regions showed an increase in average temperature at altitudes below 500 m for all seasons. At the 500-1000 m altitude, the average winter temperature pattern was positive, while the average spring, summer, and fall temperature patterns were negative. Areas above 1000 m and 1500 m in the northern and central regions, respectively, showed negative patterns for all seasons; however, in the winter, an opposite trend was found at 1000-1500 m. The average temperature trend was negative in the north but positive in the central region. At the same altitude, spring, summer, and fall all showed decreasing average temperature trends. The average temperature in the southern region was found to be positive in all seasons except summer at elevations below 500 m.
Nevertheless, at the elevation of 500-1000 m, winter was the only season showing a positive trend. Decreased average temperatures were found at altitudes of 1000-4000 m. The eastern region's increased average temperature was only observed in spring, fall, and winter at the lower elevations (below 500 m). On the other hand, altitudes above 500 m in the eastern region had decreased average temperature in all seasons. The linear trend of the maximum temperature in different seasons is shown in Figure 12. A negative trend was found in the fall at lower elevations (below 500 m); other than that, the maximum temperature decreased in all seasons at all elevations. In the central region, spring and winter showed increasing maximum temperatures at elevations below 500 m and between 1000 and 2000 m; other than that, decreased maximum temperatures were observed yearround in the central region. However, the southern region had a more complex pattern: at the lower altitudes (below 500 m), the maximum temperature trend was positive in spring, autumn, and winter.
On the other hand, a negative trend took control at altitudes of 500-3000 m for all seasons. At altitudes of 3000-3500 m, trends were positive in spring and summer, but negative in fall and winter. Areas at the highest elevation (3500-4000 m) were all positive. For the eastern region, negative trends were found from 0 to 3500 m; conversely, there was a positive spring trend and a negative trend in the other three seasons at altitudes of 3500-4000 m. As shown in Figure 13, the minimum temperature in the north was positive in all seasons at the altitudes of 0-500 m and 1000-2500 m; however, at elevations of 500-1000 m, the minimum temperature decreased in spring, but increased in fall and winter. Additionally, as the elevation rose to 2500-3500 m, the minimum temperature decreased in spring and winter, but increased in summer and fall. Thus, an overall positive minimum temperature trend was found in the central, southern, and eastern regions.
Overall, large variations in rainfall and temperature were found in different regions, elevations, and seasons. Thus, a comprehensive understanding of the spatial and temporal variations is challenging. Based on many valuable previous works, the association between the ENSO and rainfall in Taiwan is strong and complex [84,91]. As the ENSO affects an abnormal lower tropospheric anticyclonic circulation around the Philippine Sea in the El Niño phase, altering the East Asian subtropical front would influence rainfall in the meiyu season [92]. Jiang et al. (2003) [84] studied 50 years of daily rainfall station data in Taiwan to examine the relationship between spring rainfall and the ENSO; the authors found that a weak mid-latitude frontal system in the East China coastal area and an anomalous anticyclone over the Philippine Sea area may correspond to Taiwan's heavy spring rainfall events. Lin et al. (2015) [91] further investigated two types of ENSO (canonical ENSO and ENSO Modoki) in detail, and evaluated their effects on rainfall; they concluded that during canonical El Niño episodes, the rainfall increases in winter and spring but decreases in summer and fall; however, in El Niño Modoki episodes, the rainfall pattern shows the opposite response. Meanwhile, during two types of La Niña episodes, only fall rainfall increases. The effects of the Central Mountain Range (CMR), monsoons, and topography in Taiwan should also be considered.       Therefore, we adopted the definitions of the dry and wet seasons from two previous studies [93,94] in order to investigate how the dry (spring and winter) and wet (summer and fall) seasons had changed in terms of the amount of rainfall across Taiwan over the three-decade study period. As shown in Figure 14, the differences in the amounts of rainfall between the dry and wet seasons revealed increasing discrepancies over time across Taiwan, indicating that the difference in rainfall between the dry and wet seasons was becoming larger. Furthermore, the largest difference was found in the southern region, while the smallest was found in the north, revealing an unbalanced potential water resource issue in Taiwan.

Correlation and Serial Autocorrelation Analysis of NDVI and Climatic Factors
In this study, Pearson's correlation coefficient was used to analyze the correlat between NDVI and climatic factors at all altitudes in Taiwan; the results are shown Figure 15a. Most of Taiwan's NDVI and climatic factors had a significant negative cor lation (p < 0.05). The only positive relationship was found between maximum and m mum temperatures at low altitudes below 500 m. At altitudes above 3000 m, the negat effect of the maximum temperature on NDVI gradually reduced. For rainfall, insignific relationships were found at altitudes below 500 m.

Correlation and Serial Autocorrelation Analysis of NDVI and Climatic Factors
In this study, Pearson's correlation coefficient was used to analyze the correlation between NDVI and climatic factors at all altitudes in Taiwan; the results are shown in Figure 15a. Most of Taiwan's NDVI and climatic factors had a significant negative correlation (p < 0.05). The only positive relationship was found between maximum and minimum temperatures at low altitudes below 500 m. At altitudes above 3000 m, the negative effect of the maximum temperature on NDVI gradually reduced. For rainfall, insignificant relationships were found at altitudes below 500 m. On the other hand, rainfall slightly affected NDVI in the southern and eastern regions at elevations below 500 m. For the southern region, the significant correlation level of temperatures decreased when the elevation increased. The correlation level of temperature decreased in the mid-latitude area (1000-3000 m). However, the positive effects of temperatures extended to altitudes up to 2500 m. Notably, the eastern region displayed more diverse responses, attributed to its broader spatial coverage across latitudes. However, a potential lag correlation may need further investigation at a later date. Figure 16 illustrates the outcome of the Durbin-Watson test. As shown in Figure 16a, the NDVI residual component revealed no significant autocorrelation. However, for the majority of the temperature datasets, significant positive autocorrelation was observed (Figure 16b-d). On the other hand, the residual rainfall component displayed significant negative autocorrelation across Taiwan (Figure 16e). The presence of significant autocorrelation indicates that the residual components were not independent of one another. However, autocorrelation in the residual components was expected, as the h parameter in the BFAST algorithm restrains the level of detail within the linear segments. Therefore, potential variations may exist over shorter periods than the trend segment length permitted by the h parameter values used [44]. On the other hand, rainfall slightly affected NDVI in the southern and eastern regions at elevations below 500 m. For the southern region, the significant correlation level of temperatures decreased when the elevation increased. The correlation level of temperature decreased in the mid-latitude area (1000-3000 m). However, the positive effects of temperatures extended to altitudes up to 2500 m. Notably, the eastern region displayed more diverse responses, attributed to its broader spatial coverage across latitudes. However, a potential lag correlation may need further investigation at a later date. Figure 16 illustrates the outcome of the Durbin-Watson test. As shown in Figure 16a, the NDVI residual component revealed no significant autocorrelation. However, for the majority of the temperature datasets, significant positive autocorrelation was observed (Figure 16b-d). On the other hand, the residual rainfall component displayed significant negative autocorrelation across Taiwan (Figure 16e). The presence of significant autocorrelation indicates that the residual components were not independent of one another. However, autocorrelation in the residual components was expected, as the h parameter in the BFAST algorithm restrains the level of detail within the linear segments. Therefore, potential variations may exist over shorter periods than the trend segment length permitted by the h parameter values used [44].

Conclusions
This study is the first comprehensive assessment evaluating Taiwan's vegetation and climate factors for associated long-term trends and breakpoints from a vertical perspective. Based on the overall, regional, and seasonal perspectives, our major findings are summarized as follows: • The overall NDVI trend in Taiwan starts with an increasing pattern, and then follows two downward trends until 2000, followed by a gentle increase. Among climate factors, only the minimum temperature has one breakpoint, found in the late 1990s; • Based on the vertical investigation, the trends and breakpoints of NDVI across Taiwan are highly consistent, with three breakpoints for the altitudes below 3500 m (1989, 1994, and 2000) and two breakpoints above 3500 m (1989 and 1994). Rainfall gradually increases, with no breakpoints. Average and maximum temperature have similar trends, while minimum temperature shows variations-especially at higher altitudes.

•
Regarding the occurrence of breakpoints for regional assessment, the northern region displays a stable situation, while the central region shows higher variations. Higher elevations show larger variations in the NDVI and temperature datasets. NDVI, average temperature, and minimum temperature reveal their largest changes earlier in higher altitude areas. Specifically, the increasing minimum temperature direction is more prominent in the mid-to-high-altitude areas of the eastern and central regions;

•
The incidence of La Niña events is larger than that of El Niño events for the whole study period, and especially after 2000; • Seasonal variations are observed for each region. However, the difference between the dry and wet seasons is getting larger, with the smallest difference in the northern

Conclusions
This study is the first comprehensive assessment evaluating Taiwan's vegetation and climate factors for associated long-term trends and breakpoints from a vertical perspective. Based on the overall, regional, and seasonal perspectives, our major findings are summarized as follows:

•
The overall NDVI trend in Taiwan starts with an increasing pattern, and then follows two downward trends until 2000, followed by a gentle increase. Among climate factors, only the minimum temperature has one breakpoint, found in the late 1990s; • Based on the vertical investigation, the trends and breakpoints of NDVI across Taiwan are highly consistent, with three breakpoints for the altitudes below 3500 m (1989, 1994, and 2000) and two breakpoints above 3500 m (1989 and 1994). Rainfall gradually increases, with no breakpoints. Average and maximum temperature have similar trends, while minimum temperature shows variations-especially at higher altitudes.

•
Regarding the occurrence of breakpoints for regional assessment, the northern region displays a stable situation, while the central region shows higher variations. Higher elevations show larger variations in the NDVI and temperature datasets. NDVI, average temperature, and minimum temperature reveal their largest changes earlier in higher altitude areas. Specifically, the increasing minimum temperature direction is more prominent in the mid-to-high-altitude areas of the eastern and central regions; • The incidence of La Niña events is larger than that of El Niño events for the whole study period, and especially after 2000; • Seasonal variations are observed for each region. However, the difference between the dry and wet seasons is getting larger, with the smallest difference in the northern region and the largest difference in the southern region; • Taiwan's NDVI and climatic factors have a significant negative correlation (p < 0.05), but the maximum and minimum temperatures have significant positive effects at altitudes below 500 m. The northern and central regions reveal similar responses, while the south and east display different feedbacks.
Our study's results illustrate the vertical differences in climate change evidence in Taiwan, improving our overall understanding of vegetation and climatic dynamics over long periods across different regions and altitudes. Our results have potential applications in supporting ecological management and climate adaptation strategies in Taiwan to meet future climate change challenges.