A probabilistic tsunami hazard assessment for Indonesia

Probabilistic hazard assessments are a fundamental tool for assessing the threats posed by hazards to communities and are important for underpinning evidence-based decision-making regarding risk mitigation activities. Indonesia has been the focus of intense tsunami risk mitigation efforts following the 2004 Indian Ocean tsunami, but this has been largely concentrated on the Sunda Arc with little attention to other tsunami prone areas of the country such as eastern Indonesia. We present the first nationally consistent probabilistic tsunami hazard assessment (PTHA) for Indonesia. This assessment produces time-independent forecasts of tsunami hazards at the coast using data from tsunami generated by local, regional and distant earthquake sources. The methodology is based on the established monte carlo approach to probabilistic seismic hazard assessment (PSHA) and has been adapted to tsunami. We account for sources of epistemic and aleatory uncertainty in the analysis through the use of logic trees and sampling probability density functions. For short return periods (100 years) the highest tsunami hazard is the west coast of Sumatra, south coast of Java and the north coast of Papua. For longer return periods (500–2500 years), the tsunami hazard is highest along the Sunda Arc, reflecting the larger maximum magnitudes. The annual probability of experiencing a tsunami with a height of > 0.5 m at the coast is greater than 10 % for Sumatra, Java, the Sunda islands (Bali, Lombok, Flores, Sumba) and north Papua. The annual probability of experiencing a tsunami with a height of > 3.0 m, which would cause significant inundation and fatalities, is 1–10 % in Sumatra, Java, Bali, Lombok and north Papua, and 0.1–1 % for north Sulawesi, Seram and Flores. The results of this national-scale hazard assessment provide evidence for disaster managers to prioritise regions for risk mitigation activities and/or more detailed hazard or risk assessment.


Introduction
Indonesia has the third highest population exposure to tsunami in the World, with an estimated 5.5 million people at risk of once-in-500-years tsunami (UNISDR, 2013;Lovholt et al., 2014). In the past 100 years alone, there have been 24 tsunami events that have caused fatalities totalling over 231 900 (with 226 900 in the Indian Ocean tsunami alone) (NGDC, 2013). Tsunami events have caused significant runup, with many events exceeding 10 m in most parts of Indonesia (apart from Kalimantan), as can be seen by the historical run-up data in Fig. 1 (Latief et al., 2000).
Since the devastating Indian Ocean tsunami in 2004 there has been a large amount of work focused on reducing the impact of future tsunami in Indonesia. This has included the completion of detailed tsunami inundation studies for  Latief et al. (2000) and the NOAA NGDC database. Geographic locations referred to in the text are in bold black font.
high-risk locations such as Padang, Denpasar and Cilacap (Kongko and Schlurmann, 2011;Strunz et al., 2011;Taubenbock et al., 2013), which have included a large focus on preparing communities for tsunami evacuation. In addition, the creation of the Indonesian tsunami early warning system (InaTEWS) has provided a platform to disseminate warnings less than five minutes after an earthquake (Lauterjung et al., 2010;Steinmetz et al., 2010;Strunz et al., 2011), albeit with a scenario database that exists only for the Sunda Arc and not other parts of Indonesia. Despite this large amount of work, there is still no nationally consistent assessment of tsunami hazard or risk, as most of the previously mentioned work has focused on the Sunda Arc. On one hand this seems justified, as the Sunda Arc has the potential to host the largest earthquakes that may cause tsunami that impact Indonesia; however on the other hand, there has been more tsunami events in the past 100 years in eastern Indonesia than in western Indonesia (Latief et al., 2000;Lovholt et al., 2012a;NGDC, 2013), as shown in Fig. 1.
To effectively address and mitigate the hazard posed by tsunami, accurate and uniform information on tsunami hazard is required for the basis of risk studies to prioritise the implementation of mitigation measures. The probabilistic tsunami hazard assessment (PTHA) framework is one such method that allows broad scale assessments of tsunami hazard (Lin and Tun, 1982;Rikitake and Aida, 1988;Geist and Parsons, 2006;Burbidge et al., 2009;Power et al., 2012;Sorensen et al., 2012). Taking a probabilistic approach for national and regional tsunami hazard assessments is favourable because it provides a geographically consistent, long-term assessment of the tsunami hazard, includes tsunami generated from multiple earthquake sources and incorporates uncertainty in the model parameters . Furthermore, it provides the first step towards estimating probabilistic tsunami risk, which is critical for under-pinning evidence-based decision-making about tsunami risk mitigation measures.
This paper presents the first national probabilistic tsunami hazard assessment for Indonesia. The outputs are long-term forecasts of tsunami hazard at the coast and are presented as tsunami hazard curves which describe the annual probability of exceeding a range of tsunami heights for sites around the coast of Indonesia. These hazard curves are then converted to hazard maps describing the tsunami height at the coast (1 m water depth) for a range of return periods and probability maps illustrating the annual probability of exceeding a set of given tsunami heights related to the warning thresholds of the Indonesian tsunami early warning system (InaTEWS). The tsunami hazard is also disaggregated to identify unit sources that contribute most to the tsunami hazard for a given coastal location, which can subsequently be used to select scenarios for detailed inundation modelling. Together, these results provide a nationally consistent assessment of tsunami hazard for Indonesia that can be used to underpin tsunami risk mitigation activities.

