J. Eng. Technol. Sci., Vol. 52, No. 5, 2020, 707-731

707

Self-Potential Method to Assess Embankment Stability:
A Study related to the Sidoarjo Mud Flow
Sungkono1,*, Masithoh N. Wasilah1, Yekti Widyaningrum2, Wildan M.
Hidayatullah1, Fandi A. Fathoni1 & Alwi Husein3
1

Physics Department, Faculty of Sciences, Institut Teknologi Sepuluh Nopember,
Jalan Arif Rahman Hakim, Surabaya 60111, Indonesia
2
Physics Department, Universitas Bangka Belitung,
Jalan Kampus Peradaban, Bangka Belitung, Indonesia
3
Pusat Pengendalian Lumpur Sidoarjo, Jalan Gayung Kebonsari No. 50,
Surabaya, 60235, Indonesia
*E-mail: hening_1@physics.its.ac.id

Highlights:






Noise-Assisted Multivariate Empirical Mode Decomposition (NA-MEMD) can be
used for qualitative interpretation of self-potential (SP) data.
Continuous Wavelet Transform (CWT) analysis was applied to identify multiple
anomaly sources from SP data.
NA-MEMD and CWT were applied to SP data for qualitative and quantitative
interpretations, respectively.
The combination of NA-MEMD and CWT on SP data was used to identify the
location of fractures and seepage in the LUSI embankment.
The interpretation of the result was supported by direct current resistivity.

Abstract. The stability of an embankment is generally influenced by a number of
factors, such as deformation, fractures, overtopping, seepage, etc. Fractures and
seepage are commonly found in the LUSI (Sidoarjo mud flow) embankment. In
this study, analysis of self-potential (SP) data was applied to identify fractures and
seepage in the LUSI embankment. Noise-Assisted Multivariate Empirical Mode
Decomposition (NA-MEMD) and Continuous Wavelet Transform (CWT) were
applied to determine the location of seepage and fractures in the subsurface based
on SP data. The results were correlated with the 2D direct current resistivity
(DCR) method, which showed that both methods worked well and were
compatible in detecting and localizing fracture and seepage in the LUSI
embankment.
Keywords: fractures; self-potential; signal processing method; seepage; DCR.

1

Introduction

Failure of the LUSI (the Sidoarjo mud flow) embankment is mainly caused by
seepage. Seepage flows through cracks in the embankment body. The degree of
Received August 12th, 2019, Revised May 4th, 2020, Accepted for publication September 19th 2020.
Copyright ©2020 Published by ITB Institute for Research and Community Services, ISSN: 2337-5779,
DOI: 10.5614/j.eng.technol.sci.2020.52.5.8

708

Sungkono, et al.

water seepage in earth-fill embankments can be identified using the self-potential
(SP) method [1-4]. The SP method includes a passive method that measures the
natural static voltage of the ground. This method only needs voltage sensitive
measuring devices and two special electrodes to conduct measurements, making
this method quite fast and simple to execute. The spatial distribution of the
measured SP data can indicate possible anomalous water flow.
Generally, inversion is used to identify the locations (distance and depth) of
sources causing SP anomalies. Forward modeling is required in the inversion
process, where forward modeling of SP in the inversion process generates
idealized bodies (example: sphere, cylinder, inclined sheet, inclined thick sheet)
[5-9]. Several phenomena that disturb embankment stability can be approximated
by simple geometrical models [3]: piping can be represented by a cylindrical
body, internal erosion can be represented by a sphere, differential settlement can
be represented by horizontal fractures, which can be represented by an inclined
sheet. These assumptions cannot be used in general cases, because in some cases
the sources of SP anomalies cannot be represented by idealized bodies but
complex bodies are required [10]. Moreover, embankment instability is caused
by a combination of several idealized bodies [3]. Thus, using a single idealized
body model in the localization of SP sources has high ambiguitiy.
In order to locate instability in the LUSI embankment, the filtered adaptive
method Noise-Assisted Multivariate Empirical Mode Decomposition (NAMEMD) and Continuous Wavelet Transform (CWT) analysis were used. NAMEMD removes noise from the SP data for a qualitative interpretation, while
CWT determines the subsurface position of SP body anomalies for a quantitative
interpretation. In order to identify anomaly sources, the CWT method does not
need assumptions. Although seepage in the LUSI embankment generally occurs
through fractures or rock pores in the embankment body, the CWT analysis of
anomaly sources is assumption-free for determining the seepage positions. As
most geophysical methods, SP suffers from inherent ambiguity. Consequently, in
identifying the positions of SP body anomalies, uncertainty model solutions are
needed [7,11]. CWT analysis can provide the uncertainty positions of SP anomaly
sources.

2

Theory of Self-Potential

The principle of the SP method is to measure the natural potential that is formed
in the subsurface by several mechanisms, namely mineralization (geobattery),
electrochemical (liquid junction or diffusion), electrokinetic (streaming), and
thermoelectric potential. Firstly, SP anomalies associated with the presence of
ore deposits are caused by redox half-reactions [12] involving electron donors
and electron acceptors. The resulting potential is called electrochemical potential

Self-Potential Method to Assess Embankment Stability

709

and generally causes strong SP anomalies (in the order of hundreds of millivolts).
Secondly, SP anomalies related to different mobility of anions and cations in
solutions with different concentrations are LUSI embankment called diffusion
potential [13]. Thirdly, streaming potential related to the motion of ground water
in porous rocks is caused by the pore pressure gradient or hydraulic head [1].
Generally, the streaming potential can have an amplitude from some millivolts to
several hundred of millivolts. The thermoelectric potential occurs due to the
thermal gradient inside the rock, where the thermal gradient increases the energy
of the ions leading to a differential displacement between the ions, which
generates an electrical current [13].
Furthermore, streaming potential with null total current density can be defined
mathematically as follows [1,14,15]:

C

QK
Q KF

 v  v
h
b
f

(1)

