http://journal.uir.ac.id/index.php/JGEET 
 

 
 

E-ISSN : 2541-5794  
   P-ISSN :2503-216X 

Journal of Geoscience,  

Engineering, Environment, and Technology 
Vol 5 No 4 2020 

 

 
Arifianto, et al./ JGEET Vol 5 No 4/2020 175 

 

RESEARCH ARTICLE 
 

Analysis of the Surface Subsidence of Porong and Surrounding Area, 

East Java, Indonesia based on Interferometric Satellite Aperture Radar 

(InSAR) Data 

Indra Arifianto1, Rahmat C. Wibowo2* 
1 Department Geological Engineering, Universitas Gadjah Mada, Yogyakarta, Indonesia. 

2 Department Geophysical Engineering, Universitas Lampung, Bandar Lampung, Indonesia. 

 

* Corresponding author: rahmat.caturwibowo@eng.unila.ac.id 

Tel.: +6281 327 507 517 

Received: Jun 17, 2016; Accepted: Sept 21, 2020. 

DOI: 10.25299/jgeet.2020.5.4.5149 

 
Abstract 

Since 2006, the mud volcano erupted in the Porong area due to wellbore failure triggered by an earthquake (2006) epicenter in the Jogjakarta 
area. The mud volcano buried several villages with mud and continued erupted until today. Based on the InSAR data, it can be seen that the 

subsidence is still happening near the dam area and another area that is not related to mud volcano eruption such as the production of two gas fields 

in the Porong area. Moreover, the Porong area is flat and low, less than 4 meters above sea level. The analysis shows that the subsidence rate in 
this area is up to 0.5 m/yr. If this subsidence is continuing, the city can be sinking and flooding during the rainy season. The prediction result from 

this method is about 10 years more and 36 years since in 2006 based on the mudflow rate method. 

 
Keywords: surface subsidence, InSAR, Porong, East Java 
 

 

 

1. Introduction 

The study area is located in the working area of one of the 

oil companies in Indonesia, namely Lapindo Brantas Ltd. The 

Brantas Production Sharing Contract (PSC) has a concession 

area of approximately 7,250 km2, which consists of three oil and 

gas fields in this area: Wunut (WU), Carat (CA), and 

Tanggulangin (TA) (see Fig. 1). In May 2006, the company 

drilled an exploration well in this area (Banjar Panji-1 well) and 

targeted reservoir in Miocene Kujung Limestone (Sawolo et al., 

2010). However, during the drilling, the well penetrated a thick 

formation of under compacted shale of Pliocene Upper 

Kalibeng Formation (from 900 m to 1,871 m). This formation 

is characterized by a low-velocity interval (Figure 1). The 

borehole was cased down to 1,091 m (Pucangan Fm.). After the 

well had penetrated to a total depth of 2,834 m (Kujung Fm), 

this time without a protective casing and they stop the operation 

by putting a cement plug in the steel casing. Unfortunately, the 

Jogjakarta earthquake (2006) trigger and reactivated 

Watukosek fault that crosses the well location, which makes the 

mud is flowing from overpressure formation to surface through 

the existing fault, where the first eruption is about 200 m 

southwest of the Banjar Panji-1 well. 

Since it started to erupt on the 29th of May 2006, the 

remaining lifespan of the Lusi eruption is still questionable. 

However, it is in stable condition recently, and yet, no signs that 

the mudflow will stop to flow. Where initially, the mud eruption 

up to 180,000 m3/day with the subsidence rate is up to 5 cm/day 

in the central vent (Abidin et al., 2008 and Istadi et al., 2009). 

The mud covers area about 7 km2, and the subsidence area may 

cover up to 100 km2 (Davies et al., 2011). In the mid of 2011, 

the mud is remaining to expulse with a discharge rate of 10,000 

m3/day, although it was decreased significantly from previous 

years when mud is flow at 100,000 m3/day, it was expected that 

the mudflow would continue for next 25 to 30 years. This 

eruption has a devastating effect on the Porong District of 

