2-D Cross-Plot Model Analysis Using Integrated Geophysical Methods for Landslides Assessment

Featured Application: The two-dimensional (2-D) cross-plot model represents an alternative integrated method for different geophysical parameters. The outcome highlights the beneﬁts, con-sidering that the geophysical parameter in landslide studies as the cross-plot analysis was introduced to enhance the subsurface resolution where the volume of mass sliding was successfully calculated. Abstract: The large or small scale of a landslide is a natural, widespread process, resulting from the downward and outward movement of slope-forming materials, such as sculpting the landscape. Characterized landslide material and properties’ inhomogeneities conditions become a challenge as the process required the availability of a wide range of data, observations, and measurements with an evaluation of geological and hydrological conditions. Detailed investigations represent an essential component of the landslide risk mitigation process, relying on subsurface investigations, discrete subsurface sampling, and laboratory tests. To extend this approach, seismic refraction and two-dimensional (2-D) resistivity were utilized to study the landslides activities in Ulu Yam. The cross-plot analysis was introduced to integrate the geophysical results based on the criteria of the model. Velocity distributions from seismic refraction revealed the stiffness of the soil, where weak zones identiﬁed with values of V p ≤ 1200 m/s, deﬁned as threshold frequency for failure to occur. The 2-D resistivity shows that the weak zones were identiﬁed with resistivity values of <1200 Ω m. The 2-D cross-plot model gives a comprehensive interpretation where a low velocity and resistivity value represents the failure plane of materials to failure. The volume of mass sliding was calculated based on retrieved information from the model.


Introduction
Landslides are characterized by a mass movement of soil and rocks on slopes that occur across the globe due to several factors [1,2]. Landslides, known as the most destructive types of geohazards, affect human activities, and have serious socioeconomics consequences [3]. Landslides are referred to as 'mass-wasting', which occurs due to the down-slope movement of soil and rocks in a large variety of shapes and volume [3,4]. The various types of landslides are differentiated by the kinds of material involved and the type of movement. Varnes (1978) classified the landslides into six different categories, namely falls, topples, slides (rotational and translational), lateral spreads, flows, and composites (a combination of types of movements) [5]. Landslides are widespread in many regions, especially under tropical or subtropical climates with prolonged and intense rainfall [6,7]. The worldwide distribution of landslides is not uniform, with landslides occurring primarily where the requisite topographic, climatic, and environmental conditions are prevalent [8].
Landslides are often considered a secondary hazard, forming part of the problem of cascading hazards [9], as they are triggered by catastrophic events, such as storms, floods, volcanic eruptions, and earthquakes [10]. Landslides have been found to have size-frequency distributions that follow a power function for both area and volume of large events [11][12][13][14][15]. The above-referenced landslide studies have been conducted in different locations and are associated with different triggering mechanisms [16]. This phenomenon may occur in natural slopes in a variety of materials, including residual and colluvial soils [17]. Although the impact of this geohazard is inevitable, it can be significantly reduced by increasing the capacity to assess and predict the risks using different mitigation methods [3]. In past decades, numerous modeling methods have been designed to assess slope stability, predict slope response to various triggers and evaluate slope deformation [3]. Investigation of subsurface landslide features are necessary to provide preliminary input for forwarding modeling and subsequent predictions of potential failure events [18,19]. However, such models still require access to detailed information of geological, mechanical, and hydrogeological properties, as well as conditions surrounding boundaries. Therefore, integrated multidisciplinary approaches are essential to being implemented in landslide investigations to improve the reliability of the deterministic model. The increase in geophysical methods has become more widespread. The comprehensive result with proper numerical modeling helps to develop the predictive model, especially for complex landslides event. Due to complexity, classic geotechnical problem for stability assessment give verifiable information on the mechanical and hydraulic characteristics of the landslides, but in a discrete point of the subsurface [20]. Satellite and airborne methods provide useful information regarding the surface characteristics of the slope, such as geomorphological features, areas of extension of the landslide body, superficial displacement, and velocity [21,22]. However, all of this information has not included subsurface features, which plays a vital role in landslide characteristics. In-situ geophysical methods, such as resistivity and seismic refraction, are capable of evaluating physical parameters, directly or indirectly, such as lithology of the subsurface, hydrology, and geotechnical characteristics [23]. The 2-D resistivity is widely used in evaluating slope deformations and in predicting the surface movements [24][25][26][27]. The information obtained is useful to define the geological setting of the investigated subsoil, reconstruct the landslide body, and volume of the slide material [23,28]. Seismic refraction has been proven in detecting the landslide extent and material velocity of the subsurface [29,30]. These methods are less invasive and provide comprehensive information on a large volume of soil, thus overcoming the point-scale features of geotechnical measurements. Successful geophysical methods depend on the presence of a significant contrast in the physical properties of the different lithological units [31]. Therefore, 2-D resistivity and seismic refraction surveys were utilized to visualize the failure zones as a purpose to identify the subsurface characteristics with the trigged factors for the slope failure for this study. As the continuous work from both methods, a cross-plot analysis with the development of a 2-D cross plot model was introduced based on the integration of resistivity and velocity values in order to enhance the subsurface resolution besides the provided comprehensive subsurface images.

Malaysia Rainfall-Landslides Overview
Numerous landslide tragedies in Malaysia involve a large number of fatalities and economic loses. The slope failure events are mostly triggered by heavy rainfall due to the tropical rainforest climate with relatively high temperatures and receive rainfall throughout the year. In a tropical region, such as Malaysia, rainfall-induced landslides have become a common geohazard phenomenon, where the soil typically deposits by residual soils [32,33]. Figure 1a shows a map of the landslide-prone area in Malaysia [34]. The significant thick-ness of the unsaturated zone with a deep groundwater table shows the features of residual tropical soil [35]. Malaysia experiences a hot and humid climate all year around with temperature ranges from 22-32 degrees Celsius, with annual rainfall from 2000-2500 mm, and extremities during southwest (occurring during April-October) and northeast (occurring during October-February) monsoons [36]. The formation of residual soil is caused by intense chemical weathering due to the combination of climatic parameters. Under certain geological and topography effects, the rainfall becomes a major factor that triggers the occurrence of landslides, involving both natural and cut slopes. The intense precipitation, such as in highland area, causes the situation to become more hazardous. The hilly area is stable (of a certain range of water saturation); however, under the critical limit, the shear resistance is enhanced due to the effect of suction, which creates an apparent cohesion between particles [37]. Major landslides occur due to water saturation, exceeding a critical limit in certain parts of the slope [24]. The infiltration of rainfall into the soil and its effect on slope failure has been of particular interest for the past few years [38][39][40][41]. The rainfall intensity, soil properties, groundwater table, and slope geometry plays an important role in rain-induced landslide events [42]. The continuous infiltration of rainwater will increase the pore-air pressure and decrease the matrics suction, which leads to instability [43,44]. According to Collins and Znidarcic, two distinct failure mechanisms have been observed for rainfall-induced landslides [45]. For the first mechanism, a significant build-up of positive pressures is observed in a low area on the slope, or along the soil/bedrock interface. Movements along the sliding surface lead to liquefaction along with it, resulting in rapid movements, long run-out distances and, finally, complete liquefaction of the failed mass [46]. In the second mechanism, the soil in the unsaturated state, and slope failure, is due mainly to rainfall infiltration, and a loss in shear strength when soil suctions are decreased or dissipated [47].
Infiltration and near-surface flow create an increase in pore water pressure that causes the stress path to move nearly horizontally to intersect a failure envelope, initiating a failure. In unsaturated loose soils, suction decrease and coupled volumetric collapse may be involved in the failure process [48]. Landslides in weathered granitic rocks for tropical and humid regions have been well documented as the phenomena are strongly affected by types of weathering. In a tropical region, such as Malaysia, an in-depth weathering profile undergoes a thickness of up to 100 m [49]. The nature of weathered material and mass structure is highly controlled the failure mechanism. The soil or rock at the ground surface is highly weathered and tends to slide with a loose condition. The high erodibility with the existence of relict discontinuities, sharp soil-rock boundary, dynamic weathered material under the saturated condition, and the inhomogeneity of weathering structures will lead to frequent slope failure events. The presence of a fracture and structural discontinuity causes cohesion and internal friction decreases or is reduced by severing bonds between particles, creating the discrete failure surfaces [50,51]. Therefore, evolving scientific comprehension of the geological and physical process of earth's rheology is vital to assess and develop a predictive model for hazard problems. Comprehensive knowledge of the slope, such as identification and features of unstable areas, the volume of the affected slope, and others, are the pre-requisite analyses that are important in hazard assessments.

