Supershear shock front contributions to the tsunami from the 2018
$\mathbf{M}_{\mathbf{w}}$ 7.5 Palu earthquake

Faisal Amlani,1 Harsha S. Bhat,2,∗ Wim J. F. Simons,3 Alexandre Schubnel,2
Christophe Vigny2, Ares J. Rosakis4, Joni Efendi5, Ahmed Elbanna6 and
Hasanuddin Z. Abidin5,7

  1. 1.

Department of Aerospace and Mechanical Engineering, University of Southern
California, Los Angeles, USA.

  2. 2.

Laboratoire de Géologie, École Normale Supérieure, CNRS-UMR 8538, PSL Research
University, Paris, France.

  3. 3.

Faculty of Aerospace Engineering, Delft University of Technology, Delft,
Netherlands.

  4. 4.

Graduate Aerospace Laboratories, California Institute of Technology, Pasadena,
California, USA.

  5. 5.

BIG (Badan Informasi Geospasial / Geospatial Information Agency of Indonesia),
Java, Indonesia.

  6. 6.

Department of Civil and Environmental Engineering, University of Illinois at
Urbana Champaign, USA.

  7. 7.

Department of Geodesy and Geomatics Engineering, Institute of Technology
Bandung, Indonesia.

  8. *

Corresponding author<EMAIL_ADDRESS>
Hazardous tsunamis are known to be generated predominantly at subduction zones
by large earthquakes on dip (vertical)-slip faults. However, a moment
magnitude ($\mathbf{M}_{\mathbf{w}}$) 7.5 earthquake on a strike
(lateral)-slip fault in Sulawesi (Indonesia) in 2018 generated a tsunami that
devastated the city of Palu. The mechanism by which this large tsunami
originated from a strike-slip earthquake has been debated. Here we present
near-field ground motion data from a GPS station that confirms that the 2018
Palu earthquake attained supershear speed, i.e., a rupture speed greater than
the speed of shear waves in the host medium. We study the effect of this
supershear rupture on tsunami generation by coupling the ground motion to a 1D
non-linear shallow-water wave model that accounts for both the time-dependent
bathymetric displacement and velocity. With the local bathymetric profile of
the Palu bay around a tidal gauge, we find that these simulations reproduce
the tsunami motions measured by the gauge, with only minimal tuning of
parameters. We conclude that Mach (shock) fronts, generated by the supershear
speed of the earthquake, interacted with the bathymetry and contributed to the
tsunami. This suggests that rupture speed should be considered in tsunami
hazard assessments.

Tsunamis are well-known to be amongst the most destructive consequences of
earthquakes [3, 4, 2, 1], and the 2018 Sulawesi earthquake was no exception:
it generated a devastating tsunami[5, 6] in the nearby Palu bay in which
hundreds were killed and tens of thousands more displaced from their homes
(ASEAN Situation Update No.15 - Earthquake & Tsunami Sulawesi). However, this
was an unexpected event since the earthquake was associated with the
predominantly in-plane ground motion produced by strike-slip ruptures. As
these motions are not known to excite significant waves, the underlying
physical mechanisms behind the tsunami have largely remained a mystery [7].
Many studies conducted to explain the phenomenon have not arrived at
definitive conclusions [8] nor have adequately captured tidal gauge records
[9, 10]; the main concensus appears to be that some form of ground motion
(e.g., landslides [11] or the reverse-slip motion of the fault [12]),
amplified by the bay, is to blame.

However, a key notable feature of this earthquake is that it ruptured at
supershear speed [13, 14], which results in a manifestation of two shock (or
Mach) fronts carrying significant vertical velocity with relatively slow
attenuation over large distances [15]. The existence of supershear earthquakes
has been proven theoretically and experimentally since the early 1970s [16,
17, 18, 19, 20, 21, 22]. The $M_{w}$ 6.5 Imperial Valley earthquake
(California, 1979) was the first naturally observed supershear earthquake
rupture [23]. Since then, many more (although rare) earthquakes have been
recorded to propagate at supershear speeds: the $M_{w}$ 7.4 1999 Izmit in
Turkey [24], the $M_{w}$ 7.8 2001 Kunlun[25] and the $M_{w}$ 7.8 2002 Denali
[26, 27], to name a few.

