A long source area of the 1906 Colombia–Ecuador earthquake estimated from observed tsunami waveforms

The 1906 Colombia–Ecuador earthquake induced both strong seismic motions and a tsunami, the most destructive earthquake in the history of the Colombia–Ecuador subduction zone. The tsunami propagated across the Pacific Ocean, and its waveforms were observed at tide gauge stations in countries including Panama, Japan, and the USA. This study conducted slip inverse analysis for the 1906 earthquake using these waveforms. A digital dataset of observed tsunami waveforms at the Naos Island (Panama) and Honolulu (USA) tide gauge stations, where the tsunami was clearly observed, was first produced by consulting documents. Next, the two waveforms were applied in an inverse analysis as the target waveform. The results of this analysis indicated that the moment magnitude of the 1906 earthquake ranged from 8.3 to 8.6. Moreover, the dominant slip occurred in the northern part of the assumed source region near the coast of Colombia, where little significant seismicity has occurred, rather than in the southern part. The results also indicated that the source area, with significant slip, covered a long distance, including the southern, central, and northern parts of the region.


Introduction
The Colombia-Ecuador subduction zone is one of many active global subduction zones from which large earthquakes can be generated. Past earthquakes with moment magnitudes (M w ) of more than 7.0 occurred in this region in 1906in , 1942in , 1958in , 1979in (Duda 1965Kanamori 1977;Mendoza and Dewey 1984;Sennson and Beck 1996;Ye et al. 2016;Heidarzadeh et al. 2017). The earthquake properties of these events have been investigated by many previous studies; the most destructive (or largest) earthquake is considered to be the 1906 event, which had a magnitude greater than 8.0 (Gutenberg and Richter 1954;Scholz and Campos 2012). However, estimation of the detailed source process was limited because there were few quantitative datasets to specify it.
Inverse analysis of an earthquake slip distribution using observed tsunami waveforms has sometimes been conducted in cases where an earthquake triggered a tsunami (Tanioka et al. 2004). Because the 1906 earthquake induced not only strong seismic motion around the source region, but also a tsunami, sea surface fluctuations were observed by local eyewitnesses and by tide gages in some countries that faced the Pacific Ocean (Honda et al. 1908;Soloviev and Go 1975). Source models for the 1906 earthquake have been assumed/estimated by recent studies (Mayorga and Sánchez 2016;Yoshimoto et al. 2017). In particular, a nonuniform slip distribution of the earthquake was estimated by Yoshimoto et al. (2017), who conducted a slip inverse analysis using some tsunami waveforms observed at tide gauge stations. However, they excluded the tsunami waveforms observed at Naos Island (Panama) from their analysis, even though the tsunami was clearly observed there. The waveforms of Naos Island, located close to the source region, should be correctly considered in the inverse analysis because more beneficial information should be contained in waveforms obtained nearer to the observation site. Furthermore, although their source model was able to reproduce the tsunami waveforms observed in Honolulu (USA), the consistency of the model with previous studies had not been confirmed.
Therefore, the primary objective of this study was to estimate the slip distribution of the 1906 Colombia-Ecuador earthquake more accurately using the tsunami waveform data observed at appropriate tide gage stations and to compare this analysis with the results of previous studies including Yoshimoto et al. (2017) and the seismic history of the Ecuador-Colombia subduction zone.

