The Role of Citrus Groves in Rainfall-Triggered Landslide Hazards in Uwajima, Japan

: Landslides often cause deaths and severe economic losses. In general, forests play an important role in reducing landslide probability because of the stabilizing effect of the tree roots. Although fruit groves consist of trees, which are similar to forests, practical land management, such as the frequent trampling of ﬁelds by laborers and compression of the terrain, may cause such land to become prone to landslides compared with forests. Fruit groves are widely distributed in hilly regions, but few studies have examined their role in landslide initiation. This study aims at ﬁlling this gap evaluating the predisposing and triggering conditions for rainfall-triggering landslides in part of Uwajima City, Japan. A large number of landslides occurred due to a heavy rainfall event in July 2018, where citrus groves occupied about 50% of the study area. In this study, we combined geodata with a regression model to assess the landslide hazard of fruit groves in hilly regions. We developed maps for ﬁve conditioning factors: slope gradient, slope aspect, normalized difference vegetation index (NDVI), land use, and geology. Based on these ﬁve maps and a landslide inventory map, we found that the landslide area density in citrus groves was larger than in forests for the categories of slope gradient, slope aspect, NDVI, and geology. Ten logistic regression models along with different rainfall indices (i.e., 1-h, 3-h, 12-h, 24-h maximum rainfall and total rainfall) and different land use (forests or citrus groves) in addition to the other four conditioning factors were produced. The result revealed that “citrus grove” was a signiﬁcant factor with a positive coefﬁcient for all models, whereas “forest” was a negative coefﬁcient. These results suggest that citrus groves have a higher probability of landslide initiation than forests in this study area. Similar studies targeting different sites with various types of fruit groves and several rainfall events are crucial to generalize the analysis of landslide hazard in fruit groves.


Introduction
Landslides are frightening and destructive natural disasters. The Global Fatal Landslide Database recorded 55,997 fatalities from 4862 landslide events from 2004 to 2016 [1], and the number of landslides has increased over time [2,3]. For example, the Centre for Research on the Epidemiology of Disaster and United Nations Office for Disaster Risk Reduction [4] reported a 48% increase in landslide events from 1999 to 2019. Although several approaches for landslide mitigation exist, understanding the hazards in areas prone to landslides is a fundamental approach.
Geology, morphology, vegetation, and human activities are factors affecting landslide hazards [5,6], and land cover is a key determining factor. Forests and trees play important roles in reducing landslide risk mainly due to shear strength enhancement by roots on soil particles of slip surfaces [7,8]. A comparison of landslide distribution among land use types muddy turbidite, and turbidite soil. More detailed topography, geology, and land use of the study area are shown in the results section. Figure 1. Location of the study area (in Uwajima City, southwest part of Ehime prefecture) and landslide distribution in July 2018. The height of elevation means height above sea level.

The Rainfall Event in July 2018
Heavy rainfall events occurred from 1 to 8 July 2018. At the Uwajima station, the total rainfall was 578.0 mm; the maximum hourly rainfall observed from 06:00 to 07:00 on 7 July was 49.0 mm ( Figure 2). The rainfall event induced numerous landslides in several cities in Ehime prefecture, particularly in Uwajima City. Uwajima City suffered the greatest landslide disaster on 7 July 2018; in total, 51 debris flows and 102 steep slope failures affected houses and public facilities [26].

The Rainfall Event in July 2018
Heavy rainfall events occurred from 1 to 8 July 2018. At the Uwajima station, the total rainfall was 578.0 mm; the maximum hourly rainfall observed from 06:00 to 07:00 on 7 July was 49.0 mm ( Figure 2). The rainfall event induced numerous landslides in several cities in Ehime prefecture, particularly in Uwajima City. Uwajima City suffered the greatest landslide disaster on 7 July 2018; in total, 51 debris flows and 102 steep slope failures affected houses and public facilities [26].
Water 2022, 14, x FOR PEER REVIEW 3 of 16 muddy turbidite, and turbidite soil. More detailed topography, geology, and land use of the study area are shown in the results section. Figure 1. Location of the study area (in Uwajima City, southwest part of Ehime prefecture) and landslide distribution in July 2018. The height of elevation means height above sea level.