Eastern Java; The flow of mud and ground surface deformation 

has caused significant physical, environmental, economic, and 

social damage to the local and regional communities. Although 

the clay is continuing erupting, the area is also growing as a 

satellite city of Surabaya and Sidoarjo since this area is a very 

strategic location for business and located in the main 

transportation route in the East Java area. 

Since this area is experiencing deformation on the surface 

due to subsurface material eruption, we can use the 

Interferometric Sattelite Aperture Radar (InSAR) methodology 

to study the deformation rate through time in this area. By doing 

time series analysis, we also can predict when the mud volcano 

will stop to erupt since we know that the mud eruption 

decreased exponentially through time and very closely related 

to the subsidence rate as the side effect of the mud eruption. 

This result of this research will also be compared to an 

investigation of studying the probabilities calculation of the 

longevity mud eruption based on the geotechnical aspect by 

Davies et al. (2011). 

2. Geology 

The basement configuration has a NE–SW structural 

orientationand comprises a series of well-defined basement 

ridges withintervening grabens serving as depocenters, which 

contain Tertiary sediments. Clastic sediment deposition and 

carbonate buildup ofthe Ngimbang Formation took place during 

the Eocene and Early Oligocene. The Late Oligocene and 

Miocene sequence is separated from the underlying sequence 

by an unconformity  that served asthe foundation of ENE–

WSW oriented carbonate trends. This platform development, 

which is known as the Kujung limestone,occurred in the late 

Oligocene while the Prupuh and Tuban reefaldevelopment took 

http://journal.uir.ac.id/index.php/JGEET


 

 
176 Arifianto, et al./ JGEET Vol 5 No 4/2020  
 

place in the Early to Middle Miocene. The latestand most 

pronounced period of compressional tectonism began inthe 

Late Miocene time and continued episodically into the 

Pleistocene. The resulting transpressional regime caused the 

left-lateralmovement that runs East-West which resulted in an 

East-Westorientation of the anticline structures (Istadi et al., 

2009). 

Subsequent Pliocene and Pleistocene sedimentation 

consistedof an Eastward-prograding mudstone dominated 

volcaniclastic wedge of Kalibeng and Pucangan Formation, 

with a thickness of 2400–3000m. The volcaniclastic materials 

are derived from theJava volcanic arc south of the Sidoarjo area. 

Themudstone of theKalibeng Formation is over-pressured in 

most parts of the basinwhere the rapid pressure transition 

occurs. The high sedimentationrate, rapid deposition and burial 

of thick shale in the depocenters that occurred during Pliocene 

and Pleistocene haveresulted inzones of under-compacted 

shales (Willumsen and Schiller, 1994;Schiller et al., 1994). This 

highly plastic zone is over-pressured dueto the presence of 

excessive trapped water and maturing organic-rich materials. 

The many mud volcanoes found in East and CentralJava 

correlate with this thick and rapidly deposited clay-

bearingKalibeng Formation sediments. 

The active tectonics coupled with over-pressured 

maturingorganic-rich sediments makes East Java an ideal 

setting for mudvolcanism. The existence of the NE–SW 

trending Watukosek faultsystem in the LUSI area provides the 

necessary conduit for mud,fluid and gases extruding from a 

deep horizon to flow up to thesurface to form the LUSI mud 

volcano (Istadi et al., 2009). 

 

Fig.1. The study area is 10 x 10 km located in the south of Sidoarjo city (left) and the NE-SW seismic line show subsurface formation and the 

mud volcano source (right) (Istadi et al., 2009) 

3. Methods 

3.1 Data Source 

In this paper, we used 27 SAR images from ascending data, 

acquired by the Sentinel-1 satellite between January 2019 and 

April 2020. The spatial and temporal baseline ranged from 15 

to 250 m and 12 days to 48 days, respectively.  

Gamma software was applied to perform the 

interferometric process, and included the following steps: 

master image selection, topographicphase removal, and 

interferogram generation. In the study, we employed the 

