Effects of Topography on Vegetation Recovery after Shallow Landslides in the Obara and Shobara Districts, Japan

Abstract

Intense rainfall-induced shallow landslides can have severe consequences, including soil erosion and vegetation loss, making in-depth research essential for disaster risk management. However, vegetation recovery processes after shallow landslides and their influencing multivariate factors are not well known. This study aims to address this gap by investigating the vegetation recovery processes after shallow landslides and the impact of topography on this recovery. We focus on two regions in Japan: the Shobara district in Hiroshima Prefecture and the Obara district in Aichi Prefecture. The Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) derived from long-term Landsat images, as well as aerial photographs and environmental datasets, are used to measure vegetation recovery. Then, statistical analysis and the Seasonal Autoregressive Integrated Moving Averages (SARIMA) model were employed to investigate the dynamic response of vegetation under different combinations of environmental conditions using NDVI and EVI time series. Historical aerial photographs and vegetation index trend analysis suggest that vegetation in the study areas will take more than ten years to return to a stable state. The results also demonstrate the influence of atmospheric and land cover conditions when monitoring vegetation response using NDVI and EVI. In Obara, concave and convergent terrain positively influenced NDVI, while non-steep, low-elevation, and north-facing terrain positively influenced EVI. In Shobara, gentle and northwest-facing slopes were positively correlated with NDVI, and gentle and west-facing slopes were positively correlated with EVI. SARIMA modeling found that NDVI is more suitable for modeling the middle and late stages of vegetation recovery within 10–25 years after the landslide. In comparison, EVI is better for modeling the early stage of vegetation recovery within 10 years after the landslide.

1. Introduction

Rainfall-induced landslides can have severe consequences on communities [1]. Japan is prone to rainfall-induced landslides due to its geological, geomorphological, and climatic conditions [2]. Examining vegetation recovery following shallow landslides can offer valuable insights into improving future disaster recovery efforts.

Recent advances in Geographic Information Systems (GIS) and remote sensing technology have enabled researchers to collect detailed information to monitor vegetation recovery over space and time [3,4,5,6]. Vegetation indices are commonly used to extract data from multiple spectral bands to investigate vegetation changes [7,8,9,10,11]. These indices usually rely on the measurement of greenness using algebraic combinations between different spectral bands, particularly the red and near-infrared (NIR) bands [12,13]. The Normalized Difference Vegetation Index (NDVI) [14] and Enhanced Vegetation Index (EVI) [15], which are derived from NIR-based spectral indices, are frequently used to monitor vegetation recovery. The NDVI is chlorophyll sensitive, whereas the EVI is more responsive to canopy structural variations [16]. Li et al. (2010) compared and evaluated the accuracy of predicting natural vegetation cover using NDVI and EVI based on the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor by developing optimal regression equations for grassland, shrub, and forest samples. They found that NDVI has apparent advantages for predicting natural vegetation coverage over EVI [17]. Matsushita et al. (2007) presented a theoretical analysis of the topographic effect on EVI and NDVI using a non-Lambertian model. In a case study utilizing airborne-based images from a mountainous area with a dense Japanese cypress plantation, they found that EVI is more sensitive to topographic conditions than NDVI [18]. However, the comparative applicability and performance of NDVI and EVI in monitoring and simulating dynamic vegetation recovery in mountainous forest areas under various environmental conditions are still poorly understood. Thus, identifying appropriate and interpretable vegetation indices is necessary to comprehend vegetation responses after landslides under different environmental conditions.

The effectiveness of vegetation indices in representing vegetation recovery is heavily influenced by the type of landscape disturbance, soil characteristics, and topography [19]. For example, aspect is a crucial control of an area’s microclimate, localized hydrological variability, and soil hydrophobicity [20]. The slope gradient is critical in controlling soil moisture and affecting growing vegetation [21]. One or more topographic variables can significantly change the vegetation index response, and indices have different sensitivities to environmental conditions [22]. Vegetation also affects geomorphic processes through precipitation–vegetation–erosion feedback [23,24,25]. Considering the interdependence between topography and vegetation is crucial for effective vegetation recovery monitoring. Several studies have analyzed the suitability and relevance of different vegetation indices, considering several confounding factors, notably topographic variables, for vegetation monitoring [26,27,28].

Recent improvements in remote sensing technology and enhanced computing capabilities have provided researchers with unprecedented access to data, leading to a better comprehension of large-scale vegetation recovery over extended periods. These advances have created novel opportunities to investigate the complex connection between topography and vegetation recovery over time. Therefore, this study examines topography’s influence on vegetation recovery in landslide-affected areas in two study areas, the Shobara district in Hiroshima Prefecture and the Obara district in Aichi Prefecture, Japan. These areas were selected based on their susceptibility to shallow landslides owing to geological conditions and the availability of Landsat satellite images, aerial photographs of vegetation conditions, and environmental variables datasets. The study comprises four key components to achieve its objective. Firstly, we used historical aerial photographs to identify areas impacted by landslides. Secondly, we assessed vegetation recovery using vegetation indices derived from long-term satellite imagery. Thirdly, we analyzed correlations and time series data to evaluate the relationship between vegetation recovery and environmental variables. Finally, we compared the effectiveness of different vegetation indices for monitoring vegetation recovery under varying ecological conditions.

2. Materials and Methods

