Vietnam, a tropical country locating within the Southeast Asian typhoon belt, often suffered from tropical cyclones and depressions annually. Recently, abnormal weather phenomena such as torrential rain have been occurring with more frequencies and higher intensities; water-related disasters including floods, debris flows and landslides have also borne catastrophic consequences to the society. Being parallel with the requirements of economic development and population pressure in Northeast region, a large number of arteries connecting residential zones have been constructing and rehabilitating in such hilly region. Exploiting activities from mining, deforestation and infrastructure constructions in mountainous areas have led to the increase of bared land on hill-sides. High rainfall intensities in addition to human interferences in the natural slope condition (slope-cut and tree-cutting on hilly sites), and changes in land cover condition are major triggers to the landfall occurrences in this area. The expansion of the development activities could lead to the rise in slope failures especially during triggering events such as prolonged extreme precipitation, which has been influenced greatly by climate change.
According to   and  , there are two major techniques for landslide susceptibility; they are qualitative method and quantitative method. The qualitative method, including inventory and heuristic approaches, refers to expert-based method; it based on landslide inventory and historical information to evaluate and determine main parameters among them. It also identifies sites with similar conditions of geology and geomorphology to evaluate the susceptibility of failure. The inventory map is implemented using satellite images, ground survey and historical database of landslide occurrence. This method highlights the location and scales of historical landslides, and is seen as one of the most susceptibility mapping techniques. It would give spatial distribution of past landslides and might be utilized for evaluating reducing risks of slope failures on a regional scale; however it does not provide prediction of landslide susceptible areas. The heuristic, which requires long-term landside data and causative parameters, estimates potential landslide from preparatory variables. Such an expert-driven approach depends on the knowledge and experiences of scientists in deciding degree and type of landslide risks. The quantitative methods, including statistical and deterministic methods, utilized mathematical model to assess the probability of slope-failure occurrence so as to identify risk regions. The statistical/stochastic method analyses the relationship between existing landslides and instability factors; it would provide prediction of landslide in region that might not have landslide in the past based on statistical analysis of landslide ground characteristics. However it needs a large amount of data collection and is suitable for medium scale prediction. The deterministic approach requires physical processes to analyse the equilibrium mechanical of potential block, and to calculate a factor of safety  . It can be applied for detailed studies at large scale  . However, difficulties in obtaining, checking and processing data especially large spatial dataset in regional slope stability studies reveal the limitation of this approach. The application of deterministic is often in case of simple landslide types in fairy homogeneous geology condition  .
Recent studies analysed the slope instabilities in North area using the statistical method for a regional scale; the deterministic method is often applied for a local scale. In addition, the Vietnamese Government conducted a national project 2012-2020 to produce the landslide inventory maps for mountainous provinces. In fact, it is difficult to evaluate and obtain records of slope failures locating in remote areas; this project thus implemented investigations along main arteries to map landfall records that are related to the man-made slope failures in populated regions. It is noted that the Northern region is relied strongly on agriculture activities  ; it thus is more prone to changes in weather when it has experienced with the growing number of extreme weather events during the past few years. This paper therefore would like to address the probability of landslide susceptibilities through the development of physical model in a regional scale focusing the roles of rainfall in relation to conditioning factors including changes land-use under the influence of the warming weather.
2. Study Area and Data
The Cau river basin belonging to Bac Kan and Thai Nguyen provinces is the focus of this study. The catchment locates in the Northeast region, lying from 105.48˚ to 106.13˚ East longitude, and from 21.35˚ to 22.32˚ North latitude. According to the historical record of daily rainfall during 1962-2016 which was collected from the Vietnam Meteorological and Hydrological Administration, the average annual precipitation ranges from 1360 mm (Cho Ra) to 2572 mm (Diem Mac); the highest 24-hour rainfall varies from 208 to 475 mm (Figure 1(b)). Regarding to ecological stratification, there are two major sub forest-ecological regions: Low mountains (TV10) and midland (TV12)   . Major forest types are mixed closed evergreen humid forests (TV12), and mixed closed forests on limestone mountains/valleys (TV10, TV12). According to the national statistic data, within the period 2011-2016, this region illustrated significant rises of 1.24 - 1.54 percent per year in the population size. The number of residents in Bac Kan province grew from 300.4 thousand in 2011 to 319 thousand with the density of 65.64 people per square kilometre in 2016. With total people of 1227.4 thousand in 2016, Thai Nguyen province experienced a rapid growth rate of 88 thousand people when compared to population size in 2011; it also brought about the population density of 348.03 people per square kilometre in 2016. Such significant rise in population has burdened high pressure on
Figure 1. (a) The meteorological and rain-gauge stations in the Thiessen network within three provinces Tuyen Quang, Bac Kan & Thai Nguyen; (b) The highest 24-hour precipitation (Rmax) in the Cau basin acorrding to the statistical records from 1962-2016.
requirements of residential area in this rural region. Regarding to the general statistics office of Vietnam, total deforestation area in Bac Kan and Thai Nguyen during the period 2006-2016 were 329 ha and 1546 ha, respectively; whereas only 172.8 ha and 153.2 ha of forest was planted within this period. It is crucial difficult to collect real situation of deforestation in Vietnam; the official statistic data hence, were often smaller than real area of logging in mountainous region. However the statistics reveal that new constructions in relation to inadequate land-use management would lead to significant decreases in vegetation cover in these two provinces.
In this paper, rainfall data from 1960 to 2016 at rain-gauge stations and meteorological stations were adopted from the Institute of Meteorology in Vietnam. Lithology map and Geology map were extracted from seven tiles of the Geological and Mineral Resources Map of Vietnam at the scale of 1:200,000 and in Northeast region. Among 19 formations in this river basin, the north part is mainly covered by PhuNgu and Mia Le formations; whereas Ha Coi formation dominates in South region in Thai Nguyen province (Figure 2(a)). The lithology was classified into seven subclasses: 1) Quaternary deposits; 2) Sedimentary aluminosilicate rocks and sedimentary quartz rocks; 3) Sedimentary carbonate rocks; 4) Mafic-ultramafic magma rocks; 5) Acid-neutral magmatic rocks; 6) Metamorphic rocks with rich aluminosilicate components; and 7) Metamorphic rocks with rich quartz components (Figure 2(b)). The terrain condition was extracted from the Aster Global Digital Elevation Model (GDEM) version 2, which
Figure 2. (a) Formation map in the Cau basin; (b) Lithology map in the Cau basin [Group description: 1) Quaternary deposits; 2) Sedimentary aluminosilicate rocks and sedimentary quartz rocks; 3) Sedimentary carbonate rocks; 4) Mafic-ultramafic magma rocks; 5) Acid-neutral magmatic rocks; 6) Metamorphic rocks with rich aluminosilicate components; and 7) Metamorphic rocks with rich quartz components].
was released in 2011 by the Ministry of Economy, Trade and Industry, Japan (METI) and the National Aeronautics and Space Administration (NASA)  with the spatial resolution of 30 × 30 m (Figure 3(a)). The slope map was derived from DEM map as described in Figure 3(b). The land use map in North Vietnam was adopted from the Japan Aerospace Exploration Agency (JAXA)  . This high resolution land cover map was produced in 2007 and 2015 with nine categories: 1) Water; 2) Urban and built-up; 3) Rice; 4) Other crops; 5) Grasslands; 6) Orchards; 7) Barren; 8) Forest; and 9) Mangrove (Figure 4).
In this paper, the landslide susceptibility is analysed based on the deterministic method using a Transient Rainfall Infiltration and Grid-Based Region Slope-Stability (TRIGRS) model  . In this program, the pore-water pressure and the value of factor of safety (Fs) are calculated on a cell-by-cell basis and can be operated in the geographic information system  . The model evaluates changes in the transient pore pressure and the factor of safety due to rainfall infiltration; it is thus often applied for analysing rain-induced landslides. The infiltration model for initially unsaturated conditions of finite depth of the bedrock was selected for this river basin; as described Figure 5(a) the soil profile was divided into two main layers including a saturated zone beneath the water
Figure 3. (a) The Digital Elevation Model (DEM) map; (b) The slope map.
(a) (b) (c)
Figure 4. Land use map from the Japan Aerospace Exploration Agency (JAXA) regarding to (a) Landuse map in 2007 with the resolution of 100 m; (b) Landuse map in 2015 with the resolution of 100 m; (c) Landuse map in 2015 with the resolution of 50 m.
table, and an unsaturated zone expanding to the ground surface  . The factor of safety Fs at depth Z below the ground surface was derived from the balance between the downslope components of the gravitational driving stress and the resisting stress for saturated soils; Fs was defined by separating the time-variant term and the steady terms as in Equation (1)  :
Figure 5. (a) Conceptual diagram of TRIGRS model; (b) Basic input layers for TRIGRS model.
where φ is the soil friction angle (˚); δ is the slope angle (˚); c is the cohesion (kN/m2); γs is the unit weight of soil (kN/m3); γw is the unit weight of water (kN/m3); Z is the depth below the ground surface (m); t is the time (sec); ψ refers to the ground water pressure head as a function of depth Z and time t (m).
In this study, a scenario-based approach is applied for landslide risk assessment under extreme weather condition; main input layers for TRIGRS model include rainfall intensities and durations in relationship with soil physics, hydro-geology and slope conditions (Figure 5(b)). In addition to the Digital Elevation Model (DEM), information of soil and precipitation was discretized and assigned in each cell in order to simulate rain infiltration as well as vertical flow through the unsaturated zone. Regarding to the thickness of topsoil―a major parameter for the slope instability evaluation, it is difficult to obtain this information and it is not feasible for field investigation upon a large scale region. The spatial distribution of this parameter is usually predicted based on empirical relationship between topography attributes and soil thickness. In this paper, we applied the Saulnier method (1997)  to determine soil thickness in each pixel; this approach is based on the assumption of decreasing linear function between topographic slope and soil depth as described in Equation (2):
In which Zmin and Zmax refer to minimum and maximum values of effective soil depths, respectively; δmin and δmax are the minimal and maximal values of slope angle, respectively; Zi and δi represent the soil depth and slope angle at the calculated cell i.
Moreover, in order to evaluate the role of vegetation cover on the hilly slope, in this article, we considered the contribution of root systems to the shear strength (s) by additional cohesion component (Δc) following the Mohr-Coulomb formulation  (Equation (3)):
where the cohesion c in Equation (1) in vegetation cover area is the sum of soil cohesion cs (kN/m2) and apparent cohesion provided by roots (root cohesion) Δc (kN/m2); δN refers to the normal stress on the shear plane (kN/m2). For simplicity, in this study, we ignored the influence of tree-surcharge; the mechanical impact of tree cover on the shear strength of soil with vegetation cover is therefore represented only by the root cohesion (Δc).
In order to simulate shallow landslide under different extreme precipitation, it is necessary to evaluate spatial and temporal distribution of storm events. As described in Figure 1, there are 8 meteorological stations at 3 provinces in the Northeast region covering the area of the Cau basin including Bac Kan (BK), Cho Ra (CR), and Ngan Son (NS) stations in Bac Kan province, DinhHoa (DH) and Thai Nguyen (TN) stations in Thai Nguyen province, as well as ChiemHoa (CH), Ham Yen (HY), and TuyenQuang (TQ) stations in TuyenQuang province. In this paper, we analysed historical records from 1960-2016 at the 8 aforementioned meteorological stations in addition to 5 rain-gauge stations in Bac Kanprovince (Cho Don (CD) and Phu Thong (PT)) and Thai Nguyen province (Dai Tu (DT), Pho Yen (PY), and Vo Nhai (VN)) covering the entire area of this catchment. Extreme precipitation was selected from storm events that resulted in landslides. Additionally, under the impacts of climate change, the rises in temperature are believed to drive stronger downpours in many regions across the world. The magnitudes of extreme rainstorm under warmer climate were computed using the Probable maximum precipitation (PMP) relation to the moisture maximization method   . The temporal rainfall distributions were also analysed regarding to the most heavy rainstorm events in the study site.
4.1. Scenarios Setup
To evaluate the slope stability condition in the Cau river basin, we applied the scenario-based approach, in which we considered options of rainfall intensity and distribution, characteristics of surface condition, and the soil depth estimation. The actual rainstorm event was selected based on rain records from 2010 to 2016. During this period, the annual rainfall in the study area varied from 1391 mm (2011) to 2029 mm (2013); the highest values of the maximum-24-hour precipitation at most rain-stations occurred in the year 2013. Heavy rain in 2013 also resulted in a great number of landslide occurrences across the region; for example, according to the report of Bac Kan province, the road 258 from Bac Kan to Ba Be with the length of 40 km was closed in 1 month and was rehabilitated in about 3 months as a consequence of landslide occurrences triggering after the downpour in the late of May 2013. In addition, the national project also conducted the landslide inventory map in Bac Kan province in 2013. Despite the fact that the investigation was implemented through the transect walks method across traffic routes in this province, and the landfall often associated with slope-cut activities; this database of landslides still provided important references for the study of sliding issues in the study area. Therefore we utilized heavy rainstorm in 2013 (Figure 6(a)) to analyse the probability of landfall in this catchment under the impact of the historical rain event. The Thiessen polygons were constructed around stations within and near the basin to determine the areas that were represented by the stations as illustrated in Figure 1(a)  ; the Inverse Distance Weight (IDW) interpolation was also applied to estimate the spatial distribution of rainfall in the catchment as described in Figure 6(a)    .
Analysing 68 heavy storm events from 1960 to 2016 in the Northeast provinces, we found that the heavy storm events often lasted in less than 48 hours; the rain durations fluctuated from 5 hours (TuyenQuang station on 1984/6/22) to 53 hours (Bac Kan station on 1990/9/21) as illustrated in Figure 7(a). In this article, proposed extreme rainfall events were hence selected based on the probable maximum precipitation regarding to the duration of 24 hours (Figure 6(b)) and 48 hours  . The analysis of hourly rainfall regard to 68 heavy rainstorm events illustrated that the hourly distributions of precipitation varied significantly and could be grouped into three main distribution types: 1) Type 1: peak rainfall locates at the beginning of the storm; 2) Type 2: peak rain locates at the middle of the storm time; and 3) Type 3: peak precipitation appears at the end of the storm (Figure 7(b)). The values of heavy rainstorm in 2013, as well as 24-hour PMP and 48-hour PMP at 8 stations in the study site are described in Table 1.
Information of soil and hydro-geology in the catchment was gathered from recent studies  -  as well as the field investigations in 2017 and 2018, in which soil samples were tested in the geotechnics laboratory in Thuy Loi university. The average input parameters of soil properties regarding to main geology regions are summarized in Table 2.
In addition to the worst case of bared land across the region (the cohesion c in Equation (1) refers to the soil cohesion cs), the land-use data obtaining from JAXAin two time periods 2007 and 2015 was put in the model to set up scenarios; the cohesion c in Equation (1) is the sum of the soil cohesion cs and the root cohesion Δc. The root cohesion was estimated using the trial-and-error method    ; the average values varied from 0.7 KN/m2 in orchards, savannas and shrub lands to 1.0 KN/m2 in forest cover.
Applying the Saulnier method  (Equation (2)), information of average depth of top soil was gathered from recent studies from the government project of landslide inventory map, and was set up in three options including 1) case 1: upper value of soil depth (0.5 m - 4 m) based on lithology map across the whole region; 2) case 2: average value of soil depth (0.5 m - 3 m) based on lithology
Figure 6. (a) Distribution of selected rainstorm in 2013; (b) Distribution of the 24-hour PMP.
Figure 7. (a) Relationship between total amount and duration of heavy rainstorm regarding to average air temperature during the rain-event; (b) Three main distribution types of heavy rainstorm in Cau basin.
Table 1. Total rain amount regarding to the historical rainfall event as well as the 24-hour PMP and the 48-hour PMP (mm).
Table 2. Average values of soil properties in Cau basin.
Because the information of groundwater table was unavailable, it was assumed that the initial groundwater table was at the bottom of the weather soil layer in summer time without antecedent rain. Based on options of rainstorm event amount (R), rainfall distribution (D), land cover condition (L), and soil depth assumption (S), we set up scenarios (RDLS) to compute the factor of safety for the whole region of Cau river basin as described in Table 4.
4.2. Evaluation of Landslide Susceptibility
In case of actual rainstorm distribution in 2013, we calculated the factor of safety Fs regarding to three options of soil-depth as well as the changes in land-use option i.e. from land use data in 2007 (JAXA) (Ha12) to land-use data in 2015 (JAXA) (Ha21, Ha22, Ha31, Ha32 & Ha33) in addition to the worst option of land cover with bare soil (Ha02). Illustrations of factor of safety (Fs) relation to the historical rainfall event are described in Figure 9; Figure 10 shows the determination of Fs regarding to changes in the extreme-precipitation event including 24-hour PMP (scenarios 1132 and 1332), and 48-hour PMP (scenarios 2332 and 2333). With regard to historical rainfall event, case 1 in soil-depth determination (Ha21 and Ha31) provided lower Fs values than the results from case 2 and case 3. In case of bared soil across the whole region (Ha02), more unstable areas were seen in northern and south-eastern parts of the catchment when compared to results from the three options Ha12, Ha22, and Ha32. With regard to historical landslide data in 2013 from the Government project, we obtained shapes of 44 landslide points using the application of Google Earth. The summaries of the Fs values within those polygons were illustrated in Figure 11; about 60 percent of the sliding area had the factor of safety that was less than 2.0. The case 1 of soil-depth (S = 1) still obtained the highest proportion of area having low values of Fs; while the option using the land-cover data from 2007 (Ha12) indicated only 51% of the area that had values of Fs below 2.0. In addition, throughout two field trips in September 2017 and September 2018, we investigated 56 sliding points in the catchment. Since the shape and area of each landfall was not mapped, we tried to evaluate the responses of slope stability regarding to scenario-based simulations within a region of 30 meter × 30 meter at each slide point. About 80 percent of the total area obtained Fs from 0.82 to 2.
Table 3. Summary of options in soil depth range in the study area.
(a) (b) (c)
Figure 8. Soil depth determination regarding to scenarios (a) Case 1; (b) Case 2; (c) Case 3.
The case 1 of soil-depth still proved as the upper boundary in the relationship between accumulation area and factor of safety. When the rainfall amount increased to extreme value of 24-hour PMP (1332) and 48-hour PMP (2332), area of cells that could be prone to the instability was higher than that of option Ha32 (Figure 12).
In addition, in order to estimate the changes in Fs, we considered the option Ha32 as a base case and computed the difference between Fs_Ha32 and factor of safety for other scenarios throughout the basin (Figure 13). In general, the values of average ΔFs in the whole basin were 0.218 (Ha32_Ha31), 0.196 (Ha32_2333), 0.189 (Ha32_2332), 0.18 (Ha32_Ha21), 0.097 (Ha32_1332), 0.069 (Ha32_Ha02), 0.040 (Ha32_1132), 0.014 (Ha32_Ha33), −0.01 (Ha32_Ha12), and
Table 4. (a) Scenarios [RDLS] summary; (b) Description of scenarios [RDLS]*.
*.For example, scenario Ha32 refers to the option applying the historical rainfall event (R = H) having an actual distribution (D = a) using JAXA land-use in 2015 with the resolution of 50 meter (L = 3), and the soil depth case 2 (S = 2).
−0.046 (Ha32_Ha22). Within the range of Fs_Ha32 from 0 to 1.5, the option Ha12 and Ha22 obtained smaller threshold of ΔFs (negative values); using the same resolution of 100 meter in the land-use data, changes in the land cover such as from forest cover (Ha12 in 2007) to crop land and construction sites (Ha22 in 2015) revealed the reduction of the slope stability regarding to the range of Fs below 1.0. Interestingly when the whole region is assumed to be replaced by soil (Ha02), a significant reduction of Fs was shown when compared to the option of changes in land-use database (Ha12 and Ha22) and changes in soil-depth option (Ha33). The worst cases in assessing the slope stability still came from the case 1 of soil-depth option (Ha21 and Ha31). The differences between two aforementioned lines within the range of Fs_Ha32 from 1 - 1.5 illustrated the great impact of grid-cell resolution to the model simulation. When the rainstorm amount changed from actual event (Ha32) to 24-hour-PMP (1132 and 1332) and 48-hour-PMP (2332 and 2333), the factor of safety Fs also declined in the whole basin.
According to  , the stability condition was classified into four groups: 1) Unstable when Fs < 1, the stability requires stabilizing factors; 2) Quasi-stable when 1 ≤ Fs < 1.25, the instability would be resulted from minor destabilizing factors; 3) Moderately stable when 1.25 ≤ Fs < 1.5, the instability would be resulted from
(a) (b) (c) (d) (e) (f)
Figure 9. The factor of safety (Fs) in case of rainfall event in 2013 regarding to scenarios: (a) Ha02; (b) Ha12; (c) Ha22; (d) Ha31; (e) Ha32; (f) Ha33
moderate destabilizing factors; and 4) Stable when Fs ≥1.5, the instability would be resulted from major destabilizing factors. Above tested cases also indicated that about 30% - 50% of the historical landslide area obtained the Fs below 1.5. In fact, most of historical landslides were investigated along main arteries and often associated with road-cut slopes; the instability condition of such man-made slopes were more prone to changes in the environmental condition
(a) (b) (c) (d)
Figure 10. The factor of safety (Fs) in case of the 24-hour PMP and the 48-hour PMP regarding to scenarios: (a) 1132; (b) 1332; (c) 2332; (d) 2333.
than that in case of natural slope (as in the model simulation using the global DEM data). It thus reveals that the TRIGRS is sensitive to the DEM data input especially in case of changes in slope conditions. However, the simulation demonstrated the agreement with the risk of slope-failures in a relationship with destabilizing factors including slope-cutting activities.
Differences in the resolution of input data also influenced the output data as illustrated in scenarios Ha22 (JAXA 2015 using 100-meter resolution) and Ha32
Figure 12. (a) Distribution of the factor of safety (Fs) in accordance with percent of accumulation area of the regions covering the landslide points in 2017 & 2018; (b) Detail of Figure 12(a) regarding to the Fs below 2.0.
(JAXA 2015 using 50-meter resolution). As discussed in Section 4.2, changes in land-cover between two time periods resulted in differences in the value of factor of safety in the basin; within the cells under the unstable condition (Fs < 1), the factor of safety regarding to the land-cover in 2015 decreased by 0.01 - 0.02 when compared to the values of Fs at the same cells in case of the land-use 2007.
In addition, under the changes in climatic condition, extreme precipitation would result in higher probability of unstable regions. It is noted that under the same amount of the storm-rain and storm duration, for example 24-hour-PMP, the landslide prone areas would expanse in case of changing the storm distribution from type 1 (scenario 1132) to type 3 (scenario 1332). Comparing the rise of ΔFs in case of changing rain-event from actual storm in 2013 to 24-hour-PMP (126.9 mm - 332.4 mm rise in rain amount), the factor of safety still decreased significantly when the downpour total of 24-hour-PMP increased by 21.7 mm - 169.1 mm in a further 24 hours to reach the magnitude of 48-hour-PMP (Figure 13). This reveals the importance of three major components of a storm events
Figure 13. (a) Relationship between factor of safety for option Ha32 (Fs_Ha32) and Delta Fs (ΔFs) [difference between Fs_Ha32 and factor of safety for other options (Fs_RDLS)]; (b) Detail of Figure 13(a) regarding to the Fs_Ha32 below 1.5.
including rain-amount, rain-duration, and rain-distribution to the stability condition of hilly sites.
In conclusion, the model simulations obtained uncertainties raising from the estimations of soil-depth and groundwater based on geologic and hydrologic background for a regional scale, as well as the application of average values of soil properties and the same thresholds of root cohesion Δc. Subjective errors in drawing polygons from the Google Earth soft-ware, especially in case of small landslides also hindered the evaluation of the historical landfalls. However, the findings would provide important information for authorities in developing adequate land-management in the river basin considering the relationship between landslide-prone area and changes in land-cover situations.
This research was supported by the Vietnamese Government scholarship (project 911); and the Social Implementation Program on Climate Change Adaptation Technology (SI-CAT) of MEXT, Japan.
 Chalkias, C., Ferentinou, M. and Polykretis, C. (2014) GIS-Based Landslide Susceptibility Mapping on the Peloponnese Peninsula, Greece. Geosciences, 4, 176-190.
 Dieu, T.B., Owe, L., Inge, R. and Dick, O. (2011) Landslide Susceptibility Analysis in the Hoa Binh Province of Vietnam Using Statistical Index and Logistic Regression. Natural Hazards, 59, 1413-1444.
 Metternicht, G., Hurni, L. and Gogu, R. (2005) Remote Sensing of Landslides: An Analysis of the Potential Contribution to Geo-Spatial Systems for Hazard Assessment in Mountainous Environments. Remote Sensing of Environment, 98, 284-303.
 Zhou, G., Esaki, T., Mitani, Y., Xie, M. and Mori, J. (2003) Spatial Probabilistic Modeling of Slope Failure Using an Integrated GIS Monte Carlo Simulation Approach. Engineering Geology, 68, 373-386.
 van Westen, C.J. , Rengers, N., Terlien, M.T.J. and Soeters, R. (1997) Prediction of the Occurrence of Slope Instability Phenomenal through GIS-Based Hazard Zonation. Geology of Rundschau, 86, 404-414.
 D.L.M. (2015) Adjustment of the Land Use Planning up to Year 2020 and the Schedule of Land Use in the Later Term in a National Scale (2016-2020). Department of Land Management in Vietnam. (In Vietnamese)
 Le, T.T.T. and Kawagoe, S. (2017) Landslide Detection Analysis in North Vietnam Base on Satellite Image and Digital Geographical Information—Landsat 8 Satellite and Historical Data Approaches. Japan Society of Civil Engineering, 73, I_239-I_249.
 Nguyen, N.L., Do, D.S., Nguyen, X.Q., Tran, V.L., Ngo, D.Q., Tran, V.C. and Nguyen, D.K. (2011) Forest Ecological Stratification in Vietnam. Research Center on Forest Ecology and Environment. (In Vietnamese)
 ASTER GDEM Validation Team (2011) ASTER Global Digital Elevation Model Version 2—Summary of Validation Results. NASA and METI, 27.
 J.A.X.A. (2016) High-Resolution Land Use and Land Cover Map of Northern Region of Vietnam (Released in Sep.2016/Version 16.09).
 Baum, R.L., Savage, W.Z. and Godt, J.W. (2008) TRIGRS—A Fortran Program for Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis, Version 2.0. U.S. Geological Survey Open-File Report, 75.
 Baum, R.L., Godt, J.W. and Savage, W.Z. (2010) Erratum: Estimating the Timing and Location of Shallow Rainfall-Induced Landslides Using a Model for Transient, Unsaturated Infiltration. Journal of Geophysical Research, 115, 1-26.
 Pollen-Bankhead, N. and Simon, A. (2009) Enhanced Application of Root-Reinforcement Algorithms for Bank-Stability Modeling. Earth Surface Processes and Landforms, 34, 471-480.
 Le, T.T.T., Kawagoe, S. and Sarukkalige, R. (2018) Estimation of Probable Maximum Precipitation at Three Provinces in Northeast Vietnam Using Historical Data and Future Climate Change Scenarios. Journal of Hydrology: Regional Studies, unpublished.
 Duc, D.M. (2004) Numerical Simulation for Slope Stability Analysis and Shallow Landslide Prediction in Bac Kan Province, Northeast Vietnam. International Sympposium on Geoinformatics for Spatial Infrastructure Development in Earth and Allied Sciences, Hanoi, Vietnam.
 Duc, D.M., Duc, D.M. and Ngoc, D.M. (2018) Effects of Residual Soil Characteristics on Rainfall-Induced Shallow Landslides Along Transport Arteries in Bac Kan Province, Vietnam. In: Advances and Applications in Geospatial Technology and Earth Resources, Springer, Cham.
 Duc, D.M., Khang, D.Q., Binh, N.V., Loi, V.V., Sinh, P.T., Duc, P.H. and Truc, N.N. (2012) Mechanism of Sinkhole Formation in Underground Limestone Karstic Areas (A Case Study of Na Tum Area, Bac Kan Province). Journal of Mining and Earth Science, 38, 16-24. (In Vietnamese)
 Kuriakose, S.L. (2010) Physically-Based Dynamic Modelling of the Effect of Land Use Changes on Shallow Landslide Initiation in the Western Ghats of Kerala, India. Ph.D. Thesis, University of Twente, The Netherlands.
 Docker, B.B. and Hubble, T.C.T. (2008) Quantifying Root-Reinforcement of River Bank Soils by Four Australian Tree Species. Geomorphology, 100, 401-418.
 Tran, T.V., Lee, G., An, H. and Kim, M. (2017) Comparing the Performance of TRIGRS and TiVaSS in Spatial and Temporal Prediction of Rainfall-Induced Shallow Landslides. Environmental Earth Sciences, 76, 315.