Coupled, Physics-Based Modeling Reveals Earthquake Displacements are Critical to the 2018 Palu, Sulawesi Tsunami

The September 2018, Mw\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_w$$\end{document} 7.5 Sulawesi earthquake occurring on the Palu-Koro strike-slip fault system was followed by an unexpected localized tsunami. We show that direct earthquake-induced uplift and subsidence could have sourced the observed tsunami within Palu Bay. To this end, we use a physics-based, coupled earthquake–tsunami modeling framework tightly constrained by observations. The model combines rupture dynamics, seismic wave propagation, tsunami propagation and inundation. The earthquake scenario, featuring sustained supershear rupture propagation, matches key observed earthquake characteristics, including the moment magnitude, rupture duration, fault plane solution, teleseismic waveforms and inferred horizontal ground displacements. The remote stress regime reflecting regional transtension applied in the model produces a combination of up to 6 m left-lateral slip and up to 2 m normal slip on the straight fault segment dipping 65∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$65^{\circ }$$\end{document} East beneath Palu Bay. The time-dependent, 3D seafloor displacements are translated into bathymetry perturbations with a mean vertical offset of 1.5 m across the submarine fault segment. This sources a tsunami with wave amplitudes and periods that match those measured at the Pantoloan wave gauge and inundation that reproduces observations from field surveys. We conclude that a source related to earthquake displacements is probable and that landsliding may not have been the primary source of the tsunami. These results have important implications for submarine strike-slip fault systems worldwide. Physics-based modeling offers rapid response specifically in tectonic settings that are currently underrepresented in operational tsunami hazard assessment.

and periods that match those measured at the Pantoloan 23 wave gauge and inundation that reproduces observations 24 from field surveys. We conclude that a source related 25 to earthquake displacements is probable and that land-26 sliding may not have been the primary source of the 27 tsunami. Our results have important implications for 28 submarine strike-slip fault systems worldwide. Physics-29 based modeling offers rapid response specifically in tec-30 tonic settings which are currently underrepresented in 31 operational tsunami hazard assessment.

