Precise Locating of the Great 1897 Shillong Plateau Earthquake Using Teleseismic and Regional Seismic Phase Data

The collision and continued convergence of the India plate with Eurasia has caused widespread deformation along the northern part of the Indian subcontinent. These tectonic movements formed the Himalayan arc, one of the best studied continental collision orogenic belts, also producing several megathrust earthquakes in historical and modern times (e.g., Bollinger et al., 2016; Hetényi, Cattin, et al., 2016; Bilham, 2019; Dal Zilio et al., 2021). The seismic activity in the eastern Himalaya appears spatially more widely distributed, as the deformation zone extends to and includes the Shillong Plateau (e.g., Grujic et al., 2018; Fig. 1). The Shillong Plateau is the only high‐elevation terrain in the Himalayan foreland, where one of the largest and most damaging great earthquakes have occurred, on 12 June 1897.

The 1897 Shillong Plateau–Assam earthquake is one of the best documented earthquakes of the nineteenth century. A detailed description of damage, field observations, and various other data has been compiled by Oldham (1899), and his monograph has been the source for several papers describing and modeling the local and regional effects. Distinct seismic phase arrivals were recorded at various stations, and from some stations full waveforms are available (Fig. 2). Surface waves arriving on the longer great circle (R2) were historically first recorded for this event by six seismometers across Europe. The local shaking effects were considerable, including liquefaction. The macroseismic intensities and their distribution have been evaluated: several studies have built on these reports and further data, assessing the maximum intensity at X on the Rossi–Forel intensity scale and IX on the Medvedev–Sponheuer–Kárník scale. These historical studies were used recently for further evaluation and seismic hazard assessment (e.g., Ambraseys and Bilham, 2003; Hough et al., 2005). Rupture mechanism analysis using geodetic data reveals that the stress drop implied by the rupture geometry and fault slip explain that peak ground acceleration exceeded 1 g (Bilham, 2008), which is consistent with field observations (Oldham, 1899).

Regarding the earthquake’s magnitude, several estimated values have been proposed in the literature, ranging between M 8.0 and 8.4 (Table S1). The most recent estimates are of moment magnitude 8.25 ± 0.1 (England and Bilham, 2015) and 8.1 (Pettenati et al., 2017). In contrast, there is significantly more spread regarding the mainshock’s location proposed in the literature (Fig. 3a): some of them are outside the likely epicentral area contoured by Oldham (1899). The only proposed location with an uncertainty is by Gutenberg (1956): the coordinates are the same as proposed in Rudolph (1903), with a nominal uncertainty of 1° (ca. 111 km). The corresponding area is relatively large compared to the size of the region and would allow the epicenter to be located on the Himalayan Main Frontal Thrust or on the southern edge of the Shillong Plateau (Dauki or Dapsi fault). One of our motivations is therefore to better constrain the epicenter of the 1897 earthquake using instrumental recordings, especially as at the time of the earthquake the opposite approach was used: the epicenter was determined by the analysis of macroseismic records, and research focused on wavespeeds using arrival times. Today, we have the opportunity to benefit from modern velocity models to address the question of the epicenter with more precision.

Not only the epicentral location but also the origin time (Table S1) and the geometry of the causative fault of the 1897 earthquake have been debated. Earlier studies suggest that the earthquake rupture occurred on a north‐dipping Himalayan fault propagating south of Bhutan (e.g., Oldham, 1899; Molnar, 1987; Gahalaut and Chander, 1992). However, it has later been proposed that the rupture had occurred on a hidden reverse fault (named the Oldham fault), which is ca. 110 km in length and dipping steeply to the south, away from the Himalaya (Bilham and England, 2001). In addition, one of the largest average slip values was proposed, up to 25 ± 5 m (England and Bilham, 2015).

In summary, research on the 1897 earthquake has been primarily based on damage reports and local data, and the location and origin time of the mainshock are still rather uncertain (Table S1). Here, we analyze early instrumental seismic phase arrival data registered at teleseismic and regional distances to precisely constrain the location and origin time of the 1897 earthquake. We are not aware of an earlier quantitative attempt, primarily as the data were considered to be inaccurate and/or insufficient. In our calculations, however, we consider the uncertainties of arrival times, and through an iterative approach constrain the location and origin time of the event using modern, 3D velocity models. The solution we find is the first of its kind, which we interpret in terms of regional seismotectonics.