where F and C are the formation factor and the streaming potential coupling
coefficient, respectively, while  and h represent the electric potential (V) and
hydraulic head (m), respectively. Furthermore, Qv indicates the excess of
electric charge per unit pore volume [C/m3],  b and  f are the bulk and pore
water electrical conductivity [S/m], respectively, and K is the hydraulic
conductivity [m/s]. The equation has as a consequence that an upstream condition
(the hydraulic head of the measurement electrode is located above the hydraulic
head of the reference electrode) will produce negative SP anomalies on the
surface and vice versa. Consequently, there are negative and positive SP
anomalies for inflow and outflow, respectively.
Pressure is directly proportional to height. Consequently, using Eq. (1) it can be
known that the streaming potential also depends on pressure. In general, fluid
flows from high pressure to low pressure, which means that the location of water
discharge (inflow in embankment seepage) will produce negative SP anomalies,
while the direction of the fluid flow (outflow) has positive SP anomalies on the
surface.
SP data can be disturbed by several types of noise [12], namely: 1) noise
generated by different transient sources; 2) spatial noise associated with strong
heterogeneity of the resistivity distribution in the near surface; 3) other artifacts
associated with the non-polarizing electrodes themselves. Transient noise can be
removed if a fixed dipole is used to record the electrical signals during mapping,
while spatial noise can be filtered if measurements are made with high spatial

710

Sungkono, et al.

density and can be removed using several processing techniques (Fourier-,
wavelet- and empirical mode decomposition-based filters). In addition, other
artifacts can be corrected if the temperature is measured at the place where the
self-potential measurements are carried out.

3

Study Area

Geologically, the LUSI area is part of the Brantas river system with the following
rock sequence (from old to young): Ngimbang, Kujung, Tuban, Ngrayong,
Wonocolo, Ledok, and Lidah formations. At the surface, the formation consists
of alluvial deposits from the Brantas delta (in the north of the Porong river) and
Quaternary volcanic tuff deposits (in the south of the Porong river). The alluvial
delta is mainly derived from the Brantas river, which separates into two rivers,
namely the Kali Mas river in the north and the Porong river in the south, in the
form of alluvium, composed of gravel, sand, and silty clay [16].

Figure 1 Geological map and structures in the LUSI area. Points with numbers
are used to easily find positions in the LUSI embankment.

Self-Potential Method to Assess Embankment Stability

711

The LUSI embankment is the main embankment and is built on the outside of the
affected area (Figure 1). The embankment was built in expansive soil with high
resistivity. Construction was done in view of the LUSI mud eruption and the main
goal was to protect the public and avoid greater losses to the state. The LUSI
embankment was constructed over two faults, the Siring fault and the Watukosek
fault (Figure 1). Nevertheless, the embankment is well constructed and has high
standard penetration test (SPT), shear wave velocity, and resistivity [16-18]. In
order to facilitate the reference to a point or a pond in the LUSI embankment,
each position is denoted by ‘P’ followed by a number (for example: P79 for point
79).
Figure 1 shows the measurement of the SP data lines. Observation of the SP lines
was focused on the north and northeast of the LUSI embankment, because these
areas have highly collapse potential [19]. The embankment in the north and
northeast parts was designed with thicknesses of 11 m and 5 m for the width of
the top part (Figure 2) [20,21], but has been repaired each year (especially before
the rainy season) to handle unstability. Consequently, P79 to P82 of the LUSI
embankment was identified as being about 12 to 15 m thick in 2017 using
integrated vertical electrical sounding and Rayleigh wave dispersion [19].

Figure 2 Dashed lines indicate the localization of the SP and DCR electrodes for
each line at the crest of the dam.

Instability of the LUSI embankment is generally caused by three factors [16,19]:
(1) vertical and horizontal deformations, which will eventually lead to changes
of position of the embankment and cracks in the embankment; (2) the effects of
fluid eruption in large quantities causing overtopping and seepage in the
embankment; (3) the driving force of the LUSI mud flow or fluid on the
embankment, causing embankment failure due to the transformation of the
potential energy of the mudslide into kinetic energy. Moreover, the latter
generally contributes to collapse of the embankment in the rainy season
(especially when it rains a lot).

712

Sungkono, et al.

In this study, the potential difference (PD) method was used, where one electrode
is fixed to a base station and another one is moved at a regular interval.
Measurement of the self-potential signals was carried out with a digital voltmeter
and non-polarizing Cu/CuSO4 electrodes. The spacing between the SP readings
was 5 m and the measured data were corrected [12]. The NA-MEMD and CWT
methods were applied to correct the SP data for qualitative and quantitative
interpretation, respectively.
Furthermore, the direct current resistivity (DCR) method was used to measure the
same line as SP (Figure 2) in order to integrate both data. Measurements of 2D
DCR are generally used to determine the true electrical resistivity of the
subsurface. The DCR method measures the electrical resistivity distribution of
the subsurface using a DC current. The current is transmitted into the ground with
two electrodes (A and B) and the potential difference between a second pair of
electrodes (M and N) is measured. Several electrode arrays can be used to
determine subsurface resistivity, including Wenner, Schlumberger, dipole-dipole
and pole-pole arrays. The Wenner configuration was selected for measuring with
the DCR method.

4

Data Processing

4.1

Self-Potential Data

4.1.1

Noise-Assisted Multivariate Empirical Mode Decomposition
(NA-MEMD)

In this study, Noise-Assisted Multivariate Empirical Mode Decomposition (NAMEMD) [22,23] was used to filter the SP data. This method is able to eliminate
mode mixing in the multivariate empirical mode decomposition (MEMD)
algorithm. Thus, the method is suitable for filtering geophysical data [24]. It was
used to remove linear trends from the data, which were generally correlated with
electrode drift and with base slope and level. Furthermore, the method also
reduced high frequency noise that reflects the telluric effect and was associated
with the very heterogeneous nature (resistivity distribution) of the shallow
subsurface.
The NA-MEMD method operates on SP data and several white Gaussian noise
(WGN) channels with a certain amplitude was added to separate channels from
the SP data. Further, Multivariate Empirical Mode Decomposition (MEMD) [25]
was applied to decompose the SP data and several WGN channels to obtain
several multivariate intrinsic mode functions (IMFs). Several channels of the
multivariate IMFs corresponding to WGN were excluded, while the other

Self-Potential Method to Assess Embankment Stability

713

channels were IMFs of the SP data. Furthermore, some IMFs were selected for
reconstructing free noise of SP data, while a IMFs reflected a linear trend and
IMFs associated with the telluric effect and the very heterogeneous resistivity of
the shallow subsurface were excluded. The filtered SP data were used for
qualitative interpretation, the result of which was compared with the CWT and
DCR results.