35
Tsunamis occur due to abrupt perturbations to the wa-36 ter column, usually caused by the seafloor deforming 37 during earthquakes or submarine landslides. Devastating 38 tsunamis associated with submarine strike-slip earth-39 quakes are rare. While such events may trigger landslides 40 that in turn trigger tsunamis, the associated ground dis-41 placements are predominantly horizontal, not vertical, 42 which does not favor tsunami genesis. 43 However, strike-slip fault systems in complex tectonic 44 regions, such as the Palu-Koro fault zone cutting across 45 the island of Sulawesi, may produce vertical deformation. 46 Strike-slip systems may also include complicated fault 47 geometries, such as non-vertical faults, bends or en 48 echelon step-over structures. These can host complex 49 rupture dynamics and produce a variety of displacement 50 patterns when ruptured, which may promote tsunami 51 generation (Legg and Borrero, 2001;Borrero et al, 2004). 52 To mitigate the commonly under-represented hazard 53 of strike-slip induced tsunamis, it is crucial to funda-54 mentally understand the direct effect of coseismic dis-55 quake rupture model to a tsunami model, but these are 109 restricted to using the final, static seafloor displacement 110 field as the tsunami source. 111 To capture the physics of the interaction of the Palu 112 earthquake and tsunami we utilize a physics-based, cou-113 pled earthquake-tsunami model. The dynamic earth-114 quake rupture model incorporates spatial variation in 115 subsurface material properties, spontaneously develop-116 ing slip on a complex, non-planar system of 3D faults, off-117 fault plastic deformation, and the non-linear interaction 118 of frictional failure with seismic waves. The coseismic de-119 formation of the crust generates time-dependent seafloor 120 displacements, which we translate into bathymetry per-121 turbations to source the tsunami. The tsunami model 122 solves for non-linear wave propagation and inundation 123 at the coast.
The Indonesian island of Sulawesi is located at the 136 triple junction between the Sunda plate, the Australian 137 plate and the Philippine Sea plate (Bellier et al, 2006;138 Socquet et al, 2006138 Socquet et al, , 2019 (Fig. 1a). Convergence of 139 the Philippine and Australian plates toward the Sunda 140 plate is accommodated by subduction and rotation of 141 the Molucca Sea, Banda Sea and Timor plates, leading 142 to complicated patterns of faulting (Fig. 1a). 143 In central Sulawesi, the NNW-striking Palu-Koro 144 fault (PKF) and the WNW-striking Matano faults (MF) 145 (Fig. 1a) comprise the Central Sulawesi Fault System. 146 The Palu-Koro fault runs off-shore to the north of Su-147 lawesi through the narrow Palu Bay and is the fault 148 that hosted the earthquake that occurred on 28 Septem-149 ber 2018. With a relatively high slip rate of 40 mm/yr 150 inferred from recent geodetic measurements (Socquet 151 et al, 2006;Walpersdorf et al, 1998) and clear evidence 152 for Quaternary activity (Watkinson and Hall, 2017), the 153 Palu-Koro fault was presumed to pose a threat to the 154 region (Watkinson and Hall, 2017). In addition, four 155 tsunamis associated with earthquakes on the Palu-Koro 156 fault have struck the northwest coast of Sulawesi in the 157 The Palu earthquake triggered a local but powerful 205 tsunami that devastated the coastal area of the Palu 206 Bay quickly after the earthquake. Inundation depths of 207 over 6 m and run-up heights of over 9 m were recorded at 208 specific locations (e.g. Yalciner et al, 2018). At the only 209 tide gauge with available data, located at Pantoloan 210 harbor, a trough-to-peak wave amplitude of almost 4 m 211 was recorded just five minutes after the rupture (Muhari 212 et al, 2018). In Ngapa (Wani), on the northeastern shore 213 of Palu Bay, CCTV coverage show the arrival of the 214 tsunami wave after only 3 minutes.

215
Coseismic subsidence and uplift, as well as submarine 216 and coastal landsliding, have been suggested as causes of 217 the tsunami in Palu Bay (Heidarzadeh et al, 2018). Both 218 displacements and landsliding are documented on land 219 (Valkaniotis et al, 2018;Løvholt et al, 2018;Sassa and 220 Takagawa, 2019), and also at coastal slopes (Yalciner 221 et al, 2018 The computational model to generate the tsunami 308 scenario is StormFlash2D, which solves the nonlinear 309 shallow water equations using an explicit Runge-Kutta 310 discontinuous Galerkin discretization combined with a 311 sophisticated wetting and drying treatment for the inun-312 dation at the coast (Vater and Behrens, 2014;Vater et al, 313 2015Vater et al, 313 , 2017. A tsunami is triggered by a (possibly time-314 dependent) perturbation of the discrete bathymetry. 315 StormFlash2D allows for stable and accurate simula-316 tion of large-scale wave propagation in deep sea, as 317 well as small-scale wave shoaling and inundation at the 318 shore, thanks to a multi-resolution adaptive mesh refine-319 ment approach based on a triangular refinement strategy 320 (Behrens et al, 2005;Behrens and Bader, 2009). Bottom 321 friction is parameterized through Manning friction by a 322 split-implicit discretization (Liang and Marche, 2009). 323 The model's applicability for tsunami events has been 324 validated by a number of test cases (Vater et al, 2018), 325 which are standard for the evaluation of operational 326 tsunami codes (Synolakis et al, 2007).

327
Coupling between the earthquake and tsunami mod-328 els is realized through the time-dependent coseismic 3D 329 seafloor displacement field computed in the dynamic 330 earthquake rupture scenario, which is translated into 331 2D bathymetry perturbations of the tsunami model. The 3D dynamic rupture model of the Sulawesi earth-334 quake requires initial assumptions related to the struc-335 ture of the Earth, the structure of the fault system, the 336 stress state, and the frictional strength of the faults. 337 These input parameters are constrained by a variety of 338 independent near-source and far-field data sets. Most 339 importantly, we aim to ensure mechanical viability by a 340 systematic approach integrating the observed regional 341 stress state and frictional parameters and including state-342 of-the-art earthquake physics and fracture mechanics 343 concepts in the model (Ulrich et al, 2019).    We constrain the 3D structure of these faults using 391 focal mechanisms and geodetic data. We assume that is consistent with the analysis of Bellier et al (2006).

403
The southern end of the Palu segment bends towards 404 the Saluki segment and features a dip of 60 • to the 405 northeast, as constrained by the source mechanism of 406 the 2005 M w 6.3 event (see Fig. 1b). In contrast, we 407 assume that the Saluki segment is vertical. The assigned 408 dip of 90 • acknowledges the inferred ground deformation 409 of comparable amplitude and extent on both sides of 410 this fault segment (De Michele, 2018). All faults reach 411 a depth of 20 km. The fault system is subject to a laterally homogeneous 414 regional stress field, systematic constraints based on 415 seismo-tectonic observations, fault fluid pressurization 416 and the Mohr-Coulomb theory of frictional failure fol-417 lowing Ulrich et al (2019). This is motivated by the 418 fact that tractions on and strength of natural faults 419 are difficult to quantify. With this approach, only four 420 parameters must be specified to fully describe the state 421 of stress and strength governing the fault system, as fur-422 ther detailed in the appendix (Sec. 8.3). This systematic 423 approach facilitates rapid modeling of an earthquake.

424
Using static considerations and few trial dynamic 425 simulations, we identify an optimal stress configuration 426 for this scenario that simultaneously (i) maximizes the 427 ratio of shear over normal stress all across the fault 428 system; (ii) determines shear traction orientations that 429 predict surface deformation compatible with the mea-430 sured ground deformation and focal mechanisms; and 431 (iii) allows dynamic rupture across the fault system's 432 geometric complexities.

433
The resulting physical model is characterized by a 434 stress regime acknowledging transtension, high fluid 435 pressure, and relatively well oriented, apparently weak 436 faults. The effective confining stress increases with depth 437 by a gradient of 5.5 MPa/km. From 11-15 km depth, 438 we taper the deviatoric stresses to zero, to represent 439 the transition from a brittle to a ductile deformation 440 regime. The depth range is consistent with the 12 km 441 interseismic locking depth estimated by Vigny et al 442 (2002).

444
Failure is initiated within a highly overstressed circular 445 patch with a radius of 1.5 km situated at a depth of 446 10 km. This depth is at the shallow end of the range 447 of inferred hypocentral depths (Valkaniotis et al, 2018) 448 and shallower than the modeled brittle-ductile transition 449 mimicking the lower end of the seismogenic zone.    Tanioka and Satake (1996). This is motivated by the potential influence of Palu 504 Bay's steep seafloor slopes (more than 50%). The ground 505 displacement of the earthquake model is translated into 506 the tsunami generating bathymetry perturbation by is the bathymetry (increasing in the 509 upward direction). ∆b is time-dependent, since ∆x, ∆y 510 and ∆z are time-dependent (cf. Fig. S2). The tsunami is 511 sourced by adding ∆b to the initial bathymetry and to-512 pography of the tsunami model. It should be noted that 513 a comparative scenario using only ∆z as bathymetry 514 perturbation (see appendix, Sec. 8.5) did not result in 515 large deviations with regards to the preferred model.

516
The domain of the computational tsunami model 517 (latitudes ranging from −1 • to 0 • , longitudes ranging 518 from 119 • to 120 • , see Fig. 2) encompasses Palu Bay 519 and its near surroundings in the Makassar Strait, since 520 we here focus on the wave behavior within the Bay of 521 Palu. The tsunami model is initialized as an ocean at 522 rest, for which (at t = 0) the initial fluid depth is set in 523 such manner that the sea surface height (ssh, deviation 524 from mean sea level) is equal to zero everywhere in 525 the model domain. Additionally, the fluid velocity is 526 set to zero. This defined initial steady state is then 527 altered by the time-dependent bathymetry perturbation 528 throughout the simulation, which triggers the tsunami. 529 The simulation is run for 40 min (simulation time), 530 which needs 13 487 time steps.

531
The triangle-based computational grid is initially re-532 fined near the coast, where the highest resolution within 533 Palu Bay is about 3 arc seconds (or 80 m). This results 534 8 in an initial mesh of 153 346 cells, which expands to 535 more than 300 000 cells during the dynamically adaptive 536 computation. The refinement strategy is based on the 537 gradient in sea surface height (ssh).

538
The parametrization of bottom friction includes the Manning's roughness coefficient n. We assume n = 0.03,     (Fig. 4a). 579 South of this small bend, the slip magnitude increases 580 and remains mostly homogeneous, ranging between 6 581 and 8 m. Peak slip occurs on the Palu segment.

582
Over most of the fault network, the faulting mech-583 anism is predominantly strike-slip, but does include a 584 small to moderate normal slip component (Fig. 4b). This 585 dip-slip component varies as a function of fault orienta-586 tion with respect to the regional stress field. It increases 587 at the junction between the Northern and Palu segment 588 just south of Palu Bay, and at the big bend between the 589 Palu and Saluki fault segments, where dip-slip reaches 590 a maximum of approx. 4 m. Pure strike-slip faulting 591 is modeled on the southern part of the vertical Saluki 592 segment (Fig. 4b). The dip-slip component along the 593 rupture shown in Fig. 4b produces subsidence above the 594 hanging wall (east of the fault traces) and uplift above 595 the foot wall (west of the fault traces). The resulting 596 seafloor displacements are further discussed in Sec. 4.2. 597

Earthquake rupture speed 598
The earthquake scenario features an early and persis-599 tent supershear rupture velocity (Fig. 4d). This means 600 that the rupture speed exceeds the seismic shear wave 601 velocity (V s ) of 2.5 to 3.1 km/s in the vicinity of the 602 fault network from the onset of the event. This agrees 603 with the inferences for supershear rupture by Bao et al 604 (2019) from back-projection analyses and by Socquet 605 et al (2019) from satellite data analyses. However, we 606 here infer supershear propagation faster than Eshelby 607 speed ( √ 2V s ), and thus faster than Bao et al (2019), well 608 within the stable supershear rupture regime (Burridge, 609 1973). The dynamic rupture scenario satisfactorily reproduces 613 the teleseismic surface waves (Fig. 5a) and body waves 614 (Fig. 5b). Synthetics are generated at 5 teleseismic sta-615 tions around the event (Fig. 5c)

653
The patterns and magnitudes of the final horizon-654 tal surface displacements in two dimensions (black ar-655 rows in Fig. 6a) are inferred from subpixel correlation 656 of coseismic optical images acquired by the Coperni-657 cus Sentinel-2 satellites by the European Space Agency 658 (ESA) (De Michele, 2018). We use both, east-west and 659 north-south components from optical image correlation. 660 We also infer coseismic surface displacements by in-661 coherent cross correlation of synthetic aperture radar 662 (SAR) images acquired by the Japan Aerospace Ex-663 ploration Agency (JAXA) Advanced Land Observation 664 Satellite-2 (ALOS-2). SAR can measure surface displace-665 ments horizontally in the along-track direction and in the 666 slant direction between the satellite and the ground that 667 is a combination of vertical and horizontal displacement. 668 Here, we use the along-track horizontal displacements 669 (Fig. 6c) that are nearly parallel to the strike of the fault. 670 Further details about our data processing approach and 671 the dataset used can be found in appendix Sec. 8.6.

672
The use of two independent but partially coinciding 673 datasets provides additional insight on data quality. We 674 compare the SAR data and the optical data by project-675 ing the optical data into the along-track direction of the 676 SAR data. This allows for identification of the robust fea-677 tures in the imaged surface displacements. Along most 678 of the rupture, fault displacements are sharp and linear, 679 highlighting smooth and straight fault orientations with 680 some bends. Both datasets appear to be consistent to 681 first order (±1m)    The tsunami generated in this scenario is mostly 737 localized in Palu Bay, which is illustrated in snapshots 738 of the dynamically adaptive tsunami simulation after 739 20 s and 600 s simulation time in Fig. 8. This is expected 740 as the modeled fault system is offshore only within 741 the Bay. At 20 s, the seafloor displacement due to the 742 earthquake is clearly visible in the sea surface height 743 (ssh) within Palu Bay. Additionally, the effect of a small 744 uplift is visible along the coast north of the Bay.The 745 local behavior within Palu Bay is displayed in Fig. 9   with observations is shown in Fig. 12. The overall agree-798 ment is quite remarkable, with some overestimation of 799 the run-up in the northern margins of the bay and some 800 slight underestimation in the southern part near Grand-801 mall Palu City. In general, errors are lower than 10%. 802 What we can conclude is that large misfit in the run-up 803 heights are more or less randomly distributed, suggest-804 ing local amplification effects that cannot be captured in 805 the scenario due to insufficient bathymetry/topography 806 resolution. Fig. 13 shows maximum inundation depths 807 computed from the tsunami scenario near Palu City. 808 Qualitatively, the results from the scenario agree quite 809 well with observations, as the largest inundation depths 810 are close to the Grandmall area, where vast damage due 811 to the tsunami was reported.

812
In summary, the tsunami scenario sourced by coseis-813 mic displacements from the dynamic earthquake rupture 814 scenario yields results that are qualitatively compara-815 ble to available observations. Wave amplitudes match 816 well, as do the run-up distribution and the inundation 817 distances, given the limited quality of the available to-818 pography data. The Palu, Sulawesi tsunami was as unexpected as it 821 was devastating. While the Palu-Koro fault system 822 was known as a very active strike-slip plate boundary 823 tsunamis from strike-slip events are generally not antic-824 ipated. Fears arise that other regions, currently not ex-825 pected to sustain tsunami-triggering ruptures, are at risk. 826 The here presented physics-based, coupled earthquake-827 tsunami model shows that a submarine strike-slip fault 828 can produce a tsunami, if a component of dip-slip fault-829 ing occurs. In the following, we discuss advantages and 830 limitations of physics-based models of tsunamigenesis 831 as well as of the earthquake and tsunami model indi-832 vidually. We then focus on the broader implications of 833 rapid coupled scenarios for seismic hazard mitigation 834 and response. Finally, we look ahead to improving the 835 here presented coupled model in light of newly available 836 information and data. We constrain the initial conditions for our coupled model 840 according to the available earthquake data and physi-841 cal constraints provided by previous studies, including 842 those reporting regional transtension (Walpersdorf et al, 843 1998;Socquet et al, 2006;Bellier et al, 2006). A stress 844 field characterized by transtension induces a normal 845 component of slip on the dipping faults in the earth-846 quake scenario. The here assumed degree of transtension 847  magnitude. (Fig. 7c). This is sufficient for triggering 856 a realistic tsunami that reproduces the observational 857 data quite well. In particular it is enough to obtain the 858 observed wave amplitude at the Pantoloan harbor wave 859 gauge and the recorded run-up heights.

860
However, we point out that transtension is not an 861 indispensable condition to generate oblique faulting in 862 such a fault network. From static considerations, we 863 indeed infer that specific alternative stress orientations 864 can equally induce a considerable dip-slip component in 865 biaxial stress regimes (Fig. S3). Finally, incorporating the effect of landslides is likely 889 to be necessary to capture local features of the tsunami 890 wave and inundation patterns. Constraining these sources 891 is very difficult without pre-and post-event high-resolution 892 bathymetric charts. Our study suggests that these sources 893 play a secondary role in explaining the overall tsunami 894 magnitude and wave patterns, since these can be gen-895 erated by strike-slip faulting with a normal slip compo-896 nent. The speed of this earthquake is of utmost interest, al-899 though it does not provide an important contribution 900 to the tsunami generation in this scenario. We review 901 our results here and note avenues for additional mod-902 eling. The initial stress state and lithology included in 903 the physical earthquake model are areas that could be 904 improved with more in-depth study and better available 905 data.

906
The dynamic earthquake model requires supershear 907 rupture velocities to produce results that agree with the 908 teleseismic data and moment rate function. This scenario 909 also provides new perspectives on the possible timing 910 and mechanism of this supershear rupture. Bao et al 911 (2019) infer an average rupture velocity of about 4 km/s 912 from back-projection. This speed corresponds to a barely 913 stable mechanical regime, which is interpreted as being 914 promoted by a damage zone around the mature Palu-915 Koro fault that formed during previous earthquakes.

916
In contrast, our earthquake scenario features an early 917 and persistent rupture velocity of 5 km/s on average, 918 close to P-wave speed. Supershear rupture speed is en-919 abled in our model by a relatively low fault strength 920 and triggered immediately at rupture onset by a highly 921 overstressed nucleation patch. Supershear transition is 922 enabled and enhanced by high background stresses (or 923 more generally, low ratios of strength excess over stress 924 drop) (Andrews, 1976). The so called transition distance, 925 the rupture propagation distance at which supershear 926 rupture starts to occur, also depends on nucleation 927 energy (Dunham, 2007;Gabriel et al, 2012Gabriel et al, , 2013. Ob-928 servational support for the existence of a highly stressed 929 nucleation region arises from the series of foreshocks 930 that occurred nearby in the days before the mainshock, 931 including a M w 6.1 on the same day of the mainshock. 932 We conducted numerical experiments reducing the 933 level of overstress within the nucleation patch, reach-934 ing a critical overstress level at which supershear is not 935 anymore triggered immediately at rupture onset. These 936  2018). However, in light of an absence of data or models 963 justifying the introduction of complexity, we here use 964 the simplest option with a laterally homogeneous stress 965 field that honors the regional scale transtension. 966 We also note that the earthquake scenario is depen-967 dent on the subsurface structure model (  Overall, the tsunami model shows good agreement with 977 available key observations. Wave amplitudes and peri-978 ods at the only available tide gauge station in the Bay 979 match well. Run-up and inundation data from our model 980 show satisfactory agreement with the observations by 981 international survey teams (Yalciner et al, 2018).

982
Apart from the above discussed earthquake model 983 limitations that may influence the tsunami characteris-984 tics, the following additional reasons may cause devia-985 tions to tsunami observations: (a) insufficiently accurate 986 bathymetry/topography data; (b) simplified coupling 987 between earthquake rupture and tsunami scenarios; (c) 988 approximation by hydrostatic shallow water wave theory. 989 In the following we will briefly discuss these topics.

990
The insufficient resolution of the bathymetry and with the fault geometry and the regional stress field, the 1045 dynamic rupture model produces mechanically consis-1046 tent ground deformation, even in submarine areas where 1047 space borne imaging techniques are blind. These seafloor 1048 displacement time-histories, which include the influence 1049 of seismic waves, in nature contribute to source the 1050 tsunami and are utilized as such in this coupled frame-1051 work. However, the earthquake-tsunami coupling is not 1052 physically seamless. For example, as noted above, seis-1053 mic waves cannot be captured using the shallow water 1054 approach, but rather require a non-hydrostatic water 1055 body (e.g. Lotto et al, 2018). However, the coupled sys-1056 tem remains mechanically consistent to the order of 1057 the typical spatio-temporal scales governing tsunami 1058 modeling. Thus, a physics-based, coupled model is well-1059 posed to shed light on the mechanisms and competing 1060 hypotheses governing earthquake-tsunami sequences as 1061 puzzling as the Sulawesi event.

1062
The use of a dynamic rupture earthquake source 1063 has distinct contributions relative to the standard finite-1064 fault inversion source approach, which is typically used 1065 in tsunami models. The latter enables close fitting of 1066 observations through the use of a large number of free 1067 parameters. Despite recent advances (e.g., Shimizu et al, 1068 2019), kinematic models typically need to pre-define 1069 fault geometries. Naive first-order finite-fault sources 1070 are automatically determined after an earthquake and 1071 this can be done quickly (e.g. by USGS or GFZ German 1072 Research Centre for Geoscience), which is a great advan-1073 tage. Models can be improved later on by including new 1074 data and more complexity. However, kinematic models 1075 are characterized by inherent non-uniqueness and do not 1076 ensure mechanical consistency of the source (e.g., Mai 1077 et al, 2016). The physics-based model also suffers from 1078 non-uniqueness, but this is reduced, since it excludes 1079 scenarios that are not mechanically viable.

1080
These advantages and the demonstrated progress 1081 potentially make physics-based, coupled earthquake-1082 tsunami modeling an important tool for seismic haz-1083 ard mitigation and rapid earthquake response. We fa-1084 cilitate rapid modeling of the earthquake scenario by 1085 systematically defining a suitable parameterization for 1086 the regional and fault-specific characteristics. We use 1087 a pre-established, efficient algorithm, based on phys-1088 ical relationships between parameters, to assign the 1089 ill-constrained stress state and strength on the fault 1090 using a few trial simulations (Ulrich et al, 2019). This 1091 limits the required input parameters to subsurface struc-1092 ture, fault structure, and four parameters governing the 1093 stress state and fault conditions. This enables rapid response in delivering physics-driven interpretations that 1095 can be integrated synergistically with established data-1096 driven efforts within the first days and weeks after an 1097 earthquake.  The fully dynamic earthquake model captures im-1141 portant features, including the timing and speed of the 1142 rupture, 3D geometric complexities of the faults, and 1143 the influence of seismic waves on the rupture propaga-1144 tion. We find that an early-onset of supershear rupture 1145 speed, sustained for the duration of the rupture across 1146 geometric complexities, is required to match a range of 1147 far-field and near-fault observations.

1148
The modelled tsunami amplitudes and wave run-ups 1149 agree with observations within the range of modeling 1150 uncertainties dominated by the available bathymetry 1151 and topography data.We conclude that the primary 1152 tsunami source may have been coseismically generated 1153 vertical displacements. However, in a holistic approach 1154 aiming to match high-frequency tsunami features, local 1155 effects such as landsliding, non-hydrostatic wave effects, 1156 and high resolution topographical features should be 1157 included.

1158
The coupling of physics-based models, as tackled 1159 within the ASCETE framework, is specifically useful 1160 to assess tsunami hazard in tectonic settings currently 1161 underrepresented in operational hazard assessment. We 1162 demonstrate that high-performance computing empow-1163 ered dynamic rupture modeling produces well-constrained 1164 studies integrating source observations and earthquake 1165 physics very quickly after an event occurs. In the future, 1166 such physics-based earthquake-tsunami response can 1167 complement both on-going hazard mitigation and the 1168 established urgent response tool set. We thank Taufiqurrahman for helping us accessing data 1171 on Indonesian websites, and for putting us in contact 1172 with Indonesian researchers. We thank Dr. T. Yudis-1173 tira for providing their crustal velocity model of Su-1174 lawesi and Dr. Andreas Fichtner for providing us a 1175 chunk of their 'Collaborative seismic earth model'. We 1176 thank Dr. Marcello de Michele for providing his in-1177 ferred ground-deformations data and for fruitful discus-1178 sions. The ALOS-2 original data are copyright JAXA 1179 and provided under JAXA RA6 PI projects P3278 and 1180 P3360. Dr. Widodo S. Pranowo provided access to very 1181 early field survey observations. Furthermore, Dr. Abdul 1182 Muhari supported this work by providing 1-minute tide 1183 gauge data for the Pantoloan tide gauge. Finally we 1184 thank the participant of the AGU special session about 1185 the Palu earthquake for interesting discussions.

1186
The work presented in this paper was enabled by 1187 the Volkswagen Foundation (project "ASCETE", grant 1188 no. 88479). these data are not yet available. In Figure S2, we provide 1633 the displacements time histories at a few of these sites. 1634 We hope future access to this data will provide further 1635 constraints to our model. and γ. SH max is the azimuth of the maximum horizontal 1645 compressive stress; ν is a stress shape ratio balancing 1646 the principal stress amplitudes; R 0 is a ratio describing 1647 the relative strength of the faults; and γ is encapsulating 1648 fluid pressure.

1661
The fault prestress ratio R 0 describes the closeness 1662 to failure of a virtual, optimally oriented plane according 1663 to Mohr-Coulomb theory (Aochi and Madariaga, 2003). 1664 On this virtual plane, the Coulomb stress is maximized. 1665 Optimally oriented planes are critically loaded when 1666 R 0 = 1. Faults are typically not optimally oriented in 1667 reality. In a dynamic rupture scenario, only a small 1668 part of the modeled faults need to reach failure in order 1669 to nucleate sustained rupture. Other parts of the fault 1670 network can break cascadingly even if well below failure 1671 before rupture. The propagating rupture front raises the 1672 local shear tractions to match fault strength locally. 1673 We assume fluid pressure γ throughout the crust is 1674 proportional to the lithostatic stress: P f = γσ c , where 1675 γ is the fluid-pressure ratio and σ c = ρgz is the litho-1676 static pressure. A fluid pressure of γ = ρ water /ρ = 0.37 1677 indicates purely hydrostatic pressure. Higher values cor-1678 respond to overpressurized stress states. Together, R 0 1679 24 and γ control the average stress drop dτ in the dynamic 1680 rupture model as: (2)

1682
The such prescribed average stress drop dτ is a critical where τ 0 and σ n are the initial shear and normal trac-   (2019). 1735 Friction drops rapidly from a steady-state, low-velocity 1736 friction coefficient, here f 0 = 0.6, to a fully weakened 1737 friction coefficient, here f w = 0.1 (see Table S1). For computing the seafloor displacement used as source 1741 for the tsunami model, we apply the method of Tanioka 1742 and Satake (1996) to additionally account for horizontal 1743 displacements, computed from the earthquake simula-1744 tion. The final states of the three components ∆x,∆y 1745 and ∆z are given in Fig. S5. Applying the approach 1746 of Tanioka and Satake by using Eq. (1) the vertical 1747 displacement translates into ∆b, which is given in Fig. 7. 1748 The difference between ∆z and ∆b locally amounts up 1749 to 0.6m as shown in Fig. S6. Although this difference is 1750 quite remarkable and compared to the overall magnitude 1751 more than 30%, it is only very local. Due to the local 1752 bathymetry of Palu bay it also not only amplifies the 1753 displacement, but also diminishes it at some locations. 1754 The local influence of the method by 1755 Satake (1996) can be seen by comparison to the re-1756 sults section. We have run a similar simulation as de-1757 scribed in the main part of the paper, but with the 1758 computed seafloor displacement ∆z as source for the 1759 tsunami model. Snapshots of this scenario in Palu Bay 1760 can be seen in Fig. S7. Compared to the original sce-1761 nario (cf. Fig. 9) only local effects are visible, especially 1762 at points along the coast. The maximum inundation at 1763 Palu city is given for this alternative scenario in Fig. S8 S3 Magnitude and rake of prestress resolved on the fault system for a range of plausible SH max values, assuming a stress shape ratio ν = 0.5 (pure-shear). For each stress state , we show the spatial distribution of the pre-stress ratio (left) and the rake angle of the shear traction (right). Here we assume R 0 = 0.7 on the optimal plane, which results in R < R 0 for all faults since these are not optimally oriented. In blue, we label the (out-of-scale) minimum rake angle on the Palu-Saluki bend.
computation which includes horizontal displacements in 1766 the source (cf. Fig. 13). This illustrates that the method 1767 by Tanioka and Satake (1996)