2.1. Study Areas

2.1.1. Obara District

The first study area is the Obara district in Toyota City, northern Aichi Prefecture, central Japan (Figure 1a,b). The region’s topography consists of low-relief mountains with elevations ranging from 200 to 600 m. This area experienced heavy rainfall on 12–13 July 1972, with precipitation of 218 mm per 5 h, resulting in numerous shallow landslides [29]. The affected areas consist mainly of weathered granitoid, making the area prone to shallow landslides [30]. The two main types of granitoid in the study area are coarse-grained granite (CG) and medium-grained granite (MG) [31]. Vegetation is scarce on cliffs, and the other areas are mainly covered with red spruce trees [32]. The subsequent post-landslide vegetation recovery process is natural.

2.1.2. Shobara District

The Shobara district is in Shobara City, in northeastern Hiroshima Prefecture of western Japan (Figure 1c,d). On 16 July 2010, a severe rainstorm caused over a thousand shallow landslides [33]. The region’s geology primarily comprises rhyolite (RL) and andesite (AS). The topography is characterized by low-relief mountains, with elevations ranging from 250 to 650 m. Deciduous broad-leaved forests cover most slopes, with some areas consisting of conifer forests [34].

2.2. Data

2.2.1. Aerial Photographs

Historical aerial photographs were obtained from the Geospatial Information Authority of Japan for different years: 1968, 1972, 1974, 1977, 1982, 1998, and 2000 for the Obara district and 1975 and 2010 for the Shobara district. We compared the pre- and post-landslide conditions of the study areas to create maps by visual inspection to show the distribution of damaged areas affected by the landslides (Figure 2). While areas with known human activities were excluded during the extraction of damaged areas, it should be acknowledged that some areas may have been affected by human activities, such as conversion to farmland, during the subsequent recovery period. However, the proportion of affected pixels in both study areas is less than 1%, rendering their impact insignificant for our analysis.

2.2.2. Satellite Imagery

Long-term time series satellite data were used to investigate the vegetation recovery processes. For Obara, Landsat 5 surface reflectance spanning 28 years (January 1984 to May 2012) with 30 m resolution was obtained from the National Aeronautics and Space Administration (NASA). Satellite imagery before 1984 was unavailable for Obara due to significant coverage gaps in Landsat 1–3 Multispectral Scanner (MSS) data. Furthermore, the absence of a blue band in Landsat 1–3 MSS hindered the calculation of EVI. For Shobara, Landsat 7 and 8 surface reflectance data were acquired at 30 m resolution over 22 years (January 1999 to November 2021) and eight years (April 2013 to December 2021), respectively. We then conducted other processing steps on the Landsat time series through Google Earth Engine (GEE), a cloud-computing platform capable of handling satellite imagery and geospatial datasets [35].

2.2.3. Geological Type Database

We obtained a geological distribution map from the geological map display system of the Geological Survey of Japan to investigate the impact of geological type on vegetation recovery. We generated a vector shapefile layer using the map to represent the geological conditions of various regions in the two study areas. This geological layer was co-registered with satellite imagery to analyze the correlation between geological type and vegetation recovery.

2.2.4. Digital Elevation Models

Digital elevation models (DEMs) were obtained from NASA’s Shuttle Radar Topography Mission (SRTM) digital elevation 30 m and used in GEE to calculate various topographic variables (

Table 1. Summary of topographic variables used for statistical analysis. The value ranges are those of data from the sources mentioned above.

Topographic Variables	Unit	Description	Value Range
			Obara	Shobara
Elevation (E)	m	Height of terrain above sea level	0–330	90–500
Slope gradient (S)	degrees	Slope gradient	0–80	0–80
Aspect (A)	degrees	Compass direction	0–360	0–360
Horizontal curvature (HC)	m	Curvature tangent to the contour line	−0.1–0.2	−0.1–0.2
Vertical curvature (VC)	m	Curvature tangent to the slope line	−0.01–0.04	−0.1–0.1
Maximal curvature (MC)	m	Highest values of curvature	−0.01–0.2	−0.01–0.3

). The consistent spatial resolution with Landsat imagery enabled direct comparison of calculated variables to Landsat data. These variables were selected based on their known impact on vegetation growth and survival. Elevation affects temperature and moisture, creating climatic gradients as elevation changes [36]. Slope gradient affects land stability and erosion rates and has been shown to impact vegetation recovery in eroded rills and gullies [21]. Aspect, the direction of a slope, creates spatial differences in solar insolation that can influence hydrologic, ecological, and geomorphic processes [20]. Additionally, curvature, divided into horizontal curvature, vertical curvature, and maximal curvature, is a crucial factor affecting soil moisture distribution and water transport, ultimately affecting vegetation growth [37].

3. Method

The methodology used in this study consists of three main stages: preprocessing of the satellite image time series, computation of vegetation recovery trends on a pixel-by-pixel basis, and identification of significant variables via statistical analysis to be utilized in subsequent time series analysis.

3.1. Preprocessing Landsat Time Series