Study Area
Generally, the Peninsular Malaysia granite is distributed into three parallel belts, known as Western (WB), Central (CB), and Eastern belts (EB). The Bentong-Raub zone separates the Western from the Central Belt (Figure 1b). Western granite is characterized by a massive mountain range extending from Malacca in the south to Thailand in the north [52]. Two main batholith masses can be distinguished in WB known as Main Range batholith on the east part and Bintang batholith in the west [53]. The study area was located on the East of Kuala Lumpur, underlain by Kuala Lumpur granite of the Main Range batholith, which intruded into folded and regionally metamorphose clastic and calcareous Paleozoic rocks. The Kuala Lumpur granite (Figure 1c) consists of coarse to medium grain biotite granite and is commonly intruded by microgranite, coarse grained leucogranite, aplite, and pegmatites. The granite texture generally varies from coarse to medium grain, porphyritic to slightly porphyritic, with colors ranging from white to pale grey. The major minerals are K-feldspar, plagioclase, and quartz, while biotite, muscovite, and tourmaline usually occur in minor amounts, except in the late phase differentiates, where these minerals may be dominant [54]. The overburden materials consist mainly of completely weathered residual soils, which are derived from the weathering of granitic rock and generally consist of silty sand and sandy silt with traces of gravel [55]. The Ulu Yam landslide was considered as a rotational landslide type, where the crown and the toe of the landslide was identified in the study area. Figure 1e shows the schematic diagram of the landslide.
Appl. Sci. 2020, 10, x FOR PROOFREADING 4 of 24 grain biotite granite and is commonly intruded by microgranite, coarse grained leucogranite, aplite, and pegmatites. The granite texture generally varies from coarse to medium grain, porphyritic to slightly porphyritic, with colors ranging from white to pale grey. The major minerals are K-feldspar, plagioclase, and quartz, while biotite, muscovite, and tourmaline usually occur in minor amounts, except in the late phase differentiates, where these minerals may be dominant [54]. The overburden materials consist mainly of completely weathered residual soils, which are derived from the weathering of granitic rock and generally consist of silty sand and sandy silt with traces of gravel [55]. The Ulu Yam landslide was considered as a rotational landslide type, where the crown and the toe of the landslide was identified in the study area. Figure 1e shows the schematic diagram of the landslide.

Field Procedure
In this study, 2-D resistivity and seismic refraction are designed in line. Six 2-D resistivity survey lines (HR1-HR6) were designed with a length of 200 m (m) for HR1 and HR4 and 300 m for HR2, HR3-HR6. HR1-HR3 was oriented parallel to the slope failure, while HR4-HR6 was in a perpendicular direction to the slope (Figure 1d). The Aktiebolaget Elektrisk Malmletning-Signal Averaging System (ABEM SAS 4000) was utilized to acquire the data with an electrode spacing of 5 m. The pole-dipole array was used for data protocol acquisition in providing comprehensive data resolutions. In seismic refraction, three survey lines (HS1-HS3) were designed with a length of 230 m and 115 m for HS1 and HS2-HS3, respectively, with geophone spacing of 5 m. HS1 and HS2 were designed parallel to the slope, while HS3 was oriented perpendicular to the slope. ABEM Terraloc Mark 8 (Mk8) was used as a seismograph to record the seismic waves. Table 1 shows the detailed of the survey lines. In this survey, four boreholes (BH) were provided, with two boreholes located inside the landslide area and another two outside of the landslide boundary ( Figure 1d). The minimum electrodes spacing of 2-D resistivity and the length of the survey were chosen based on the research objectives and accessibility of the area. There are several factors considered in deciding the length of the survey, included the site constraint and depth of investigation. The poledipole array was chosen based on the deepest subsurface penetration, as well as the ability to provide a good resolution. Otherwise, the accessibility and site constraints were counted by selecting the best array to use ( Figure 2). In addition, several factors that needed to be taken into account included the sensitivity of array to vertical and horizontal changes in subsurface resistivity, depth of investigation,  [56,57]); (c) geology map (modified from [58]); (d) study area (the survey lines of two-dimensional (2-D) resistivity and seismic refraction was shown with red and yellow line respectively); (e) schematic diagram of the rotational-landslide of Ulu Yam.

Field Procedure
In this study, 2-D resistivity and seismic refraction are designed in line. Six 2-D resistivity survey lines (HR1-HR6) were designed with a length of 200 m (m) for HR1 and HR4 and 300 m for HR2, HR3-HR6. HR1-HR3 was oriented parallel to the slope failure, while HR4-HR6 was in a perpendicular direction to the slope (Figure 1d). The Aktiebolaget Elektrisk Malmletning-Signal Averaging System (ABEM SAS 4000) was utilized to acquire the data with an electrode spacing of 5 m. The pole-dipole array was used for data protocol acquisition in providing comprehensive data resolutions. In seismic refraction, three survey lines (HS1-HS3) were designed with a length of 230 m and 115 m for HS1 and HS2-HS3, respectively, with geophone spacing of 5 m. HS1 and HS2 were designed parallel to the slope, while HS3 was oriented perpendicular to the slope. ABEM Terraloc Mark 8 (Mk8) was used as a seismograph to record the seismic waves. Table 1 shows the detailed of the survey lines. In this survey, four boreholes (BH) were provided, with two boreholes located inside the landslide area and another two outside of the landslide boundary ( Figure 1d). The minimum electrodes spacing of 2-D resistivity and the length of the survey were chosen based on the research objectives and accessibility of the area. There are several factors considered in deciding the length of the survey, included the site constraint and depth of investigation. The pole-dipole array was chosen based on the deepest subsurface penetration, as well as the ability to provide a good resolution. Otherwise, the accessibility and site constraints were counted by selecting the best array to use ( Figure 2). In addition, several factors that needed to be taken into account included the sensitivity of array to vertical and horizontal changes in subsurface resistivity, depth of investigation, horizontal data coverage, and the signal strength [59]. The pole-dipole array shows good horizontal coverage with significantly higher signal strength. The sensitivity section shows higher for large n factor as this array is more sensitive to vertical structure with good horizontal coverage [59]. Res2Dinv software was used to process and interpret the data by using a mathematical algorithm to generate a 2-D resistivity inversion model. The inversion routine used by the program is based on the smoothness-constrained least-square method [60-63].
Appl. Sci. 2020, 10, x FOR PROOFREADING 6 of 24 horizontal data coverage, and the signal strength [59]. The pole-dipole array shows good horizontal coverage with significantly higher signal strength. The sensitivity section shows higher for large n factor as this array is more sensitive to vertical structure with good horizontal coverage [59]. Res2Dinv software was used to process and interpret the data by using a mathematical algorithm to generate a 2-D resistivity inversion model. The inversion routine used by the program is based on the smoothness-constrained least-square method [60-63]. The seismic refraction method was employed in this survey with a 5 kg sledgehammer and was used as a source of impulsive waves with 28 Hz geophones as a detector. A remote trigger geophone was used to trigger the seismograph at instant impact as the sledgehammer smacked on a steel striker plate. Several stacking is needed to increase the signal-to-noise ratio by repeatedly smacking the striker plate to produce the energy. In this survey, seven shot points were carried out, as five shot points were inline, while the others were the positive and negative offsets. The offset was important in delineating the last layer of the media as it affected the depth of penetration. Figure 3 shows the schematic diagram of seismic refraction data acquisition. FIRSTPIX v4.21 software is used to pick the first arrival time, and the 2-D model was generated by calculating the ray path tracing of the arrival time of every shot point using SeisOpt@2D. Generally, SeisOpt@2D uses a nonlinear optimization technique called adaptive simulated annealing involves forward modeling. Test velocity models were generated through which travel times were calculated. The calculated travel time was compared with observed data, and the error that exits between the data will be optimized and generate a velocity model with minimum travel time error.