4.1.2

Continuous Wavelet Transform (CWT)

Continuous Wavelet Transform (CWT) can be applied to accurately calculate the
position (depth and distance) of sources bodies of potential field anomalies
(gravity, magnetic, and self-potential data) [26,27]. The wavelet transform is a
signal analysis method that uses simultaneous characterization of a signal in the
time and frequency or the space and wave number domain. CWT is obtained by
convolution or cross-correlation between a signal (in this study SP measured
data) f  x  and multiple scales (multiple  ) of wavelets g n with different
translation parameters x ' . The wavelet transform is described in detail in
[28,29]. Mathematically, the wavelet transforms results, Wf ( x ', ) , can be
expressed as follows:


W  g n , f  ( x ', )    n  g n  x  x ' /   f  x  dx

(2)



where g n  x  x ' /   denotes a wavelet that has spatial support centered at x '
and is proportional to  . Using the properties of the wavelets, the scale
parameter (  ) plays the role of the physical dimension depth (Z), while the
translation parameter is the equivalent of position x [30]. The scale parameter is
proportional to the dilation parameter [30]. Furthermore, the scale is described as
the window size of the wavelet. A small scale means a detailed view, which will
increase the resolution in the same time, and vice versa.
In order to increase efficiency, the CWT is processed in the wavenumber u
domain (frequency domain in the spatial function). Thus, this process needs fast
Fourier transforms. Generally, the derivative of order n of the Poisson kernel
family and their Hilbert transforms are used as wavelets, horizontal and vertical,
respectively. The wavelet equation is as follows:
H n  u    2 u  exp  2 u 
n

Vn  u   2 u  2 ui 

n 1

exp  2 u 

where u denotes the wavenumber of the spatial variable x .

(3)

714

Sungkono, et al.

Using the convolution concept, the wavelet transform of a signal using several
scales has the following properties:
1. High frequency (wavenumber) in low scales and vice versa. This has two
consequences, namely: a high frequency signal is clearly visible in the low
scales, and vice versa, and the wavelet transform can act as a filtering process.
2. It contains properties of the signal and wavelet is used in this process.
Based on the properties of the wavelet transforms, the SP and potential data
indicate subsurface anomaly source that are visible in high and low scales. The
wavelet transform property is different from the noise property disturbing the SP
and potential data [27,31].
Furthermore, CWT can identify singularities very well [32]. A singularity is
described by extrema lines in the matrix wavelet transform results (wavelet
transform modulus maxima lines, WTML). For identifying the source of an
anomaly, two or more extrema (minima and maxima) lines are needed, where at
least one line is of maxima and one line of minima [27]. Generally, each body
source of self-potential anomalies will generate a singularity in the CWT results.
Mauri et al. [26] have shown that the dilation value is considered as the space
above the ground surface (z > 0) and as projection below the ground surface (z <
0). The lines of extrema converge at z < 0 in a cone-shaped structure toward the
singularity. The convergence position of the extrema represents the position
(depth and distance) of the SP anomaly source. A more clear description of this
methodology can be found in References [26,27,33].
In order to investigate the accuracy of CWT to predict positions of seepage (fluid
flow) disturbing embankment stability, this study used a combination of real and
imaginary wavelets based on the Poisson kernel function. The real wavelets of
the Poisson family using this process are the second horizontal derivative and the
third vertical derivative, while the imaginary wavelets are the third horizontal
derivative and the second vertical derivative, as discussed in [26].
Variations and errors in depth calculations are caused by three main sources [27].
Firstly, noise in the data (due to a heterogeneous medium) will distort the lines of
the extrema for small values of dilation. Consequently, this will affect the point
of convergence. Secondly, the error due to the heterogeneous medium, which is
calculated based on the quality of the best fit of the cone-shaped structure.
Thirdly, the sampling step has a significant impact on the measured signal and
the signal may characterize the source. The error due to measurement noise can
be further limited by processing the signal with several wavelets (more than one)
in order to statistically constrain the depth.

Self-Potential Method to Assess Embankment Stability

4.2

715

Resistivity Data

The DCR data were analyzed using the L1-norm optimization method. L1-norm’s
error approach is less sensitive to noise levels and produces fewer artifacts [34],
while L1-norm’s modeling commonly gives better imaging results for the
contrast of the model parameter than smoothness constrained least-squares (L2norm) inversion [35,36]. Inversion using L1-norm is done to minimize the
objective function, as follows:
  R d  dobs -dcal   λi R m Wri

(4)

where R d and R m denote the weighted misfit and data roughness, respectively,
J indicates the Jacobian matrix of the partial derivatives, while W describes the
roughness filter. Further, d obs and d cal are the observed and calculated data, and

ri indicates the changed model parameters for the ith iteration.
The Lagrange multiplier for each iteration ( λ i ) is used to minimize model
roughness and to stabilize the inversion process. A simple approach to solving
L1-norm based optimization is the iteratively reweighted least-squares method
(IRLS) [34,35]. In the IRLS algorithm, both R d and R m are iteratively
calculated [36].

5

Result and Discussion

SP and DCR data were acquired in the rainy seasons. The pond condition for each
line is presented in Figure 3. The figure shows that generally the pond contained
water, except at line P76-P77. This condition was used to help in the
interpretation of the processed SP and DCR results.
In order to locate anomalies, the NA-MEMD approach was applied to the SP data.
Figure 4 shows SP data decomposition into 5 IMFs, where the last IMF indicates
the trend of the data. This trend correlates with electrode drift and topography
(Figure 4). Thus, the last IMF must be excluded from the SP data reconstruction.
The first and second IMFs show transient data, which may be caused by strong
heterogeneity of the resistivity variation in the shallow subsurface and may also
have been disturbed by transient sources, such as the telluric effect due to
transient current flow in the ionosphere, lightning, large cumulus clouds, etc.

716

Sungkono, et al.

Figure 3 Pond (point) condition of each line, a) P79-P82; b) P78-P79; c) P76P77; d) P75-P73; e) P75A-P75. Each pond contained fluid, except P76-P77.

Figure 4 SP data decomposed using NA-MEMD into five IMFs.

