Assessment of vegetation dynamic and its effects in a large-scale landslide in Central Taiwan with multitemporal Landsat images

1.   Introduction

Earthquake-induced landslides are frequent natural hazards in Taiwan due to its unique geographical location and climatic conditions. One of the most severe landslides was triggered by the 921 Earthquake in 1999, which was classified as a large-scale landslide ^[1,2,3]. In Taiwan, large-scale landslides are defined as those with a collapsed area exceeding 10 hectares, an earth volume greater than 100,000 cubic meters, or a collapse depth of more than 10 meters ^[4]. The exposed soil and rock slopes caused by landslides are prone to sediment-related disasters, especially during heavy rainfall when vegetation cover is insufficient. Therefore, monitoring and understanding vegetation recovery in landslide-affected areas are critical for the effective management and mitigation of landslide hazards ^[5,6,7,8,9]. Furthermore, monitoring vegetation restoration provides valuable insights into the restoration cycle within landslide-affected areas. This information serves as a critical reference for developing and improving future vegetation restoration strategies in such regions. Although traditional vegetation survey methods offer important information on flora, succession pathways, and biodiversity, their application in large-scale landslide areas is challenging. These methods require long-term monitoring and assessment, which are both time-consuming and labor-intensive ^[10,11]. Furthermore, their implementation is often impractical in remote or inaccessible regions.

In response to these challenges, remote sensing technology has been used as an effective alternative for monitoring vegetation changes. The Normalized Difference Vegetation Index (NDVI) is commonly used to track dynamic changes in vegetation ^{[12,13,14,15,16,17,18,19,20]}. Researchers have used multi-temporal NDVI data to detect changes in vegetation cover after landslides and assess the extent of vegetation recovery ^{[21,22,23,24,25,26,27]}. However, these methods are limited to examining past vegetation changes and do not facilitate the prediction or simulation of future vegetation restoration scenarios. To address these limitations, landscape change models have been developed to simulate the functional and dynamic changes in land use systems, which offers a more comprehensive approach to forecasting future vegetation recovery ^[28,29,30]. These models can explore the interactions of natural processes and evaluate proposed management treatments ^[31,32,33]. Furthermore, landscape change models are regarded as an effective tool for post-disaster vegetation restoration and simulating future vegetation succession. However, many models overlook the sequential nature of vegetation succession. This may lead to simulation results that do not accurately reflect real-world conditions. This gap highlights the necessity of incorporating additional succession processes in the modeling approach.

We leveraged multi-temporal Landsat remote sensing imagery, captured before and after the 921 Earthquake in the Chiufenershan landslide area, in combination with a land use change model and the principles of vegetation succession priority. For this integrated approach, we aim to simulate and assess the dynamic restoration of vegetation and potential future landscape changes within the affected area. The findings are expected to contribute to the development of effective landslide management strategies, mitigate future disaster risks, and improve the efficacy of vegetation restoration efforts.

2.   Materials and methods

2.1.   Study area

The study area is in Nangang Village, Nantou County in the central part of Taiwan (R.O.C.). The Chiufenershan landslide occurred between 23°58′08″N and 23°56′52″N and between 120°49′36″E and 120°51′01″E (Figure 1). The altitude of the study area varies from 500 to 1000 m above sea level. The Shizikeng and Jiucaihu Rivers transformed into barrier lakes due to the 921 Earthquake that collapsed the area and blocked both rivers. The collapsed area is 102.5 ha, with the depth collapse ranging from 30 to 50 m. The collapsed volume was 32.85 million m³. The disaster is classified in a large-scale landslide category. The area is divided into three parts based on a top-down view of the location: A large collapsed area, a deposition area, and a conservation park. There is an almost intact semi-natural area in the western part outside the landslide area with some orchards and betel nut plantations. In contrast, the eastern part is steep areas resulting in minimal disturbance with complex categories of vegetation ^[24].

2.2.   Data acquisition

2.2.1.   Multitemporal satellite imagery