We applied various filtering and processing techniques to preprocess the reflectance dataset acquired from Landsat 5, 7, and 8 satellite imagery. We used the damaged areas distribution map extracted from the pre- and post-landslide aerial photographs to constrain the Landsat data to the boundary of the damaged areas. The Landsat surface reflectance images have already been calibrated, converted to top-of-atmosphere reflectance, and atmospherically corrected [38] to ensure high quality images suitable for time series analysis. Data with less than 80% cloud cover over land and geometric root mean squared error of less than 10% were filtered, and terrain correction was applied to minimize topographic shading. Given the smaller size of the study areas compared to the Landsat footprint and the acquisition of multiple satellite images per month or year, an 80% cloud cover filter is unlikely to significantly impact the data extracted from study areas because image composition was performed per month or year. After terrain correction, the images were harmonized [39], classified by year and month, and then composited separately yearly and monthly, resulting in datasets of yearly and monthly median multispectral images. Annual values were employed to assess long-term trends and changes in vegetation conditions. Monthly values were used to evaluate seasonal trends in vegetation conditions.

3.2. Computing Vegetation Indices and Deriving Recovery Trends

This study analyzed vegetation recovery in landslide areas using two commonly used vegetation indices, the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI). EVI is preferred for its improved vegetation monitoring capability [40]. On the other hand, NDVI is chosen because it provides reliable vegetation cover estimates across various vegetation types [41] and different soil and lithological backgrounds [42]. Its ability to normalize minimizes the effect of background variations [43]. Furthermore, NDVI is frequently preferred for long-term observations compared to the Soil Adjusted Vegetation Indices (SAVI), which are more suitable for sparsely vegetated areas [44]. Although NDVI has drawbacks such as topographic illumination and shading effects issues [18,45], terrain correction during the preprocessing of Landsat time series data can generally address these limitations caused by the topographic relief [46]. To calculate average annual vegetation index values after rainfall-induced landslides, the study period was set from 1985 to 2011 for Obara and from 2011 to 2020 for Shobara. The study periods were chosen to ensure complete data coverage during the summer months of each year (June, July, and August), thus providing an accurate representation of the annual vegetation index. Annual NDVI and EVI values were obtained using linear regression for Obara time series data and log-linear regression for Shobara, reflecting the distinct vegetation recovery patterns in each study area. These regression methods were chosen based on approaches used in previous studies [47,48] to provide a simplified representation of the vegetation recovery patterns observed in each study area.

NDVI is calculated as

$$NDVI=\frac{NIR-RED}{NIR+RED}$$

(1)

where NIR is the near-infrared band intensity and RED is the red band intensity.

EVI corrects for some atmospheric conditions and canopy background noise [15]. It is calculated as

$$EVI=2.5\times \frac{(NIR-RED)}{(NIR+6RED-7.5BLUE+1)}$$

(2)

where BLUE is the blue band intensity, the value of 2.5 represents the gain factor, 6 and 7.5 represent the coefficients of the aerosol resistance term, and adding the value of 1 adjusts the canopy background [49]. All parameters were defined according to the Landsat EVI description from USGS [50].

The differences in $d_{NDVI}$ and $d_{EVI}$ calculated from two images after landslides were used to assess the net change in vegetation index value for vegetation change distribution calculation and further statistical analysis:

$$d_{NDVI}={NDVI}_{year2}-{NDVI}_{year1}$$

(3)

$$d_{EVI}={EVI}_{year2}-{EVI}_{year1}$$

(4)

In this study, $year1$ and $year2$ represent the start and end years of the study period for each area. We used linear regression for Obara and log-linear regression for Shobara to determine the values of ${NDVI}_{year1}$, ${NDVI}_{year2}$, ${EVI}_{year1}$, and ${EVI}_{year2}$. These regression models were applied to the NDVI and EVI time series data, enabling us to derive the indicators for the start and end years of the study period.

3.3. Statistical Analysis

To investigate the relationship between environmental variables and vegetation recovery, locations of damaged areas, Landsat time series, and topographic variables were converted into a vector shapefile layer and georeferenced to the EPSG:4326 coordinate system using GEE. Then, we performed multivariate analysis, including analysis of variance (ANOVA) and Tukey Honestly Significant Difference (Tukey HSD) tests. One-way, two-way, and three-way ANOVAs were conducted to understand environmental variables’ effects on vegetation recovery comprehensively. One-way ANOVA reflects the impact of each environmental variable on vegetation recovery. Two-way ANOVA shows the combined effects of two environmental variables and their interactions. Three-way ANOVA allows us to investigate the simultaneous effects of three environmental variables and their interactions. In addition, Tukey HSD [51] analysis was conducted following ANOVA to understand which specific groups of environmental variables differ by determining the honest significant differences between group means. Furthermore, Spearman’s correlation coefficient [52] was employed to assess the positive or negative relationships between the identified significantly different groups and vegetation indices.

Categorizing the data for conducting ANOVA is necessary to determine if there are statistically significant differences between the categories. The categories presented in

Table 2. Summary of environmental variables’ categories used for statistical analysis.