SBAStechnology to acquire land subsidence information. The 

topographic phase was removed by using the 90 m resolution 

Shuttle Radar Topography Mission (SRTM) DEM data from 

the United States Geological Survey (USGS). The spatial 

adaptive filter and temporal filter were applied to separate the 

atmospheric phase, noise phase, and nonlinear deformation 

phase. Then, two-dimensional spatial high-pass filtering was 

used to remove the orbital linear ramp and other long-

wavelength noise from the differential interferograms. A total 

of 27 interferogram pairs were built with spatial and temporal 

baselines smaller than 250 m and 24 days, respectively (Fig. 2). 

3.2 SB Technique 

In SB methods (see, e.g. Berardino et al., 2002; Schmidt et 

al., 2003), a network of inter ferogram pairs is created with 

small temporal and geometrical baselines to limit decorrelation 

noise. Noise is further reduced by applying range and azimuth 

filters (Just et al., 1994) and spatial multi-looking. However, 

filtering reduces spatial resolution, which can fail to detect 

some stable isolated pixels. We use the full-resolution SB 

method of Hooper (2008), in which pixels are identified among 

the candidate pixels, based on phase analysis. For SB 

processing, we form interferograms with an estimated mean 

coherence above a threshold of 0.5. We also add additional 

connections to ensure that there are no isolated clusters of 

images (see Fig.2), creating a total of 27 interferograms. 

 
Fig.2. A network of interferogram pairs obtained from Sentinel-1 used 

in small baseline subset(SBAS) InSAR, circles represent images and 
lines represent the interferograms formed. InSAR Network Build set 

the spatial 250m and temporal baselines 24 days 

4. Results and Discussion 

4.1 Results 

The first image processing result is the interferogram SAR 

image. This image clearly shows the effect of surface condition 

on the coherence and decorrelation. From satellite image, city 

or village have good coherence in the InSAR phase image, 

while vegetation or plantation area show scattered temporal 

decorrelation. This scattered could happen due to the 

randomness of the phase by the signal reflection of the trees. 

The interferometric SAR image generates the deformation 

phase, noise phase, and the atmospheric signal. Later on, this 

Atmospheric signal will be reduced and neglected for small 

scale area or deformation. Furthermore, some decorrelation and 

speckle noise can be removed by applying some filtering 

parameters. Since we have several temporal decorrelations in 

our area, we can use adaptive filtering to reduce the noise on 



 

 
Arifianto, et al./ JGEET Vol 5 No 4/2020 177 

 

our data without changing the pixel size. To do the adaptive 

filter is to adapt the local slope before averaging the value. In 

Gamma, the adaptive program filter is used for complex-valued 

images, such as interferograms. This operation helps to reduce 

decorrelation noise and improves the appearance of the InSAR 

phase and helps in the unwrapping phase. 

For the incoherent area, we can apply masking in the 

software to remove lousy phase data. In this study area, bad 

coherence area exists inside the LUSI dam, where the study 

target area. Therefore, we don't use the masking image to 

analyze displacement near the vent of the mud volcano. Image 

result from the processing that we use for analysis is mainly an 

unwrapping image because we will investigate the 

displacement in a metric unit in LOS (line of sight) view. In this 

study, we use Geocode, LOS, filtered unwrapping image to do 

time series analysis.  

First of all, we need to export the image data to .tiff in order 

can be analyzed in Matlab or other software. In this case, we 

analyze the image of the GIS software. The GIS software allows 

us to load the other geographic data from online, make a 

section, and easier to add a feature such as the border of the 

study area, dam location, georeferenced image from 

publications. Moreover, the GIS software has a simple 

mathematical calculator that can be used for image stacking and 

time series analysis. 

Before calculating the time series analysis, we need to do 

calibration of the input of tiff data (Geo_Disp_LOS_m_filt.tiff). 

As we know, there is several InSAR image that has noise from 

the atmospheric signal, which make the undeformed area in the 

study area as if experiencing deformation but in the same rate 