Satellite imagery is commonly adopted for environmental monitoring, landscape change, and vegetation restoration assessments due to its wide detection range, fixed period, and multi-temporal-spectral properties. Landsat imagery has been widely used for landscape and vegetation monitoring due to its long-term data availability and moderate spatial resolution ^[34]. Similarly, Sentinel-2 provides high-resolution optical data with a 5-day revisit time, suitable for detailed environmental analyses ^[35]. MODIS (Moderate Resolution Imaging Spectroradiometer), with its daily global coverage and multi-spectral capabilities, supports large-scale monitoring of vegetation dynamics and land-use changes ^[36]. These satellite platforms, among others, enable comprehensive monitoring of environmental changes over time through the integration of multi-temporal and multi-spectral data. We used satellite imageries from Landsat 5 and Landsat 8 due to the long-term availability ^[37], spatial resolution (30 m), and free application ^[38] for vegetation dynamics assessment and landscape change simulation. The Chiufenershan landslide occurred in a mountainous area, where imageries are often affected by clouds and shadows. Therefore, for the period 1990 to 2020, 10 satellite images were selected from April 1, 1999 (pre-earthquake), September 24, 1999 (post-earthquake), 2000, 2002, 2006, 2009, 2013, 2015, 2018, and 2020. The images from 1999 were applied for landslide mapping, and the remaining images were used for land use classification (Figure 2).

2.2.2.   Long-term rainfall data

Rainfall is a critical factor for vegetation growth; however, excessive precipitation can lead to runoff and surface erosion. We utilized rainfall as a key driving factor and considered its impact on spatial distribution. Long-term annual rainfall data from 2000 to 2020 were obtained from nearby weather stations. The data was employed to model the spatial distribution of rainfall. The rainfall data were collected by the Water Resources Department of the Ministry of Economic Affairs and the Central Weather Bureau of the Ministry of Transportation and Communications of Taiwan (Table 1).

2.3.   Analytical methods

2.3.1.   NDVI

Vegetation detection relies on the unique spectral properties of plants, which absorb blue and red light while reflecting near-infrared radiation ^[39,40]. The Normalized Difference Vegetation Index (NDVI), proposed by ^[41], represents this spectral difference by calculating the ratio of the difference between the near-infrared and red bands to their sum ^[42,43]. NDVI is widely utilized in remote sensing for evaluating vegetation restoration, classifying land use, and modeling vegetation changes. The values of NDVI range from −1 to 1, with negative values indicating non-vegetated areas and positive values reflecting varying degrees of vegetation cover. The formula for calculating NDVI is as follows:

where R is the red band and NIR is the near-infrared band.

2.3.2.   Landslide mapping

The NDVI difference values were calculated for images taken before and after the landslide event. The most severely affected area was identified on the map and selected as the seed point to determine the initial threshold for the landslide area. The landslide triggered by the 921 Earthquake was mapped through a comparison of images captured before and after the event ^[44]. The stratified sampling was used to select 250 random samples of collapsed and non-collapsed points from the map to evaluate the assessment accuracy. The accuracy of landslide mapping was evaluated based on the overall accuracy (OA) and the kappa coefficient ^[45,46].

2.3.3.   Land use classification

Image classifier—Support vector machine

In the absence of multi-temporal land use maps for the study area, remote sensing image classification techniques are necessary to generate land use maps and assess land cover changes over time. Support Vector Machine (SVM) is a supervised learning model based on statistical theory. SMV identifies an optimal hyperplane within the input space to separate the best classes in the data. SVM can handle both linearly separable and non-linearly separable data by employing different types of kernels, such as linear, polynomial, or radial basis function kernels, which effectively map the input space into a higher-dimensional feature space where the classes become more easily separable. It is used to map an inseparable sample from a low-dimensional space into a higher-dimensional space where it identifies the optimal straight line, or hyperplane, that separates the sample sets within that space. The optimal hyperplane is defined as the one that maximizes the margin, or the greatest possible distance, between the sample sets, particularly in the context of binary classification (Figure 3).

The sample set x with n records, where ${x}_{i}$ represents the feature vector of the ith record, and the records belong to two categories, ${w}_{1} = 1$ or ${w}_{2} = 1$. The hyperplane formula is as follows:

where w is the normal vector of the hyper-plane and b is the threshold.

Classification is performed based on the distance between the sample and the optimal hyperplane. If $f\left({x}_{i}\right) > 0$, the sample ${x}_{i}$ is classified into ${w}_{1}$. Conversely, if $f\left({x}_{i}\right) < 0$, the sample is classified into ${w}_{2}$. The classification constraints are expressed as follows:

SVM is commonly applied for remote sensing image classification. In this study, the segmentation and classification module of ArcGIS 10.5 (ESRI) was used for land use classification.

Accuracy assessment

