Methodological reconstruction of historical seismic events from anecdotal accounts of destructive tsunamis: a case study for the great 1852 Banda arc mega-thrust earthquake and tsunami
Hayden Ringer, Jared P. Whitehead, Justin Krometis, Ronald A. Harris, Nathan Glatt-Holtz, Spencer Giddens, Claire Ashcraft, Garret Carver, Adam Robertson, McKay Harward, Joshua Fullwood, Kameron Lightheart, Ryan Hilton, Ashley Avery, Cody Kesler, Martha Morrise, Michael Hunter Klein
MMETHODOLOGICAL RECONSTRUCTION OF HISTORICALSEISMIC EVENTS FROM ANECDOTAL ACCOUNTS OFDESTRUCTIVE TSUNAMIS: A CASE STUDY FOR THE GREAT1852 BANDA ARC MEGA-THRUST EARTHQUAKE ANDTSUNAMI
H. RINGER, J. P. WHITEHEAD, J. KROMETIS, R. A. HARRIS, N. GLATT-HOLTZ, S.GIDDENS, C. ASHCRAFT, G. CARVER, A. ROBERTSON, M. HARWARD, J. FULLWOOD,K. LIGHTHEART, R. HILTON, A. AVERY, C. KESLER, M. MORRISE, AND M. H. KLEIN
Abstract.
We demonstrate the efficacy of a Bayesian statistical inversionframework for reconstructing the likely characteristics of large pre-instrumentationearthquakes from historical records of tsunami observations. Our frameworkis designed and implemented for the estimation of the location and magnitudeof seismic events from anecdotal accounts of tsunamis including shoreline wavearrival times, heights, and inundation lengths over a variety of spatially sep-arated observation locations. As an initial test case we use our framework toreconstruct the great 1852 earthquake and tsunami of eastern Indonesia. Rely-ing on the assumption that these observations were produced by a subductingthrust event, the posterior distribution indicates that the observables werethe result of a massive mega-thrust event with magnitude near 8.8 Mw and alikely rupture zone in the north-eastern Banda arc. The distribution of pre-dicted epicentral locations overlaps with the largest major seismic gap in theregion as indicated by instrumentally recorded seismic events. These resultsprovide a geologic and seismic context for hazard risk assessment in coastalcommunities experiencing growing population and urbanization in Indonesia.In addition, the methodology demonstrated here highlights the potential forapplying a Bayesian approach to enhance understanding of the seismic historyof other subduction zones around the world. Introduction
Indonesia is one of the most tectonically active and densely populated places onEarth. It is surrounded by subduction zones that accommodate the convergence ofthree of Earth’s largest plates. Some of the largest earthquakes, tsunamis and vol-canic eruptions known in world history happened in Indonesia [44, 27]. Since theseevents, population and urbanization has increased exponentially in areas formerlydestroyed by past geophysical hazards. Recurrence of some of these large eventsduring the past two decades have claimed a quarter million lives [2].Most casualties from natural disasters in Indonesia are caused by tsunamis,which, over the past 400 years, occur on average every 3 years (e.g. [23]). Manypotential tsunami source areas, such as the eastern Sunda [49] and Banda [26]subduction zones have no recorded mega-thrust earthquakes [53]. However, somehistorical accounts of earthquakes and tsunamis in Indonesia provide enough de-tail about wave arrival times and wave heights from multiple locations to verify ifmega-thrust events have happened in apparently quiescent regions, and assess the a r X i v : . [ phy s i c s . g e o - ph ] S e p H. RINGER ET. AL. potential consequence of a similar event occurring in the future. Indeed relianceon modern instrumental records of earthquake events to determine seismic risk se-verely biases hazard assessments, as the relevant temporal scales are hundreds orthousands of years on a given fault zone. To improve risk estimates, it is imperativeto draw from historical records of damaging earthquakes, which reach beyond thefifty to seventy year horizon provided by modern instrumentation.To this end, there has been substantial effort invested in the quantification ofthe characteristics of pre-instrumental earthquakes and tsunamis; see e.g. [49, 61,46, 45, 47, 31, 48, 8, 9, 21, 56, 39, 27, 15, 20, 41]. As noted in these references, thehistorical and prehistorical data sources are sparse in details and laced with highlevels of uncertainty. To improve the usage of these imprecise accounts, we developa systematic framework to estimate earthquake parameters along with quantita-tive bounds on the uncertainty of these parameter estimates. We do this usinga Bayesian statistical inversion approach already leveraged in a variety of disci-plines in the physical, social and engineering sciences, (see [65, 34, 13] as well as[40, 16, 62, 63]), to reconstruct large seismic events from historical accounts of theresulting tsunamis.Our focus here is on an initial case study concerning the reconstruction of the1852 Banda arc earthquake and tsunami in Indonesia detailed in the recently trans-lated Wichmann catalog of earthquakes [27, 71] and from contemporary newspaperaccounts [64]. To proceed with the Bayesian description of this inverse problem,we describe uncertainties in the noisy anecdotal observations of the 1852 tsunamivia probability distributions. We next supplement this historical data with a priorprobability distribution for the seismic parameters calibrated using modern instru-mental seismic data. Finally, we develop a forward model mapping seismic param-eters to shoreline observations using the Geoclaw software package [36, 37, 19, 5] tonumerically integrate the shallow water equations, predicting the evolution of thetsunami initiated by seafloor deformation due to the earthquake itself. These threeelements are then combined with Bayes theorem to produce a posterior distributionon the location, magnitude, and geometry of the most likely mega-thrust source forthe 1852 tsunami.Detailed information concerning the Bayesian posterior distribution, the outputof our framework, is drawn from large scale computational simulations using Markovchain Monte Carlo sampling techniques [38]. Note that while the structure of theprior distribution or the interpretation of the data may be a subject of debate, thesolution of the inverse problem detailed here is reproducible from the describedassumptions. Any of the assumptions can be modified by changing a few lines ofcode using a Python based software package available to the public upon request,and accessible via GitHub: https://github.com/jpw37/tsunamibayes . Allowingfor a mega-thrust event as justified below, we find that the most likely cause of theobservations for this 1852 account was an event with magnitude near . Mw andepicenter to the south east of the island of Maluku.The rest of this article is organized as follows. The next section includes a reviewof previous efforts related to this specific historical event, and a discussion of thesource of the 1852 tsunami (mega-thrust earthquake or submarine slump). Section3 describes the tectonic setting of the region in consideration. Section 4 givesa very brief overview of the Bayesian methodology, a description of the differentassumptions and parameterizations used for this particular event as well as an
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 3 overview of the relevant historical observations and the forward tsunami modelused here. Section 5 discusses the results of the inference and describes in somedetail the posterior distribution that yields information on the possible earthquakesthat may have resulted in the observed tsunami. Finally Section 6 discusses theimplications that can be derived from the posterior distribution, and a discussionof future work.2.
Available observational data, and previous modeling for the 1852Banda arc event
For the Banda arc in particular, we note that there were two major tsunamispossibly connected to mega-thrust earthquakes in eastern Indonesia, witnessed in1629 and 1852 [70]. Numerical models of these events [39, 15] show they werelikely sourced from the shallow subduction interface of the Banda arc subductionzone [26]. This subduction zone is largely ignored as a potential source for mega-thrust earthquakes because it includes the large, continent asperity of northernAustralia [29, 12]. However, the Banda subduction zone is similar in many waysto the western Sunda arc near Sumatra, which has sourced several mega-thrustearthquakes documented in instrumental, historical and geological records. One ofthe largest mega-thrust earthquakes in modern history happened along the northernSunda arc in 2004. Geological records of tsunamigenic events from this regionindicate average repeat times over the past 7400 years of around 450 years [59].This statement is somewhat misleading though, as the temporal interval betweenevents ranges from more than 2,000 years down to half a century. Such a variancein the temporal scales of seismic activity for a single region indicates that a lackof instrumental or even historical accounts of mega-thrust earthquakes on a givensegment of a subduction zone should not be interpreted as unlikely, but rather asinevitable. [42] argues that no low angle convergent plate boundary with over 20mm/a of convergence should be ruled out for producing mega-thrust earthquakes.Plate convergence across the Tanimbar and Seram Troughs varies from 30-70 mm/a[7].The only giant earthquake recorded instrumentally near the Tanimbar and SeramTroughs is a Mw = > m [3]. Apparently the earthquake also did not cause a landslide induced tsunamieven though there was intense shaking for several minutes. On the other hand,historical accounts of the 1852 event are much more characteristic of a mega-thrustearthquake. Like the 1938 event, it was widely felt with an estimated MMI IIIminimum diameter > > H. RINGER ET. AL.
First, we note that it is well established that landslide induced tsunamis, whilelocally more forceful and with much higher run-up wave heights, typically havea significantly shorter wave-length from the initial seafloor disturbance and hencedissipate much more rapidly than a seismically induced tsunami. With this inmind, [54] use experimental data to develop a quantitative heuristic to assist indetermining the relative likelihood of a given tsunami being induced by a landslideor directly by the earthquake. They define I as the ratio between the maximumrun-up wave height, and the horizontal extent of the wave run-up, and show that I < − for tsunamis induced by a pure seismic event, and I > − for tsunamisinduced by a landslide. Although [54] assumes a single, straight shoreline per-pendicular to the wave, the heuristic is shown to be remarkably accurate even formore realistic geometries. For the current setting, even if we presume a maximalrun-up wave height of m based on the historical account at Banda Neira, then I = m km ≈ . × − which is right in line with the values calculated by [54] forseismically induced tsunamis, but significantly less than that for landslide inducedones. In summary the lateral reach of the tsunami itself, evident from the historicalaccount for the 1852 event, is far too broad to warrant a landslide induced source.Second, [12] rightly emphasize that the eastward convergence of a thrust earth-quake along the west-dipping subduction interface would produce a negative wavein the Banda Sea prior to the arrival of the first positive wave, something that isnot recorded in any of the historical records, but shown in numerical models [15]including those utilized below. While this argument is certainly accurate, we notethat the absence of such information in the historical account does not necessar-ily preclude the occurrence of such a negative wave. The morning of the event inquestion was a spring tide which caused an extremely low tide level so that a smallnegative wave would be far less noticeable. In addition, in the 1850’s the BandaSea region observed Moluccan time, which in November would indicate a sunriseat approximately 7:30AM, only minutes prior to the earthquake. It is very likelythat the Dutch officers (the primary source of the historical record) would not takenotice of a weak negative wave so near sunrise at a spring low-tide, particularlyfollowing a devastating earthquake.In addition to these two primary reasons, we note that the 1938 Mw = 8.4earthquake [53] was very nearly in the same location as the proposed landslidesource in [12], yet there is no evidence of a submarine slump occurring as a result.This makes it less likely that the earthquake source proposed in [12] may haveinduced a slump sufficient to yield the recorded tsunami.These arguments do not eliminate the possibility that the tsunami in questionwas caused by submarine slumping along the edge of the Weber Deep as proposed by[12], but it does indicate that the potential for a mega-thrust event as the primarysource is more likely and therefore can not be disregarded. This paper addresses themega-thrust hypothesis directly and systematically shows that a mega-thrust sourcealong the Tanimbar and Seram troughs can produce a tsunami that matches thehistorical account. In addition, the systematic approach taken here clearly describesthe most likely location, strength and geometric layout of the 1852 earthquake if itwas indeed a mega-thrust event. ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 5 Tectonic Setting
The Banda arc subduction zone is the eastward extension of the Sunda arc whereboth arcs rise above where the Australian Plate subducts beneath the EurasianPlate. The composition of the lower and upper plate reverses at the boundarybetween the Sunda and Banda arcs: the lower plate changes eastward from oceanicto continental and the reverse happens to the upper plate. Both arcs are partialsubduction zones, and partial collision zones. The western Sunda arc is a collisionbetween deep sea fan deposits riding on an oceanic lower plate and a continentalarc upper plate, while the Banda arc is a collision between passive margin depositsof NW Australia and an oceanic arc [22]. At both subduction zones the Australianplate moves NNE relative to the Eurasian plate. However, the rates of motionacross the subduction zone in the Banda arc (70 mm/a, [50]) are nearly twice thoseof the western Sunda arc [7].Like the western Sunda arc, the Banda arc consists of two chains of mountains, aninner volcanic arc and an offshore chain of rising islands associated with offscrapingand accretion of thick layers of mostly sedimentary rock riding on the subductingplate [22]. In the Banda arc the accreted layers are part of the distal Australianpassive continental margin [11, 25], which according to multiple sources of data,has subducted to at least a depth of 300 km [26, 66]. The Banda volcanic arc isstill active though it is largely contaminated by subducted Australian continentalcrust [69, 30]. What was a deep trench before the buoyant continental lithospherearrived at the plate boundary is now a ‘trough’ known from SW to NE as the Timor,Tanimbar and Seram troughs, see Figure 1. Both subduction zone interfaces areactive, and rupture through sedimentary layers that drape over the subductioninterface. Earthquakes along this interface yield fault-plane solutions with low-angle thrusts and mostly dip-slip slip vectors. Shortening throughout the islandsrising above these active thrusts also verge perpendicular to the plate boundary[26].Another similarity between the Sumatra section of the Sunda arc and the Bandaarc is that each changes considerably in strike, which causes highly oblique plateconvergence in some areas. Notwithstanding this obliquity, several mega-thrustearthquakes are recorded in the western Sunda arc. Slip rakes of these events aremostly parallel to one another, but largely perpendicular to the trench at the epi-center, which is nearly 90 degrees from the Australian plate convergence direction.We challenge the idea based on geological, geophysical and historical data thatthe eastern Sunda arc and Banda arc can be dismissed as potential sources of mega-thrust earthquakes and tsunamis. In addition, the primary results of this study,illustrated in Figure 9 and Figure 10, imply that the 1852 event was most probablylocated along a narrow region in the eastern portion of the Banda arc and was amassive mega-thrust earthquake on the same scale as the December 2004 Sumatraevent. 4.
Methodology
This section lays out the complete details of the methodology we use to developour Bayesian statistical model of the 1852 event. We start with an overview of thecomponents of the model.
H. RINGER ET. AL.
Figure 1.
The seismic and geologic setting for this event. Theconvergent plate boundaries are in red and the transvergent bound-aries are in black. The arrows indicate motion of the Pacific andAustralian plates relative to the Eurasian plate. The Australianplate is converging at a rate of 70mm/yr and the Pacific plate isconverging at 110mm/yr. The nine observation locations from theWichmann catalog for the 1852 Banda arc earthquake and tsunamiare also labeled. The green rectangle indicates the region that isdepicted in Figure 9 and Figure 13. The yellow rectangle is theregion depicted in Figure 11.4.1.
Overview of the Bayesian Statistical Model.
The Bayesian approach tostatistical inversion, cf. [17, 34, 13, 65], provides a methodology for convertinguncertain outputs of a physical model into probabilistic estimates of model param-eters. This framework is perfectly suited for the anecdotal, uncertain nature of thehistorical accounts utilized here. Using this Bayesian methodology provides a lot offlexibility and will be adjusted to treat a variety of other tsunami producing seismicevents from historical and pre-historical data sets in future studies.The primary inputs required for Bayesian statistical inversion, particularly ap-plied to the determination of historical earthquakes are:(1) A prior probability distribution describing the best guess of a set of earth-quake parameters without considering the observations. For the 1852 event,we formulate the prior distribution via independent distributions on thedepth (and hence location) and the magnitude (and hence length and widthof the rupture zone). Further details as well as a description of the geomet-ric layout of the ruptures are provided in Section 4.2.
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 7 (2) A likelihood probability distribution describing measured data and observa-tional uncertainties, which for the 1852 event are observations of tsunamiarrival time, height, and coastal inundation taken from the Wichmann cat-alog. The associated uncertainties are estimated from a direct textual anal-ysis combined with other information about the shoreline locations wherethe event was recorded. Section 4.3 describes the selection of this dataset.(3) A forward model describing the relationship between model parameters andobservations. For the 1852 event, we use the Geoclaw software package topropagate off-shore tsunami waves, supplemented with a heuristic modelthat maps on shore wave heights and shoreline geometry to inundationlength. See Section 4.4 for additional details.With these inputs specified, Bayes’ Theorem gives the posterior probability dis-tribution , which describes in probabilistic terms the seismic parameters that bestmatch both our understanding of reasonable parameter values based on the priordistribution and observations extracted from historical records associated with the1852 Banda arc event. In order to extract quantitative information from our model,including a variety of marginal distributions correlating variables of interest, wemake use of Markov chain Monte Carlo (MCMC) statistical sampling techniques[34, 38]. See Section 4.5 for details.4.2.
Calibrating the parameter space and the prior distribution.
To makeefficient use of Bayesian methods, it is necessary to consider the dimensionality ofthe parameter space. As the number of parameters to be estimated increases, sodoes the difficulty of the sampling problem. This ‘curse of dimensionality’ appearsin this setting because Bayesian inference boils down to the computation of highdimensional integrals. It is known that random walk MCMC methods convergearbitrarily slowly as the dimension of the parameter space increases [57, 58, 6].Thus, we need to ensure that the dimension of the parameter space we use todescribe the earthquakes remains relatively ‘small’.A zeroth order approach to parameterizing an earthquake, is to consider the9-dimensional parameter space for the Okada model [51, 52]. However, it is un-reasonable to model the source earthquake as a single rectangular rupture alongthe Banda arc where the strike changes quite rapidly with latitude and longitude.Naively, an alternative would be to creat an N -subfault rupture zone made up of N rectangular subfaults that follow the geometry of the subduction zone, allowingfor realistic changes in the strike. Such a model would require a N -dimensionalparameter space, which produces an intractable sampling problem for any usefulvalue of N (to capture the curvature of the Banda arc, we need a minimum of N ≥ ).To reduce the dimensionality of the parameter space, we make a distinctionbetween the (forward) model parameters and the sample parameters . The modelparameters are the direct inputs to the forward model: in this case, the N Okadaparameters for an N -subfault rupture. On the other hand the sample parametersare a reduced representation of the model parameters. As such this lower dimen-sional subset of sample parameters is where we define our posterior distributionand therefore determines the dimension of the space where the MCMC is carriedout.For reducing the number of sample parameters, a good starting point is to con-sider model parameters that may be assumed to take constant values. Among the H. RINGER ET. AL. nine Okada parameters, the rake angle can be reasonably fixed to 90 ◦ . This corre-sponds to pure thrust motion, which acts perpendicular to the strike of the fault.While strike-slip motion is certainly present in real mega-thrust earthquakes, thrustmotion is the primary driver of seafloor deformation, and thus tsunami formation.Within the Okada model, rake angles other than θ = 90 ◦ are roughly equivalent toa reduction in slip by a factor of sin( θ ) . This coupling between the slip and rakeangle make it nearly impossible to distinguish the effect of one over the other inthe inference, so we leave the rake fixed throughout, and infer the slip instead.Another avenue to reduce the dimensionality is to seek model parameters thatcan be determined from other model parameters in the context of prior informa-tion. For the 1852 Banda arc earthquake, a detailed model of the subduction zonegeometry is available from the USGS Slab2 dataset [28]. The Slab2 data for theBanda arc is depicted in Figure 2. Depth, dip angle, and strike angle can be deter-mined from latitude and longitude, although each of these parameters has inherentuncertainties that are associated with each variable. To demonstrate this, we alsoshow the depth uncertainty in the upper right plot of Figure 2.Fixing the rake angle, and determining depth, strike and dip from latitude-longitude, we are left with five of the Okada parameters: latitude, longitude, length,width, and slip. These could be chosen as the sample parameter space. However, aproblem arises in choosing the triple of (length, width, slip) as sample parameters,due to their relationship with earthquake magnitude. The scalar seismic moment M of an earthquake of length L , width W , and average slip S is defined as(1) M = µLW S, where µ is the shear modulus of the rock (or Earth’s crust), with dimensions of forceper unit area. The scalar seismic moment was introduced by H. Kanamori in hisdefinition of moment magnitude M w [24]. Moment magnitude is an improvementover the classical Richter magnitude scale, and is now the standard magnitude scaleused by the U.S. Geological Survey [4]. Moment magnitude is defined as(2) M w = 23 (log M − . . It is clear that the empirical frequency of earthquakes of a given magnitude followsan exponential distribution [32], i.e. smaller earthquakes are exponentially morelikely to occur than large magnitude earthquakes. In order to ensure that mag-nitude follows an exponential prior distribution, we remove slip from the sampleparameters and replace it with moment magnitude. Given values of magnitude,length, and width, slip can be back-calculated via Equations 1 and 2.Equations 1 and 2 also highlight a challenge when using a random walk pro-posal kernel (something we are restricted to because the forward model is far tocomplicated for any type of gradient based method and the non-Gaussianity of theprior eliminates all other options) with these parameters. Since magnitude growswith the logarithm of length and width, any fixed choice of variance for length andwidth in the Gaussian proposal kernel will be inappropriate for all but a limitedrange of magnitudes. Therefore, we introduce magnitude-normalized substitutesfor length and width as sample parameters. Using the Wells-Coppersmith dataset[68] (augmented with additional data collected for more recent earthquakes), wecomputed linear least squares fits for log L and log W against magnitude. Thesefits are displayed in Figure 3. Our magnitude-normalized substitutes for lengthand width are ∆ log L and ∆ log W , the “residuals” compared to the linear best fit. ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 9
Depth (km)
Depth uncertainty (km)
127 128 129 130 131 132 1339876543
Strike angle (degrees)
Dip angle (degrees)
Figure 2.
Slab2 depth, depth uncertainty (as provided in theSlab2 dataset), strike angle, and dip angle maps for the Banda arcsubduction interface. Although there is some uncertainty associ-ated with the strike and dip angles, it is relatively minor comparedwith the relative uncertainty in depth.In other words: given values for M w , ∆ log L , and ∆ log W , length and width arecomputed as: log L = aM w + b + ∆ log L log W = cM w + d + ∆ log W where a, b, c, d are the coefficients of the linear best fits as shown in Figure 3. Oncethe magnitude M w , and length and width are found, then the slip can be computedfrom (1) and (2).To the five sample parameters (latitude, longitude, magnitude, ∆ log L , ∆ log W ),we add a sixth parameter: depth offset . The Slab2 data includes estimates ofuncertainty in the subduction interface depth (see Figure 2). Depth offset accountsfor this uncertainty by allowing for earthquakes that are situated somewhat deeperor shallower than is specified in the Slab2 depth map. By the above approach, we Length vs magnitude
Best fit+/ L Wells-Coppersmith 6 7 8 9 1010 Width vs magnitude
Best fit+/ W Wells-Coppersmith
Figure 3.
Wells-Coppersmith and modern data plotted alongsidelinear best fits for log L and log W against M w .reduce dozens or hundreds of Okada parameters to a set of six – latitude, longitude,magnitude, ∆ log L , ∆ log W , and depth offset – that is both low enough in dimensionto be computationally tractable and sufficiently independent to well-represent thepossible earthquakes, as we describe next.4.2.1. Computing Subfault Model Parameters.
As mentioned above, it is necessaryto model the earthquake as a collection of rectangular subfaults that conform tothe subduction interface geometry. Here we describe our approach for “decompress-ing” the six sample parameters introduced above into the Okada parameters for Nrectangular subfaults that follow the interface geometry provided by Slab2.The basic approach is to “break” a single rectangular rupture zone into an m × n grid of identical subrectangles, which are then oriented to conform to the interfacegeometry. Each of these subrectangles has length L/m and width
W/n , where L and W are the length and width of the full rupture zone. The difficulty thenlies in choosing the number of subfaults, and identifying the dip angle and strikeorientation for each one so as to best match the fault geometry.For ease of implementation, we use odd values of m and n . We found that m = 11 and n = 3 made an appropriate fit for the 1852 Banda arc earthquake, adequatelycapturing the geometry without over-fitting the parameter space. To determine theorientation and location of each subfault, we place a single point at the latitudeand longitude of the centroid of the full rupture zone. Using the Slab2 map ofstrike angle, we move in opposite directions, staying parallel to strike. Every L/m kilometers, we place another point. This continues until m points are placed. Foreach point, we then move in opposite directions, perpendicular to the strike angle,placing points every W/n kilometers, until all mn points have been placed. Thesepoints are the latitude/longitude coordinates for the centers of the subrectangles.This procedure is displayed in Figure 4.Having specified the latitude, longtitude, length, and width for each subrectan-gle, the remaining Okada parameters are determined as follows. Each subrectangleis given the same slip value as determined by (1) and (2). The strike and dip angles ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 11
130 131 132654
Rupture zone center
130 131 132654
Strikeward centers
130 131 132654
All centers
Fault planes
Figure 4.
Placing subrectangles contoured to interface geometry.First, a point is placed at the center of the rupture zone. Pointsare then placed forwards and backwards following the strike angle(essentially following level curves of depth). Additional points areplaced up-dip and down-dip. Using Slab2 depth, dip, and strikedata, Okada parameters for rectangles centered at each point arecomputed.are determined by the Slab2 strike and dip maps using the epicenter of each sub-rectangle as the reference point. The depth is determined by the Slab2 depth map,plus the value of the depth offset sample parameter (universal to all subfaults). Asdiscussed above, all subfaults are assigned a rake angle of 90 ◦ .4.2.2. Specifying the full prior distributions.
Selection of appropriate prior distribu-tions is a key step in good Bayesian inference. An over-specified prior can overwhelmthe data, and an under-specified prior may allow for parameter values that are non-physical. Having discussed the map from sample parameters to model parameters,we now discuss our choice of prior distributions for latitude, longitude, magnitude, ∆ log L , ∆ log W , and depth offset for the 1852 Banda arc event.Prior constraints on earthquake latitude and longitude are derived from the sub-duction interface geometry. Large earthquakes can only be supported in a certainrange of depth: too deep, and the crust is too plastic to store the strain energynecessary for a large earthquake [60], to shallow, and the rupture interface wouldextend above the surface. We take the approach that, a priori, depth is the primary Table 1.
Prior distributions for all sample parameters of the 1852Banda arc earthquake
Parameter name(s) Kind Distribution ParametersLatitude & longitude pre-image of truncated nor-mal via depth (in km) µ = 30 , σ = 5 , ( a, b ) =(2 . , λ = . , ( a, b ) = (6.5,9.5)∆ log L normal µ = 0 , σ = . W normal µ = 0 , σ = . µ = 0 , σ = 5 constraint on earthquake location. Since the Slab2 dataset gives a depth map for thesubduction zone along the entire Banda arc, any probability distribution on depthproduces an implied distribution on latitude and longitude, at least for mega-thrustevents like those considered here. Based on the augmented Wells-Coppersmithdataset, we chose a truncated normal distribution for depth. This distribution issupported on [2 . , kilometers, with a mean of 30km and a standard deviation of5km. Evalutating the pdf of this distribution at each latitude/longtitude coordi-nate, via the Slab2 depth map, gives a non-negative continuous function. Althoughthis function does not integrate to unity, the normalizing constant cancels out in theevaluation of the Metropolis-Hastings acceptance parameter α , i.e. this functionprovides an adequate weight for the latitude and longitude sample parameters. Theunnormalized logpdf of the latitude/longitude prior is displayed in the left panel ofFigure 9.As discussed above, earthquake magnitude is observed to approximately followan exponential distribution. Clearly however, the exponential scaling cannot con-tinue indefinitely in the large magnitude regime, and a number of approaches havebeen used to address this (see [32]). We take the simple approach of right-truncatingthe exponential distribution at magnitude 9.5, i.e. we do not allow for an eventof magnitude greater than . Mw. A consensus estimate for the parameter of theexponential distribution is λ = . [32], which we use here.Since ∆ log L and ∆ log W are magnitude-normalized length and width, definedas residuals against a linear best-fit, we chose Gaussian prior distributions withmean zero for each of these sample parameters. The standard deviations for thesedistributions are determined from the sample variances for the residuals in theaugmented Wells-Coppersmith dataset against the linear fit (see Figure 3). Thesevalues are σ ∆ log L = 0 . and σ ∆ log W = 0 . .The prior for depth offset was chosen based on the Slab2 depth uncertaintydata. The average reported uncertainty is roughly 5km, so a mean-zero normaldistribution with standard deviation of 5 was selected.4.3. The Historical Dataset, the relevant uncertainties, and assigned like-lihood distributions.
Overview of historical account and potential observations.
Observations areselected from the historical accounts in the Wichmann catalog [70, 71] based ontwo key criteria. First, the account has to provide an identifiable location (latitude-longitude) that can be incorporated into the modeling. In other words, the details
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 13
Table 2.
Likelihood distributions for all tsunami observations ofthe 1852 Banda arc earthquake. Height and inundation are inmeters, and arrival is in minutes.
Observation location & type Kind Distribution ParametersPaulau Ai (height) normal µ = 3 , σ = 0 . µ = 1 . , σ = 0 . µ = 15 , σ = 5 , a =2Banda Neira (height) normal µ = 6 . , σ = 1 . µ = 185 , σ = 65Pulau Buru (height) χ µ = 0 . , σ = 1 . , df = 1 . χ µ = 0 . , σ = 2 , df = 1 . µ = 45 , σ = 5Saparua (height) normal µ = 5 , σ = 1Saparua (inundation) normal µ = 125 , σ = 40Kulur (height) normal µ = 3 , σ = 1Ameth (height) normal µ = 3 , σ = 1Amahai (height) normal µ = 3 . , σ = 1 provided in the historical account must be sufficiently accurate to yield a preciselocation via modern-day maps and information. Second, the account has to besufficiently detailed that some level of confidence can be placed on the observablein question. Note that drawing from a catalog of this kind introduces unavoidableambiguities that do not apply to modern instrumental data. For example, wespecify the wave height based on passages of the form “[t]he water rose to the roofsof the storehouses and homes,” as described in more detail below.Thirteen different observations for the 1852 Banda arc tsunami meet these crite-ria spread across nine locations, which are shown in Figure 1. These include threetypes of observations:(1) Arrival time.
The arrival of the first significant wave after the shakingstopped. We assume that the arrival time refers to the first wave, not themaximal one.(2)
Maximum wave height.
This is the most frequent observable, and is identi-fied at every location.(3)
Inundation length.
This refers to the distance inland that the wave traveledonshore, and is actually interpreted for our purposes as a deterministicfunction of the wave height. This essentially places a double amount ofweight on those locations that have observations of both wave height andinundation.Based on the text of each account, a probability distribution is developed de-scribing the probabilistic likelihood that each observation took a given value. Thesedistributions, which are assumed to be independent, are shown in Figure 5 (exactspecification as given in Table 2). Rather than explain the reasoning behind allthirteen of these likelihood distributions for each of the nine locations, we onlyprovide a detailed discussion of the likelihood for a single location: Banda Neira.
Banda Neira: a sample likelihood distribution.
Observations at this locationare primarily taken from page 242 in the Wichmann catalog which states in part:“Barely had the ground been calm for a quarter of an hour when the flood wavecrashed in...The water rose to the roofs of the storehouses and homes...[the wave]reached the base of the hill on which Fort Belgica is built on Banda Neira”. Weexpect the wave height observation to be near the boat dock on Banda Neira which isjust east of Fort Nassau. For the available bathymetry data, we seek a location nearthat point that will maintain a sizable wave for a reasonably initiated tsunami. Withthis in mind, we select − . ◦ latitude and . ◦ longitude as the observationlocation.Using 15 minutes as the anticipated arrival time of the wave at Banda Neira istoo simplistic for these circumstances. In particular it is noted in other locationsthat the shaking lasted for at least 5 minutes, but the modified Okada model used inGeoclaw here assumes an instantaneous rupture. Hence we build into the likelihood,a skew toward longer times with a mean of 15 minutes. This is done with a skew-normal distribution with a mean of 15 minutes, standard deviation of 5 minutes,and skew parameter 2.Assuming standard construction for the time period for the homes (and store-houses) we can assume the water rose at least 4 meters above standard flood levelsas most buildings of the time were built on stilts and had steep vaulted roofs. Basedon the regular storm activity in the region we can expect that with high tide, andnormal seasonal storm surge, the standard flood level is also approximately 2 me-ters in this region. This leads us to select a normally distributed likelihood for waveheight with a mean of . m and standard deviation of . m , allowing for reasonablelikelihood values for wave heights in the range from m to m .To quantify the wave reaching the base of the hill, we measured the distance from20 randomly selected points along the beach to the edge of said hill in ARCGIS. Themean of these measurements was 185 meters, with a standard deviation of roughly65 meters. Thus we choose a normal distribution with those parameters. Withoutmore detailed information about the coastline, and a direct idea of the direction ofthe traveling wave, we can not be more precise with regard to the inundation, butthis is sufficient for the model we use (as described below).4.3.3. Overview of all likelihoods.
The likelihood distributions for the other 8 loca-tions are constructed in a very similar manner to that described above for BandaNeira with exact parameters given in Table 2 and shown graphically in Figure 5.The total likelihood of a given event is then computed as the product of theseindividual observational likelihoods (we rely heavily on the assumption that eachobservable is independent of the others). The assumption of independence of thedifferent observations is certainly questionable, but there is also no reason to sup-pose that a more complicated construction of the total likelihood is preferable, i.e.we have chosen to take the most simplified approach without making additionalunjustifiable assumptions about the structure of the likelihood.
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 15
First wave arrival time (minutes)
Maximum shoreline wave height (m)
Maximum inundation length (m)
Pulu AiAmbonBanda NeiraBuruHulaliuSaparuaKulurAmethAmahai
Figure 5.
The Forward Model.
Calculation of the forward model that maps seismicevents to quantitative observables is the most complicated and computationallyexpensive part of the inversion process. We compute the tsunami observationsresulting from seismic events by numerically integrating the shallow water equationsin a restricted region surrounding the Banda Sea.4.4.1.
Geoclaw Integration and ‘high’ resolution bathymetry.
The propagation ofthe tsunami waves is computed via the nonlinear shallow water equations supple-mented with the appropriate initial and boundary conditions dictated by the spec-ified Okada parameters and bathymetry of the region. We simulate the tsunamigenerated by each Monte Carlo sample using the Geoclaw software package, [36,37, 19, 5] which employs an adaptively-generated mesh for a finite volume basedscheme. For bathymetry (sea-floor topography) we use the 1-arcminute etopodatasets available from the open access NOAA database referred to hereafter asNOAA bathymetry, and for the coastline near each observational point we utilizehigher resolution Digital Elevation Models (DEM) from the Consortiom for SpatialInformation (CGIAR-CSI) referred to below as DEM coastlines. These higherresolution topographical files yield a 3-arcsecond resolution on land, but give no http://srtm.csi.cgiar.org/srtmdata/ Table 3.
Specification of the extent of the DEM coastline filesused near each of the historically observed accounts.Observation location Latitude extent Longitude extentBanda Neira & Pulau Ai [ − . , − . . , . Ambon, Saparua,Haruku, & Nusa Laut [ − . , − . . , . Pulau Buru [ − . , − . . , . Amahai [ − . , − . . , . additional information on the sub-surface bathymetry. The extent of each of thesefiles is provided in Table 3.In addition to these DEM coastline datasets and the NOAA bathymetry, wealso took advantage of detailed sounding maps available from the Badan NasionalPenanggulangan Bencana (BNPB or Indonesian National Agency of Disaster Coun-termeasure, see http://inarisk.bnpb.go.id ). To convert this data into digitallyaccessible information, contours were taken from images exported from the websiteand then traced and interpolated in arcGIS to produce approximate depths in thesame regions as specified in Table 3. For example, the bathymetric readings basedon this data are shown in Figure 6 for the bay of Amahai. The upper left panelin Figure 6 depicts the bathymetry data that is gleaned from the BNPB and digi-tized by interpolating across contours of constant depth in arcGIS. The upper rightpanel of Figure 6 depicts the bathymetry/topography from the NOAA bathymetrydataset. Using the built in interpolative methods in Geoclaw’s topotools package( topotools.interp unstructured with the cubic interpolant, and a proximity radius of1000), we interpolate the coastline and coarse bathymetry from the NOAA datasetto match the bathymetric contours from the upper right panel to produce the lowerleft panel. This lower left panel does not accurately capture any of the topographi-cal features of the coastline and suffers significantly from interpolant error onshoreas there are no bathymetric readings there. The actual shoreline and onshore to-pography is then overlaid from the DEM coastlines on top of the bottom left panelof Figure 6 to create the final product which is seen in the bottom right panel ofthe same Figure. This retains the improved bathymetric contours, and yields anaccurate coastline and near-shore topographical profile.This same process is repeated for Palau Buru, and the coastline near the islandsof Ambon, Saparua, Haruku, and Nasu Luat. The resultant final bathymetric filesare shown in Figure 7. Finally, all of these high resolution bathymetric files areused by Geoclaw when the wave approaches these locations onshore.For the region near Banda Neira and Palau Ai, the bathymetric data was stillquite rough, particularly for the narrow channels between Banda Neira, Banda Api,and Lonthor. We obtained a set of soundings for this region from a map publishedby the Kepala Dinas Hidro-Oseanografi (the Indonesian Navy Hydrography andOceanography Center) from data collected primarily in 1928/1929 [1]. Using thesame approach as described above for the bay of Amahai, these discrete soundingsare interpolated for the entire region surrounding the Banda islands (except that alinear interpolant is used instead of cubic due to the sparsity of the measurements) ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 17
Figure 6.
Combining all of the bathymetric and topographicalsources into a single file for the bay near Amahai. The upper leftfigure demonstrates the bathymetry drawn from the level curvesexported from http://inarisk.bnpb.go.id . The upper right fig-ure shows the level of resolution for the NOAA bathymetry data.The lower left figure shows the interpolation of these two datasets (omitting the interior of the coast, i.e. all grid points formthe NOAA bathymetry that are not below sea level, or border agrid below sea level). The lower right figure is the final product,combining the improved bathymetric data with the DEM coastlinedataset.and overlayed with the DEM coastlines. The resultant bathymetry files for PalauAi and the Banda islands are shown in Figure 8.For the forward simulations of the tsunami wave, we employ an adjoint-basedadaptive mesh strategy [14]. This entails solving a linearized adjoint equationbackward in time with sources centered at each gauge location The solution of theadjoint equation produces waves that propagate backward in time from the desiredobservation locations to indicate what part of the forward wave will eventuallyinfluence the tsunami at those locations (see [14] for details). To initialize theadjoint solver, we place a smoothed Gaussian perturbation h ( x, y ) to the waveheight at each gauge location given by:(3) h ( x, y ) = (cid:88) k exp( − r k / , Figure 7.
The final bathymetry/topography renderings for theremaining historical observations sites excluding Banda Neira andPalau Ai. The figure on the left includes the bay of Ambon,Haruku, and all of Palau Saparua, and the right hand figure depictsthe refined bathymetry for Kayeli Bay in Palau Buru.
Figure 8.
The refined bathymetric data for the Banda islands andPalau Ai. The upper left panel depicts the discrete sounding datathat was available. This is interpolated with the NOAA bathymet-ric data to produce the top right panel, and then overlayed with theDEM coastlines and cropped for the final bathymetry files shownon the bottom left panel (Banda islands) and bottom right (PalauAi).
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 19 where r k is the distance from the point ( x, y ) to the gauge location ( x k , y k ) . Thesolution of the linearized adjoint problem guides the choice of refinement regionsof the fully nonlinear forward model, indicating where the wave that will reach theobserved locations will be at specific times. The benefit of using this approachas noted in [14] is that only those parts of the wave that will reach the desiredlocations are refined, i.e. the mesh refinement is restricted to those parts of thedomain (in both space and time) that will most influence the final wave at thedesired location. In addition, for the application at hand, we only need to run thebackward adjoint solver once, and then the generated output can be used for everysample so long as the gauge locations are not changed. This saves a substantialamount of computational cost, allowing us to use a much finer mesh near theobservational locations than a standard adaptive mesh would have allowed.We use an adaptive mesh with 6 levels, starting with 6 arcminute resolutionin the open water with no motion, and then going through × , × , × , × , and × grid refinements to those regions where the adjoint indicates the wave will be,resulting in the finest grid of 3 arcseconds which matches the fine resolution of theDEM coastline files. This means that the mesh levels are given by 6 arcminute, 3arcminute, 1.5 arcminute, 45 arcsecond, 15 arcsecond, and 3 arcsecond resolutionrespectively. In addition to this dynamic adaptation of the mesh, we staticallyfix regions near each gauge at the highest mesh resolution (3 arcseconds) for theentirety of the simulation, thus accurately capturing the wave characteristics nearthe observed locations. These regions are explicitly specified in Table 4. Imple-mentation of such a highly refined grid for the region in question required someminor modification of the default list lengths in the fortran code as described inthe code repository. The backward adjoint solver is run on a 15 arcsecond grid andthe output files are saved every 5 minutes to ensure adequate spatial and temporalresolution for the dynamic grid refinement. Geoclaw interpolates these output filestemporally to determine the wave location throughout the entire simulation.All other settings in Geoclaw are set to their default values. An adaptive timestep is adjusted according to the Courant-Friedrichs-Lewy (CFL) condition witha desired CFL of . . The spatial discretization in Geoclaw is a second orderscheme with the MC limiter [35] employed to avoid the development of un-physicalshocks. All simulations are run for a physical time window of . hours to ensurethat the wave has reached all of the relevant locations (for this event the longesthistorically recorded time between the earthquake and the arrival of the wave wasapproximately 40-45 minutes as shown in Figure 5). Each simulation of Geoclawgenerates wave heights and arrival times at the locations shown in Figure 1.4.4.2. Wave inundation calculation.
The observations of wave inundation at BandaNeira and Saparua are very precise, and seem to be important to infer the earth-quake. Despite the precision of this measured distance, it is unclear what specificpart of the shoreline these observations were recorded for, and even if this was clear,the highest resolution topography available isn’t sufficient to completely trust a sim-ulated wave inundation. For these reasons, we opted to use a simplified model ofwave inundation that is a deterministic function of the wave height on the shoreline.The model that dictates the mapping from on-shore wave height to wave inun-dation length is taken from [10] (equation (2.15)), yielding a relationship betweenthe maximum inundation length, maximal wave height on-shore, and the average
Table 4.
Specification of the statically refined regions labeled ac-cording to the historically observed data. To maintain computa-tional tractability, some choices were necessarily made regardingwhich regions in the computational domain were needed at thehighest mesh refinement level. For instance, the narrow channelbetween the islands of Haruku and Palau Saparua is captured viatwo distinct refined grids to avoid having to much spatial refine-ment unnecessarily placed over land.Observation location Latitude extent Longitude extentBanda Neira [ − . , − .
49] [129 . , . Pulau Ai [ − . , − . . , . Ambon [ − . , − .
66] [127 . , . Hulaliu [ − . , − . . , . Pulau Buru [ − . , − .
27] [127 . , . Saparua (near port) [ − . , − . . , . Saparua (main bay) [ − . , − . . , . Nusa Laut [ − . , − . . , . Amahai [ − . , − . . , . Channel between Haruku& Saparua [ − . , − . . , . Channel between Haruku& Saparua [ − . , − .
54] [128 . , . slope of the shoreline. This is given by: I β,n ( H ) = k · H . cos( β ) n (4)where I is the maximum landward inundation distance as a function of the max-imum shoreline wave height H . Here n is Manning’s coefficient of friction [67](selected using arcGIS and Google images near the coastline in question, β is theslope of the coastline and k is an empirically determined constant equal to . .We computed β as the average slope taken from a series of 1-dimensional verticalprofiles taken from arcGIS perpendicular to the coastline near each gauge. Thismodel was used to convert on-shore wave heights into on-shore inundation whichis then evaluated against the likelihood probabilities of on-shore inundation fromFigure 5.4.5. Sampling and Convergence.
To quantify the results of the Bayesian infer-ence, Markov Chain Monte Carlo (MCMC) was used to generate samples from theposterior measure. Because we did not have an adjoint solver for this PDE-basedforward map, gradient-based methods like Hamiltonian Monte Carlo were not avail-able. We therefore employed random walk-style Metropolis-Hastings MCMC withperiodic Sequential Monte Carlo-style resampling according to posterior probability(see [13, Section 5.3]). A diagonal covariance structure was used for the proposal
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 21
Table 5.
Random walk standard deviations by parameter.Parameter Standard DeviationLatitude . Longitude . Magnitude . log L . log W . Depth offset 0.5kernel, with the step size in each of the six parameters tuned to approximate theoptimal acceptance rate of roughly . [17, Section 12.2]. The final standard devi-ations for the random walk proposal kernel are given in the table below: To ensurethat all viable seismic events were considered, we initialized fourteen (14) MCMCchains at locations around the Banda arc with initial magnitudes of either . or . Mw. Additional chains were tried at . Mw; however, these were quickly discardedas they consistently failed to generate a sufficiently large wave to reach each of theobservation points (Figure 1) and therefore produced likelihoods of zero probability.Each of the 14 remaining chains was run for 24,000 samples, including a “burn in”of 6,000 samples, with resampling according to posterior probability every 6,000samples. Our approximation to the posterior therefore is made up of a total of252,000 samples. To ensure accurate statistics, chains were run well beyond whenthey appeared to have converged; for example,
Gelman-Rubin diagnostic R [18, 17]for all parameters fell below 1.1 (a common convergence criterion) after 8,000 sam-ples. Samples were computed using the compute resources available through BYU’sOffice of Research Computing, consuming a total of nearly 200,000 core-hours inall. 5. Results
Summary.
The results of the inference are summarized in Figures 9 and 10and show stark differences between the prior and posterior distributions. Whereasthe prior encompassed all parts of the Banda arc with a reasonable depth, theposterior for the epicenter is narrowly concentrated in a region near 4.5 ◦ S, 131.5 ◦ E,which is situated in a shallow part of the subduction interface. That is, the modelingimplies that the epicenter must have been located within this region to best matchthe observations summarized in Figure 5. Also notable is the marginal posterior formagnitude: despite a prior that heavily preferred lower magnitudes, the posterioris concentrated around earthquakes of magnitude 8.8. This is because the tsunamisimulations indicated that a large event was required to produce the wave heightsdescribed in the historical accounts at several of the observation locations.More subtle inference is seen in magnitude-normalized length and width. Theposterior favors rupture zones that are relatively narrow for their magnitude, i.e.,the length and width are smaller than is typical, given the magnitude of the event.This is likely because the inference is trying to balance observations of wave heightand arrival time. Simulations indicate that an earthquake needs to be quite large in
127 128 129 130 131 132 1338765432
Prior logpdf
127 128 129 130 131 132 133
Samples
Figure 9.
Sample epicenters from the posterior compared withthe prior in latitude/longitude (the same region in the green rec-tangle of Figure 1). The posterior is concentrated in a small regionin the northeast of the arc. Note that the color gradient is not thesame quantitatively between these two plots, but the general char-acterization is accurate, i.e. the prior is distributed evenly over theentire arc (via the depth calculation) while the posterior is heavilyconcentrated in the northeastern section of the arc.
Magnitude
10 0 10 20
Depth offset log L log W
10 20 slip depth
Figure 10.
Magnitude, depth offset, ∆ log L , and ∆ log W pos-terior histograms, compared to the associated prior distributiondensities (plotted in green) The bottom two plots are histogramsof the slip and depth (in meters for both) drawn from the posterior.order to produce the observed wave heights in, for instance, Banda Neira. However,larger earthquakes, all else being equal, have rupture zones that are closer to Banda ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 23
Neira, thus reducing the arrival time of the wave. The posterior therefore favorsa smallish rupture zone given the magnitude (keeping the arrival time in check)balanced by a larger slip to achieve the observed wave heights.5.2.
Fault Characteristics by Epicenter Location.
We now describe how thefault parameters change with geographic location according to the computed pos-terior. Figure 11 displays the approximate conditional expectation for the sampleparameters magnitude, depth offset, ∆ log L , and ∆ log W , conditioned on latitudeand longitude as well as the conditional expectation for the model parameters slipand depth; these figures therefore show the expected value of each parameter if wewere to assume a given epicenter location. Several trends are apparent: • The farther outside the arc, the higher the expected magnitude.
This is not surprising, as higher magnitudes would be required to producelarge enough waves at that distance from the observation points. • The farther outside the arc, the greater the value of depth offset.
This appears to counteract the shallowing of the fault towards the outside ofthe arc, ultimately producing earthquakes at nearly constant depth amongaccepted samples. This can be seen even more clearly in the conditionalexpectation of depth, where there are variations in the depth, but they arerelatively mild (most of the points in this plot are on the shallower end ofthe colorbar). • The closer the center of the rupture is to the coast of Seram,the shorter the length (smaller ∆ log L ), and higher the slip of therupture. This is likely due to the rupture extending underneath SeramIsland, which leads to a smaller tsunami (as only some of the rupture occursbeneath the ocean). Thus, a shorter rupture zone increases the slip as seenin the conditional expectation of slip (and thus wave height), counteractingthe influence of Seram Island on the tsunami generation. On the otherhand, the earthquakes to the south or west have lower slip values as thetsunami source is then closer to the observation locations.5.3.
Forward Model and Output.
It is also worth considering the implied ob-servation distributions, known as the posterior predictive distributions . This hastwo primary interpretations, first in understanding the drivers of the results, andsecond as an assessment of hazard risk. We describe each of these in turn.
Understanding the Results.
Comparing the posterior predictive distributionsand the likelihood distributions, shown in Figure 12, helps describe what drove themodel’s conclusions. Banda Neira and Saparua provided the largest contributionto the likelihood, and we see that the posterior samples broadly matched our in-terpretation of the observations there. The posterior samples at Kulur and Amethstand out as different from the proposed likelihood, with waves smaller than ourinterpretation of the accounts. However this is acceptable and actually validatesthe current approach, given that these accounts were not specific, and we assignedwide distributions to them, indicating that the posterior distribution matched thelikelihood from the other locations but was unable to match the proposed distri-butions at these locations. Overall, the posterior distribution is consistent withthe observations recorded in our sources, with sufficient differences to make theinferred posterior distribution distinct from both the likelihood density and theprior distribution. This indicates that the posterior is not just a re-sampling of the
Magnitude
Depth offset logL log W Slip
Depth
Figure 11.
Posterior conditional expectation of sample param-eters magnitude (Mw), depth offset (km), ∆ log L , and ∆ log W as well as the slip (m) and depth (km) conditioned on lati-tude/longitude, shown for the same region as the yellow rectanglein Figure 1.likelihood and/or prior but has instead found earthquakes that best match the dataand physical limitations of the situation.Additional insights into the limitations of our approach can be gleaned from thedifferences between the posterior predictive distribution and the likelihood. Forinstance in addition to the differences for the wave height at certain locations, wenotice a significant difference between the distributions for arrival time at bothBanda Neira and Saparua. In fact, it appears that our assigned likelihood distribu-tion over estimated the actual arrival time. Such a discrepancy is not unexpectedhowever, because the current forward model assumes an instantaneous rupturewhereas the historical account implies that the earthquake lasted approximately5 minutes which would adequately account for the shifted distributions shown inFigure 12. Hazard Risk.
The posterior predictive distribution can also be interpreted asan estimate of the hazard risk for the specified observation locations. The his-tograms in Figure 12 represent what communities in these locations might be ex-pected to experience – the wave heights, arrival times, and inundations – shoulda similar event happen in the future. For instance, if an event of this magnitudeoccurred in the same location on the Banda arc, we anticipate a wave of approxi-mately . m to reach the populous city of Ambon (approximately 300,000 people).For those living in the bay of Ambon, this gives a probabilistic hazard assessmentthat can be coupled with detailed topographical information to assess potentialflood levels as well as economic and societal impact from a similar future event. ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 25
Pulu Ai height
Ambon height
Banda Neira arrival
Banda Neira height
Banda Neira inundation
Buru height
Hulaliu height
30 40 50 600.000.050.100.150.200.25
Saparua arrival
Saparua height
Saparua inundation
Kulur height
Ameth height
Amahai height
Figure 12.
Model output compared to likelihood densities (plot-ted in green). Arrival times are in minutes, wave height and inun-dation length are in meters.5.4.
Claim & Corroborating Evidence.
The implied claim of our posteriordistribution is this: if the 1852 Banda arc tsunami was caused by a mega-thrustsubduction zone earthquake, it was a magnitude ∼ near 4.5 ◦ S, 131.5 ◦ E, and this type of event matches the historical account quitewell. During the analysis, we discovered an item of corroborating evidence for thisclaim, in the form of the Slab2 depth uncertainty data. The Slab2 model of thesubduction zone is based on seismically collected instrumental data that can be usedto infer the interface geometry. The more earthquakes that have occurred recentlyon a particular segment of a fault, the more certain we can be of the geometry.Regions of uncertainty correspond to “seismic gaps”; fault segments that have beenrelatively silent during the modern period of instrumental data collection. A seismicgap may represent a location where hundreds of years of stress has accumulated,which eventually results in a large earthquake when the fault slips and the stressis released [43]. While not all seismic gaps turn out to be dangerous [33], they arestill important to consider as possible sources for an event such as the 1852 Bandaarc earthquake.
127 128 129 130 131 132 1338765432
Samples vs Banda Arc Seismicity
Figure 13.
Banda arc Slab2 seismicity compared to the locationof samples from the posterior distribution. The seismicity is com-puted using a Gaussian kernel density estimator from locationsof instrumentally recorded earthquakes (as provided in the Slab2dataset) and is displayed on a blue-yellow color gradient with thedark blue regions referring to areas where little to no seismic eventswere recorded, and yellow being the highest concentration of seis-mic happenings. The logarithm of the expected posterior distribu-tion conditioned on the latitude-longitude location is shown in theorange color scheme. Clearly the posterior distribution is concen-trated in a region of low instrumentally recorded seismicity. Thisis the same region as the green rectangle in Figure 1.
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 27
Both the Slab2 depth uncertainty, and the underlying seismic dataset, demon-strate the presence of a seismic gap in the region where our posterior distribution isconcentrated (see Figure 13). This can be viewed as evidence, distinct and separatefrom our usage of Slab2, that supports the results of our analysis.6.
Discussion and Future Work
This paper has presented a robust approach to determining the strength andlocation of historically observed mega-thrust earthquakes via observations of theresultant tsunami. The Bayesian nature of this investigation not only provides anunderstanding of the earthquake parameters that may have caused this event, butalso yields information regarding the uncertainty of these post-dicted parameters,and the correlations between them. We note that this 1852 Banda arc event iscertainly not the only event for which this methodology is feasible. There are severalother events both in Indonesia and elsewhere for which there is either historical orgeological evidence for significant seismic events, but little or no instrumental datato draw from. In addition, although this article focuses on mega-thrust events, themethodology is sufficiently robust to work for a submarine slump induced tsunami,or any other type of hazard for which a forward model is available. We will lookinto these possibilities in future studies, including analyzing the potential for the1852 event to have been caused by a submarine slump as suggested by [12, 55].We note that since the 1852 earthquake investigated here, the Banda arc regionhas seen an order of magnitude increase in population, increasing the seismic andtsunami disaster potential for the inhabitants. However, risk assessments basedon the short (geologically speaking) period of instrumental records for this andother densely populated regions (i.e. Java and Bali) underestimate the seismic andtsunami disaster potential. Some of these regions have not experienced mega-thrustearthquakes for several hundred years. At the current rates of seismic loading onthese subduction zones, enough elastic strain has accumulated to cause anotherMw 8.5-9.0 event akin to that described here in the Banda Sea. This fact, coupledwith the evidence provided here, indicate that at least in the case of the Banda arc,mega-thrust events are not only possible, but highly likely to have occurred in thepast and thus likely to recur in the future.
Acknowledgements.
Data supporting the conclusions of the work (in particularthe construction of the prior and likelihoods, as well as csv files containing all datafor the posterior) are available through the GitHub repository.The authors acknowledge Advanced Research Computing at Virginia Tech andResearch Computing at BYU for providing computational resources and technicalsupport that have contributed to the results reported within this paper. JPW andRH would like to thank the Office or Research and Creative Activities at BYUfor supporting several of the students’ efforts on this project through a MentoringEnvironment Grant, as well as generous support from the College of Physical andMathematical Sciences and the Mathematics and Geology Departments. We alsoacknowledge the visionary support of Geoscientists Without Borders. JPW wouldlike to thank Tulane University for supporting him as a visiting scholar during theSpring 2018 semester where a portion of this work was carried out. http://rc.byu.edu We would like to thank our colleagues W. F. Christensen, J. Borggaard, E.J. Evans, R. LeVeque, V. Martinez, S. McKinley, C. Mondaini, C. S. Reese, G.Simpson & F. Viens for valuable feedback on this work. In addition to the authors,there were several students at BYU that provided some feedback and assistance tovarious parts of this project including: A. (Oveson) Bagley, R. Howell, J. Lapicola,C. Carter, I. Sorenson, B. Berrett, C. Ringer, J. Voorhees, and Y. Zhao.
References Detailed soundings of the banda islands , 2011, Based on data collected by Dinaa HidrografiNegara Belanda (The Hydrographic Service of the Royal Netherlands Navy) in 1928/1929.2.
Tsunami Sources 1610 B.C. to A.D. 2017 from Earthquakes, Volcanic Eruptions,Landslides, and Other Causes , , 2017, Accessed: 2020-09-01.3. Anonymous, Aardbevingen in den Oost Indischen Archipel waargenomen gedurende het jaar1938 , Natuurk. Tijdschr. Ned. Indie (1940), 38–74.4. William Bakun, USGS Earthquake Magnitude Working Group , 2002 (accessed July 19, 2020).5. Marsha J Berger, David L George, Randall J LeVeque, and Kyle T Mandli,
The GeoClawsoftware for depth-averaged flows with adaptive refinement , Advances in Water Resources (2011), no. 9, 1195–1206.6. Alexandros Beskos, Gareth Roberts, and Andrew Stuart, Optimal scalings for localMetropolis–Hastings chains on nonproduct targets in high dimensions , The Annals of Ap-plied Probability (2009), no. 3, 863–898.7. Yehuda Bock, L Prawirodirdjo, J F Genrich, C W Stevens, R McCaffrey, C Subarya, S S OPuntodewo, and E Calais, Crustal motion in Indonesia from global positioning system mea-surements , Journal of Geophysical Research: Solid Earth (2003), no. B8.8. Stein Bondevik,
Earth science: The sands of tsunami time , Nature (2008), no. 7217,1183.9. E Bryant, Grant Walsh, and Dallas Abbott,
Cosmogenic mega-tsunami in the Australia re-gion: are they supported by Aboriginal and Maori legends? , Geological Society, London, Spe-cial Publications (2007), no. 1, 203–214.10. Edward Bryant,
Tsunami: the underrated hazard , Springer, 2014.11. David J Carter, Michael G Audley-Charles, and AJ Barber,
Stratigraphical analysis of islandarccontinental margin collision in eastern Indonesia , Journal of the Geological Society (1976), no. 2, 179–198.12. Phil R Cummins, Ignatius R Pranantyo, Jonathan M Pownall, Jonathan D Griffin, IrwanMeilano, and Siyuan Zhao,
Earthquakes and tsunamis caused by low-angle normal faulting inthe banda sea, indonesia , Nature Geoscience (2020), no. 4, 312–318.13. Masoumeh Dashti and Andrew M Stuart, The Bayesian approach to inverse problems , Hand-book of Uncertainty Quantification (2017), 311–428.14. B Davis and R LeVeque,
Adjoint methods for guiding adaptive mesh refinement in tsunamimodeling. , Pure & Applied Geophysics (2016), no. 12.15. TszMan L Fisher and Ron A Harris,
Reconstruction of 1852 Banda Arc megathrust earthquakeand tsunami , Natural Hazards (2016), no. 1, 667–689.16. Junichi Fukuda and Kaj M Johnson, A fully bayesian inversion for spatial distribution of faultslip with objective smoothing , Bulletin of the Seismological Society of America (2008), no. 3,1128–1146.17. Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin, Bayesian data analysis ,vol. 2, Taylor & Francis, 2014.18. Andrew Gelman, Donald B Rubin, et al.,
Inference from iterative simulation using multiplesequences , Statistical Science (1992), no. 4, 457–472.19. Frank I Gonz´alez, Randall J LeVeque, Paul Chamberlain, Bryant Hirai, Jonathan Varkovitzky,and David L George, Validation of the geoclaw model , NTHMP MMS Tsunami InundationModel Validation Workshop. GeoClaw Tsunami Modeling Group, 2011.20. Jonathan Griffin, Ngoc Nguyen, Phil Cummins, and Athanasius Cipta,
Historical Earthquakesof the Eastern Sunda Arc: Source Mechanisms and Intensity-Based Testing of Indonesia’s
ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 29
National Seismic Hazard Assessment , Bulletin of the Seismological Society of America (2018).21. Barbara Dix Grimes, . Mapping Buru: The Politics of Territory and Settlement on an EasternIndonesian Island , Sharing the earth, dividing the land: land and territory in the Austronesianworld (2006), 135–155.22. Warren Bell Hamilton,
Tectonics of the indonesian region , vol. 1078, US Government PrintingOffice, 1979.23. Latief Hamzah, Nanang T Puspito, and Fumihiko Imamura,
Tsunami catalog and zones inindonesia , Journal of Natural Disaster Science (2000), no. 1, 25–43.24. Thomas C. Hanks and Hiroo Kanamori, A moment magnitude scale , Journal of GeophysicalResearch: Solid Earth (1979), no. B5, 2348–2350.25. RA Harris, Temporal distribution of strain in the active Banda orogen: a reconciliation ofrival hypotheses , Journal of Southeast Asian Earth Sciences (1991), no. 3-4, 373–386.26. Ron Harris, The nature of the Banda Arc–continent collision in the Timor region , Arc-Continent Collision, Springer, 2011, pp. 163–211.27. Ron Harris and Jonathan Major,
Waves of destruction in the East Indies: the Wichmanncatalogue of earthquakes and tsunami in the Indonesian region from 1538 to 1877 , GeologicalSociety, London, Special Publications (2016), SP441–2.28. Gavin P. Hayes, Ginevra L. Moore, Daniel E. Portner, Mike Hearne, Hanna Flamme, MariaFurtney, and Gregory M. Smoczyk,
Slab2, a comprehensive subduction zone geometry model ,Science (2018), no. 6410, 58–61, Slab2 GitHub: https://github.com/usgs/slab2 , datadownload: .29. Arnauld Heuret, Serge Lallemand, Francesca Funiciello, Claudia Piromallo, and Claudio Fac-cenna,
Physical characteristics of subduction interface type seismogenic zones revisited , Geo-chemistry, Geophysics, Geosystems (2011), no. 1.30. DR Hilton and H Craig, A helium isotope transect along the Indonesian archipelago , Nature (1989), no. 6252, 906–908.31. Kruawun Jankaew, Brian F Atwater, Yuki Sawai, Montri Choowong, Thasinee Charoentitirat,Maria E Martin, and Amy Prendergast,
Medieval forewarning of the 2004 Indian Oceantsunami in Thailand , Nature (2008), no. 7217, 1228.32. Yan Y. Kagan,
Seismic moment distribution revisited: I. Statistical results , Geophysical Jour-nal International (2002), no. 3, 520–541.33. Yan Y. Kagan and David D. Jackson,
Seismic gap hypothesis: Ten years after , Journal ofGeophysical Research: Solid Earth (1991), no. B13, 21419–21431.34. Jari Kaipio and Erkki Somersalo, Statistical and computational inverse problems , AppliedMathematical Sciences, vol. 160, Springer Science & Business Media, 2005.35. Randall J LeVeque,
Finite volume methods for hyperbolic problems , vol. 31, Cambridge uni-versity press, 2002.36. Randall J LeVeque and David L George,
High-resolution finite volume methods for the shallowwater equations with bathymetry and dry states , Advanced numerical models for simulatingtsunami waves and runup, World Scientific, 2008, pp. 43–73.37. Randall J. LeVeque, David L. George, and Marsha J. Berger,
Tsunami modelling with adap-tively refined finite volume methods , Acta Numerica (2011), 211–289.38. Jun S Liu, Monte Carlo strategies in scientific computing , Springer Science & Business Media,2008.39. Z. Y. C. Liu and R. A. Harris,
Discovery of possible mega-thrust earthquake along the SeramTrough from records of 1629 tsunami in eastern Indonesian region , Natural Hazards (2014), 1311–1328.40. Alberto Malinverno, Parsimonious bayesian markov chain monte carlo inversion in a non-linear geophysical problem , Geophysical Journal International (2002), no. 3, 675–688.41. Stacey Servito Martin, Linlin Li, Emile A. Okal, Julie Morin, Alexander E. G. Tetteroo,Adam D. Switzer, and Kerry E. Sieh,
Reassessment of the 1907 Sumatra ”Tsunami Earth-quake” Based on Macroseismic, Seismological, and Tsunami Observations, and Modeling ,Pure and Applied Geophysics (2019).42. Robert McCaffrey,
The next great earthquake , SCIENCE-NEW YORK THENWASHINGTON- (2007), no. 5819, 1675.
43. W. R. McCann, S. P. Nishenko, L. R. Sykes, and J. Krause,
Seismic gaps and plate tectonics:Seismic potential for major boundaries , Pure and Applied Geophysics (1979), no. 6,1082–1147.44. K. McCue,
Seismic hazard mapping in Australia, the Southwest Pacific and Southeast Asia ,Annals of Geophysics (1999), 1191–1198.45. Aron J Meltzner, Kerry Sieh, Hong-Wei Chiang, Chuan-Chou Shen, Bambang W Suwargadi,Danny H Natawidjaja, Belle Philibosian, and Richard W Briggs, Persistent termini of 2004-and 2005-like ruptures of the Sunda megathrust , Journal of Geophysical Research: Solid Earth (2012), no. B4.46. Aron J Meltzner, Kerry Sieh, Hong-Wei Chiang, Chuan-Chou Shen, Bambang W Suwargadi,Danny H Natawidjaja, Belle E Philibosian, Richard W Briggs, and John Galetzka,
Coralevidence for earthquake recurrence and an AD 1390–1455 cluster at the south end of the2004 Aceh–Andaman rupture , Journal of Geophysical Research: Solid Earth (2010),no. B10.47. Aron J Meltzner, Kerry Sieh, Hong-Wei Chiang, Chung-Che Wu, Louisa LH Tsang, Chuan-Chou Shen, Emma M Hill, Bambang W Suwargadi, Danny H Natawidjaja, Belle Philibosian,et al.,
Time-varying interseismic strain rates and similar seismic ruptures on the Nias–Simeulue patch of the Sunda megathrust , Quaternary Science Reviews (2015), 258–281.48. Katrin Monecke, Willi Finger, David Klarer, Widjo Kongko, Brian G McAdoo, Andrew LMoore, and Sam U Sudrajat,
A 1,000-year sediment record of tsunami recurrence in northernSumatra , Nature (2008), no. 7217, 1232.49. KR Newcomb and WR McCann,
Seismic history and seismotectonics of the Sunda Arc ,Journal of Geophysical Research: Solid Earth (1987), no. B1, 421–439.50. Hendro Nugroho, Ron Harris, Amin W Lestariya, and Bilal Maruf, Plate boundary reorga-nization in the active banda arc–continent collision: Insights from new gps measurements ,Tectonophysics (2009), no. 1-2, 52–65.51. Yoshimitsu Okada,
Surface deformation due to shear and tensile faults in a half-space , Bul-letin of the seismological society of America (1985), no. 4, 1135–1154.52. , Internal deformation due to shear and tensile faults in a half-space , Bulletin of theSeismological Society of America (1992), no. 2, 1018–1040.53. Emile A Okal and Dominique Reymond, The mechanism of great banda sea earthquake of 1february 1938: applying the method of preliminary determination of focal mechanism to ahistorical event , earth and planetary science letters (2003), no. 1-2, 1–15.54. Emile A Okal and Costas E Synolakis,
Source discriminants for near-field tsunamis , Geo-physical Journal International (2004), no. 3, 899–912.55. Ignatius Ryan Pranantyo and Phil R Cummins,
The 1674 ambon tsunami: Extreme run-upcaused by an earthquake-triggered landslide , Pure and Applied Geophysics (2020), no. 3,1639–1657.56. Anthony Reid,
Two hitherto unknown Indonesian tsunamis of the seventeenth century: Prob-abilities and context , Journal of Southeast Asian Studies (2016), no. 1, 88–108.57. Gareth O Roberts, Andrew Gelman, Walter R Gilks, et al., Weak convergence and optimalscaling of random walk Metropolis algorithms , The annals of applied probability (1997),no. 1, 110–120.58. Gareth O Roberts and Jeffrey S Rosenthal, Optimal scaling of discrete approximations toLangevin diffusions , Journal of the Royal Statistical Society: Series B (Statistical Methodol-ogy) (1998), no. 1, 255–268.59. Charles M Rubin, Benjamin P Horton, Kerry Sieh, Jessica E Pilarczyk, Patrick Daly, NazliIsmail, and Andrew C Parnell, Highly variable recurrence of tsunamis in the 7,400 yearsbefore the 2004 indian ocean tsunami , Nature communications (2017), 16019.60. Valent´ı Sallar`es and C´esar R. Ranero, Upper-plate rigidity determines depth-varying rupturebehaviour of megathrust earthquakes , Nature (2019), no. 7785, 96–101.61. Kerry Sieh, Danny H Natawidjaja, Aron J Meltzner, Chuan-Chou Shen, Hai Cheng, Kuei-Shu Li, Bambang W Suwargadi, John Galetzka, Belle Philibosian, and R Lawrence Edwards,
Earthquake supercycles inferred from sea-level changes recorded in the corals of west Sumatra ,Science (2008), no. 5908, 1674–1678.62. Ihab Sraj, Kyle T Mandli, Omar M Knio, Clint N Dawson, and Ibrahim Hoteit,
Uncertaintyquantification and inference of mannings friction coefficients using dart buoy data during thet¯ohoku tsunami , Ocean Modelling (2014), 82–97. ECONSTRUCTION OF 1852 BANDA ARC EARTHQUAKE AND TSUNAMI 31
63. ,
Quantifying uncertainties in fault slip distribution during the t¯ohoku tsunami usingpolynomial chaos , Ocean Dynamics (2017), no. 12, 1535–1551.64. Jacob Swart, Verhandelingen en Berigten Betrekkelijk het Zeewegen en de Zeevaartkunde(English: Treatises and Reports Related to the Seaways and Nautical Sciences) , (1853),257–274.65. Albert Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation ,siam, 2005.66. Garrett W Tate, Nadine McQuarrie, Douwe JJ van Hinsbergen, Richard R Bakker, RonHarris, and Haishui Jiang,
Australia going down under: Quantifying continental subductionduring arc-continent accretion in Timor-Leste , Geosphere (2015), no. 6, 1860–1883.67. Donald L Turcotte and Gerald Schubert, Geodynamics , Cambridge university press, 2002.68. Donald L Wells and Kevin J Coppersmith,
New empirical relationships among magnitude,rupture length, rupture width, rupture area, and surface displacement , Bulletin of the seismo-logical Society of America (1994), no. 4, 974–1002.69. DJ Whitford, W Compston, IA Nicholls, and MJ Abbott, Geochemistry of late Cenozoic lavasfrom eastern Indonesia: Role of subducted sediments in petrogenesis , Geology (1977), no. 9,571–575.70. A Wichmann, The earthquakes of the Indian archipelago to 1857 , Verhandl. Koninkl. Akad.van Wetenschappen, 2nd sec (1918), no. 4, 1–193.71. , The earthquakes of the Indian archipelago from 1858 to 1877 , Verhandl. Koninkl.Akad. van Wetenschappen, 2nd sec (1922), no. 5, 1–209.(1922), no. 5, 1–209.