The Rainfall Event in July 2018
Heavy rainfall events occurred from 1 to 8 July 2018. At the Uwajima station, the total rainfall was 578.0 mm; the maximum hourly rainfall observed from 06:00 to 07:00 on 7 July was 49.0 mm ( Figure 2). The rainfall event induced numerous landslides in several cities in Ehime prefecture, particularly in Uwajima City. Uwajima City suffered the greatest landslide disaster on 7 July 2018; in total, 51 debris flows and 102 steep slope failures affected houses and public facilities [26].

Data Preparation
We used a LIM, distributions of slope gradient, slope aspect, NDVI, land use, and geology as conditioning factors, and rainfall for analysis. Data of the conditioning factors were generally selected based on the abundance, accessibility, region characteristics, and research purposes [27,28]. Because we aimed to evaluate the effects of citrus groves on landslide initiation, we considered land use as a primary conditioning factor. We also included morphological factors such as slope gradient and aspect, as well as NDVI and geology, as conditioning factors. The number of conditioning factors may not be sufficient for a landslide susceptibility study. However, they were appropriate for the hazard assessment in the area according to the study purposes. Table 1 summarizes the datasets used. We used a LIM provided by the Geospatial Information Authority of Japan (GSI) [29], which was produced using aerial photographs taken both before and shortly after the rainfall event. The LIM included 162 landslides and occupied 0.53% of the study area ( Figure 1). The slope gradient and aspect were calculated using a digital elevation model (DEM). ASTER Global DEM (GDEM) was used for the calculations [30], and NDVI was calculated from the reflectance of the near-infrared wavelength (NIR) and red wavelength (RED): Herein, we used the NIR and RED of Landsat 8, collection 2, level-2, measured on 20 April 2018 [31], which was a clear day just before the landslide events. Landsat series collection 2, level-2, was generated from collection 2, level-1 inputs. The level-1 consisted of quantized and calibrated scaled digital number values representing multispectral image data, whereas level-2 data consisted of pixel values estimating the spectral reflectance measured at the Earth's surface without atmospheric scattering or absorption effects [32]. Land use and geological maps were provided by the Ministry of Land, Infrastructure, Transport, and Tourism in vector format (polygon coverage) [33]; these were converted into raster data with a resolution of 100 m. The land use map included six categories (i.e., "forest," "other agricultural land (land where vegetables and fruits are cultivated)," "field (agriculture)," "building site," "road," and "water"). Within the study area, the category "other agricultural land" mainly included citrus groves. Although the original classification was used when examining the landslide density for each category, these were further classified as "forest," "citrus grove," and "the others", which represented the field, building site, road, and water categories. The purpose of the new classification was to highlight the difference in landslide initiation between "forest" and "citrus grove." The geological map in this area included five geological units: "sandy turbidite," "muddy turbidite," "turbidite," "limestone block," and "chert block." For rainfall analysis, rainfall distributions based on the X-band MP radar network were used [34]. We downloaded the data that were recorded every 10 min on 1-8 July, and we calculated the 1-h, 3-h, 12-h, and 24-h maximum rainfall and total rainfall during the study period.