Accurate land use classification is crucial for ensuring the reliability of the results. Statistical methods are commonly applied to assess the quality of classification. One widely used approach is the confusion matrix, which evaluates the performance of land use classification ^[48,49]. Among key metrics, the OA and kappa coefficient show the validity of classification results. These metrics are crucial for measuring the comprehensive accuracy between simulated and observed maps ^{[50,51,52,53]}. The kappa coefficient, which ranges from 0 to 1, provides a quantitative measure of classification accuracy. The higher values indicate greater accuracy. In remote sensing applications, accuracy assessments are vital for determining the suitability of classification results for specific purposes. A kappa coefficient greater than 0.7 is generally considered indicative of valid classification results ^[54,55,56]. The OA, calculated from the confusion matrix, represents the proportion of correctly classified sample points, weighted by the number of samples in the i-th row and j-th column. It offers an objective measure of classification performance, with higher values reflecting higher accuracy. The formula for calculating OA is as follows:

where ${\mathrm{X}}_{\mathrm{i}\mathrm{i}}$ represents the number of sample points on the diagonal of the confusion matrix indicating correctly classified instances, and ${\mathrm{X}}_{\mathrm{i}\mathrm{j}}$ refers to the number of sample points in the i-th row and j-th column of the confusion matrix, reflecting misclassifications between categories.

The Kappa coefficient is calculated by comparing the classified data with actual ground truth information. This comparison is used to assess the accuracy of the classification, providing a more reliable measure of agreement between the predicted and observed values. The formula is as follows:

where n represents the number of rows in the confusion matrix and N is the total number of samples. ${\mathrm{X}}_{\mathrm{i}\mathrm{i}}$ is the number of samples on the diagonal of the confusion matrix, ${\mathrm{X}}_{\mathrm{i}+}$ is the number of samples in each row, and ${\mathrm{X}}_{+\mathrm{i}}$ is the number of samples in each column of the confusion matrix.

2.3.4.   Landscape change analysis

Determining driving variables of post-landslide landscape changes

Landscape patterns are shaped by the interactions between natural processes and human activities ^[57,58]. In this study, human activities were not considered effective driving factors since there was no residential or agricultural activity in the study area after the earthquake disaster. Consequently, the landscape changes were primarily caused by environmental stressors. Factors such as altitude, slope, topographic relief, topographic position, and topographic wetness significantly influence temperature, transportation, and distribution of moisture and nutrients across the landscape. Therefore, the distribution of landscape vegetation is strongly related to topographic patterns. In addition, solar radiation is a key factor driving various physical and biological processes on the Earth's surface. It is essential for photosynthetic plant growth and is influenced by topography, surface features, and seasonal variations. Therefore, solar radiation serves as another critical determinant of landscape change. Furthermore, road development disrupts vegetation and contributes to surface runoff, which accelerates soil erosion, impeding plant growth and hindering vegetation succession. Another important factor is the presence of rivers. When a river flows through the surface soil, capillary action enables the soil to retain water, which is beneficial for plant growth. Therefore, we utilized a combination of meteorological data, digital elevation models (DEMs), road maps, and river network maps to derive the driving factors for landscape change simulation using the spatial analysis module in ArcGIS 10.5. Meteorological data, including annual rainfall records, were obtained from the Water Resources Agency, Ministry of Economic Affairs (WRA), and the Central Weather Administration, Ministry of Transportation and Communications (CWA). The DEM, provided by the Ministry of the Interior (MOI), was utilized to calculate terrain-related indices such as slope, topographic relief, topographic position, and topographic wetness. Road and river network maps, sourced from the National Land Surveying and Mapping Center, Ministry of the Interior (NLSC), were used to compute proximity variables, including the distance to roads and rivers (Table 2).

Landscape change prediction

The artificial neural network (ANN) is a nonlinear statistical method used to model relationships between influencing factors and output variables. ANNs can predict or classify variables more accurately than traditional statistical methods by leveraging mathematical and statistical learning techniques. ANNs identify patterns from large data sets more accurately compared to traditional statistical methods. Given the complexity of factors affecting plant growth, identifying correlations between these variables enhances the accuracy of simulations and predictions. We measured the relationship between each variable and each land use category. The selected variables included annual rainfall, solar radiation, altitude, slope, topographic relief index, topographic position index, topographic wetness index, distance to roads, and distance to rivers. Then, the dynamic landscape change simulation was run, with particular attention paid to calibrating the model and setting appropriate hyperparameters to achieve optimal results. Calibration was performed using images from 2003 and 2009, as well as from 2016 and 2018. The trial-and-error method was applied to determine the optimal model.