We performed an extensive search for original data and resources to collect information about seismological records from the 1897 earthquake, including books, reports, articles in various languages, querying early seismological station lists with operational sites at the end of the nineteenth century. We also contacted observatories, libraries, and other archives in search for arrival time information and waveforms. Ultimately, we collected arrival time data from 32 stations at teleseismic distances, with 29 P‐wave and 19 S‐wave arrival times. Search for data at further 28 stations (from 15 present‐day countries on 5 continents) have not returned data, either because the operation of stations started shortly after the 1897 event, or as the records were lost, or because the recording paper was being changed right at the time of the earthquake. Original waveforms from five seismological instruments are available (see Fig. 2) but did not provide additional information beyond the reported arrival times.

The main challenge with the arrival times was their precision and accuracy. Some of the records came with information only on the minutes but without the seconds; these could not be used in the location procedure. The arrival times of all wave phases recorded at each location have been converted into Coordinated Universal Time (UTC). This was straightforward for teleseismic data, but the accuracy of the analog clocks at the time was mostly unknown; only few records provided uncertainties (from a few to 15 s). All station coordinates could be identified with sufficient precision (Tables S2, S3) so that their uncertainty contributed insignificantly to the calculations compared with the arrival time uncertainty.

A complementary set of body wave arrival time information came from regional distances: as reported by Oldham (1899), train station masters across India have noted the arrival of clearly felt waves or the stopping of clocks due to the earthquake. Oldham has already compiled this data, first by readjusting the diverse local time zones to a common reference and then by assessing the quality of data subsets from different train lines. We adopted arrival time information only from the East Indian Railway running between Calcutta and Delhi, “the busiest, and by common sense the best managed line in India,” for which Oldham (1899) already showed the approximate wave propagation speed, matching S waves (see hodograph on his plate XXXIX). We discarded stations located at <500 km distance from the Shillong Plateau to avoid any near‐field effects such as Pg–Pn crossover or crustal thickness change. The available time information is typically given in minutes (Table S3), which we consider to be the full and complete minutes on the clock; hence we assume that time uncertainties up to +60 s are possible, but no negative deviations. Finally, local uplift data (Fig. 3b) were taken from Oldham (1899).

To constrain the area of the best fitting epicenter and the related origin time, we start from an initial location and time, and calculate misfits from the observed data, as described in the following. This includes selecting phases that are coherent with that solution within uncertainty and discarding others. The location and time of origin are then updated, and the calculations are repeated. After three iterations we find convergence, meaning that the epicenter location and the origin time coincide with the computed minimum misfit map (for location) and curve (for time). Details of the individual steps were the following.

Initial location: We adopt that recently proposed by Pettenati et al. (2017) at 25.61°N, 91.42°E. This is the easternmost location proposed in the literature (Fig. 3a), close to the eastern end of the geodetically proposed subsurface Oldham fault.

Theoretical travel times: We compute for all possible locations over the area of interest—the broader area of the Shillong Plateau, 88°–94°E and 23°–29°N, at 0.1° $×$ 0.1° grid spacing—the theoretical travel time of each phase to all recording station locations (Fig. 4a). The theoretical travel times for P and S phases were computed using LLNL‐Earth3D, a recent, global, 3D velocity model (Simmons et al., 2015) that accounts for crustal thickness changes and mantle velocity anomalies along the raypaths (Fig. 4a).

Origin time misfit considerations: The travel times are subtracted from the observed arrival times for each phase. This heterogeneous set of possible origin times from all stations are compared with a range of plausible origin times based on the literature (from 11:05:00 to 11:08:00 UTC), and from the absolute value differences a cumulative L1‐misfit is computed for every assumed time at 1 s spacing at the location of the current epicenter. The updated origin time is at the minimum of this cumulative misfit curve.

Selection of coherent phases: Theoretical time versus distance curves are plotted for the current location and origin time for the distances corresponding to the stations with data. The observed arrival times are compared to these curves, considering the uncertainty of each phase. Based on our reading of the arrival time information, we assign ±15 s uncertainty for teleseismic P and S waves, and +60 s uncertainty for phases at regional distances. Phase arrivals lying inside these ranges are kept for further calculation.

Earthquake epicenter: To update the epicenter solution, we back‐projected the observed arrival times to the area of interest and considered the same uncertainties as mentioned earlier. Two maps were produced: a hit‐count map, on which all locations within uncertainty of each phase are treated equally (+1 on the grid node), and a more classical misfit map, on which absolute time misfit values are added within the uncertainty band. In the first map, we seek for the maximum value area (the highest number of phases coherent with that solution), whereas in the latter map we look for the lowest cumulative misfit area.

In summary, cumulative travel‐time misfit curves and misfit maps are used jointly to evaluate the results. Three iterations were required to find a coherent solution. The selection of fitting phases was the strongest in the first round, with only few changes in subsequent iterations.

