Nowcasting earthquakes in Sulawesi Island, Indonesia

Large devastating events such as earthquakes often display frequency–magnitude statistics that exhibit power-law distribution. In this study, we implement a recently developed method called earthquake nowcasting (Rundle et al. in Earth Space Sci 3: 480–486, 2016) to evaluate the current state of earthquake hazards in the seismic prone Sulawesi province, Indonesia. The nowcasting technique considers statistical behavior of small event counts between successive large earthquakes, known as natural times, to infer the seismic progression of large earthquake cycles in a defined region. To develop natural-time statistics in the Sulawesi Island, we employ four probability models, namely exponential, exponentiated exponential, gamma, and Weibull distribution. Statistical inference of natural times reveals that (i) exponential distribution has the best representation to the observed data; (ii) estimated nowcast scores (%) corresponding to M ≥ 6.5 events for 21 cities are Bau-bau (41), Bitung (70), Bone (44), Buton (39), Donggala (63), Gorontalo (49), Kendari (27), Kolaka (30), Luwuk (56), Makassar (52), Mamuju (58), Manado (70), Morowali (37), Palopo (34), Palu (62), Pare-pare (82), Polewali (61), Poso (42), Taliabu (55), Toli-toli (58), and Watampone (55); and (iii) the results are broadly stable against the changes of magnitude threshold and area of local regions. The presently revealed stationary Poissonian nature of the underlying natural-time statistics in Sulawesi brings out a key conclusion that the seismic risk is the same for all city regions despite their different levels of cycle progression realized through nowcast scores. In addition, though the earthquake potential scores of the city regions will be updated with the occurrence of each small earthquake in the respective region, the seismic risk remains the same throughout the Sulawesi Island.


