Reconstruction of the slip distributions in historical earthquakes on the Sunda megathrust, W. Sumatra

Abstract

Spatially heterogeneous inversions for the slip in major earthquakes are typically only available for modern, instrumentally recorded events. Stress reconstructions on active faults, which are potentially vital elements of earthquake hazard assessment, are usually based on one or two instrumental inversions and some semi-quantitative information concerning historical seismic activity. Here we develop a new Bayesian Monte Carlo inversion method for sparse, large uncertainty data, such as is available for many historical earthquakes. We use it to reconstruct the slip in two great historical earthquakes on the Sunda megathrust, from 200 yr of paleogeodetic data recorded in the stratigraphy of coral microatolls on the forearc islands. The technique is based on the stochastic forward modelling of many slip distributions, which are constrained by the observed fractal scaling in the slip field. A set of static elastic Green's functions is used to estimate the vertical displacement at the coral locations resulting from each trial slip distribution. The posterior expected value of the slip and its variance are obtained from the average of the trial slips, weighted by their ability to reproduce the coral displacements. We successfully test the method on a synthetic slip distribution, before applying it to reconstruct the slip in a recent, instrumentally recorded event. We then invert for the great 1797 and 1833 megathrust events under the Mentawai Islands west of Sumatra. When compared with the slip distributions in recent earthquakes, our results unambiguously indicate that the sequence is not consistent with the classic characteristic earthquake model; earthquakes tend to incrementally tessellate the active fault plane rather than repeatedly breaking a segment of it. The results indicate that homogeneous loading of a fault with heterogeneous initial stress is enough to explain the observations on the Mentawai segment of the Sunda megathrust; no time dependence of material properties, for example, is necessary. Furthermore, we speculate that even small amounts of nonlinearity in the rupture process would ensure that any tessellating sequence will not be repeated, even over very long times. This method could be adapted to incorporate other historical information, such as records of tsunami run up and Mercalli Intensity estimates. Finally, we suggest that reconstruction of stress based on long-term slip histories may prove useful in identifying the possible locations of future earthquakes.

1 INTRODUCTION

The sequence of great earthquakes on the Sunda megathrust since 2004, in particular the 2004 M9.2 Sumatra–Andaman and 2005 M8.7 Simeulue–Nias events on the northern part of the fault zone, have been responsible for massive loss of life. While the large moment release in 2004 and 2005 might suggest that these sections of the fault zone are not likely to produce another great earthquake in the near future, the most recent events under the southern Mentawai islands in September 2007 (two events of M8.4 and M7.9, referred to here as 2007a and 2007b respectively), September 2009 (M7.9) and October 2010 (M7.7) have not significantly reduced the moment deficit there (Konca et al.2008; McCloskey et al.2010). Additionally, a 300–400 km long section under Siberut island in the northern Mentawai islands (just off-shore from the city of Padang, population 1 million) has not produced a major event in this sequence, the last being more than 200 yr ago (Newcomb & McCann 1987). The risk of another great, tsunamigenic earthquake in Sumatra has not diminished and, by transferring stress onto the Mentawai section of the fault, previous events in the sequence may have increased this risk on the short to medium-term (Nalbant et al.2005; McCloskey et al.2010).

It has been shown that predicted tsunami wave heights along the west coast of Sumatra are sensitive to the distribution of moment release in the causative earthquake (McCloskey et al.2008). In assessing the risk to populations in Sumatra, we would therefore ideally like to gain information not only on the possible magnitude and rupture extent of an earthquake but also on the distribution of slip.

The amount of seismic slip in any event is constrained to a large extent by the accumulated strain on the plate interface, which is described to first order by the interseismic loading rate and the distribution of interplate coupling (Chlieh et al.2008; see Fig. 1). The maximum strain is accumulated where there is near complete coupling at the plate interface.

In the classic ‘characteristic’ model of earthquake rupture, where interseismic strain accumulates until it exceeds a limiting friction and the segment ruptures in a quasi-periodic sequence of similar events, we would expect great earthquakes to correlate both with the coupling and with each other. However, none of the slip distributions in 2005, 2007a or 2007b seem to correlate well with the coupling (Fig. 2; Briggs et al.2006; Konca et al.2008, see also Nalbant et al.2013). While it is true that the maximum slip in all cases occurred on strongly coupled areas of the megathrust, secondary slip-asperities seem to have occurred on low to intermediate coupling (under southern Simeulue in 2005 and south of Pagai in 2007a). Note that here we use ‘slip-asperity’ to refer to a patch of high slip in an earthquake and do not imply any causation due to a property of the fault surface.

This simplified, characteristic model relies on the assumption that all of the interseismically accumulated moment is released in each great earthquake. In reality, earthquake moment release is spatially heterogeneous and, as shown by recent observations (e.g. Lorito et al.2011 and Simons et al.2011), approximately co-located great earthquakes can produce very different slip distributions. The strain therefore probably retains a memory of an event for longer than what would be conventionally considered the length of the earthquake cycle. This heterogeneous pre-strain is therefore the result of both interseismic strain accumulation and ‘imperfect’ co- and post-seismic strain release in previous earthquakes, which have only partially reduced the moment deficit. This is a model where the rupture physics need not be complex, but the initial conditions of accumulated strain necessarily are. In this case, reconstructing the slip history on the megathrust may improve our picture of the total accumulated strain, thereby improving forecasts for the spatial distribution of moment release in future events. Another implication of this is that we should not expect one earthquake slip distribution, but the total cumulated history of co- and post-seismic slip to correlate with the coupling distribution, and that this correlation should improve as our knowledge of the slip history is extended to earlier great earthquakes on the interface.

