Mapping of Landslide Potential in Payung, Batu City, Indonesia, Using Global Gravity Model Plus (GGMplus) Data as Landslide Mitigation

Abstract


Introduction
Landslides can cause damage to economic assets such as buildings and productive land (Priyono et al., 2020).One of the important buildings that became a means of access to transportation is the highway.On highways, areas that experience shifts, landslides, or subsidence are often found; some of these are forms of ground movement.Ground movement can occur due to the displacement of soil or rock masses in a vertical, flat, or oblique direction from their original position.This movement includes the movement of soil material in the form of rock, embankment material, soil, or mixed material.The most common problem is landslides when entering the rainy season; in Indonesia, landslides are spread in almost all provinces with high risk (Hasan et al., 2022;Karimah et al., 2022).
Several areas in Batu City often experience landslides when entering the rainy season; landslides occur at several central points of Batu City.Based on the data we obtained from the Regional Disaster Management Agency (BPBD) Batu City, one of the spots that experience landslides every year is the Payung area on Trunojoyo Street, Songgokerto Village, Batu City (BPBD Batu City, 2022); throughout 2022, 11 landslide events were recorded in this village.This street connects several areas in East Java, such as Malang Regency, Kediri City, and Jombang Regency.If this area experiences a landslide, the access road will be closed, forcing road users to look for another route that is further away.If viewed directly, the condition of the road is on a mountain slope, so the potential for landslides is also relatively high (Duc-Nghiem and Viet-Phuong 2021).When the slope's surface is exposed to heavy and prolonged rains, the potential for landslides will be even more significant (Lias et al., 2022).
The solution to this problem is that it is necessary to carry out studies related to areas that have a high level of landslide vulnerability.One of the methods commonly used for preliminary studies of landslides is the geophysical method, such as the gravity method.Based on previous research, Jacob et al. (2018) used the gravity method to map clay basins on cliffs for landslide hazard assessment and research.Other researchers, such as Legowo et al. (2022), use the gravity method to model fault structures in determining landslide potential.The gravity method is effective for detecting bedrock structures (Al-Khafaji et al., 2021); gravity data information can be used efficiently in various geological problems related to the exploration of the earth's crust, such as in natural disaster mitigation (Suprianto et al., 2021).The gravity method is advantageous for initial surveys because it can provide detailed information about the geological structure and rock density contrast.Measurement of the gravitational field aims to find out local density information, including the formation process, and study the irregularities of the earth (Telford et al., 1990).The working principle of the gravity method is to measure gravity anomalies caused by density differences between the subsurface of rocks (Xue et al., 2022).Gravity anomalies reflect lateral density fluctuations in the earth's crust, which is essential in investigating subsurface layers (Chen et al., 2022).Therefore, the gravitational method is commonly used to map the density contrast under the surface (Dobrin and Savit, 1960;Peng and Liu, 2021).Differences in rock layers that are quite contrasting will become slip planes that will trigger landslides.The bedrock or lava rock with a large density tends to remain fixed and not experience movement.When above the bedrock, there is a layer of rock with a lower density and higher porosity; it will become a place for water to accumulate.The impermeable layer will become a sliding zone so that the layer above it will move in the direction of the slope due to water or hydrostatic pressure (Ling et al., 2016).
Generally, the density contrast between clay and bedrock is quite large, where clay rock has a density of 1.63 mg/m 3 while bedrock is above 2.35 mg/m 3 (Telford et al., 1990); this can be detected based on the measured gravity anomaly.Gravity anomaly is defined as the difference between the observed gravity value and the value predicted by the model (Rambat et al., 2021).The gravity anomaly value is related to the rock density value; the higher the gravity anomaly value, the higher the rock density value.So, a low density value is indicated as a cavity or rock that can hold water.In the case of landslides, the difference in rock density is a reference in investigations using the gravity method.Low density areas are indicated as rocks with high porosity to hold water.If the rock layers are saturated with water, it will accelerate the process of landslides.This study aims to map areas potentially prone to landslides on the Payung tourism road in Batu City, Indonesia.The results of this research can be an initial reference for further research, such as using another geophysical method and geotechnical methods.