2-D Cross-Plot Analysis Model Development
Cross-plot analysis widely implements in the oil and gas industry for formation evaluations and lithological delineations [64][65][66][67]. This method applies multiple input parameters, such as electrical resistivity, nuclear, and acoustic logging in determining lithological reservoir characteristics. Inazaki, 2011, used cross-plot analysis to integrate the S-wave velocity and resistivity result as qualitative vulnerability assessment in producing a map of permeability distribution and N-value [68]. Hayashi and Suzuki, with Immamura et al., have shown a similar application of cross-plot analysis where degrees of compaction (loose or dense) are divided using correlations between S-velocity values and N-values, and soil types by comparing to borehole data [69,70]. Cross-plot analysis between seismic refraction and 2-D resistivity is an analytical approach where different soil type subsurface conditions are classified based on their seismic velocity and resistivity values [71]. Instead of individually analyzing the result from each method, four-quadrant criteria, based on the ranges of seismic velocity and electrical resistivity, are introduced. Based on the measured seismic velocity and resistivity at a given location, and where that point falls within four quadrants (Q1-Q4), the soil types and vulnerability of the subsurface is estimated. A new model of cross-plot analysis was introduced in integrating two different geophysical parameters as a goal to enhance the subsurface resolution using The seismic refraction method was employed in this survey with a 5 kg sledgehammer and was used as a source of impulsive waves with 28 Hz geophones as a detector. A remote trigger geophone was used to trigger the seismograph at instant impact as the sledgehammer smacked on a steel striker plate. Several stacking is needed to increase the signal-to-noise ratio by repeatedly smacking the striker plate to produce the energy. In this survey, seven shot points were carried out, as five shot points were inline, while the others were the positive and negative offsets. The offset was important in delineating the last layer of the media as it affected the depth of penetration. Figure 3 shows the schematic diagram of seismic refraction data acquisition. FIRSTPIX v4.21 software is used to pick the first arrival time, and the 2-D model was generated by calculating the ray path tracing of the arrival time of every shot point using SeisOpt@2D. Generally, SeisOpt@2D uses a nonlinear optimization technique called adaptive simulated annealing involves forward modeling. Test velocity models were generated through which travel times were calculated. The calculated travel time was compared with observed data, and the error that exits between the data will be optimized and generate a velocity model with minimum travel time error.
Appl. Sci. 2020, 10, x FOR PROOFREADING 6 of 24 horizontal data coverage, and the signal strength [59]. The pole-dipole array shows good horizontal coverage with significantly higher signal strength. The sensitivity section shows higher for large n factor as this array is more sensitive to vertical structure with good horizontal coverage [59]. Res2Dinv software was used to process and interpret the data by using a mathematical algorithm to generate a 2-D resistivity inversion model. The inversion routine used by the program is based on the smoothness-constrained least-square method [60-63]. The seismic refraction method was employed in this survey with a 5 kg sledgehammer and was used as a source of impulsive waves with 28 Hz geophones as a detector. A remote trigger geophone was used to trigger the seismograph at instant impact as the sledgehammer smacked on a steel striker plate. Several stacking is needed to increase the signal-to-noise ratio by repeatedly smacking the striker plate to produce the energy. In this survey, seven shot points were carried out, as five shot points were inline, while the others were the positive and negative offsets. The offset was important in delineating the last layer of the media as it affected the depth of penetration. Figure 3 shows the schematic diagram of seismic refraction data acquisition. FIRSTPIX v4.21 software is used to pick the first arrival time, and the 2-D model was generated by calculating the ray path tracing of the arrival time of every shot point using SeisOpt@2D. Generally, SeisOpt@2D uses a nonlinear optimization technique called adaptive simulated annealing involves forward modeling. Test velocity models were generated through which travel times were calculated. The calculated travel time was compared with observed data, and the error that exits between the data will be optimized and generate a velocity model with minimum travel time error.

2-D Cross-Plot Analysis Model Development
Cross-plot analysis widely implements in the oil and gas industry for formation evaluations and lithological delineations [64][65][66][67]. This method applies multiple input parameters, such as electrical resistivity, nuclear, and acoustic logging in determining lithological reservoir characteristics. Inazaki, 2011, used cross-plot analysis to integrate the S-wave velocity and resistivity result as qualitative vulnerability assessment in producing a map of permeability distribution and N-value [68]. Hayashi and Suzuki, with Immamura et al., have shown a similar application of cross-plot analysis where degrees of compaction (loose or dense) are divided using correlations between S-velocity values and N-values, and soil types by comparing to borehole data [69,70]. Cross-plot analysis between seismic refraction and 2-D resistivity is an analytical approach where different soil type subsurface conditions are classified based on their seismic velocity and resistivity values [71]. Instead of individually analyzing the result from each method, four-quadrant criteria, based on the ranges of seismic velocity and electrical resistivity, are introduced. Based on the measured seismic velocity and resistivity at a given location, and where that point falls within four quadrants (Q1-Q4), the soil types and vulnerability of the subsurface is estimated. A new model of cross-plot analysis was introduced in integrating two different geophysical parameters as a goal to enhance the subsurface resolution using