The Markov chain and cellular automata models are widely used to simulate the spatial distribution of landscape patterns. However, these models often overlook the sequential nature of vegetation succession, which follows a defined progression rather than a random sequence (e.g., from bare land to grassland and then to forest; Figure 4; ^[67]). This sequential progression was carefully considered in this study, which aimed to explore vegetation restoration and dynamically simulate landscape changes across different periods. TerrSet, developed by Clark Labs at Clark University, is an integrated geospatial software system designed for monitoring and modeling Earth systems to support sustainable development. We used the Land Change Model (LCM) within TerrSet to analyze potential future vegetation changes in study areas. The modeling framework combines historical data, spatial variables, and transition potential models to simulate and predict spatial phenomena. Machine learning techniques in TerrSet, specifically ANNs, were employed to create correlations between explanatory factors and observed landscape changes. The input data consisted of two land use maps and key driving factors (Table 2). ANN-based modeling iteratively adjusted weights and biases through backpropagation to minimize error, thereby generating predictions for potential future land use patterns, including forest, grassland, bare land, and water. The optimal ANN parameters used in this study included a maximum number of iterations of 20,000, a learning rate of 0.001, a momentum of 0.5, and hidden layers of 12.

3.   Results

3.1.   Chiufenershan landslide mapping

The landslide triggered by the 921 Earthquake was mapped using the NDVI image subtraction and threshold change approach. The total landslide area was found to be 223.74 ha. OA and Kappa coefficient statistical tests were used to measure the accuracy of the landslide delineation. OA represents the proportion of correctly classified pixels relative to the total number of reference pixels. The Kappa coefficient measures the agreement between the mapped results and reference data, accounting for chance-level agreement. The results showed an OA of 86.4% and a Kappa coefficient of 0.728, as presented in Table 3 and Figure 5. These values indicate good performance of the model based on the classification accuracy standards proposed by ^[54,68]. This suggested that the methodology provided more reliable results for assessing landslide impacts caused by seismic activity. The achieved accuracy highlighted the effectiveness of NDVI-based image subtraction in detecting landslides, particularly in complex mountainous terrain, before and after earthquakes.

3.2.   Accuracy assessment of land use classification by period

After coupling the NDVI data with the original image bands, Support Vector Machine (SVM) classification was applied to categorize land use. The accuracy evaluation results showed OA values exceeding 80% and Kappa coefficients greater than 0.7 for each image (Table 4). These results highlighted high classification accuracy ^[54,55]. Therefore, the classification results were considered reliable for the subsequent dynamic simulation of future landscape changes. Image inspection revealed that most error points were along the edge of the landslide or at the boundaries of land use categories (Figure 6). These errors are likely attributable to the spatial resolution (30 m) and the presence of mixed cells, where a single grid contains multiple land use types. The confusion matrix for each period indicated that most errors occurred between grassland and forest areas. This is primarily because grassland and forest have similar spectral reflectance values, making them difficult to distinguish in classification.

3.3.   Dynamic simulation of landscape change in the future

3.3.1.   Landscape change model calibration and validation

The relationship between each driving factor and landscape category was identified through ANN. The Markov chain and cellular automata models were integrated to complete a dynamic simulation of future landscape changes. ANNs are commonly used to identify relationships between factors. However, optimizing parameters such as the maximum number of iterations, learning rate, momentum, and the number of hidden layers is crucial for achieving acceptable results. The Markov chain model simulates the probability distribution of landscape change for the next period (t+∆t) based on the differences observed between the two periods (∆t). For this study, data from 2003 and 2009, as well as from 2016 and 2018, were used to calibrate and validate the landscape change for the years 2015 and 2020. In addition, the vegetation succession process was analyzed in terms of landscape patterns under restricted conversion conditions. Through the trial-and-error method, the optimal parameters were determined: A maximum of 20,000 iterations, a learning rate of 0.001, a momentum of 0.5, and 12 hidden layers. The evaluation results of the landscape change simulation showed OAs of 72.80% for 2015 and 82.26% for 2020, with Kappa coefficients of 0.53 and 0.72, respectively. These findings indicated medium to high accuracy for both simulations. The results demonstrated that optimized parameters effectively contributed to the dynamic simulation of subsequent landscape changes.

3.3.2.   Landscape change prediction

The optimized parameters derived from the 2015 and 2020 simulations were applied to predict dynamic landscape changes. The area percentages of each land cover category were then calculated based on the dynamic simulation results (Figure 7 and Table 5). The findings suggested that, in the absence of external disturbances, the land cover categories will stabilize by 2075. Forests comprised the largest proportion of land cover, approximately 60%, followed by grasslands at around 27%. Bare land and water bodies were 13% of the total land cover. These results revealed that vegetation in the Chiufenershan landslide area has smoothly recovered, with forests emerging as the dominant land cover type.