Materials and Methods
The data collection location was in the Payung tourism area in Songgokerto Village, Batu City, East Java, Indonesia.The research survey design is shown in Fig. 1; it shows a major road that connects several (red colour) locations in East Java, such as Malang Regency, Kediri City, and Jombang Regency.The gravitational field can be measured with a gravimeter, magnetometer, or satellite observation data (Rambat et al., 2021).In this research, gravity anomaly measurements were carried out using Global Gravity Model Plus (GGMplus) satellite data with a distance between data points of 220 m.GGMPlus refers to a combination of data from the Gravity Recovery and Climate Experiment (GRACE), Gravityfield and steady-state Ocean Circulation Explorer (GOCE), Earth Gravity Model (EGM 2008) combined with forward modeling of topographical gravity based on spherical harmonic expansion to improve resolution (Hirt et al., 2013;Akhadi et al., 2021;Legowo et al., 2022).The interpretation of gravity data needs to be done by an inversion process, commonly by performing repeated forward modeling calculations (Biswas, 2020;Wang et al., 2021).According to (Suprianto et al., 2021), data from GGMPlus generally has perfect resolution.The data obtained are the position of latitude and longitude, altitude, and Free Air Anomaly (FAA) (Sutasoma et al., 2022).Measurement points are designed in a grid form so that they can represent each measurement point (Fig. 1).
The gravity method is measured on the value of the vertical component of the acceleration due to gravity (vector field) at a location.However, the earth's shape is not spherical, so there are variations in the acceleration value due to gravity for each place.Gravity data from GGMPlus can be obtained using the approximative method and then analyzed spectrally using the discrete Fourier technique to get a model of the degree of variation; the discrete Fourier technique is formulated by the following equation (Rexer and Hirt, 2015): (1) Based on equation 1, C is the spherical harmonic representation coefficient (SH), GM is the product of the gravitational constant times the mass of the earth, and (l) is the spherical harmonic degree.The GGMPlus equations are used to decay the gravitational signal down to 90,000-degree harmonics based on 3 billion GGMPlus data points.The analysis that has been carried out produces a 2D model with an error rate of 10% (Rexer and Hirt, 2015).The gravity data processing stage begins with converting readings into milligal units.The data is then processed to produce a Bouguer Anomaly value by calculating the average density (Akhadi et al., 2021).Data processing uses Microsoft Excel to calculate Free Air Correction, Bouguer Correction, and Terrain Correction.Things that can affect the value of the acceleration due to gravity are differences in degrees of latitude, differences in topographical conditions, the position of the earth in the solar system, variations in the density of rock masses below the earth's surface, and differences in elevation where measurements are made.The acceleration value due to gravity depends on the place of measurement on the distance to the center of the earth; therefore, it is necessary to make corrections if there are differences in the conditions of the point in the observation of these theoretical conditions.Several factors generally influence gravity readings, so it is necessary to make some corrections to remove all unnecessary gravitational effects to obtain a complete Bouguer anomaly (Jacob et al., 2018).
The results of the gravity data processing will obtain complete Bouguer anomaly data, which will be transformed into a flat plane and a Butterworth filter to separate regional and residual anomalies.Reduction to a flat plane is the process of bringing Bouguer anomaly values that are still affected by topography to a flat plane at a certain height (Sutasoma et al., 2022).Reduction of the regional anomaly from the Bouguer anomaly produces a residual anomaly representing the density distribution below the surface (Jacob et al., 2018).The residual anomaly data obtained are interpreted by considering the geological information of the study area and other supporting data.Residual anomaly data obtained aims to identify low-density anomalies.Anomalies with low density and geological information in the study area are used as a preliminary study to determine landslide-prone areas.The next interpretation process, the residual anomaly data, is combined with the landslide vulnerability map on the InaRISK portal, which was previously compiled by the National Agency for Disaster Countermeasure (BNPB).

