Geophys. J. Int.
(2022)
230,
2089–2097
https://doi.org/10.1093/gji/ggac162
GJI Seismology
Supershear shock front contribution to the tsunami from the 2018
M
w
7.5 Palu, Indonesia earthquake
Faisal Amlani
,
1
,
*
Harsha S. Bhat
,
2
Wim J.F. Simons,
3
Alexandre Schubnel,
2
Christophe Vigny,
2
Ares J. Rosakis,
4
Joni Efendi,
5
Ahmed E. Elbanna,
6
Pierpaolo Dubernet
2
and Hasanuddin Z. Abidin
5
,
7
1
Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, CA 90089, USA
2
Laboratoire de G
́
eologie,
́
Ecole Normale Sup
́
erieure, CNRS-UMR
8538
, PSL Research University, 75231 Paris, France. E-mail:
harsha.bhat@ens.fr
3
Faculty of Aerospace Engineering, Delft University of Technology, 2629 HS Delft, Netherlands
4
Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA
5
BIG (Badan Informasi Geospasial / Geospatial Information Agency of Indonesia), Java,16911, Indonesia
6
Department of Civil and Environmental Engineering, University of Illinois at Urbana Champaign, Urbana, IL 61801, USA
7
Department of Geodesy and Geomatics Engineering, Institute of Technology Bandung, Kota Bandung, Jawa Barat 40132, Indonesia
Accepted 2022 April 15. Received 2022 January 10; in original form 2022 March 1
SUMMARY
Hazardous tsunamis are known to be generated predominantly at subduction zones. However,
the 2018
M
w
7.5 Palu (Indonesia) earthquake on a strike-slip fault generated a tsunami that
devastated the city of Palu. The mechanism by which this tsunami originated from such an
earthquake is being debated. Here we present near-field ground motion (GPS) data confirming
that the earthquake attained supershear speed, i.e. a rupture speed greater than the shear wave
speed of the host medium. We subsequently study the effect of this supershear rupture on
tsunami generation by coupling the ground motion to a 1-D non-linear shallow-water wave
model accounting for both time-dependent bathymetric displacement and velocity. With the
local bathymetric profile of Palu bay around a tidal station, our simulations reproduce the
tsunami arrival and motions observed by CCTV cameras. We conclude that Mach (shock)
fronts, generated by the supershear speed, interacted with the bathymetry and contributed to
the tsunami.
Key words:
Earthquake dynamics; Geodetic instrumentation; Tsunamis; Numerical mod-
elling; Numerical approximations and analysis.
1 INTRODUCTION
Tsunamis are well-known to be amongst the most destructive con-
sequences of earthquakes (Synolakis & Okal
2005
;Bryant
2008
;
Pugh & Woodworth
2014
;R
̈
obke & V
̈
ott
2017
), and the 2018 Palu
earthquake was no exception: it generated a devastating tsunami
(Fritz
et al.
2018
;Mai
2019
) in the nearby Palu bay in which hun-
dreds were killed and tens of thousands more displaced from their
homes (ASEAN
2018
). However, this was a very unexpected occur-
rence since the earthquake was associated with the predominantly
in-plane ground motion produced by strike-slip ruptures.
∗
Now at: Universit
́
e Paris-Saclay, CentraleSup
́
elec, ENS Paris-Saclay, CNRS,
LMPS - Laboratoire de M
́
ecanique Paris-Saclay, Gif-sur-Yvette, France.
These motions are not known to excite significant waves and
hence the underlying physical mechanisms behind the tsunami have
largely remained a mystery (Umar
et al.
2019
).
Studies conducted to explain the phenomenon have not arrived
at definitive conclusions (Muhari
et al.
2018
) nor have adequately
captured observed records (Jamelot
et al.
2019
; Heidarzadeh
et al.
2019
; Ulrich
et al.
2019
); the main consensus appears to be that
some form of ground motion [e.g. landslides (Sassa & Takagawa
2019
) or the reverse-slip motion of the fault (He
et al.
2019
)],
amplified by the bay, is to blame.
However, a key notable feature of this earthquake is that it rup-
tured at supershear speed (Bao
et al.
2019
; Socquet
et al.
2019
),
which results in a manifestation of shear and Rayleigh Mach fronts
carrying substantial vertical velocity with relatively slow attenu-
ation over large distances (Bernard & Baumont
2005
; Dunham
&Bhat
2008
). The existence of supershear earthquakes has been
C
The Author(s) 2022. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access
article distributed under the terms of the Creative Commons Attribution License (
https://creativecommons.org/licenses/by/4.0/
), which
permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
2089
Downloaded from https://academic.oup.com/gji/article/230/3/2089/6604081 by California Institute of Technology user on 23 June 2022
2090
F. Amlani et al.
proven theoretically and experimentally since the early 1970s (Wu
et al.
1972
; Burridge
1973
; Andrews
1976
;Das&Aki
1977
;
Rosakis
et al.
1999
;Xia
et al.
2004
;Passel
`
egue
et al.
2013
). The
1979
M
w
6.5 Imperial Valley (California) earthquake was the first
naturally observed supershear earthquake rupture (Archuleta
1984
).
Since then, several more (although rare) earthquakes have been
recorded to propagate at supershear speeds: the
M
w
7.4 1999 Izmit
in Turkey (Bouchon
et al.
2001
), the
M
w
7.8 2001 Kunlun (Robin-
son
et al.
2006
) and the
M
w
7.8 2002 Denali (Ellsworth
et al.
2004
;
Mello
et al.
2014
), to name a few.
Although the overall tsunami behaviour at Palu is likely a combi-
nation of several effects that include these supershear dynamics as
well as landslides, recent studies (Ulrich
et al.
2019
;Jamelot
et al.
2019
;Oral
et al.
2020
) suggest that the influence from phenomena
such as the latter may be secondary: the rupture itself may have
created adequate seafloor movement to excite the tsunami, which
was subsequently amplified by the shallow and narrow 2-D/3-D ge-
ometric features of the Palu bay. Indeed, high-frequency waveform
observations (1 Hz) from carefully calibrated analysis of CCTV and
social media camera footage near the Pantoloan (PANT) station sug-
gest a near instantaneous, high-frequency, tsunami arrival (Carvajal
et al.
2019
)—consistent with a coseismic source near the coast. This
arrival is not captured by observations that were made by the one
working acoustic sensor at the PANT tidal gauge (Sep
́
ulveda
et al.
2020
), whose resolution (0.02 Hz, or 1 measurement per minute)
is too coarse to have captured the much shorter wavelength [a 1–2
min period (Carvajal
et al.
2019
)] of the tsunami as observed by the
high-resolution camera analysis.
Hence the primary objective of this work is to explain the near
instantaneous arrival of the tsunami by elucidating the tsunami gen-
eration process of the supershear strike-slip Palu earthquake in order
to more fully understand the role played by the corresponding rup-
ture dynamics on the observed timing and first motions of the sub-
sequent tsunami. In particular, we incorporate a feature neglected
in previous modelling studies on Palu (Ulrich
et al.
2019
;Jamelot
et al.
2019
) that is a defining characteristic of supershear earth-
quakes: the
velocity
of the ground motion (Bernard & Baumont
2005
; Dunham & Bhat
2008
). Using a model validated by the first
near-field evidence (also presented in this paper) of supershear at
Palu, our results imply that ground velocities, which better represent
the intricacies of the Mach fronts, may further explain the observed
motions of the tsunami. Since other studies (including those investi-
gating landslides and liquefaction) have adequately captured much
of the observed run-up amplitudes and some local inundations, the
scope of this paper is to focus on the arrival, first motions and
phases inferred from CCTV camera records near the PANT station
(Carvajal
et al.
2019
;Sep
́
ulveda
et al.
2020
).
This manuscript is organized as follows. Section 2 describes
the overall methods and data used in this study, including earth-
quake displacements/velocities simulated by a supershear rupture
model (Section 2.1); a corresponding tsunami model (Section 2.2)
that accounts for such dynamic displacements/velocities (numeri-
cally simulated via a novel pseudo-spectral methodology for solving
the shallow water wave equations); and GPS ground displacement
data recorded at the PALP station during the Palu earthquake (Sec-
tion 2.3). The results and discussion of Section 3 provide the afore-
mentioned evidence of supershear observed directly from those GPS
records (Section 3.1), where the corresponding rupture dynamics
are then numerically modelled and subsequently incorporated into
the tsunami equations for comparison with observed waveforms
acquired from the PANT observations (Section 3.2). Concluding
remarks are provided in Section 4.
2 METHODS AND DATA
2.1 Supershear modelling
For the considered supershear earthquake dynamics and the cor-
responding rupture modelling (mutually validated by GPS data
in Section 3.1 and subsequently used to source the Palu tsunami
configuration in Section 3.2), we use existing numerical simula-
tions conducted by Dunham and Bhat (Dunham & Bhat
2008
).
Such simulations have been produced by a staggered-grid finite-
difference (FD) code (Favreau
et al.
2002
) with the fault boundary
conditions implemented using a staggered-grid split-node (SGSN)
method (Dalguer & Day
2007
). Since Dunham & Bhat (
2008
)have
provided non-dimensionalized solutions, we simply dimensionalize
their results for the Palu earthquake by using a shear modulus of 30
GPa, stress drop of 20 MPa and a shear wave speed of 3.5 km s
–1
.
The depth of the rupture is assumed to be 7.5 km. These parameters,
reasonable for crustal earthquakes, have been chosen to best fit the
observations. The resulting particle velocities and displacements
are presented in Section 3.
2.2 Tsunami modelling
2.2.1 Governing shallow water wave equations with dynamic
ground displacement and velocity
Using the synthetic particle motions generated by the 3-D supers-
hear earthquake model described above (which, as later discussed
in Section 3.1, agree with PALP GPS records and are reasonably
assumed to sweep past the bay near Pantoloan), a 1-D non-linear
shallow water wave model incorporating time-dependent ground
movements of velocity and displacement (Dutykh & Clamond
2016
)
is utilized to simulate the generation and propagation of the tsunami.
Such a model employs the depth-averaged shallow water approxi-
mation of the Euler equations, which can be written as a system of
coupled hyperbolic partial differential equations given by
⎧
⎪
⎪
⎨
⎪
⎪
⎩
∂
H
∂
t
+
∂
(
Hu
)
∂
y
=
0
,
∂
(
Hu
)
∂
t
+
∂
(
Hu
2
)
∂
y
+
gH
∂η
∂
y
=
0
,
0
≤
y
≤
L
,
t
≥
0
.
(1)
Here,
u
(
y
,
t
) is the fluid velocity,
η
(
y
,
t
) is the sea surface height
and
H
(
y
,
t
)
=
η
(
y
,
t
)
+
h
0
(
y
)
−
h
(
y
,
t
) is the absolute height from
the bed-level to the water surface for an initial at-rest bathymetry
tpo
h
0
(
y
). The constant
g
is the acceleration due to gravity. The
entire domain of length
L
is subjected to a time-dependent ground
perturbation
h
(
y
,
t
) which—together with the corresponding ground
velocity
∂
h
(
y
,
t
)/
∂
t
included in eq. (
1
)—sources the subsequent
tsunami dynamics. In the specific Palu bay configuration considered
in this work (Section 3.2), these values are determined from the 3-D
supershear earthquake model as discussed in Section 3.1.
2.2.2 Pseudospectral numerical analysis based on Fourier
continuation
The complete non-linear system given by (
1
) is solved using a
numerical scheme based on an accelerated Fourier continuation
(FC) methodology for accurate Fourier expansions of non-periodic
functions (Lyon & Bruno
2010
; Albin & Bruno
2011
;Amlani&
Bruno
2016
). Considering an equispaced Cartesian spatial grid on,
for example, the unit interval [0,1] (given by the discrete points
y
i
=
i
/(
N
−
1),
i
=
0,...,
N
−
1), Fourier continuation algorithms append
Downloaded from https://academic.oup.com/gji/article/230/3/2089/6604081 by California Institute of Technology user on 23 June 2022
Supershear contribution to the 2018 Palu tsunami
2091
a small number of points to the discretized function values
η
(
y
i
),
u
(
y
i
) in order to form (1
+
d
)-periodic trigonometric polynomials
η
cont
(
y
),
u
cont
(
y
) that are of the form
η
cont
(
y
)
=
M
k
=−
M
a
k
e
2
π
iky
1
+
d
,
u
cont
(
y
)
=
M
k
=−
M
b
k
e
2
π
iky
1
+
d
(2)
and that match the given discrete values of
η
(
y
i
),
u
(
y
i
), i.e.
η
cont
(
y
i
)
=
η
(
y
i
),
u
cont
(
y
i
)
=
u
(
y
i
)for
i
=
0,...,
N
−
1. Spatial derivatives for
the shallow water system given by (
1
) are then computed by exact
term-wise differentiation of (
2
)as
∂η
∂
y
(
y
i
)
=
∂η
cont
∂
y
(
y
i
)
=
M
k
=−
M
2
π
ik
1
+
d
a
k
e
2
π
iky
i
1
+
d
,
∂
u
∂
y
(
y
i
)
=
∂
u
cont
∂
y
(
y
i
)
=
M
k
=−
M
2
π
ik
1
+
d
b
k
e
2
π
iky
i
L
+
d
.
(3)
In essence, FC algorithms add a (fixed) handful of additional
values to the original discretized function in order to form a pe-
riodic extension in [1, 1
+
d
] that transitions smoothly from
η
(1)
back to
η
(0) (similarly for
u
). The resulting continued functions
can be viewed as sets of discrete values of periodic and smooth
functions that can be approximated to high-order on slightly larger
intervals by a trigonometric polynomial. Once these discrete peri-
odic continuation functions have been constructed, corresponding
Fourier coefficients
a
k
,
b
k
in eq. (
2
) can be obtained rapidly from
applications of the Fast Fourier Transform (FFT). The adopted FC
parameters employed in this work as well as a detailed presenta-
tion on the accelerated construction of FC functions can be found
in Amlani & Bruno (
2016
).
Using these discrete continuations in order to evaluate spatial
function values and derivatives on the discretized physical domain
modelled by the shallow water wave equations, the algorithm is com-
pleted by using the explicit fourth-order Adams-Bashforth scheme
(Amlani & Bruno
2016
; Amlani & Pahlevan
2020
;Amlani
et al.
2021
) to integrate the corresponding ordinary differential equa-
tions in time from the given initial conditions
η
(
y
i
,
t
)
=
u
(
y
i
,
t
)
=
0 up to a final given time. The final full solver enables high-
order accuracy and nearly dispersionless resolution of propagating
waves with mild, linear Courant–Friedrichs–Lewy constraints on the
temporal discretization—properties that are important for adequate
resolution of the different spatial and temporal scales involved be-
tween the supershear source dynamics and the subsequent tsunami
dynamics. Both implicit and explicit FC-based partial differential
equation solvers have been successfully constructed and utilized for
a variety of physical problems including those governed by radia-
tive transfer equations (Gaggioli
et al.
2019
), classical wave and
diffusion equations (Lyon & Bruno
2010
; Bruno & Prieto
2014
),
Euler equations (Shahbazi
et al.
2013
), convection-diffusion equa-
tions (Amlani
et al.
2021
), Navier–Cauchy elastodynamics equa-
tions (Amlani & Bruno
2016
;Amlani
et al.
2019
), Navier–Stokes
fluid equations (Albin & Bruno
2011
; Bruno
et al.
2019
; Fontana
et al.
2020
) and fluid-structure hemodynamics equations (Amlani
& Pahlevan
2020
).
2.3 GPS data from the Palu earthquake
The dual-frequency GPS is processed using the scientific GIPSY-
OASIS II software version 6.4 (Webb
1997
). The (post-processing)
Precise Point Positioning (PPP) method (Zumberge
et al.
1997
)
is used in kinematic (1 s) mode to derive precise absolute coor-
dinates for the PALP station. Precise ephemeris of GPS satellites
(non-fiducial style, using high-rate 30 s clocks) along with Earth
rotation parameters (ERP) in the IGS14 framework (Rebischung
& Schmid
2016
) are obtained from the Jet Propulsion Laboratory
(JPL). A satellite elevation mask angle of 7
◦
and absolute Inter-
national GNSS Service (IGS) antenna phase centre corrections are
applied. The Vienna tropospheric Mapping Functions (VMF1) are
used in estimating both zenith delay and gradients, downloaded from
the Global Geodetic Observing System website (re3data.org: VMF
Data Server
2021
). The global ocean tide model applied in the GPS
data processing is FES2014b, and the ocean loading parameters have
been retrieved from the Onsala Space Observatory website (M.S.
Bos and H.-G. Scherneck,
http://holt.oso.chalmers.se/loading/
). To
enhance the coordinate solutions, the daily global wide lane phase
bias (wlpb) files from JPL are used to resolve the phase cycle am-
biguities (Bertiger
et al.
2010
). Although each kinematic position
has a higher uncertainty and is affected by biases which usually
cancel out over long periods of measurements, the instantaneous
coseismic displacements at PALP are much higher than the high-
frequency noise of around 1 cm and 2–3 cm for, respectively, the
horizontal and vertical positions. Finally, the GPS time tags are cor-
rected to UTC time by subtracting 18 s. The coseismic displacement
of the station simply follows from epoch-to-epoch coordinate differ-
ences. The standard available script has been modified to properly
weigh the phase/code measurements of the stations and also to out-
put the correlations. The
X
-
Y
-
Z
Cartesian component positions are
then converted to the north-east-up positions along with their for-
mal standard deviations. They are scaled using the weighted-root-
mean-square of all the positions up to the time of the earthquake
and generally reach a relative precision (3
σ
) of about 30 mm on
the horizontal components. The resulting displacement field is then
differentiated by computing adaptive linear fits adapted to satisfy
an error to fit criteria. The slope of the linear fit then gives the local
velocity. The resulting data is then resampled again at 1 Hz by lin-
ear interpolation. The corresponding velocity data is presented in
Section 3.
3 RESULTS AND DISCUSSION
3.1 Direct evidence of a supershear rupture
In this section we provide the first-ever observation of supershear
by a high-rate GPS station, accomplished by considering the most
unmistakable signature of a supershear rupture: that the fault parallel
particle velocity dominates over the fault normal velocity (Dunham
& Archuleta
2005
; Mello
et al.
2014,
when the rupture velocity
v
is greater than
√
2
c
s
for a shear wave speed
c
s
). The opposite
signature is expected for a subshear rupture. Fig.
1
(a) shows the
Palu-Koro fault system (comprising of three segments) with the
location of the high-rate, 1 Hz, PALP GPS station. Figs
1
(b) and (c)
show the particle velocities recorded during the Palu earthquake,
clearly demonstrating a fault parallel particle velocity greater than
the fault normal velocity (
∼
1.0 m s
–1
versus
∼
0.7 m s
–1
). This proves
that the rupture, as it passed by the PALP station, definitively went
supershear and hence attained a speed between
√
2
c
s
and the
P
-wave
speed,
c
p
, of the medium (the absolute limiting speed of the rupture).
Socquet
et al.
(
2019
)andBao
et al.
(
2019
) have also inferred
that this earthquake went supershear, but mainly through far-field
observations employing geodetic and teleseismic data, respectively.
The only other near-field evidence of a supershear earthquake was
obtained using an accelerometer (250 Hz) at Pump Station 10 (PS10)
during the 2002
M
w
7.9 Denali earthquake (Ellsworth
et al.
2004
;
Downloaded from https://academic.oup.com/gji/article/230/3/2089/6604081 by California Institute of Technology user on 23 June 2022
2092
F. Amlani et al.
Figure 1.
The earthquake rupture and near-field evidence of supershear. (a) The Palu-Koro fault system, where the Pantoloan tidal gauge and the PALP GPS
station are marked. The green line of dots represents the slice of the bay considered for the tsunami model used in this study. (b) Comparison between th
e
fault parallel particle velocities recorded at the PALP station with those generated by the numerical supershear rupture model (Dunham & Bhat
2008
). (c)
Comparison between the corresponding fault normal particle velocities. (d, e) Same as (b, c) but for a subshear rupture.
Mello
et al.
2014
). We emphasize here that we have not performed
any kinematic inversion of the GPS data; we instead have employed
well-known unique signatures of near source ground velocity for
supershear ruptures (Dunham & Archuleta
2005
; Mello
et al.
2014
)
that indubitably confirm that the rupture, at least as it passed by the
PALP station, was supershear.
We can further compare the PALP records against a 3-D supers-
hear earthquake simulation (Section 2.1) whose rupture propagates
at a speed of
v
=
1.6
c
s
and whose corresponding particle velocities
are computed at 100 Hz and then decimated to match the 1 Hz
sampling rate of the GPS observations. The synthetic data and the
GPS records are in excellent agreement for the main rupture pulse
(Figs
1
b and c). Subsequent arrivals are not as well-captured since
the numerical model does not account for local velocity structure
nor detailed fault geometry. A similar comparison with synthetic
velocities computed for a subshear rupture (
v
=
0.8
c
s
) finds that
they are in poor agreement with GPS data (Figs
1
d and e). This
clearly suggests that the supershear rupture speed was 1.6
c
s
(around
5.3 km s
–1
) when it passed by PALP (Ulrich
et al.
2019
,alsofind
a speed greater than
√
2
c
s
). We have thus provided the definitive
first near-field high-rate GPS-based proof that the earthquake rup-
ture actually did go supershear as claimed and, further, have vali-
dated the numerical data used to source the tsunami model in what
follows.
3.2 Capturing the arrival and first motions at Pantoloan
The specific Palu bay configuration is outlined in Fig.
2
along the
horizontal
y
-axis, where
z
=
η
(
y
,
t
) represents the water height
relative to the background sea level. The bathymetry shape closely
approximates that of the segment demarcated by the dotted green
line near the Pantoloan tidal gauge in Fig.
1
(a) [basin width 9.2 km,
maximum depth 710 m and an average slope of 7
◦
to the east and
27
◦
to the west of the bay (Weatherall
et al.
2015
)]. The shallowest
part is taken to be 1 m, and the distance between the virtual gauge
and the fault is 4.3 km. The complete computational domain is taken
to be twice the basin width (
L
=
18.4 km).
Fig.
3
presents a temporal snapshot in the (
x
,
y
)-plane (the ground
surface) illustrating the dynamic vertical velocity field (and associ-
ated Mach fronts) which is input as a synthetic source in conjunction
with its corresponding time-dependent displacement field. The fault
and the sense of slip (left-lateral) are indicated in red, and the data
applied to perturb the bathymetry is taken along the green dotted
line (whose locations correspond to the same markers indicated in
Fig.
2
). For an example point located at (
x
0
,
y
0
) and highlighted in
a larger light green circle, Fig.
3
additionally presents the temporal
evolution of both the vertical velocity (which can reach
∼
1ms
–1
along the domain) as well as its corresponding ground displacement
(which, in the 1-D setting, can reach
∼
40 cm). As already noted,
the shapes and the maximum values of these profiles remain fairly
unattenuated at large distances from the original earthquake—a
hallmark of the energy carried by supershear shock fronts (Bernard
& Baumont
2005
; Dunham & Bhat
2008
).
For the results that follow, Fig.
4
additionally presents the anal-
ogous inputs for classical modelling of seismogenic tsunamis. In a
classical setting (Pedlosky
2013
), the source is often modelled as a
static displacement perturbation applied to the bathymetry (rather
than dynamic ground motion), that is a static
h
(
y
,
t
)
=
h
(
y
)that
neither accounts for the time-dependence nor the velocity of the
seafloor (other simple approximations to more complicated sources
are also standard (Kajiura
1963
; Tanioka & Satake
1996
). From the
supershear earthquake results, this corresponds to the final, perma-
nent ground displacement at the end of the profiles in Fig.
3
and is
expectedly on the order of a few centimetres.
Downloaded from https://academic.oup.com/gji/article/230/3/2089/6604081 by California Institute of Technology user on 23 June 2022