Probabilistic Tsunami Hazard Assessment in The Lesser Sunda Islands, Indonesia
ShofiaKarima2✉Emailshof002@brin.go.id
HamzahLatief3,5
IyanE.Mulia4
1Programme Study Earth Sciences, Faculty of Earth Science and EngineeringBandung Institute of TechnologyJalan Ganesha No. 10100190BandungIndonesia
2Research Center for Geological DisasterNational Research and Innovation Agency40293BandungIndonesia
3Environmental and Applied Oceanography Research Group, Faculty of Earth Science and TechnologyBandung Institute of Technology40132BandungIndonesia
4Hydrography Research Group, Faculty of Earth Sciences and TechnologyBandung Institute of Technology40132BandungIndonesia
5Center for Marine and Coastal DevelopmentBandung Institute of TechnologyJalan Ganesha No. 10100190BandungIndonesia
Shofia Karima1,2*, Hamzah Latief3,5, Iyan E. Mulia4
1* Programme Study Earth Sciences, Faculty of Earth Science and Engineering, Bandung Institute of Technology, Jalan Ganesha No. 10, Bandung, 100190, Indonesia.
2Research Center for Geological Disaster, National Research and Innovation Agency, Bandung 40293, Indonesia
3Environmental and Applied Oceanography Research Group, Faculty of Earth Science and Technology, Bandung Institute of Technology, Bandung 40132, Indonesia
4 Hydrography Research Group, Faculty of Earth Sciences and Technology, Bandung Institute of Technology, Bandung 40132, Indonesia
5Center for Marine and Coastal Development, Bandung Institute of Technology, Jalan Ganesha No. 10, Bandung, 100190, Indonesia.
*Corresponding E-mail: shof002@brin.go.id
Abstract
The Lesser Sunda Islands (comprising Bali, West Nusa Tenggara, and East Nusa Tenggara) are situated between three seismic sources: the subduction zone, the fold-thrust back (FTB), and the back-arc thrust fault system in the north. The combined activity of these three sources makes the region highly vulnerable to tsunami hazards. The high population density in coastal areas and the presence of critical infrastructure in this region increase the urgency of conducting a tsunami hazard assessment. This study aims to analyze the tsunami hazard using a Probabilistic Tsunami Hazard Assessment (PTHA) based on a Green’s function approach, derived from simulations of the Shallow Water Equation (SWE). A total of 10,120 tsunami scenarios were developed from 9 fault sources using the stochastic slip distribution based on the Von Karman correlation function. All faults were divided into 108 subfaults to estimate the response of the tsunami wave at 255 observation points. The locations with the highest frequency per year for waves exceeding 1 m include Jembrana and Buleleng, northern and southern Sumbawa, as well as Sikka and west Manggarai. The tsunami hazard maps indicate that for a return period (RP) of 2475 years, more than 57% of the locations experience tsunami waves exceeding 3 m, with maximum heights reaching over 11 m. For the return period of 475 years, the rate is more than 0.4%, and for the return period of 72 years, there’s no location where the tsunami wave height exceeds 3 meters. The faults from the Bali subduction zone, the Sumbawa thrust, the Flores thrust, and the Timor FTB are the dominant contributors to the tsunami hazard in 20 strategic locations.
Keywords:
Green’s function
hazard curves
hazard maps
stochastic slip method
Shallow Water Equation model
and tsunami
1 Introduction
The Lesser Sunda Islands (including Bali, West Nusa Tenggara (NTB), and East Nusa Tenggara (NTT)are situated in a tectonically complex and highly seismically active region of Southeast Asia (Hasanah et al., 2024). This archipelago is bounded by three major earthquake-generating systems: the Sunda subduction zone to the south, the Timor Fold-Thrust Belt (FTB) to the east, and the Flores back-arc thrust system to the north (Zhao et al., 2025). Interactions among these tectonic features significantly increase the likelihood of large-magnitude earthquakes and tsunamis (Newcomb & McCann, 1987). Historical records document at least 29 tsunami events since 1815, including the devastating 1992 Mw 7.9 Flores earthquake and tsunami, which claimed more than 1,700 lives (Pararas-Carayannis, 2018). The presence of shallow faulting near the coastline, particularly along the Timor FTB and back-arc thrust systems, further heightens the risk of rapidly arriving and highly localized tsunamis (Hasanah et al., 2024).
The region’s vulnerability is amplified by dense coastal populations and critical infrastructure. According to Badan Pusat Statistik (BPS) (2013), Bali has approximately 4.3 million inhabitants, NTB 5.5 million, and NTT 5.6 million. Many of these communities reside in low-lying coastal areas that function as economic and social hubs, particularly for tourism, fisheries, and maritime transport. This combination of high exposure and elevated hazard potential underscores the need for a scientifically robust, probabilistic approach to tsunami hazard assessment. Unlike deterministic models, which may overlook variability in earthquake rupture processes and tsunami propagation, a Probabilistic Tsunami Hazard Assessment (PTHA) framework allows for long-term hazard estimation while explicitly accounting for uncertainties in seismic source parameters and slip distributions.
Previous studies have provided important insights into tsunami hazards in Indonesia. Mulia et al. (2020) introduced PTHA research using historical data and numerical modeling to estimate tsunami probabilities along Japanese coastlines, offering methodological advances later adapted to Indonesian contexts. In the Lesser Sunda region, Felix et al. (2022) modeled tsunamis in Bali and NTB using historical events and deterministic scenarios of Flores back-arc thrust sources. Wibowo et al. (2021) assessed tsunami hazards from Flores back-arc sources affecting Bali and Lombok, but their analysis incorporated limited scenario variability and excluded stochastic slip distributions. Horspool et al. (2013) developed a national-scale PTHA for Indonesia—including Bali, NTB, and NTT—yet their models did not incorporate random slip variability. Similarly, Pranantyo et al. (2013) conducted a PTHA along Bali’s northern coast, focusing on the Java subduction zone, and reported relatively low tsunami heights (< 1 m) for 100- and 200-year return periods; however, their work excluded NTB and NTT and did not consider back-arc or fold-thrust sources.
To address these gaps, this study applies a PTHA framework that integrates stochastic slip modeling with the Green’s function method and numerical solutions of the Shallow Water Equation (SWE). A total of 10,120 tsunami scenarios were simulated from nine earthquake fault sources, each subdivided into 108 subfaults, generating tsunami responses at 255 offshore observation points. The stochastic slip model, based on the von Karman autocorrelation function, enables realistic spatial variation in fault dislocation, while the Green’s function method efficiently computes tsunami propagation via linear superposition. Offshore wave heights are subsequently extrapolated to the shoreline using Green’s Law to capture shoaling effects.
The stochastic slip modeling framework introduced by Mai & Beroza (2002) provided a key methodological foundation for representing the heterogeneity of earthquake rupture processes. Mulia et al. (2020) later implemented this approach in Japan, coupling stochastic slip models with Green’s function computations to enhance PTHA accuracy and efficiency. Building on these advances, the present study combines stochastic slip, the Green’s function method, and Monte Carlo simulations to generate a large ensemble of realistic tsunami scenarios for the Lesser Sunda Islands. Unlike earlier work, this study incorporates spatially variable earthquake sources from subduction, back-arc, and fold-thrust-belt systems surrounding Bali, NTB, and NTT, and provides detailed hazard characterizations across 20 strategic coastal locations. Those locations are strategic and considerable as observation points due to several reasons: the location of high activity of people, such as airports and seaports, residential areas, and tourism areas.
By employing a comprehensive probabilistic framework and state-of-the-art stochastic modeling techniques, this research aims to advance the scientific understanding of tsunami hazards in one of Indonesia’s most exposed and tectonically complex regions. The findings will inform long-term risk reduction strategies, coastal zoning, early warning systems, and infrastructure planning in the Lesser Sunda Islands.
2. Data and Method
2.1 Active faults in the Lesser Sunda Islands
The Lesser Sunda Islands are bounded by three major fault systems: the Sunda subduction zone south of Java and Bali, the Flores back-arc thrust to the north, and the Timor Fold-Thrust Belt (FTB) near Timor Island (Fig. 1; Suarbawa et al., 2021). These tectonic structures are capable of generating large earthquakes and locally destructive tsunamis (Gunawan et al., 2022; Suarbawa et al., 2021).
Fig. 1
Location of major fault systems in the Lesser Sunda Islands.
The Flores back-arc thrust extends north of Alor Island, from Flores to Sumbawa (Suarbawa et al., 2021). The Sunda subduction zone, located south of Bali and East Java, is driven by subduction of the Indian Ocean plate beneath the Sunda plate (Nugraha et al., 2014). Meanwhile, the Timor FTB results from plate convergence, in which marine sedimentary layers are compressed and folded northward, forming thrust faults and intense deformation structures (Audley-Charles, 1986; Harris et al., 2024).
In this study, each fault was discretized into subfaults (Fig. 2) with variable slip distributions, expressed in meters. Stochastic slip modeling was employed to generate a wide range of slip variability consistent with seismic activity (Goda et al., 2016; Mai & Beroza, 2002). A total of nine earthquake sources were considered, subdivided into 108 subfaults. Subfault sizes were adapted to fault type: 40 × 40 km for subduction zones, 50 × 50 km for Timor FTB, and 30 × 30 km for back-arc thrusts. Earthquake parameters for each source are summarized in Table 1 (Goda et al., 2016; Mai & Beroza, 2002).
Fig. 2
Faults and their subfault discretization in the study area (red box).
No | Source | L (km) | W (km) | Subfault size (km) | Dip (deg) | Depth (m) | Max Mw | Sample | Scenario |
|---|---|---|---|---|---|---|---|---|---|
1 | East Java Subduction | 440 | 120 | 40 | 15 | 20000 | 8.9 | 110 | 2200 |
2 | Bali Subduction | 560 | 120 | 40 | 5.56 | 18200 | 9 | 110 | 2310 |
3 | Timor FTB 1 | 100 | 50 | 50 | 45 | 10000 | 7.5 | 110 | 660 |
4 | Timor FTB 2 | 200 | 50 | 50 | 45 | 9000 | 7.5 | 110 | 660 |
5 | Timor FTB 3 | 350 | 50 | 50 | 45 | 5600 | 8 | 110 | 1210 |
6 | Bali Thrust | 180 | 30 | 30 | 49 | 33000 | 7.4 | 110 | 550 |
7 | Sumbawa Thrust | 180 | 30 | 30 | 53.5 | 32850 | 8 | 110 | 1210 |
8 | Flores Thrust | 180 | 30 | 30 | 50.8 | 33000 | 7.5 | 110 | 660 |
9 | Flores Thrust | 120 | 30 | 30 | 55.6 | 33500 | 7.5 | 110 | 660 |
Total | 10120 | ||||||||
2.2 Green’s Function Method and Shallow Water Equation Model
The maximum tsunami height at observation points was computed using the Green’s function approach, in which waveforms from the Shallow Water Equation (SWE) are linearly combined with variable slip distributions:
1
2
where
is the tsunami waveform at the observation points,
is Green’s Function (synthetic waveforms from subfaults),
is the slip distribution in a fault,
denotes subfault locations,
is the observation points, and
is time. The number of subfault columns is identified by nr, and their rows are identified by nc.
To construct the
, each subfault was treated as an independent tsunami source, with initial seafloor deformation computed using Okada’s model (Okada, 1985), assuming 1 m slip (Fig. 3). Resulting waveforms were recorded at 255 offshore observation points and 20 strategic coastal sites at − 50 m depth.. The location of both the points of interest and strategic locations is provided in the Supplementary Fig. 1.
Fig. 3
Schematic of Green’s function computation in PTHA.
In this study, we utilized the Cornell Multi-grid Coupled Tsunami (COMCOT) to solve the linear SWE. COMCOT has been used by several previous studies and has been validated against historical data (Syamsidik et al., 2015; Tri Laksono et al., 2020). The total simulation time was 4 hours with an output interval time of one minute. We considered a two-layer nested grid system with bathymetry obtained from BATNAS (https://sibatnas.big.go.id/) (Fig. 4). The first domain has a grid size of 0.5 arc minutes, and the smallest one is 0.167 arc minutes. To obtain the stability condition (CFL), the numerical time step needed was one second.
Fig. 4
Nested computational domains for COMCOT simulations..
An amplification factor using Green’s law was used to transform the tsunami waveforms at 50 m to 1 m water depth (Green, 1838).
3
where
and
are the tsunami heights at coastal and offshore points, and
and
are the corresponding water depths.
2.3 Stochastic slip model and variance analysis
The RuptGen model developed by Mai & Beroza (2002) was simulated in this study to generate various slip distributions in the 9 faults. This method applied a von Karman autocorrelation function, a spatial random field of anisotropic wavenumber spectra, which can be used to describe the complexity of seismic slip. Initially, the slip distribution value was calculated on a grid spacing with a size of 1 x 1 km and later interpolated into the subfault size with bilinear interpolation, as previously applied in Mulia et al. (2020). One of the examples is shown in Fig. 5, one of the slip distributions for the East Java subfault with a maximum magnitude of Mw 8.7. The minimum magnitude was Mw 7 for all faults, the magnitude bin was 0.1, and the maximum magnitude was different for each fault, as shown in Table 1.
Fig. 5
One of the slip distributions for East Java Subduction (Mw 8.9), and. The location of East Java Subduction (black rectangle)
To estimate the number of samples that are sufficient for PTHA calculation, a statistical analysis will be conducted using variance analysis methods as earlier applied for PTHA in the Sea of Japan (Mulia et al, 2020). This method measures the variability of samples around the mean by using the equation
. Delta (
) is the standard deviation, and µ is the mean of the maximum tsunami wave height at the coastal points, resulting from the Green’s function multiplication with the stochastic slip. The calculation is applied for all the magnitudes; in this case, the samples range from 2 to 110. The variance analysis result for the West Java Subduction fault is presented in Fig. 6. The black dashed line marks the minimum number of samples, beyond which the samples no longer contribute to the variability of tsunami heights indicated by the saturated curve of the coefficient of variation. We applied the same calculation for all the faults (Supplementary Fig. 2 – Fig. 10). The minimum samples required for this project are about 40 to 100 samples, depending on the magnitudes. However, for conciseness, we opted to use 110 samples for all faults, as it does not affect the computing time; consequently, the total number of scenarios considered in this study is 10.120 samples (Table 1).
Fig. 6
Monte Carlo calculation for East Java subduction.
2.4 Comparison with historical events and tsunami hazard assessment
According to historical data (NGDC, 2025; Hamzah and Puspito, 2000), the tsunami event that occurred in the area of interest has the same source location as this study, is the tsunami Flores (1992) with the magnitude of Mw 7.8, and the tsunami in Bali Sea with the magnitude of Mw 6.4 (2018) and Mw 6.9 (2018). Those three events will be compared. The source locations are in the Bali Back-arc Thrust for the 2018 tsunami events, and the Flores Back-arc Thrust for the 1992 tsunami event. 12 locations record tsunami wave height according to historic data for the 1992 tsunami event, and 3 locations for the 2018 tsunami event (Fig. 7). The simulation results from both events are compared to the historical data, and the Normalized Root Mean Square Error (NRMSE) value is around 0.13 (1.53 m for the RMSE). Additionally, the k value and κ value from the Aida calculation are also estimated in order to compare the data and get the β value for the PTHA calculation. The results for k value are 0.56, and for κ or β are 1.37. In the PTHA, the value of κ will be utilized to account for the aleatory uncertainty related to our model's parameterization and resolution.
Fig. 7
Comparison between historical observations and modeled tsunami heights.
A PTHA is a systematic approach used to evaluate tsunami hazard by estimating the probability of a tsunami intensity, such as the maximum wave height exceeding a certain threshold at a specific coastal location within a specified time period (Geist & Lynett, 2014). These estimates are typically presented in the form of a tsunami hazard curve, which shows the occurrence rate of exceeding a given tsunami intensity as a function of probability (Griffin et al., 2019; Suppasri et al., 2015; Thio & Somerville, 2012).
The tsunami hazard calculation is based on the tsunami exceedance rate (
) at a given location, where the exceedance rate is defined by the following equation (Cornell, 1968):
4
where
is the total number of sources or faults considered,
is the number of magnitude levels considered in a given interval, and
is the occurrence rate of earthquakes with magnitudes greater than or equal to M min, calculated using the Gutenberg-Richter distribution (Ishimoto and Iida, 1939; Gutenberg & Richter, 1944):
5
The values of
and
are constants determined based on seismicity data in the study area. This data was obtained from PusGen (2017) and Yuliatmoko et al. (2022) based on historical data at each fault location (Table 2). The
value is a constant describing the level of seismic activity. The
value is a constant describing the ratio of large earthquakes to small earthquakes.
Value | East Java Subduction | Bali Subduction | Timor FTB 1 | Timor FTB 2 | Timor FTB 3 | Bali Thrust | Sumbawa Thrust | Flores Thrust 1 | Flores Thrust 2 |
|---|---|---|---|---|---|---|---|---|---|
a | 5.60 | 5.63 | 5.14 | 4.59 | 7.67 | 3.19 | 3.19 | 3.19 | 3.19 |
b | 1.08 | 1.11 | 0.76 | 0.56 | 1.18 | 0.76 | 0.76 | 0.76 | 0.76 |
The probability of each magnitude level is calculated using the Gutenberg-Richter distribution with the formula:
6
The probability of a particular magnitude
is obtained from:
7
The probability of a tsunami height
exceeding a tsunami level
for a particular magnitude m is calculated as:
8
where
is the cumulative normal distribution function,
is the median logarithm of the tsunami height for a given source and magnitude. Thio (2009) proposed a total
value of 0.519, which was also used by Horspool et al. (2014). In this study, the value of β was sourced from the comparison with historical data, resulting in the
value of 1.37.
The probability value of lambda (
) can be calculated using the following equation:
9
The results of the annual frequency of occurrence (
) calculation can be used for disaggregation analysis or the identification of the dominant sources that contribute to the tsunami hazard. We applied the deaggregation analysis only for the 20 strategic locations. This calculation will produce a percentage contribution from each fault source, which can then be used to determine which source has the most dominant influence at that location.
3. Result
3.1 Tsunami Hazard Curves
The tsunami hazard curve describes the relationship between tsunami height and the annual rate of exceedance on a logarithmic scale. Hazard curves for all scenarios (Fig. 10) show a decreasing annual exceedance rate with increasing wave height, which is natural.
Fig. 10
Hazard curves for 20 strategic locations
Figure 10a depicts that the highest level of hazard in Bali is in Buleleng, Gianyar, Tanjung Benoa, and Jembrana, in which the annual rate of exceeding 1 m tsunami is approximately 10− 1.77, corresponding to a ~ 60-year return period (RP). In contrast, Badung demonstrates the lowest level of tsunami hazard among other cities in Bali. In NTB (Fig. 10.b), the highest rates (> 1 m) are in the southern and northern Sumbawa due to its locations that are exposed to the open seas; the lowest annual exceedance rates are in Bima, Mataram, and Central Lombok because of its position that blocked by an island (Bima) and in the middle of the strait (Mataram and Central Lombok). Lastly, for NTT (Fig. 10.c), the highest rates (> 1 m) are in Sikka because it is a bay near the open water and West Manggarai, which is directly influenced by the seas without any barrier. The lowest rates are in Ende, which may be because the tsunami will hit the cape first, and East Flores, as its nearby islands act as natural barriers.
The probability of experiencing a tsunami in percentage can be derived from the
value (Eq. 9). The first tsunami warning BMKG (2019) and Horspool et al. (2014) has the threshold of 0.5 m (“alert”), 0.5 m up to 3 m is “caution”, while the most dangerous wave height is considered more than 3 m (“watch out”). Figure 11 pictures the probability of a tsunami event for more than 0.5 m, 1 m, 2 m, and 3 m along the coast. Highest values (> 1.4% for > 0.5 m) occur in northern/southern Sumbawa, western Bali, and northern Flores. Extreme tsunami probabilities (> 3 m) are < 0.215%, with hotspots in the northern/southern Sumbawa and Bima.
Fig. 11
Maps exceedance probabilities for > 0.5, 1, 2, and 3 m wave heights.
According to ASCE (2017), BSN (2019), and USGS (2020), there are three tsunami probabilities considered for coastal infrastructures: 50%, 10% and 2% in 50 years. Therefore, the RP for those probabilities is 72, 475, and 2475 years. The tsunami wave height for such RPs for the Bali region is presented in Fig. 12, while for other regions it is shown in the Supplementary Figs. 11 and 12. Tsunami heights increase with the years of RP, consistent with probabilistic hazard principles. Spatial patterns or tsunami hazard map (Fig. 13) show substantial hazard in northern Sumbawa, southern Bali, and southern Sumbawa for all RPs. Tsunami > 3 m affects 0.4% of sites for 475-year RP and 57.3% for 2475-year RP. Heights range from 0.02–0.72 m (72-year), 0.02–3.3 m (475-year), and 0.12–11 m (2475-year).
Fig. 12
Hazard curves for three different return periods
Fig. 13
Hazard map for return periods of 72 years, 475 years, and 2475 years.
3.2 Deaggregation
The deaggregation analysis was applied in the 20 strategic locations along Bali, NTB, and NTT with the RPs of 72 years, 475 years, and 2475 years. The deaggregation results for 475 years and 2475 years are presented in the Supplementary Fig. 13–15. Figure 14 depicts the percentage of tsunami sources for the observation points in Bali. For Badung, Jembrana, and Gianyar, the Bali Subduction has the highest influence on tsunami hazard, while for Tanjung Benoa and Buleleng, the Flores Fault 2 has the highest influence. Not only does the distance matter in this analysis, but also the value of a, b, and the condition of bathymetry, and the propagation of the tsunami. Even the Flores Fault is situated farther than other Bali subduction or Bali back thrust arc, but the influence is greater. Three different RPs show the same dominant fault for Bali, East Nusa Tenggara, and West Nusa Tenggara, which are dominated by the Bali subduction zone. As for West Nusa Tenggara, the highest influence comes from Bali Subduction, Flores Thrust, Timor FTB 2, and Timor FTB 3. Bali subduction and the Sumbawa Thrust are dominant for the tsunami hazard in the East Nusa Tenggara. A previous study from Prasetya et al (2015) stated that the subduction zone usually contributes to a higher tsunami hazard than other types of faults. It is proven in this study that the Bali Subduction zone is the dominant source across the region.
Fig. 14
Deaggregation map for the Bali region (return period 72 years)
Location | Dominant Fault | |
|---|---|---|
Bali | Badung | Bali Subduction |
Tj Benoa | Flores Thrust 2 | |
Buleleng | Bali Subduction | |
Gianyar | Bali Subduction | |
Jembrana | Bali Subduction | |
West Nusa Tenggara | Bima | Timor FTB 3 |
Sumbawa | Flores Thrust | |
Sumbawa (North) | Flores Thrust | |
Central Lombok | Timor FTB 2 | |
Mataram | Bali Subduction | |
East Nusa Tenggara | Sikka | Bali Subduction |
East Flores | Bali Subduction | |
Kupang | Bali Subduction | |
Lembata | Bali Subduction | |
East Sumba | Bali Subduction | |
Ende | Sumbawa Thrust | |
Belu | Bali Subduction | |
Alor | Bali Subduction | |
West Manggarai | Bali Subduction | |
Manggarai | Bali Subduction | |
4. Discussion
The extreme wave heights in Klungkung and southern Sumbawa emphasize the potential for catastrophic local impacts from nearby large-magnitude subduction earthquakes. The mean wave heights, however, are much lower, suggesting that while extreme scenarios are possible, they are relatively rare. The hazard curve patterns indicate clear spatial variability in tsunami risk across the Lesser Sunda Islands. Coastal regions facing subduction zones (e.g., Bali’s southeast coast, southern Sumbawa) show higher probabilities for significant wave heights, consistent with their proximity to major tsunami-generating faults. This spatial heterogeneity matches findings from previous regional studies (e.g., Grezio et al., 2024; Suppasri et al., 2015)
The annual exceedance probability threshold of 10⁻² or 100-year RP for 1 m waves highlights priority zones for mitigation infrastructure, notably in Sumbawa, western Bali, and northern Flores. The exceedance probability maps also indicate that smaller tsunamis (< 1 m) have a much broader distribution, which, while less destructive, may still disrupt maritime activities and coastal operations. Return period maps reinforce the role of southern Bali, southern Sumbawa, and northern Sumbawa as persistent high-hazard zones, regardless of return period. This consistency across time horizons underlines the importance of long-term, location-specific tsunami preparedness.
Deaggregation results reveal that the Bali subduction zone is the most critical source for the majority of sites, with significant contributions from other thrust and fold-and-thrust belt sources depending on location. This reinforces the need for source-specific scenario planning and emphasizes that regional tsunami hazard cannot be attributed to a single fault system.
The probabilistic tsunami hazard assessment shows that tsunami contributions vary systematically across the islands, reflecting the influence of different fault systems (Supplementary Fig. 16). In Bali, the Bali Subduction Fault is the dominant contributor, with relative contributions ranging from 34.62% at Jembrana to 53.14% at Karangasem. Smaller but notable inputs come from the Flores Thrust (up to 12.48% at Karangasem) and the Lombok Fault (up to 18.55% at Karangasem). This pattern indicates that while Bali is primarily exposed to subduction-driven tsunamis, secondary sources such as the Flores and Lombok faults remain relevant.
In Lombok, the tsunami hazard is distributed more evenly among sources. The Bali Subduction contributes 20.74–27.46%, while the Sumbawa Fault provides 19.93–27.22%. Additional influence is observed from the Flores Thrust (11.45–18.71%) and Lombok Fault (13.47–15.29%). This balanced distribution suggests that tsunami hazard in Lombok arises from a combination of regional subduction and local crustal faults, rather than being dominated by a single source.
The Sumbawa region shows a different pattern, with the Flores Thrust emerging as the principal contributor. Percentages range from 28.53% in Dompu to 34.57% in Bima. The Sumbawa Fault also provides substantial contributions, particularly at Sumbawa Besar (30.64%) and Dompu (29.29%). In contrast, contributions from the Bali Subduction decrease to below 10% across most Sumbawa sites. These results highlight the importance of back-arc thrust activity in governing tsunami hazard along northern Sumbawa.
For Flores, Komodo, and Padar, the Flores Thrust is the dominant source, contributing 35.34% at Labuan Bajo, 38.27% at Komodo, and 34.71% at Padar. Secondary contributions vary by location: the Timor Fault accounts for 16.38% at Padar and 15.65% at Komodo, while the Sumba Fault contributes up to 15.63% at Komodo. These percentages suggest that although the Flores Thrust is the primary driver of hazard, multiple regional sources, including the Timor and Sumba faults, play a significant role in shaping local hazard potential.
Overall, Fig. 15 reveals a clear eastward shift in dominant fault contributions. In Bali, hazard is primarily subduction-driven, while in Lombok, hazard reflects a mixture of subduction and thrust activity. By Sumbawa, the Flores Thrust becomes the leading source, and this dominance intensifies further east in Flores, Komodo, and Padar. These findings emphasize that tsunami hazard in the Lesser Sunda Islands is not uniform, but instead tightly coupled to local tectonic settings. For risk management, this implies that mitigation strategies should account for island-specific dominant sources: subduction-driven hazards in Bali, combined subduction and thrust in Lombok, and Flores Thrust–dominated hazards in the eastern islands.
Fig. 15
Relative contributions of individual fault sources to tsunami hazard at selected sites across the Lesser Sunda Islands. Percentages represent the proportion of tsunami hazard attributable to each fault system.
5. Conclusion
Our PTHA for the Bali, West Nusa Tenggara (NTB), and East Nusa Tenggara (NTT) regions, considering both subduction and back-arc fault sources, incorporated stochastic slip distributions for nine faults, resulting in 10,120 earthquake scenarios generated using a 0.1 magnitude bin and 110 samples per magnitude. The estimated annual occurrence frequency ranged from less than 10⁻¹ (once in 10 years) to 10⁻⁴ (once in 10,000 years), with the highest frequencies observed in Gianyar, Jembrana, and Buleleng in Bali, Sumbawa in NTB, and Sikka and West Manggarai in NTT. Based on wave height classifications, hazardous tsunamis (> 3 m) were projected to affect 0.4% of locations for a 475-year (RP) and 57.3% for a 2,475-year RP, with maximum heights of 0.72 m, 3.3 m, and 11 m for RPs of 72, 475, and 2,475 years, respectively. Overall, the greatest contributions to tsunami wave heights at 20 strategic locations were attributed to the Bali subduction zone, the Flores Thrust, the Sumbawa Thrust, and Timor FTB 2 and 3. The analysis further reveals a distinct eastward shift in dominant tsunami sources: in Bali, hazard is primarily driven by the Bali Subduction, while in Lombok it reflects a balance between subduction and local thrust faults. In Sumbawa, the Flores Thrust becomes the main contributor, and its dominance continues into Flores, Komodo, and Padar, where additional input from the Timor and Sumba faults is also observed. These findings underscore the spatial heterogeneity of tsunami hazard across the Lesser Sunda Islands and highlight the need for island-specific mitigation strategies tailored to their dominant seismic sources.
Acknowledgement
We thank Prof. Dr. Nining Sari Ningsih and Dr. Astyka Pamumpuni for the inputs regarding the methods and results, and Gandhi Napitupulu, M.Si for the additional analysis.
A
Author Contribution
S. Karima: Performed data processing, analysis, drafted, and conducted a thorough review of the manuscript. Hamzah Latief: Concept the research, revised the analysis and drafted. Iyan E. Mulia: Provided a thorough review of the manuscript and revised the substance of the content.
Declarations
Conflicts of Interest
The authors declare no conflict of interest.
A
Funding
The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.
Ethical Approval
Not Applicable.
Electronic Supplementary Material
Below is the link to the electronic supplementary material
References
A
American Society of Civil Engineers (ASCE) (2017) Minimum design loads and associated criteria for buildings and other structures. Minimum Design Loads and Associated Criteria for Buildings and Other Structures. American Society of Civil Engineers (ASCE). https://doi.org/10.1061/9780784414248Audley-Charles MG (1986) Rates of Neogene and Quaternary tectonic movements in the southern Banda Arc based on micropalaeontology. J Geol Soc 143(1):161–175. https://doi.org/10.1144/GSJGS.143.1.0161
A
Badan Meteorologi (2019) Klimatologi, dan Geofisika (BMKG). Katalog gempabumi signifikan dan merusak 1821–2018 [Catalog of significant and damaging earthquakes 1821–2018] [in Indonesian]. BMKGA
Badan Pusat Statistik (BPS) (2023) Statistik jumlah penduduk Nusa Tenggara Barat tahun 2023 [Statistical data on the population of West Nusa Tenggara in 2023] [in Indonesian]. Badan Pusat Statistik. https://ntb.bps.go.id/idA
Badan Standardisasi Nasional (2019) SNI 1726:2019 Tata cara perencanaan ketahanan gempa untuk struktur bangunan gedung dan nongedung [SNI 1726:2019 Procedures for earthquake resistance planning for building and non-building structures] [in Indonesian]. Badan Standardisasi Nasional. https://pesta.bsn.go.id/produk/detail/12762-sni17262019Cornell CA (1968) Engineering Seismic Risk Analysis. In Bulletin of the Seismological Society of America (Vol. 58, Issue 5)
Felix RP, Hubbard JA, Bradley KE, Lythgoe KH, Li L, Switzer AD (2022) Tsunami hazard in Lombok and Bali, Indonesia, due to the Flores back-arc thrust. Nat Hazards Earth Syst Sci 22(5):1665–1682. https://doi.org/10.5194/NHESS-22-1665-2022
Geist EL, Lynett PJ (2014) Source Processes for the Probabilistic Assessment of Tsunami Hazards. Source: Oceanogr 27(2):86–93. https://doi.org/10.2307/24862158
Goda K, Yasuda T, Mori N, Maruyama T (2016) New scaling relationships of earthquake source parameters for stochastic tsunami simulation. Coastal Eng J 58(3):1650010–1650011
Grezio A, Anzidei M, Baglione E, Brizuela B, Di Manna P, Selva J, Taroni M, Tonini R, Vecchio A (2024) Including sea-level rise and vertical land movements in probabilistic tsunami hazard assessment for the Mediterranean Sea. Sci Rep 14(1). https://doi.org/10.1038/s41598-024-79770-9
Griffin J, Nguyen N, Cummins P, Cipta A (2019) Historical earthquakes of the eastern sunda arc: Source mechanisms and intensity-based testing of Indonesia’s national seismic hazard assessment. Bull Seismol Soc Am 109(1):43–65. https://doi.org/10.1785/0120180085
Gunawan E, Kongko W, Kholil M, Widyantoro BT, Widiyantoro S, Supendi P, Hanifa NR, Anjasmara IM, Pratama C, Gusman AR (2022) The 2019 Mw 7.0 Banten, Indonesia, intraslab earthquake: investigation of the coseismic slip, tsunami modelling and Coulomb stress change. Geoenvironmental Disasters 9(1):1–12. https://doi.org/10.1186/S40677-022-00215-4/FIGURES/11
Gutenberg B, Richter CF (1944) Frequency of earthquakes in California. Bull Seismol Soc Am 34(4):185–188
Green G (1838) On the motion of waves in a variable canal of small depth and width. Trans Camb Philosophical Soc 6:457
A
Hamzah L, Puspito NT, Imamura F (2000) Tsunami catalog and zones in Indonesia. J Nat Disaster Sci 22(1):25–43Harris R, Meservy W, Sulaeman H, Bunds M, Andreini J, Sharp B, Berrett B, Whitehead J, Carver G, Setiadi G, Hapsoro S, Prasetyadi C (2024) Discovery of imbricated beachrock deposits adjacent to the Java trench, Indonesia: influence of tsunami and storm waves, and implications for mega-thrust earthquakes. Nat Hazards: J Int Soc Prev Mitigation Nat Hazards 120(9):8209–8238. https://doi.org/10.1007/S11069-023-06327-W
Hasanah MU, Supendi P, Nugraha AD, Widiyantoro S, Syaifuddin F (2024) The new insight of tectonic setting in Sunda-Banda transition zone using tomography seismic. Case study: 7.1 M deep earthquake 29 August 2023. Open Geosci 16(1). https://doi.org/10.1515/GEO-2022-0710/ASSET. /GRAPHIC/J_GEO-2022-0710_FIG_007.JPG
Horspool N, Pranantyo IR, Griffin J, Latief H, Natawidjaja D, Kongko W, Thio HK (2013) A national tsunami hazard assessment for Indonesia. Australia–Indonesia Facility for Disaster Risk Reduction, Australian Government Department of Foreign Affairs and Trade, Canberra
Horspool N, Pranantyo I, Griffin J, Latief H, Natawidjaja DH, Kongko W, Cipta A, Bustaman B, Anugrah SD, Thio HK (2014) A probabilistic tsunami hazard assessment for Indonesia. Nat Hazards Earth Syst Sci 14(11):3105–3122. https://doi.org/10.5194/NHESS-14-3105-2014
Mai PM, Beroza GC (2002) A spatial random field model to characterize complexity in earthquake slip. J Geophys Research: Solid Earth 107(B11) 10 – 1. https://doi.org/10.1029/2001JB000588
Mulia IE, Ishibe T, Satake K, Gusman AR, Murotani S (2020) Regional probabilistic tsunami hazard assessment associated with active faults along the eastern margin of the Sea of Japan. Earth Planets and Space 72(1):1–15. https://doi.org/10.1186/S40623-020-01256-5/FIGURES/10
Newcomb KR, McCann WR (1987) Seismic history and seismotectonics of the Sunda Arc. J Geophys Research: Solid Earth 92(B1):421–439
Nugraha J, Pasau G, Sunardi B, Widiyantoro S (2014) Analisis Hazard Gempa dan Isoseismal Untuk Wilayah Jawa-Bali-NTB. Jurnal Meteorologi Dan Geofisika 15(1). https://doi.org/10.31172/JMG.V15I1.168
Okada Y (1985) Surface Deformation Due to Shear And Tensile Faults in A Half-Space. Bull Seismol Soc Am 75(4):1135–1154
Pararas-Carayannis G (2018) Brief History of Early Pioneering Tsunami Research – Part A. Science of Tsunami Hazards, 37(1), 49–129. https://doaj.org/article/5c7b49cac0794266a66ab3eb13795451
A
Pranantyo IR, Horspool N, Latief H, dan, Meilano I (2012) Studi Pendahuluan Analisis Probabilistik Bahaya Tsunami (PTHA) untuk Pulau Bali. 37th HAGI Annual Convention and Exhibition, 37Suarbawa KN, Sukarasa K, I., Riyono E (2021) Identification of Bali Island Deformation Based on Recorded GPS Data, Using GAMIT/GLOBK Software 10.6. Buletin Fisika 22(1):47–52
Suppasri A, Goto K, Muhari A, Ranasinghe P, Riyaz M, Affan M, Mas E, Yasuda M, Imamura F (2015) A Decade After the 2004 Indian Ocean Tsunami: The Progress in Disaster Preparedness and Future Challenges in Indonesia, Sri Lanka, Thailand and the Maldives. Pure appl Geophys 172(12):3313–3341. https://doi.org/10.1007/s00024-015-1134-6
Syamsidik S, Al M, Rasyif M, T., Fahmi M (2015) Banda Aceh, 23111 Indonesia. 2 Post graduate student at Civil Engineering Master Program. Jl. Syeh Abd. Rauff, 7. https://www.researchgate.net/publication/280134748
A
Thio HK (2012) Treatment of Uncertainties in Probabilistic Tsunami Hazard. In AGU Fall Meeting Abstracts (Vol. 2012, pp. NH24A-05)Thio HK, Somerville P (2009) A probabilistic tsunami hazard analysis of California. In tclee 2009: Lifeline earthquake engineering in a multihazard environment (pp. 1–12)
Tri Laksono FA, Aditama MR, Setijadi R, Ramadhan G (2020) Run-up Height and Flow Depth Simulation of the 2006 South Java Tsunami Using COMCOT on Widarapayung Beach. IOP Conference Series: Materials Science and Engineering, 982(1). https://doi.org/10.1088/1757-899X/982/1/012047
U.S. Geological Survey (USGS) (2020) Seismic design maps and hazard tools. U.S. Department of the Interior. https://www.usgs.gov/programs/earthquake-hazards/design-ground-motions
Wibowo M, Kongko W, Hendriyono W, Karima S (2021) Tsunami Hazard Potential Modeling as Tourism Development Considerations in the North of Lombok Strait. IOP Conference Series: Earth and Environmental Science, 832(1). https://doi.org/10.1088/1755-1315/832/1/012047
A
Yuliatmoko RS, Perdana YH, Martha AA (2021) Distribusi Frekuensi Gempabumi dan Dimensi Fraktal pada Seismik Gap di Indonesia. Jurnal Meteorologi dan Geofisika 22(2):55–56Zhao S, Cummins PR, McClusky S, Miller MS (2025) Interseismic Coupling Along the Java-Timor Subduction-Collision Zone at East Indonesia. Geophys Res Lett 52(6). https://doi.org/10.1029/2024GL112563;CTYPE:STRING:JOURNAL e2024GL112563