Abstract

Landslides are major natural hazards in mountainous regions, frequently triggered by earthquakes and causing significant vegetation loss and long-term ecological degradation. This study evaluates the impact of the 2018 earthquake-triggered landslides on vegetation dynamics in Mount Rinjani National Park (MRNP), Indonesia, using multi-temporal Landsat 8 and 9 Normalized Difference Vegetation Index (NDVI) data from April 2018 to February 2026. NDVI values were extracted from two landslide-affected sites, the Central Rinjani site and the East Rinjani site, to analyze vegetation degradation and recovery across growing and non-growing seasons. A total of 341 statistically significant correlations were identified before multiple-comparison correction, with 281 retained after Benjamini-Hochberg false discovery rate adjustment (q < 0.05). Results indicate a pronounced decline in NDVI during the 2018-2019 disturbance period, followed by a gradual recovery. The Central Rinjani site exhibited positive but statistically non-significant recovery trends (+0.0041 NDVI·yr^-1 growing season; +0.0036 NDVI·yr^-1 non-growing season), while the East Rinjani site showed higher growing-season NDVI (mean = 0.592 vs. 0.540), likely due to more favourable local moisture conditions. Pre-disturbance NDVI was positively associated with recovery trajectories, indicating its value as a predictor of post-disturbance regeneration. However, NDVI values remained below pre-disturbance levels by February 2026, suggesting incomplete ecological recovery. These findings highlight the effectiveness of Landsat NDVI for monitoring post-landslide vegetation dynamics and underscore the need for continued long-term observation.

Keywords

NDVI, Landslide, Vegetation Recovery, Landsat, Remote Sensing, Earthquake Disturbance

Share and Cite:

1. Introduction

Landslides are natural disasters that constitute one of the biggest environmental problems in mountainous areas worldwide [1]. Landslides may strip the ecosystem, infrastructure, and human settlements of vast territory and lead to significant changes in the natural landscape [2]. A multiplicity of factors are usually responsible for these events, including severe rainfall, seismic activity, volcanic processes, and anthropogenic changes in land use [3]. The topography of mountainous areas, rough slopes, loose soils, and high levels of precipitation make the terrain vulnerable to slope failure [4]. Consequently, vegetation cover is often eroded, soil is washed away, habitats are destroyed, and ecological imbalances are caused by landslides in the long term [5]. Knowledge of the environmental effects of landslides thus plays a crucial role in effective environmental monitoring, ecosystem restoration, and disaster risk management [6].

Vegetation plays a central role in restoring slope stability following landslide disturbances, as the return of tree root systems helps anchor soil and reduce erosion [7]. But this recovery process might be prolonged, and not all sites, including those located in harsh conditions, will ever recover [8]. While it is stated that human-assisted revegetation can have positive outcomes in certain studies, mountainous terrain frequently prevents human involvement, so recovery in those locations is likely to be facilitated by natural processes [9].

Mount Rinjani is one of the most prominent volcanic mountains in Indonesia and represents a critical ecological and environmental system within the region [10]. Mount Rinjani is rugged in terms of topography, with steep slopes and intricate geological features due to its elevation of about 3726 meters above sea level [11]. Rainfall is very high, and volcanic activity is dynamic in the region, which collectively increases susceptibility to landslides and other natural disasters [12]. The presence of landslides in this region may also cause considerable changes in vegetation distribution, land cover degradation, and the disruption of ecological processes [13]. Other adverse impacts of such environmental disruptions can include effects on biodiversity, watershed stability, and soil fertility [14]. As a result, vegetation changes in landslide-prone areas like Mount Rinjani should be monitored to understand the ecosystem’s response to natural disasters and to aid in sustainable land management strategies [15].

The development of remote sensing technology has provided a potent means of monitoring the environment and assessing damage in the event of disasters, especially in recent decades [16]. The large land-surface areas and long observation durations of satellite-based observations enable researchers to analyze land-surface conditions [17]. The satellite data obtained via the Landsat program are long-term and consistent, and they are free, which is very common in environmental research [18]. Such satellite observations help scientists trace changes in land cover, the dynamics of vegetation, and other environmental disturbances at various times [19].