2-D Cross-Plot Analysis Model Development
Cross-plot analysis widely implements in the oil and gas industry for formation evaluations and lithological delineations [64][65][66][67]. This method applies multiple input parameters, such as electrical resistivity, nuclear, and acoustic logging in determining lithological reservoir characteristics. Inazaki, 2011, used cross-plot analysis to integrate the S-wave velocity and resistivity result as qualitative vulnerability assessment in producing a map of permeability distribution and N-value [68]. Hayashi and Suzuki, with Immamura et al., have shown a similar application of cross-plot analysis where degrees of compaction (loose or dense) are divided using correlations between S-velocity values and N-values, and soil types by comparing to borehole data [69,70]. Cross-plot analysis between seismic refraction and 2-D resistivity is an analytical approach where different soil type subsurface conditions are classified based on their seismic velocity and resistivity values [71]. Instead of individually analyzing the result from each method, four-quadrant criteria, based on the ranges of seismic velocity and electrical resistivity, are introduced. Based on the measured seismic velocity and resistivity at a given location, and where that point falls within four quadrants (Q1-Q4), the soil types and vulnerability of the subsurface is estimated. A new model of cross-plot analysis was introduced in integrating two different geophysical parameters as a goal to enhance the subsurface resolution using P-wave velocity and resistivity data. Based on the previous Hayashi model, four quadrants used to divide the parameters. Figure 4a shows the schematic relationship between P-wave velocity and resistivity data.
Appl. Sci. 2020, 10, x FOR PROOFREADING 7 of 24 P-wave velocity and resistivity data. Based on the previous Hayashi model, four quadrants used to divide the parameters. Figure 4a shows the schematic relationship between P-wave velocity and resistivity data. In this research, the cross-plot analysis was implemented to integrate both resistivity and seismic refraction results, which were in line with each other. The goal of this analysis is to provide detailed subsurface images so as to enhance the visualization of the features below since determining the failure zone was a very important part in slope behaviors. The end results provided the 2-D cross plot model for slope area, where the subsurface characteristics, such as boulders, saturated zones, weakness zones, bedrock, and other features, were well visualized. In general, three stages of processing were performed: geometry correction, graphical, and datum point analysis and clustering. In geometry correction, the datum point for each result was re-arranged based on the position of x and y that fall in the same position. The datum point, which does not fall at the same position of x and y, were removed from the data. After all the data were ensured to be precisely (or nearly) the same with the position of x and y locations for both data, the graphical and datum point analyses were performed to plot the resistivity-velocity graph. At this stage, the data, with a value much bigger than average range data, were removed. Clustering was performed to sort the data based on the ranged values of the quadrant from the scattered plot of the resistivity-velocity graph. Four quadrants (Q1-Q4) were assigned and have been sorted accordingly based on the quadrant limit of the resistivity and velocity. Each quadrant represents the different values with different characteristics, as shown previously in Figure 4a. In the clustering method, the sorted option in Microsoft Excel was utilized to categorize all of the data based on the quadrant ranged. Figure 4b shows the resistivity-velocity graph based on the quadrant range. All of the data for each quadrant were output as a (.xlsx) file and plotted as a 2-D model, using a post-map function in Surfer8 tools. The model was plotted accordingly to positions x and y of the data. Figure 5 shows the workflow of the cross-plot analysis. In this research, the cross-plot analysis was implemented to integrate both resistivity and seismic refraction results, which were in line with each other. The goal of this analysis is to provide detailed subsurface images so as to enhance the visualization of the features below since determining the failure zone was a very important part in slope behaviors. The end results provided the 2-D cross plot model for slope area, where the subsurface characteristics, such as boulders, saturated zones, weakness zones, bedrock, and other features, were well visualized. In general, three stages of processing were performed: geometry correction, graphical, and datum point analysis and clustering. In geometry correction, the datum point for each result was re-arranged based on the position of x and y that fall in the same position. The datum point, which does not fall at the same position of x and y, were removed from the data. After all the data were ensured to be precisely (or nearly) the same with the position of x and y locations for both data, the graphical and datum point analyses were performed to plot the resistivity-velocity graph. At this stage, the data, with a value much bigger than average range data, were removed. Clustering was performed to sort the data based on the ranged values of the quadrant from the scattered plot of the resistivity-velocity graph. Four quadrants (Q1-Q4) were assigned and have been sorted accordingly based on the quadrant limit of the resistivity and velocity. Each quadrant represents the different values with different characteristics, as shown previously in Figure 4a. In the clustering method, the sorted option in Microsoft Excel was utilized to categorize all of the data based on the quadrant ranged. Figure 4b shows the resistivity-velocity graph based on the quadrant range. All of the data for each quadrant were output as a (.xlsx) file and plotted as a 2-D model, using a post-map function in Surfer8 tools. The model was plotted accordingly to positions x and y of the data. Figure  5 shows the workflow of the cross-plot analysis.

Volumetric Calculation of Data
The volume of mass movement was estimated based on the integration of 2-D resistivity and seismic refraction methods. In this part, the location of the potential slip surface identified in 2-D cross-plot models was digitized for further analysis of volumetric calculation. The mass movement after the landslide was calculated based on three methods; Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule, in determining the volume in Surfer8 as follows in Equation (1). Mathematically, the volume under the function f (x,y) is defined as; In Surfer, the software computed the values by first integrating the X (column) to produce an area under the individual rows with following the integrating of Y (rows) for the final volume. The software uses three classical numerical integration algorithms, which is Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule [72]. The estimated volume of the mass movement was calculated with Equation (2) based on the calculation of the volume of total area subtracted with the volume under the failure plane. Figure 6 shows the graphical illustration of volumetric mass movement calculations.
VT = volume of total area VS = volume below the failure plane/sliding plane VMM = volume of mass movement

Results
In this section, the result of 2-D resistivity and seismic refraction were discussed with the integration of the data using cross-plot analysis. The output model of 2-D cross-plot revealed the vulnerability of the subsurface and the information from the model were used in estimated mass sliding. Figure 7a shows the inversion model of 2-D resistivity for HR1 and seismic refraction of HS1.

Volumetric Calculation of Data
The volume of mass movement was estimated based on the integration of 2-D resistivity and seismic refraction methods. In this part, the location of the potential slip surface identified in 2-D cross-plot models was digitized for further analysis of volumetric calculation. The mass movement after the landslide was calculated based on three methods; Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule, in determining the volume in Surfer8 as follows in Equation (1). Mathematically, the volume under the function f (x,y) is defined as; In Surfer, the software computed the values by first integrating the X (column) to produce an area under the individual rows with following the integrating of Y (rows) for the final volume. The software uses three classical numerical integration algorithms, which is Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule [72]. The estimated volume of the mass movement was calculated with Equation (2) based on the calculation of the volume of total area subtracted with the volume under the failure plane. Figure 6 shows the graphical illustration of volumetric mass movement calculations.

Volumetric Calculation of Data
The volume of mass movement was estimated based on the integration of 2-D resistivity and seismic refraction methods. In this part, the location of the potential slip surface identified in 2-D cross-plot models was digitized for further analysis of volumetric calculation. The mass movement after the landslide was calculated based on three methods; Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule, in determining the volume in Surfer8 as follows in Equation (1). Mathematically, the volume under the function f (x,y) is defined as; In Surfer, the software computed the values by first integrating the X (column) to produce an area under the individual rows with following the integrating of Y (rows) for the final volume. The software uses three classical numerical integration algorithms, which is Trapezoidal Rule, Simpson's Rule, and Simpson's 3/8 Rule [72]. The estimated volume of the mass movement was calculated with Equation (2) based on the calculation of the volume of total area subtracted with the volume under the failure plane. Figure 6 shows the graphical illustration of volumetric mass movement calculations.
VT = volume of total area VS = volume below the failure plane/sliding plane VMM = volume of mass movement

Results
In this section, the result of 2-D resistivity and seismic refraction were discussed with the integration of the data using cross-plot analysis. The output model of 2-D cross-plot revealed the vulnerability of the subsurface and the information from the model were used in estimated mass sliding. Figure 7a shows the inversion model of 2-D resistivity for HR1 and seismic refraction of HS1.