Newcomb & McCann (1987) have compiled and summarized approximate rupture extents for a number of historical earthquakes under the Mentawai islands, based on historical intensity reports and tsunami impact on the Sumatran coast. The two largest of these, in 1797 and 1833, resulted in significant sea-floor deformation (generating tsunamis in both cases), which resulted in a drop in relative sea-level at the Mentawai islands. A record of the highest level of survival of corals microatolls (Taylor et al.1987), which can only grow up to the level of the lowest tide, is preserved in their stratigraphy. The corals may therefore be used to reconstruct deformation associated with historical earthquakes.

In combination with historical intensity records, damage and tsunami reports, which give information on the rupture extent of the events, this has allowed the broad pattern of slip in 1797 and 1833 to be estimated. Natawidjaja et al. (2006) have estimated slip distributions composed of several rectangular, uniform-slip rupture segments, which were validated by forward modelling coral displacements. However, the slip distributions are clearly simplified compared to the complexity of slip observed in modern inversions.

It has been shown (Mai & Beroza 2002) that the power spectrum of slip in earthquakes, as determined from seismic finite source inversions, is well described by |$\omega _k = {\text{c} /\text{{k^D }}}$|⁠, where ω_k is the power spectral density for wave number k and D is a fractal dimension (Turcotte 1997), which remains invariant at D ≈ 2.2 over all scales. This suggests that there is significant high-frequency variability in the magnitude of the slip, at least up to the spatial resolution of the inversion method. Even geodetic inversions observe slip heterogeneity (see, e.g. Johanson et al.2006; Subarya et al.2006), with a higher low spatial frequency limit.

Although the coral displacement data are sparser than modern GPS, InSAR and seismic data and contain only the vertical component of deformation, we suggest that the corals may contain more information about the distribution of slip than the composite box models suggest. Reproducing more of the slip heterogeneity will allow us to make a more detailed comparison between historical and modern slip distributions.

As an extension of the forward model validation approach of Natawidjaja et al. (2006), we propose that trial slip models which contain more realistic slip variability may produce historical slip distributions that are consistent both with the coral data and with the degree of slip heterogeneity expected from modern earthquake inversions.

In this paper, we describe a new technique: trial slip models are generated from fractal distributions; vertical displacements from the trial slip distribution are modelled at the location of published coral measurements for the historical events; and a solution is produced from an average of trial slip distributions, weighted by their ability to reproduce the coral uplift measurements. We investigate the quality of the solutions for 1797 and 1833 by applying the technique to recover the slip in a synthetically generated earthquake and the 2007 events, for which we have modern instrumental inversions (Konca et al.2008; Lorito et al.2008) for comparison. Finally, we calculate the total accumulated slip since 1797 in great earthquakes under the Mentawai islands and make recommendations for further work.

Displacements estimated from corals for the 1797 and 1833 earthquakes are taken from Natawidjaja et al. (2006), and we have also used the detailed summary of the historical information provided in that paper to constrain the trial slip models. The coral displacement measurements for 2007 were taken from Konca et al. (2008).

2 METHODOLOGY

2.1 Generating trial coral displacements

Chlieh et al. (2008) described a curved megathrust geometry at the Mentawai islands (based on locations and focal mechanisms of plate interface seismicity): δ = 5–7° at the trench, 15–20° under the islands and 30° below the coast of Sumatra. This was simplified to a planar geometry, with a best fit δ = 13°, for the formal inversion for the coupling, which is consistent with dip inferred from a recent tomographic study (Hayes et al.2012, SLAB1.0) under the Mentawai islands. Closer examination of the SLAB1.0 model shows that, while this dip might be appropriate for the source region of the 1797 earthquake, the interface is steeper (about 16°) to the south in the source region of the 1833 event. We will therefore leave the dip as a free parameter in our inversions. While depth-dependent dip would represent the geometry of the plate interface more accurately, displacements are least sensitive to slip where the curvature is largest at the trench and coast. Results are therefore unlikely to be sensitive to the curvature of the interface under the islands where the slip is well resolved.

In summary, we approximate the megathrust geometry using a planar fault surface (Fig. 3), defined by the latitude and longitude of its origin and the depth of the trench [ϕ₀, λ₀, d₀] (d₀ is allowed to take values between 5 and 7 km), its strike φ = 325°, maximum rupture dimensions along strike L and down dip W (see Table 1 for the range of allowed values for each earthquake) and its dip, which is allowed to take values between 5° and 20°.

Geodetic modelling has shown that the locked part of the megathrust, which can support large seismic slip, extends to 35–57 km depth (Simoes et al.2004) and seismicity hypocentres extend to at least 50 km at the Mentawai islands (Collings et al.2012). Taking the seismogenic thickness to be 50 km, we obtain a maximum down-dip rupture extent (or seismogenic width) of 220 km for δ = 13°. In order to remain consistent with the coupling grid produced by Chlieh et al. (2008), we therefore fix the maximum seismogenic width to 240 km.

We constrain the range of allowed magnitudes and rupture lengths for trial slip models according to available historical records and the distribution of coral data, based broadly on empirical relationships (see, e.g. Henry & Das 2001) for the scaling between L and M_w, when L is greater than the seismogenic width. Constraints are described in detail in Section 2.6 and Table 1.

1. Sample the origin, length and width of |${\boldsymbol \xi }_n $| uniformly from some range, so that it fits within |${\boldsymbol \xi }$|⁠.

2. Sample M uniformly from its allowed range and calculate the seismic moment M₀.

3. Generate a 2-D fractal distribution |${\boldsymbol \Psi } = \left[ { \cdots ,{\rm \Psi }_{kl} \in \left[ {0,1} \right], \cdots } \right]$| over the sub-area |${\boldsymbol \xi }_n $| (using the procedure in Turcotte 1997). The fractal dimension D is uniformly sampled from [1.5,2.5].

4. Set |${\bf u}_n = u^* {\boldsymbol \Psi }$|⁠, where