Self-Potential Method to Assess Embankment Stability

717

Furthermore, the last, the first and the second IMF must be removed.
Consequently, in order to reconstruct the SP anomalies (filtered SP), IMF 3-IMF4
can be summed as in Figure 5. This figure shows that compared with SP
measurement, the SP anomalies were improved by removing the last IMF and
transient data compression by excluding the first and the second IMF. These steps
were performed for all SP measurement lines. Further, the filtered SP data were
interpreted qualitatively (based on the properties of the SP anomalies) and
quantitatively (wavelet multiscale processing).

Figure 5 Comparison of measured SP and filtered SP using the NA-MEMD
approach.

As described in Eq. (1), the qualitative interpretation of self-potential anomalies
of streaming potentials for identifying seepage is easy, where seepage is generally
shown by negative SP anomalies [14,15,37,38]. In addition, the relative intensity
of the discharge is demonstrated by the amplitude of the anomalies. The SP data
filtered using the NA-MEMD algorithm were applied in the qualitative
interpretation for identifying seepage positions.
In the quantitative interpretation of the SP anomalies, the measured SP data were
analyzed using wavelet multiscale analysis to estimate the depth of the source
responsible for the measured SP anomalous field [26,27]. In this study, the second
and third orders of the vertical and horizontal derivatives of the Poisson family
were used to provide the uncertainty positions of the subsurface anomalies. For

718

Sungkono, et al.

this process, the MWTmat software [26] was used to analyze the measured selfpotential field. In the software, the SP data were processed using CWT to identify
the minimum and maximum singularities. Further analysis of both singularities
was done to identify the sources of the SP anomalies. As a quality check of the
results for the source positions, only positions with at least two of the four wavelet
analyses were considered significant. For each wavelet and profile a position was
obtained by estimating the best fit of the cone-line structure [26].
The convergence of the line’s extrema through intersection point z < 0 gives the
depth of the SP source. As shown in Figure 6, CWT analysis was applied to the
SP data measured at P79-P82 (Figure 5) using the second order of the vertical
derivative of the Poisson family.
The figure shows that sources of SP anomalies can be identified by singularity
analysis, where negative and positive singularities (dashed lines) will meet at a
point (diamond) and form a cone. The points correlate to the positions of the SP
anomaly sources. The results were divided into four groups of SP anomaly
sources, where each group was identified through two or more singularities (see
Figure 6). Furthermore, the data were also analyzed using another wavelet, which
is not presented here.

Figure 6 CWT analysis of measured SP to determine anomaly sources.

For each estimated depth, the associated error is based on the quality of the best
fit. The final localization of each source (horizontal and vertical) was based on
the depths and distances of all solutions from the analyzed results and their
associated uncertainty (σ). The analysis results for all of the measured SP data
are presented in Table 1. It shows the CWT analysis result consisting of the depth
and position of each SP source. The WTML in the CWT analysis of the SP data
cannot determine the nature of the object associated with an SP anomaly. In order
to describe the nature of the SP anomaly object, dipolar tomography (advanced
complex CWT analysis of the SP data) and SP data inversion using the finite
element approach can be used [31,39-41].

Self-Potential Method to Assess Embankment Stability

719

Table 1 Source depths (Z) and positions (X) Estimated by CWT of self-potential
data in the LUSI embankment for several lines.
Structure Number of Wavelets
Line P79-P82
S1
4
S2
4
S3
4
S4
4
Line P78-P79
SS1
4
SS2
4
SS3
4
SS4
4
SS5
3
SS6
3
SS7
4
SS8
4
Line P76-P77
WK1
3
WK2
3
WK3
3
WK4
3
WK5
2
WK6
2
Line P75A-P75
KD1
3
KD2
4
KD3
4
KD4
4
KD5
3
Line P75-P73
KT1
3
KT2
4
KT3
4
KT4
3
KT5
3
KT6
2
KT7
2

X

σX

Z

σZ

67.3967
190.4354
138.3866
273.5462

6.5865
9.231651
4.136642
20.5845

-10.2437
-65.1711
-15.6224
-32.2226

9.550729
15.99827
5.483529
9.073569

194.3154
85.26599
32.88858
239.0813
281.9847
352.4498
325.5904
146.5408

0.93842
2.699981
0.72523
0.765973
0
30.30011
5.60797
3.638915

-20.2692
-23.3347
-4.91201
-22.8812
-7.94234
-8.24386
-44.0136
-36.7046

4.851243
8.868472
3.794171
3.88584
7.9434
7.547819
7.640628
19.04271

81.54776
108.4507
197.961
256.6205
191.0173
36.26095

15.03222
3.917961
14.14983
3.27331
11.11475
2.56701

-7.55414
-67.2183
-5.76768
-39.1937
-29.5334
-3.97696

5.48696
21.66722
2.884647
12.449
7.198908
4.493881

222.7747
44.20035
160.4179
290.9196
102.5827

0.345972
0.679158
3.036411
7.335863
1.167678

-8.06587
-10.3666
-31.574
-68.47
-28.991

0.807549
1.859948
12.09849
15.65895
2.241301

88.6605
275.451
506.921
607.0037
592.8443
731.802
420.1123

2.28848
6.781233
1.632278
6.992968
32.07147
18.39118
3.609299

-4.64522
-61.3532
-14.9403
-50.4533
-4.23121
-1.40786
-19.1776

2.188764
25.06461
6.094641
18.28238
4.086568
2.230054
6.172374

Note: Column 2 contains the number of wavelets used in the CWT calculation to locate
the source depths and positions for each line. Depth and position are in meters. σ is the
standard deviation to describe the result uncertainty.

As calibration, the analysis result of SP was compared with 2D imaging from the
DCR method measured over the same lines. The DCR method was analyzed to
minimize the L1-norm model and errors (Eq. (4)) using the IRLS method. The
method is good to minimize tail errors and properly reconstruct the model of the
embankment (the difference in resistivity in the embankment and under the
embankment was assumed to be large).

720

5.1

Sungkono, et al.

Line P79-P82

