# Supershear shock front contribution to the tsunami from the 2018 Mw 7.5 Palu, Indonesia earthquake

Faisal Amlani, Harsha S Bhat, Wim J F Simons, Alexandre Schubnel, Christophe Vigny, Ares J Rosakis, Joni Efendi, Ahmed E Elbanna, Pierpaolo Dubernet, Hasanuddin Z Abidin

► To cite this version:

Faisal Amlani, Harsha S Bhat, Wim J F Simons, Alexandre Schubnel, Christophe Vigny, et al. Supershear shock front contribution to the tsunami from the 2018 Mw 7.5 Palu, Indonesia earthquake. Geophysical Journal International, 2022, 230, pp.2089 - 2097. 10.1093/gji/ggac162. hal-03869566

HAL Id: hal-03869566

https://ens.hal.science/hal-03869566

Submitted on 24 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_1.png" /></figure>

# Supershear shock front contribution to the tsunami from the 2018 Mw 7.5 Palu, Indonesia earthquake

Faisal Amlani<sup>1,*</sup> Harsha S. Bhat<sup>2</sup> Wim J.F. Simons,<sup>3</sup> Alexandre Schubnel,<sup>2</sup> Christophe Vigny,<sup>2</sup> Ares J. Rosakis,<sup>4</sup> Joni Efendi,<sup>5</sup> Ahmed E. Elbanna,<sup>6</sup> Pierpaolo Dubernet<sup>2</sup> and Hasanuddin Z. Abidin<sup>5,7</sup>

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<sub>w</sub> 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 modelling; Numerical approximations and analysis.

## 1 INTRODUCTION

Tsunamis are well-known to be amongst the most destructive consequences of earthquakes (Synolakis & Okal 2005; Bryant 2008; Pugh & Woodworth 2014; Röbke & Vött 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 hundreds were killed and tens of thousands more displaced from their homes (ASEAN 2018). However, this was a very unexpected occurrence since the earthquake was associated with the predominantly in-plane ground motion produced by strike-slip ruptures.

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 ruptured 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 attenuation over large distances (Bernard & Baumont 2005; Dunham & Bhat 2008). The existence of supershear earthquakes has been

<footnote n="1"><sup>1</sup>Department of Aerospace and Mechanical Engineering, University of Southern California, Los Angeles, CA 90089, USA</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Géologie, École Normale Supérieure, CNRS-UMR 8538, PSL Research University, 75231 Paris, France. E-mail: harsha.bhat@ens.fr</footnote>

<footnote n="3"><sup>3</sup>Faculty of Aerospace Engineering, Delft University of Technology, 2629 HS Delft, Netherlands</footnote>

<footnote n="4"><sup>4</sup>Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="5"><sup>5</sup>BIG (Badan Informasi Geospasial / Geospatial Information Agency of Indonesia), Java, 16911, Indonesia</footnote>

<footnote n="6"><sup>6</sup>Department of Civil and Environmental Engineering, University of Illinois at Urbana Champaign, Urbana, IL 61801, USA</footnote>

<footnote n="7"><sup>7</sup>Department of Geodesy and Geomatics Engineering, Institute of Technology Bandung, Kota Bandung, Jawa Barat 40132, Indonesia</footnote>

<footnote>*Now at: Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de Mécanique Paris-Saclay, Gif-sur-Yvette, France.</footnote>
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ègue et al. 2013). The 1979 M<sub>w</sub> 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<sub>w</sub> 7.4 1999 Izmit in Turkey (Bouchon et al. 2001), the M<sub>w</sub> 7.8 2001 Kunlun (Robinson et al. 2006) and the M<sub>w</sub> 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 combination 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 geometric 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 suggest 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úlveda 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 generation process of the supershear strike-slip Palu earthquake in order to more fully understand the role played by the corresponding rupture dynamics on the observed timing and first motions of the subsequent 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 earthquakes: 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 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 arrival, first motions and phases inferred from CCTV camera records near the PANT station (Carvajal et al. 2019; Sepúlveda et al. 2020).