Topographic Variables		Categories
Obara	Shobara
E (m)
0–80	90–190	Lower level
80–160	190–290	Lower-middle level
160–240	290–390	Upper-middle level
240–330	390–500	Upper level
S (degrees)
0–20	0–20	Gentler slopes
20–40	20–40	Gentler-moderate slopes
40–60	40–60	Steeper-moderate slopes
60–80	60–80	Steeper slopes
A (degrees)
0–45	0–45	N
45–90	45–90	NE
90–135	90–135	E
135–180	135–180	SE
180–225	180–225	S
225–270	225–270	SW
270–315	270–315	W
315–360	315–360	NW
HC (m)
−0.1–0	−0.1–0	Convergent terrain
0–0.1	0–0.1	Slightly divergent terrain
0.1–0.2	0.1–0.2	Divergent terrain
VC (m)
−0.01–0	−0.1–0	Concave terrain
0–0.02	0–0.02	Slightly convex terrain
0.02–0.04	0.02–0.1	Convex terrain
MC (m)
−0.01–0	−0.01–0	Convergent and concave terrain
0–0.1	0–0.1	Slightly divergent and convex terrain
0.1–0.2	0.1–0.3	Divergent and convex terrain
G
CG	RL	-
MG	AS	-

were defined based on the range of values observed in the data (Figure 3), with varying bin widths chosen to suit the data distribution. Elevation, slope gradient, and aspect were divided into several groups based on the range of values. Horizontal, vertical, and maximal curvature were split into three groups centered around zero and boundary values. Geology (G) was categorized based on geological type. These groups were considered as independent variables, while $d_{NDVI}$ and $d_{EVI}$ were dependent variables. A p-value of less than 0.01 was used to determine statistical significance.

3.4. Modeling Vegetation Recovery

The vegetation response was modeled using the Seasonal Autoregressive Integrated Moving Averages (SARIMA) Box–Jenkins model [53]. A comparative analysis of NDVI and EVI responses and modeling performance was conducted under different environmental conditions in the two study cases. Specifically, we considered Shobara as a case representing the early stage of vegetation recovery, while Obara represented the middle and late stages.

The SARIMA model helps model stochastic processes considering the past behavior of a time series. It can be composed of an autoregressive (AR) model, a moving average (MA) model, or a combination of both (ARMA). The ARMA model is appropriate for stationary series where values vary uniformly around the mean without trend. However, for non-stationary series, such as those represented by the vegetation indices in this study, it is necessary to achieve stationarity by identifying and removing all factors contributing to non-stationarity, resulting in an Autoregressive Integrated Moving Averages (ARIMA) model. Moreover, vegetation index time series often exhibit marked seasonality, which can be incorporated into the Seasonal ARIMA (SARIMA) model through seasonal differences. The SARIMA model is represented as

SARIMA(p,d,q)(P,D,Q)_S

(5)

where p and q are the orders of the AR and MA models, respectively, and d is the order of the difference applied to the series to achieve stationarity. P, D, and Q are the orders of the seasonal differences used concerning a seasonal period S [54]. The simplified equation for the SARIMA model is

$${\varphi }_p\left(B\right){\Phi }_p\left(B^S\right){\nabla }^d{\nabla }_S^D{Y}_t={\theta }_q\left(B\right){\Theta }_Q\left(B^S\right){\epsilon }_t$$

(6)

where $Y_t$ is the observed value of the time series in period t with a seasonal period of S, ${\epsilon }_t$ is a white noise process in period t, and $B$ is a simple or seasonal backshift operator where for each ${\epsilon }_t$: $B{\epsilon }_t={\epsilon }_{t-1}$ and $B^S{\epsilon }_t={\epsilon }_{t-S}$. In this way, the operators of the simple, seasonal AR and MA models can be written concerning B. Furthermore, the $\nabla$ operator of the simple and seasonal differences is also written in terms of B for each ${\epsilon }_t$.

Building a SARIMA model involves three main steps: identification, parameter training and test number estimation, and model assessment. In the identification step, the orders of differencing d and D were selected to achieve stationarity, and the order of the autoregressive (AR), moving average (MA), and seasonal components (P and Q) of the model were determined based on the analysis of autocorrelation functions (ACFs) and partial autocorrelation functions (PACFs) of the differenced series. After the identification step, we estimated all parameters in Equation (5) and the number of training and test data points. Three statistical indices were used to evaluate the model performance: the determination coefficient (R²), the Mean Absolute Percentage Error (MAPE), and the Root Mean Square Error (RMSE). The best model was selected based on the highest overall performance for each time series. Furthermore, we evaluated the performance of the best model to determine the vegetation indices suitable for modeling vegetation recovery trends under different environmental conditions in the two study cases.

4. Results

4.1. NDVI and EVI Responses to Vegetation Recovery

The temporal changes in NDVI and EVI exhibit a distinct seasonal pattern. High NDVI and EVI characterize the summer season from July to October, while lower NDVI and EVI were observed during the withering season from December to April (Figure 4b,d). There is a negligible difference between NDVI and EVI temporal changes in both study areas and seasonal cycles. The annual mean NDVI and EVI on the damaged slopes ranged from 0.2 to 0.6 in the Obara district (Figure 4a). In Shobara, the mean annual NDVI and EVI in the damaged areas ranged from 0.1 to 0.6 (Figure 4c). Before the heavy rainfall in 2011, the NDVI values at Shobara ranged from 0.25 to 0.45. However, in 2011, a year after the rainfall-induced landslide event, the mean NDVI and EVI decreased to 0.2 and 0.3, respectively (Figure 4c). For Obara, although satellite imagery was unavailable before 1985, the increasing trend in annual mean NDVI and EVI between 1985 and 2011 indicates continuous recovery even after landslides occurred approximately ten years earlier (Figure 4a). Aerial photographs showed rainfall-induced landslides in Obara caused vegetation loss in July 1972 (Figure 5b). Post-rainfall aerial photographs from 1972 to 1998 showed increased vegetation coverage in the damaged areas over the years, suggesting vegetation recovery (Figure 5b–f).