Data and methods
Some tide gauge records, including tsunami fluctuations, were obtained from sites near and far from the source region, such as in Panama, Japan, and the USA. For example, tide gauge stations at Honolulu captured the tsunami with heights of several tens of centimeters (National Geophysical Data Center/World Data Service). However, the details of the tsunami arrival and fluctuations at some of these tide gauge stations were not clearly visible because chaotic waveforms with noise comparable to the observed wave height were mixed into the observed waveforms. Thus, we decided to use only the tide gauge data from Honolulu and Naos Island ( Fig. 1), where the tsunami was clearly visible in the data, as transcribed by Honda et al. (1908) and Soloviev and Go (1975).
First, these tsunami waveforms were estimated from the brightness matrix of images obtained through the digitization and pixelization of the documents of Honda et al. (1908) and Soloviev and Go (1975). The time resolutions of the estimated data from Honolulu and Naos Island were expected to be 1.9 and 1.6 min, respectively, and height resolutions were expected to be 0.8 and 1.5 cm, respectively, in view of the relationship between the data scale and pixel resolution of the matrix. Time series datasets with an interval of 0.5 min were made based on the estimated data. Fourier analysis with a band pass filter of 10-120 min was applied to the datasets at both stations to estimate the water surface elevation change due to the tsunami. The waveforms estimated through the above process were then determined to be the observed waveforms of the tsunami.
The inverse analysis of this study estimated the constrained slip distribution in space using a smoothing factor, as described in Tanioka et al. (2004). The intensity of the smoothing factor for the analysis was determined with a trial-and-error step. Furthermore, the weight coefficients between the data at the Honolulu and Naos Island stations were set to 1.5 and 1 for the analysis because the observed tsunami height at Honolulu was significantly smaller than that at Naos Island.
The earthquake was first assumed to consist of 45 subfault zones (each with a length and width of 50 km) that covered the entire source area, as suggested by Kelleher (1972) (Fig. 1). Strike/dip angles and the depth of each subfault were assumed based on slab depths and their coordinates (Bird 2003;Hayes et al. 2012). A rake angle was assumed and fixed for all subfaults based on a subducting angle to the trench (Chlieh et al. 2014). The vertical displacement of the ocean floor due to each subfault was first estimated based on Okada (1985), and sea surface deformation was assumed to be equal to that displacement. Furthermore, the effect of horizontal displacement within the bathymetry on the sea surface deformation was considered based on Tanioka and Satake (1996) using computation bathymetry data with a 90 arc-s spatial resolution obtained by resampling the General Bathymetric Chart of the Oceans data with a 30 arc-s resolution (Smith and Sandwell 1997;Becker et al. 2009). The propagation of each sea surface deformation generated by each subfault was then simulated with a time step of 5 s to specify the tsunami waveforms at each tide gauge station. The occurrence time of the earthquake was assumed to be 15:30 (UTC) on January 31, 1906, based on Soloviev and Go (1975), and the rupture/rise times of the earthquake were not taken into consideration. For the Honolulu and Naos Island stations, the propagation computations were conducted based on the linear long wave model, which considers no dispersion, and the linear Boussinesq model, which considers high frequency dispersion. Both the high and the low frequency dispersion effects are expected to significantly influence the waveforms at Honolulu. Thus, all dispersion effects are introduced into the waveforms simulated by the linear long wave model based on the phase correction method developed by Watada et al. (2014). In contrast, only the high frequency dispersion effect is expected to influence the waveforms of Naos Island with respect to distance relations. Furthermore, although the phase collection method can yield all dispersion effects, some assumptions are required, such as propagation distance. Therefore, the linear Boussinesq model is more suitable and useful than the phase collection method for estimating the waveforms of Naos Island. The governing equations for the propagation computations are as follows (Goto 1991): where η is the water surface elevation, M (= uh) and N (= vh) are the flux in the λ and θ directions, u and v are the velocity in the λ and θ directions, h is the still water depth, g is the acceleration due to gravity, R is the radius of the Earth, f is the Coriolis parameter, λ is the latitude, θ is the longitude, and t is the time. The linear long wave model is obtained using the above equations with an assumption that F1 and F2 are zero; the linear Boussinesq model is obtained when the above equations are fully solved. The linear long wave and the Boussinesq models were discretized based on the concepts of Goto et al. (1997) and Saito et al. (2014); however, the latitude dependence of grid size should be taken into consideration.
Here, the phase correction method requires a propagation distance and average water depth during tsunami propagation. The distance was assumed to be the arc distance along the great-circle path through the subfault location and the tide gauge station, as defined by Watada et al. (2014), who used a depth of 4 km as the average water depth during the propagation. However, arrival times of the waveforms simulated by the linear long wave model at the Honolulu station assumed a range of approximately 12-13 h. When the tsunami propagated following the linear long wave theory over distance and time, the equal average depth was comparable to 1.5-1.7 km with respect to the dispersion relation in the long wave condition. Therefore, the equal average depth for the phase correction was estimated and applied based on the relationship between the simulated arrival times of waveforms and the propagation distance using the dispersion relation in the long wave condition. Furthermore, modified dispersion relation curves corresponding to the equal average water depth were estimated by linear interpolation based on the dispersion relation curves with depths of 2, 4, and 6 km shown in Watada et al. (2014). Table 1 and Fig. 1 display the assumed/estimated subfault parameters and the estimated slip distribution, while Fig. 2 displays comparisons of the tsunami waveforms between the observations and estimations. As shown in Fig. 2, the estimated first waveforms show good agreement with the observed waveforms at both stations. The observed waveforms at Naos Island include two pulses during the first tsunami, but this was not reproduced in the estimated waveforms. It may be presumed that the assumed size of each subfault was overly large for reproduction of these pulses or that evolution of the slip distribution with time was necessary. A smaller subfault size, such as 20 or 30 km, may be useful for reproducing the pulses. After the first tsunami ended, the estimated waveforms did not match the observed waveforms. The fluctuation after the first tsunami could not be estimated sufficiently well because the nonlinearity of the tsunami and the microtopography near the coasts significantly influenced the waveforms. In particular, Naos Island is located in the Gulf of Panama where the water depth is very shallow, less than 100 m. Therefore, the waveforms estimated by the linear theory tended to significantly mismatch the observed waveforms following the first tsunami. Figure 2 includes the waveforms due to the subfault zones A01, A08, and A15 located near the southernmost, central, and northernmost sites along the trench (Fig. 1), respectively. As shown in Fig. 2b, the observed tsunami waveforms at Naos Island are very sensitive to the source location; the arrivals of the waveforms from the central and southern parts are significantly delayed compared with the observed arrival times. These results indicate the necessity of slip occurring in the northern region and support the results of the estimated slip distribution. On the other hand, an initial rise of the waveform in Honolulu seems to be insufficiently reproduced without the slip in the southern region, as shown in Fig. 2a. Thus, both northern and southern slips are necessary for proper interpretation of this event.