The set of trial coral displacements for the nth slip distribution, z_n = [⋅⋅⋅z_{k, n}⋅⋅⋅], is computed as follows: we define a set of Green's functions G_{k, ij}(δ, d₀, μ, ν) (where ν is Poisson's ratio) for the vertical displacement at coral k as a result of unit reverse slip on cell ij, with area δξ².

Green's functions are calculated using a homogeneous half space elastic model (Okada 1992). Elastic parameters are typically assigned values of μ ∼ 3 × 10⁵ bars and ν = 0.25 for homogeneous media. Previous work has shown that the correction to surface displacements required to account for the full heterogeneity of elastic behaviour would be on the order of 5 per cent (McCloskey et al.2008), which, particularly given the other sources of uncertainty in the solution slip distribution (chiefly poor data coverage), is unlikely to make a substantial difference to the results. However, in order to investigate the nature of the systematic errors that could be introduced into the solution slip distributions from the assumption of fixed, homogeneous elastic parameters, we allow both μ and ν to vary.

2.2 Weighting solutions

2.3 The Bayesian Monte Carlo method

The procedure outlined above has an interpretation in terms of the Bayesian Monte Carlo (BMC) method (see, e.g. Dilks et al.1992; Qian et al.2003). Rather than applying a prior probability distribution for the slip at each cell, where Monte Carlo simulations are weighted by both the prior probability and the likelihood, we apply the prior distribution of slip at a single location via the autocorrelation properties of the whole slip distribution. A prior distribution is therefore implicit in the Monte Carlo sampling of trial slip models. We can reconstruct the prior and posterior distributions during the inversion process, using the following procedure.

2.4 Notation

The posterior expected slip distribution is denoted U_Σ and has χ = χ_Σ. We distinguish between the magnitude and maximum slip of this solution, M_Σ and U_{max, Σ}, and the posterior expected values of the magnitude and maximum slip, E(M|D) ≡ 〈M〉 and E (U_max|D) ≡ 〈U_max〉. We refer also to the posterior expected dip and fractal dimension E(δ|D) ≡ 〈δ〉 and E(D|D) ≡ 〈D〉. We discuss the best trial slip distribution U₁ and its parameters, for example, its magnitude M₁, maximum slip U_{max, 1}, dip δ₁ and fractal dimension D₁.

2.5 Convergence

As observed by Qian et al. (2003), the Monte Carlo method in general requires very large numbers of simulations to avoid the risk of undersampling the correct solution space. This is particularly the case for large data sets, which is what makes the BMC approach possible in this application: where the data are sparse and the uncertainties are large, solutions converge relatively quickly.

However, in this context, problems of undersampling can also be offset to some extent by restricting the input ranges for parameters, such as the magnitude and the rupture dimensions (i.e. the BMC prior is again adjusted based on some semi-qualitative assumptions). If these constraints are poorly chosen, there is the risk of not sampling the correct space at all, so appropriate choices for the range of magnitude and rupture length are crucial to obtaining a solution that converges. The parameter input ranges we have chosen are discussed in Sections 2.1 and 2.6. A further constraint is provided by the choice of a fractal model for the slip distribution, which also greatly reduces the phase space to be searched.

Given the limitations of the BMC method, it is necessary to establish that the solution space has been adequately sampled and the slip solution has converged in each inversion. We discuss these issues with reference to the following:

1. Separate realizations of solutions for the same event should result in essentially the same posterior slip distribution, so for each earthquake we present results from two independent inversions.

2. We show the convergence of the values of χ, 〈M〉, 〈U_max〉, 〈D〉 and 〈δ〉.

3. We show posterior density functions for M, U_max, D, δ and also for the slip at four test locations A–D on the fault surface.

2.6 Coral data and historical constraints on the 1797 and 1833 Mentawai earthquakes

We have constrained magnitudes and rupture areas, summarized in Table 1, based on both historical records, as detailed in Natawidjaja et al. (2006), and the distribution of coral data.

The 1797 coral displacement set consists of 12 displacement measurements, extending from the south end of Pagai island to the Batu islands in the north of the study area. (There are 14 published displacements, but we find, as did Natawidjaja et al. (2006), that two of these cannot be reconciled with the other data. See the discussion in Section 6.) Damage reports from shaking in Padang, along with reports of a tsunami that had 5–10 m run-up at Padang and which also reached the Batu islands, suggest that there was significant slip as far north as Siberut island (Newcomb & McCann 1987; Natawidjaja et al.2006). We have allowed the rupture to extend northwards to the Batu islands, under the northernmost coral displacement measurement. However, unfortunately, the data in the northern Mentawais are limited to one displacement measurement on Siberut (with a large error) and one on the Batu islands; the slip is likely to be poorly resolved in this area, which will have implications for our ability to define the earthquake and tsunami hazard in Padang. Towards the south, the data density is much better: there are three measurements on Sipora and seven on Pagai. Since there are no records of damage as far south as Bengkulu, which was a major settlement at the time, we have assumed the rupture did not extend much further south than Pagai island. The fault origin is therefore at the trench, just south of Pagai, and we allow values for the rupture length of between 300 and 480 km, so that the slip extends to just south of the Batu islands. For consistency with empirical relations between the magnitude and rupture length, we constrain the magnitude to lie in the range M = 8.3–9.0, which includes the suggested M = 8.5–8.7 of Natawidjaja et al. (2006).