4.   Discussion

4.1.   Assessment of long-term vegetative restoration

The transition probability matrix for landscape change from 2000 to 2020 indicated a significant negative net change of −136.89 ha in bare lands. However, forest areas experienced the largest positive net change of 90.18 ha (Table 6). The results showed that the landslide area has largely reverted to forest since 2000, with forest cover at approximately 52.29% of the total landslide area, followed by grassland at 28.80% (Table 7 and Figure 8). The spatial distribution analysis of the landscape pattern within the landslide demonstrated that water areas were predominantly concentrated in the two barrier lakes (Figure 6, Table 7). Following the landslide event, the water areas of these lakes showed only minor changes. Bare land has been steadily decreasing, with the primary fragmentation occurring in the deposition area beneath the landslide. The long-term monitoring and assessment of the Chiufenershan landslide following the earthquake by ^[69] showed that vegetation restoration was closely related to terrain conditions. Their study indicated that the deposited area, characterized by loose soil, provides appropriate conditions for vegetation invasion and the natural renewal of residual vegetation. This deposited area serves as a primary zone for natural vegetation recovery. In contrast, restoring vegetation in the collapsed area of the landslide is challenging due to exposed rock and unsuitable environmental conditions. Furthermore, ^[70] highlighted that the topography of the landslide significantly impacts the growth potential of plants. The collapsed area, which features a steep slope, had low vegetation coverage and was dominated by herbaceous plants. In contrast, the deposited area was covered by more vegetation due to its gradual slope, which is more conducive to the growth of woody plants. Our findings align with these studies, confirming the role of topography in influencing vegetation recovery. Additionally, a comparative analysis across various time points revealed a gradual decline in the rate of forest succession, followed by an increase in grassland succession rate after 2006. This trend suggested that grassland has emerged as the predominant form of vegetation succession during this period. Moreover, an examination of spatial vegetation changes within the landscape indicated vegetation restoration process has progressively shifted from the deposited area to the collapsed area over time.

The dynamic simulation of future landscape changes indicated that the deposited area has essentially stabilized and forest covers become the dominant land cover (Figure 7). However, the collapsed area (upper landslide), characterized by thin soil layers, steep slopes, and ongoing terrain instability, remains less conducive to the establishment of deep-rooted tree species. Therefore, vegetation recovery in the collapsed area is currently dominated by grassland. It is expected by improvements in soil formation and slope stability over time, more favorable conditions would emerge for tree development. From the perspective of landscape fragmentation, forest regeneration is primarily occurring along the edges of fragmented patches, where environmental conditions are more stable. The results showed a decreasing trend in bare land areas, which are in the central area. On the other hand, grassland continues to expand inward from the margins. From a successional perspective, the findings suggested the collapsed area will require a significant period to reach a climax community. The results showed the stability of the deposition area and the ongoing potential for recovery in the collapsed area for current conditions in the study area. Therefore, vegetation restoration efforts between 2025 and 2075 will focus on enhancing recovery in the collapsed area.

4.2.   Vegetation restoration process and effects of the landslide

The simulated spatial distribution showed vegetation succession progresses from the bottom upward (Figure 6 and Figure 9c). This pattern was created throughout the landslide event. The soil and rock fall downward and accumulate, which results in a thick soil layer and a moderate slope. These conditions provide an appropriate environment for plant growth ^[69,71,72]. Consequently, the vegetation recovery rate in this area was relatively rapid and led to the highest levels of forest coverage. In contrast, the upper part of the slope, where the collapse was deeper and bedrock is exposed, provides minimal topsoil for plant growth. Therefore, vegetation recovery in this area was slower and relies on the process of soil genesis. In addition, gramineous plants have been identified as the pioneer species in this region due to the steep and varied environmental conditions (Figure 9a and Figure 9c). The simulation indicated that the right side of the collapsed area, which is closer to the epicenter and characterized by exposed bedrock, will likely remain bare once the landscape pattern stabilizes. Vegetation restoration in this area is expected to take more time due to these challenging conditions.

4.3.   Field investigation and verification of Chiufenershan landslide vegetation