Introduction
Sulawesi is one of the four Greater Sunda Islands in Indonesia which has high seismicity. Earthquake sources in these regions come from tectonic processes on land and sea, fault systems in the middle, and subduction zone in the north (Cardwell et al. 1980;Mc Caffrey & Sutardjo 1982;Mc Caffrey 1982;Walpersdorf et al. 1998;Stevens et al. 1999;Gomez et al. 2000;Hall 2009). Many destructive events in the Sulawesi province have occurred in the vicinity of the major cities that often knit together to form a global hub of economic and industrial activity (Cipta et al. 2017). In recent years, the entire island has observed a rapid growth of urban population and infrastructure development as an outflow of worldwide economic empowerment and industrial revolution. However, the fear of large damaging earthquakes looms over the Sulawesi region which is surrounded by a network of active fault lines and has many seismic nonresistant traditional structures (Widiyanto et al. 2019;Omira et al. 2019).
The seismicity in the Sulawesi area is so high that it has experienced 25 major earthquakes (M w ≥ 7.0) during a span of about 50 years . The latest of them is the 28th September 2018 M w 7.5 "supershear" Sulawesi earthquake occurring in the northern part of the Palu-Koro fault that has not been well studied before (Bao et al. 2019). This earthquake was categorized as one Open Access *Correspondence: sumanta.pasari@pilani.bits-pilani.ac.in; 86.sumanta@gmail. com 1 Birla Institute of Technology and Science Pilani, Pilani Campus, Jhunjhunu 333031, Rajasthan, India Full list of author information is available at the end of the article of the deadliest earthquakes in 2018 as it was accompanied by an unusual tsunami caused by the unexpectedly large vertical component of the fault slip (Bao et al. 2019;Ulrich et al. 2019). Liquefaction was observed in the Palu city which sits at the end of the bay surrounded by a river delta (Bradley et al. 2019). Similar to the September 2018 Sulawesi earthquake, other disastrous events may also occur along the less studied geological faults that pass through the major cities (Cipta et al. 2017). Therefore, it is imperative to evaluate the current level of seismic hazards for all populous cities in the Sulawesi province for a nationwide disaster preparation and mitigation. In light of this, the present study develops a statistical, datadriven nowcasting approach (viz. Rundle et al. 2016) to estimate the contemporary state of earthquake progression in several regions of the Sulawesi Island. As in financial nowcasting, the earthquake nowcasting technique seeks to define the current state of hazard, rather than to produce a forecast of future years. Since the basic stress variables that define the current state of the earth's crust are fundamentally unobservable (Scholz 2019), one seeks to find proxy variables that can be used to estimate the current state of hazard, thus to provide some form of useful information for policymaking (Rundle et al. , 2018(Rundle et al. , 2019ab).
The nowcasting method in seismology has evolved from the understanding that in a dynamic, non-linear, power-law behaving earthquake system (similar to typhoons, landslides, market economy crashes), the distribution of events can be broadly aligned in a frequency-magnitude spectrum where the large magnitude events are usually preceded by several smaller events (Rundle et al. , 2018(Rundle et al. , 2019ab). This method has found significant applications in economics (e.g., growth nowcast, inflation nowcast, market-stress nowcast), meteorology (e.g., thunderstorm nowcast), and political sciences (e.g., electoral campaign nowcast) where the interest is to assess the uncertain current state in view of the data of immediate past, present, or very near future. The nowcasting technique in seismology utilizes the key idea of natural times (e.g., Varotsos et al. 2011), the interevent small earthquake counts braced by successive large earthquakes in a specified region, rather than the usual calendar or clock times to elucidate the present state of earthquake system through earthquake cycle. It is built upon the assumption that the underlying seismicity statistics in local regions is embedded in a homogeneous larger region from which the seismicity statistics will be developed. This assumption seems reasonable in light of the ergodic dynamics in the statistical physics of seismicity, in which the statistics of smaller regions over longer times are considered to be similar to the statistics of broader regions over large spatial domains and longer periods (Tiampo et al. 2003;Holliday et al. 2016). Under this setup, the nowcasting method describes the seismic progression of a region in terms of earthquake potential score (EPS) measured through the cumulative probability for the current event counts (Rundle et al. , 2018Pasari, 2019b). The method has been successfully applied to several seismogenic cities such as Ankara, Dhaka, Islamabad, Kolkata, Lima, Manila, New Delhi, Taipei, Tokyo, Los Angeles and San Francisco Pasari 2019bPasari , 2020Pasari et al. 2021a, b, c;Pasari and Sharma 2020;Bhatia et al. 2018).
The method of earthquake nowcasting is conceptually different from earthquake forecasting which looks ahead in time. Nowcasting technique, unlike forecasting, is applicable for a dataset involving dependent as well as independent events . Nowcasting often enables a systematic ranking of cities as to their current exposure to the earthquake hazard, whereas forecasting deals with several probability measures for longterm planning and preparation (Rundle et al. , 2018. The nowcasting method may provide a fast and sophisticated alternative means to understand the current state of progress of a regional fault system, which is otherwise traditionally determined from direct or indirect, physical or remote observation of a tectonic stress regime (Scholz 2019).

Geology and tectonics
Split by the equator, the study area encompasses a quadrangle bounded by the geographical limits 115°-130° E and 10° S-5° N (Fig. 1). Due to the inherent setting from the ongoing trio-convergence between the north moving Australian plate, south-southeast moving Eurasia plate and the west moving Pacific plate, the Sulawesi island and its contiguous areas exhibit an exceedingly complex seismotectonic pattern ). The tectonic process, mainly the drifted collision of the three plates, started since the Mesozoic era when it experienced a spreading exposure in the northwestern part of several microcontinent fragments derived from the Australian continent (Audley-Charles et al. 1988). The relentless tectonic process in a number of episodes eventually constitutes a triple junction-a unique, unusual K-shape of the Sulawesi island with four characteristic arms (Monnier et al. 1995;Katili 1978). The convergence zone of this tri-junction accommodates a composite domain of accretionary complexes, microcontinental fragments, melange terrains, island arcs and ophiolites (Hall 2012). The stratigraphic development of Sulawesi is mainly controlled by the successive accretion from the east of oceanic and microcontinental material (Wilson et al. 2000;Bosence and Wilson 2003). As a consequence, the K-shaped Sulawesi Island is considered to be a response to the post-collision rotation of the curvatures of four arms originally being convex to being concave (Monnier et al. 1995).
Sitting in the central part of Indonesia archipelago, Sulawesi Island is one of the most enigmatic areas with its geological formation, process and structure Monier et al. 1995;Wakita et al. 1996;Villeneuve et al. 2002;Hall 1996Hall , 1998Hall ,2002Hall ,2009. Under the influence of Australia-Philippine plate convergence, the Sulawesi domain persistently collides with the Eurasian plate. The tectonic collision is accommodated by subduction at the north Sulawesi trench and by the motions along the Palu-Koro and Matano faults at the southwestern and southern domains (Walpersdrof et al. 1998;Socquet et al. 2006;Bellier et al. 2006). One of them forms the southern margin of the Celebes Sea (Walpersdorf et al. 1998). Associated with the tectonic collisions, there was development of NNW-SSE trending Palu-Koro left-lateral transcurrent fault, along which the part of eastern Sulawesi has moved northwards with regard to Fig. 1 Topography of the Sulawesi island including general tectonics and major cities western Sulawesi. Opening of pull-apart basins, such as lakes of Poso, Matano and Towuti, Palu depression and other tectonic features are evident as a result of recent transtensional movements during the Quaternary period continuing to the present time. Several disastrous events along the Palu-Koro and nearby major faults reveal that the system is probably active (Hall 2009).
Tectonically, the Sulawesi Island, due to enormous external pressure, deforms continuously. Such pressure comes from the Flores Sea in the southern part and activates Palu-Koro, Walanae and Banggai-Sula fault (Jaya and Nishikawa 2013). External pressure from Banda Sea in the eastern to central part activates Matano, Batui, Lawanoppo and Kolaka fault, whereas the pressure from Sulawesi Sea in the northern part activates North Sulawesi subduction and Gorontalo fault. Pressure from North Maluku subduction plate from northeast often causes large earthquakes and volcanoes in North Sulawesi area. These ongoing tectonic phenomena in the Sulawesi island produce four major transcurrent faults, namely the Sorong-Matano fault, Palu-Koro fault, Walanae sinistral fault, and the Gorontalo dextral fault (Watkinson et al. 2011;Watkinson and Hall 2017). The pathways of these active faults exhibit high earthquake and tsunami potential in the study region.
Several earthquakes located in the sea floor have tsunami potential. The extent of the tsunami size is often influenced by steep subsurface topography of the seabed and coast (Hamzah et al. 2000). Both southern and northern Sulawesi have witnessed disastrous tsunamis (Baeda 2011). In the southern part, stretched from Majene to Mamuju or from Palu to Toli-toli, tsunamis occur due to earthquakes in the intersection of Paternoster fault and Makassar Strait normal fault or in the intersection of Palu-Koro fault and Makassar Strait fault. In the northern part, tsunamis occur in Gorontalo, Luwuk-Banggai and Kendari-Wawoni-Buton areas due to events located in the intersections of Gorontalo fault and North Sulawesi subduction, Gorontalo fault and North Maluku subduction, and Lawanoppo fault and Wawoni thrust (Watkinson et al. 2011;Watkinson and Hall 2017). Areas that have been hit by tsunami since 1967 are Majene-Pinrang in 1967, Palu in 1968, Mamuju in 1969, Donggala in 1996, Toli-toli in 2000and Luwuk-Banggai in 1999. Exploring the nowcasting technique, which provides a surrogate way to determine the current progress of an earthquake cycle in a fault system thus appears to be an interesting task to indirectly assess earthquake/tsunami hazards in the Sulawesi Island ). However, it may be noted that in a nowcasting analysis, the entire large geographic region is assumed to be a single driven threshold system to develop nowcast scores for several local regions embedded in the large region.