There are 15 coral displacement measurements for the 1833 event, including the corals used for 1797 and several more on southern Pagai. Records of damage, shaking intensity and tsunami reports suggest that the earthquake had significant slip as far south as Bengkulu, while the duration of the shaking was shorter and the damage less severe in Padang, suggesting there was not significant slip under Siberut Island. Without data south of Pagai, it is not possible to resolve slip as far south as Bengkulu, so we have chosen to restrict the southern extent of the rupture to approximately 150 km southeast of Pagai island; allowing slip models to extend further increases the computational time before the solution converges, without gaining any additional information about the slip beyond the prior input information. The origin of the fault surface is at the trench, about 150 km southeastwards along strike from south Pagai. Slip in the north under Siberut will, again, be difficult to resolve. However, as before, we have allowed the models to extend as far north as the Batu islands. As we will see, the 1833 displacements taper to near zero on Siberut and the Batu islands and their uncertainties are smaller than for 1797, which will rule out most slip distributions that place large slip here. Trial rupture lengths are constrained to lie in the range 400–600 km, so that the slip extends as far north as the Batu islands. We constrain the magnitude to lie in the range M = 8.3–9.2, which includes the suggested M = 8.9 of Natawidjaja et al. (2006).

3 RECOVERY OF A SYNTHETIC SLIP DISTRIBUTION

The ability of this method to resolve slip on the fault surface will be dependent on the sensitivity of the coral network to slip, as a function of its location. The along-strike distribution of data is good for all study earthquakes, particularly on Pagai island, but the fact that they are confined to the islands in a band around 20–50 km wide in the trench-coast direction, which thins to only one coral on Siberut, suggests that while it may be possible to resolve relatively strong gradients in slip in the strike direction, down-dip resolution is likely to be poor.

In order to investigate the implications of this for the quality of inversions for the historical events, we have generated a synthetic test fractal slip distribution (Fig. 4a), with M = 8.6 and D = 2.5, broadly based on the parameters of the 1797 earthquake (see Table 2).

We forward model synthetic vertical displacements from this earthquake at a test coral network based on the 1797 coral locations (see Fig. 6 for the displacements). We assign a random error to each displacement, by sampling from a Gaussian distribution with a standard deviation equal to the uncertainty in the corresponding 1797 displacement measurement. Using the algorithm in Section 2 and the prior inputs associated with the 1797 event (see Table 2 for the allowed ranges for the input parameters and Fig. 4b for U_Prior), we invert to recover the test synthetic slip distribution.

Two realizations, containing 1 × 10⁸ trial slip models, for the posterior expected values for the slip in the synthetic earthquake are shown in Figs 5(a) and (b), along with the posterior estimate of the standard deviation (Fig. 5c).

The coral displacements (Fig. 6) are reproduced well by the solution slip distribution and the χ_Σ-values (Fig. 7a) converge rapidly.

The inversion reproduces the large slip-asperity under S. Pagai island and is able to resolve the tapering slip under N. Pagai, towards Sipora. The slip-asperity 60–80 km down dip of N. Pagai and Sipora seems to be moved up-dip, to less than 40 km to the NE of the islands, suggesting that although there is some resolution down-dip, the solutions are not reliable far from the islands. This is particularly the case towards the trench: the standard deviation of the slip is very large here, suggesting that individual trial slip models can have arbitrarily large slip, with no penalty in terms of fit to the data. The best trial slip model (Fig. 5d) has patches of large slip both at the trench and under the Sumatran coast line.

Additionally, this set of coral locations does not have good enough coverage for the solution to be able to resolve the large slip under Siberut island, confirming that we will not be able to determine the slip under Siberut in either of the historical events (the data coverage is not much better for 1833) using the corals alone.

We find as a result that although the values of the magnitude and maximum slip (Figs 7b and d) have converged well and their posterior density functions are peaked (Figs 8a and c), their values are underestimated by the solution: we obtain M_Σ = 8.50 and U_{max, Σ} = 4.4 m. We suggest, however, that this could also be due to the smoothing of the slip distribution in the final solution: we find that the moment magnitude of the best trial slip model is closer to what we expect at M₁ = 8.63 (Fig. 7b) and the fractal dimension of the best model is large at D₁ = 2.4 (Fig. 7c).

For data sets of this quality, we do not expect roughness in the slip to be clearly resolvable (the roughness of the final solution will be determined by the data quality), but the posterior density function for the fractal dimension does show a slight tendency for the solution to favour larger fractal dimensions (see Fig. 8b). 〈D〉 = 2.1 (for both realizations of the solution), rather than 2.0, and D_KL(D) is very slightly positive at 0.04 and 0.05 for the two solutions.

We find that there does seem to be some information on the fault geometry contained in the coral displacement data. δ for the interface converges to a value of 11.6° (with a corresponding D_KL(δ) = 0.24), which is smaller than the prior expected value of 15° and closer to the input value of 10° for the synthetic. This result is more likely to be the result of the varying modelled depth of the interface under the islands (which has the strongest influence on the solution) than direct resolution of the dip itself.

However, in addition, we would expect the posterior density function for the dip to have a peak around the input value, which is not what we find (Fig. 8d). We find when the inversion for the synthetic slip is repeated with the observation errors set to zero, the posterior density for the dip has a better peak around the expected value for the dip 13.5°, but this value is also closer to the prior than the input value and has a smaller D_KL(δ) = 0.05. We therefore suggest that posterior values for the dip (which are only weakly resolved) may be dominated by errors in the displacements. We note that these results are reproduced exactly between the two realizations of the solution and the parameter values have converged to a very small tolerance (e.g. <0.01° for the dip). The apparent information gain is therefore not the result of inadequate convergence.

With respect to the material parameters and the depth to the trench, we find no evidence in the solutions that the coral displacements contain information on either μ or d₀. However, there is a weak preference for Poisson ratios of between 0.2 and 0.4. The slight trade-off between the shear modulus and the maximum slip suggests that the magnitude of the posterior slip distribution will depend in detail on the assumptions made regarding elastic properties. However, we find no significant difference in the geometries of the solution slip distributions with different elastic constants or with d₀; the conclusions of this paper are robust to uncertainties in these parameters.