Geological Condition
Based on the geological map of the Kediri Quadrangle (Fig. 2 a), the research area consists of 3 different formations, namely the Old Anjasmara Volcanics Formation (Qpat), Kawi-Butak Volcanics Formation (Qpkb), and Upper Quaternary Volcanics Formation (Qv (pandermen)).The Old Anjasmara Volcanics Formation (Qpat) was formed during the Quaternary period in the Early Pleistocene.The rock formation consists of volcanic breccia, lava, tuff, and dikes.The Kawi-Butak Volcanics Formation (Qpkb) was formed in the Quaternary Period in the Middle into Late Pleistocene age, consisting of volcanic breccias, lava, tuff, and lahars.Then, the Upper Quaternary Volcanics Formation (Qv (pandermen)) was formed in the Holocene that consists of volcanic breccia, tuff breccia, lava, and tuff (Santosa and Atmawinata 1992).The geological column of the research area can be seen in Fig. 2b.In the current conditions, around the research location, there are several volcanoes.To the west, there are Mount Kawi and Anjasmoro, which are no longer active, while to the north to the northeast, there is Mount Arjuno-Welirang, which was recorded as having experienced its last eruption in 1952.Based on observations by the Indonesian Ministry of Energy and Mineral Resources from June to August 2023, the measured seismicity was dominated by gusts and distant tectonics earthquakes.The gusts earthquake with an amplitude of 6 mm with a duration of 17.52 seconds and two distant tectonic earthquakes with an amplitude of 15-20 mm, S-P 32.29-38.28seconds and the duration was 70.01-122.36seconds (Kuswano, 2023).Most of the earthquakes or seismic activity produced by Mount Arjuna-Welirang cannot be felt.Hence, several volcanoes' activity has so far not influenced landslide events in the research location.There are several active faults at the research location, as can be seen in Fig. 3.In the south, there is the Blitar thrust fault which extends west-east at a distance of 40 km from the research location, to the southeast, there is the Turen fault which extends northwest-southeast at a distance of 50 km, to the north there is a fault on Mount Arjuno which extends northwest-southeast for a distance of 16 km, and to the northeast there is a fault trending southwest-northeast for a distance of 19 km.

Results
Fig. 4 is the result of the complete Bouguer Anomaly, which has an anomaly value range from 369.8 to 475.6 mgal for the highest to the lowest color scale indicated by blue-pink colors.The anomaly value shows the density in an area, where the lower density value is in the northeast with a value range 369.8-415.2mgal, which is marked in blue, and conversely, the higher density is in the southwest in pink, which has a value range 415.2-475.6 mgal.It is necessary to separate regional and residual anomalies using the Butterworth filter in the Oasis Montaj software to obtain regional and residual anomalies.The complete Bouguer map anomaly has a large gravity anomaly, probably because volcanic rocks dominate the study area.Based on Table 1, the density values in volcanic rocks have a fairly high average density when compared to alluvial areas (Telford et al., 1990).5a is a regional anomaly ranging from 370 to 475.4 mgal, and Fig. 5b is a residual anomaly ranging from -13.4 to 12 mgal.Based on the two results above, the difference between the two is that the regional anomaly is almost the same as the complete Bouguer anomaly, and this is because the regional anomaly is broad or regional.Whereas the residual anomaly looks more varied, this is because the anomaly is local.In addition, the residual anomaly is also shallow, so we prefer to use residual anomalies in the case of landslides because the observed landslides occurred at the ground surface.Based on Fig. 5a, the regional anomaly map has high anomaly values in the southwest and decreases towards the northeast with a value range 370-415 mGal.This anomaly pattern is in accordance with the Bouguer anomaly map, which has the same pattern where high anomalies are in the southwest, and low anomalies are in the northeast.Based on these findings, it can be interpreted that the anomaly that influences the regional anomaly trend originates from the Bouguer anomaly.Meanwhile, on the residual omaly map (Fig. 5b), the anomaly pattern looks more varied and different from that found on the complete and regional Bouguer anomaly map.This is because the residual anomaly is local, shallow, and unaffected by regional anomalies.In the case of landslides, we use the residual anomaly because the observed landslides occurred on a shallow ground surface.
Based on the results of research using the resistivity geoelectric method carried out by Agung et al. (2023), it is known that the lava rocks that become the bedrock in the study area are at a depth of 25 m at Old Anjasmara Volcanics Formation and 50 m at Upper Quaternary Volcanics Formation.Lava rock has a higher density value than clay rock, so the meeting between two types of rock layers with contrasting densities will become a slip surface and be the main trigger for ground movement.(Jacob et al., 2018).Based on the residual anomaly map, it is known that the density of the study area is divided into three zones, namely low, medium, and high.The division of this zone is determined from the residual anomaly value where the low gravity ranges from -13.4 to -5 mgal, medium -5 to 5 mgal, and high from 5 to 12 mgal.The low value shows the area with the lowest density, marked in blue-green, while the medium-density value is green-orange.The high value with an orange-purple color shows high rock density.It can also be used as a reference in determining rock lithology, rock age, and rock density.

