Slip distribution of the 2005 Nias earthquake ( M w 8.6) inferred from geodetic and far-ﬁeld tsunami data

geodetic and far-ﬁeld tsunami data, which is also similar to slip distributions from previous studies based on local geodetic data. This demonstrates that far-ﬁeld tsunami waveforms, once corrected for propagation effects, can be used to estimate the slip distribution of large submarine earthquakes leading to results that are similar to those obtained using sparse local geodetic data.


I N T RO D U C T I O N
In Indonesia, large interplate earthquakes have repeatedly occurred along the Sunda and Java trenches (Fig. 1). The slip distributions on the fault of some significant earthquakes have been estimated by tsunami waveform inversions: the 2004 Sumatra-Andaman earthquake , the 2006 West Java tsunami earthquake (Fujii & Satake 2006), the 2007 Bengkulu earthquake (Fujii & Satake 2008) and the 2010 Mentawai tsunami earthquake (Satake et al. 2013). The source models of other historical earthquakes in 1797 and 1833 were estimated from the sea level changes reconstructed from coral microatolls (Natawidjaja et al. 2006). Plate coupling models have been constructed on the basis of the slip distribution of past earthquakes and current geodetic observations in the region west of Java (Hanifa et al. 2014) and off Sumatra (Chlieh et al. 2008;McCloskey et al. 2008). They are used to assess the future activities of the interplate earthquakes along the Sunda and Java trenches.
The Nias earthquake (2.085 • N, 97.108 • E, depth 30.0 km, M w 8.6 at 16:09:36 UTC according to USGS) occurred on 2005 March 28, three months after the 2004 Sumatra-Andaman earthquake (M w 9.1). It was one of the largest earthquakes in the 21st century. Tsunamis generated by the 2005 Nias earthquake caused less severe damage compared to the 2004 earthquake, with maximum heights of up to 4.2 m on Simeulue and Nias islands (Borrero et al. 2011). These heights are substantially smaller than those from the 2004 earthquake, and even smaller than expected from its large earthquake magnitude (Kerr 2005). The 2005 tsunamis were recorded not only at tide gauges (TGs) in and around the Indian Ocean, but also at far-field ocean bottom pressure gauges (OBPGs) near the Syowa Base in Antarctica. A tsunami source model, however, has not yet been constructed from the tsunami waveform data. 1162 C The Author 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.  , the 2006 West Java tsunami earthquake (Fujii & Satake 2006), the 2007 Bengkulu earthquake (Fujii & Satake 2008) and the 2010 Mentawai tsunami earthquake (Satake et al. 2013). The red and blue stars are the epicentres of the 2004 Sumatra earthquake and the 2005 Nias earthquake, respectively. The black rectangles show the source models of historical earthquakes in 1797 and 1833 by Natawidjaja et al. (2006). The region surrounded by a black line is the plate coupling model by Hanifa et al. (2014).
The USGS finite-fault model (https://earthquake.usgs.gov/earthq uakes/eventpage/of ficial20050328160936530 30/finite-fault, G38) inferred from teleseismic body waves shows deep slips of up to 9 m. Two large slip patches, with maximum slips of 15 and 9 m, were estimated beneath Simeulue and Nias islands from seismic and geodetic data (Konca et al. 2007). These maximum slip amounts are relatively small for an earthquake with M w 8.6 and consistent with relatively small tsunami runup heights. Slip distribution models were also derived from local geodetic data, for example, coral microatoll level change measurements (Briggs et al. 2006) and GPS measurements (Kreemer et al. 2006). In the source region of the 2005 earthquake, a large earthquake (M ∼8.5) occurred in 1861 and a somewhat smaller (M 7.6) earthquake occurred in 1907 (Kanamori et al. 2010;Fujino et al. 2014). A persistent barrier was inferred beneath Simeulue Island, which controlled the source extent and slip distributions of the 2004 and 2005 events, as well as previous earthquakes (Meltzner et al. 2012;Morgan et al. 2017).
The tsunami waveforms recorded in the far-field have been known to have problems of delayed arrivals and reversed polarity of the first wave relative to the simulated waves based on linear long-wave theory and were not used for tsunami waveform inversions. Recently, Watada et al. (2014) developed a method to correct the phase of tsunami waveforms computed by conventional linear long waves to take into consideration the effects of elasticity of the solid earth, the compressibility of sea water and gravitational potential variation. The phase-correction method was applied to the inversions of subfault slips of the 2010 Chile (Maule) earthquake (Yoshimoto et al. 2016), the 2011 Tohoku earthquake (Ho et al. 2017) and the 1960 Chile (Valdivia) earthquake (Ho et al. 2019), and successfully revealed slip distributions which are consistent with those estimated with other geodetic/geophysical data.
In this study, we invert the tsunami waveform data recorded at TGs and far-field OBPGs to estimate the slip distribution on the fault for the 2005 Nias earthquake and compare the far-field synthetic tsunami waveforms with observed data recorded on OBPGs near the Syowa Base. We also invert local geodetic data, coastal land level changes from ASTER images, sea level changes inferred from coral microatolls (Briggs et al. 2006) and GPS data (Kreemer et al. 2006) separately and jointly with tsunami data, and compare the slip distribution results.

T S U N A M I WAV E F O R M A N D G E O D E T I C DATA
The tsunamis generated from the 2005 Nias earthquake propagated in the Indian Ocean and were recorded at TGs and OBPGs (Fig. 2). Most of the TG data were available from an open website of the University of Hawaii Sea-Level Center (UHSLC), with a digital sampling rate of 2 or 4 min. The Cocos TG data (1-min sampling) were provided by the National Tide Center, Australia. For Sibolga TG data, which were the only TG data available at the west coast of Sumatra, we digitized the tide-signal waveform from Aydan et al. (2005) with varying sampling intervals of ∼3-6 min. We refer to the station list from UNESCO/IOC (http://www.ioc-sealevelmoni toring.org/) for the TG locations.
The 2005 tsunami was also recorded at two far-field OBPG stations, near (BPG1) and off-shore (BPG2) the Syowa Base in Antarctica, which also recorded the 2004 Sumatra-Andaman earthquake tsunami (Nawa et al. 2007). The digital data of BPG1 were downloaded from the Japan Coast Guard website, and the BPG2 data were provided by the National Institute of Polar Research. The sampling intervals were 30 s and 1 min for BPG1 and BPG2 data, respectively.
To remove the ocean tide signals from the TG and OBPG data, we fitted a polynomial function to each data using the GMT command 'trend1d'. The residuals of the original records from the fitted functions were used as observed tsunami waveforms. The observed waveforms (Fig. 3a) showed maximum amplitudes of 0.1 to 0.3 m at most of the TG stations. At Sibolga TG, the trough and crest heights were approximately −1 and 1 m, respectively. At OBPGs in Antarctica (Fig. 3c), the amplitudes were significantly smaller than those at TGs, with the maximum values of 0.1 m at BPG1 and 0.02 m at BPG2. We used the tsunami waveforms observed at TGs and OBPGs for a tsunami waveform inversion.
Local geodetic data were obtained from ASTER satellite images, microatoll level change measurements (Briggs et al. 2006) and GPS measurements (Kreemer et al. 2006) before and after the main shock. We used level changes at 18 sites from ASTER, 96 sites from microatoll and 54 displacement data (three components at 18 stations) from GPS measurements (Figs 4a and e).

Subfault model and tsunami computations
To estimate the slip distribution on the fault, we divided the fault plane of the 2005 Nias earthquake into eight subfaults (Table 1, and Figs 4 and 5). The subfaults were placed to cover the 1-d aftershock area (Fig. 5). We used the subfault size and configuration similar to those of the 22-subfault model of the 2004 Sumatra-Andaman earthquake . Each subfault has a length of 100 km along strike and a width of 100 km along dip. The southernmost subfaults of the 2004 earthquake (1 and 2 in ) and the northernmost subfaults (1A and 1B) of the 2005 earthquake are overlapped. Subfaults 2A-4A are aligned with the trench axis, changing the strike angles. Deeper subfaults 2B-4B are set following the shallower subfaults with the same strike angles. We referred to the Slab1 and Slab2 models (Hayes et al. , 2018 in order to set top depths and dip angles (Table 1) for the shallow and deep subfaults. The bottom edge depths are approximately 20 and 54 km for the shallow and deep subfaults, respectively, which are consistent with the depth contours of the Slab models. The rake angle of 107 • was taken from the focal mechanism of the W-phase solution (M w 8.6) by USGS.
We computed horizontal and vertical displacements on the ground surface at the geodetic observation points from each rectangular subfault with a unit slip of 1 m using equations of Okada (1985). For the tsunami initial conditions, the displacements at seafloor were first computed for a coarser grid data of 2 arcmin. We included the effects of horizontal displacement on a steep bathymetric slope (Tanioka & Satake 1996) to compute the sea surface displacement. We then resampled the 2 arcmin sea surface height distribution into the 24 arcsec grid data for the initial condition of the tsunami computation.
We assumed the simultaneous rupture of all the subfaults. As discussed in , the far-field tsunami waveform data are not sensitive to the rupture velocity. The total fault length is 400 km and the main-shock epicentre is located almost at the centre of the source area (Fig. 5). The rupture delay time from the epicentre to the edges of each subfault would be less than 100 s, and thus, the finite rupture effect can be negligible because the sampling interval is 60 s for the tsunami waveforms.
We calculated the tsunami propagation from each subfault to the TG and OBPG stations. The tsunami computation area ranges 5 • E-125 • E and 71 • S-25 • N (Fig. 2). The 30 arcsec bathymetry grid data from GEBCO 2014 (Weatherall et al. 2015) were resampled with a 24 arcsec grid interval, and thus, the numbers of computation grid points were 18 000 and 14 400 in the longitude and latitude directions, respectively. We also used the latest global bathymetry data, GEBCO 2019 (GEBCO Compilation Group 2019), and compared the results in the Supporting Information. We numerically solved the linear shallow-water equations (Satake 1995) in the spherical coordinate system for 20-hr tsunami propagation and set up a time interval of 1 s to satisfy the stability condition. The tsunami initial condition was given with a rise time of 60 s. The computation time was approximately 2 hr, with GPGPU (TESLA P100 16GB and CUDA 8.0), as used in Satake et al. (2017).

Phase correction of green's functions
We applied the phase-correction method of Watada et al. (2014) to account for the tsunami arrival time delay and the initial-phase reversal due to the elastic and gravitational coupling effect between the solid earth and the ocean. Before applying the phase correction, the calculated tsunami waveforms at nearby stations such as Sibolga were shifted vertically to remove the offset due to crustal deformation. This pre-correction of Green's function was negligible for the  (Briggs et al. 2006) and those calculated from (b) tsunami waveform inversion, (c) geodetic inversion and (d) joint inversion. Bottom row figures are the (e) observed data (red arrows) from GPS measurements (Kreemer et al. 2006) and the calculated displacements (blue arrows) from (f) tsunami inversion, (g) geodetic inversion and (h) joint inversion. The red and blue stars are the epicentres of the 2004 Sumatra-Andaman earthquake and the 2005 Nias earthquake, respectively. other far-field stations. The time-domain tsunami waveforms were then converted to the frequency domain by using the fast Fourier transform and the phases were corrected for the differences corresponding to the epicentral distance from the subfault location to the stations. We set the reference ocean depth to 4 km following Watada et al. (2014) to calculate the phase differences. The spectra were finally transformed back to the time domain. We used the phase-difference table of Ho et al. (2017) to also include the effect of stratified ocean water layers. The phase-corrected Green's functions were resampled as 1-min data, the same sampling rate as the observed tsunami waveforms.

Tsunami and geodetic inversions
For the tsunami waveform inversion, the method of Fujii & Satake (2007) was applied to simultaneously estimate the slip amount for each subfault and its error by using the non-negative leastsquares method (Lawson & Hanson 1974) and the delete-half jackknife method (Tichelaar & Ruff 1989), respectively. For the tsunami waveform data, both TG and OBPG data with different weights were used. The OBPG data weights were 10 times larger than the TG data weights because the amplitudes of the OBPG data were smaller than those of the TG data by an order of magnitude. We found the following sensitivity of different weights for the tsunami data. When we give larger weights to TG (e.g. 2:10 for TG:OBPG), the estimated slip in the northeastern subfault (1B) becomes larger, while, when we give larger weights to OBPG (e.g. 1:20 for TG:OBPG), the northeastern slip on subfault (1B) becomes smaller and the large deep slip on subfault (3B) becomes larger. We also performed geodetic data inversion and joint inversion of the tsunami and geodetic data. For the geodetic and joint inversions, we set weight parameters between the tsunami waveform and geodetic data following the method described by Satake (1993). First, we calculated the norm of both data sets, then set the weights to make the amplitude scales of both data similar. For the geodetic inversion, we set three different weights for ASTER data, coral data (Briggs et al. 2006) and GPS data (Kreemer et al. 2006). We started to test the weights from 17:1:2 (inverse number of data norm ratio for ASTER, coral and GPS data, respectively) successively reducing the weight of the ASTER data and eventually found that the weights of 2:1:2 for ASTER, coral and GPS data lead to better results reproducing the observed geodetic data, which show the better values of chi-square and variance reduction than the other weights. When we give a larger weight to ASTER data (e.g. 10:1:2 for ASTER, coral and GPS data) the amount of slip in the deep subfault (2B) becomes smaller and makes the reproducibility of the observed geodetic data worse. This is because most of the ASTER data are located in a narrow area inside the subfault (2B) whereas the coral and GPS data are widely located on subfaults (1B, 2B, 3B and 4B). For the joint inversion, we calculated the norm of tsunami data (TG and OBPG) and geodetic data (ASTER, coral and GPS) in the same way, and found that the ratio of geodetic/tsunami data norm was approximately 5. Consequently, we weighted each data set for the joint inversion as follows: TG data: 5, OBPG data: 50, ASTER data: 2, coral data: 1 and GPS data: 2. Before applying the inversions described above to the actual observed data, we performed inversion tests for different data sets synthesized from a checkerboard slip distribution (Supporting Information Fig. S1). We confirmed that tsunami and geodetic data have good resolutions for shallow and deep subfaults, respectively, and that by jointly inverting these data types the solutions achieved good resolution at both shallow and deep subfaults. We did not employ any smoothing parameters in all inversions. Tsunami data have good enough resolution in this size of subfault (100 km × 100 km), especially for shallower regions as shown in the checkerboard inversion tests.

The 2005 tsunami source model
The slip distributions estimated from different data sets, that is, from the tsunami waveform inversion, geodetic inversion and joint inversion, all showed a common feature: large slips concentrate in the deeper subfaults (Figs 5a-c). This is also consistent with the USGS finite-fault model. The largest slip of 7.3-8.6 m was located in the southeast subfault (3B) of the main shock. Slips smaller than the deep slip were estimated in the shallower subfaults, which may correspond to the source area of the previous 1907 earthquake (Kanamori et al. 2010). This complementary relationship of slips of two earthquakes between the overlapping faults is similar to the relationship of the 2007 Bengkulu earthquake (deep slips) and the following 2010 Mentawai tsunami earthquake (shallow slips) as discussed in Satake et al. (2013; Fig. 1). The range of the slip on the northernmost shallow subfault (1A) was 0.2-0.9 m ( Table 1). The slip on this subfault was also estimated to be zero or very small from the tsunami and joint inversions of the 2004 earthquake . For the deeper subfault (1B), the 2005 slip was 1.5-2.2 m, larger than the 1A, but significantly smaller than the joint inversion results (4.3-19.6 m) using TG and satellite altimetry data of the 2004 earthquake.
The slip distribution estimated in this study is very similar to the slip models of previous studies from local geodetic data (Briggs et al. 2006;Kreemer et al. 2006;Qiu et al. 2019). The calculated displacements using the estimated fault slips obtained by the tsunami waveform inversion and the joint inversion reproduced the measurement data very well and the uplift and subsidence patterns were reasonably well explained (Figs 4 and 5d-f, Supporting Information Fig. S2). The variance reductions between observed and modelled displacements yield 0.38, 0.48 and 0.46 from the inversion results of the tsunami, geodetic and joint inversions, respectively. Since we adopted a larger size of subfaults, our tsunami and geodetic data inversions do not have enough resolution to reveal local slip and deformation patterns. For example, the northern patch of peak co-seismic slip seen in the other slip distribution models (Briggs et al. 2006;Kreemer et al. 2006;Qiu et al. 2019) and zero-to-small uplifts located at the northern tip of Simeulue Island, which are well resolved by Briggs et al. (2006), are not clear features in our models.
For the tsunami waveforms at most of the TG stations (Figs 3  and 6), the first wave cycles were well reproduced, although the maximum amplitude was underestimated in Hanimadhoo. We did not use the Colombo and Port Louis data for the inversion because of the discrepancy in arrival time and lack of recorded signal of the first waves. For Sibolga, we only used the first negative or receding waves because of the large discrepancy in the later positive waves. While the phase correction did not produce a significant effect for the TG stations in the Indian Ocean (Figs 3a and b), the synthetic tsunami waveforms reproduced the observation slightly better: The variance reductions between observed and calculated TG data with the phase correction improved by 5 and 10 per cent in the tsunami waveform and joint inversions, respectively, compared with those without the phase correction. The synthetic tsunami waveforms from the geodetic inversion almost explained the observed waveforms (Fig. 6), although they overestimated those observed at some stations and also showed faster arrivals of the first positive waves due to the large slip on the shallow subfault (3A).

Inversion using different bathymetry data
Bathymetry data are essential for tsunami propagation simulations which provide Green's functions used for the tsunami data inversions. We also performed inversions using the latest global bathymetry data from GEBCO 2019 (GEBCO Compilation Group 2019) to calculate Green's functions for the tsunami (Supporting Information Table S1 and Figs S3 and S4). The tsunami and joint inversion results showed slip distributions similar to the previous results (Figs 5a and c, and Supporting Information Fig. S3). We found that the largest slip on the deep subfault (3B) was robustly estimated, despite the use of different bathymetry data. At Hanimaadhoo station, the reproducibility of the tsunami waveform became better than using the GEBCO 2014 data, while at Pointe La Rue it worsened (Figs 3 and 6, and Supporting Information Fig.  S4). The GEBCO 2019 bathymetry depth contours around Hanimaadhoo show shallower depth and substantially differ from those of the GEBCO 2014. Around Pointe La Rue we do not see much difference in the depth contours between the two sets of bathymetry data, and therefore the faster tsunami arrivals must be caused by the difference in other wider area including the tsunami propagation path to the station from the source.

Far-field tsunami waveforms
We compared the far-field tsunami waveforms calculated from the estimated slip distributions (Figs 3c and d, and 6b) to evaluate the source models. The observed tsunami waveforms were well reproduced at far offshore OBPG (BPG2) that is located at a depth of approximately 4.5 km, suggesting that our slip models of the 2005 earthquake are reasonably good. The negative wave with a reversed polarity before the first positive wave is not as clear in the observed waveform, but is clearly seen in the calculated waveforms with phase correction. At Syowa BPG1, located approximately 5 km offshore at a depth of approximately 200 m, the calculated waveforms did not reproduce the observed waveform well, especially for the short period. This may be due to the poor quality of the bathymetry data and coastal shape data around the OBPG station.

C O N C L U S I O N S
We estimated the slip distribution of the 2005 Nias earthquake (M w 8.6) using the phase-corrected Green's functions (linear long wave), including the elastic and gravitational coupling effect between the solid earth and the ocean. While the phase correction did not produce a significant effect for the TG stations in the Indian Ocean, the synthetic tsunami waveforms with the phase correction performed slightly better in terms of reproducing the observations. For the 2005 earthquake, deeper slips on the plate boundary were revealed from the observed tsunami waveforms at TGs and far-field OBPGs. The deep slips contributed to the generation of smaller tsunamis expected from such a large magnitude (M w 8.6) subduction zone earthquake. The obtained slip distribution was nearly identical to that obtained from local GPS and coastal uplift/subsidence data, demonstrating that the near-and far-field tsunami waveform analysis can reliably estimate the fault slip distribution of large earthquakes with a spatial resolution comparable to that of the local geodetic analysis. The obtained slip distribution reproduced well the OBPG tsunami waveforms at the Syowa Base in Antarctica, which is located at distances of more than 8800 km.

A C K N O W L E D G E M E N T S
We thank the editor and three reviewers for their valuable comments to improve the manuscript. All of the figures were generated using Generic Mapping Tools (Wessel & Smith 1998). Observed data at an OBPG off-shore the Syowa Base were provided by the National Institute of Polar Research (NiPR). This work was supported by JST J-RAPID grant number JPMJJR1805 and by JSPS KAKENHI grant numbers JP16H01838 and 19K04034.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
Table S1. Subfault parameters for the 2005 Nias earthquake, and the estimated slip amounts, seismic moment (Mo) and moment magnitude (M w ) from different inversions using GEBCO 2019 bathymetry data to calculate tsunami waveforms. Figure S1. Results of inversion test with a checkerboard slip distribution. (a) Target slip distribution to be inverted. The slip amounts are 10 m (red) and 0 m (white). (b) Slips inverted from synthetic tsunami waveform data. (c) Slips inverted from synthetic geodetic data. (d) Slips inverted jointly from tsunami waveform and geodetic data. Figure S2. Comparison of the observed geodetic data and calculated ones from different inversions using (a) tsunami data, (b) geodetic data and (c) joint data. Site numbers correspond to the location points of the original measurements listed in tables S1 and S2 of Briggs et al. (2006) for ASTER and coral data, and table 1 of Kreemer et al. (2006) for GPS three-component data. Red and blue lines indicate the observed and calculated displacement data, respectively. Black (positive values) and white (negative values) bars are residuals of the observed and calculated data. Figure S3. Inversion results using GEBCO 2019 bathymetry data to calculate tsunami waveforms. (a) Slip distribution of tsunami inversion. (b) Slip distribution of joint inversion of tsunami and geodetic data. Note that the slip distribution inverted from geodetic data should be the same as Fig. 5(b) because bathymetry data affect only tsunami waveforms. Figure S4. Comparison of observed and calculated tsunami waveforms using GEBCO 2019 bathymetry data to calculate tsunami waveforms at TGs and OBPGs (a) from tsunami waveform inversion, (b) geodetic inversion and (c) joint inversion. Thick lines in black and grey are observed waveforms, used and not used for the inversions. Red and blue lines are calculated waveforms with and without the phase correction, respectively. Thick lines in red are the calculated waveforms used for the inversions. The grey bars under the waveforms show the time windows used for inversions.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.