Finally, in Fig. 9 we have plotted the prior and posterior density functions for the slip at four test locations on the fault plane (see Fig. 4a and Table 3 for their coordinates). The prior density functions have large densities at zero slip, since the trial slip distributions have, in general, a smaller area than the full model fault plane. All four posterior distributions are well converged and are reproduced exactly between the two realizations. The changes in the posterior densities with respect to the prior are reflected in the values of D_KL(u) of the solution at each location. We would expect D_KL(u) to be relatively small for the locations off the islands, which is what we find (see Table 3): the standard deviation of the slip towards the trench is large, which is consistent with the small D_KL value at location C. Although the standard deviation of the slip towards the Sumatran coastline is smaller than at the trench, location D is quite far from patches of large slip, which would also contribute to a small D_KL(u) value. By comparison, locations A and B have high D_KL(u) values (particularly location A which is beneath the best data on Pagai island).

It is clear from these results that χ_Σ is not sufficient to demonstrate convergence to the correct slip distribution: for sparse, poorly distributed data, χ may easily approach 1 if the displacements are not sensitive to the slip. Additionally, while χ₁ may be larger than χ_Σ, U_Σ is the better solution in general, since individual trial slip models may contain arbitrarily large slip in regions where the data are sparse, whereas the posterior expected slip will converge to the prior expected value in these locations: comparing the posterior expected slip both with the prior and with the standard deviation will allow the identification of areas which are poorly resolved by the data.

4 THE 2007 EARTHQUAKES

We apply this method to recover the published slip distributions from the modern instrumentally recorded earthquakes in 2007 (Konca et al.2008), shown in Fig. 2. The instrumental solution is the result of joint inversions of InSAR, GPS, seismic and coral-uplift data: vertical coral displacements account for a relatively small amount of the data used, so although the coral data will account for some part of the solution, it is likely to be dominated by the other data. We also have a slip distribution for the 2007a event from inversion of tsunami data (Lorito et al.2008), for comparison. As a result, this provides an opportunity not only to test the adequacy of our approach for this inversion problem, but also to test the general consistency of the coral data solutions (including our set of physical assumptions) with those based on other data sets.

4.1 Data and constraints

Common features of the instrumental inversions suggest that the first M8.4 earthquake in 2007 had an epicentre about 70 km south of Mega island (the southernmost coral measurement), and ruptured unilaterally northwest. According to Konca et al. (2008), there were two major slip-asperities, one just north of the epicentre, under Mega, and another under the southern end of Pagai. The M7.9 event nucleated down dip of this asperity, about 12 hr later and ruptured in two distinct patches; one down dip of Pagai island and one down-dip of Sipora. The total moment magnitude of the combined ruptures was approximately 8.5.

We have 18 coral data points for 2007, distributed, as before, between Sipora, Pagai and also on two small islands to the south. The southernmost asperity is in a region of sparse data coverage, with two measurements on Mega, which are however very close together. Slip under Pagai is well covered, as before. While there are displacement measurements on north Pagai and Sipora, the data seem to have a very narrow trench-coast distribution and are therefore unlikely to resolve the small asperity under Sipora. We have allowed trial slip models to extend south to the epicentre of the M8.4 event, and north as far as Siberut, with a range of magnitudes of 8.0–8.8.

4.2 Results and discussion

Slip in 2007 is broadly reproduced in solution U_Σ (Figs 10a and b) after 8 × 10⁹ simulations. The slip distribution under Mega and S. Pagai is reproduced particularly well and in both solutions the low slip between N. Pagai and Sipora is very clearly resolved. As expected, the low slip patch towards the Sumatran coast from Sipora does not appear and the down-dip extent of slip at Sipora is not reproduced. We note that although the slip solutions in Figs 10a and b have in common all the major features of the published slip distribution, the solution in Fig. 10b has some artefacts, most notably under Siberut island, suggesting that the slip solution has not fully converged.

The coral displacements (Fig. 11) are, in general, reproduced well and the χ_Σ values appear to be close to convergence (Fig. 12a), although the χ_Σ after 8 × 10⁹ is small. The posterior expected values of 〈M〉 = 8.70 and 〈U_max〉 = 14.2 m have converged well (Figs 12b and d). Both the magnitude and the maximum slip show large D_KL values of 2.11 and 4.47 bits, respectively. We find 〈δ〉 = 9.2° and 〈D〉 = 2.3, with D_KL(δ) = 2.19 bits and D_KL(D) = 2.07 bits. However, although these values appear to have converged and the density functions (Fig. 13) show peaks around the expectation values, the densities do not look well sampled, suggesting again that there is a problem with convergence. Similarly, although slip appears to be very well resolved at all but the last of the test locations, the posterior density functions are not properly sampled (see Fig. 14 for the posterior density functions and Table 3 for D_KL(u)). The magnitude is considerably larger than the M8.5 magnitude of Konca et al. (2008) and the posterior density function (Fig. 13a) for the magnitude, which is strongly peaked around the expectation value, suggests that an M8.5 event would not explain the coral displacements well. One explanation for this could be that we must assume the displacements were produced by one event rather than two subevents; it is possible that the magnitude–length scaling is not correct under these circumstances. Alternately, the magnitude of Konca et al. (2008) might be dominated by other data sources in the joint data set.

With respect to the convergence of the solutions, the 2007 data set contains more displacements with much smaller uncertainties than the historical data sets, which would lead us to expect large computation times to be required for convergence. However, we find that the convergence does not improve even after 30 × 10⁹ trial slip models. Other causes for the relatively poor convergence are suggested by the coral displacement data: there is a large gradient in displacement between the two corals at around −3.5°S latitude, which is not reproduced by any of the solutions. This suggests either that there is some error in the displacements that is not well reflected in the published uncertainties or that simple double-couple failure on the planar fault model we have assumed here is not adequate. Local anelastic deformation of the accretionary prism, activity on a splay fault or post-seismic movement of a coral head are just a few possible causes. Any of these could explain the failure of the solution to converge completely.