Results
In this section, the result of 2-D resistivity and seismic refraction were discussed with the integration of the data using cross-plot analysis. The output model of 2-D cross-plot Appl. Sci. 2021, 11, 747 9 of 23 revealed the vulnerability of the subsurface and the information from the model were used in estimated mass sliding. Figure 7a shows the inversion model of 2-D resistivity for HR1 and seismic refraction of HS1. The 2-D resistivity shows resistivity values ranging from 1-10,000 Ωm with depth coverage of about 80 m. The orientation of these lines is parallel to the landslide area, which is oriented from the East to West (E-W) direction. The resistivity result indicated two different contrast zones of the resistivity layer. The resistivity zones between 100 and 1200 Ωm indicate as a weathered zone, with a depth 20 m from the surface. In this profile, this layer is categorized as highly weathered zones, with the presence of boulders with resistivity values ranging between 1200 and 3500 Ωm. A saturated zone with a resistivity of <100 Ωm was identified at depth >25 m. The fractures zone was identified due to contrast in resistivity values, as shown in the dotted line. The seismic refraction profile of HS1 shows three layer cases, where the first layer is indicated as a highly weathered zone with velocity values of <1200 m/s. The second layer with a velocity ranging of 1000-2000 m/s indicates the weathered zone, which represents the loose material or highly weathered bedrock. The third layer was indicated as bedrock at velocity values of >3800 m/s. contrast in resistivity values, as shown in the dotted line. The seismic refraction profile of HS1 shows three layer cases, where the first layer is indicated as a highly weathered zone with velocity values of <1200 m/s. The second layer with a velocity ranging of 1000-2000 m/s indicates the weathered zone, which represents the loose material or highly weathered bedrock. The third layer was indicated as bedrock at velocity values of >3800 m/s. Figure 7b shows the inversion model of 2-D resistivity for HR2 and seismic refraction for HS2. In 2-D resistivity, two distinct zones were identified based on the contrast values. Highly weathered colluviums zones ranged from 100-1200 Ωm with a depth of 30-70 m from the surface. The contrast resistivity zone at a distance of 100 m was suspected as a fracture zone for this profile (dotted line). The saturated zone with a resistivity of <100 Ωm was identified at a depth >20 m. For the seismic refraction profile, in line with the resistivity, showed the first layer dominate with velocity values ranging <1200 m/s at a depth of 5-20 m. The second layer with granitic bedrock identified velocity values of >3800 m/s. Figure 7c shows the inversion model of 2-D resistivity for HR4 and seismic refraction of HS3. For HR4, two distinct zones were identified as same as the lines HR1 and HR2, where colluvium was identified at resistivity values of 100-1200 Ωm with a depth of 20-70 m, while the saturated zone identified at values of <100 Ωm. The bedrock was identified at resistivity values >7000 Ωm with the fracture zones mark with a dotted line in this profile. The seismic profile of HS3 shows the three layer cases, where the first layer identified at a velocity of <600 m/s while the second layer at 1250-1500 m/s, with a depth of 5-10 m. The third layer, with velocity values ranging from 2900-3300 m/s, identified as the weathered bedrock for this profile, with a depth of 20-30 m. Figure 7d-f shows the inversion model of 2-D resistivity of lines, HR3, HR5, and HR6. In HR3, the highly weathered zone was identified at resistivity values of <300 Ωm with fracture as indicted the contrast (low and high) of resistivity region at a distance of 130 m. In a profile of HR5 and HR6, two distinct zones were classified as colluviums at resistivity values of 100-1200 Ωm at a depth of 20-70 m while secondary zones were indicated as granitic bedrock, with resistivity values of >7000 Ωm. Bedrock was identified with resistivity values >7000 Ωm at depth >40 m.

Discussion
The result of 2-D resistivity and seismic profiles show the characteristic of the subsurface material for the study area. The failure zone was identified for this area as the existence of the fracture zone with the boulder, and may trigger the occurrence of the mass movement. The low resistivity region, indicated as a highly saturated zone, contributed to the factor of the occurrence mass movement as the soil profile dominated by the sand and sandy silt. The continuous infiltration of surface water and heavy rain as it is seeping through to the fracture area, and is accumulated, may  Highly weathered colluviums zones ranged from 100-1200 Ωm with a depth of 30-70 m from the surface. The contrast resistivity zone at a distance of 100 m was suspected as a fracture zone for this profile (dotted line). The saturated zone with a resistivity of <100 Ωm was identified at a depth >20 m. For the seismic refraction profile, in line with the resistivity, showed the first layer dominate with velocity values ranging <1200 m/s at a depth of 5-20 m. The second layer with granitic bedrock identified velocity values of >3800 m/s. Figure 7c shows the inversion model of 2-D resistivity for HR4 and seismic refraction of HS3. For HR4, two distinct zones were identified as same as the lines HR1 and HR2, where colluvium was identified at resistivity values of 100-1200 Ωm with a depth of 20-70 m, while the saturated zone identified at values of <100 Ωm. The bedrock was identified at resistivity values >7000 Ωm with the fracture zones mark with a dotted line in this profile. The seismic profile of HS3 shows the three layer cases, where the first layer identified at a velocity of <600 m/s while the second layer at 1250-1500 m/s, with a depth of 5-10 m. The third layer, with velocity values ranging from 2900-3300 m/s, identified as the weathered bedrock for this profile, with a depth of 20-30 m. Figure 7d-f shows the inversion model of 2-D resistivity of lines, HR3, HR5, and HR6. In HR3, the highly weathered zone was identified at resistivity values of <300 Ωm with fracture as indicted the contrast (low and high) of resistivity region at a distance of 130 m. In a profile of HR5 and HR6, two distinct zones were classified as colluviums at resistivity values of 100-1200 Ωm at a depth of 20-70 m while secondary zones were indicated as granitic bedrock, with resistivity values of >7000 Ωm. Bedrock was identified with resistivity values >7000 Ωm at depth >40 m.