Figure 7 demonstrates the correlation between SP anomalies, positions sources
and their uncertainty estimated by CWT analysis, and the resistivity of the
subsurface at line P79-P82. At this line, the thickness of the embankment is
estimated around 10 to 15 m, as identified by the first resistivity layer and
supported by Rayleigh wave analysis [18,19]. The low resistivity in the first layer
indicates seepage in the embankment, which is correlated to negative SP
anomalies. Thus, the SP anomalies can be used to determine seepage positions.

Figure 7 (a) Filtered SP anomalies (mV) and (b) 2D resistivity model estimated
by DCR method and source positions of SP anomalies (diamonds) for P79-P82 of
the LUSI embankment, Sidoarjo, Indonesia. Lines indicate a fracture or fault
interpreted by DCR imagery.

The CWT of the SP data at this line identified four groups (Figure 7 and Table 1)
that characterize four sources. Each anomaly controls a different area with
uncertainty less than 21 m for position and 16 m for depth. Furthermore, the data
were inverted using global optimization methods [4,42,43], where the results
suggest that all the anomaly bodies at this line are horizontal cylinders. The
anomalies can be interpreted as piping or seepage [3]. Moreover, seepage sources
were also identified through CWT processing of SP anomalies, namely S1, S3,
and S4.
Furthermore, the negative SP anomalies at 0-70 m, 80-110 m, and 170-250 m
were considered near-surface fractures correlated with S1, S2, and S4,
respectively. Fractures at S1 and S2 were clearly mapped in the 2D DCR

Self-Potential Method to Assess Embankment Stability

721

inversion result, while for S4 the fracture could not be identified from the 2D
DCR analysis, but it was clearly revealed by inversion of the VLF-EM data [19]
and by the GPR method [44]. This may be caused by the measurement point at
around 170-250 m distance being limited (the observed data are based on the
length of the embankment). The fracture network at this line is the result of
deformation. This line has a high rate of vertical deformation as was demonstrated
by the total station data analysis [19].

5.2

Line P78-P79

In order to know the cause of potential instability of the embankment at this line,
qualitative and quantitative interpretations of the SP anomalies were applied. As
described above, the qualitative interpretation is applied to presume seepage
positions from the SP anomalies. During the measurements of the SP and DCR
methods, the embankment was being repaired to increase the height and strength
of the dikes. The thickness of the embankment at this line was around 12 m.
Consequently, the first layer of the embankment, which is made of compacted
earth fill, normally has higher resistivity than the second layer. Nevertheless, if
the first layer is disturbed by fluid (seepage), the embankment will have low
resistivity.
Figure 8 shows negative SP anomalies in several locations, such as at 0-50,
around 100 m and 200 m, and at 175-215 m distance. If we compare the SP
anomalies with 2D resistivity, seepage not always correlates to negative
anomalies, as it can be seen around 250 m distance. In that location the SP
anomalies were not only an effect of seepage but also of fracture (complex
sources). Thus, the seepage position is difficult to identify through qualitative
analysis. Furthermore, CWT analysis should be applied to identify the source of
the SP anomalies.
A CWT analysis was conducted to determine the minimum and maximum
singularities to estimate the sources of the SP anomalies. At this line, the
anomalies indicated structures in 8 different groups. These groups reflect 8
different sources (Figure 8 and Table 1).
The sources of SP anomalies that can be seen in the near surface (in the
embankment) were at SS5 (281.98 m distance) and SS6 (352.45 ± 30.30 m
distance), where pores were the media of seepage. Furthermore, those at SS2
(85.27 ± 2.70 m distance), SS8 (146.54 ± 3.64 m distance), SS1 (194.32 ± 0.94
m distance), and SS4 (239.08 ± 0.77 m distance) can be considered faults. Faults
also disturb embankment stability because water can seep in through the
embankment. Fractures SS2 and SS8 were also revealed in the VLF-EM imaging

722

Sungkono, et al.

[19]. Thus, the seepage at this line may be caused by fractures and pores in the
embankment.

Figure 8 The same description as Figure 7 but for line 78-79.

5.3

Line P76-P77

Figure 1 shows that the measurements of the SP and DCR methods at line P76P77 were taken over the Watukosek fault, where the fault is close to P76. The
thickness of the embankment at this line was around 11-12 m. During the
measurements, the pond at line P76-P77 was not filled with fluid. Consequently,
in the rainy season, the fault or fracture will contain fluid and have low resistivity.
The CWT analysis of the SP data indicated structures in six different groups,
which characterize six sources of SP anomalies (Figure 9 and Table 1). Figure 9
is a comparison between the SP and the DCR analysis results. Negative SP
anomalies at a distance of 75-190 m correlate with low resistivity anomalies in
the near surface (embankment), where the positions contain some structures. The
structures are assumedly fractures or faults breaking through the embankment
material. The structure is considered a water-filled fracture.
In addition, if we compare Figure 1 with Figure 9, the Watukosek fault system is
correlated to WK2 (108.45 ± 3.92 m distance) in the CWT result of the SP
anomalies and has low resistivity around both [45]. Fractures were also found at
WK5 (191.02 ± 11.11 m distance) and WK6 (36.26 ± 2.57 m distance) consistent
with low resistivity. WK5 and WK2 were also identified using VLF-EM imaging
[19].

Self-Potential Method to Assess Embankment Stability

723

Another anomaly, WK1 (81.55 ± 15.03 m distance) is located in a ‘near’ fracture,
which may be an uncertainty effect in the CWT analysis. The last anomaly, WK3
(197.96 ± 14.15 m distance) is considered porous rock containing fluid and has
relatively low resistivity. Furthermore, the low resistivity anomaly from DCR
located around 150 m distance was not identified using CWT, but the NA-MEMD
result shows a small peak.

Figure 9 The same description as Figure 7 but for line 76-77.

5.4

Line P75A-P75

The measurements of the SP and DCR methods were done while the pond at this
line was filled with fluid. Thus, fractures in the embankment would have been
filled with fluid. Consequently, the fracture had lower resistivity than the
embankment. The embankment at this line was estimated around 10 m (at a
different height than P76), where during the measurement the embankment was
under construction.
Figure 10 shows that the amplitude of the negative anomalies of the SP data found
around 200-265 m distance was high. The negative SP zone values indicate
probable paths of seepage in the downstream direction. If we compare this with
DCR resistivity, seepage filling the fracture was shown by the low resistivity
zone. The other seepage position could not be identified through qualitative
analysis. This may be caused by the seepage position being close to other
anomaly sources (for example fractures).