In the final iteration, 24 phases have been found to be coherent with the earthquake’s epicenter and origin time within uncertainties. Here, we present the results using these phases: 8 P and 6 S at teleseismic, and 10 stations at regional distances—an acceptably balanced contribution of the three sets of available information. The arrival times incoherent with the solution represent inaccurate time records; these being distributed on both sides beyond the uncertainty thresholds (Fig. 5) suggests that there was no systematic bias in the selection.

The cumulative misfit curves for various origin times have been computed for two focal depths: 10 and 40 km (Fig. 4b). However, the assumed focal depths make little difference to the solution, and we cannot distinguish between them based on this analysis, especially considering the width of the minima. Overall, the minimum of cumulative misfit between observed and computed arrival times is at 11:06:46 UTC, at a cumulative misfit value of 217 s for 24 phases (ca. 9 s per phase). The origin time uncertainty is ca. 15 s, considering the width of the misfit curve at 240 s (10 s per phase; Fig. 4b). This time uncertainty is comparable to the source time function of similar size earthquakes that typically exceeds 40 s for M ∼8 events that occurred during 1976–2021 at < 50 km depth (see Data and Resources). Finally, the obtained origin time also fits four additional teleseismic observations for which no seconds information were available (Fig. 5a,b).

With the best estimate for the origin time, the arrival times of all the selected phases were back projected to verify and constrain the earthquake’s epicenter. The first approach, in which equal weight is given to all points within the time uncertainty, highlights a narrow band in the north‐central part of Oldham’s proposed epicenter limit, where all the 24 phases fit (Fig. 4c). The surrounding, still well‐fitting area is centered on Oldham’s contour and includes the Oldham fault.

The second approach shows cumulative time misfit over the map with the uncertainty bounds described earlier (Fig. 4d). The point of minimum misfit falls very close (24 km precision) to the assumed location in this iteration and thus confirms this epicenter, located at 26.0°N, 90.7°E, at the western end of the Oldham fault, at its junction with the Chedrang fault (see pink star in Fig. 4d). The accuracy is less than 54 km, considering the distance between the assumed location in the previous iteration and the final one. An area representing the origin time uncertainty is projected onto the map, and shown at 228 and 240 s cumulative (on average 9.5 and 10 s per phase) misfit values (see black dashed contours in Fig. 4d). This area is consistent with the damage reports, and its elongation axis is perpendicular to that of Oldham. The area within 228 s resembles an ellipse of ca. 0.3° semi‐minor and 0.7° semi‐major axis in longitude and latitude, comparable to the ca. 50 km location accuracy. The area at 240 s reaches the Main Frontal Thrust in the north and Dapsi and Dauki faults in the south; however, given their peripheral location and earlier studies, these faults are unlikely to be the source of the 1897 earthquake. Most previous epicenters’ estimates did not provide an uncertainty range; Gutenberg (1956) estimated it at a nominal ±1° (the best category in that paper), our result has a better accuracy (< 50 km).

The most robust results from this study are the origin time and epicenter of the 1897 Shillong Plateau mainshock (Fig. 4), constrained by joint inversion of phase arrivals at teleseismic and regional distances.

The earthquake’s origin time is not discussed extensively in earlier studies, except Oldham (1899) (see Table S1): Our origin time matches well with his origin time estimate, with only 4 s difference, which is remarkably close. Two more studies are close to our origin time but without seconds information (11:06 UTC; Ambraseys, 2000; Ambraseys and Douglas, 2004). Other previously reported origin times are either too early (e.g., 11:05:01 UTC; Milne, 1899) or too late (e.g., 11:09 UTC; Ambraseys and Bilham, 2003), and are clearly beyond the time–uncertainty range determined in our calculation.

Regarding the earthquake depth, most prior results point toward a fault reaching deep (30–50 km) beneath the Shillong Plateau (e.g., Chen and Molnar, 1990; Pettenati et al., 2017) and without surface rupture. The dataset presented here does not allow to constrain the focal depth, as the time difference for the 10 and 40 km depth solutions is only 5 s, whereas the teleseismic phase arrival uncertainty is ±15 s.

Our epicenter result lies inside Oldham’s probable epicenter boundary, which was plotted based on local damage reports. Most of the published earthquake epicenters are to the south of our solution (Fig. 3a). Yet, epicenters reported by Rudolph (1903), Gutenberg (1956), and one by Oldham (1899) are relatively close, at <40 km distance from our result. It could be that previously reported epicenters were influenced by earlier publications, by the position and orientation of the Oldham fault, and by the fact that Oldham observed the most severe surface deformation east of 91°E (Fig. 3b). Oldham thus concluded that the epicenter was required to be on the northern side of the Shillong Plateau, and our solution corresponds to this. Our epicenter is located 84 km to the west of the solution of Pettenati et al. (2017) based on macroseismic data: We believe the two solutions can coexist, as our solution points to the place from which the first seismic waves were emitted, whereas the macroseismic solution reflects the highest damage (Fig. 3).