Results
The depth source parameters shown in Table 1 vary within a range of about 0-40 km. Taking these depths into account, the average rigidity of the source area that induced the earthquake would assume a range from about 2 × 10 10 to 5 × 10 10 N/m 2 , based on the Preliminary Reference Earth Model (PREM) by Dziewonski and Anderson (1981), and the M w of the source model was assumed to be 8.3-8.6. Many previous studies on the 1906 earthquake have attempted to estimate its magnitude based on various data. For example, Gutenberg and Richter (1954), Kelleher (1972), Abe (1979), Kanamori (1977), and Yoshimoto et al. (2017) mentioned/estimated that the magnitude of the 1906 earthquake was 8. 6, 8.9, 8.7, 8.8, and 8.4, respectively. Okal (1992) pointed out that the seismic moment for the earthquake would not be larger than 6 × 10 21 Nm (equivalent to M w < 8.4-8.5). The estimated magnitude of this study was consistent with that of Gutenberg and Richter (1954), Okal (1992), and Yoshimoto et al. (2017). However, further studies should be conducted to investigate why the estimations by Kelleher (1972), Kanamori (1977), and Abe (1979) were higher than a magnitude of 8.7.

Discussion
First, in order to examine the resolution of the inverse analysis as the result, an inverse analysis, with same conditions as described in the previous section using waveforms produced by a slip distribution as shown in Fig. 3a, was conducted and reflected the result shown in Fig. 3b. As shown in Fig. 3, the slips in the northern part were sufficiently reproduced, but those in southern part were underestimated. This is because both the Honolulu and Naos Island stations are located north of the source area, and the local slips occurring in the southern part had a relatively low sensitivity in the analysis. However, Table 1 Assumed/estimated source parameters of each subfault, whose numbers correspond to those in Fig. 1 The latitude and longitude indicate the left corner and top dip of each subfault. The depth is the distance from the seafloor to the top dip although the local slips in the southern part were estimated to be smaller than the assumed slips, the estimated source length along the strike direction was almost same with the assumed source length. Therefore, the slips in the southern part shown in Fig. 1 and Table 1 could have been underestimated, but the slip area should be estimated sufficiently.