724

Sungkono, et al.

Figure 10 The same description as Figure 7 but for line 75A-75.

Figure 11 Internal erosion on 5 October 2018 is at around 150 to 250 m distance
from P75A [46]. This means internal erosion caused by seepage.

CWT was also applied to the measured SP data to identify the position of each
SP anomaly source. The CWT result indicated five groups, which characterize
five sources of SP anomalies (Table 1 and Figure 10). The first source, KD1
(222.77 ± 0.35 m distance), may correspond to seepage sources, while the second

Self-Potential Method to Assess Embankment Stability

725

anomaly, KD2 (44.20 ± 0.67 m in distance), is a fluid outlet. Both sources were
supported by the 2D resistivity data, with low resistivity in those positions. The
anomalies were probably caused by structures (fractures) as denoted by KD1.
Thus, the water saturated in this position passed through the fracture.
Furthermore, Figure 11 shows slide embankment or internal erosion in this area
on October 5, 2018 [46], which is correlated to negative SP anomalies (200-265
m distance) and may be caused by seepage in the embankment.
Furthermore, KD3 (160.42 ± 3.04 m distance) and KD5 (102.58 ± 1.17 m
distance) as sources of SP anomalies were predicted as fractures based on the
slightly different resistivity around both positions. Both fractures probably are
not affected seepage, because they are relatively deep. Similar to KD3 and KD5,
KD4 (290.92 ± 7.34 m in distance) may be considered a fracture, although the
DCR inversion result did not reach to this position. The interpretation was based
on the geological condition, where the area is controlled by a high rate of
deformation. The deformation is mainly caused by mud loading and collapse of
the overburden due to the removal of mud from the subsurface.

5.5

Line P75-P73

During the SP and DCR measurements, the embankment was being reinforced,
where the embankment part consisted of two parts. Firstly, P75-P74 is relatively
new, with a thickness of around 8-10 m. Secondly, P74-P73 is relatively old, with
a thickness of around 11 m. Consequently, the second part of the embankment is
more compact than the first. This condition is reinforced by the 2D resistivity as
the inversion result of the DCR data (Figure 12). The pond at both parts of the
embankment contained fluid from mud eruption (75% fluid and 25% solid) and
rainwater.
Seepage assessment is indicated by the position of negative anomalies values. At
line P75-P73 there are negative anomalies with potential values between -28 to 40 mV at 0-150 m and 595-795 m distance. The potential value was clearly
measured in the electrical field, possibly caused by seepage flow. Furthermore,
in order to know the depth of the SP anomalies, the CWT method was applied,
which is resulted in five groups (Figure 12 and Table 1), indicates five sources.
The source of the negative anomaly at 0-150 m distance is denoted and interpreted
as KT1 and seepage, respectively. This result was supported by the qualitative
interpretation and low resistivity anomalies. Using the DCR inversion result,
seepage through porous rock in this position was assumed. The deep source of
SP anomalies denoted as KT2, probably indicates a fracture or a fault, where the
source was demonstrated by contrast resistivity in the inversion result of the DCR
data. The fracture assumedly are not affected by seepage in the embankment
because its position is relatively deep. Furthermore, the sources of SP anomalies

726

Sungkono, et al.

denoted as KT3-KT7 are indicated as seepage through the near surface of a
fracture (low resistivity in the near surface).

Figure 12 The same description as Figure 7 but for line 75-73.

6

Conclusions

Embankment stability can be disturbed by fluid. In this paper, the cause of
embankment instability was identified using the SP method. Interpretation of
measured SP data is difficult, because the data represent several different sources,
which can potentially distort subsurface anomalies. NA-MEMD was applied to
improve the SP data by decomposing the data and selecting several IMFs to
represent the SP anomalies. This process is called SP data filtering using NAMEMD. The filtering result is easy to use to identify seepage positions. Seepage
is indicated by negative SP anomalies. This result is supported by low resistivity
in the embankment. Furthermore, CWT was also applied to the SP data to
determine the positions of SP anomaly sources. The result showed that the
anomalies from the SP analysis were similar to those from the 2D resistivity
analysis, where the anomalies were assumed to be faults or fractures and seepage
through fractures and pores. The SP method is easy to implement, simple to
interpret and fast in measurement. In the future, the SP method will be used to
assess the stability of the LUSI embankment.

Acknowledgements
We would like to express our gratitude to the BPLS and PPLS team who helped
with data acquisition in the field. The first author was supported by the Institute

Self-Potential Method to Assess Embankment Stability

727

of Research and Public Services, Institut Teknologi Sepuluh Nopember,
Indonesia, through Project No. 31467/IT2.11/PN.08/2016.

References
[1]