Because of the lack of clear field evidence of surface rupture over considerable lengths, the location and geometry of the responsible fault have been debated. Two main proposals were made: a hidden, south‐dipping Oldham fault along the northern margin of the Shillong Plateau (e.g., Bilham and England, 2001) and the Dauki fault along the southern relief of the plateau (e.g., Morino et al., 2014). However, surface uplift (Fig. 3b) and strain decreasing from north to south across the plateau (England and Bilham, 2015) refute the latter proposal. Our epicenter solution is at the western end of the proposed Oldham fault—at the Chedrang fault, where up to 35 ft (10.7 m) of normal faulting has been reported (Oldham, 1899). The epicenter uncertainty assessment excludes that this earthquake occurred on another major fault and corroborates the model solution based on geodetic data (England and Bilham, 2015), with the fault plane’s northwest and upper corner matching our seismological solution within uncertainty and tightly constraining the epicenter (Fig. 4d). The location of the epicenter at the western end of the Oldham fault suggests that the earthquake have likely ruptured unilaterally toward the east along the Shillong Plateau.

Thirty‐three years later and less than one fault length to the west of the Oldham fault, the 2 July 1930 M 7.1 Dhubri earthquake (Gee, 1934; Fig. 4c,d) may have been induced by stress changes accompanying the 1897 earthquake. The 1930 event occurred on the Dhubri–Chungthang fault zone (DCF)—a dextral, mid‐ to deep‐crustal fault that causes segmentation of the Indian basement and extends beneath the Himalaya (Hetényi, Cattin, et al., 2016; Diehl et al., 2017). The Bhutan Himalaya and Shillong Plateau stress regimes are weakly connected, but large earthquake in both the regions can stress load the DCF in the west and the Kopili strike‐slip fault zone in the east (Grujic et al., 2018). The 1897 earthquake could have loaded sufficient stress on the DCF to ultimately trigger the 1930 event. These two earthquakes attest to the complex deformation regime of the region, though the absence of surfaces rupture makes this more difficult to investigate.

To locate the 1897 M 8+ earthquake, we have made an extensive search for original records of wave phase arrival times from this early‐instrumental era, leading to the teleseismic data from Europe and regional distance data from India. These were combined with 3D velocity models to iteratively constrain the origin time and epicenter of the earthquake. Arrival time uncertainty was used both to select reliable data and to assess the result’s precision.

Our results constrain the origin time at 11:06:46 UTC, with an uncertainty of ca. ±15 s. The epicenter at 26.0°N, 90.7°E is located on the northwest Shillong Plateau, at the western end of the Oldham fault at its junction with the Chedrang fault. The location’s accuracy is much better than that of previous studies, and it matches the geodetically inferred thrust plane. This corroborates the scenario for an event on the buried Oldham fault and suggests that rupture propagated eastward. The 1930 M 7.1 Dhubri earthquake on the DCF may have been triggered by the 1897 event. Further studies on deep geophysical structure and paleoseismicity could further advance the seismotectonic understanding of this area, the nucleation of and the interaction between large earthquakes, and contribute to better seismic hazard models.

This study did not use any new data. All data used in this article came from published sources listed in the references. The LLNL‐G3D velocity model is available from https://www-gs.llnl.gov/nuclear-threat-reduction/nuclear-explosion-monitoring/global-3d-seismic-tomography. The source time function data of 1976–2021 earthquakes are available at https://www.globalcmt.org/CMTcite.html. International Seismological Centre–Global Earthquake Model (ISC‐GEM) earthquake data are available at http://www.isc.ac.uk/iscgem/. All websites were last accessed in November 2021. The supplemental material for this article includes summary table of body waves’ arrival time, and previously published epicenters and origin times of the 1897 earthquake.

The authors acknowledge that there are no conflicts of interest recorded.

The authors are grateful to Remo Grolimund for his help in locating and retrieving original information on seismological stations and arrival times and for revising a draft of this article; to Graziano Ferrari and Luis Rivera for discussions on the early days of seismology; and to two anonymous reviewers for their constructive comments that improved this work. The authors acknowledge the Geosciences library of the University of Lausanne for keeping a well‐preserved copy of Oldham (1899), including pages we opened for the first time. They also kindly acknowledge support of the Swiss National Science Foundation (Grants PP00P2_157627, PP00P2_187199).