The regression models used to derive ${NDVI}_{year1}$, ${NDVI}_{year2}$, ${EVI}_{year1}$, and ${EVI}_{year2}$ for $d_{NDVI}$ and $d_{EVI}$ calculations demonstrated the following best fits: linear regression in Obara (NDVI: R² = 0.36; EVI: R² = 0.38) and log-linear regression in Shobara (NDVI: R² = 0.14; EVI: R² = 0. 34). Changes in vegetation indices were categorized using the threshold values of $d_{NDVI}$ and $d_{EVI}$, as described in

Table 3. Threshold values of $d_{NDVI}$ and $d_{EVI}$ and the corresponding frequency distribution in the Obara district and the Shobara district. “Poor” to “Excellent” categories indicate increasing positive changes in vegetation index values. “Very poor” means a negative change in the vegetation index value.

Category	Obara			Shobara
	$d_{NDVI}$ and $d_{EVI}$	$\mathrm{Distribution} (\%) \mathrm{of} d_{NDVI}$	$\mathrm{Distribution} (\%) \mathrm{of} d_{EVI}$	$d_{NDVI}$ and $d_{EVI}$	$\mathrm{Distribution} (\%) \mathrm{of} d_{NDVI}$	$\mathrm{Distribution} (\%) \mathrm{of} d_{EVI}$
Excellent	>0.3	33.8	5.1	>0.2	18.6	1.9
Very good	0.2–0.3	28.1	6.9	0.1–0.2	39.4	12.2
Good	0.1–0.2	21.3	37.8	0.05–0.1	22.9	19.9
Poor	0–0.1	16.4	31.8	0–0.05	19.1	29.3
Very poor	<0	0.4	18.4	<0	0	36.7

. The same threshold values for $d_{NDVI}$ and $d_{EVI}$ were used for each area to compare the differences in NDVI and EVI changes over time after landslides. Overall, $d_{NDVI}$ exhibited greater changes than $d_{EVI}$ in both study areas.

4.2. Statistical Analyses on Environmental Variables Influencing Vegetation Recovery

The results of the one-way ANOVA analysis (

Table 4. Statistical analysis results of one-way ANOVA.

Variable	Obara		Shobara
	$d_{NDVI}$	$d_{EVI}$	$d_{NDVI}$	$d_{EVI}$
E	0.18	0.05	0.02	0.06
S	0.02	0.03	0.02	0.02
A	0.03	0.02	0.01 *	0.009 *
HC	0.002 *	0.06	0.07	0.19
VC	0.049	0.003 *	0.085	0.271
MC	0.05	0.003 *	0.044	0.128
G	0.005 *	0.02	0.07	0.14

* indicates statistical significance equal to or lower than 1%.

) reveal that in Obara, horizontal curvature and geological type significantly affect $d_{NDVI}$, whereas vertical curvature and maximum curvature affect $d_{EVI}$. In contrast, in Shobara, the aspect significantly affects both $d_{NDVI}$ and $d_{EVI}$. The two-way ANOVA results (

Table 5. Statistical analysis results of two-way ANOVA.

Variable1	Variable2	Obara		Shobara
		$d_{NDVI}$	$d_{EVI}$	$d_{NDVI}$	$d_{EVI}$
E	S	0.1	0.1	0.139	0.1
	A	0.001 *	0.006 *	0.028	0.559
	HC	0.089	0.017	0.084	0.076
	VC	0.001 *	0.001 *	0.075	0.001 *
	MC	0.045	0.001 *	0.01 *	0.071
	G	0.023	0.2	0.078	0.34
S	A	0.001 *	0.1	0.01 *	0.006 *
	HC	0.081	0.001 *	0.103	0.766
	VC	0.104	0.301	0.01 *	0.001 *
	MC	0.503	0.143	0.023	0.607
	G	0.100	0.100	0.56	0.67
A	HC	0.309	0.592	0.305	0.107
	VC	0.11	0.022	0.337	0.09
	MC	0.001 *	0.001 *	0.003 *	0.640
	G	0.02	0.05	0.56	0.41
HC	VC	0.488	0.107	0.65	0.572
	MC	0.100	0.200	0.01 *	0.021
	G	0.146	0.166	0.27	0.34
VC	MC	0.002 *	0.236	0.2	0.276
	G	0.2	0.14	0.1	0.1
MC	G	0.18	0.20	0.5	0.67

* indicates statistical significance equal to or lower than 1%.

) indicate that in Obara, the combination of maximum curvature or aspect with elevation, slope gradient, and vertical curvature significantly affects both $d_{NDVI}$ and $d_{EVI}$. Similarly, in Shobara, the combination of maximum curvature or aspect with elevation, slope gradient, and vertical curvature significantly affects $d_{NDVI}$ and $d_{EVI}$. Interestingly, the aspect effect is significant in vegetation recovery from one-way and two-way ANOVA. Furthermore, geological type does not show effective results in the two-way ANOVA analysis in either study area, indicating its independence from other topographic variables. The three-way ANOVA analysis (

Table 6. Statistical analysis results of three-way ANOVA.

Variable1	Variable2	Variable3	Obara		Shobara
			$d_{NDVI}$	$d_{EVI}$	$d_{NDVI}$	$d_{EVI}$
E	S	A	0.005 *	0.001 *	0.001 *	0.005 *
		HC	0.01 *	0.001 *	0.007 *	0.057
		VC	0.01 *	0.001 *	0.068	0.018
		MC	0.001 *	0.001 *	0.685	0.117
		G	0.278	0.346	0.510	0.690
S	A	HC	0.01 *	0.003 *	0.019	0.02
		VC	0.007 *	0.001 *	0.004 *	0.002
		MC	0.001 *	0.001 *	0.009 *	0.326
		G	0.402	0.153	0.013	0.307
A	HC	VC	0.002 *	0.003 *	0.313	0.025
		MC	0.001 *	0.001 *	0.718	0.654
		G	0.02	0.05	0.56	0.41
HC	VC	MC	0.01 *	0.01 *	0.001 *	0.029
		G	0.12	0.24	0.267	0.646
VC	MC	E	0.01 *	0.001 *	0.966	0.803
		G	0.584	0.301	0.560	0.170
MC	G	E	0.100	0.200	0.02	0.021

* indicates statistical significance equal to or lower than 1%.

) reveals that in Obara, almost all the combinations significantly affect both $d_{NDVI}$ and $d_{EVI}$. In Shobara, the significant combinations identified from two-way ANOVA results (maximum curvature or slope gradient and aspect) continued to show substantial effects when combined with other factors, such as elevation and vertical curvature, on $d_{NDVI}$. Only the combination of elevation, slope gradient, and aspect significantly affects $d_{EVI}$. Notably, the number of significant results is larger for $d_{NDVI}$ than for $d_{EVI}$ in both the two-way and three-way ANOVA analyses conducted in both study areas.