Ikard, S.J., Revil, A., Schmutz, M., Karaoulis, M., Jardani, A. & Mooney
Mooney, M., Characterization of Focused Seepage through an Earthfill
Dam Using Geoelectrical Methods, Groundwater, 52(6), pp. 952-965,
Nov. 2014. DOI: 10.1111/gwat.12151.
[2] Minsley, B., Burton, B., Ikard, S. & Powers, M., Hydrogeophysical
Investigations at Hidden Dam, Raymond, California, USGS Staff –
Published Research, Jan. 2011. http://digitalcommons.unl.edu/
usgsstaffpub/517.
[3] Rozycki, A., Ruiz, Fonticiella, J.M. & Cuadra, A., Detection and
Evaluation of Horizontal Fractures in Earth Dams Using the Self-Potential
Method, Engineering Geology, 82(3), pp. 145-153, Jan. 2006. DOI:
10.1016/j.enggeo.2005.09.013.
[4] Sungkono & Warnana, D.D., Black Hole Algorithm for Determining Model
Parameter in Self-Potential Data, Journal of Applied Geophysics, 148, pp.
189-200, Jan. 2018. DOI: 10.1016/j.jappgeo.2017.11.015.
[5] Biswas, A. & Sharma, S.P., Resolution of Multiple Sheet-Type Structures
in Self-Potential Measurement, J. Earth Syst. Sci., 123(4), pp. 809-825,
Jun. 2014. DOI: 10.1007/s12040-014-0432-1.
[6] Biswas, A. & Sharma, S.P., Optimization of Self-Potential Interpretation
of 2-D Inclined Sheet-Type Structures Based on Very Fast Simulated
Annealing and Analysis of Ambiguity, Journal of Applied Geophysics, 105,
pp. 235-247, Jun. 2014. DOI: 10.1016/j.jappgeo.2014.03.023.
[7] Biswas, A. & Sharma, S.P., Interpretation of Self-Potential Anomaly Over
Idealized Bodies and Analysis of Ambiguity Using Very Fast Simulated
Annealing Global Optimization Technique, Near Surface Geophysics,
13(2), pp. 179-195, 2015. DOI: 10.3997/1873-0604.2015005.
[8] Biswas, A. & Sharma, S.P., Interpretation of Self-Potential Anomaly Over
2-D Inclined Thick Sheet Structures and Analysis of Uncertainty Using
Very Fast Simulated Annealing Global Optimization, Acta Geod Geophys,
52(4), pp. 439-455, Dec. 2017. DOI: 10.1007/s40328-016-0176-2.
[9] Monteiro Santos, F.A., Inversion of Self-potential of Idealized Bodies’
Anomalies Using Particle Swarm Optimization, Computers &
Geosciences, 36(9), pp. 1185-1190, Sep. 2010. DOI: 10.1016/j.cageo.
2010.01.011.
[10] Patella, D., Introduction to Ground Surface Self-Potential Tomography,
Geophysical Prospecting, 45(4), pp. 653-681, Jul. 1997. DOI: 10.1046/
j.1365-2478.1997.430277.x.

728

Sungkono, et al.

[11] Fernández-Martínez, J.L., García-Gonzalo, E. & Naudet, V., Particle
Swarm Optimization Applied to Solving and Appraising the StreamingPotential Inverse Problem, Geophysics, 75(4), pp. WA3-WA15, Jul. 2010.
DOI: 10.1190/1.3460842.
[12] Revil, A. & Jardani, A., The Self-potential Method: Theory and
Applications in Environmental Geosciences, Cambridge: Cambridge
University Press, 2013.
[13] Sharma, P.V., Environmental and Engineering Geophysics, Cambridge:
Cambridge Univ. Press, 1997.
[14] Bolève, A., Vandemeulebrouck, J. & Grangeon, J., Dyke Leakage
Localization and Hydraulic Permeability Estimation through SelfPotential and Hydro-Acoustic Measurements: Self-Potential ‘Abacus’
Diagram for Hydraulic Permeability Estimation and Uncertainty
Computation, Journal of Applied Geophysics, 86, pp. 17-28, Nov. 2012.
DOI: 10.1016/j.jappgeo.2012.07.007.
[15] Grangeon, J.R., Boleve, A., Sanders, J.W. & Glaser, S.D., Self-Potential
Investigation of Moraine Dam Seepage, Journal of Applied Geophysics,
74(4), pp. 277-286, Aug. 2011. DOI: 10.1016/j.jappgeo.2011.06.014.
[16] Sungkono, Husein, A., Prasetyo, H., Bahri, A.S., Monteiro Santos, F.A. &
Santosa, B.J., The VLF-EM Imaging of Potential Collapse on the LUSI
Embankment, Journal of Applied Geophysics, 109, pp. 218-232, 2014.
DOI: 10.1016/j.jappgeo.2014.08.004.
[17] Husein, A., Santosa, B.J. & Bahri, A.S., Seepage Monitoring of an
Embankment Dam Using Resistivity Method: A Case Study of LUSI Mud
Volcano P.79 - P.82 Embankment, Applied Mechanics and Materials, 771,
pp. 213-217, Instrumentation and Measurement Systems, pp. 213-217,
2015.
[18] Laby, D.A., Sungkono, Santosa, B.J. & Bahri, A.S., RR-PSO: Fast and
Robust Algorithm to Invert Rayleigh Waves Dispersion, Contemporary
Engineering Sciences, 9, pp. 735-741, 2016. DOI: 10.12988/ces.
2016.6685.
[19] Sungkono, Feriadi, Y., Husein, A., Prasetyo, H., Charis, M., Irawan, D.,
Rochman, J.P.G.N., Bahri, A.S. & Santosa, B.J., Assessment of Sidoarjo
Mud Flow Embankment Stability Using Very Low Frequency
Electromagnetic Method, Environmental Earth Sciences, 77(196), 2018.
DOI: 10.1007/s12665-018-7333-6.
[20] Team Work of Operational Planning, TTG and BM Monitoring Activity,
Sidoarjo Mud Volcano Agency, Reports of Year-end, 2013. (Text in
Indonesian)
[21] Team work of Operational Planning, TTG and BM Monitoring Activity,
Center of Mud Volcano Control, Reports of Year-end, 2015. (Text in
Indonesian)

Self-Potential Method to Assess Embankment Stability

729