(from 0.01m-0.04 m) as we can see in the following Fig.3. After 

checking all data, we pull up or pull down the undeformed 

reference point to zero levels of deformation. 

The next step is time series analysis. The time series 

analysis itself helps us to track deformations trough time in our 

study area. Furthermore, we can make a prediction of the mud 

eruption in our study area based on the exponential of 

subsidence. The processed atmospheric signal data above is 

then processed by computing or by applying time-series linear 

equation d = Gm. 

𝑑 = 𝐺𝑚 ->[

0
𝜑1(𝑃)
𝜑2(𝑃)

] =  [
1 0 0
1 −1 0
0 1 −1

] [

∅1(𝑃)

∅2(𝑃)
∅3(𝑃)

]     (1) 

 

Where d is data vector of InSAR unwrapped phase data 

(Geo_Disp_LOS_m_filt), m is the deformation of each SAR 

image in a particular time (mt0, mt1,mt2, mt3,….mtn), and G is 

the interferogram baseline network model (in which every 

subset is connected to). We assumed that mt0 is zero 

deformation. In this project, we doing time series in thirteen 

data set with the temporal baseline is 24 days, where we only 

choose a data set that has a continued timeline (Fig.4). Thus, it 

will be easier to create time series because we only need to add 

the InSAR deformation data of m(t0)+d(t1) to get mt1, add 

m(t1)+d(t2) to get mt2, etc. see the Time series in the (Fig. 4). In 

this study area, the identified deformation value is negative 

corresponding to subsidence since the direction azimuth is 

descending from north to south, and the line of sight (range) is 

to the east. While calculating the subsidence rate, we just divide 

the subsidence (in meter) with the period (in day or year) of the 

two-deformation time or from the slop of the plot between 

subsidence and time. 

 

Fig.3. Interferometric SAR image after applying adaptive filter 0.5 (study area in the black box) and Reducing the effect of the atmosphere in 
Geo_Disp_LOS_m_filt.tiff (red dot is the reference point) 

 

Fig.4. Selected Interferogram 

4.2 Discussion As mentioned in the introduction that this area is part of the 

concession area of one oil company (Lapindo Brantas Ltd). 



 

 
178 Arifianto, et al./ JGEET Vol 5 No 4/2020  
 

Lapindo Brantas has four oil fields in this block area, namely; 

Tanggulangin, Wunut, Carat, and Watu Dakon oil field. Two of 

them are located near the LUSI eruption center. They are 

Tanggulangin in the East of the study area and the Wunut field 

in the West of the study area (see Fig. 1). Based on the time 

series InSAR image, we can see several deformations in the 

area. The deformation in this area corresponds to the LUSI 

eruption center, and both of the gas field (Fig.5). As we can see 

the significant deformation is still happening in the Northeast 

and the West of the levee area until now, where the deformation 

also happened in the past in the North East levee wall and the 

west region in the railroad (Istadi et al. 2012). Therefore, the 

government or operator should pay attention to this area highly 

deformed area to prevent dam failure. 

From the time series map, we choose and plot three points 

reference in time series graph to know the subsidence rate in 

those three depocenters. The first point is corresponding to the 

LUSI, the second point, which is located in the East of the study 

area is corresponding to the Tanggulangin gas field, and the 

third point is located in the West of study area corresponding to 

Wunut gas field (Fig.6). This subsidence rate also becomes 

necessary in this area because this area is a flat and low 

elevation area; therefore, a significant subsidence rate within a 

year can lead to sinking and flooding when entering the rainy 

season. 

 

Fig.5. Average subsidence velocity in the study area, in the middle of the study area corresponding to Lusi, in the Northeast area is the 
Tanggulangin gas field and in the West area is the Wunut gas field.

 

Fig.6. Time series analysis with temporal baseline maximum of 24 days

From the graph, we can see that there is a change in the 

subsidence rate or displacement velocity of each point (Fig.7). 

The subsidence rate change in the gas field must be caused by 