Discussion
The gravity anomaly values obtained from GGMplus need to be evaluated with the available field data (Yahaya et al., 2022).The correlation of the gravity anomaly map results with the landslide hazard map aims to compare the potential for landslides from the two maps with the support of local geological maps so that differences and similarities between the two maps can be identified.As a comparison, we use the landslide hazard map we accessed from the INarisk website issued by BNPB.InaRISK is a portal for risk assessment results using ArcGIS as a data service that describes the coverage of disaster hazard areas, affected populations, potential losses, and potential environmental damage (BNPB, 2016)  Rock density describes the bond between the constituent materials and is inversely proportional to the rock's porosity.Porosity is the ratio between the pore volume in a rock mass or the part of the rock volume that is not filled with solid matter.Landslides can occur when two rocks with contrasting density differences meet to become a slip zone.Rocks with low density indicate large porosity values.Rocks with large porosity and steep slopes have a high level of landslide susceptibility, and this is because water can enter through the rock pores, thereby increasing the rock load.If we look at the geological map, the blue color on the residual anomaly map (Fig. 6b) is in the Upper Quaternary Volcanics Formation (Qv (pandermen)) where the age of the rock formation is in the Holocene or the youngest, so the rock density is the lowest of those three formations.Meanwhile, based on the data we collected from the BPBD, the points that most often experience landslides in the Payung area are at the same location, namely in Fig. 5a (red dot), with a relatively high level of damage.
If it is correlated with the Residual Anomaly map in Fig. 6b, the lowest density value (blue color) is in the same area as the highest landslide hazard in Fig. 6a (red color).Trunojoyo Street, marked with a black line, indicates that this road passes through areas with high landslide rates or low density.Areas with a high level of landslide hazard can be caused by a steep slope, so the landslide hazard map (Fig. 6a) is classified as high.The correlation between the gravity anomaly map and the landslide hazard map shows that the Trunojoyo road, which connects the cities of Batu and Kediri, passes through areas prone to landslides, thus endangering road users and potentially cutting off access between areas when a landslide occurs.Based on these conditions, it is necessary to anticipate handling landslides when entering the rainy season to avoid casualties and losses.This research is preliminary study, so it is necessary to conduct further research to determine the subsurface layers of the underground in detail.This research can be carried out using the geoelectrical resistivity method find slip planes, the microtremor method to find seismic vulnerability indices and peak ground acceleration, or the geotechnical method to find maximum soil shear strength.In addition, rainfall study also needs to be identified to calculate the probability of landslides based on rainfall.While the presence of the active Arjuna-Welirang volcano and active faults around the research area can have an impact on slope stability in the research location.However, because the distance is quite far, the impact of seismic activity of Mount Arjuna-Welirang and active faults cannot be concluded as a triggering factor for landslides in the study area.To prove this theory, further research is needed regarding the influence of the seismic activity of the Arjuno-Welirang Volcano and active faults on landslide events at the research location.

Conclusions
The residual anomaly map represents the shallower effect of rock densities, so to analyze the landslide potential based on gravity data, we use the residual anomaly map because landslide potential generally occurs in the surface layer of the soil.Based on the residual anomaly map, the measured values are between -13.4 to -5 mgal (low), -5 to 5 mgal (moderate), and 5 to 12 mgal (high).The low zone shows the area that has the lowest density, which is marked in blue-green; low density values are indicated as cavities or rocks that can hold water and large porosity values so that they have the potential to become slip planes that trigger landslides.The correlation between the gravity anomaly map and the landslide hazard map made by BNPB shows that the Trunojoyo road, which connects the city of Batu with several surrounding areas, is an area prone to landslides.The study written in this article is a preliminary study, so further research is needed to obtain detailed subsurface information.

Fig. 4 .
Fig.4.The map of the complete bouguer anomaly

Fig. 5 .
Fig.5.Modeling of anomaly map gravity result (a).Regional Anomaly (b).Residual Anomaly Fig. 6a is a landslide hazard map prepared by the National Agency for Disaster Countermeasure (BNPB) of the Republic of Indonesia.The figure shows landslide hazard levels consisting of 3 color scales: low, shown in green; medium, in orange; and high, in red.Areas in red indicate a high level of landslide hazard, where the area covers the Trunojoyo road linking Malang district, Batu City, Kediri City, and Jombang Regency.