This manuscript is organized as follows. Section 2 describes the overall methods and data used in this study, including earthquake 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 (numerically 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 (Section 2.3). The results and discussion of Section 3 provide the aforementioned 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 corresponding 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 simulations 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⁻¹. 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 supershear 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 approximation of the Euler equations, which can be written as a system of coupled hyperbolic partial differential equations given by

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

Here, $u(y, t)$ is the fluid velocity, $\eta(y, t)$ is the sea surface height and $H(y, t) = \eta(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 $tph_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 $\partial h(y, t)/\partial 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, \dots, N-1$), Fourier continuation algorithms append
a small number of points to the discretized 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) = \sum_{k=-M}^{M} a_k e^{\frac{2\pi ik y}{1+d}}, \quad u_{\text{cont}}(y) = \sum_{k=-M}^{M} b_k e^{\frac{2\pi ik y}{1+d}} \quad (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, \dots, N - 1$. Spatial derivatives for the shallow water system given by (1) are then computed by exact term-wise differentiation of (2) as

$$ \begin{aligned} \frac{\partial \eta}{\partial y}(y_i) &= \frac{\partial \eta_{\text{cont}}}{\partial y}(y_i) = \sum_{k=-M}^{M} \left( \frac{2\pi ik}{1+d} \right) a_k e^{\frac{2\pi ik y_i}{1+d}}, \\ \frac{\partial u}{\partial y}(y_i) &= \frac{\partial u_{\text{cont}}}{\partial y}(y_i) = \sum_{k=-M}^{M} \left( \frac{2\pi ik}{1+d} \right) b_k e^{\frac{2\pi ik y_i}{1+d}}. \end{aligned} \quad (3) $$

In essence, FC algorithms add a (fixed) 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 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 presentation 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 completed 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 equations 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 constraints on the temporal discretization—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 partial differential equation solvers have been successfully constructed and utilized for a variety of physical problems including those governed by radiative 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 equations (Amlani et al. 2021), Navier-Cauchy elastodynamics equations (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 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 (Rebischung & Schmid 2016) are obtained from the Jet Propulsion Laboratory (JPL). A satellite elevation mask angle of 7° and absolute International 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 ambiguities (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 corrected to UTC time by subtracting 18 s. The coseismic displacement of the station simply follows from epoch-to-epoch coordinate differences. The standard available script has been modified to properly weigh the phase/code measurements of the stations and also to output the correlations. The X-Y-Z Cartesian component positions are then converted to the north-east-up positions along with their formal 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 linear 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 $\sqrt{2c_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⁻¹ versus ~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{2c_s}$ and the P-wave speed, $c_p$, of the medium (the absolute limiting speed of the rupture). Socquet et al. (2019) and Bao 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;
<figure><img src="image_2.png" /><figcaption>**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 the 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.</figcaption></figure>

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 supershear earthquake simulation (Section 2.1) whose rupture propagates at a speed of $v = 1.6c_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 1b 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.8c_s$) finds that they are in poor agreement with GPS data (Figs 1d and e). This clearly suggests that the supershear rupture speed was $1.6c_s$ (around $5.3 \text{ km s}^{-1}$) when it passed by PALP (Ulrich *et al. 2019*, also find a speed greater than $\sqrt{2}c_s$). We have thus provided the definitive first near-field high-rate GPS-based proof that the earthquake rupture actually did go supershear as claimed and, further, have validated 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 = \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. 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 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. 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 ~1 m s⁻¹ 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 analogous 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, permanent ground displacement at the end of the profiles in Fig. 3 and is expectedly on the order of a few centimetres.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Diagram of the non-linear shallow water wave system for tsunami height $\eta(y, t)$, initial bathymetry $h_0(y)$ (basin width 9.2 km, maximum depth 710 m) and bathymetry perturbation $h(y, t)$ (i.e. the source).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Snapshot of the dynamic vertical velocity (top panel) from a supershear earthquake and the temporal evolution (bottom panel) of both velocity and displacement at an example point $(x_0, y_0)$ (indicated by the light green circle). The dark green dots correspond to the source locations used to perturb the bathymetry domain illustrated in Fig. 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The static displacement field (top panel) due to a supershear earthquake and the corresponding spatial profile in $y$ with its temporal evolution (bottom panel) at an example point $(x_0, y_0)$ (indicated by the light green circle). The dark green dots correspond to the source locations used to perturb the bathymetry domain in Fig. 2.</figcaption></figure>

Using such inputs with the FC-based tsunami model described in Section 2.2, Figs 5(a) and (b) present the corresponding results of the simulated water height $z = \eta^*(y, t)$, normalized by the absolute maximum from the dynamic case (i.e. $\eta^*(y, t) = \eta(y, t)/\max_i|\eta_{\text{dynamic}}(y, t)|$), at various synthetic stations (whose locations are indicated in Fig. 5c) simulated by both the dynamic and static (classical) sources generated from the same supershear earthquake simulation. Fig. 6

additionally presents the complete spatiotemporal evolution. The numerical modelling has been conducted at a much higher temporal resolution (a time-step of $\Delta t = 2.62 \times 10^{-3}$ s) but plotted at 10 Hz. The effects of the dynamic source, which are on the order of seconds, clearly produce high-frequency and high-amplitude waves in contrast with the static source (see Fig. 7 for a comparison of the spectral content between the two). These high-energy waves are
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Simulated tsunamis generated by dynamic and static (classical) sources. (a, b) The time histories (sampled at 10 Hz) of normalized water heights $z = \eta^*(y, t)$ predicted at various synthetic stations [located at L1–L4 and R1–R3 in (c)] along the Palu bay for tsunamis generated by a supershear earthquake due to (a) dynamic and (b) static sources. The dashed line in (b) corresponds to the duration of the earthquake. (c) The computational domain overlaid with the locations of the synthetic stations L1–L4 and R1–R3.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** The complete solution, over the first 5 min, of the normalized water height $\eta^*(y, t)$ due to (left-hand panel) a static source and (right-hand panel) a dynamic source (both generated from the same supershear earthquake).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Magnitude scalogram of the spectral contributions at the synthetic tidal gauge station R3 (whose location is indicated in Fig. 5) due to static (left-hand panel) and dynamic (right-hand panel) sources.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Simulated snapshots at various times of normalized tsunami waves along the entire Palu bay generated by static source and dynamic source models. Here, $\eta^*$ is the tsunami height normalized by the absolute maximum from the dynamic case, i.e. $\eta^*(y, t) = \eta(y, t)/\max_t|\eta_{\text{dynamic}}(y, t)|$.</figcaption></figure>

generated earlier than those of the static case but start shedding their high energy content as they slow down in their progress towards the coastline; the two begin to resemble one another in shape (Fig. 8 presents an alternate visualization in the form of snapshots-in-time across the bay). We note that, for comparison throughout, we have presented normalized water heights: since more energy of the Mach fronts is carried *along* the fault (Bernard & Baumont 2005; Dunham & Bhat 2008) running in the direction *x* (Fig. 3), the 1-D model in *y* will naturally generate lower amplitudes (on the order of half a metre). However, similar tsunami signatures can still be expected and, indeed, Elbanna et al. (2021) have demonstrated that, by incorporating horizontal motions in generic 2-D/3-D bay-like bathymetry, similar behaviour can be observed but with amplitudes on the order of metres.

Although the final waveforms are similar, a notable feature of Fig. 8 is the earlier arrival at the coastline for the dynamic case. This is more clearly illustrated in Fig. 9(a), which presents the corresponding simulated time histories at the PANT station (whose geographic location is indicated in Fig. 1a) and, more importantly, presents a comparison between the waveforms of these models with those generated at 1 Hz by carefully calibrated, and timed, CCTV and other video sources in the vicinity of the PANT tidal gauge (Carvajal et al. 2019; Sepúlveda et al. 2020). The simulations and camera records indicate an arrival that is expectedly (Sepúlveda et al. 2020) not captured by the coarse (0.02 Hz) tidal gauge at PANT, whose corresponding comparison to simulations is presented in Fig. 9(b). Remarkably, the arrival and first motions observed from the camera records in Fig. 9(a) are in excellent agreement with the 1-D approximation generated by excitation from the dynamic source. 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 localized effects of the 2-D/3-D bathymetric
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparisons of model predictions and observations at the PANT station. (a) The time histories of normalized water heights predicted by simulations and those observed by the high-resolution (1 Hz) PANT video record waveforms obtained from the author data provided in (Carvajal *et al.* 2019). Here, $\hat{\eta} \equiv \eta^* = \eta/\max_t|\eta_{\text{dynamic}}|$ for the simulations and $\hat{\eta} = \eta_{\text{CCTV}}/\max_t|\eta_{\text{CCTV}}|$ for the CCTV video-generated waveforms. (b) Corresponding normalized comparison with the low-resolution (0.02 Hz) tidal gauge.</figcaption></figure>

profile. Nevertheless, the tsunami arrival and primary dynamics are correctly reproduced.

By contrast, the static source model predicts a much later arrival. However, this is to be expected since we have used a common approach of the classical modelling community where static results are shifted by the earthquake duration [i.e. the time taken to establish the final vertical displacement that is used for the static source, or about 42 s at Palu (USGS 2018)]. Some static models do account for the finite duration of the rupture by gradually increasing the static offset over this duration, but this does not account for the full wave dynamics of the source (Satake *et al.* 2013). This is a reasonable assumption for far-field tsunamis, but it is not clear that this is justifiable for a near-field source like at Palu bay, nor is it clear how much of a shift should be introduced (Lotto *et al.* 2017). Indeed, for a fair comparison between static and dynamic models, one should wait until the end of the rupture to obtain the final static offset because of secondary slip pulses and various reflections from the surrounding medium. This can be seen in the GPS records in Fig. 1, where non-negligible ground motion is still being recorded after the passage of the main rupture pulse. We also note that the correct timing prediction is only possible through simulations informed by the full supershear dynamics (which need not make any such assumptions), where the corresponding comparisons in Fig. 9 suggest an essentially Occam’s razor explanation for the arrival observed by the PANT video waveforms: when the Palu rupture went supershear, the high-frequency ground velocities carried by the shock fronts initiated a tsunami in Palu bay at the instance when the rupture swept past the station at $t \approx 13$ s (see also Figs 6 and 8).

## 4 CONCLUSIONS

Hence we confirm that the Palu, Indonesia earthquake went supershear (via the first direct observation of such a rupture in a GPS station, accomplished here using the unique near-source signatures

of supershear (Dunham & Archuleta 2005; Mello *et al.* 2014)) and conclude that, by modelling the effects of supershear on the generation of tsunamis in a shallow geometry, the corresponding ground motion resulting from the associated Mach fronts (which carry minimally attenuated velocities to large distances) may well have contributed to the initiation of the Palu tsunami. This work provides a robust proof-of-concept, albeit in 1-D, on the contribution of shock fronts in tsunami generation. In order to gain further insight into this process, more detailed modelling in 3-D is needed to account for, e.g. geometrical spreading, attenuation and detailed 3-D velocity structure from tomographic studies. Regardless, since nothing geologically specific about the bay, except its geometry, 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. Well-known strike-slip faults that cut through various gulfs and bays (Robinson *et al.* 2010) include:

(i) Tomales bay in California, which is crossed offshore by the San Andreas fault system (Johnson & Beeson 2019).

(ii) Izmit bay in Turkey, which is crossed by the North Anatolian fault (Altinok *et al.* 2001).

(iii) Gulf of Tonkin in Vietnam, which is intersected by the Red river fault system (Tapponnier *et al.* 1986).

(iv) Gulf of Martaban in Burma, which is cut by the Sagaing fault (Vigny 2003).

(v) Gulf of Aqaba in the northern tip of the Red Sea, which is crossed by the Dead Sea fault system (Ambraseys *et al.* 1994).

(vi) Several bays and straits in the Philippines that are cut through by the Luzon fault system (Yumul *et al.* 2003).

Some of 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 su-
pershear 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
re-emphasize that both the supershear rupture *and* a shallow bay
are key to generate contributions to tsunami motions. We thus sug-
gest 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 mecha-
nism, the depth *and* the speed of the rupture all contribute towards
the generation of tsunamis.

ACKNOWLEDGMENTS

FA would like to thank N.M. Pahlevan at USC for encouraging this
work. HSB acknowledges the European Research Council grant
PERSISMO (grant 865411) for partial support of this work. AS
would like to acknowledge the European Research Council grant
REALISM (2016-grant 681346). AJR was supported by the Cal-
tech/MCE Big Ideas Fund (BIF). AE acknowledges support by the
National Science Foundation (CAREER Award Number 1753249).
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). GPS data acquisition and research was 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 for hosting the
GPS station and for being available 24/7 to assist us with optimal
operation of the equipment. A special thanks to our local survey
staff, B.R. Umar and A. Urif, for their continued support (including
directly after the earthquake) in inspecting the GPS stations.

DATA AVAILABILITY

All codes are available upon reasonable request to the corresponding
author. All relevant data relating to the tsunami modelling and the
PALP GPS are available on Zenodo: https://doi.org/10.5281/zeno
do.5018003

REFERENCES

Albin, N. & Bruno, O.P., 2011. A spectral FC solver for the compressible
Navier-Stokes equations in general domains. I: explicit time-stepping, *J*
*Comput. Phys.*, **230**(16), 6248–6270.

Altinok, Y., Tinti, S., Alpar, B., Yalciner, A., Ersoy, Ş., Bortolucci, E. & Armigliato, A., 2001. The tsunami of August 17, 1999 in Izmit bay, *Turkey, Natural Hazards*, **24**(2), 133–146.

Ambraseys, N.N., Melville, C.P. & Adams, R.D., 1994. *The Seismicity of Egypt, Arabia and the Red Sea*, Cambridge Univ. Press.

Amlani, F. & Bruno, O.P., 2016. An FC-based spectral solver for elastody-
namic problems in general three-dimensional domains, *J. Comput. Phys.*,
**307**, 333–354.

Amlani, F., Bruno, O.P., López-Vázquez, J.C., Trillo, C., Doval, Á.F., Fernández, J.L. & Rodríguez-Gómez, P., 2019. Transient propagation and scattering of quasi-Rayleigh waves in plates: quantitative comparison between pulsed TV-holography measurements and FC (Gram) elastody-namic simulations, preprint (arXiv:1905.05289).

Amlani, F. & Pahlevan, N.M., 2020. A stable high-order FC-based method-
ology for hemodynamic wave propagation, *J. Comput. Phys.*, **405**,
doi:10.1016/j.jcp.2019.109130.

Amlani, F., Wei, H. & Pahlevan, N.M., 2021. A new pseudo-spectral method-
ology without numerical diffusion for conducting dye simulations and
particle residence time calculations, preprint (arXiv:2112.05257).

Andrews, D., 1976. Rupture velocity of plane strain shear cracks, *J. Geophys. Res.*, **81**(32), 5679–5687.

Archuleta, R.J., 1984. A faulting model for the 1979 Imperial Valley earthquake, *J. Geophys. Res.*, **89**(B6), 4559–4585.

ASEAN, 2018. Situation Update No. 15 - FINAL: M 7.4 Earthquake and Tsunami, Sulawesi, Indonesia.

Bao, H., Ampuero, J.-P., Meng, L., Fielding, E.J., Liang, C., Milliner, C.W., Feng, T. & Huang, H., 2019. Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu earthquake, *Nat. Geosci.*, **12**(3), 200–205.

Bernard, P. & Baumont, D., 2005. Shear mach wave characterization for kinematic fault rupture models with constant supershear rupture velocity, *Geophys. J. Int.*, **162**, 431–447.

Bertiger, W., Desai, S.D., Haines, B., Harvey, N., Moore, A.W., Owen, S. & Weiss, J.P., 2010. Single receiver phase ambiguity resolution with GPS data, *J. Geod.*, **84**(5), 327–337.

Bouchon, M., Bouin, M.-P., Karabulut, H., Toksöz, M.N., Dietrich, M. & Rosakis, A.J., 2001. How fast is rupture during an earthquake? New insights from the 1999 Turkey earthquakes, *Geophys. Res. Lett.*, **28**(14), 2723–2726.

Bruno, O.P., Cubillos, M. & Jimenez, E., 2019. Higher-order implicit-
explicit multi-domain compressible Navier-Stokes solvers, *J. Comput.
Phys.*, **391**, 322–346.

Bruno, O.P. & Prieto, A., 2014. Spatially dispersionless, unconditionally stable FC- AD solvers for variable-coefficient PDEs, *J. Sci. Comput.*, **58**(2), 331–366.

Bryant, E., 2008. *Tsunami: The Underrated Hazard*, Springer-Verlag.

Burridge, R., 1973. Admissible speeds for plane-strain self-similar shear cracks with friction but lacking cohesion, *Geophys. J. Int.*, **35**(4), 439–455.

Carvajal, M., Araya-Cornejo, C., Sepúlveda, I., Melnick, D. & Haase, J.S., 2019. Nearly instantaneous tsunamis following the Mw 7.5 2018 Palu earthquake, *Geophys. Res. Lett.*, **46**(10), 5117–5126.

Dalguer, L.A. & Day, S.M., 2007. Staggered-grid split-node method for spontaneous rupture simulation, *J. Geophys. Res.*, **112**(B2), doi:10.1029/2006JB004467.

Das, S. & Aki, K., 1977. A numerical study of two-dimensional spontaneous rupture propagation, *Geophys. J. Int.*, **50**(3), 643–668.

Dunham, E.M. & Archuleta, R.J., 2005. Near-source ground motion from steady state dynamic rupture pulses, *Geophys. Res. Lett.*, **32**(3), doi:10.1029/2004GL021793.

Dunham, E.M. & Bhat, H.S., 2008. Attenuation of radiated ground motion and stresses from three-dimensional supershear ruptures, *J. Geophys. Res.*, **113**(B8), doi:10.1029/2007JB005182.

Dutykh, D. & Clamond, D., 2016. Modified shallow water equations for significantly varying seabeds, *Appl. Math. Modell.*, **40**(23–24), 9767–9787.

Elbanna, A., Abdelmeguid, M., Ma, X., Amlani, F., Bhat, H.S., Synolakis, C. & Rosakis, A.J., 2021. Anatomy of strike-slip fault tsunami genesis, *Proc. Natl. Acad. Sci.*, **118**(19), doi:10.1073/pnas.2025632118.

Ellsworth, W. et al., 2004. Near-field ground motion of the 2002 Denali fault, Alaska, earthquake recorded at pump station 10, *Earthq. Spectra*, **20**(3), 597–615.

Favreau, P., Campillo, M. & Ionescu, I.R., 2002. Initiation of shear instability in three-dimensional elastodynamics, *J. Geophys. Res.*, **107**(B7), ESE-4.

Fontana, M., Bruno, O.P., Mininni, P.D. & Dmitruk, P., 2020. Fourier con-
tinuation method for incompressible fluids with boundaries, *Comp. Phys.
Commun.*, **256**, doi:10.1016/j.cpc.2020.107482.

Fritz, H.M. et al., 2018. Field survey of the 28 September 2018 Sulawesi tsunami, in *AGU Fall Meeting Abstracts*, Vol. **2018**, pp. NH22B-04.

Gaggioli, E.L., Bruno, O.P. & Mitnik, D.M., 2019. Light transport with the equation of radiative transfer: the Fourier continuation–discrete ordinates (FC–DOM) method, *J. Quant. Spectrosc. Radiat. Transf.*, **236**, 106589.
He, L., Feng, G., Li, Z., Feng, Z., Gao, H. & Wu, X., 2019. Source parameters and slip distribution of the 2018 Mw 7.5 Palu, Indonesia earthquake estimated from space-based geodesy, *Tectonophysics*, **772**, doi:10.1016/j.tecto.2019.228216.

Heidarzadeh, M., Muhari, A. & Wijanarto, A.B., 2019. Insights on the source of the 28 September 2018 Sulawesi tsunami, Indonesia based on spectral analyses and numerical simulations, *Pure Appl. Geophys.*, **176**(1), 25–43.

Jamelot, A., Gailler, A., Heinrich, P., Vallage, A. & Champenois, J., 2019. 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**(8), 3351–3376.

Johnson, S.Y. & Beeson, J.W., 2019. Shallow structure and geomorphology along the offshore Northern San Andreas Fault, Tomales point to Fort Ross, California, *Bull. Seism. Soc. Am.*, **109**(3), 833–854.

Kajiura, K., 1963. The leading wave of a tsunami, *Bull. Earthq. Res. Inst., Univ. Tokyo*, **41**(3), 535–571.

Lotto, G.C., Nava, G. & Dunham, E.M., 2017. Should tsunami simulations include a nonzero initial horizontal velocity?, *Earth, Planets Space*, **69**(1), 1–14.

Lyon, M. & Bruno, O.P., 2010. High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations, *J. Comput. Phys.*, **229**(9), 3358–3381.

Mai, P.M., 2019. Supershear tsunami disaster, *Nat. Geosci.*, **12**(3), 150–151.

Mello, M., Bhat, H., Rosakis, A. & Kanamori, H., 2014. Reproducing the supershear portion of the 2002 Denali earthquake rupture in laboratory, *Earth Planet. Sci. Lett.*, **387**, 89–96.

Muhari, A., Imamura, F., Arikawa, T., Hakim, A.R. & Afriyanto, B., 2018. 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** (Scientific Communication),.

Oral, E., Weng, H. & Ampuero, J.P., 2020. Does a damaged-fault zone mitigate the near-field impact of supershear earthquakes? Application to the 2018 7.5 Palu, Indonesia, earthquake, *Geophys. Res. Lett.*, **47**(1), e2019GL085649.

Passelègue, F.X., Schubnel, A., Nielsen, S., Bhat, H.S. & Madariaga, R., 2013. From sub-Rayleigh to supershear ruptures during stick-slip experiments on crustal rocks, *Science*, **340**(6137), 1208–1211.

Pedlosky, J., 2013. *Geophysical Fluid Dynamics*, Vol. **710**, Springer-Verlag.

Pugh, D. & Woodworth, P., 2014. *Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes*, Cambridge Univ. Press.

re3data.org: VMF Data Server; editing status 2021-08-24; re3data.org - Registry of Research Data Repositories. http://doi.org/10.17616/R3RD2H last accessed: 2022-05-05.

Rebischung, P. & Schmid, R., 2016. IGS14/igs14. atx: a new framework for the IGS products, in *Proceedings of the AGU Fall Meeting 2016*, San Francisco, CA.

Robinson, D., Brough, C. & Das, S., 2006. The Mw 7.8, 2001 Kunlunshan earthquake: extreme rupture speed variability and effect of fault geometry, *J. Geophys. Res.*, **111**(B8), doi:10.1029/2005JB004137.

Robinson, D., Das, S. & Searle, M., 2010. Earthquake fault superhighways, *Tectonophysics*, **493**(3–4), 236–243.

Röbke, B. & Vött, A., 2017. The tsunami phenomenon, *Prog. Oceanogr.*, **159**, 296–322.

Rosakis, A., Samudrala, O. & Coker, D., 1999. Cracks faster than the shear wave speed, *Science*, **284**(5418), 1337–1340.

Sassa, S. & Takagawa, T., 2019. Liquefied gravity flow-induced tsunami: first evidence and comparison from the 2018 Indonesia Sulawesi earthquake and tsunami disasters, *Landslides*, **16**(1), 195–200.

Satake, K., Fujii, Y., Harada, T. & Namegaya, Y., 2013. Time and space distribution of coseismic slip of the 2011 Tohoku earthquake as inferred from tsunami waveform data, *Bull. Seism. Soc. Am.*, **103**(2B), 1473–1492.

Sepúlveda, I., Haase, J.S., Carvajal, M., Xu, X. & Liu, P.L., 2020. Modeling the sources of the 2018 Palu, Indonesia, tsunami using videos from social media, *J. Geophys. Res.*, **125**(3), doi:10.1029/2019JB018675.

Shahbazi, K., Hesthaven, J.S. & Zhu, X., 2013. Multi-dimensional hybrid Fourier continuation–WENO solvers for conservation laws, *J. Comput. Phys.*, **253**, 209–225.

Socquet, A., Hollingsworth, J., Pathier, E. & Bouchon, M., 2019. Evidence of supershear during the 2018 magnitude 7.5 Palu earthquake from space geodesy, *Nat. Geosci.*, **12**(3), 192–199.

Synolakis, C.E. & Okal, E.A., 2005. 1992–2002: perspective on a decade of post-tsunami surveys, in *Tsunamis: Case Studies and Recent Developments*, pp. 1–29, Springer.

Tanioka, Y. & Satake, K., 1996. Tsunami generation by horizontal displacement of ocean bottom, *Geophys. Res. Lett.*, **23**(8), 861–864.

Tapponnier, P., Peltzer, G. & Armijo, R., 1986. On the mechanics of the collision between India and Asia, *Geol. Soc., Lond., Spec. Publ.*, **19**(1), 113–157.

Ulrich, T. et al., 2019. Coupled, physics-based modeling reveals earthquake displacements are critical to the 2018 Palu, Sulawesi tsunami, *Pure Appl. Geophys.*, **176**(10), 4069–4109.

Umar, M., Margaglio, G., Fitrayansyah, A., et al., 2019. Post-tsunami survey of the 28 September 2018 tsunami near Palu Bay in Central Sulawesi, Indonesia: impacts and challenges to coastal communities, *Int. J. Disaster Risk Reduct.*, **38**.

USGS, 2018. M 7.5-72 km N of Palu, Indonesia.

Vigny, C., 2003. Present-day crustal deformation around Sagaing Fault, Myanmar, *J. Geophys. Res.*, **108**(B11), doi:10.1029/2002JB001999.

Weatherall, P. et al., 2015. A new digital bathymetric model of the world's oceans, *Earth Space Sci.*, **2**(8), 331–345.

Webb, F.H. & Zumberge, J. F., 1997. *An Introduction to GIPsy/oasIs-II*, Jet Propulsion Laboratory, Rep. JPL D-11088.

Wu, F.T., Thomson, K. & Kuenzler, H., 1972. Stick-slip propagation velocity and seismic source mechanism, *Bull. Seism. Soc. Am.*, **62**(6), 1621–1628.

Xia, K., Rosakis, A.J. & Kanamori, H., 2004. Laboratory earthquakes: the sub-Rayleigh-to-supershear rupture transition, *Science*, **303**(5665), 1859–1861.

Yumul, G.P., Dimalanta, C.B., Tamayo, R.A. & Maury, R.C., 2003. Collision, subduction and accretion events in the Philippines: a synthesis, *Island Arc*, **12**(2), 77–91.

Zumberge, J., Heflin, M., Jefferson, D., Watkins, M. & Webb, F., 1997. Precise point positioning for the efficient and robust analysis of GPS data from large networks, *J. Geophys. Res.*, **102**(B3), 5005–5017.