The relatively large value for 〈D〉 suggests that slowly varying slip distributions do not explain the 2007 data well. This might be expected, given that the data quality is very much better for this event than for 1797, so any high-frequency variability in the slip would be better resolved. However, a preference for slip models with high spatial frequencies could also be caused by data errors or some complexity in the fault geometry (e.g. activation of a splay fault) not accounted for here. Finally, the fact that for the second realization of the solution, we find 〈δ〉 = 9.7° and, again, D_KL(δ) is large, at 3.16 bits, might again point to data or model error, rather than noise resulting from too small a number of trial slip solutions.

5 SLIP SOLUTIONS FOR THE HISTORICAL EARTHQUAKES

5.1 1797

Solution slip distributions for 1797 are shown in Fig. 15. The two realizations of the solutions are indistinguishable (Figs 15a and b), coral displacements are well reproduced by the posterior model (Fig. 16) and χ_Σ converges rapidly (Fig. 17a), suggesting the solution is well converged after 1 × 10⁹ trial slip models. The standard deviation (Fig. 15c) of the slip is less than 3 m over much of the fault. At the trench and under Siberut, however, it can be as large as 5 m.

We find |$M_\Sigma = \langle M\rangle = 8.67$| and 〈M〉 converges very quickly (Fig. 17b). However, these values are considerably smaller than the magnitude of the best trial model M₁ = 8.82 and the magnitude is not very well resolved (see Fig. 18a for the posterior density function), suggesting that there may be large unconstrained slip in the trial slip solutions; under Siberut, in particular, the slip values are the same as the prior expected slip. The maximum slip is, similarly, poorly resolved, with a wide posterior density function (Fig. 18c). 〈U_max〉 is significantly larger at 12.8 m than U_{max, Σ} = 8.5 m, which suggests that there is considerable smoothing during averaging of the slip distributions. That there is only a weak preference (see the posterior density function in Fig. 18b) for large fractal dimension trial slip models is consistent with this result. 〈D〉 = 2.1, the best trial model has D₁ = 1.6 and D_KL(D) is small at 0.08 (see Fig. 18b for the posterior density function).

The results here are very much consistent with what was expected from the inversion of the synthetic earthquake: resolution is relatively good under Pagai and Sipora islands and less reliable elsewhere, particularly towards the trench.

One interesting feature occurs under south Pagai, where the slip decreases rapidly (over 20–40 km) from around 8 to 2–3 m. This feature of the slip is extremely robust and occurs in a region of quite small standard deviation. Additionally, test locations A and B, which lie under the islands in this region, have positive D_KL(u) values of 0.46 and 0.49 (see Table 3 and Fig. 19). These values are relatively small, but they are extremely robust (i.e. they do not vary at all between different realizations of the solution) so we suggest that the low slip here may be a real feature of the historical earthquake.

We note that 〈δ〉 = 10.7°, which is less than the prior expected value of 15°; also D_KL(δ) = 0.40, which is small but positive. We discuss this in Section 5.2, by comparison with the corresponding result for 1833.

Finally, the slip distribution is not materially affected by performing the inversion on a coarse 40 km grid (Fig. 15d). The location and extent of the slip are broadly the same as with the 20 km grid, although U_{max, Σ} increases to 10.2 m and M_Σ = 8.71.

5.2 1833

The solution slip distributions for 1833 are shown in Fig. 20. The solution converges well in 1 × 10⁹ trial slip models: the two realizations of the solution are identical, the coral displacements are very well reproduced (Fig. 21) and the parameters of the inversion converge within very small tolerances (Fig. 22). Again, the standard deviation of the slip is low over most of the fault surface, but increases to between 6 and 8 m at the trench in the South.

The moment magnitude is larger than for 1797, at M_Σ = 〈M〉 = 8.77, as expected, with U_{max, Σ} = 15.2 m. Although M₁ is not significantly larger than M_Σ, we find that 〈U_max〉 is significantly larger at 20.1 m than U_{max, Σ}, suggesting, as for 1797, that there is considerable unresolved roughness in the slip distribution. Again, the fact that there is only very weak preference for high fractal dimension slip models (〈D〉 = 2.1, with D_KL(D) = 0.05, see Fig. 23b for the posterior density function) is consistent with this result.

The posterior density functions at the test locations are well sampled (Fig. 24) and are the same for the two realizations. They show similar D_KL(u) values to the 1797 results (see Table 4 for the coordinates of the 1833 test locations and D_KL(u)). The slip under S. Pagai shows tapering towards the north, with a decrease from almost 15 to 8 m in 40–50 km. This is in the best resolved part of the fault so, again, we suggest that the detail of the slip is quite reliable under S. Pagai.

A final point to note is that the posterior expected value for the dip for the 1833 solution is significantly larger at 14.7° than for the 1797 solution. D_KL(δ) almost zero at 0.02, but there is a slight peak in the posterior density around the prior expected value for the dip (Fig. 23d), where the prior density is uniform, and this is reproduced in the second realization of the solution (〈δ〉 = 14.6° and D_KL(δ) = 0.03). Although we have shown that posterior values for the dip can be sensitive to random errors in the observations, this result is also consistent with seismological models for the geometry of the interface at this margin (SLAB1.0, Hayes et al.2012), which show an increase in dip of approximately 12–16° between Sipora island and S. Pagai. Since the location of the maximum slip in 1797 and 1833 coincide with regions of the interface with quite different dip, this could suggest the inversion is sensitive to the model dip, at least in as far as it determines the differing depth to the interface under Sipora and Pagai islands.