As ANOVA only provides a general test of the differences among variables, it is necessary to perform Tukey HSD tests to determine which specific categories are significantly different. It is important to note that the Tukey HSD results may differ from the ANOVA results, as a variable that shows a significant difference in ANOVA may not yield a significant difference in specific categories according to Tukey HSD. While ANOVA examines relationships among all levels, Tukey HSD compares individual levels.

Table 7. Summary of environmental variables’ significantly different group(s) derived from Tukey HSD analysis. The corresponding effects, derived from Spearman’s correlation coefficient, are listed in the right column for each group. A negative effect is denoted by ‘−’, a positive effect by ‘+’, and empty spaces indicate variables that did not yield significantly different groups or distinct effects.

Variables	Obara				Shobara
	NDVI		EVI		NDVI		EVI
E (m)	80–160, 240–330	−	80–160, 160–240	−	190–290, 290–390	−	120–250, 250–300	−
			0–80	+
S (degrees)	0–20	+	0–20, 20–40	+	0–20	+	0–20	+
	60–80	−			60–80	−
A (degrees)			0–45	+	315–360	+	270–315	+
			135–180	−
VC (m)			−0.01–0	+
HC (m)	0–0.1	−	0–0.1	−	−0.1–0	−
MC (m)	−0.01–0	+	−0.01–0	+
G

presents the Tukey HSD analysis results, categorizing the significantly different groups of environmental variables based on their positive or negative impact on vegetation recovery. The findings in Table 7 suggest that the middle elevation levels in the local damaged areas differ considerably from other groups in both study areas. For slope gradient, the flat and steep slopes differ substantially from other groups in both study areas. For aspect, the north and southeast groups in Obara and the west and northwest groups in Shobara show significant differences. For horizontal curvature, the group with a positive and higher value shows significant variations in Obara, whereas the group that approaches a zero value shows substantially different results in Shobara. Finally, the groups with negative vertical and maximum curvature values show significant differences.

4.3. Trend Analyses and SARIMA Modeling of Environmental Variables

This study analyzed the effects of various environmental variables on vegetation regrowth, categorizing them as positive, neutral, and negative conditions represented by ‘+’, ‘blank’, and ‘−’, respectively, as shown in Table 7. The results are illustrated in Figure 6, where every value represents an average of the pixels in a particular group: the green dotted line represents vegetation recovery under neutral conditions, and the blue and purple lines represent vegetation recovery under positive and negative situations, respectively. The effects of different combinations of environmental variables on NDVI were not evident in Obara (Figure 6a). However, the impact of various combinations of environmental variables on EVI was more apparent than on NDVI, especially under the positive environmental condition (blue line in Figure 6b). In Shobara, the impact of different environmental factors was notable for NDVI and EVI time series trends after landslides. Specifically, the effect was more pronounced on the NDVI time series trend, whereas the unfavorable condition had a distinct impact on the EVI time series trend, different from the neutral and positive conditions (Figure 6c,d).

SARIMA was used to model the NDVI and EVI time series trends for the two study areas based on different combinations of environmental variables. The optimal order combination for each model was determined by feeding all the time series data into the SARIMA model and referring to the Akaike Information Criterion (AIC) value [55], commonly used to assess the relative quality of different statistical models. Three statistical indices, R², MAPE, and RMSE, were used to evaluate the performance of each model. The results are presented in

Table 8. Optimal order combination of the SARIMA model and the corresponding performance indices for each environmental condition in the Obara and Shobara districts. The subscript number in the Order line indicates a seasonal period S in Equation (5).

	NDVI			EVI
	Neutral	Positive	Negative	Neutral	Positive	Negative