Earthquake data
To develop nowcasting technique, we prepare an earthquake catalog  of instrumental data mainly from the compilation of two global catalogs, namely the Advanced National Seismic System (ANSS) comprehensive catalog (http:// www. ncedc. org/ anss/ catal og-search. html) and the International Seismological Centre (ISC) catalog (http:// www. isc. ac. uk/iscbulletin/search/catalogue/). Both ANSS and ISC (Reviewed ISC Bulletin) are reliable sources of earthquake data. To maintain the consistency of data sources, we used "magnitude author" NEIC (National Earthquake Information Center) in the ISC search. For the recent data (January to November, 2020), we use the regional Meteorological, Climatological, Geophysical Agency of Indonesia (BMKG) catalog. As the BMKG network comprising 162 broadband seismometers is locally maintained, recent events have been directly adopted from this catalog, without any possible magnitude alignment between BMKG and others. A pictorial summary of earthquake data is provided in Fig. 2, whereas the entire dataset is provided as Additional file 1.
In the analysis, we consider magnitude 4.5 as the threshold of small earthquake events (also, the magnitude completeness threshold) and 200 km as the threshold of maximum focal depth, since deep-seated oceanic earthquakes often exhibit weaker signal and are generally less damaging Pasari 2019b). During January 01, 1969 to November 21, 2020, a total of 17,082 events (4.5 ≤ M ≤ 7.9) have occurred in the study region providing 87 earthquake cycles of M ≥ 6.5 events. The associated earthquake interevent counts (natural times) will be used to develop the natural-time statistics. Recall that nowcasting analysis via natural-time statistics does not require declustering of dependent events, such as foreshocks and aftershocks ).