One of the most widely applied techniques for assessing vegetation conditions from satellite imagery is the Normalized Difference Vegetation Index (NDVI), a spectral index derived from the red and near-infrared bands that reflects vegetation health, density, and productivity [20] [21]. NDVI provides a reliable indicator for detecting vegetation loss and monitoring recovery, making it well-suited for evaluating ecosystem changes following disturbance events [22] [23]. However, few investigations have combined precise landslide mapping with rigorous multi-temporal NDVI analysis over extended periods exceeding five years, particularly within Indonesian seismic regions. This gap is especially significant in seismically active volcanic mountain environments such as Mount Rinjani, where earthquake-triggered mass movements are frequent, yet post-disaster ecological assessments remain rare, hindered by logistical constraints, inadequate monitoring infrastructure, and the remoteness of the affected terrain.

Therefore, the main objective of this study is to analyze vegetation dynamics in landslide-affected areas of Mount Rinjani National Park over the period from 2018 to 2026 using multi-temporal Landsat satellite imagery. Specifically, the study aims to: 1) quantify the extent of vegetation loss caused by the 2018 earthquake-induced landslides; 2) assess post-disturbance vegetation recovery trajectories across growing and non-growing seasons; and 3) evaluate inter-annual NDVI relationships to identify spatial and temporal patterns of vegetation change within the affected gullies.

2. Method

2.1. Study Area

The study was conducted within Mount Rinjani National Park (MRNP), located on Lombok Island, West Nusa Tenggara Province, Indonesia. The park encompasses approximately 41,330 hectares, with terrain spanning elevations from 600 to 3726 m above sea level. Two study sites were selected to represent the landslide-affected zones generated by the 2018 seismic sequence: the Central Rinjani site and the East Rinjani site. Both sites are situated at elevations of approximately 1700 - 2200 m.a.s.l., on slopes ranging from near-flat to exceeding 40˚, with predominantly south, southeast, and east-facing aspects.

A total of 50 fixed observation points were established across the two study sites: 26 points within the Central Rinjani Site and 24 points within the East Rinjani Site, as shown in Table 1. Points were distributed using a systematic grid approach in Google Earth Pro (version 7.3.7.1094), with a minimum spacing of 30 m between adjacent points to ensure spatial independence while capturing terrain variability across slope, aspect, and distance from the main gully channel. All points were located within mapped landslide scar boundaries delineated from the 2018 post-earthquake imagery, ensuring that all observations directly represent disturbed vegetation. The same fixed geographic coordinates (WGS84 decimal degrees) were used for NDVI extraction at every image date throughout the 2018-2026 study period, enabling a true time-series comparison at each location.

Table 1. Presents the geographic coordinates and distribution of sampling points across the two study sites.

The East Rinjani site represents a severely disturbed area subject to direct earthquake-induced mass movements, characterized by steep volcanic slopes and initially sparse post-disturbance vegetation cover. The Central Rinjani site served as a comparative location with similar geomorphological and climatic conditions but distinct recovery dynamics. The regional climate is monsoonal, with a wet season from September to March and annual rainfall of 2000 - 4000 mm. The park’s geology is predominantly Quaternary volcanic rock, contributing to inherent slope instability.

2.2. Satellite Data Acquisition and Pre-Processing