Analysis
First, maps for the one triggering and five conditioning factors were developed. For rainfall, maps for the 1-h, 3-h, 12-h, and 24-h maximum rainfall and the total rainfall during the study period were developed. Second, we calculated the landslide area proportional to the entire study area, and the LAD for each category was calculated. The LAD was calculated as the ratio of the area with landslides to the total area of each category [35,36]. We established eight categories for slope aspect, and six categories with the same data range were set for the slope gradient and NDVI. For the five rainfall indices, four categories with the same data range were set. Third, the LAD was compared in the new land use categories for each category of slope gradient, slope aspect, NDVI, and geology. Finally, multivariate analyses were performed to comprehensively examine the effects of each factor on landslide initiation. For this analysis, we applied logistic regression although random forest has also been frequently used recently. However, according to a previous review [37], logistic regression has been the most commonly used for landslide susceptibility modeling. The logistic regression can be expressed as follows [38]: i is the number of independent variables, and n (= 6 in this study) is the number of independent variables. Landslides and non-landslide points are required to predict the probability of landslide occurrences in the logistic regression model and evaluating the model accuracy by trueor false-positive and negative rates. Landslide points at pixels with the highest elevation for each landslide were set [39,40]. Non-landslide points were randomly selected in areas without landslides using the "random points inside polygons" function in QGIS. We set the minimum distance to 10 m, and the observed mean distance was 262.3 m as a result. Following the previous literature [38,41], the number of non-landslide and landslide points was the same. Both landslide and non-landslide points were divided into training sets for model development and testing sets for model validation. According to the suggestions in previous studies [38,[42][43][44], the ratios of the training and testing sets were 80% and 20%, respectively.
The slope gradient, slope aspect, NDVI, land use, geology, and rainfall were considered as independent variables. Dummy variables were used for categorical factors of land use and geology [45] (see Appendix A Tables A1 and A2). We developed ten models with different rainfall indices, namely 1-h, 3-h, 12-h, and 24-h maximum rainfall and total rainfall and different land use categories, namely "forest" and "citrus grove" (Appendix A  Table A3). To investigate the linearity among independent variables [46], the variance inflation factor (VIF) was used, which is widely used for multicollinearity analysis. VIF indicates the effect of the determination coefficients on the variance of the estimated regression coefficients [45]. A high VIF value indicates that the model is not feasible due to multicollinearity among independent variables [45,46]. We set an acceptable VIF of 2.5, because several studies stated that VIF ≥ 2.5 indicates considerable collinearity [47,48], and VIF > 5 or 10 indicates a serious collinearity problem. The two land use types were separately calculated because the VIF values were >2.5 without separation.
The independent variables for each model were reduced stepwise based on the akaike information criterion (AIC). A stepwise regression automatically assigned the best subset of independent variables throughout several stages [49]. Each conditioning factor was evaluated one-by-one at every step using t statistics to determine the variable coefficients. To ensure acceptable model accuracy, the final models were evaluated using the area under the curve (AUC), which was obtained by constructing a receiver operating characteristic (ROC) curve comparing true-and false-positive rates, an F1 score by Equation (3) [50], and the overall accuracy by Equation (4) [51]. The all-accuracy score takes a value between 0% and 100%, and a higher score indicates higher model accuracy [46,51,52]. All geoprocessing was performed with QGIS (ver. 3.10), and logistic regression analysis was performed using R (ver. 4.1.0) programming language.
Overall Accuracy = TP + TN where TP is the true positive, TN is the true negative, FP is the false positive, and FN is the false negative.
In the evaluations, we set the threshold of 0.5 between the occurrence and nonoccurrence of landslides. Figure 3 shows the maps of the five conditioning factors. The proportional area and LAD for each category in all conditioning factors are shown in Figure 4. The center of the LAD distribution was a slope gradient of 16-24 • . The LAD was smaller when the slope gradient was smaller or larger than 16-24 • (Figure 4a). There were no considerable differences in the proportional area among the eight slope aspects. The LADs in the southern-and northern-facing slopes tended to be larger than those in the western-and eastern-facing slopes (Figure 4b). The NDVI for the study area was 0.05 and 0.83; over 70% of the area had an NDVI of >0.68 (Figure 4c). For land use, "citrus grove" was the dominant (52.74%), followed by "forest" (35.80%) and "building site" (7.31%) (Figure 4d). The highest LAD was found on "citrus grove" (0.41%), followed by "forest" (0.11%). The LADs of the other four land use types were 0.01% (Figure 4d). Limestones and chert blocks occupied <2% (Figure 4e). The LAD of sandy turbidites was more than twice as large as that of muddy turbidites and turbidite.

Landslide Area Density Based on Rainfall Indices
Areas with a high 1-h maximum rainfall (35-63 mm) were distributed in the western part of the study area (Figure 5a). In contrast, areas with a high 3-h maximum rainfall were distributed in the northern part of the study area (Figure 5b). The distributions of the 12h ( Figure 5c) and 24-h ( Figure 5d) maximum rainfall were similar to that of the total rainfall (Figure 5e); these were higher in the northeast to southwest areas and lower in the southeastern and southwestern areas. The LAD increased with an increase in the 1-h maximum rainfall (Figure 5f). The pattern of the 3-h maximum rainfall was different from that of the 1-h maximum rainfall. The middle range had the highest LAD (Figure 5g). The LADs for the 12-h (Figure 5h) and 24-h (Figure 5i) maximum rainfall and the total rainfall ( Figure 5j) were comparable due to the similarity in rainfall distributions. The LAD steeply increased from the first to the second categories, and it gradually increased from the second to the fourth categories.

Landslide Area Density Based on Rainfall Indices
Areas with a high 1-h maximum rainfall (35-63 mm) were distributed in the western part of the study area (Figure 5a). In contrast, areas with a high 3-h maximum rainfall were distributed in the northern part of the study area (Figure 5b). The distributions of the 12-h ( Figure 5c) and 24-h ( Figure 5d) maximum rainfall were similar to that of the total rainfall (Figure 5e); these were higher in the northeast to southwest areas and lower in the southeastern and southwestern areas. The LAD increased with an increase in the 1-h maximum rainfall (Figure 5f). The pattern of the 3-h maximum rainfall was different from that of the 1-h maximum rainfall. The middle range had the highest LAD (Figure 5g). The LADs for the 12-h (Figure 5h) and 24-h (Figure 5i) maximum rainfall and the total rainfall ( Figure 5j) were comparable due to the similarity in rainfall distributions. The LAD steeply increased from the first to the second categories, and it gradually increased from the second to the fourth categories.   Figure 6 compared the LADs of "forest," "citrus grove," and "the others" for each category of slope gradient, slope aspect, NDVI, and geology. For the slope gradient (Figure 6a), the LAD for citrus groves was larger than that for forests on gentle slopes (<24 • ); in contrast, it was an inconsiderable difference on steep slopes (>32 • ). For the slope aspect (Figure 6b), NDVI (Figure 6c), and geology (Figure 6d), the LAD for "citrus grove" was larger than that for "forest" in all categories. Figure 6 compared the LADs of "forest," "citrus grove," and "the others" for each category of slope gradient, slope aspect, NDVI, and geology. For the slope gradient (Figure 6a), the LAD for citrus groves was larger than that for forests on gentle slopes (<24°); in contrast, it was an inconsiderable difference on steep slopes (>32°). For the slope aspect (Figure 6b), NDVI (Figure 6c), and geology (Figure 6d), the LAD for "citrus grove" was larger than that for "forest" in all categories.

Significant Effects of Landslide Conditioning Factors
All the independent variables for the ten models had a VIF of <2 (Appendix A Tables A4 and A5), indicating that all the variables could be simulated together. Table 2 summarized the selected factors; after the stepwise processes of logistic regression analysis, NDVI, "forest," and "citrus grove" were selected as significant factors. Because rainfall was not selected for models 2-5 and 7-10, the same result was obtained for the eight models. For models 1-10, NDVI and slope gradient were selected with negative coefficients, and "sandy turbidites" (Dummy 3) was selected with positive coefficients. In models 1-5, "forest" (Dummy 1) was selected with negative coefficient, and "citrus grove" (Dummy 2) was selected with positive coefficients in models 6-10. Additionally, rainfall was selected with a positive coefficient for models 1 and 6. The factors with a positive coefficient indicated an increasing value (or the presence) of factors that increase the landslide hazard. The AIC value for model 1 and 6 (321.1 and 311.7, respectively) were lower than those for models 2-5 and 7-10 (323.9 and 314.2, respectively). Figure 7 shows the ROC of all models; the AUC for models 1 and 6 were larger than that for models 2-5 and 7-10 as well as the F1 scores and the overall accuracy for models 1-5. However, the same accuracy values for models 6-10 were found in the F1 scores and the overall accuracy (Table 3).

Significant Effects of Landslide Conditioning Factors
All the independent variables for the ten models had a VIF of <2 (Appendix A Tables A4 and A5), indicating that all the variables could be simulated together. Table 2 summarized the selected factors; after the stepwise processes of logistic regression analysis, NDVI, "forest," and "citrus grove" were selected as significant factors. Because rainfall was not selected for models 2-5 and 7-10, the same result was obtained for the eight models. For models 1-10, NDVI and slope gradient were selected with negative coefficients, and "sandy turbidites" (Dummy 3) was selected with positive coefficients. In models 1-5, "forest" (Dummy 1) was selected with negative coefficient, and "citrus grove" (Dummy 2) was selected with positive coefficients in models 6-10. Additionally, rainfall was selected with a positive coefficient for models 1 and 6. The factors with a positive coefficient indicated an increasing value (or the presence) of factors that increase the landslide hazard. The AIC value for model 1 and 6 (321.1 and 311.7, respectively) were lower than those for models 2-5 and 7-10 (323.9 and 314.2, respectively). Figure 7 shows the ROC of all models; the AUC for models 1 and 6 were larger than that for models 2-5 and 7-10 as well as the F1 scores and the overall accuracy for models 1-5. However, the same accuracy values for models 6-10 were found in the F1 scores and the overall accuracy (Table 3).

Discussion and Conclusions
The land use category "citrus grove" occupied >50% of the study area ( Figure 4). The LAD for "citrus grove" was the highest among the six land use categories (Figure 4). When LAD was compared with "forest", "citrus grove," and "the others" for each category of conditioning factors, the LAD for "citrus grove" was larger than that for "forest" and "the others" in most cases ( Figure 6). Specifically, the LAD of citrus groves on gentle terrain was higher than that of forests on steeper terrains (32-40° of gradient) ( Figure 6). This result may be attributed to the study location, which is predominantly covered by gentle terrain (78.46%), and "citrus grove" was more distributed on gentle terrain than steep terrain with 42.82%, and 9.92% of the area, respectively. Furthermore, "citrus grove" was selected as a significant factor with positive coefficients in the logistic regression analysis

Discussion and Conclusions
The land use category "citrus grove" occupied >50% of the study area ( Figure 4). The LAD for "citrus grove" was the highest among the six land use categories (Figure 4). When LAD was compared with "forest", "citrus grove," and "the others" for each category of conditioning factors, the LAD for "citrus grove" was larger than that for "forest" and "the others" in most cases ( Figure 6). Specifically, the LAD of citrus groves on gentle terrain was higher than that of forests on steeper terrains (32-40 • of gradient) ( Figure 6). This result may be attributed to the study location, which is predominantly covered by gentle terrain (78.46%), and "citrus grove" was more distributed on gentle terrain than steep terrain with 42.82%, and 9.92% of the area, respectively. Furthermore, "citrus grove" was selected as a significant factor with positive coefficients in the logistic regression analysis for all models.
Thus, citrus groves contributes more to increasing the landslide hazard than forests in the study area.
Some studies [11,14] also reported a higher landslide density of fruit groves to landslides compared with forests, although there are few extensive studies examining the differences in landslide hazard between forest lands and fruit groves. However, there are several possible reasons for fruit groves contributing more to the landslide hazard. First, fruit groves are often developed by cutting the upper side of the slope and leveling the lower side of the slope to build pits for transplanting the seedlings. Generally, manmade slopes are more susceptible to landslides than natural slopes [53]. Manmade slopes and leveled natural slopes tend to change the physical properties of the soil. The main cause of hazard in manmade slopes is excavation [54,55]; geomechanical studies of soft rock slopes state that large deformation of the slope may occur on layered excavations without retaining the structure, due to the unloading of stress, which can lead to slope instability [55,56]. Second, irrigation in fruit grove areas can maintain a high soil water content. Third, fruit groves are composed of a single species. Areas with low biodiversity may have a weak root network system in terms of soil-binding capacity [57]. In addition, soil compaction by trampling may also alter soil physical properties to influence landslide initiation.
This study used radar-based rainfall data with high spatial and temporal resolutions. The result showed 1-h maximum rainfall significantly affected the landslide occurrence ( Table 2). The AIC and AUC of model 1 and 6, including 1-h maximum rainfall, were lower and larger than models 2-5 and models 7-10, respectively, where rainfall indices were excluded. In a previous study [58] that compared rainfall data from XRAIN radar with those based on rain gauges, the reported differences were 2% for normalized bias, 22% for normalized error, 0.96 for correlation coefficient, and 3.5 mm for root mean square error. High-quality radar-based rainfall measurement can improve the evaluations of landslide hazard.
The results of this study suggest that fruit groves should not be developed on slopes near houses or public facilities, where people can be harmed and assets destroyed. This may apply to other areas with similar topographic and geological conditions. The results of this study (i.e., higher contribution in landslide initiation and distribution in fruit groves than forests) are identical to previous studies [11,14]. However, previous studies have only examined landslide initiation among land use types without considering variations in rainfall distributions; this study examined landslide hazard in several steps while considering several rainfall distributions, which expands upon previous work. This is an improvement over the previous studies. However, this study had some limitations. First, this study only included one rainfall event, several events would have been desirable for generalization. Second, historical land use changes and landslides were not considered, which may affect the current landslide hazard [59,60]. Third, the mechanism that leads to higher landslide initiation in fruit groves compared with forests could not be clarified. Examining and comparing the root system as well as hydrological aspects of fruit groves with forests will further increase the understanding of these mechanisms. Finally, other data inputs (i.e., natural vs. manmade slopes, soil characteristics, and root types) for the logistic regression model will help to better understand the landslide hazard of fruit groves. Our final goal is to be able to generalize landslide hazard in fruit groves, and the study presented in this paper is a step in this direction.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Rule of dummy variables on land use.