Procedure and results
Nowcasting method comprises three major steps in succession: tabulation of natural times (interevent counts between large events), performing statistical inference, and computation of earthquake potential score (nowcast scores). While tabulation of natural time essentially requires the description of "large" and "small" events in a specified homogenous region, the statistical inference involves probability model description, parameter estimation, and model performance analysis. Knowing the data-derived seismicity statistics of natural times, we use mathematical cumulative distribution function (CDF) of the best-fit probability distribution to calculate earthquake potential score (EPS) for a number of selected Pasari et al. Geosci. Lett. (2021) 8:27 Fig. 2 Pictorial summary of earthquake data  that are used for the present nowcasting analysis; subplots highlight several characteristics of the catalog, including epicentral distributions on map, cross section views of hypocentral depth, magnitude of completeness, and occurrence time of large earthquakes in the study region cities. These EPS scores, as "thermometer" readings, provide a snapshot of the current progression of a city in its earthquake cycle of large events.
Let A be the geographical area of a broad seismically active region and A 1 , A 2 , A 3 , · · · , A n be the geographical areas of the n number of selected cities in a way that A 1 ⊂ A, A 2 ⊂ A, A 3 ⊂ A, · · · , A n ⊂ A . Also, let M σ and M , respectively, denote the magnitude of the small and large earthquake events in a sense that "large" events have the potential to cause human death or notable destruction to society. The value of M σ is usually guided by the data-completeness threshold in the frequency-magnitude relationship (Scholz 1990;Rundle et al. 2016), whereas the value of M is decided based on the extent of damages in a city. Once M σ and M are decided, we define natural time (say, X) as the interspersed small event counts between two successive large events. It is obvious that natural time (X) exhibits randomness.
For the present study with 1969-2020 catalog, we consider M σ = 4.5 and M = 6.5 . This yields a random sample of natural times {X 1 , X 2 , · · · , X 87 } to infer underlying seismicity statistics in the Sulawesi region. While the sample range of interevent counts varies from 1 to 1034, other descriptive statistics such as mean, median, standard deviation and skewness turn out to be 192, 145, 189 and 2, respectively. Thus, the observed dataset reveals asymmetry with a skew towards right, having sample mean higher than sample median. The histogram of the observed natural time counts is shown in Fig. 3.
The nowcasting approach or the natural-time statistics is independent to the background seismicity rate (a) as long as the Gutenberg-Richter b-value (slope in frequency-magnitude relation) is close to a constant. To illustrate, let the cumulative number of small earthquake counts be N σ = 10 a−bM σ and cumulative number of large event counts be N = 10 a−bM . Then, using the The number of small earthquakes, N , that occur between two large earthquakes can then be obtained by setting N = 1 (Luginbuhl et al. 2018b). Thus, N = 10 b(M −M σ ) . As a consequence, N is independent to the background seismicity rate (a) and it scales exponentially with the difference of magnitudes M and M σ . As a consequence, it is reasonable to consider exponential distribution and its primary variants gamma, Weibull, and exponentiated exponential in developing the mathematical cumulative distribution function (CDF) and associated earthquake potential score (EPS) computation.
These time-dependent, time-independent, and exponentiated class of probability distributions (see Tables 1 and 2) have noteworthy applications in statistical seismology (e.g., Pasari and Dikshit 2014a, b, 2015a, b, 2018Pasari 2015Pasari , 2018Pasari , 2019a. For estimating model parameters, we use the method of maximum likelihood that involves maximizing a likelihood function for a given data. We also employ a surrogate Fisher Information Matrix (FIM) based approach coupled with Cramer-Rao lower bound theorem to compute parametric uncertainties (see, Pasari 2015; Pasari and Dikshit 2015a, b). After parameter estimation, the performance of the studied distributions is analyzed on the basis of two popular goodness-of-fit measures: Akaike information criterion (AIC) value and Kolmogorov-Smirnov (K-S) distance. While AIC is founded on information theory and is designed to prioritize a number of distributions based on the relative amount of discrepancy (lost information) from the true distribution, the K-S test is a non-parametric procedure that prioritizes candidate distributions based on the vertical distances between the empirical cumulative distribution function and the mathematical CDF of the assumed distribution (Johnson et al. 1995). Parameter estimation and model selection results are presented in Tables 1 and 2. Measures of K-S and AIC statistics suggest that exponential distribution has the best fit to the observed natural times in the Sulawesi region. This means that the underlying natural-time seismicity statistics in Sulawesi Island exhibits a stationary Poissonian process having the memoryless property.

Earthquake potential score: illustration for cities
Once the natural-time statistics of the larger region (A) is derived, we define earthquake potential score (EPS) for a circular city region (A i , i = 1, 2, · · · , n) as EPS ≡ P{X ≤ n(t)}; n(t) is the number of small earthquake counts at (calendar) time t since the last large event within the city circle. Notice that EPS is a monotonically non-decreasing function which retunes to 0 immediately after the completion of an earthquake cycle.
For the present analysis, we compute nowcast values or EPS scores for 21 cities in the Sulawesi region using the best-fit mathematical (exponential) CDF formula. These cities are Palu, Makassar, Manado, Kendari, Bitung, Gorontalo, Palopo, Bau-bau, Luwuk, Toli-toli, Pare-pare, Poso, Bone, Kolaka, Buton, Watampone, Polewali, Donggala, Taliabu, Mamuju, and Morowali. Centered at the city center, a circle of radius 300 km defines the city region (A i , i = 1, 2, · · · , 21) (Fig. 3). Assuming that the natural-time statistics of A and A i are the same, the EPS scores are nothing but the exponential CDF value evaluated for the current number (as on 21 November, 2020) of small event counts n(t) . A demonstration of the EPS score for the Palu city is provided in Fig. 3. The computed nowcast scores, as on Fig. 3 Histogram of natural time counts is shown along with a demonstration of EPS computation. The dark-green bars are the histograms of the number of small earthquakes (4.5 ≤ M < 6.5) between the large earthquakes (M ≥ 6.5). The purple step-function is the empirical distribution function derived from the histogram values. The dashed blue curve is the best-fit exponential distribution for the present natural time data. Just for demonstration, we show the black solid circle that corresponds to the number of small earthquakes as on November 21, 2020 since the last large earthquake in Palu region. The thermometer-like red bar represents the Earthquake Potential Score (EPS) showing the progression of the earthquake cycle in Palu as on November 21, 2020. However, as the seismicity in Sulawesi exhibits stationary Poisson process, the current level of seismic hazard (risk) is the same for the entire region  For the exponentiated exponential model, the parametric uncertainties cannot be computed as the trigamma function ψ ′ (β − 1) is not defined for β < 1 (see Pasari and Dikshit 2015, 2015a, 2015bPasari 2018) 2 For the computation of asymptotic standard deviation and 90% confidence interval of the model parameters, we followed a surrogate approach of Fisher Information Matrix and associated Cramer-Rao lower bound theorem (see Pasari and Dikshit 2015, 2015a, 2015bPasari 2018) Table 3 Nowcast values for M ≥ 6.5 earthquakes as on November 21, 2020 in the Sulawesi area with M σ = 4.5 and R = 300 km These EPS scores represent the current level of seismic progression of each city region in its cycle of large earthquakes. However, as mentioned before, due to the stationary Poisson process in Sulawesi, these EPS scores do not have any relevance to the associated seismic risk in the study region is about 82% of the way through its cycle. The estimated 90% CI of the EPS values is also mentioned in Table 3. However, as mentioned above, the information of different levels of cycle progression has no relevance to the current state of seismic hazards in the Sulawesi Island, as the underlying stationary Poisson process is memoryless. Therefore, the seismic risk is the same in time and space throughout the entire Sulawesi province.

Sensitivity testing
Although nowcasting technique inevitably provides a surrogate tool to examine the current progression of earthquake cycle of large events, numerical EPS values are often influenced by the selection of magnitude thresholds (M σ , M ) and geographical areas (A, A i ) . Therefore to perform a sensitivity testing, we vary small earthquake magnitude M σ from 4.0 to 5.0 and circular city radius (R) from 250 to 350 km. The associated results are pictorially depicted in Fig. 4. Although the change in small magnitude threshold M σ affects the EPS values by 5%-15% in the study region, the results are fairly stable against the changes of magnitude threshold and area of local regions (Pasari and Sharma 2020).

Discussion and summary
The present paper provides an application of a relatively new technique, known as earthquake nowcasting , to evaluate the current progression of earthquake cycles at different cities of Sulawesi region through the concepts of natural time. The nowcast scores may characterize the contemporary state of earthquake hazard in a geographical region. In this method, the key idea is to derive the seismicity statistics of a large defined area and then apply the same to the regional areas. The process requires constructing cumulative distribution function for natural times obtained from a sequence of earthquake cycles in the defined geographic region. In effect, natural times, counts between major events, act as a measure of stress-strain accumulation between large earthquakes in a specific region. There are at least three advantages of this approach. First, the method inherently accounts for the possible contribution from dependent events as natural time count is unaffected when aftershocks dominate, when background seismicity dominates, and when both contribute. Second, statistics of natural times remains unaltered with the changes of background seismicity rate as long as b-value remains constant. Third, the nowcast score-a measure of the current progression of earthquake cycle, often enables a systematic ranking of cities based on their current exposure to the seismic risk in a rapid, efficient and reproducible way. However, we re-emphasize that in the Sulawesi region, as the underlying natural-time statistics exhibits a stationary Poissonian process, the EPS values are totally irrelevant to the current level of earthquake hazards that is the same throughout the study region. As noted earlier, the method of earthquake nowcasting requires an assumption that the b-value is homogenous in space and time (e.g., Rundle et al. 2016Rundle et al. , 2018. To verify this basic assumption, we estimated b-values (using a least-squares approach) for 21 sub-catalogs corresponding to the circular city regions, as summarized in Table 4. The frequency-magnitude statistics of these city-circle sub-catalogs demonstrate that the b-value estimates are largely homogenous in space and time.
In addition to computational advances (e.g., Rundle et al. 2018), improvements of the nowcasting idea in seismology have taken place in two predominant directions-to strengthen the foundation of nowcasting in light of the ergodicity dynamics of statistical physics coupled with information theory (e.g., Holliday et al. 2016;Rundle et al. 2016Rundle et al. , 2018Rundle et al. , 2019a and to explore applicability of nowcasting for anticipating the occurrence of induced seismicity (e.g., Luginbuhl et al. 2018a, b), global seismicity (e.g., Luginbuhl et al. 2018c;Rundle et al. 2019a, b) and global tsunami risk (e.g., Rundle et al. 2019a, b). In order to support the presumptions in nowcasting analysis, Rundle et al. (2003Rundle et al. ( , 2012Rundle et al. ( , 2018 have used several concepts like percolation, ergodic property and phase transitions of statistical mechanics of complex dynamical threshold systems (Tiampo et al. 2003). Efforts are also made to generalize nowcasting into forecasting through a Weibull projection (e.g., Holliday et al. 2016) or based on the knowledge of the current state of the system and the assumptions of the probability distributions (e.g., Perez-Oregon et al. 2020). Recently, Rundle & Donnellan (2020), through some machine learning algorithms of the seismicity clustering patterns of aftershock sequences and seismic swarms, have further extended the application of nowcasting method for regional fault characterization in Southern California.
To summarize, the nowcast scores of 21 cities in Sulawesi provide a snapshot evaluation of the extent of cycle progression in a nice 0-100% scale of extremeness. However, the presently revealed stationary Poissonian nature of the underlying natural-time statistics in the study region indicates that the current state of earthquake hazard (risk) is the same for all city regions despite their different levels of cycle progression at current time. In addition, though earthquake potential scores of the city regions will be updated with the occurrence of each small earthquake in the respective city region, the seismic risk remains the same throughout the Sulawesi Island. Therefore, the present nowcasting analysis in Sulawesi fails to produce a space-time variant risk assessment as sought by engineers, city planners or policymakers. Fig. 4 Sensitivity testing of nowcast scores for different cities in Sulawesi region. The first subplot demonstrates the circular city region for the Palu city corresponding to a radius of 250 km, 300 km, and 350 km. Current number of small event counts since the last major earthquake in each circular city region are shown in the lowermost panel. These EPS values, however, bear no information of the current level of earthquake hazards due to the stationary Poissonian nature of the underlying natural-time statistics Data and resources Earthquake data were downloaded from public and regional catalogs: advanced National Seismic System comprehensive catalog (http:// www. ncedc. org/ anss/ catal og-search. html), Meteorological Geophysical Agency of Indonesia catalog (http:// repog empa. bmkg. go. id/ repo_ new/), and International Seismological Centre catalog (http:// www. isc. ac. uk/ iscbu lletin/ search/ catal ogue/). All websites were last accessed in November, 2020. The compiled dataset is provided as a additional file.
Authors'contributions SP formulated the problem, carried out preliminary analysis and drafted the manuscript; AS compiled earthquake data, prepared figures and participated in drafting; Neha and YS carried out statistical analysis and prepared some of the figures. All authors read and approved the final manuscript.

Funding
Currently no funding is available for this research.

Competing interests
No potential conflict of interest was reported by the authors.