Obara
Order	(1,1,0) (0,0,0)12	(1,1,0) (0,0,0)12	(1,1,0) (0,0,0)12	(1,1,0) (0,0,0)12	(1,1,0) (0,0,0)12	(1,1,0) (0,0,0)12
Train, Test	(238, 80)	(238, 80)	(238, 80)	(238, 80)	(238, 80)	(238, 80)
MAPE (%)	0.90	0.65	0.63	1.10	0.97	1.12
RMSE	0.0095	0.0070	0.0070	0.0114	0.0125	0.0118
R2	0.9738	0.9867	0.9872	0.9615	0.9650	0.9614
Shobara
Order	(1,0,1) (1,0,1)12	(2,1,1) (1,0,1)12	(2,1,1) (1,0,1)12	(1,0,1) (1,0,1)12	(2,1,1) (1,0,1)12	(1,0,1) (1,0,1)12
Train, Test	(155, 97)	(155, 97)	(155, 97)	(155, 97)	(155, 97)	(155, 97)
MAPE (%)	0.62	0.77	0.74	0.56	0.57	0.35
RMSE	0.0061	0.0080	0.0071	0.0068	0.0068	0.0041
R2	0.9739	0.9584	0.9600	0.9677	0.9678	0.9865

. Furthermore, we plotted the modeled values (orange line) of the test dataset for each model and compared them with the observed values (blue line) using 99% confidence intervals (gray area) in Figure 7 and Figure 8.

5. Discussion

5.1. Prolonged Vegetation Recovery

The NDVI and EVI trend analysis revealed that the vegetation recovery processes in the study areas exhibit a protracted duration, potentially extending beyond ten years (Figure 4 and Figure 5). This timeframe is longer than the average revegetation timescale documented in previous studies on Japan’s earthquake- and rainfall-induced landslide-affected areas [56,57,58]. Determining complete vegetation recovery in this study involves a comprehensive approach combining aerial photograph interpretation with the long-term stability of vegetation index trends. Estimating the time required for the vegetation recovery processes in Obara can be inferred by synthesizing two critical aspects of the findings. Firstly, historical aerial photograph analysis reveals that the vegetation coverage within the damaged regions remained incomplete by 1982 (Figure 5b,e). Secondly, trend analysis of NDVI and EVI calculated from Landsat imagery from 1985 to 2011 shows that both NDVI and EVI exhibited a consistent upward trend during this timeframe, influenced by seasonal fluctuations (Figure 4a,b). For Shobara, from 2011 to 2020, the observed trend revealed an upward trajectory of seasonal fluctuations in NDVI and EVI values without reaching a stable plateau of no further increase (Figure 5c,d).

The slower vegetation recovery observed in these study areas can be attributed to several factors. One is the severity of the landslide disturbance compared to the other areas. The distribution of damaged areas in both study areas is dense and extensive, and previous studies have shown a positive correlation between the disturbance magnitude and the time required for vegetation to recover fully [59,60]. The density of landslide-affected areas is around 10 and 15 patches km⁻² for Obara and Shobara, respectively, which is denser than 8–10 patches km⁻² for rainfall-induced landslides in Chiba Prefecture [56] and 5.1 patches km⁻² in the Kiso Mountains [57], central Japan. In addition, artificial vegetation was carried out on severely damaged slopes and landslide patches in the south Kiso Mountains [57]. By contrast, both study areas were left to recover naturally, extending vegetation recovery time. Therefore, the prolonged vegetation recovery can be attributed to the severity and density of the landslide disturbance and the absence of artificial vegetation interventions.

5.2. Environmental Variables Controlling Vegetation Recovery

Our analysis revealed that topographic variables (elevation, slope gradient, aspect, and curvature) significantly influenced vegetation recovery and that geological type exerted an independent influence. Previous studies have highlighted the independent impact of geological type on vegetation dynamics [61,62,63]. Different geological contexts form soil with different physical, chemical, and mineralogical properties [64,65]. Therefore, geology’s independent effect on vegetation recovery can be attributed to its influence on soil properties. However, our dataset does not enable the differentiation of soil properties, although NDVI may indicate soil reflectance to some extent [66].

For topographic variables, in Obara, rapid vegetation recovery was observed in concave and convergent terrain according to NDVI and in non-steep, low-elevation, north-facing terrain according to EVI. In Shobara, fast vegetation recovery was observed on gentler and northwest-facing slopes in terms of NDVI, and on gentler and west-facing slopes in terms of EVI. Rapid vegetation recovery in the north- to west-facing and non-steep regions may reflect higher soil moisture because of lower solar radiation (Figure 3c) [67]. Additionally, steep slope gradients may limit vegetation recovery due to limited water availability and difficulties in seed settlement and root anchoring [68]. While the elevation of the study areas did not induce significant variation in water or thermal conditions, the ANOVA results suggest that when combined with aspect or curvature, elevation may interact with vegetation recovery. Curvature shows substantial results from the three-way ANOVA (Table 6), especially in groups with slope gradients. Slightly convergent and concave terrain enhances vegetation recovery (Table 7), likely because concave depressions accumulate more water and nutrients [69,70].