Resolving whether this result is a reflection of the data errors or a true reflection of the megathrust geometry is beyond the scope of this work, but we suggest that repeating the inversions

1. using a 3-D interface geometry and

2. applying random errors to the coral displacement observations (to estimate the possible effect on the solution for the dip)

could provide further insight into this issue.

5.3 Towards recreating the slip history under the Mentawai islands

See Fig. 25 for a summary of the history of earthquake slip under the Mentawai earthquakes since 1797. The maximum slip-asperity for 1797 occurs under north Pagai, where for 1833 it occurs under south Pagai. To a large extent these slip distributions are complimentary, in that they sequentially tile the regions of relatively high interface coupling (Fig. 1), particularly where data coverage is best at Pagai (notwithstanding the problems with resolution in the down-dip direction). They clearly do not correlate with each other, which would be expected if they were consistent with the characteristic earthquake hypothesis. We note that these findings are entirely consistent with the conclusions of Simons et al. (2011), who found that the Giant Tohoku-Oki earthquake completed the rupture of the Pacific-Japan megathrust in 2011.

The northern slip-asperity of the 2007 M8.4 event is concentrated under and slightly up-dip of south Pagai, overlapping the maximum slip in the 1833 earthquake. Additionally, the two slip-asperities in the M7.9 earthquake seem to be either side of the largest asperity in the 1797 earthquake, between Pagai and Sipora, which would to some extent explain the pattern of moment release there in 2007. These results would suggest that the complex pattern of slip in the 2007 earthquakes was at least in part constrained by the slip in 1797 and 1833.

These results support a model where complex conditions of initial stress dictate the geometry of slip in an earthquake, without any requirement for balance between loading and moment release on the timescale of one or even a number of great earthquakes. Although this loading-moment release balance must, of course, be satisfied overall, the characteristic timescale for this to occur is a matter for further research.

The greatest seismic and tsunami hazard facing the population of west Sumatra is from the megathrust under Siberut island. A tsunami from such an event would pose a serious threat to the city of Padang; the death toll in the event of a great earthquake there could be very large. Unfortunately, since the data coverage is so poor on Siberut (see Section 3), we are unable to accurately estimate the cumulated slip under the island. However, the data do indicate the (poorly constrained) possibility that slip under Siberut for the last 200 yr has been modest and that there is therefore significant accumulation of moment deficit there. This issue is obviously of great social as well as scientific consequence. As discussed in Section 5 below, incorporating historical data more explicitly into the modelling may help to improve resolution and thereby clarify the hazard for Padang. This study is not able to rule out a near-future great earthquake under Siberut island.

6 DISCUSSION

While our information on the slip in the 1797 and 1833 historical earthquakes on the Sunda megathrust is still far from complete, this method has contributed to resolving the pattern of moment release in both events.

Advantages of the technique are the ability to quantify the resolution of slip as a function of location on the fault surface (particularly important for sparse, large error data sets), and the ability to introduce prior assumptions about the slip heterogeneity numerically. However, it takes long computation times to get convergence for large, high-quality data sets. This technique is therefore not, at present, an alternative to standard inversion methods in most modern applications. Rather, it is well suited to historical inversion problems, where the information is limited by imprecise and incomplete records.

Based on GPS observations of time-dependent deformation in the aftermath of the 2005 and 2007 events (Hsu et al.2005; Konca et al.2008), post-seismic aseismic moment release can be as much as 30 per cent of the co-seismic moment. While modern coral measurements were made within several weeks of the earthquakes, historical displacements record a year of total deformation between growth rings in the coral colonies, which could contribute significant adjustment to the observed vertical displacements. For the 2005 and 2007 events, most of this displacement occurred up-dip of the islands, towards the trench, suggesting that the co-seismic displacements for 1797 and 1833 could be underestimated. Given the poor resolution on this part of the fault, it is not clear how far this would influence our results, though it is unlikely to make a significant difference to the main conclusions of this work. However, the post-seismic component of the slip on the seismogenic part of the interface also contributes to the reduction of moment deficit, so for these purposes, the cumulative moment release from both co- and post-seismic slip may be the important feature, rather than the co-seismic alone.

The physical model of the fault zone we have used includes two important simplifications: the assumption that deformation is the result of failure on a planar fault interface and the assumption of material elastic homogeneity. The stochastic approach we have used allows us to quantify errors introduced by these assumptions only partially and, while the uncertainties associated with poor data coverage may dominate these errors in general, more complex models for the medium and the fault geometry will undoubtedly improve solutions, if they can be well enough constrained. In particular, the fact that a number of coral measurements on the southern end of Pagai island in both the 1797 and 2007 earthquakes cannot be reconciled with the other displacement measurements may suggest either that there was some failure on a subsidiary structure such as a splay fault or that there is a problem with the observation (e.g. some movement of the coral head unrelated to co-seismic uplift). Although it is unlikely that resolving problems with a small number of corals will change the broad geometry of the slip distributions or their moment, these issues are important in the general context of inversion of great subduction earthquakes.

It may be difficult, for a number of reasons, to generalize this approach to other subduction zones. In Chile and Japan, for example, there are no islands along most of the forearc, so we are unlikely to be able to obtain deformation measurements that have reasonable coverage near the main asperities in historical earthquakes. The importance of good data coverage for resolving heterogeneity in slip has further support from modern inversions for the 2010 Tohoku-Oki earthquake in Japan: sea-floor geodetic instruments make a significant difference to the solution slip distribution, compared to inversions based only on deformation at the coast (Romano et al.2012). Similarly, joint inversions that include observations of tsunami wave heights, since they are sensitive to the sea-floor displacement, show significant increase in resolution towards the trench.