Discussion
The result of 2-D resistivity and seismic profiles show the characteristic of the subsurface material for the study area. The failure zone was identified for this area as the existence of the fracture zone with the boulder, and may trigger the occurrence of the mass movement. The low resistivity region, indicated as a highly saturated zone, contributed to the factor of the occurrence mass movement as the soil profile dominated by the sand and sandy silt. The continuous infiltration of surface water and heavy rain as it is seeping through to the fracture area, and is accumulated, may trigger the instability of the slope, as the permeability, and the porosity of the soil is high. The disturbance of the soil properties with the variation of porosity and permeability increased the effective stress of the soil, which is able to trigger the mass movement for the instability slope.
Previous research conducted by Ghazali et al., in residual soil for the granitic area, shows the same characteristic in this area where the fracture and highly weathered mass contribute a major causative factor for slope failure [73]. These features lead to the movement of the seepage and underground water through it. The groundwater flow was determined as a plane of weakness for that research. In this case, the hillslope soil condition, with loose material dominant for all profiles, which able to trigger the slope failure under the continuous infiltration of the surface water above the water table. The occurrence of the fracture zones provided a wider zone for surface water to seep through to the soil and accumulate. The unstable condition, with loose materials and increased accumulation of water, triggered the occurrence of slope failure. The failure zone was identified for this area as the existence of the fracture zone with the boulder, and may trigger the occurrence of the mass movement. The presence of the boulder can be a key factor for the instability due to change and removing the rock that holds the slope in place [74,75]. Other than that, highly fractured zone and water-conducting zones are indicated by the presence of low resistivity zones, within the granite body. The thin layer of residual soil due to domination of the granite body at the slope face triggered the failure or landslide events. Figure 8 shows the water circulation pathway of the Gapis River associated with sediments in this survey area. Low resistivity values (HR6-HR2-HR1) are observed as the water circulation quite deep beneath this area associated with water from the Gapis River directed towards northeast to southwest (NE-SW). Other than, in HR3-HR4-HR5, the water circulation may be relatively shallow at this location at a depth of 5-20 m. The survey enables to characterize the unstable region along with the profiles as the existence of fracture and scattered boulders. The resistivity value of <1200 Ωm at a depth of <20 m represents the weathered bedrock. Relatively, the unaltered bedrock generally has a high resistivity value of >7000 Ωm and the elongate zone of lower resistivity mostly of the range of 500-1600 Ωm parallel to topography interpreted as altered zones. The seismic refraction survey shows the bedrock was identified at velocity values of >3800 m/s. The first layer dominates with velocity values of <1200 m/s represented as highly weathered zones. The lower velocity values show that this zone had a high potential for instability to occur. The broad zones of low resistivity values (50-300 Ωm) with the appearance of high resistivity values of (1200-2500 Ωm) are associated with the fracture or failure zone along with profiles. Low resistivity region (<100 Ωm), as the saturated zones, may correlate with high permeability of water-bearing fractures [76,77].
Appl. Sci. 2020, 10, x FOR PROOFREADING 11 of 24 body. The thin layer of residual soil due to domination of the granite body at the slope face triggered the failure or landslide events. Figure 8 shows the water circulation pathway of the Gapis River associated with sediments in this survey area. Low resistivity values (HR6-HR2-HR1) are observed as the water circulation quite deep beneath this area associated with water from the Gapis River directed towards northeast to southwest (NE-SW). Other than, in HR3-HR4-HR5, the water circulation may be relatively shallow at this location at a depth of 5-20 m. The survey enables to characterize the unstable region along with the profiles as the existence of fracture and scattered boulders. The resistivity value of <1200 Ωm at a depth of <20 m represents the weathered bedrock. Relatively, the unaltered bedrock generally has a high resistivity value of >7000 Ωm and the elongate zone of lower resistivity mostly of the range of 500-1600 Ωm parallel to topography interpreted as altered zones. The seismic refraction survey shows the bedrock was identified at velocity values of >3800 m/s. The first layer dominates with velocity values of <1200 m/s represented as highly weathered zones. The lower velocity values show that this zone had a high potential for instability to occur. The broad zones of low resistivity values (50-300 Ωm) with the appearance of high resistivity values of (1200-2500 Ωm) are associated with the fracture or failure zone along with profiles. Low resistivity region (<100 Ωm), as the saturated zones, may correlate with high permeability of water-bearing fractures [76,77]. In order to validate the geophysical results, a borehole record was used as supported data in providing realistic interpretations. In this study area, only two boreholes were in line with the geophysical survey: BH9 and BH10. Figure 9 shows the correlation of borehole and geophysical results for line HR1 and HS1 with BH10 and HR2 and HS2 with BH9 at a distance of 225 m and 175 m, respectively. In Figure 9a, the borehole record, BH10, shows the record with investigation depth about 8.5 m from the surface. The result showed a much lower N-value at a depth range of 0-2.5 m that shows a low stiffness of soil for this particular point. At this depth range, both 2-D resistivity and seismic refraction show no data availability due to the blanking effect during processing. As the borehole reaches at depth >3 m, the record shows the drastic increases of the trend where N-value of 50 blows were recorded until the end of the borehole. From this, N-values of 50 blows show that the compacted or stiffer/very dense soil was reached as the rock head was recorded at a depth of 5.5 m from the surface. In addition, both 2-D resistivity and seismic refraction shows convincing correlation at a depth of >3 m, the resistivity shows the value of 1600 Ωm consistently until the depth of 5 m, and increases >3500 Ωm at depth >5 m. The proportional increases of the resistivity values indicate the In order to validate the geophysical results, a borehole record was used as supported data in providing realistic interpretations. In this study area, only two boreholes were in line with the geophysical survey: BH9 and BH10. Figure 9 shows the correlation of borehole and geophysical results for line HR1 and HS1 with BH10 and HR2 and HS2 with BH9 at a distance of 225 m and 175 m, respectively. In Figure 9a, the borehole record, BH10, shows the record with investigation depth about 8.5 m from the surface. The result showed a much lower N-value at a depth range of 0-2.5 m that shows a low stiffness of soil for this particular point. At this depth range, both 2-D resistivity and seismic refraction show no data availability due to the blanking effect during processing. As the borehole reaches at depth >3 m, the record shows the drastic increases of the trend where N-value of 50 blows were recorded until the end of the borehole. From this, N-values of 50 blows show that the compacted or stiffer/very dense soil was reached as the rock head was recorded at a depth of 5.5 m from the surface. In addition, both 2-D resistivity and seismic refraction shows convincing correlation at a depth of >3 m, the resistivity shows the value of 1600 Ωm consistently until the depth of 5 m, and increases >3500 Ωm at depth >5 m. The proportional increases of the resistivity values indicate the changes of soil stiffness as the values >3500 Ωm indicate the weathered bedrock/boulder. In the seismic profile, at depth ranges of 3-5 m, the velocity is represented by the values of 800-1400 m/s and increases >1800 m/s at depth >5 m. The velocity of 800-1400 m/s indicates the medium stiffness of the soil, which probably shows the weathered layer for the material as the velocity increases >1800 m/s corresponding with the N-value of 50 blows, as rock head was found at a depth of 5.5 m. The velocity shows the increase in trends at a depth of >8.5 m as investigation depth cover ups to 25 m with velocity ranges of 1800-4200 m/s. The wide ranges, with proportional trends of the velocity, show the increased strength of soil material, as the bedrock was interpreted with values of >3400 m/s.
In other profiles of borehole records of BH9, only line HR2 was in line at a distance of 175 m as HS1 showed limited survey distance, as shown in Figure 9b. The borehole record shows a deeper depth of investigation, up to 33 m from the surface. However, at this particular distance, the resistivity shows no data availability due to signal coverage at a depth >15 m. Only the depth range of 1.5-15 m was used to correlate with N-values data. From the record of BH9, it was clearly seen the variation of N-values that represents the heterogeneous condition with different soil stiffness at depth <10 m. The N-values show a lower number of blows from the range 1-10 at this depth, which indicates low stiffness or represents loose material for this depth range, as the resistivity shows the variation of values of 150-800 Ωm. As the depth increases, at depth 12 m the borehole record shows the N-value of 50 blows for this particular depth, and decreases at a depth of 13-15 m, where N-values were recorded <5 number of blows. The upsurge decreases of N-values from 50, which indicate this layer consists of very poor conditions that shows very loose material was underlined. However, the N-values show the gradual increases from 18-50 of N-values, up to a depth of 27 m from the surface. In resistivity profiles, the values show the gradual changes of resistivity distribution from a range of 500-1200 Ωm up to the depth of 15 m, as the data were not available above that depth. The borehole record shows a consistent number of blow at depth 27-33 m with values of 50, which indicate the hard layer (medium dense) for these profiles as no indicator was marked as bedrock compared to others boreholes. The schematic diagram of soil classification, [78] based on the distribution of N-value for BH7-BH10 are shown in Figure 9c. Noted that BH7-BH8 was located in the landslide area. Generally, the profile shows very loose-loose material at a depth range of 3-10 m. This shows the soil profile was considered as highly weathered with the presence of the boulder.
Appl. Sci. 2020, 10, x FOR PROOFREADING 12 of 24 increases >3500 Ωm at depth >5 m. The proportional increases of the resistivity values indicate the changes of soil stiffness as the values >3500 Ωm indicate the weathered bedrock/boulder. In the seismic profile, at depth ranges of 3-5 m, the velocity is represented by the values of 800-1400 m/s and increases >1800 m/s at depth >5 m. The velocity of 800-1400 m/s indicates the medium stiffness of the soil, which probably shows the weathered layer for the material as the velocity increases >1800 m/s corresponding with the N-value of 50 blows, as rock head was found at a depth of 5.5 m. The velocity shows the increase in trends at a depth of >8.5 m as investigation depth cover ups to 25 m with velocity ranges of 1800-4200 m/s. The wide ranges, with proportional trends of the velocity, show the increased strength of soil material, as the bedrock was interpreted with values of >3400 m/s. In other profiles of borehole records of BH9, only line HR2 was in line at a distance of 175 m as HS1 showed limited survey distance, as shown in Figure 9b. The borehole record shows a deeper depth of investigation, up to 33 m from the surface. However, at this particular distance, the resistivity shows no data availability due to signal coverage at a depth >15 m. Only the depth range of 1.5-15 m was used to correlate with N-values data. From the record of BH9, it was clearly seen the variation of N-values that represents the heterogeneous condition with different soil stiffness at depth <10 m. The N-values show a lower number of blows from the range 1-10 at this depth, which indicates low stiffness or represents loose material for this depth range, as the resistivity shows the variation of values of 150-800 Ωm. As the depth increases, at depth 12 m the borehole record shows the Nvalue of 50 blows for this particular depth, and decreases at a depth of 13-15 m, where N-values were recorded <5 number of blows. The upsurge decreases of N-values from 50, which indicate this layer consists of very poor conditions that shows very loose material was underlined. However, the Nvalues show the gradual increases from 18-50 of N-values, up to a depth of 27 m from the surface. In resistivity profiles, the values show the gradual changes of resistivity distribution from a range of 500-1200 Ωm up to the depth of 15 m, as the data were not available above that depth. The borehole record shows a consistent number of blow at depth 27-33 m with values of 50, which indicate the hard layer (medium dense) for these profiles as no indicator was marked as bedrock compared to others boreholes. The schematic diagram of soil classification, [78] based on the distribution of Nvalue for BH7-BH10 are shown in Figure 9c. Noted that BH7-BH8 was located in the landslide area. Generally, the profile shows very loose-loose material at a depth range of 3-10 m. This shows the soil profile was considered as highly weathered with the presence of the boulder. The correlation of overview hydrogeological parameter and landslide initiation was presented using available rainfall data distribution based on the nearest rainfall station. As mentioned before, rainfall distribution in Malaysia is mainly influenced by two monsoons seasons, namely southeast and northeast monsoons, from May-September and November-March. However, the highest precipitation occurs during the transition period between the monsoon seasons, which are known as inter-monsoon of March-April and October-November, as shown in Figure 10a. The inter-monsoon seasons show the highest domination of rainfall distribution of 361-242 mm and 273.2-267 mm for Mach-April and October-November, respectively. This prolonged rainfall shows as a factor in building up to the mechanism of landslides. However, in a rainfall-induced mechanism, it is widely debated regarding the antecedent rainfall, referring to the rainfall in the days immediately preceding landslide events, while triggering rainfall represented as the rain that falls in a day of landslide occurrence [79,80]. Previous researcher has well accepted the effects of antecedent rainfall on slope instability; however, the duration needed for cumulative rainfall needs to be considered [81]. In this study, 15-days rainfall data were used, suggested by Zezera et al., as antecedent rainfall for shallow landslide [82]. Figure 10b shows the overview of rainfall distribution for three months before the landslide event. The bar chart shows the rainfall amount upsurge trend from September to October and gradually increases from October to November. Based on the chart given, the rainfall is intense in October-November with a rainfall amount of 30-90 mm per day. This shows this area exhibits daily rainfall and leads to the progressive weathering of the soil. The cumulative rainfall data (Figure 10c

2-D Cross-Plot Model of Resistivity-Velocity
From the previous discussion, the resistivity of <1200 Ωm and velocity of <1200 m/s was identified as the threshold value for the material to undergo the mass movement. Note that this survey was acquired after the landslide event, where line HR2 and HS2 and HR4 and HS3 was crossing the landslide boundary while HR1 and HS1 outside formed the landslide boundary. From the graph ( Figure 11) the resistivity-velocity ranged for line HR1 and HS1 is much larger as compared to the other two lines as this line is located outside from the boundary line of a landslide, which The variation of weak zones indicated by the geophysical results show the presence of subsurface materials, which exhibit the additional weathering process, which tends to produce poor rock mass quality. The extreme weather climate leads to the progressive weathering processes due to the abundance of rainfall and moistures as accelerating the transformations of homogenous subsurface materials in heterogeneous conditions [83]. This process allows the water to flow and weaken the subsurface geo-materials according to the current topography of surface direction. The intense precipitation will contribute to surface runoff, as the increasing groundwater level in saturated conditions, as well as increasing soil mass become the additional trigger factors for failure. The intense and continuous precipitation causes the pore pressure to increase as the groundwater rises rapidly. The continuous process reduces the shearing resistance of geo-materials and contributes to slope failure. The reducing shear strength and stiffness of the materials (due to new formations of fractures, or due to the extension as incipient fractures lose the tensile strength), causes the discontinuity of rock to become weak. The presence of this situation causes the infiltration of water to flow through the surface crack and failure structure easily. The process intensively weakens or reduces the geo-materials of the subsurface, especially by chemical alteration due to the weathering process as it takes place via water movement throughout mass materials [84]. The weathered fracture materials can disintegrate and decompose into fine-grained materials, such as sand, silt, and silty sand, which later fill the existing joint or fractures. Thus, the process leads to variations in subsurface velocity, as infilling material shows in high velocity, as compared to loose material.

2-D Cross-Plot Model of Resistivity-Velocity
From the previous discussion, the resistivity of <1200 Ωm and velocity of <1200 m/s was identified as the threshold value for the material to undergo the mass movement. Note that this survey was acquired after the landslide event, where line HR2 and HS2 and HR4 and HS3 was crossing the landslide boundary while HR1 and HS1 outside formed the landslide boundary. From the graph ( Figure 11) the resistivity-velocity ranged for line HR1 and HS1 is much larger as compared to the other two lines as this line is located outside from the boundary line of a landslide, which represents the originality of subsurface condition before the occurrence of the landslide. This explains the high resistivity-velocity range of the line.

2-D Cross-Plot Model of Resistivity-Velocity
From the previous discussion, the resistivity of <1200 Ωm and velocity of <1200 m/s was identified as the threshold value for the material to undergo the mass movement. Note that this survey was acquired after the landslide event, where line HR2 and HS2 and HR4 and HS3 was crossing the landslide boundary while HR1 and HS1 outside formed the landslide boundary. From the graph ( Figure 11) the resistivity-velocity ranged for line HR1 and HS1 is much larger as compared to the other two lines as this line is located outside from the boundary line of a landslide, which represents the originality of subsurface condition before the occurrence of the landslide. This explains the high resistivity-velocity range of the line. The 2-D cross-plot model was an enhancement work for cross-plot analysis in integrating different geophysical parameters. The goal of this work is to visualize and provide comprehensive subsurface images as well as to integrated both parameters of resistivity and velocity to be more precisely. In this study, the method was applied, and the final model of the 2-D cross plot was presented in Figure 12. Based on the previous discussion, this work was divided into three different stages of processing until the final data. Four quadrants were divided, representing four different characteristics of the subsurface, which are Q1 (≤1200 Ωm, ≤1200 m/s), Q2, (≥1200 Ωm, ≤1200 m/s), Figure 11. Graph correlation of resistivity-velocity for the same X-Y location for line: (a) HR1 and HS1; (b) HR2 and HS2; (c) HR4 and HS3. The 2-D cross-plot model was an enhancement work for cross-plot analysis in integrating different geophysical parameters. The goal of this work is to visualize and provide comprehensive subsurface images as well as to integrated both parameters of resistivity and velocity to be more precisely. In this study, the method was applied, and the final model of the 2-D cross plot was presented in Figure 12. Based on the previous discussion, this work was divided into three different stages of processing until the final data. Four quadrants were divided, representing four different characteristics of the subsurface, which are Q1 (≤1200 Ωm, ≤1200 m/s), Q2, (≥1200 Ωm, ≤1200 m/s), Q3, (≥1200 Ωm, ≥1200 m/s), and Q4 (≤1200 Ωm, ≥1200 m/s). The sorting analysis is applying to extract the values according to its depth and position that inline each other. Quadrant 1 (Q1) compromise of low resistivity and velocity values of ≤1200 Ωm and ≤1200 m/s respectively. From the preliminary results of 2-D resistivity and seismic refraction, this value plays as a threshold factor for the occurrence of landslides event. From the 2-D cross-plot model, the information of the highly weathered zone was retrieved and used to estimate the volume of mass sliding. The estimated volume of the mass sliding was calculated using Equation (2), where the calculation of the volume of total area is subtracted from the volume under the failure plane. The Kriging interpolation methods were used in generating the model. From the previous figure of the study area (Figure 1), the boundary of the landslide area was digitized and plotted, as shown in Figure 13a.
In order to identify the volume, m 3 of mass sliding the Universal Transverse Mercator (UTM) coordinate in meter, m unit was used with the elevation map (m unit) of the study area. It is very important to ensure that the unit used for all parameters is the same to avoid wrong interpretations. From the elevation map, the elevation map of the landslide area was digitized and plotted, as shown in Figure 13b. The elevation map under the failure plane was plotted (Figure 13c) after being digitized from geophysical results. This map represents the surface under the mass sliding. From this data, the volume was computed using the volume function in Surfer tools. The volume in the landslide area was extracted as the volume of the surface (Figure 14a) and the volume under the failure plane (Figure 14b). The volume difference between the two figures of 15a and b represent the volume of mass sliding. Table 2 shows a summary of the calculation.  , marked with a solid black line. In Q1, the region represents the highly weathered zone with lower resistivity and velocity values. This zone indicates the potential zone for mass movements or triggered zone for slope failure. The Q2 region at depth <5 m (solid rectangle box), shows low velocity and high resistivity values, which indicates the presence of moderate weathered material for this profile. The highest region represents by Q3 may indicate the presence of hard materials/boulder, which is aided as an additional factor for the occurrence of landslide, besides the movement of water flow due to saturated zones. The Q4 region dominates with a resistivity value of 1-1200 Ωm and a velocity of >1200 m/s.
In Figure 12b,c, it shows the 2-D resistivity, seismic refraction, and 2-D cross-plot model for HR2 and HS2 and HR4 and HS3, respectively. As previously mentioned, this line was crossing the landslide boundary, which is generally the top surface cover with landslide material. From 2-D resistivity results, the saturated and highly weathered zones were identified at a resistivity value of <100 Ωm and 100-1200 Ωm, respectively. The seismic refraction profile shows the same characteristics where a low velocity zone identifies with a velocity of <1200 m/s. This zone mostly probably covers with loose material, such as sand with depth <5 m, as shown in the 2-D cross-plot model. The model reveals, at a depth of <5 m, the layer dominated with values from Q1 region, which classified as low velocity and resistivity values, indicating poorly compacted soil or the existence of higher porosity material, such as sand. The Q4 region looks dominant for this profile, which indicates the weathered zone; the region of Q2 at a depth of <5 m (solid rectangle box) shows low velocity but high resistivity, and may indicate the presence of a boulder at shallow subsurface, or represent as moderate weathered material for this profile. The Q3 region, with high resistivity and velocity, may represent the hard material at a depth of >15 m. The suspected failure plane was marked with a solid black line.
From the 2-D cross-plot model, the information of the highly weathered zone was retrieved and used to estimate the volume of mass sliding. The estimated volume of the mass sliding was calculated using Equation (2), where the calculation of the volume of total area is subtracted from the volume under the failure plane. The Kriging interpolation methods were used in generating the model. From the previous figure of the study area (Figure 1), the boundary of the landslide area was digitized and plotted, as shown in Figure 13a.  In order to identify the volume, m 3 of mass sliding the Universal Transverse Mercator (UTM) coordinate in meter, m unit was used with the elevation map (m unit) of the study area. It is very important to ensure that the unit used for all parameters is the same to avoid wrong interpretations. From the elevation map, the elevation map of the landslide area was digitized and plotted, as shown in Figure 13b. The elevation map under the failure plane was plotted (Figure 13c) after being digitized from geophysical results. This map represents the surface under the mass sliding. From this data, the volume was computed using the volume function in Surfer tools. The volume in the landslide area was extracted as the volume of the surface (Figure 14a) and the volume under the failure plane (Figure 14b). The volume difference between the two figures of 15a and b represent the volume of mass sliding. Table 2 shows a summary of the calculation.    The parameters above were generated by the Surfer volume function as three methods of volume calculation were presented. The final volume is represented by positive volume, which is 167,904.81 m 3 , and indicates the volume of mass sliding. The negative and net volume is neglected since the input data of upper and lower surfaces (Figure 14c) did not overlap each other.

Conclusions
The geophysical survey was successfully applied to identify the subsurface characteristics of slope failure. Through the study, weathering and distressing have weakened the soil profiles as velocity, <1200 m/s, and resistivity, <1200 Ωm, defined as the threshold for failure to occur. The weak zones, suspected to be landslides with sliding, were determined from the 2-D cross plot model. From the results, it reveals that the failure plane occurred at a depth of <5 m (Quadrant 1, Q1).
The factor of failure in this area was suspected of loose material (sand), the existence of boulder, and accumulated saturated zone at the bottom part of the slope, and may trigger the occurrence of the landslide event. Major landslides at local slope failure occur due to water saturation exceeding a critical limit in certain parts of the slope [47]. The continuous infiltration of rainwater increases the pore-air pressure and decreases the matric suction, which compromises the slope stability [44,45]. The process causes a change in effective stress, which produces the variations in porosity [45]. The infiltration process of near-surface flow increases the pore water pressure as created the stress path to shift horizontally to intersect at the failure envelope. The continuous infiltrated process will trigger the occurrence of slope failure. Other than that, for unsaturated/loose materials, the decreasing of suction force and collapsing of coupled volumetric may be involved in the failure process [48]. The sliding mass can also be estimated based on the integrated values from both results. The identification of the failure plane is vital in identifying the failure zone for all profiles. This survey gave valuable information regarding the slope instability investigation, and it was concluded that it was successfully integrated to evaluate the slope failure for this area.
In a nutshell, the cross-plot analysis shows a qualitative method in integrating the different geophysical parameters. The geophysical approaches give comprehensive images to characterize the landslide properties, with validation from geotechnical and hydrogeological parameters providing substantiate interpretation of geophysical results. The integrated method of cross-plot analysis with the 2-D cross-plot model development shows an informative method in characterization and enhances subsurface images.