the increase of gas production in the gas field. As reported by 

the online media, that Lapindo Brantas will increase twice their 

production by the end of this year. That is mentioned in the 

news that the production of the gas field in 2018 to 2019 is 

around 1.1 million m3 gas/day, and they will increase the 

production twice to 2.2 million m3 gas/day total from two gas 

fields. This gas production account for the significant 

subsidence of the Tanggulangin area and the Wunut area from 

0.22 m/yr to 0.6 m/yr and from 0.15 m/yr to 0.3 m/yr, 

respectively. From the graph, we also can see that the 

production increase in the Tanggulangin area starts around 

early October, while in Wunut area begins at the end of 

December 2019 (Fig. 7). 



 

 
Arifianto, et al./ JGEET Vol 5 No 4/2020 179 

 

Inside the LUSI area, the subsidence rate is decreased 

nonlinearly through time, which is agree with Davies et al., 

(2011) flow rate model. Davies et al. (2011) suggest that the 

flowrate of the LUSI will decrease to less than 100 liters/day in 

a minimum of 26 days (Fig.8). The other estimation research is 

by Istadi et al. (2009) predict that the source of mud will be 

depleted in 23-35 days. However, according to Davies, this 

analysis is flawed because the eruption will continue after the 

mud source depleted, and they use a constant eruption rate 

rather than reducing through time. 

Then in this study, we plot the subsidence rate from our 

research with the initial subsidence rate acquired from Abidin 

et al. (2008). Abidin analyzes the data taken from a GPS 

database that has an initial subsidence rate is ranging from 14 

meters/year to 6.5 meters/year (exponentially decrease) near the 

vent center of mud volcano from May 2006 to early 2017 

(Figure 8). This research is also supported by Deguchi et al. 

(2007), which show the same result from InSAR data, where 

from October to November 2006, the maximum subsidence rate 

was 12.4 meters/year and from November 2006 to January 2007 

was 7.3 meters/year (Fig.8). However there also a publication 

from Shirzae et al. (2015) suggest that average deformation 

velocity from InSAR data during early eruption phase (2006-

2007) is only 80 cm/year while during while the average 

deformation velocity in the Lusi and Wunut during 2006 to 

early 2011 is about 5 cm/year to 8cm/year, respectively. This 

research from Shirzae et al. (2015) and one of the results from 

Deguchi et al. (2007) shows meager subsidence rate in the early 

time of the eruption. This mistake is because they use a longer 

temporal baseline that will result in the phase unwrapping 

problem due to a very high deformation rate in a small area. 

Finally, from the plot of the subsidence rate in the 

logarithmic scale from a few works of literature mentioned 

above, we can predict when is the subsidence rate decrease, 

which is linear with the volume of mud erupted, or flow rate of 

eruption is not significant anymore. In this case, we determine 

the subsidence rate is less than 1 cm per year is achieved in 2036 

or after 30 years since the first eruption day. The result of this 

project is still linear with the other method of prediction of the 

longevity of mud eruption.

 

Fig.7. Time series of surface subsidence in three reference point (TA = Tanggulangin; WU = Wunut) 

 

Fig.8. Mud eruption longevity prediction from Davies et al., 2011 (left) and from this research (right) 
 

5. Conclusion 

In this study, we were successfully processing and applied 

InSAR to detect and calculate deformation (Time Series) within 

the study area in Porong, Sidoarjo. During 2019-2020, the 

subsidence rate in the Porong area due to the mud volcanoes is 

exponentially decreasing in one year from 0.9 m/yr – 0.3 m/yr. 

Due to Field production, the Porong area is experiencing more 

than twice the increase of subsidence rate from 0.22 m/yr to 0.6 

m/yr in the Tanggulangin area and twice subsidence rate in the 

Wunut area (0.15 m/yr – 0.3 m/yr). This study and the previous 

study from Davies et al. (2011) suggest that LUSI will stop to 

erupt the mud in 2036 (after 30 years) from the subsidence rate 

and mud discharge rate, respectively. 