This may suggest a way of improving the data set for historical earthquakes. Historical records of tsunamis and damage from shaking, while they have large errors compared to modern measurements and are again very sparsely distributed, may provide an additional source of information, which can be added in via the prior, or, more quantitatively, used to adjust the weights of trial slip model. This would require a detailed review of the available historical data in order to obtain as much quantitative information from the records as possible, but also to estimate the likely errors associated with this data. In addition, the stratigraphy and geochronology of cores through tsunami deposits along the coast adjacent to the megathrust might provide further important data. At the very least this approach might be expected to provide better constraints on the rupture extent and magnitude of historical earthquakes, but it may also contribute to resolving the slip in areas where there is poor coral coverage.

The tsunami records at Padang may say something more about the slip under Siberut island in 1797, for example. Historical records seem to indicate that there was significant slip here in 1797, but not in 1833: tsunami modelling may determine whether large slip is required under Siberut to produce the recorded 5–10 m run-up at Padang, or whether large slip further south is sufficient. This question is particularly important since, undoubtedly, significant moment deficit has accumulated there since the historical events: as noted by Konca et al. (2008), while the 2007 earthquakes occurred on a region that had significant slip in 1833, they released only a fraction of the moment accumulated since then. These results are consistent with the view that there is potential for a great M > 8.5 earthquake under the Mentawai islands and that the possibility of significant slip under Siberut cannot be discounted.

Finally, the interpretation of the cumulated slip in terms of slip deficit on the megathrust, and therefore its application to forecasting, relies on the assumption that the coupling distribution is stable over time. There is some evidence that this is the case on the Mentawai section of the fault zone (Natawidjaja et al.2007), but it is less clear in northern Sumatra in the rupture area of the 2005 earthquake (Meltzner et al.2012). In other subduction zones, while we may be able to estimate the present coupling distributions using modern geodetic measurements, without information on historical deformation rates, it will be difficult to demonstrate that the distribution is a stable feature of the subduction zone.

7 CONCLUSIONS

These results indicate strongly that the total moment deficit on the Sunda megathrust is likely to be a strong control on the magnitude and slip distribution of the next great earthquakes there. While, to first order, the moment deficit is controlled by (1) the distribution of coupling and (2) the co- and post-seismic slip history at the plate interface, our ability to estimate it accurately is limited by the relatively short period for which there has been instrumentally recorded data for the great earthquakes. The quality and quantity of modern data allow inversions for slip that, at some level, contain realistic heterogeneity. However, understanding the accumulation and release of moment deficit on the subduction interface is likely to require longer term knowledge about the detail of slip history on the megathrust than modern records can provide.

We have developed a procedure, based on the BMC methodology, to invert for the slip in historical megathrust earthquakes from paleogeodetic data, where we assume slip heterogeneity a priori. The technique allows posterior and prior probability distributions for the slip to be estimated as a function of location on the fault surface: understanding the resolution is crucial for this application, given the sparse distribution and relatively large errors compared to modern data sets.

With the aim of reconstructing the slip history in great earthquakes on the Sunda megathrust, we have applied this technique to obtain solutions for the distribution of slip in the 1797 and 1833 earthquakes under the Mentawai islands. We obtain moment magnitudes of 8.8 and 9.0 for the 1797 and 1833 events. However, since the rupture lengths of the historical events are constrained only by historical records, and neither event has good data coverage around the estimated north and south termination of the ruptures, the magnitudes of the events are not well constrained beyond the prior.

That we are able broadly to recover the published slip distributions for the modern events in 2007, given data which are of comparable density and coverage to the historical measurements, confirms that the method can resolve slip heterogeneity, even with very sparse data. However, synthetic tests demonstrate that the resolution of the slip is strongly dependent on the data coverage. On one hand, the locations of measurements for the 1797 and 1833 deformation give very good resolution under Pagai island. On the other hand, very small values of KL divergence indicate poor resolution under Siberut to the north, and also at the coast and trench, where the data coverage is sparse: on locations of the megathrust where coral displacements are not sensitive to the slip, the posterior distributions for the slip at these locations are very similar to the prior input distributions.

The good constraints under Pagai island for the 1797 and 1833 inversions allow us to identify tessellation of the slip between these two earthquakes, suggesting that the distribution of moment release in 1833 was controlled by the regions of maximum slip in 1797. Similarly, the complex pattern of slip exhibited in the two 2007 ruptures is explained to a large extent by the pattern of cumulated slip in the 1797 and 1833 events. This suggests, first, that the characteristic earthquake model of a quasi-periodic sequence of similar earthquakes is inadequate to explain the heterogeneity visible in both historical and modern inversions for slip on the megathrust and, second, that a ‘slip-deficit’ model, where a heterogeneous distribution of pre-strain results from both heterogeneous interseismic strain accumulation (due to the interplate coupling) and heterogeneous slip in previous events on the segment, has real potential to give constraints on the distribution of slip in future earthquakes.

Since, with present data constraints, the slip history under the Mentawai islands is poorly resolved under Siberut island, we are unable to determine the cumulated slip deficit on this segment since 1797. We suggest that incorporating other historical information, such as damage and tsunami reports, more explicitly in the modelling may allow us to improve the technique further by adding constraints where coral measurements are sparse. However, even if, as historical information seems to suggest, this segment experienced significant slip in 1797, it is likely that significant slip deficit has accumulated there since: the hazard for Padang as the result of a great M > 8.5 event is still high.

This work was funded by NERC under grant number NE/H008519/1. MNB is currently funded under the Griffith Geoscience Research Award from the Department of Communications, Energy and Natural Resources under the National Geoscience Programme (Ireland). The authors wish to thank the editor for her careful reading and correction of the paper and for constructive criticism of our presentation. We also want to express our sincere thanks to two anonymous reviewers whose perceptive and knowledgeable critique improved this contribution significantly.

REFERENCES