PTHA framework
The aim of a PTHA is to calculate the probability of exceeding a set of tsunami heights at the coast or near shore. This is accomplished by superposition of results from simulations of unit source tsunami to increase computation efficiency.
The PTHA framework can be summarised as follows: In this assessment the hazard is defined as the tsunami wave amplitude, where the term amplitude is identical to wave height, keeping with common convention in the tsunami hazard literature.

Seismic sources
For this assessment only tsunami generated by submarine earthquakes are considered. While tsunami can be generated from other mechanisms, such as submarine mass failures, submarine and sub-aerial volcanic eruptions and meteorite impacts, it is difficult to associate probabilities of such events, which prohibits including them in a probabilistic framework (Harbitz et al., 2014;ten Brink et al., 2014). Furthermore, the return period for significant tsunami generated by these alternative mechanisms is large (e.g. > 5000 years, Brune et al., 2010), so they are unlikely to influence the hazard estimates for short return periods (100-2500 years), which are the main interest for this study. In addition, 92 % of tsunami observed in Indonesia in the past 100 years have been caused by submarine earthquakes. Consequently, including only earthquake sources captures most of the tsunami hazard (Latief et al., 2000).
Seismic sources of tsunami that could impact Indonesia are divided into local, regional and distant sources. While it is assumed the hazard for Indonesia will be dominated by local sources, regional and distant sources are included to provide an all-encompassing assessment of the tsunami hazard. Regional sources are those likely to provide short warning times (less than a few hours), while distant sources are those that will have long warning times (many hours). Furthermore, this also enables a database of earthquake sources and tsunami events to be generated, which can be used for future tsunami inundation modelling or integration into the existing InaTEWS (Steinmetz et al., 2010).  Table 1) that represent subduction megathrusts (e.g. Sunda Arc, North Papua and North Sulawesi) and crustal faults are included in the local seismic source model. Crustal faults include submarine thrust faults (e.g. Flores back arc fault) and large submarine strike-slip faults that have an oblique rake component that may generate tsunami (e.g. Sorong fault). The recurrence parameters for the local fault sources (excluding the maximum magnitude, M max ) are taken from the recent national probabilistic seismic hazard map for Indonesia (Irsyam et al., 2008(Irsyam et al., , 2010 which represents the current state of knowledge about earthquake recurrence in Indonesia. Furthermore, this maintains consistency in earthquake recurrence rates between the current national seismic and tsunami hazard maps. The reason for not using the M max of Irsyam et al. (2010) is that the seismic hazard assessment for Indonesia has a less conservative estimation of M max based on instrumental catalogues. These may underestimate the maximum magnitude event due to their short records (e.g. 50 years) relative to actual recurrence intervals of M max , which may be 500-1000 years. Tsunami hazard assessments are much more sensitive to M max than probabilistic seismic hazard assessment (PSHA) because tsunami heights do not saturate with increasing magnitude as seismic ground motions do . For the PTHA, a conservative approach was taken to assigning M max , where the M max for each source was calculated using magnitude scaling laws from Blaser et al. (2010). First the length of the fault was determined, then the corresponding magnitude for this length was calculated from the scaling laws. In order to incorporate the uncertainty in scaling laws (e.g. the spread in measured values about the mean regression line), the M max was increased and decreased by 0.2 magnitude units and included in a logic tree (see following section on epistemic uncertainty). This accounts for events with smaller lengths compared to their magnitude; for example, the 11 March 2011 M w Tohokua earthquake was much smaller than the mean length for a M w 9.0 event de- rived from the scaling laws by Simons et al. (2011). The approach used here of assigning M max based on fault length is less subjective than using expert judgement. However, it is sensitive to the scaling law adopted.
The dip of seismic sources is a critical parameter in any tsunami modelling, since the resulting tsunami heights are sensitive to the fault dip (Geist, 1999(Geist, , 2002Lovholt et al., 2012b). The dip of local crustal faults was sourced from published estimates obtained from seismic reflection data (Rigg and Hall, 2012;Watkinson et al., 2011Watkinson et al., , 2012 and from the global CMT catalogue (Ekstrom et al., 2012). For the Sunda Arc, the fault geometry used by the InaTEWS was adopted to allow compatibility between the two tsunami databases Strunz et al., 2011). The InaTEWS fault model for the Sunda Arc is based on the regionalised upper mantle seismic model (RUM) of Gudmundsson and Sambridge (1998), which has variable dip and subfault dimensions (with an average of 40 × 25 km) along the length of the trench. The depth of the seismogenic zone for the sources was set between 3 and 50 km for subduction megathrusts, and 3 and 20 or 30 km for crustal faults. These parameters were taken from the Indonesian seismic hazard map by Irsyam et al. (2010). For subduction megathrusts and crustal thrust faults a rake of 90 • (pure thrust) was assumed; for large strike-slip faults (e.g. Sorong Fault) the rake was estimated by taking the average rake from large events (M w > 7.0) that have occurred on these faults. The focal mechanism data are sourced from the global Centroid Moment Tensor (GCMT) catalogue (Ekstrom et al., 2012).  Irsyam et al. (2008Irsyam et al. ( , 2010 and regional and distant sources are from Burbidge et al. (2008

Regional and distant sources
The regional and distant source model is represented by twelve subduction megathrusts around the western and northern Pacific Ocean and the Makran subduction megathrust in the northwestern Indian Ocean. The subduction zone geometry and recurrence rates were taken from the PTHA for Australia (Burbidge et al., , 2009, which uses seismicity data to assign annual seismic moment rates to the plate interface. Sources in the eastern Pacific Ocean are assumed not to influence tsunami hazard for Indonesia, based on historical data and preliminary modelling.

Unit sources
Each source is discretised into rectangular subfaults prior to the tsunami summation calculation. The dimensions of the subfault represent the smallest magnitude considered in the synthetic catalogue. For local subduction sources, each subfault measures 40 km along strike by 25 km down dip, which is approximately a M w 7.0 earthquake (Blaser et al., 2010). Local crustal sources have dimensions of 20 km × 10 km, which corresponds to an M w 6.4 event (Blaser et al., 2010). Regional and distant subfaults are discretised into larger subfaults (100 km × 50 km) since the minimum magnitude that could generate a tsunami that would impact Indonesia is much larger. This subfault size is equivalent to a M w 7.6 event (Blaser et al., 2010). A total of 4220 subfaults were generated for the 32 seismic sources, with the Sunda Arc alone divided into 1800 subfaults (Fig. 3). The vertical deformation is calculated using the layered earth model of Wang et al. (2003) and the crustal parameters from the Crust 2.0 model (Bassin et al., 2000).

3110
N. Horspool et al.: A probabilistic tsunami hazard assessment for Indonesia

Unit source method
The use of unit sources have become increasingly popular for PTHA (Geist and Parsons, 2006;Burbidge et al., 2009;Power et al., 2012) and tsunami early warning systems (Titov et al., 2005;Greenslade et al., 2011) due to their computational efficiency at calculating the tsunami hazard for a large number of sites over a large area. The fundamental assumption in using a unit source approach is that tsunami waves behave linearly in deep to moderate water depths . This assumption means that any tsunami can be constructed by the summation and scaling of smaller individual tsunami generated by each of the subfaults that define the larger earthquake rupture (Fig. 3). This means that the computationally intensive tsunami propagation models need only be computed once for each subfault. From these unit sources any number of earthquakes and resulting tsunami can be calculated efficiently by summing and scaling the tsunami time series from individual unit source tsunami.

Tsunami propagation simulations
Tsunami propagation modelling was carried out using the finite difference code of Satake (1995), which solves the linear shallow water wave equation over a regular grid in spherical coordinates (WGS84 projection). This code has been used extensively for modelling tsunami scenarios Fujii and Satake, 2006;Baba et al., 2008) and for PTHA Burbidge et al., 2008Burbidge et al., , 2009). Tsunami propagation simulations were performed for 1 m of slip on each of the 4220 subfaults and run for a model time of 6 h for local sources and 24 h for regional and distant sources. The model allows a variable computational grid with a coarse master grid, and finer nested grids around areas of interest. The coarse grid (Fig. 2) has a spatial resolution of 3 arcmin and was derived from the 30 arcsecond general bathymetric chart of the oceans (GEBCO) model. For a region surrounding Indonesia (Fig. 1), a 1 arcmin grid was used that was derived from a commercial 90 m resolution bathymetric data set from TCarta Marine. This is a commercial bathymetry data set combining Navy Charts, multibeam surveys and GEBCO data. A reflective boundary at the coast was used for the propagation modelling. For each propagation simulation, tsunami time series were recorded at 3050 hazard points located at the 100 m water depth contour at approximately 5-10 km intervals around the Indonesia coastline. In water depths shallower than 100 m, non-linear terms in the shallow water wave equation can become important and the unit source approach can not be applied.
Since the goal of the study is to estimate the probabilistic tsunami hazard at the coast, we extrapolate the tsunami height from the 100 m water depth contour to the coast (1 m water depth) using Green's Law. This approach has been used in Japan (Kamigaichi, 2009) and the Mediterranean (Sorensen et al., 2012) for similar studies. Kamigaichi (2009) demonstrated that extrapolation using Green's Law was consistent with the mean offshore amplitude obtained when using detailed non-linear modelling on high-resolution grids for a case study in Japan.

Treatment of uncertainty
Incorporating uncertainties is a critical component of any probabilistic assessment. For this study there are two types of uncertainty that are dealt with: aleatory uncertainty and epistemic uncertainty. Aleatory uncertainty is the variability from randomness in a natural process and cannot be reduced with increasing knowledge. In PTHA, an example would be the variability of slip of a future tsunamigenic earthquake. Aleatory uncertainty can be accounted for by sampling from a probability density function (Bommer and Abrahamson, 2006). Epistemic uncertainty is the uncertainty stemming from our lack of knowledge about a natural process. Epistemic uncertainty is often referred to as scientific uncertainty. With increased data and knowledge the epistemic uncertainty can theoretically be reduced to zero. An example of epistemic uncertainty in PTHA is whether a subduction zone interface will only rupture in single segments or is able to rupture as a multi-segment rupture. Logic trees are a tool used to capture and quantify the epistemic uncertainty and are used widely in PSHA Bommer et al., 2005;Bommer and Abrahamson, 2006).

Epistemic uncertainties
Sources of epistemic uncertainty that are included in the PTHA are fault segmentation, slip rate, magnitude-frequency distribution (MFD) type and maximum magnitude. Variations in the values for these parameters are included in logic trees for each source zone and represent different possible models for the input parameters (Fig. 4). For most source zones, the logic trees include variation in the MFD type, maximum magnitude and slip rate if different published estimates exist ( Table 1). The weighting and structure of the logic tree was developed during a consensus style  technical working group meeting of key tsunami and earthquake researchers from Indonesia, and as such the weights assigned to the logic tree are subjective and biased by the workshop participants' points of view and opinions. All weights assigned in the following discussion were agreed upon by the participants based on current understanding of the seismic behaviour of submarine fault sources in Indonesia.  The Sunda Arc has a complex logic tree ( Fig. 4 and Table 1) which reflects the different interpretations of the earthquake mechanics of the megathrust. Two main branches represent whether the Sunda Arc could (a) rupture in a single large rupture anywhere along the length of the megathrust (non-segmented model) assuming we do not have enough data to suggest this could not occur, or (b) rupture in segments (segmented model), as observed by the recent Sumatran megathrust earthquake sequence from 2004 to present (Natawidjaja et al., 2007;Chlieh et al., 2008;Sieh et al., 2008), where we assume ruptures are stopped by structural barriers in the subducting plate. For subsequent branches the maximum magnitudes are constrained by scaling laws (Blaser et al., 2010). For all logic trees, the maximum magnitude from scaling laws is given a weighting of 0.6; two alternative maximum magnitudes that are ±0.2 magnitude units from the mean of the scaling law are given a weighting of 0.2. This allows uncertainty in the scaling law prediction error to be incorporated into the PTHA. Table 1 shows the source parameters used in the logic trees with their associated weighting. To incorporate different expert opinion on the magnitude frequency distribution (MFD) shape for each source, two different MFDs were used. For each source, a truncated Gutenberg-Richter MFD (Gutenberg and Richter, 1994) was given a weighting of 0.66 and a characteristic earthquake MFD (Schwartz and Coppersmith, 1984) was given a weighting of 0.34. A b value of 1.0 was used for both MFDs.

Aleatory uncertainties
There are three main sources of aleatory uncertainty in the PTHA: computational modelling uncertainty (σ M ), uncertainty in the geometry of the sources (σ D ) and uncertainty due to random slip distributions (σ S ) . Uncertainty due to the aforementioned sources is described as the standard deviation of a log-normal distribution with a zero mean.
Modelling uncertainty can be described as errors in numerical models that are used to model seafloor displacement and tsunami propagation. By comparing tsunami waveforms from historical events to modelled waveforms of the event, an estimate of this uncertainty can be quantified. Thio (2012) estimated this uncertainty from validation of several wellconstrained tsunami on the California coast and estimated a standard deviation (σ M ) of 0.345.
Uncertainty in the dip of the fault would normally be treated as an epistemic uncertainty, since increased knowledge of the fault geometry, for example from seismic reflection data, could reduce the uncertainty of the dip. However, this would require recomputing the tsunami unit sources for each variation in dip which would make this approach more computationally intensive and may eliminate the benefit of using tsunami unit sources. Instead, the dip uncertainty is included as an aleatory uncertainty. Thio (2012) modelled a range of scenarios for a range of dip angles with a mean of 10 • and a standard deviation of 5 • , which resulted in a standard deviation of the tsunami height (σ D ) of 0.292.
The distribution of earthquake slip in future earthquakes is assumed to be a random process. When the tsunami heights are calculated from the events in the synthetic earthquake catalogue the slip is modelled as uniform across all subfaults. To include this uncertainty in the tsunami unit source summation would require a large number of slip distributions to be modelled for every event in our catalogue (> 100 000 events). A more efficient way is to treat this uncertainty similar to the dip, in which a large number of slip distributions are modelled for a single scenario event and a suite of resulting tsunami heights are estimated. The spread of the tsunami heights from different slip distributions can be used to estimate uncertainty resulting from this randomness and be represented as a probability density function against which the final wave heights can be evaluated to include this uncertainty. The sensitivity testing of slip distribution results in a standard deviation of tsunami height of 0.256 (Thio, 2012).

N. Horspool et al.: A probabilistic tsunami hazard assessment for Indonesia
The three sources of aleatory uncertainty are combined to compute a total sigma term (σ total ): which results in a σ total of 0.52.

Synthetic earthquake catalogue
Combining all the information from the source models and logic trees, a synthetic catalogue is generated which represents the range of possible earthquake magnitudes, locations and sources for every logic tree branch. The catalogue was generated by iterating through each magnitude in the MFD, from the minimum magnitude of M w 7.0 to M max , and calculating the rupture dimensions using the scaling laws for subduction zone earthquakes from Blaser et al. (2010) and crustal fault earthquakes from Wells and Coppersmith (1994). The rupture is then iteratively moved across the fault one subfault at a time until that magnitude has occurred on every possible location within the fault dimensions. For the megathrust sources, the subfault dimensions are 40 km × 25 km for all depth ranges, which is approximately a M w 7.0 earthquake; therefore the number of ruptures would be equivalent to the number of subfaults. The maximum magnitude earthquake would occur once and rupture the whole fault if scaling laws have been used to constrain the maximum magnitude. This iterative process ensures that all magnitudes could occur at any possible location on the fault plane. For each event, the probability of that magnitude was then weighted by 1/n, where n is the number of earthquakes represented by that magnitude. This ensures that the sum of the probabilities of all events of the same magnitude equals the annual probability of one event of that magnitude in the MFD. For each event the average slip is calculated using the fault scaling laws of Blaser et al. (2010) for subduction zone earthquakes and from Wells and Coppersmith (1994) for crustal fault earthquakes. A total of 104 000 synthetic earthquakes were generated for the catalogue across the 32 fault sources. This implementation of a PTHA assumes that earthquakes follow a Poisson process; that is, they are time-independent. In a Poisson process, events vary by a mean recurrence time. This does not consider the process of time-space clustering of earthquakes such as aftershocks. In PSHA, Poisson models have been shown to be valid for intermediate to long-term hazard forecasts (Helmstetter et al., 2006;Petersen et al., 2007), and can generally apply to PTHA because the source is also from earthquakes.

Probabilistic calculations
For each event in the synthetic catalogue, the tsunami hazard is calculated at each hazard point along the coast by summing the tsunami unit source contributions from the subfaults that make up that event and by scaling the tsunami height by the average event slip. For each site, this results in a list of tsunami heights (H i ) and associated annual probabilities (P i ) for each event (i). For a range of given tsunami heights (e.g. x = 0.1, 0.2, 0.3, 0.5, 1.0 . . . 30.0) we then compute the probability of exceeding that height by incorporating the modelled tsunami heights and the aleatory uncertainty term, σ total . The probability of exceeding a given tsunami height is as follows: where f total (u) is the PDF of the total aleatory uncertainty and P i is the earthquake probability for event i. Next we can calculate the rate at which H > x for all events (i) across all sources (n): By plotting the probability of exceeding a range of tsunami heights (λ(H > x)), a tsunami hazard curve can be constructed (Fig. 4) from which tsunami hazard maps and probability maps can be generated by iterating over each hazard curve. This model assumes a memory-less Poisson process where earthquakes are time-independent and occur at a fixed probability over time.

Results
The results from the PTHA are developed into hazard curves, hazard maps, probability maps and disaggregation maps. The maps are all derivative products of the hazard curves and can be used to provide information on different aspects of the tsunami hazard. For instance, the hazard maps illustrate the tsunami height expected at the coast over a given return period and are useful for understanding the size of tsunami expected for a community over a fixed time. Probability maps define the probability of exceeding a given tsunami height at the coast, which is linked to threat thresholds in the InaTEWS system. This is useful for assessing the annual probability of experiencing a tsunami warning, particularly for prioritising areas most likely to experience tsunami that are a threat to life safety. Disaggregation maps define the relative contribution of each source to the hazard at a particular site for a given return period or tsunami height. This information is extremely useful for identifying sources and selecting scenario tsunami events for detailed tsunami inundation modelling. tal output from the PTHA. The hazard curves are grouped together according to the tsunami zones in Indonesia as defined by Latief et al. (2000) and are shown in Fig. 4.

Tsunami hazard curves
The hazard curves in Western Indonesia show that the mean hazard for Sumatra and Java is of a similar magnitude, however the spread of the hazard curves is larger for Sumatra. This reflects the site location of the hazard curves in the Sumatra zone, as some are located on the eastern coast of the Mentawai and Nias islands that are protected from tsunami originating to the west of these islands. Furthermore, the hazard curves on the west coast of the offshore islands have higher hazard than the west coast of Sumatra and the south coast of Java, as they are located adjacent to the Sunda Arc source. The Java hazard curves have a distinct bulge above the 0.001 annual probability (∼ 1000-year return pe-riod). This is thought to be due to the influence of the maximum magnitude events in the characteristic MFD which has a high probability for large characteristic earthquakes.
In Eastern Indonesia, the tsunami hazard curves for the Banda, Papua and Sulawesi zones are similar. The mean probability of exceeding 1 m coastal tsunami height is similar (∼ 0.01) for these zones and for Java and Sumatra. However for larger coastal tsunami heights (> 1 m), Java and Sumatra have much higher probabilities. This is due to the size of M max for the sources that contribute to the hazard. Java and Sumatra have much higher M max (M w 9.5) compared to sources in eastern Indonesia which have M max of M w 8.5-9.0 (Table 1). There is a larger spread in the hazard curves in eastern Indonesia due to the presence of thousands of islands which add complexity to tsunami propagation paths.  Latief et al. (2000). Upper panel shows the location of the zones (black lines) and faults used in the PTHA (red lines). The grey curves show all the hazard curves (one for each site) in each zone, while the shading is a function of the density of the hazard curves. The red line represents the mean hazard curve for the zone weighted by the sample spacing between each coastal hazard curve location. The Kalimantan zone (including the north coast of Java and the east coast of Sumatra) has the lowest tsunami hazard of any zone due to being remote from any major tsunami sources besides the north Sulawesi megathrust and the Makassar thrust in Sulawesi. These areas are also adjacent to shallow (< 200 m) epicontinental seas that will attenuate tsunami more rapidly compared with deeper areas.
The hazard curves can be used to prioritise which areas in Indonesia have the highest tsunami hazard. The results show that eastern Indonesia has a similar tsunami hazard to western Indonesia at higher annual exceedance probabilities, whereas western Indonesia has a higher tsunami hazard at lower probabilities of exceedances.

Tsunami hazard maps
Hazard maps are generated for return periods of 100, 500 and 2500 years. Figure 5 shows the hazard map for the 100-year return period, indicating that the highest tsunami hazard is on the west coast of the Mentawai islands, offshore from the west coast of Sumatra. The maximum expected tsunami height at the coast of the Mentawai islands is 9 m. On the west coast of Sumatra the tsunami hazard ranges from 2 to 7 m, which is similar to the south coast of Java and the islands of Bali and Lombok. In eastern Indonesia, the north coast of Papua has the highest hazard at around 4 m for this return period; parts of the northern Molluccas and north Sulawesi have hazards of ∼ 2 m.
At the 500-year return period (Fig. 6), the tsunami hazard is still the highest along the west coast of the Mentawai and Nias islands. The hazard along the north Papua coast, northern Sulawesi and northern Molluccas is 10-12 m for this return period, similar to that of the south coast of Java and the northwest coast of Sumatra, near Lampung and Aceh. Areas in the Banda Arc have a tsunami hazard ranging from 2 to 6 m, with sites adjacent to tsunami sources having a higher hazard.
The hazard for the 2500-year return period (Fig. 7) is clearly the highest along the Sunda Arc and is due to the strong contribution to the hazard from low probability events approaching M max , such as the 2004 Indian Ocean tsunami. In eastern Indonesia, particularly north Papua and Sulawesi, the hazard has saturated because the maximum magnitude is only M w 8.5, compared to M w 9.5 along the Sunda Arc. The variation in hazards for different return periods is a reflection of different activity rates and maximum magnitudes of the seismic sources.
In Sumatra it is evident that the hazard varies significantly along the west coast of Sumatra at all return periods due to the presence or absence of offshore islands. Where offshore islands are present, the coastal areas of Sumatra have a lower hazard, approximately half that of the islands. This is due to two factors. Firstly the islands are located directly above the megathrust; if islands are present during the coseismic displacement, the displaced water volume will be less, resulting in a smaller tsunami on the west coast of Sumatra directly east of the islands (e.g. Padang city). Secondly, the islands act as a barrier for tsunami that have slip near the trench, which results in the west coast of the offshore islands receiving the largest tsunami waves and the west coast of Sumatra being more protected. With very short arrival times to the Mentawai and Nias islands and the very high tsunami hazard, this area clearly is of concern for tsunami mitigation efforts. In comparison, the south coast of Java has less variation along the coast due to the absence of offshore islands.

Tsunami probability maps
Tsunami probability maps, which illustrate the probability of exceeding a given tsunami height, are less common in PTHA but are useful for prioritising tsunami mitigation activities based on the likelihood of having a significant tsunami. The InaTEWS sets its warning thresholds at 0.5 m for a tsunami warning and 3.0 m for a major tsunami warning (Strunz et al., 2011). A major tsunami warning corresponds to an event that may cause land inundation and pose a serious threat to life and property safety. Therefore, tsunami probability maps are generated for these thresholds and are shown in Fig. 8. Figure 8a shows that the probability of experiencing a tsunami with a height at the coast of 0.5 m or greater in any given year is over 10 % for many parts of Indonesia which have a seismic source nearby. The probability is the highest along the Sunda Arc and slightly less in northern Sulawesi and north Papua. The area around the Banda Sea has a 2 % probability of experiencing a 0.5 m tsunami in any given year.
The probability of experiencing a tsunami with a height at the coast greater than 3.0 m, which would trigger a major tsunami warning, varies along the Sunda Arc depending on whether offshore islands are present or not. The probability that the Mentawai and Nias islands would experience a 3.0 m or greater tsunami in any given year is 5 %, compared to less than 1 % for the west coast of Sumatra, including the city of Padang. However, it is noted that the tsunami hazard in Padang is considered to be much higher in the short-term if time dependence is captured (Taubenbock et al., 2013). Most parts of eastern Indonesia have a probability of around 1 % of experiencing a major tsunami warning in a given year, with the highest probability of 2-5 % occurring along the north coast of Papua and the northern Molluccas.
These results highlight that areas in eastern Indonesia such as north Papua and north Sulawesi have similar probabilities of experiencing a major tsunami warning as western Indonesia; however, east of Lombok no tsunami scenarios exist in the InaTEWS database. The results presented here warrant  an extension of the InaTEWS to develop a scenario database for eastern Indonesia.

Hazard disaggregation
The PTHA methodology aggregates all possible tsunami events to generate tsunami hazard curves. The tsunami hazard can also be disaggregated by source and magnitude. Tsunami hazard disaggregation maps are developed to quantify how much each source contributes to the hazard for a single site and a single return period. They are useful for scenario event selection when detailed inundation modelling is undertaken. Figure 9 shows the tsunami hazard disaggregation for two populous coastal cities, Jayapura in north Papua and Mataram in Lombok, for a return period of 500 years. The tsunami hazard is disaggregated for each subfault to show which contribute the most to hazard. This is useful to define where scenario events should be located for detailed inundation modelling. Mataram, the capital of Lombok, is situated on the west coast of the island. It is located between the active Java Trench to the south and the Flores back arc thrust fault system to the north (Fig. 1). The 500-year return period hazard is dominated by the Flores back arc thrust, with a slightly lower contribution from the Java Trench. Jayapura is situated on the north coast of Papua and is located directly above the active New Guinea Trench which has produced tsunamigenic earthquakes in the past 20 years (Henry and Das, 2002;Tregoning and Gorbatov, 2004). However the disaggregation in Fig. 9 shows that tsunami hazard is dominated not only by the adjacent New Guinea Trench, but also the Philippines subduction zone and to a lesser extent by the Mariana and Ryukyu trenches. This highlights how tsunami hazard disaggregation is not only a function of the proximity to tsunami sources, but also to their activity rate. As Fig. 9 shows, the Philippines subduction zone, while further away 3118 N. Horspool et al.: A probabilistic tsunami hazard assessment for Indonesia has a high slip rate, meaning that over the long time frames considered in PTHA, it will contribute as much to the hazard for Jayapura as the local New Guinea Trench. For sites along the Sunda Arc the contribution is solely from the Sunda Arc megathrust. This contrasts with sites in eastern Indonesia where there are more local and regional sources, and the tsunami hazard comes from more than one source. Any detailed inundation modelling should use disaggregation information for event selection to ensure that all significant source zones are covered in the analysis.

Discussion
This study represents the first attempt at developing a nationally consistent probabilistic tsunami hazard assessment for Indonesia, one of the most tsunami-prone countries in the World (UNISDR, 2013). This study represents the first iteration of a national tsunami hazard map, but there are many aspects that need to be considered and improved upon for future versions. Nonetheless, this study can be used to facilitate evidence-based decision-making on tsunami risk mitigation activities and priorities in Indonesia.

Underpinning tsunami risk mitigation
Until now, the majority of tsunami risk mitigation activities have been located along the Sunda Arc, due to a wealth of data on the earthquake history of the Sumatra megathrust (Natawidjaja et al., 2007;Brune et al., 2010;Kongko and Schlurmann, 2011;Strunz et al., 2011;Taubenbock et al., 2013). However, the tsunami event database from the National Geophysical Data Centre (NGDC), shows that in the past 100 years, removing the 2004 Indian Ocean tsunami event, there have been 3651 fatalities from tsunami in eastern Indonesia (from 15 events at a rate of 1 fatal tsunami every 6.5 years) compared to 1759 fatalities along the Sunda Arc (from 10 events at a rate of 1 fatal tsunami every 10 years). These data, and the results from the PTHA, provide evidence that the tsunami hazard in eastern Indonesia is as high as that along the Sunda Arc for return periods less than 1000 years. At return periods greater than 1000 years, the M max events along the Sunda Arc dramatically increase the hazard for Java, Sumatra, Bali and Lombok. This highlights a clear need to include tsunami scenarios for eastern Indonesia within the decision support system of InaTEWS, which would strengthen the early warning system in eastern Indonesia, similar to that already established along the Sunda Arc (Steinmetz et al., 2010;Strunz et al., 2011), so that communities can receive warnings for local and regional tsunami events that may impact them. Furthermore, the results from the PTHA can be used as a prioritisation tool to allow disaster managers to prioritise focus areas for tsunami risk mitigation activities.

Poisson model
The PTHA assumes earthquakes that generate tsunami follow a Poisson process; that is, they are time-independent. This assumption is fundamental in PSHA and PTHA and over the long-term is usually a valid approach. However, earthquake faults have been shown to interact on short time scales, and earthquake sequences are now being recognised . For example, in Indonesia there is evidence of this through the recent Sumatra megathrust sequence , which started with the 26 December 2004 M w 9.15 Andaman earthquake and was followed by the 28 March 2005 M w 8.6 Nias earthquake, the 12 September 2007 M w 8.5 Bengkulu earthquake, the 25 October 2010 M w 7.7 Mentawai earthquake and the 11 April 2012 M w 8.6 earthquakes off the coast of Sumatra. This clustering in space and time highlights that there may be a time dependence to this earthquake sequence. Over the short term (i.e. the next 10-50 years), time-dependent models may predict the tsunami hazard better; however, time-dependent tsunami models have yet to be developed. In PSHA, timedependent models have been applied widely over the past 10 years (Wiemer, 2000;Gerstenberger et al., 2005). Results suggest they perform better in short-term forecasting than time-independent models in California, where rigorous statistical testing has occurred (Helmstetter et al., 2006;Petersen et al., 2007). However, the purpose of this PTHA is to provide a picture of the long-term tsunami hazard across Indonesia for the purpose of prioritising long-term risk mitigation activities. In this case, assuming a time-independent (Poisson) model is a valid assumption. The recent sequence of Sumatran earthquakes could provide an opportunity to test the application of time-dependent models for tsunami hazard assessments.

Limitations
A number of limitations are evident from our PTHA which will be improved upon in future model updates. The assessment presently only considers earthquake sources of tsunami and not other sources, such as submarine landslides (earthquake triggered or independent of earthquakes), volcanic collapse, or volcanic products such as pyroclastic flows. All of these sources have generated significant tsunami in Indonesia in the past 200 years. These events include the 1815 Tambora eruption that caused large pyroclastic flows that generated large tsunami (Self et al., 1984), the 1883 Krakatau eruption that caused a huge tsunami over 10 m in height in western Java and killed over 36 000 people (Mader and Gittings, 2006), and the 1979 Lomblen island tsunami, that was generated by a submarine landslide not associated with an earthquake and killed over 500 people (Brune et al., 2010). Although these alternative sources account for ∼ 8 % of tsunami in the past 100 years, they have the potential for large local impact. Presently it is difficult to assign probabil-ities to such events which precludes incorporating them in a probabilistic analysis (Harbitz et al., 2014;ten Brink et al., 2014). However, in other areas around the world, landslides are now being assigned probabilities, for example in California (Geist and Parsons, 2010). This work is promising and may allow landslides to be included in PTHA in the near future; however, this will require dealing with large uncertainties and the use of dispersive propagation models.
Since the 2004 Indian Ocean tsunami, the majority of tsunami hazard studies have focussed on Sumatra and to a lesser extent Java. A result of this focus is that the recurrence rate of previous earthquakes has been well-studied, using coral uplift data above the trench (Natawidjaja et al., 2007;Sieh et al., 2008), current seismic slip deficit, and locking constrained from extensive GPS networks (Chlieh et al., 2008). Meanwhile, in Java and eastern Indonesia, similar studies are rare, and the recurrence rate and geometry of submarine faults is poorly constrained. As a result, the uncertainties are large for key parameters that tsunami hazard is sensitive to, such as slip rate and dip. More work is needed in this tectonically complex area and, as new data are available, it will improve the accuracy of tsunami and seismic hazard forecasts.
Uncertainty arising from random slip distributions is incorporated into the PTHA by means of sampling a probability density function. A better approach may be to include heterogeneous slip by assigning different slip values to subfaults that are used for each event. While this approach would be more computationally intensive than the current implementation, it would produce a more direct means of incorporating slip uncertainty.
Finally, an area of further research is the inclusion of "tsunami earthquakes" in PTHA. Tsunami earthquakes are earthquakes with large tsunami relative to the magnitude of the event. At present tsunami earthquakes are not treated any differently in the PTHA. However, there are ways of including these in future PTHA studies, for example by allowing increase slip for earthquakes located on the shallow portion of the plate interface in areas where tsunami earthquakes are known to occur (e.g. Java Trench and Mentawai islands).

Testing of the PTHA
Since a PTHA is a statistical forecast of the future tsunami hazard, testing should be applied to rate the performance of these models (Geist and Parsons, 2006). In PSHA, statistical testing has become common, with specific model testing facilities set up around the world to test different seismic hazard forecasts. In particular, the Collaboratory for the Study of Earthquake Predictability (CSEP) (Schorlemmer et al., 2007) is leading the way in setting up routine testing of seismic hazard models. The most important testing methods compare seismic hazard forecasts to empirical ground motion data, which ultimately tests the final PSHA model outputs rather than single components of the PSHA (Albarello and D'Amico, 2008;Stirling and Gerstenberger, 2010). Future work should focus on rigorously testing PTHA models with observational data. However, given the limited amount of observational data in most regions, it may be necessary to integrate over large areas to develop a large enough sample of observations for testing purposes.

Conclusions
This study developed the first national probabilistic tsunami hazard assessment for Indonesia and provides forecasts of long-term tsunami hazard at the coastline from earthquake sources. Results show that the annual probability of experiencing a tsunami with a height of > 0.5 m at the coast is greater than 10 % for Sumatra, Java, Sunda islands (Bali, Lombok, Flores, Sumba) and north Papua. The annual probability of experiencing a tsunami with a height of > 3.0 m is 1-10 % in Sumatra, Java and north Papua, and 0.1-1 % for north Sulawesi, Seram and Flores. For long return period forecasts (2500-year return period) the tsunami hazard is the highest along the coast exposed to the Sunda Arc (i.e. Mentawai and Nias islands, Sumatra and Java), where tsunami heights of 25-30 m are expected. Comparatively, eastern Indonesia has lower maximum tsunami hazard at the 2500-year return period, on the order of 12-20 m in Sulawesi, Papua and Seram island. These results highlight that although the maximum tsunami hazard is higher in the Sunda Arc, the probability of experiencing a tsunami that could cause inundation and fatalities (coastal height of > 3.0 m) is similar along the Sunda Arc and parts of eastern Indonesia such as north Papua and north Sulawesi. This warrants the extension of the InaTEWS decision support system scenario database to include eastern Indonesia. Furthermore, the results from this study can be used to underpin tsunami risk assessments by highlighting high hazard areas where tsunami inundation modelling should be conducted.