No. Latitude (°N) Longitude (°W) Length (km) Width (km) Depth (km) Strike (°) Dip (°) Rake (°) Slip (m)
Next, we compared the source model of this study with that of Yoshimoto et al. (2017), who obtained a discontinuous and large local slip (comparable to 30 m) in the confined northern part of the source region when the observed tsunami waveforms of Naos Island were used for their inverse analysis. However, as shown in Fig. 1, this study's analysis estimated a continuous slip Fig. 2 Comparison of the estimated and observed waveforms at a Honolulu and b Naos Island. A01, A08, and A15 are slip locations defined by the grid shown in Fig. 1  Fig. 3 Result of the resolution test. a Assumed slip distribution that consists of 3-and 6-m slips and b estimated slip distribution distribution over the whole source region, including the northern part, even when data in Naos Island were used. The most important difference between the two models is the area where the dominant slip occurred. This study showed that the dominant slip occurred in the northern part of the source region, but Yoshimoto et al. (2017) concluded that it occurred in the southern part and that no slip occurred in the northern part. Kelleher (1972) provided important information based on a report that permanent coastal uplift occurred at Manta (− 0.59°N) and Buenaventura (3.54°N) and that a submarine cable between Buenaventura and Panama was broken. This cable would not have been destroyed without a large displacement due to an earthquake or  Ekström et al. 2012); mechanisms are shown only when M w is greater than 7.0. The blue circles indicate the epicenter locations of the 1942 and 1958 earthquakes (surface magnitudes of 7.9 and 7.8, respectively) based on Mendoza and Dewey (1984). The vectors indicate the assumed rupture length of the major historical earthquakes based on Kelleher (1972) and Kanamori and McNally (1982) large tsunami in the northern part of the source region. Soloviev and Go (1975) mentioned that the tsunami induced by the 1906 earthquake was most destructive along the low-lying coasts of Colombia from Rioverde to Micay (Fig. 1), where this study estimated the most significant slips to be. Furthermore, as mentioned in the previous section, the waveforms of Naos Island cannot be reproduced well without a dominant slip in the northern part. Together, these factors support the results of this study's source estimation in contrast to those of Yoshimoto et al. (2017). Figure 4 displays the estimated slip distribution with its epicenter, the assumed source length for the historical major earthquakes in the area, and the seismicity around the Colombia-Ecuador subduction zone. According to previous studies, the Nazca Plate is subducting beneath the South American Plate along the Colombia-Ecuador subduction zone with a convergence rate of about 5-8 cm/year (Chlieh et al. 2014;Scholz and Campos 2012;Kanamori and McNally 1982). The plate interface of the 1979 earthquake source area around its epicenter, which was considered to be ruptured by the 1906 earthquake, has been coupled with 60-80% rates (Chlieh et al. 2014). If we assume that all accumulated slip deficits in the 1979 source area were released by the 1906 earthquake, the plate interface of the 1979 earthquake source region should have been able to accumulate a slip deficit of approximately 2.4-3.2 m based on the elapsed time between the 1906 and 1979 earthquakes, assuming a convergence rate of the Nazca Plate of 5.5 cm/year (Chlieh et al. 2014). The average slip amount of the 1979 earthquake was estimated as 2.1-5.2 m based on the rupture area and the seismic moment, as suggested by Kanamori and McNally (1982) and Kanamori and Given (1981), when the average rigidity assumed a range of 2.0 × 10 10 -5.0 × 10 10 N/m 2 . These estimations support the idea that the 1906 earthquake ruptured the source area of the 1979 earthquake because the source area of the latter could accumulate enough slip deficits for the resulting earthquake. These results also imply that the source area of the 1979 earthquake could accumulate a slip deficit of at least 2 m. As shown in Fig. 4, the estimated slip of the 1906 earthquake around the source area of the 1979 earthquake was on the same order of magnitude of that of the 1979 earthquake. Therefore, the results of this study are also consistent with other earthquake activity in this region.
As shown in Fig. 4, the 1979 earthquake re-ruptured the plate interface originally ruptured by the 1906 earthquake with significant slip amounts. However, the dominant slip area of the 1906 earthquake, located in the northern part of the region, was not ruptured by the 1979 earthquake. Furthermore, few significant earthquakes with M w > 4.0 occurred after 1971 in the large slip area of the 1906 earthquake. According to the results of this study, therefore, the 1906 earthquake ruptured the plate interface with a maximum slip of approximately 6 m; thus, the plate interface should be able to accumulate a slip deficit of 6 m. In other words, this area is capable of causing a significant earthquake with a slip of 6 m, but the details of plate coupling and the seismicity in the northern part are still unclear.
As shown in Fig. 4, the epicenter of the 1906 earthquake was located in the southern part of the source region, suggesting that the rupture started there. Because this study estimated the dominant slip area to have been in the northern part of the region, the rupture must have propagated from southwest to northeast, as suggested by Kelleher (1972). Furthermore, the slip area of the 1906 earthquake seems to cover not only the 1979 earthquake area, but also the 1942 and 1958 earthquake areas, located in the central and southern parts of the region, respectively. These results imply that the rupture area of the 1906 earthquake may have involved those of the previous three earthquakes; in other words, the 1906 earthquake might have been a consolidated-type earthquake that simultaneously ruptured all these areas, similar to the case of the Nankai Trough Earthquake in Japan (Furumura et al. 2011).

Conclusions
This study estimated the slip distribution of the 1906 Colombia-Ecuador earthquake based on an inverse analysis using observed tsunami data and considered the results in the context of previous/expected earthquake activity. The main conclusions are summarized below: 1. A slip distribution model of the 1906 earthquake in space was newly developed. It indicated that the M w of the earthquake ranged from 8.3 to 8.6. 2. The newly defined presence of a local large slip area in the northern part of the source region is supported by the historical seismicity and the earthquake/ tsunami damage that occurred along the coast of Colombia. 3. Although the physical characteristics of the subduction zone, such as the seismicity and the interseismic plate coupling, had been well studied (e.g., Beck and Ruff 1984;White et al. 2003;Chlieh et al. 2014), further detailed investigations of the earthquake activity in the whole subduction zone (especially in the northern parts) are necessary. Further development of the slip distribution model with time dependence also remains as a goal for future work.