Although the overall tsunami behaviour at Palu is likely a combination of
several effects that include these supershear dynamics as well as landslides,
recent studies [29, 28, 9] 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 2D/3D geometric features of the Palu bay. Jamelot et
al. [9] capture amplitudes recorded by the Pantoloan tidal gauge, but not the
first phases and motions—ultimately conceding that a dynamic study should be
conducted to better understand the influence of supershear. The model-based
study of Ulrich et al. [29], which incorporates some dynamics of supershear in
the form of time-dependent ground displacement, better captures tidal gauge
records but, again, there remains a mismatch in the first phases and arrival.

Hence the primary objective of this work is to provide the missing link in
explaining these discrepancies and more fully understand the role played by
supershear rupture dynamics on the generation of the Palu tsunami. In
particular, we incorporate a neglected feature in the above-cited works that
is a defining characteristic of supershear earthquakes: the _velocity_ of the
ground motion [15]. Using a model validated by the first near-field evidence
of supershear at Palu, our results imply that ground speeds, which better
represent the intricacies of the Mach fronts, may further explain the observed
motions of the tsunami. Since other studies (including those investigating
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 first phases and arrival in the Pantoloan records.

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 employed in this study. b, Comparison between
the fault parallel particle velocities recorded at the PALP station with those
generated by the numerical supershear rupture model [15]. c, Comparison
between the corresponding fault normal particle velocities. d,e, Same as (b,c)
but for a subshear rupture.

### Evidence of a supershear rupture

The most unmistakable signature of a supershear rupture is that the fault
parallel particle velocity dominates over the fault normal velocity [30, 27]
(when the rupture velocity $v$ is greater than $\sqrt{2}c_{s}$ for a shear
wave speed $c_{s}$). The opposite signature is expected for a subshear
rupture. Fig. 1a shows the Palu-Koro fault system with the location of the
high-rate, $1$ Hz, PALP GPS station. Figs. 1b-c show the particle velocities
recorded during the Sulawesi earthquake, clearly demonstrating a fault
parallel particle velocity greater than the fault normal velocity
($\sim\\!\\!1.0$ m/s versus $\sim\\!\\!0.7$ m/s). This proves that the
rupture, as it passed by the PALP station, definitively went supershear and
hence attained a speed between $\sqrt{2}c_{s}$ and the P-wave speed, $c_{p}$,
of the medium (the absolute limiting speed of the rupture). This represents
the first-ever observation of a supershear rupture by a high-rate GPS station.
Socquet et al. [14] and Bao et al. [13] 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 Mw 7.9 Denali earthquake
[26, 27].

We can further compare the PALP records against a 3D supershear earthquake
simulation [15] whose rupture propagates at a speed of $v=1.6c_{s}$ and whose
corresponding particle velocites are computed at $100$ Hz and then decimated
to match the $1$ Hz sampling rate of the GPS observations (see Methods). The
synthetic data and the GPS records are in excellent agreement for the main
rupture pulse (Figs. 1b-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.8c_{s}$) finds that they are in poor agreement with GPS
data (Figs. 1d-e). This clearly suggests that the supershear rupture speed was
$1.6c_{s}$ (around $5.3$ km/s) when it passed by PALP (Ulrich et al. [29] also
find a speed greater than $\sqrt{2}c_{s}$). We have thus provided the first
near-field high-rate GPS-based evidence that the Sulawesi earthquake rupture
actually did go supershear as claimed and, further, have validated the
numerical data employed to source the tsunami model in what follows.

Figure 2: Dynamic rupture modelling. Representative synthetic particle
displacement (left) and velocity (right) histories from Dunham and Bhat[15].
Red lines correspond to the location of the Palu bay.

### Modelling the effect of supershear velocity on tsunami generation

Using the synthetic particle motions (which, again, agree with PALP GPS
records and are generated by the 3D supershear earthquake model), a 1D non-
linear shallow water wave model incorporating time-dependent ground movements
of velocity and displacement [31] has been utilised to simulate the generation
and propagation of the tsunami. This employs the depth-averaged shallow water
approximation of the Euler equations that can be written as a system of
coupled hyperbolic partial differential equations given by

$\begin{cases}\dfrac{\partial H}{\partial t}+\dfrac{\partial(Hu)}{\partial
y}=0\\\ \dfrac{\partial(Hu)}{\partial t}+\dfrac{\partial(Hu^{2})}{\partial
y}+gH\dfrac{\partial\eta}{\partial y}=0\end{cases},\quad 0\leq y\leq L,\quad
t\geq 0.$ (1)

Figure 3: The non-linear tsunami model setup that incorporates displacement
and velocity ground dynamics. The strike-slip fault and its sense of motion
are indicated in red in all the panels. a, A diagram of the non-linear shallow
water wave system for tsunami height $\eta$, initial bathymetry $h_{0}$ and
bathymetry perturbation (source) $h$. b, Snapshot of the dynamic vertical
velocity from a supershear earthquake with, c, its temporal evolution at an
example point $(x_{0},y_{0})$ (light green). d, The static displacement field
due to a supershear earthquake. The dark green dots on the supershear
earthquake data in (b,d) correspond to the source locations used to perturb
the bathymetry domain in (a). e, The spatial profile in $y$ of the static
displacement field due to a supershear earthquake.

Here, $u(y,t)$ is the fluid velocity, $\eta(y,t)$ is the sea surface height
and $H(y,t)=h_{0}(y)+h(y,t)+\eta(y,t)$ is the absolute height from the bed-
level to the water surface for an initial at-rest bathymetry $h_{0}(y)$. The
entire domain of length $L$ is subjected to a time-dependent ground
perturbation $h(y,t)$ which—together with the corresponding ground velocity
$\partial h(y,t)/\partial t$ included in system (1)—sources the subsequent
tsunami dynamics. In what follows, these values are determined from the 3D
supershear earthquake model [15] and illustrated in the seismograms of Fig. 2.
The constant $g$ is the acceleration due to gravity.

The specific Palu bay configuration is outlined in Fig. 3a along the
horizontal $y$-axis, where $z=\eta(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. 1a (basin width 9.2 km, maximum depth 710 m and an average slope of
$7^{\circ}$ to the east and $27^{\circ}$ to the west of the bay). 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. 3b presents a temporal snapshot
in the $(x,y)$-plane (the ground surface) illustrating the dynamic vertical
velocity field (and associated 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. 3a).
For an example point located at $(x_{0},y_{0})$ and highlighted in a larger
light green circle, Fig. 3c additionally presents the temporal evolution of
both the vertical velocity (which can reach $\sim\\!\\!1$ m/s along the
domain) as well as its corresponding ground displacement (which, in the 1D
setting, can reach $\sim\\!\\!40$ cm). As already noted (and further
illustrated in Fig. 2), 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 [15]. For the
results that follow, Figs. 3d-e additionally present the analogous inputs for
classical modelling of seismogenic tsunamis. In a classical setting [32], the
source is often modelled as a static displacement perturbation applied to the
bathymetry (rather than dynamic ground motion), i.e., a static $h(y,t)=h(y)$
that neither accounts for the time-dependence nor the velocity of the sea
floor (other simple approximations to more complicated sources are also
standard [33, 34]). From the supershear earthquake results, this corresponds
to the final, permanent ground displacement at the end of the profiles in Fig.
3c and is expectedly on the order of a few centimeters.

### Capturing first motions and arrival recorded at Pantoloan

Numerical solution of the non-linear shallow water wave equations has been
facilitated by a spectral Fourier continuation (FC) methodology [35, 36]
employing the bathymetry described above (see also Figs. 1,3). Such a solver
enables high-order accuracy, mild constraints on the temporal discretisation
and nearly dispersionless resolution of propagating waves over large distances
[35] (see Methods). Fig. 4a presents simulation results in the $(t,y)$-plane
(time and space) of the water height $z=\eta(y,t)$ and, more importantly,
presents a comparison between observations recorded every minute by the
Pantoloan tidal gauge (whose geographic location is indicated in Fig. 1a) with
that predicted by the dynamic source model (plotted in Fig. 4a at the same
frequency). The numerical modelling has been conducted at a much higher
temporal resolution (on the order of 10 milliseconds), and the Pantoloan
records have been obtained and processed by the Agency for Geospatial
Information (BIG), Indonesia (http://tides.big.go.id). Remarkably, both the
first motions and phases from the observation records are in excellent
agreement with the 1D approximation generated by excitation from the dynamic
supershear earthquake. Later phases, which can be attributed to wave
reflections within the bay, are not as well-captured since our model does not
fully account for the localised effects of the 2D/3D bathymetric profile.
Nevertheless, the tsunami arrival and primary dynamics are correctly
reproduced.

Figure 4: Simulated tsunamis from dynamic and classical modelling. a,b,
Spatiotemporal evolution $z=\eta(y,t)$ of the tsunami generated by a
supershear earthquake along the Palu bay due to a, dynamic bathymetry velocity
with dynamic displacement and b, static displacement. Normalised Pantoloan
tidal gauge records during the event (grey) are overlaid with model
predictions (blue)

For comparison, Fig. 4b presents the corresponding simulations that are
classically-sourced by the final (static) vertical displacement given by the
same supershear earthquake simulation. Although a tsunami is generated, there
are large contrasts in phase, width and particularly the sign (Fig. 4b) when
compared to the Pantoloan tidal gauge observations. This implies that the
complex ground motion (which is not accounted for by classical tsunami models)
must be incorporated in order to correctly predict the tsunami dynamics as in
Fig. 4a (see Supplementary Figs. 1-2 for comparisons between the two model
sources). We note that we have normalized water heights by their corresponding
maxima for comparison throughout: since more energy of the Mach fronts is
carried _along_ the fault [15] running in the direction $x$ (Fig. 3), the 1D
model in $y$ will naturally generate lower amplitudes (of the order of half a
meter). However, similar tsunami signatures can still be expected and, indeed,
Jamelot et al. [9] incorporated a more comprehensive (classical) model
treating localised effects of the 2D/3D bathymetry but generated the same
mismatches in phase as those in Fig. 4b.

Hence we confirm that the Sulawesi earthquake went supershear via the first
near-field high-rate GPS-based evidence of such a rupture and that, by
modelling the corresponding effect on the generation of tsunamis in a shallow
geometry, we conclude that the ground motion resulting from the associated
Mach fronts (which carry minimally attenuated velocities to large distances)
may well have contributed to the Palu tsunami. This work provides a robust
proof of concept, albeit in 1D, on the contribution of shock fronts in tsunami
generation. In order to gain further insight into this process, more detailed
modeling in 3D is needed to account for, e.g., geometrical spreading,
attenuation and geological structure. Regardless, since nothing geologically
specific about the bay has been introduced, our results signify the importance
of such configurations for tsunami hazard assessment due to strike-slip
earthquakes. The same physical ingredients (supershear rupture and a shallow
bay) may combine to produce similar effects elsewhere, including the Tomales
bay in California (which is crossed offshore by the San Andreas fault system
[37]) and the Izmit bay in Turkey (which is crossed by the North Anatolian
fault [38]). Both these regions, as well as the Palu bay, have suffered from
historical tsunamis. On the contrary, the 2012 Off Northern Sumatra earthquake
and the 2013 Craig, Alaska earthquake both went supershear but caused
negligible (or no) tsunamis since they occurred in deep ocean without any
shallow bay near them. Additionally, the 1999 Izmit earthquake was subshear as
it passed through the Izmit bay and thus generated only a negligible tsunami.
Hence we reemphasise that supershear rupture and a shallow bay are key to
generate contributions to tsunami motions. We thus suggest that any rapid
assessment of tsunami hazard after a strike-slip earthquake should also
involve a rapid assessment of the earthquake rupture velocity as we have shown
that ultimately the focal mechanism, the depth and the speed of the rupture
all contribute towards the generation of tsunamis.

## Methods

#### GPS data.

The dual-frequency GPS has been processed using the scientific GIPSY-OASIS II
software version 6.4 (Jet Propulsion Laboratory, 2017). The (post-processing)
Precise Point Positioning (PPP) method [39] was used in kinematic ($1$ s) mode
to derive precise absolute coordinates 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[40] were obtained
from the Jet Propulsion Laboratory (JPL). A satellite elevation mask angle of
$7$ degrees and absolute IGS antenna phase centre corrections were applied.
The Vienna tropospheric Mapping Functions (VMF1) was used (in estimating both
zenith delay and gradients) and downloaded from the Global Geodetic Observing
System website (http://vmf.geo.tuwien.ac.at/). The global ocean tide model
applied in the GPS data processing was (FES2014b), and the ocean loading
parameters were retrieved from the Onsala Space Observatory website
(http://holt.oso.chalmers.se/). To enhance the coordinate solutions, the daily
global wide lane phase bias (wlpb) files from JPL were used to resolve the
phase cycle ambiguities[41]. Although each kinematic position has a higher
uncertainty and is affected by biases which usually cancel out over long
periods of measurements, the instantaneous co-seismic displacements at PALP
were much higher than the high-frequency noise of around $1$ cm and $2$ to $3$
cm for, respectively, the horizontal and vertical positions. Finally the GPS
time tags were corrected to UTC time by subtracting $18$ s. The co-seismic
displacement of the station simply follows from epoch-to-epoch coordinate
differences. The standard available script was modified to properly weigh the
phase/code measurements of the stations and also to output the correlations.
The XYZ positions were then converted to the NEU positions along with their
formal standard deviations. They are scaled using the WRMS of all the
positions up to the time of the earthquake and generally reach a relative
precision ($3\sigma$) of about $30$ mm on the horizontal components. The
resulting displacement field was 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 was then
resampled again at $1$ Hz by linear interpolation.

#### Supershear earthquake dynamics and rupture modelling.

We have used the existing numerical simulations conducted by Dunham and
Bhat[15]. Additional details can be found in the cited manuscript. Numerical
simulations were conducted using a staggered-grid finite-difference (FD) code
[42] with the fault boundary conditions implemented using the staggered-grid
split-node (SGSN) method [43]. As the authors had provided non-dimensionalised
results we simply dimensionalized their results for the Palu earthquake by
using a shear modulus of $30$ GPa, stress drop of $10$ MPa and the shear wave
speed of $3.5$ km/s. The depth of the rupture was assumed to be $10$ km. The
resulting particle velocities and displacements are shown in Fig. 2.

#### Numerical analysis of shallow water wave equations.

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 [36, 48, 35]. 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,\dots,N-1$), Fourier
continuation algorithms append a small number of points to the discretised
function values $\eta(y_{i}),u(y_{i})$ in order to form ($1+d$)-periodic
trigonometric polynomials $\eta_{\text{cont}}(y),u_{\text{cont}}(y)$ that are
of the form

$\eta_{\text{cont}}(y)=\displaystyle\sum_{k=-M}^{M}a_{k}e^{\frac{2\pi
iky}{1+d}}~{}~{},~{}~{}u_{\text{cont}}(y)=\displaystyle\sum_{k=-M}^{M}b_{k}e^{\frac{2\pi
iky}{1+d}}$ (2)

and that match the given discrete values of $\eta(y_{i}),u(y_{i})$, i.e.,
$\eta_{\text{cont}}(y_{i})=\eta(y_{i}),u_{\text{cont}}(y_{i})=u(y_{i})$ for
$i=0,...,N-1.$ Spatial derivatives in equation (1) are then computed by exact
term-wise differentiation of (2) as

$\displaystyle\dfrac{\partial\eta}{\partial y}(y_{i})$
$\displaystyle=\dfrac{\partial\eta_{\text{cont}}}{\partial y}(y_{i})$
$\displaystyle=\sum_{k=-M}^{M}\left(\frac{2\pi
ik}{1+d}\right)a_{k}e^{\frac{2\pi iky_{i}}{1+d}},$ (3)
$\displaystyle\dfrac{\partial u}{\partial y}(y_{i})$
$\displaystyle=\dfrac{\partial u_{\text{cont}}}{\partial y}(y_{i})$
$\displaystyle=\sum_{k=-M}^{M}\left(\frac{2\pi
ik}{1+d}\right)b_{k}e^{\frac{2\pi iky_{i}}{L+d}}.$

In essence, FC algorithms add a handful of additional values to the original
discretized function in order to form a periodic extension in $[1,1+d]$ that
transitions smoothly from $\eta(1)$ back to $\eta(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
periodic continuation functions have been constructed, corresponding Fourier
coefficients $a_{k},b_{k}$ in equation (2) can be obtained rapidly from
applications of the Fast Fourier Transform (FFT). A detailed presentation on
accelerated construction of Fourier continuations can be found in, e.g.,
Amlani and Bruno [35].

Employing these discrete continuations in order to evaluate spatial function
values and derivatives on the discretised physical domain modelled by equation
(1), the algorithm is completed by employing the explicit fourth-order Adams-
Bashforth scheme [35] to integrate the corresponding ODEs in time from the
given initial conditions $\eta(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 (CFL) constraints on the temporal discretisation—properties
that are important for adequate resolution of the different spatial and
temporal scales involved between the supershear source dynamics and the
subsequent tsunami dynamics. Both implicit and explicit FC-based PDE solvers
have been successfully constructed and utilised for a variety of physical
problems including those governed by classical wave and diffusion equations
[36, 44], radiative transfer equations [45], non-linear Burgers systems [46],
Euler equations [47], Navier-Stokes fluid equations [48, 49, 50], Navier-
Cauchy elastodynamics equations [35] and fluid-structure hemodynamics
equations [51].

### Code availability

All codes are available upon request to the corresponding author.

### Data availability

All relevant data relating to the tsunami modelling and the PALP GPS station
are available on Zenodo (DOI:10.5281/zenodo.4066297).

## All References

  * [1] Röbke, B. R. & Vött, A.  The tsunami phenomenon.  _Prog. Oceanogr._ 159, 296–322 (2017).
  * [2] Pugh, D. & Woodworth, P.  _Sea-Level Science: Understanding Tides, Surges, Tsunamis And Mean Sea-Level Changes_ (Cambridge Univ. Press, Cambridge, 2014).
  * [3] Bryant, E.  _Tsunami: The Underrated Hazard_ (Springer-Verlag, Berlin Heidelberg, 2008).
  * [4] Synolakis, C. E. & Okal, E. A.  1992–2002: Perspective on a decade of post-tsunami surveys.  In _Tsunamis_ , 1–29 (Springer, Dordrecht, 2005).
  * [5] Mai, P. M.  Supershear tsunami disaster.  _Nat. Geosci._ 12, 150 (2019).
  * [6] Fritz, H. et al.  Field survey of the 28 September 2018 Sulawesi Tsunami.  In _AGU Fall Meeting Abstracts 2018_ NH22B-04 (2018).
  * [7] Syamsidik, Benazir, Umar, M., Margaglio, G. & Fitrayansyah, A.  Post-tsunami survey of the 28 September 2018 tsunami near Palu Bay in Central Sulawesi, Indonesia: Impacts and challenges to coastal communities.  _Int. J. Disast R. Re._ 38, 101229 (2019).
  * [8] Muhari, A. et al.  Solving the puzzle of the September 2018 Palu, Indonesia, tsunami mystery: Clues from the tsunami waveform and the initial field survey data.  _J. Disaster Res._ 13 Sci. Comm., sc20181108 (2018).
  * [9] Jamelot, A., Gailler, A., Heinrich, P., Vallage, A. & Champenois, J.  Tsunami simulations of the Sulawesi Mw 7.5 Event: Comparison of seismic sources issued from a tsunami warning context versus post-event finite source.  _Pure Appl. Geophys._ 176, 3351–3376 (2019).
  * [10] Heidarzadeh, M., Muhari, A. & Wijanarto, A. B.  Insights on the source of the 28 September 2018 Sulawesi tsunami, Indonesia based on spectral analyses and numerical simulations.  _Pure Appl. Geophys._ 176, 25–43 (2019).
  * [11] Sassa, S. & Takagawa, T.  Liquefied gravity flow-induced tsunami: first evidence and comparison from the 2018 Indonesia Sulawesi earthquake and tsunami disasters.  _Landslides_ 16, 195–200 (2019).
  * [12] He, L. et al.  Source parameters and slip distribution of the 2018 Mw 7.5 Palu, Indonesia earthquake estimated from space-based geodesy.  _Tectonophysics_ 772, 228216 (2019).
  * [13] Bao, H. et al.  Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake.  _Nat. Geosci._ 12, 200–205 (2019).
  * [14] Socquet, A., Hollingsworth, J., Pathier, E. & Bouchon, M.  Evidence of supershear during the 2018 magnitude 7.5 Palu earthquake from space geodesy.  _Nat. Geosci._ 12, 192–199 (2019).
  * [15] Dunham, E. M. & Bhat, H. S.  Attenuation of radiated ground motion and stresses from three-dimensional supershear ruptures.  _J. Geophys. Res._ 113, 1–17 (2008).
  * [16] Burridge, R.  Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion.  _Geophys. J. Roy. Astr. S._ 35, 439–455 (1973).
  * [17] Andrews, D. J.  Rupture velocity of plane strain shear cracks.  _J. Geophys. Res._ 81, 5679–5689 (1976).
  * [18] Das, S. & Aki, K.  A numerical study of two-dimensional spontaneous rupture propagation.  _Geophys. J. Roy. Astr. S._ 50, 643–668 (1977).
  * [19] Wu, F. T., Thomson, K. & Kuenzler, H.  Stick-slip propagation velocity and seismic source mechanism.  _Bull. Seism. Soc. Am._ 62, 1621–1628 (1972).
  * [20] Rosakis, A. J., Samudrala, O. & Coker, D.  Cracks faster than the shear wave speed.  _Science_ 284, 1337–1340 (1999).
  * [21] Xia, K. W., Rosakis, A. J. & Kanamori, H.  Laboratory earthquakes: The sub-Rayleigh-to-supershear rupture transition.  _Science_ 303, 1859–1861 (2004).
  * [22] Passelègue, F. X., Schubnel, A., Nielsen, S., Bhat, H. S. & Madariaga, R.  From sub-Rayleigh to supershear ruptures during stick-slip experiments on crustal rocks.  _Science_ 340, 1208–1211 (2013).
  * [23] Archuleta, R. J.  A faulting model for the 1979 Imperial Valley earthquake.  _J. Geophys. Res._ 89, 4559–4586 (1984).
  * [24] Bouchon, M. et al.  How fast is rupture during an earthquake? New insights from the 1999 Turkey earthquakes.  _Geophys. Res. Lett._ 28, 2723–2726 (2001).
  * [25] Robinson, D. P., Brough, C. & Das, S.  The ${M}_{w}$ 7.8, 2001 Kunlunshan earthquake: Extreme rupture speed variability and effect of fault geometry.  _J. Geophys. Res._ 111, B08303 (2006).
  * [26] Ellsworth, W. L. et al.  Near-field ground motion of the 2002 Denali fault, Alaska, earthquake recorded at pump station 10.  _Earthq. Spectra_ 20, 597–615 (2004).
  * [27] Mello, M., Bhat, H. S., Rosakis, A. J. & Kanamori, H.  Reproducing the supershear portion of the 2002 Denali earthquake rupture in laboratory.  _Earth Planet. Sc. Lett._ 387, 89–96 (2014).
  * [28] Oral, E., Weng, H. & Ampuero, J. P.  Does a damaged-fault zone mitigate the near‐field impact of supershear earthquakes?–Application to the 2018 $M_{W}$ 7.5 Palu, Indonesia, earthquake.  _Geophys. Res. Lett._ 47, e2019GL085649 (2020).
  * [29] Ulrich, T. et al.  Coupled, physics-based modelling reveals earthquake displacements are critical to the 2018 Palu, Sulawesi tsunami.  _Pure Appl. Geophys._ 176, 4069–4109 (2019).
  * [30] Dunham, E. M. & Archuleta, R. J.  Near-source ground motion from steady state dynamic rupture pulses.  _Geophys. Res. Lett._ 32, L03302 (2005).
  * [31] Dutykh, D. & Clamond, D.  Modified shallow water equations for significantly varying seabeds.  _Appl. Math. Model._ 40, 9767–9787 (2016).
  * [32] Pedlosky, J.  _Geophysical Fluid Dynamics_ (Springer-Verlag, New York, 2013).
  * [33] Tanioka, Y. & Satake, K.  Tsunami generation by horizontal displacement of ocean bottom.  _Geophys. Res. Lett._ 23, 861–864 (1996).
  * [34] Kajiura, K.  The leading wave of a tsunami.  _Bull. Earthquake Res. Inst._ 141, 535-571 (1963).
  * [35] Amlani, F. & Bruno, O. P.  An FC-based spectral solver for elastodynamic problems in general three-dimensional domains.  _J. Comput. Phys._ 307, 333–354 (2016).
  * [36] Lyon, M. & Bruno, O. P.  High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations.  _J. Comput. Phys._ 229, 3358–3381 (2010).
  * [37] Johnson, S. Y. & Beeson, J. W.  Shallow Structure and Geomorphology along the Offshore Northern San Andreas Fault, Tomales Point to Fort Ross, California.  _Seismol. Soc. Am., Bull._ 109, 833–854 (2019).
  * [38] Altinok, Y. et al.  The tsunami of August 17, 1999 in Izmit bay, Turkey.  _Nat. Hazards_ 24, 133–146 (2001).
  * [39] Zumberge, J. F., Heflin, M. B., Jefferson, D. C., Watkins, M. M. & Webb, F. H.  Precise point positioning for the efficient and robust analysis of GPS data from large networks.  _J. Geophys. Res._ 102, 5005–5017 (1997).
  * [40] Rebischung, P. & Schmid, R.  IGS14/igs14.atx: a new framework for the IGS products.  In _AGU Fall Meeting Abstracts 2016_ G41A-0998 (2016).
  * [41] Bertiger, W. et al.  Single receiver phase ambiguity resolution with GPS data.  _J. Geod._ 84, 327–337 (2010).
  * [42] Favreau, P., Campillo, M. & Ionescu, I. R.  Initiation of shear instability in three-dimensional elastodynamics.  _J. Geophys. Res._ 107, 2147–2164 (2002).
  * [43] Dalguer, L. A. & Day, S.  Staggered-grid split-node method for spontaneous rupture simulation.  _J. Geophys. Res._ 112, B02302 (2007).
  * [44] Bruno, O. P. & Prieto, A.  Spatially dispersionless, unconditionally stable FC-AD solvers for variable-coefficient PDEs.  _J. Sci. Comput._ 58, 1–36 (2014).
  * [45] Gaggioli, E. L., Bruno, O. P. & Mitnik, D. M.  Light transport with the equation of radiative transfer: The Fourier Continuation - Discrete Ordinates (FC-DOM) Method.  _J. Quant. Spectrosc. Radiat. Transf._ 236, 106589 (2019).
  * [46] Bruno, O. P. & Jimenez, E.  Higher-order linear-time unconditionally stable alternating direction implicit methods for nonlinear convection-diffusion partial differential equation systems.  _J. Fluid. Eng._ 136, 060904 (2014).
  * [47] Shahbazi, K., Hesthaven, J. S. & Zhu, X.  Multi-dimensional hybrid Fourier continuation–WENO solvers for conservation laws.  _J. Comput. Phys._ 253, 209–225 (2013).
  * [48] Albin, N. & Bruno, O. P.  A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping.  _J. Comput. Phys._ 230, 6248–6270 (2011).
  * [49] Bruno, O. P., Cubillos, M. & Jimenez, E.  Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers.  _J. Comput. Phys._ 391, 322–346 (2019).
  * [50] Fontana, M., Bruno, O. P., Mininni, P. D. & Dmitruk, P.  Fourier continuation method for incompressible fluids with boundaries.  _Comput. Phys. Commun._ 256, 107482 (2020).
  * [51] Amlani, F. & Pahlevan, N.M.  A stable high-order FC-based methodology for hemodynamic wave propagation.  _J. Comput. Phys._ 405, 109130 (2020).

## Corresponding author

Correspondance and requests for materials should be delivered to
Harsha S. Bhat<EMAIL_ADDRESS>
## Acknowledgements

F. A. acknowledges Prof. Oscar Bruno at Caltech for inspiring the numerical
technique used in this study. A. S. would like to acknowledge the European
Research Council grant REALISM (2016‐grant 681346). A.J.R. was supported by
the Caltech/MCE Big Ideas Fund (BIF). The continued (long-term) operation of
the GPS stations in Central Sulawesi, Indonesia data has been co-facilitated
by the EU-ASEAN SEAMERGES (2004-2006) and GEO2TECDI-1/2 projects (2009-2013)
in cooperation with the Geospatial Information Agency of Indonesia (BIG). The
GPS data acquisition and research was also partly funded by grants from the
Dutch NWO User Support Programme Space Research (2007-2018). We would like to
express our special thanks to the local staff of the Indonesian Meteorology,
Climatology, and Geophysical Agency (BMKG) offices in Palu, both for hosting
the GPS station and for being available 24/7 to assist us with optimal
operation of the GPS equipment. A special thanks to our local survey staff,
Bpak Rahman Umar and Urif, for their continued support, including directly
after the earthquake for inspecting the GPS stations. H. S. B. would like to
acknowledge Kiara and Aziliz for inspiring this work.

## Author contributions

F. A. and H. S. B. designed and conceived the project and conducted numerical
modelling studies. W. J. F. S., J. E. and C. V. have maintained the (high-
rate) PALP station and obtained all GPS data. W. J. F. S. did the kinematic
GPS data analysis. F. A. and H. S. B. wrote the manuscript. All authors
contributed to analysis, interpretation and manuscript preparation.