Acknowledgements 

The research was supported by the ErSE 3900 InSAR 

Course project. The authors thank the King Abdullah 

Universityof Science and Technology for facilities support. 

References 

Abidin H.Z., Kusuma M.A., Andreas H., Gamal M., and 

Sumintadireja P., 2008. GPS-Based Monitoring of 

Surface Displacements in the Mud Volcano Area, 

Sidoarjo, East Java. Observing our Changing Earth: 



 

 
180 Arifianto, et al./ JGEET Vol 5 No 4/2020  
 

Proceedings of the 2007 IAG General Assembly, 

Perugia, Italy, July 2 - 13, 2007. 

Berardino, P., G. Fornaro, R. Lanari, and E. Sansosti., 2002. “A 

New Algorithm for SurfaceDeformation Monitoring 

Based on Small Baseline DifferentialSAR 

Interferograms.” IEEETransactions on Geoscience and 

Remote Sensing 40 (11): 2375–2383. 

Davies., R.J., Mathias S.A., Swarbrick R.E., and Tingay M.J., 

2011. Probabilistic longevity estimates for the LUSI mud 

volcano, East Java. Journal of the Geological Society, 

London, Vol. 168, 2011, pp. 517–523. doi: 

10.1144/0016-76492010-129. 

Deguchi, T., Y. Maruyama, and C. Kobayashi., 2007. 

Monitoring of Deformation Caused by the Development 

of Oil and Gas Field Using PALSAR and ASTER Data, 

Powerpoint presentation file, 14 January, ERSDAC 

(Earth Remote Sensing Data Analysis Center), Japan. 

Istadi B.P., Pramono G.H., Sumintadireja P., and Alam S., 

2009. Modeling study of growth and potential geohazard 

for LUSI mud volcano: East Java, Indonesia. Marine and 

Petroleum Geology 26 (2009) 1724-1739. 

Istadi B.P., Wibowo H.T., Sunardi E., Hadi S., and Sawolo N., 

2012. Mud Volcano and Its Evolution. Earth Science 

book chapter, Intech Open, DOI:10.5772/24944 

Just, D., and R. Bamler. 1994. “Phase Statistics of 

Interferograms with Applications to SyntheticAperture 

Radar.” Applied Optics 33 (20): 4361–4368. 

Sawolo N., Sutriono E., Istadi B.P., Darmoyo A.B., 2010. Was 

LUSI caused by drilling? – Authors reply to discussion. 

Marine and Petroleum Geology 27 (2010), 1658-1575. 

Schiller, D.M., Seubert, B.W., Musliki, S., Abdullah, M., 1994. 

The reservoirpotential of globigerinid sands in 

Indonesia. In: Proceedings of the 23rdAnnual 

Convention, Proceedings of the Indonesian Petroleum 

Association,vol. 1, p. 189ff. 

Schmidt, D. A., and R. Burgmann. 2003. “Time-Dependent 

Land Uplift and Subsidence in theSanta Clara Valley, 

California, from a Large Interferometric Synthetic 

Aperture Radar Data Set.”Journal of Geophysical 

Research 108 (B9): 2416–2428. 

Shirzae M., Rudolph M.L., and Manga M., 2015. Deep and 

shallow sources for the Lusi mud eruption revealed by 

surface deformation. Geophys. Res. Lett., 42, 5274–

5281, doi:10.1002/2015GL064576. 

Willumsen, P. and Schiller, D.M., 1994. High quality 

volcaniclastic sandstone reservoirsin East Java, 

Indonesia. In: 23rd Annual Convention, vol. I. IPA, pp. 

101–111. 

 
© 2020 Journal of Geoscience, Engineering, 

Environment and Technology. All rights reserved. This 

is an open access article distributed under the terms of 

the CC BY-SA License (http://creativecommons.org/licenses/by-sa/4.0/). 

 

http://creativecommons.org/licenses/by-sa/4.0/
http://creativecommons.org/licenses/by-sa/4.0/