[22] Rehman, N.U., Park, C., Huang, N.E. & Mandic, D.P., EMD Via MEMD:
Multivariate Noise-Aided Computation of Standard EMD, Advances in
Adaptive Data Analysis, 5(2), 1350007, Apr. 2013. DOI: 10.1142/S179
3536913500076.
[23] Rehman, N.U. & Mandic, D.P., Filter Bank Property of Multivariate
Empirical Mode Decomposition, IEEE Transactions on Signal Processing,
59(5), pp. 2421-2426, 2011. DOI: 10.1109/TSP.2011.2106779.
[24] Sungkono, Bahri, A.S., Warnana, D.D., Monteiro Santos, F.A. & Santosa,
B.J., Fast, Simultaneous and Robust VLF-EM Data Denoising and
Reconstruction via Multivariate Empirical Mode Decomposition,
Computers & Geosciences, 67, pp. 125-137, 2014. DOI: 10.1016/
j.cageo.2014.03.007.
[25] Rehman, N.U. & Mandic, D.P., Multivariate Empirical Mode
Decomposition, Proc. R. Soc. A, 466(2117), pp. 1291-1302, 2010. DOI:
10.1098/rspa.2009.0502.
[26] Mauri, G., Williams-Jones, G. & Saracco, G., Mwtmat-Application of
Multiscale Wavelet Tomography on Potential Fields, Computers &
Geosciences, 37(11), pp. 1825-1835, Nov. 2011. DOI: 10.1016/
j.cageo.2011.04.005.
[27] Mauri, G., Williams-Jones, G. & Saracco, G., Depth Determinations of
Shallow Hydrothermal Systems By Self-potential and Multi-scale Wavelet
Tomography, Journal of Volcanology and Geothermal Research, 191(3-4),
pp. 233-244, Apr. 2010. DOI: 10.1016/j.jvolgeores.2010.02.004.
[28] Burrus, C.S., Gopinath, R.A. & Guo, H., Introduction to Wavelets and
Wavelet Transforms: A Primer, 1st edition, Upper Saddle River, N.J.:
Pearson, 1997.
[29] Mallat, S., A Wavelet Tour of Signal Processing, Third Edition: The Sparse
Way, Amsterdam; Boston: Academic Press, 2008.
[30] Moreau, F., Gibert, D., Holschneider, M. & Saracco, G., Wavelet Analysis
of Potential Fields, Inverse Problems, 13(1), p. 165, 1997. DOI: 10.1088/
0266-5611/13/1/013.
[31] Saracco, G., Labazuy, P. & Moreau, F., Localization of Self-Potential
Sources in Volcano-Electric Effect with Complex Continuous Wavelet
Transform and Electrical Tomography Methods for an Active Volcano,
Geophys. Res. Lett., 31(12), p. L12610, Jun. 2004. DOI: 10.1029/
2004GL019554.
[32] Mallat, S. & Hwang, W.L., Singularity Detection and Processing with
Wavelets, IEEE Transactions on Information Theory, 38(2), pp. 617-643,
Mar. 1992. DOI: 10.1109/18.119727.
[33] Saracco, G., Arneodo, A., Beylkin, G. , Fedi, M., Cella, F., Quarta, T. &
Villani, A.V., Special Issue on Continuous Wavelet Transform in Memory
of Jean Morlet, Part II2D Continuous Wavelet Transform of Potential
Fields Due to Extended Source Distributions, Applied and Computational

730

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41]

[42]

[43]

[44]

Sungkono, et al.

Harmonic Analysis, 28(3), pp. 320-337, May 2010. DOI: 10.1016/
j.acha.2010.03.002.
Loke, M.H., Acworth, I. & Dahlin, T., A Comparison of Smooth and
Blocky Inversion Methods in 2D Electrical Imaging Surveys, Exploration
Geophysics, 34, pp. 182-187, Jan. 2003. DOI: 10.1071/EG03182.
Kim, J.H., Supper, R., Tsourlos, P., & Yi, M.J., 4D Inversion of Resistivity
Monitoring Data through Lp Norm Minimizations, Geophysical Journal
International, 195(3), pp. 1640-1656, Dec. 2013. DOI: 10.1093/gji/ggt324
Borsic, A. & Adler, A., A Primal–Dual Interior-Point Framework for
Using the L1 Or L2 Norm on the Data and Regularization Terms of Inverse
Problems, Inverse Problems, 28(9), 095011, Sep. 2012. DOI: 10.1088/
0266-5611/28/9/095011.
Al-Saigh, N.H., Mohammed, Z.S. & Dahham, M.S., Detection of Water
Leakage from Dams by Self-Potential Method, Engineering Geology,
37(2), pp. 115-121, Jun. 1994. DOI: 10.1016/0013-7952(94)90046-9.
Bolève, A., Janod, F., Revil, A., Lafon, A. & Fry, J.J., Localization and
Quantification of Leakages in Dams Using Time-Lapse Self-Potential
Measurements Associated with Salt Tracer Injection, Journal of
Hydrology, 403(3-4), pp. 242-252, Jun. 2011. DOI: 10.1016/j.jhydrol.
2011.04.008.
Jardani, A., Revil, A., Bolève, A. & Dupont, J.P., Three-Dimensional
Inversion of Self-Potential Data Used to Constrain the Pattern of
Groundwater Flow in Geothermal Fields, J. Geophys. Res., 113(B9), pp.
B09204, Sep. 2008. DOI: 10.1029/2007JB005302.
Saracco, G., Mathé, P.É., Moreau, F. & Hermitte, D., Localization and
Characterization of Buried Magnetic Structures Using a Multi-Scale
Tomography. Application to Archaeological Structures on Fox-Amphoux
Site, ArcheoSciences. Revue d’archéométrie, 33(suppl.), pp. 339-343, Oct.
2009. DOI: 10.4000/archeosciences.1806.
Soueid Ahmed, A., Jardani, A., Revil, A. & Dupont, J.P., SP2DINV: A 2D
Forward and Inverse Code for Streaming Potential Problems, Computers
& Geosciences, 59, pp. 9-16, Sep. 2013. DOI: 10.1016/
j.cageo.2013.05.008.
Sungkono, Robust Interpretation of Single and Multiple Self-Potential
Anomalies via Flower Pollination Algorithm, Arabian Journal of
Geosciences, 13(100), 2020. DOI: 10.1007/s12517-020-5079-4.
Abdelazeem, M., Gobashy, M., Khalil, M.H. & Abdrabou, M., A Complete
Model Parameter Optimization from Self-Potential Data Using Whale
Algorithm, Journal of Applied Geophysics, 170, 103825, Nov. 2019. DOI:
10.1016/j.jappgeo.2019.103825.
Husein, A., Sungkono, Wijaya, A. & Hadi, S., Subsurface monitoring of
P.79 - P.82 LUSI Embankment using GPR Method to Locate Subsidence
and Possible Failure, in 15th International Conference on Ground

Self-Potential Method to Assess Embankment Stability

731

Penetrating Radar (GPR), Brussels, Belgium, pp. 268-273, Jul. 2014. DOI:
10.1109/ICGPR.2014.6970427.
[45] Husein, A., Mazzini, A., Lupi, M., Mauri, G., Kemna, A., Santosa, B.J. &
Hadi, S., Investigating the Watukosek Fault System Using Combined
Geophysical Methods Around Lusi Eruption Site, 19, 12266, Apr. 2017.
[46] Team Work of Operational Planning, TTG and BM Monitoring Activity,
Center of Mud Volcano Control, Year-end Report, 2018. (Text in
Indonesia)