Multi-temporal satellite imagery was obtained from the Landsat program via the United States Geological Survey (USGS) Earth Explorer platform (https://earthexplorer.usgs.gov/). Images were acquired across multiple dates spanning April 2018 to February 2026, corresponding to the period before, during, and after the 2018 earthquake-triggered landslide events. A total of 60 temporal NDVI datasets were generated from the downloaded imagery, covering both the NDVI growing season (October to March) and the NDVI non-growing season (April to September). This seasonal classification is consistent with the monsoonal climate of Lombok, where the wet season from September to March drives peak vegetation greenness, while the dry season from April to September corresponds to reduced vegetation productivity [24]-[26]. Images acquired during the non-growing season (April to September) were prioritised in the time series, as cloud cover during the wet growing season frequently obscured imagery, consistent with the approach adopted for the same study area by [27] and broader remote sensing practices in tropical regions [28]. All images were pre-processed to surface reflectance level using the USGS Collection 2 Level 2 product, which applies atmospheric correction procedures to enable consistent multi-temporal comparison. Images were included if: 1) the acquisition date fell within the study period (April 2018 - February 2026); 2) cloud cover over the study area was less than 10%; and 3) the image was available as a USGS Collection 2 Level 2 surface reflectance product. Images were excluded if cloud/shadow masking revealed that more than 20% of the sampling points were contaminated, as shown in Figure 1.

2.3. NDVI Calculation and Extraction

NDVI was computed using the standard formula on the atmospherically corrected surface reflectance data. Radiometric calibration was applied using the Landsat Surface Reflectance scaling equation, SR = (DN × 0.0000275) − 0.2, to convert raw digital numbers into surface reflectance values from the Landsat 8/9 under the Landsat Program and NDVI = (NIR − Red)/(NIR + Red), where NIR and Red are Band 5 and Band 4, respectively, of the Landsat 8/9 OLI/TIRS Collection 2 Level-2 imagery. NDVI ranges from −1 to +1, where higher values indicate denser vegetation. Zonal statistics were applied to each point in the image using a 3 × 3-pixel kernel centered on the coordinate to extract NDVI values and eliminate single-pixel noise. All available dates for all images were extracted during the study period, resulting in complete multi-temporal data for both study locations. All processes were conducted using the latest version of QGIS (3.44.7-Solothurn).

2.4. Validation

The NDVI dataset was validated using spatial, comparative, and temporal approaches. Landslide-affected areas consistently showed lower NDVI values than stable forest control sites, confirming NDVI’s ability to distinguish disturbance from intact vegetation. Temporal trends aligned with the 2018 earthquake events, showing a clear decline followed by a gradual recovery. A 3 × 3-pixel kernel was applied to reduce noise and improve data reliability. Despite the absence of field data, the combined spatial agreement, control comparison, temporal consistency, and virtual ground truthing provide strong indirect validation of the results.

2.5. Statistical Analysis

All possible pairs of the multi-temporal NDVI variables were evaluated using the Pearson correlation coefficient. The Pearson correlation formula was applied to each pair of NDVI variables (n = 50 observations per pair) to quantify the strength and direction of temporal relationships in vegetation dynamics across the 2018–2026 period. From the initial 60 NDVI datasets, one duplicate variable was removed, resulting in 59 variables used in the correlation analysis and a total of 1,711 pairwise comparisons. For each correlation, the unit of analysis was the 50 observation points, with NDVI values extracted at those points for each image date (n = 50 per pair). To control for multiple comparisons, the Benjamini–Hochberg false discovery rate (FDR) procedure was applied at a target FDR of 5% (q < 0.05). This method reduces the likelihood of false positives by adjusting p-values, ensuring that only statistically reliable correlations are retained for subsequent analysis. The total number of pairwise comparisons among NDVI variables was determined using the combination formula:

Figure 1. Methodological workflow for multi-temporal NDVI analysis of vegetation recovery, including surface reflectance conversion, extraction of Landsat 8/9 OLI Red and NIR bands, NDVI computation, virtual ground truth sampling, and statistical accuracy assessment.

$\left(\begin{array}{c}n\\ 2\end{array}\right)=\frac{n\left(n-1\right)}{2}$

where n represents the total number of NDVI variables. Following the removal of one duplicate variable, a total of 59 NDVI variables were retained, resulting in 1,711 pairwise comparisons. Correlation strength was classified as: Very Strong (|r| ≥ 0.80), Strong (0.60 ≤ |r| < 0.80), Moderate (0.40 ≤ |r| < 0.60), and Weak (|r| < 0.40). All analyses were conducted in Python (Version 3.11.7) using SciPy and pandas’ libraries.

Pairwise comparison formula:

$$ r = \frac{\sum \left(x_i-\stackrel{`}{x}\right)\left(y_i-\stackrel{`}{y}\right)}{\sqrt{\surd \sum (x_i-\stackrel{`}{x}{)}^2\sum (y_i-\stackrel{`}{y}{)}^2}}$$

(FDR) Benjamini-Hochberg Formula:

$p\left(1\right)\le p\left(2\right)\le \cdots \le p\left(m\right)$

Then find the largest k such that:

$p_{\left(k\right)}\le \frac{k}{m}\cdot Q$

3. Results

3.1. Characteristics of the Study Area

Mount Rinjani National Park (MRNP) is located on Lombok Island at approximately 8.42˚S and 116.45˚E (Figure 3). The park exhibits a tropical monsoonal climate strongly influenced by elevation, with temperatures decreasing from approximately 20˚C in lowland areas to below 10˚C at higher elevations. However, mean temperatures across the region generally remain above 26˚C. Rainfall is highly seasonal, occurring predominantly between September and March, with annual precipitation ranging from 2000 to 4000 mm, consistent with regional climate patterns in eastern Indonesia [29] (Figure 2). The park covers around 41,330 hectares, with terrain rising from 490 to 3726 meters above sea level [30]. Landslide-affected sites within the study area were distributed between 1700 and 2200 meters in

Source: Adopted from [27].

Figure 2. Average temperature (a) and annual rainfall index (b) in Mt. Rinjani National Park from 2018 to 2024 (Indonesian Agency for Meteorology 2025).

elevation, on slopes ranging from nearly flat to over 40˚, and were predominantly oriented toward the south, southeast, and east [31] (Figure 4). The park’s geology is largely composed of Quaternary volcanic rock formations, which contribute to the area’s inherent slope instability [32].

Mount Rinjani National Park has a well-documented history of mass movements driven by high-magnitude seismic events and intense monsoonal rainfall. The 2018 earthquake sequence alone generated at least 10,521 coseismic landslides across the park’s steep volcanic terrain. Research has consistently identified steep volcanic morphology, unconsolidated geological material, prevailing land cover conditions, rainfall intensity, and tectonic activity as the primary factors governing landslide occurrence in the region [33]. This combination of recurring disturbance drivers and extreme triggering events makes Mount Rinjani National Park (MRNP) a particularly suitable location for studying post-earthquake vegetation recovery dynamics [34].

Source: Adopted from [27].

Figure 3. Study site: (a) Republic of Indonesia, (b) Lombok Islands, (c) Mt. Rinjani National Park.

The 2018 seismic sequence comprised five major events occurring on July 29, August 5, August 9, August 18, and August 19, registering magnitudes of 6.4, 7.0, 5.9, 6.4, and 6.9, respectively [27]. These earthquakes produced a range of mass movement types, including rockfalls, debris flows, debris slides, and lateral spreads across the park’s slopes [11].

Source: Adopted from [27].

Figure 4. Study site characteristics in elevation (a), slope gradient (b), and slope aspect (c).

In terms of biodiversity, Mount Rinjani National Park supports a diverse assemblage of plant species distributed across five distinct elevational zones. The lowland zone, below 1000 meters, is characterized by Shorea spp. and Ficus spp. Between 1000 and 2000 meters, the submontane zone is dominated by Schima wallichii, Podocarpus spp., Casuarina junghuhniana, and Imperata cylindrica. The montane zone, spanning 2000 to 3000 meters, is defined by Anaphalis javanica and Vaccinium varingiaefolium. At subalpine elevations between 3000 and 3726 meters, vegetation shifts to Anaphalis javanica, Vaccinium varingiaefolium (Blume) Miq., alpine grasses, and mosses. Above 3726 meters, the summit zone supports virtually no vegetation owing to the severity of environmental conditions at that altitude [27].

3.2. Normalized Difference Vegetation Index (NDVI) Temporal Trends (2018-2026)

The temporal NDVI of the Central and East Rinjani sites showed a similar pattern of vegetation degradation during the disturbance period of 2018-2019, after which an upward recovery pattern took hold (Figures 5-10). The minimum NDVI values at the two sites occurred during the shared disturbance period (2018-2019), which was related to the 2018 earthquake sequence and mass-movement events.

At the Central Rinjani site, there was a marked decline in NDVI values during both the growing and non-growing seasons following the 2018 disturbance events. The linear trend lines fitted to the post-disturbance data show a consistently positive slope at both sites across all seasons, indicating a directional recovery trajectory; however, neither trend reached statistical significance at p = 0.05 (Central Rinjani site: r = 0.103, p = 0.576; non-growing season: r = 0.112, p = 0.703), likely reflecting the limited number of observations and high inter-annual variability characteristic of early successional recovery.

Figure 5. Combined NDVI overlay for the Central Rinjani site and the East Rinjani site, 2018-2026, illustrating synchronous disturbance responses and divergent post-disturbance recovery rates.

Figure 6. Temporal NDVI trends at the Central Rinjani Site, growing and non-growing seasons, 2018-2026. The shaded area indicates the 2018-2019 disturbance period.

A comparable initial degradation had been observed at the East Rinjani Site; the post-disturbance recovery trend was relatively steeper, and the NDVI values started to approach the magnitude of the pre-disturbance level much more quickly. This difference might reflect site-specific variations in residual seed banks, soil moisture, or microclimate conditions. Growing-season NDVI values were slightly higher than at the Central Rinjani Site, with a mean of 0.592 (standard deviation = 0.094), ranging from 0.484 to 0.711. However, only one non-growing season observation was available for the East Rinjani Site (NDVI = 0.486), which precluded robust seasonal trend analysis and limited direct comparability with the Central Rinjani Site. The increased vegetation productivity implied by the higher growing-season NDVI at the East Rinjani Site could be attributed to more favorable local moisture conditions associated with its topographic position.

Figure 7. Temporal NDVI trends at the East Rinjani Site, growing and non-growing seasons, 2018-2026. The shaded area indicates the 2018-2019 disturbance period.

Figure 8. Seasonal NDVI distributions at the central and east sites across the 2018-2026 study period. Blue boxes represent the growing season (Wet, Nov-Apr) and orange boxes represent the non-growing season (Dry, May-Oct). Median, interquartile range, whiskers, and outliers are shown. Growing season NDVI was consistently higher at both sites, reflecting the prioritisation of wet-season imagery and the influence of monsoonal precipitation on vegetation productivity.

Seasonal change was observed in the record at the two sites, where growing-season NDVI was always greater than non-growing-season NDVI, indicating the effect of monsoon precipitation on vegetation productivity.

The linear trend lines exhibited positive but statistically non-significant slopes, indicating a directional recovery trend without statistical significance. The combined NDVI overlay of the two sites, as shown in Figure 5, indicates a similarity in the disturbance responses and a dissimilarity in the recovery rates.

3.3. Pre- and Post-Disturbance (Earthquake-Induced) NDVI Comparison

NDVI declined during the 2018-2019 disturbance period, with reductions becoming more pronounced in the months following the initial seismic events. While mean NDVI across all 50 observation points decreased slightly from 0.781 ± 0.136 (April 2018) to 0.748 ± 0.220 (May 2018), this initial change was not statistically significant (paired t-test: t (49) = 1.05, p = 0.297). It should be noted that the values presented in Figure 8 represent aggregated temporal trends, whereas the statistical analysis is based on point-level NDVI extraction, which may result in slight differences in magnitude. However, both approaches consistently indicate that the most substantial decline occurred later within the disturbance period, with NDVI reaching its lowest values between late 2019 and early 2020. These reductions were most pronounced at observation points located on steeper slopes within the direct pathways of mass movements.

Figure 9. Pre- and post-eruption NDVI comparison at both study sites, quantifying the magnitude of vegetation loss associated with the 2018 earthquake sequence.

3.4. Pearson Correlation Analysis

The multi-temporal NDVI dataset analysis showed, through Pearson correlation analysis, that 341 statistically significant pairs had been identified before multiple-comparison correction, with most of the positive correlations being Very Strong and Strong across adjacent seasonal NDVI values.

Figure 10. Inter-site NDVI correlation between central and east Rinjani sites across all matched observation dates (2018-2026). Each point represents a paired site-mean NDVI observation, colour-coded by year (red = 2018, dark green = 2025). The dashed red line indicates the OLS regression fit (ŷ = 1.08x + 0.01), and the dotted grey line represents the 1:1 reference. Pearson r = 0.899, p < 0.001.

Table 2 presents the top 15 correlation pairs ranked by absolute correlation coefficient.

The perfect correlation observed between NDVI202208 and NDVI2022_1 (r = 1.000, p < 0.001) warrants clarification. These two variables were derived from the same underlying Landsat scene but were retained as separate columns due to differences in naming during data processing. As both variables represent identical NDVI values extracted from the same image at the same 50 observation points, the observed perfect correlation reflects duplicate information rather than independent measurements. Accordingly, NDVI2022_1 was excluded from subsequent analyses to avoid redundancy and potential inflation of correlation statistics.

The most significant correlations (r > 0.95) were found between the value of NDVI in the same season across adjacent years (as shown in Table 2), especially between NDVI 202309 and NDVI 202304 (r = 0.987; NDVI 202308 and NDVI 202309, r = 0.959), and spatially consistent relative vegetation productivity trends at monthly time scales. A significant negative correlation between dry-season and wet-season NDVI values (e.g., NDVI202508 and NDVI202510, r = −0.778) indicates natural phenological differences. In this growth, the wet period, which is moisture dependent, is not directly proportional to a sustained relative amount of dry-season greenness at the same sites.

Table 2. Top 15 Pearson correlation pairs from the multi-temporal NDVI dataset (ranked by absolute correlation coefficient; statistical significance is based on pre-correction p-values, with FDR-adjusted results discussed in the text).

The correlation over the years between NDVI before the disturbance (e.g., NDVI201804) and the values of the recovery period (e.g., NDVI202408, r = 0.821; NDVI202510, r = 0.609) shows that the sites with higher pre-disturbance vegetation cover were more likely to recover better, which fits the ecological succession theory and the importance of residual biomass in post-disturbance regeneration. These findings confirm that pre-disturbance vegetation conditions are a significant predictor of long-term recovery patterns in the MRNP landslide landscape. This recovery trajectory is consistent with findings from comparable studies of post-earthquake vegetation dynamics in tropical volcanic environments, where natural regeneration processes typically produce measurable but slow increases in canopy cover over multi-year to decadal time scales [33] [35].

4. Discussions

4.1. Post-Disturbance NDVI Decline and Comparison with Other Earthquake-Affected Regions

The significant reduction in NDVI at both the Central and East Rinjani Sites during the disturbance period between 2018 and 2019 is consistent with a growing body of literature documenting vegetation loss after major seismic events in volcanic mountainous regions. Studies reporting on the 2008 Wenchuan earthquake in Sichuan, China, found that the maximum vegetation damage (minimum NDVI 0.41) lasted up to one year after the main earthquake, as subsequent hazards, such as debris flows and slope reactivation, persisted in clearing vegetation cover outside the narrow disturbance window [36]. The present study also observed a similar lag in peak damage, as NDVI values at both sites were lowest throughout the 2018-2019 interval, but not immediately after the July-August 2018 earthquakes. Such a trend indicates that further geomorphic instability, coupled with extreme monsoonal precipitation in the months after the earthquakes, lengthened the duration of active vegetation clearing at both study sites.

The post-earthquake NDVI decline and the landslide scar distribution are also spatially related in the present study, as in the study by [37] that used Landsat imagery to map vegetation damage and recovery conditions following the earthquake in Wenchuan, which also showed that regions with the lowest NDVI were spatially clustered along the coseismic landslide boundaries. Their approach of calculating a Measure of Vegetation Damage Assessment (VDA), which is an offset between pre- and post-earthquake NDVI differencing, has a conceptual similarity to the pre- and post-eruption NDVI differences in the current study, providing methodological similarity across independent research environments.

4.2. Rate and Completeness of Vegetation Recovery

Although statistically insignificant, the NDVI recovery slopes over the study period of 2018-2026 were consistently positive at the Central Rinjani Site in both growing and non-growing seasons (+0.0041 NDVI/yr and +0.0036 NDVI/yr, respectively), indicating a directional trend of gradual vegetation recovery. These are lower rates compared to those reported for similar post-earthquake volcanic environments [38]. Studies of post-seismic vegetation succession in the Wenchuan area with ten years of MODIS data showed that the recovery rates of regional NDVI were about 0.006 NDVI units per year in the first recovery period, but the recovery rates were not similar across elevations and slope aspects. The reduced pace of recovery at Mount Rinjani is probably due to the increased altitude of the study sites (1700 - 2200 m.a.s.l.), the occurrence of a steep Quaternary volcanic substrate with minimal soil development, and the continuation of geomorphic instability due to frequent slope movements triggered by rainfall after the major seismic disturbance.

The observation of a failed NDVI recovery to pre-disturbance levels by 2026, eight years after the 2018 earthquakes, is not isolated but rather relates to the larger body of literature concerning recovery timescales in high-elevation tropical volcanic settings. In a study of post-landslide vegetation recovery in the southern region of Mount Rinjani National Park between 2018 and 2024, using Sentinel-2 imagery, [27] established that it took several years before full vegetation recovery on landslide scars occurred and that recovery remained incomplete across many of the sites during the six-year period (especially on steep slopes and at higher elevations). The agreement between the results of the study on the same topic using Sentinel-2 and the current one using Landsat leads to the conclusion that recovery at MRNP is a long-term process that occurs largely due to the specifics of the terrain.

On the volcanic landscape level, [39] assessed the post-eruption vegetation recovery at the Unzen volcano in Japan with Landsat time series and discovered that the rates of recovery and the extent of its comprehensiveness were heavily controlled by substrate heterogeneity, proximity to remnants of vegetation cover after the eruption, and occurrences of stochastic disturbances. Their finding of spatially inhomogeneous recovery in volcanic substrates, with certain regions being characterized by fast pioneer colonization and other regions remaining bare over decadal timescales, can be directly related to the differences in recovery processes that have been witnessed between the Central Rinjani Site and the East Rinjani Site in this study. This steeper recovery curve at the East Rinjani Site is probably the result of better local conditions in terms of propagule recruitment and soil moisture retention compared to the Central Rinjani Site.

4.3. Topographic Controls on Recovery and the Role of Aspect

The observed difference in the recovery rates and the values of NDVI in the growing season is in line with the literature on the topographic control of vegetation regeneration during post-landslide periods in volcanic environments, where the NDVI values are higher in the East Rinjani Site (mean = 0.592) compared to the Central Rinjani Site (mean = 0.540) [40]. A study of rainfall-induced and coseismic landslide scars at the Aso volcano in Japan did find a higher rate of grass vegetation recovery on north-to-west facing slopes than on south- and east-facing slopes, which the authors attribute to different levels of sunlight and, therefore, higher soil moisture levels on the sheltered slopes. The fact that the East Rinjani Site has experienced more rapid recovery can be attributed partly to comparable moisture-related topographic controls, since the geomorphological location of the valley in a gully system would favor higher convergence of moisture compared to more exposed slope locations, such as those found in the Central Rinjani Site.

The aspect of slope gradient as an obstacle to recovery also needs to be addressed. The slope at both study sites was steeper than 40˚, and the literature indicates that steep gradients are a major barrier to vegetation recovery following disturbances [41]. Analyzing the frequency and vegetation recovery across Taiwan region using Landsat time series from 1990 to 2022, the study discovered that lower elevations and gentle slopes on the island had considerably shorter vegetation recovery durations, and high-frequency disturbance sites on steep slopes had a growth impediment that exceeded their study period. The discovery that over 50% of the vegetation at the post-typhoon landslide sites in Taiwan region had recovered to pre-event levels during the period observed, in contrast to the partial recovery found at MRNP after eight years, indicates that the combination of elevation, volcanic substrate, and compound seismic-triggered events at Rinjani imposes a stricter and more enduring recovery limit than the events triggered by rainfall at typical lower-elevation sites.

4.4. Seasonal NDVI Variability and Phenological Dynamics

The observed seasonal variation in NDVI reflects monsoonal climatic control, with higher vegetation productivity during the wet season and reduced greenness during the dry season. Similar impacts due to monsoonal precipitation have also been described in the seasonal NDVI variations throughout the Indonesian vegetation system by [42]. In karst savannas in East Nusa Tenggara, Indonesia, the input of wet-season precipitation generated strong positive NDVI responses at all study locations, whereas the lack of moisture in dry seasons had a consistent negative impact on vegetation productivity. Equally, [43] reported globally that grassland and shrubland biomes in the tropics and subtropics have a growing seasonal NDVI amplitude, with greening in wet seasons surpassing browning in dry seasons, a trend that is directly similar to seasonal differences at MRNP. This is also consistent with spectral seasonality by leaf phenology that remains in the recovering or degraded vegetation of tropical seasonal biomes by [44], and the active recovery of photosynthetically active canopy that is observed in other tropical volcanic regions when earthquakes take place.

4.5. Pre-Disturbance NDVI as a Predictor of Recovery

The strong positive cross-year correlations observed between pre-disturbance NDVI (e.g., NDVI201804) and post-disturbance NDVI values (e.g., NDVI202408, r = 0.821; NDVI202510, r = 0.609) provide clear evidence that the initial vegetation condition is a key determinant of recovery trajectories. Sites with higher pre-disturbance NDVI values consistently exhibited stronger post-disturbance recovery signals, indicating greater resilience and regenerative capacity.

This finding is consistent with the biological legacy concept in disturbance ecology, where areas with higher pre-disturbance vegetation tend to recover more rapidly due to the presence of residual biomass, intact root systems, and persistent soil seed banks [37] [45]. These ecological legacies facilitate faster recolonization and reduce the time required for vegetation re-establishment following disturbance events. Similar relationships between pre-disturbance vegetation conditions and recovery rates have been widely reported across different disturbance contexts. For example, studies of post-fire vegetation recovery in Mediterranean ecosystems have demonstrated that areas with higher pre-fire NDVI return to baseline conditions more rapidly than sparsely vegetated sites, highlighting the importance of initial vegetation capital in controlling recovery dynamics [45]. Within the framework of Mount Rinjani National Park, the relatively lower pre-disturbance NDVI observed at the Central Rinjani Site compared to the East Rinjani Site may partly explain the slower, less pronounced recovery trajectory at that site. This suggests that spatial variability in initial vegetation conditions plays a critical role in shaping long-term recovery patterns within landslide-affected landscapes.

From a management perspective, these findings imply that areas with low pre-disturbance NDVI are less likely to recover rapidly through natural regeneration alone and may therefore benefit from targeted restoration interventions. Identifying such areas using pre-disturbance satellite data provides a practical and cost-effective approach for prioritizing ecological restoration efforts in mountainous and remote environments.

4.6. Methodological Considerations and Limitations

The present study employed Landsat 8 and 9 OLI/TIRS Collection 2 Level 2 imagery at 30 m spatial resolution, which imposes a fundamental constraint on the detection of fine-scale vegetation recovery processes, particularly in the early post-disturbance period when recovering vegetation consists of small herbaceous plants and seedlings with low canopy cover. This limitation is well-recognized in the remote sensing literature. Studies using higher-resolution platforms, such as [40], who employed RapidEye (5 m) and PlanetScope (3 m) imagery at Aso volcano, or UAV-based photogrammetry approaches, have demonstrated that vegetation recovery in the early succession stages is often underestimated by Landsat-resolution sensors due to the dominance of bare substrate within mixed pixels. The moderate NDVI values recorded at both study sites (growing-season means of 0.540 and 0.592) likely underestimate actual plant cover at finer spatial scales, and future studies should consider integrating Sentinel-2 (10 m) or higher-resolution data to resolve intra-gully variability in recovery status. An additional limitation of this study is the uneven temporal distribution of seasonal observations between the two study sites. Specifically, the East Rinjani Site had only one non-growing season observation available throughout the study period. This limited the ability to perform robust seasonal comparisons and reduced the statistical reliability of inter-site comparisons for non-growing season dynamics. As a result, interpretations of seasonal variability and recovery patterns between the two gullies should be considered with caution. Future studies should prioritize consistent seasonal sampling across all study sites to improve comparability and analytical robustness.

5. Conclusion

This paper assessed the post-disturbance vegetation dynamics at two sites with coseismic landslides in Mount Rinjani National Park through multi-temporal Landsat NDVI analysis between 2018 and 2026. This sequence of earthquakes caused quantifiable vegetation loss at both the Central and East Rinjani Sites, and NDVI decreased significantly during the 2018-2019 disturbance period and then began to slow down at various times due to natural succession. At the Central Rinjani Site, NDVI trends were positive but statistically non-significant in both the growing (+0.0041 NDVI/yr) and non-growing seasons (+0.0036 NDVI/yr). However, at the East Rinjani Site, there were relatively higher growing-season NDVI values (mean = 0.592 vs. 0.540), probably due to more favourable local hydrological conditions. The correlation between pre-disturbance NDVI and recovery trajectory was assessed using Pearson correlation on 281 statistically significant pairs (after Benjamini-Hochberg FDR correction, q < 0.05), demonstrating temporal coherence of the Landsat-derived signal and supporting the validity of pre-disturbance NDVI as a useful predictor of recovery trajectory. A positive slope confirms that vegetation is recovering over time (directional improvement), while NDVI values remaining below pre-disturbance baselines validate that recovery is incomplete as of February 2026. These findings demonstrate the utility of multi-temporal Landsat NDVI in monitoring post-landslide recovery following the 2018 earthquake at Mount Rinjani, Indonesia, in remote volcanic environments and highlight the need for extended observation beyond 2026 to determine whether natural regeneration alone is sufficient to restore the ecological integrity of the affected catchments.

Ethics Statement

This study is a retrospective spatiotemporal analysis conducted entirely using freely available, publicly distributed satellite imagery obtained from the United States Geological Survey (USGS) Earth Explorer platform. The research does not involve human or animal subjects, personally identifiable information, or any form of field intervention. Accordingly, no institutional ethics committee review or informed consent was required. The authors confirm that no human participants, vulnerable populations, or protected species were involved in any aspect of this research.

Acknowledgements

The authors would like to express sincere gratitude to Bu Novia Lusiana, Ph.D., for her exceptional guidance and skill-based instruction throughout the Spatial Technology Course. Her dedication to imparting practical knowledge has greatly contributed to the development of the competencies applied in this research.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References