The results of the field investigation revealed that a significant portion of the collapsed area remains bare, with a weak process of vegetation restoration (Figure 9b). This is mainly due to the steep terrain and the vulnerability of soil to erosion. In areas where soil from the upper slopes has been washed away by rainfall and formed a thick layer of soil in the upslope of the roads, vegetation has begun to regenerate and recovery has been facilitated (Figure 9a). However, grassland remained the dominant cover, with Pennisetum purpureum Schumach identified as the pioneer species in this area. Furthermore, there were significant differences in plant distribution at the intersection of the deposited and collapsed sites (Figure 9c). The collapsed area was predominantly covered by grassland, while the deposited area was characterized by secondary forest (Figure 9d). The dominant species in the secondary forest included Machilus zuihoensis Hayata and Schefflera octophylla (Lour.) Harms ^[73]. The findings from the field investigation verified that the vegetation restoration assessment results through the model were aligned with the actual vegetation succession observed in the study area. This suggested that the methodology proposed in this study addressed previously overlooked mechanisms, thereby enhancing the ability of the model to predict vegetation dynamics. As a result, this model offers a more accurate and practical reference for future land management and restoration strategies in areas affected by landslides.

4.4.   Limitations and uncertainties

We successfully incorporated vegetation succession patterns and integrated support vector machine (SVM) classification, the Markov chain model, and remote sensing imagery, yielding effective results in land use classification and the simulation of dynamic landscape changes. However, several uncertainties affected model performance and prediction reliability.

4.4.1.   Spatial resolution of remote sensing imagery

We utilized 30-meter resolution Landsat satellite imagery, which faced challenges such as mixed-pixel effects and classification ambiguities, particularly along land cover boundaries. These limitations were notable in the classification of grassland and forest due to similar spectral reflectance, which led to misclassification.

4.4.2.   Spectral similarity between land cover classes

The spectral similarity between grassland and forest further complicated spectral-based classification methods. This overlap in spectral characteristics reduced the accuracy of classification. In addition, it was difficult to distinguish these two land cover types with high precision.

4.4.3.   Assumption of stationarity in the Markov Chain Model

The Markov chain model assumes stationarity in transition probabilities, which oversimplifies the complex interactions between land use categories. This assumption fails to account for abrupt changes in land use that may be triggered by extreme weather events or human activities. It increases the predictive uncertainties, particularly in long-term simulations.

4.4.4.   Limited consideration of external factors

The model does not fully account for dynamic external factors such as human intervention, climate change, or sudden environmental disturbances. These factors may significantly alter land use patterns and reduce the ability of the model to accurately predict long-term landscape changes.

It is suggested to address these limitations and enhance the model performance, future focus on utilizing higher-resolution imagery to reduce mixed-pixel effects, incorporate additional auxiliary data (such as soil properties, land management data, or more detailed topographic features), and apply more advanced modeling techniques that can better capture non-stationary dynamics in land use transitions. These enhancements will improve classification accuracy, better reflect the complexities of landscape change, and increase the overall reliability and predictive power of the model.

5.   Conclusions

We aimed to explore dynamic vegetation restoration in the Chiufenershan area, where the large-scale landslide occurred more than 25 years ago. A combination of multi-temporal remote sensing imagery, image classification, and land use change techniques were utilized to investigate vegetation succession sequences and land cover changes. The results suggested that the secondary forest stage can be reached approximately seven years after the landslide in the deposited area. In contrast, the recovery of the collapsed area may require a longer period due to the influence of environmental stressors. It is anticipated vegetation succession to remain unstable until 2075. Although the study used a significant volume of remote sensing data and landscape change models to simulate vegetation restoration, there were uncertainties resulting from the resolution of remote sensing images, the similarity of object reflection spectra, and the foundational assumptions of model theory. This may influence the reliability of the model. It is recommended that future research integrate higher-resolution remote sensing images with pertinent auxiliary information to enhance the accuracy of vegetation succession simulations.

Author contributions

Chih-Wei Chuang contributed to the conceptualization and design of the study, supervised the research activities, conducted the investigation and methodology development, and was responsible for the original draft preparation, as well as the review and editing of the manuscript. He also administered the overall project. Hao-Yu Huang was responsible for formal data analysis and visualization. He also contributed to the critical review and editing of the manuscript. Chun-Wei Tseng participated in the investigation process and contributed to the review and editing of the manuscript.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was supported by the 2024 Agricultural Science-7.1.2-Forestry-01 Science and Technology Project of Taiwan Forestry Research Institute, MOA, R. O. C. (Taiwan). The authors would like to thank the anonymous reviewers for their insightful comments and suggestions on an earlier version of this manuscript.