Unfavorable conditions for vegetation recovery in the study areas, such as middle to upper elevations, steeper and southwest-facing slopes, and slight-divergent terrain, were widespread in the landslide-affected areas (Figure 3) and contributed to the delayed vegetation recovery observed in Section 5.1. The cumulative effect of these conditions exacerbates challenges for revegetation, including limited resources, heightened erosion risk, and reduced microclimatic suitability.

5.3. Comparisons of NDVI and EVI

This study utilized NDVI and EVI for monitoring vegetation recovery in damaged areas. NDVI is sensitive to atmospheric effects, such as cloud-induced variations [71]. Although the Landsat cloud mask helps mitigate cloud properties in per-pixel reflectance, it does not entirely eliminate bias in time series analysis [72]. On the other hand, EVI tends to show higher values in sparsely vegetated areas than NDVI [71], as NDVI is sensitive to the spectral effects of soil moisture and texture in sparsely vegetated areas [73,74]. Therefore, the high frequency of cloud cover and low vegetation coverage in Shobara has led to more pronounced differences between EVI and NDVI than in Obara (Figure 4). The greater variation in the distribution of NDVI changes (Table 3) may also reflect the high sensitivity of NDVI to soil moisture and texture [73,74].

The Shobara and Obara districts represent the early and middle-late stages of vegetation recovery, respectively. As vegetation recovery progresses, ecosystems exhibit increased resilience, diversity, and stability, leading to a decreased sensitivity to variations in environmental conditions [75] and a more consistent trend in the NDVI (Figure 6a). Although EVI was proposed to reduce atmospheric and soil background noise simultaneously, it was found to be sensitive to different environmental conditions in the two study areas (Figure 6). This may result from the soil adjustment factor in EVI that makes it more susceptible to topographic effects than NDVI [18]. Table 8 indicates that EVI outperforms NDVI in Shobara in modeling vegetation recovery trends under different environmental conditions, as the smaller MAPE value suggests. EVI incorporates the blue band, which minimizes atmospheric influences and improves sensitivity to changes in vegetation structure and density [18]. This makes EVI particularly suitable for capturing changes in early-stage vegetation conditions characterized by low structure and density, visible soil, and sparse vegetation cover. On the other hand, in Obara, NDVI demonstrates superior performance, as indicated by the lower MAPE value and higher R² value. This observation suggests that NDVI is better suited for modeling the middle and late stages of vegetation recovery within 10–25 years after the landslide. Previous studies have also reported the effectiveness of NDVI in assessing long-term vegetation dynamics in landslide-affected areas under different environmental conditions [76,77]. In addition, the Landsat time series we obtained for Obara has less cloud cover than the Shobara time series, leading to the superior performance of NDVI in Obara. However, cloud cover and its impact on NDVI do not undermine the overall effectiveness of NDVI in modeling longer-term vegetation recovery. It highlights the need to consider environmental conditions when interpreting NDVI data. Therefore, choosing appropriate and even different vegetation indices is necessary for modeling specific stages of vegetation recovery.

6. Conclusions

This study investigated vegetation recovery following rainfall-induced landslides in two areas of Japan, the Shobara and Obara districts, using Landsat data to derive NDVI and EVI time series. We found that vegetation recovery was gradual and consistent, taking over a decade without reaching a near-stable state, exceeding the previously reported vegetation recovery time for landslide-affected areas in Japan. This was caused by the severity of landslide disturbance, lack of restoration work, and the widespread occurrence of unfavorable environmental conditions. Statistical analysis of environmental variables showed that elevation, slope gradient, aspect, and curvature significantly correlated with NDVI and EVI responses in both study areas. The SARIMA model performed well in modeling NDVI and EVI under different environmental conditions. We also found that NDVI is more suitable for modeling the middle and late stages of vegetation recovery within 10–25 years after the landslide for various environmental conditions. In contrast, EVI is more suitable for modeling the early stages within 10 years after the landslides. NDVI provides insights into the long-term health and maturity of the recovering vegetation, while EVI helps capture the initial vegetation growth patterns. These results successfully addressed the primary motivation of investigating vegetation recovery and identifying suitable vegetation indices to understand post-landslide vegetation dynamics while considering the complex relationship between topography and vegetation over time. The study’s findings allow for targeted and nuanced monitoring and assessment of vegetation dynamics, enabling effective management and decision making in post-landslide ecosystems.

This study acknowledges several limitations. Firstly, Obara’s lack of certain data restricted the analysis to satellite imagery and aerial photographs, limiting the comprehensive assessment. Future studies should consider incorporating higher-resolution imagery or Light Detection and Ranging (LiDAR) data to enhance correlations with vegetation characteristics. To enhance data quality, integrating Sentinel data with Landsat imagery [78,79] can improve the accuracy of future vegetation dynamics estimates. Additionally, the contrasting results between NDVI and EVI in modeling vegetation recovery under different conditions require further investigation through additional data collection, such as soil moisture and local climate data, to better understand the underlying mechanisms. Meanwhile, further research is needed to investigate the reciprocal influences of geomorphic processes on vegetation dynamics by simulating topographic changes, enhancing the comprehension of the intricate relationship between vegetation recovery and topography. Finally, although regression models have helped to capture general trends in vegetation recovery in Obara and Shobara, future studies could consider exploring advanced approaches, such as long short-term memory (LSTM) networks, to better understand the complex nature of vegetation recovery in these regions.