Strong historical earthquakes and their relationships with the Tan-Lu fault system and modern seismicity in eastern China

The largest earthquake (Ms 8½) in eastern China occurred near Tancheng in 1668. Six M ≥ 6 earthquakes occurred in this region in historical time. We re-examine the primary historical archives from this region, adding independent constraints such as maps of active faults, local topography, intensity attenuation laws for China, and relevant administrative practices in archiving data, to reassess source parameters (epicenter, magnitude, or maximum intensity). We use the seismic density index—a measure of the clustering of seismicity weighted by source magnitude—to identify candidate zones of present-day clustering of seismicity than may represent a long-term memory of past historical events in an intraplate setting. We update at least one of the source parameters for four of the seven events. The main shock and two aftershocks of the 1668–1672 Ms 8½ Tancheng earthquake sequence locate along the strike of the Tan-Lu fault zone, each associated with significant and contiguous seismic density anomalies whose elliptical contours are elongated along its strike. A third major aftershock (Ms ≈ 6) in this sequence locates some 200 km away from the fault zone. Six of the seven historical events have seismic density anomalies that cannot be explained by aftershocks of modern seismicity (ML  ≤  5.3), and hence are candidate zones for long-lived aftershock sequences. In contrast, there are no significant seismic density anomalies along the section of Tan-Lu fault south of Suqian, despite palaeo-seismic evidence for a Ms 8½ earthquake around 6280 B.C., providing an upper bound to long-term memory of great earthquakes in this region.


Introduction
One of the outstanding problems in seismic hazard analysis is to construct a reliable seismic source model and the likelihood of recurrence as a function of frequency, magnitude and position on a map (Reiter 1991). It is common to start with a spatial seismic source model consisting of geometric zones within which seismicity is considered likely to be uniform, and/or to characterize the spatial variability with a smoothly-varying seismicity based on the local density of events. Modern instrumental earthquake catalogs can place bounds on this, but are commonly restricted to several decades. The catalog can be extended by consulting the historical record, revealing large events in the past that may or may not occur in areas where the modern seismicity is concentrated. On longer timescales still, the strain rate, the pattern of active faulting mapped at the surface, and evidence of past earthquakes from palaeoseismology, can provide independent constraints on the seismic potential to fill in gaps where modern seismicity may not be representative of the longterm recurrence rates (e.g., Field et al. 2014;Bayliss et al. 2020).
China has a particularly long historical record, with formal archives recording unusual geophysical, astronomical and meteorological events dating back to the twenty-third century B.C., and an increasing effort in palaeoseismic investigation of past earthquakes (Min et al. 1995) to place bounds on the location, magnitude and time of occurrence of past major earthquakes, defined in China as being M s ≥ 4¾. These archives are used to define 'the maximum intensity' area defined by a threshold intensity contour (usually taken to be 'the most heavily damaged'), and used after suitable calibration with an attenuation model to estimate the historical magnitude. The rupture length can be constrained by the longest dimension of the maximum intensity area, and its strike by its orientation (Musson 1998;Wang and Lin 2020), albeit with a significantly larger error than estimates of these parameters for modern strong earthquakes.
Individually, these data sources have their own limitations. The modern seismic record is limited by its short duration of up to one hundred years or so-a snapshot in time compared to the longer-term recurrence intervals for large events, especially in slowly-deforming intraplate areas. Historical data can extend this record to more representative time windows, but there is often not enough data to delineate a suitable maximum intensity zone and hence provide an accurate estimate of source magnitude or location (Musson 1998;Wang and Lin 2020). Active faults can preserve evidence of past large events in sediments laid down across the fault, but not all faults reach the surface, and recent sedimentary deposits may be reworked after deposition to wipe the record clean. Hence the palaeo-seismic record is also not complete, especially in intraplate regions (Hough and Hough 2000). While strain rate can be a very effective constraint in inter-plate areas, many earthquakes in intraplate areas occur on zones with low strain rate (e.g., the 1811-1812 New Madrid sequence in eastern the United States, Newman et al. 1999). Our research region is located in the North China Plain, where the strain rates are generally less than 3.5 ± 0.3 mm/yr. Sinistral motions across NNE trending Tan-Lu faults at a slip rate of less than ~ 2 mm/yr Wang and Shen 2020).
In some cases, there is a very close relationship between large historical earthquakes and clustering of modern seismicity, particularly in intraplate areas where the background seismicity rates are low to intermediate, allowing long-lived aftershock sequences to be detected above the background event rate more easily. There are many examples (Mueller et al. 2004;Stein and Liu 2009;Castro et al. 2010;Page and Hough 2014;Fereidoni and Atkinson 2014), including the classic case in New Madrid, eastern 1 3 US (Zoback et al. 1980) and also in China (Wang 1985). There is no uniform definition of a long-term "aftershock" sequence, but we take this to mean a cluster of modern seismicity associated with a historical event occurring after an elapsed time much longer than the few months or years of similar sequences in inter-plate areas, and which cannot be explained by short-term aftershock sequences in the modern record alone.
The existence of long-lived "aftershock" sequences implies their presence could be used to indicate the location of unrecorded paleo-seismic events or historical earthquakes (Ebel et al. 2000;Wang et al. 2017). Generally, aftershocks are strongly clustered within a finite area close to the source rupture, and hence can be used to estimate the rupture area and its orientation. This area changes only very slowly in time, consistent with anomalously slow stress diffusion (Huc and Main 2003), which allows the possibility of retaining a high degree of clustering and memory of the source rupture over long time scales. These properties imply long-term "aftershocks" can also be used to constrain the source rupture plane in intraplate areas, at least on a timescale of a few hundred years, and where seismicity clusters are elongated along strike (e.g., Zoback et al. 1980 for the New Madrid zone in the eastern US). In previous work, we confirmed the presence of modern-day seismicity clusters that may represent either long-lived "aftershock" sequences, or zones of persistent weakness in the lithosphere, near Beijing, China (Wang et al. 2017).
Here we test the hypothesis that modern-day clustering of seismicity contains long-lived "aftershocks" of historical events in eastern China. The largest earthquake in this region occurred near Tancheng in 1668, during the Qing Dynasty, with a magnitude of M s 8½ estimated from observations of damage recorded in the rich historical archive. In Chinese historical earthquake catalogs, the magnitude was expressed to the nearest quarter unit to indicate the error. The earthquake was associated with area of macrosesimic intensity VIII of 16,800 km 2 , with an elliptical shape whose long axis follows the strike of Tan-Lu fault-a large-scale fault cutting through the crust in eastern China. Some six other major earthquakes M s ≥ 6 have occurred in this region in historical time. The area also has a good instrumental record of moderate and small earthquakes and contains the Tan-Lu fault, a major active feature which crosses the region, making it an ideal place to examine the relationship between modern-day seismicity, strong historical earthquakes and large scale faulting in an intraplate setting.
First, we re-evaluate the historical record, consulting the original archives in the light of independent constraints including the map of active faults, a topographic map and clustering of modern-day seismicity, and reassign source parameters to four of the seven major historical earthquakes. We use a magnitude weighted metric of clustering known as the seismic density index to identify candidate zones containing long-lived "aftershock" sequences, and establish a close relationship between the density contours and the event location, fault location and strike for the great Tancheng earthquake mainshock and two of its major aftershocks. We find the clustering of modern seismicity for six of the seven events is not caused by short-lived aftershock sequences of intermediate magnitude seismicity (M L < 5.3) in the modern catalog, instead having a more stationary character consistent with the tail of a long-term "aftershock" sequence of the historical events or more persistent zones of weakness. Conversely, there are no significant seismic density anomalies along the section of Tan-Lu fault south of Suqian, despite palaeo-seismic evidence that this segment suffered a paleo-earthquake similar to the M8½ Tancheng earthquake around 6280 B.C. ). This places a conservative upper bound for the duration of long-lived aftershock sequences of great earthquakes of up to several thousands of years in the south section of the Tan-Lu fault.

3 2 Input data
Our research region is chosen as a quadrangle (116.5-120.5°E, 33.0-38.0°N, shown in Fig. 1) to cover the area of intensity VIII of the Tancheng M s 8½ earthquake. The map also shows the location of the seven major (M s ≥ 6) historical earthquakes, the locations of  Table 1 (black circles, dates and magnitudes as shown) and local magnitude for instrumentally-recorded earthquakes of M L ≥ 1 since 1970 (gray circles), and the distribution of mapped faults (thicker black lines for the Tan-Lu fault zone, thinner gray lines for the others). A Xian (county) is the lowest level of administration. A dot with a circle indicates the position of the capital town of a Xian. Nowadays, some Xian have been re-designated as a City 1 3 modern earthquakes, and the map of active faults, with the Tan-Lu fault zone highlighted in darker lines. In this paper, the magnitude of strong historical earthquakes is of surface wave magnitude (M s ); the magnitude of small and moderate earthquakes is of local earthquake magnitude (M L ).

Historical earthquake data
In China, there are four formal editions of the historical earthquake catalog. "The Catalog of Chinese Historical Strong Earthquakes (2300 B.C.-A.D. 1911)" is the latest edition, known as 'C4' (Min et al. 1995). This catalog was compiled using common conventions that include estimates of the uncertainties in the parameters. For example, magnitudes are expressed in multiples of 1/4, to emphasize these are more uncertain than modern magnitudes, which are commonly reported to one decimal place. Latitude and longitude of the epicenters are expressed to one decimal place for historical earthquakes, compared to at least two for modern epicenters, (Li et al. 1960(Li et al. , 1971Gu et al. 1983;Min et al. 1995;Wang and Lin 2020). The location accuracy has five categories: ± 10, 25, 50, 100 and over 100 km respectively (Min et al. 1995). The parameters and uncertainties were determined by the state of knowledge at the time the catalogs were compiled, and there are often no explanations for the reasons. The C4 catalog contains 7 strong historical earthquakes (M s ≥ 6) in the research region (listed in Table 1, along with their parameters and uncertainties, where recorded), including the Tancheng earthquake and two of its large aftershocks in 1668. The origin time of each earthquake in Table 1 is reported on the Gregorian calendar, after conversion from the Chinese lunar calendar reported in the primary archives.
Our research region is located in the core of the cradle of Chinese culture. The first historically recorded earthquake in this region occurred in 70 B.C. (indicated by a negative date in Table 1), with a large epicentral uncertainty of ± 100 km due to the greater scarcity of primary data. Generally, the descriptions of damage from early earthquakes are relative simple, and the number of locations with damage records is limited. For example, the 1668 A.D. Tancheng M s 8½ earthquake is the only event with sufficient recordings to estimate contours of macro-seismic intensity directly from the map of local intensity values (Fig. 1). In addition to local accounts, which are often not preserved in archives accessible today, it is also common for contemporary reports of the damage suffered by local towns and villages to be summarized in the central archives of a 'Xian' or a larger regional capital. A Xian was the lowest tier of local government, which governed towns and villages. If the earthquake was felt at a local village or town, the record for that location might consist of only two words "earth quake", with no any other description of the damage, unless to a famous temple, pagoda, or other notable building at that location (Wang and Lin 2020). Thus, historical seismologists must interpret simplified records where descriptions of damage are stamped primarily with the name and level of a Xian and/or regional capital, and estimate uncertainties accordingly. This loss of location data and associated bias compared to modern practice is a source of systematic error (Min et al. 1995;Wang and Lin 2020), and often makes it impossible to draw a reliable isoseismal map from such records without independent constraints. In such cases (for events Label D and E) the parameters of the earthquakes were estimated from the felt area only, so there were no estimates of the maximum intensity and its uncertainty. In this paper, we use the Chinese intensity scale, which has 12 categories. It was created in 1950s for application to China, informed by and referring to alternates such as the modified Mercalli intensity scale and the intensity scale of the former Soviet Union (Wang and Lin 2020).

Instrumental earthquake data
Seismic networks have been deployed since 1966 in eastern China. In our research region, a total of 5711 earthquakes with M L ≥ 1 are listed in the catalog between 1970 and 2019 (shown in Fig. 1). Since 1974, earthquakes with M L ≥ 2 are considered to be completely reported (Lin and Wang 1991). The largest earthquake in this period occurred in 1995 with a magnitude of M L 5.3. This maximum magnitude from the instrumental period is much smaller than the estimated magnitude of the major historical events listed in Table 1. The frequency-magnitude distribution (FMD) from the instrumental catalog for the period 1974-2019 is shown in Fig. 2. Utilizing the software of ZMAP, the magnitude of completeness (M c ) was estimated (Wiemer 2001). The maximum curvature technique (Wiemer and Wyss 2000) gives M c = 1.7, but this is known to produce underestimates when the roll-over at small magnitudes is gradual, as in Fig. 2 (Roberts et al. 2015). The maximum likelihood estimates (Aki 1965;Woessner and Wiemer 2005) of the FMD parameters a and b for M c = 1.9, are a = 4.96 ± 0.03, b = 0.804 ± 0.015, where uncertainties are quoted at one standard deviation. Only 3191 of the 5711 earthquakes have a reported depth, some 56% of total. The distribution of recorded focal depth is listed Table 2, more than 95% earthquakes with a reported depth have a depth in the range 1-25 km. This is not atypical for an intra-plate continental setting. The frequency distribution of epicenter location uncertainty is listed in Table 3, expressed as a percentage of the total of those entries which report uncertainties (3610 events). A large majority (90%) of estimated epicenter location uncertainties are less than 10 km, and predominantly less than 5 km. Some 2101 earthquakes not included in Table 3 had no recorded uncertainty, about 37% of the total of all events. The operational or scientific reasons for the omission are not recorded.

Map of active faults and topography
We use the map of active faults published by . The major feature is the Tan-Lu fault zone, first discovered by an associated aero-magnetic anomaly in 1957. It has played a key role in the tectonics of eastern Asian from the Mesozoic to the present day (Gilder et al. 1999;Zhang et al. 2003;Zhu et al. 2010;Liu et al. 2015). Its strike is NNE-SSW and its mapped length is over 3000 km. The fault can be divided into three segments: the northern segment in northeastern China, a middle segment in Bohai Bay and Shandong Province, and the southern segment in Jiangsu and Anhui provinces. Here we address its influence in the middle segment.
The historical earthquake recordings came from Xie et al. (1983Xie et al. ( , 1987. In many cases, the topography provides useful constraints on the location of damage in the historical archive, for example some felt reports are assigned to locations 'near the mountains'. Here we use the open access topography map at web:https:// map. bmcx. com/ as a reference for interpreting such observations.

Earthquake clustering and the seismic density index
The distribution of small and moderate modern earthquakes in Fig. 1 is very complicated, including some obvious seismic clusters where there are many events, events are closer together, and where there is more energy release, as indicated by the event magnitude. All of these are indicators of ongoing stress release in the area, despite the fact that no earthquakes of M L > 5.3 have been recorded since 1970. In this paper, we use a magnitudeweighted estimate of clustering known as the 'seismic density index' to map out the intensity and directional properties of the clustering of seismicity for these small and moderate earthquakes, and to identify anomalies that may be related to past historical earthquakes. 1 3

The seismic density index
The seismic density index was introduced by (Wang et al. 2017), and applied to identify evidence of long-term clustering of seismicity associated with historical events near Beijing, China, so we introduce this only briefly here. Firstly, the research region is gridded into a square array of nodes separated by a distance Δ in longitude and latitude. Only events above the magnitude threshold for complete reporting are used. In 2.1, the M c is estimated as 1.9, here M min = 2.0, slightly larger than the M c to allow for some uncertainty in this estimate. For the jth node at time t, the seismic density index is defined by: where n is the number of events in the t'th time interval; M i is the magnitude of the ith event; r ij is the distance between the ith event and the jth node; Δm is a normalization factor, defined by the dynamic range or the difference between the threshold magnitude for complete reporting and the maximum magnitude. The finite minimum distance R min avoids problems with the singularity in 1/(ln r ij ) at r ij = 0. The maximum R is determined by a trade-off between having sufficient sampling and the spatial resolution. From (1) the seismic density for a node j increases with respect to M i and n, and decreases inversely with respect to r ij . Once the calculation has been applied to every node in the spatial grid, contours of the density index are constructed by interpolation to form a map where anomalies representing zones of concentrated seismic energy release can be identified. The fixed parameters in formula (1), such as grid size Δ , search radius R and resolution R min are determined with respect to the accuracy of epicentre locations listed in Table 3. Here, the grid size Δ is taken as 0.05° of longitude and latitude, roughly equivalent to uncertainty level I (i.e., the lowest uncertainty) in Table 3 (± 5 km). R min is fixed by the ultimate resolution of epicenter location at about 2-3 km (Lin and Wang 1991), so R min = e (e ≈ 2.71828) is taken as a mathematically convenient number to represent this. The maximum radius R should be not too small, so there are sufficient numbers of events for a representative statistical sample, and not contain too much overlapping data to avoid oversmoothing the contours. Here we take R = 10 km. The dynamic range Δm = (M max − M min ) is 3.3 for M min = 2.0 and M max = 5.3. The resulting contours of seismic density and the epicenter distribution for the magnitude thresholds M min = 2.0 are drawn in Fig. 3 The seismic density maps in Fig. 3 reflect the clustering properties of earthquake weighted by magnitude, and are more representative of the concentration of ongoing seismic energy release than would be the case from the unweighted density of smaller instrumentally-recorded earthquakes alone (Fig. 1). Major features include an elongate series of anomalies stretched out along the location and strike of the Tan-Lu fault zone, more isotropic anomalies near Kenli and Laizhou to the north, near the Bo sea, and near Yazhou to the west. One anomaly near Cangshan is associated with the largest event in the modern catalog.

Preliminary analysis of the relationship between seismic density and strong historical earthquakes
If strong historical events are associated with a long-term anomaly reflected in the clustering of present day seismicity, then we can use the seismic density map as a (1) reference in re-evaluating the location of the epicenter, or to indicate the possible presence of larger events that are not yet recorded in the archive (Wang et al. 2017). We now analyze the relationship between the seismic density map of Fig. 3 and the strong historical earthquakes listed in Table 1, illustrated for comparison in Fig. 4. The contours in Fig. 4 start from a threshold of I j,t = 2 and increase in unit steps. It is immediately apparent that the elongate density anomalies on the Tan-Lu fault are associated with the largest event and two of its major aftershocks. It is also evident that the more isotropic anomalies near the high amplitude anomalies at Fengxian and Yanzhou are  Table 1 appear to be associated with modern density anomalies, for example, the first recorded earthquake on 1 June 70 B.C., M s ≥ 7. From Table 1 (indicated by its label '−70.06.01 M s ≥ 7' on the map). However, the uncertainty in epicenter location of this earthquake was ± 100 km, leaving open the possibility of a relationship. This is examined in reference to the updated information below.  Table 1. The values of contours start from I j,t = 2 and increase in unit steps

Re-evaluation of the parameters of strong historical earthquakes
The parameters and uncertainties of the C4 catalog were determined by methods and practices developed and used by experts over 20 years ago, using their experience and personal judgment, and often without documenting the reasoning. In interpreting the catalog, it is therefore necessary to be aware of this lack of documentation (Min et al. 1995). Here we address some issues and opportunities that have developed since 1995, notably: (a) Improved understanding of the historical context, notably the practices of local administration of data collection, documenting and zoning (Tan 1996a, b), (b) Better models for the attenuation of macroseismic intensity (Wang et al. 1998(Wang et al. , 2000, (c) Use of topographic maps to interpret the historical records for relative spatial information (Wang 2011;Wang and Lin 2020). For example, in some cases, there were the descriptions in original records such as "heavy damage near a mountain or a river", (d) The recognition that seismic density index (defined in Sect. 3.1) can retain a long memory of historical events in China (Wang et al. 2017) as well as elsewhere.
Unfortunately, no Intensity Data Points are reported in C4, because (a) the location of archiving in the Xian is not a point but a county (Xian) and (b) the published isoseismal maps do not retain the original intensity data -only the contours themselves (Min et al. 1995;Wang and Lin 2020). In this section, we apply these four post-C4 developments to reassess the locations of historical earthquakes in the study area, within their stated uncertainties. As a result, we have re-assigned epicenters and/or magnitudes to four historical events, as indicated in bold in Table 5, but could find no reason to do so for the other two listed in Table 1. Of those with re-assigned parameters, we start with the 1829 earthquake, because it has the greatest potential for re-examination of local historical data with a large and more detailed archive of contemporary accounts of historical earthquakes. The original records were written with Chinese (Xie et al. 1987;Min et al. 1995). Here we just check and translate relevant excerpts from the contents of these documents, which are the same as those used in compiling the National C4 catalog.

The 19 November 1829 M s 6¼ earthquake
This earthquake occurred in the Qing Dynasty, and has an epicenter in the C4 catalog at the mid-point between the regional capitals of Yidu and Linqu (blue circle in Fig. 5), with an estimated error in location of ± 25 km and an estimated magnitude of 6¼ (Table 1, Min et al. 1995). Noting the practice of assigning intensity data at local capitals, we tested the hypothesis that the C4 estimate was optimal. First, we scrutinized the original records for critical observations of damage that could be assigned with confidence to a particular location. For example, "The damage is relatively heavy in Yidu and Linqu, especially for some villages where houses are built close to the mountains. Yidu: total 28 persons died and 7047 houses collapsed or were destroyed in 27 villages. Linqu: total 89 persons died and 9544 houses collapsed or destroyed. Another key constraint is "Boshan: cracks occurred in the mountain. Zouping: some old houses were damaged" (Xie et al. 1987, pp. 811-818).
In Yidu, the archives report some 27 villages that suffered damage, as well as the total number of fatalities and damaged houses, including those in the 27 villages. While it is possible to estimate a macro-seismic intensity from the records, it is not possible to assign these to a particular village, because neither their names nor their precise locations are reported in the original records. Hence, the compilers of C4 assigned the maximum observed intensity values of VIII to the capital towns of Yidu and Linqu. Meanwhile, a lower intensity value of VI was assigned to Boshan and Zouping (Min et al. 1995, p. 339).
In the absence of other information, the epicenter was then estimated in C4 at the midpoint between Yidu and Linqu (blue circle in Fig. 5).
The primary records highlight the importance of topography, also illustrated in Fig. 5, as a reference point. The mountain ranges lie to the south-west of the C4 epicenter. According to "The Historical Atlas of China, Vol. 8-Qing Dynasty", there was a town name as "Wujing" located at the boundary of the mountain (Tan 1996b) whose reported damage would likely have been included in the records for the regional capitals. There is also a present-day seismic density anomaly just west of Wujing in Fig. 5, also near an active fault as shown on the map. If we assume the seismic density index represents a memory of past events, then it may be better to use the center of the modern day seismic density anomaly near Wujing as a candidate alternate epicenter. This location is more consistent with the Fig. 5 Topography around the epicenters and the model isoseismals of the 1829 earthquake using the updated best estimate for the location (red star) and magnitude of the event. The C4 epicenter is shown for reference (blue circle). The four sites metioned in original records: Yidu (118.476,36.685),Linqu (118.541,36.508),Boshan (117.836,36.448) and Zouping (117.727,36.886) available information from the re-examined historical data and the topography than the C4 epicenter, located in the mountains south-west of the C4 epicentre, but within its uncertainty of ± 25 km, as shown by the red star in Fig. 5.
There is a statistical relationship between magnitude (expressed in a bin of given range) and the isoseismal area for each intensity, with a significant scater due to a combination of source focal mechanism, stress drop, and site response. The estimated isoseismal area was not always accurate, owing to insufficient detailed, credible accounts of the effects of the earthquakes (Bakun and Hopper 2004), so the estimated magnitude is also uncertain. The relationship between the magnitude bin and the isoseismal area for each intensity contour can be used as a calibration to judge the veracity of the estimated intensity and reporting site location and the. In practice the isoseismal area has an irregular outline, but in the absence of other information, we assume a regular circle with an 'eqivalent radius' which encloses the same area, calibrated on events where the data is sufficient to define a reliable outline. The equivalent radius in then the value used in the ground motion attenuation law. Here, we estimate the equivalent radius from the available isoseismal area data in China including historical earthquakes, which may or may not also include information on intensity data points. For a given candidate epicentre, the macro-seismic intensity contours can then be estimated based on the available historical data and a calibrated (isotropic) attenuation law for a given magnitude bin (e.g., Wang et al. 1998). This unavoidably leads to concentric circular intensity contours on Fig. 5 rather than the more complex structures that might result with a greater density and spatial coverage of primary intensity data points. The parameters of the specific calibrated attenuation law used are listed in Table 4 for different candidate source magnitudes (after Wang and Lin 2020).
By trial and error, a magnitude bin of 6.5-6.9 is found to be more consistent with the historical data and the independent constraints mentioned above than the other magnitude bins shown in Fig. 5, notably the reported the maximum intensity of VIII assigned in C4 to Yidu and Linqu. The preferred model has an intensity VIII covering the eastern periphery of the mountain, and intensity VI assigned to Boshan and Zouping. Given the new information from the attenuation law, we can rule out a magnitude bin of 6.0-6.4, which would not be consistent with the VIII intensity observed on the plain, where 27 villages suffered heavily damage over a wide area near the mountains. Similarly, a magnitude bin of 7.0-7.4 would be too high to explain the damage observed at Boshan and Zouping. Accordingly, we re-assigned the epicenter of 1829 earthquake to the center of the seismic density anomaly near Wujing, estimated its magnitude as 6¾, and reduced the estimated intensity at the regional capitals Yidu and Linqu to VII (Fig. 5). The re-assigned magnitude and maximum recorded intensity are both one category higher than that assigned in C4 (which had M s 6¼, I 0 = VIII, Table 1).
The parameter uncertainties in C4 were determined by a previous generation of experts and there are often no explanations for the reasons. However, our understanding is that the distances between the individual intensity data points, or equivalently the density of locations recording ground shaking, is the main factor in controlling the uncertainty. In this case, the most important intensity data points include Linqu and Yidu, a distance of about 20 km apart. The half distance is then 10 km. If there were only the two points, the uncertainty would have been given as ± 10 km in C4. However, the intensities from other intensity data points, such as at Boshan and Zouping, require a higher epicenter uncertainty of ± 25 km for this event.

The 1 June 70 B.C. M s ≥7 earthquake
This earthquake occurred in 70 B.C. in the Han Dynasty, about 2090 years before present. In C4, the epicenter was put at the mid-point between the capital towns of Yingling and Dongwu (blue circle in Fig. 6), with a reported uncertainty of ± 100 km. The maximum recorded intensity was estimated to be VIII implying an epicentral intensity I 0 ≥ IX,and a corresponding magnitude estimated as the constraint M s ≥ 7 in C4 (Table 1, Min et al. 1995). The inequality here reflects the fact that the best estimate of epicentral intensity may be higher than the maximum recorded intensity, depending on the attenuation law, and the location and the distance to the epicenter of the reporting point. We now follow the same procedure as in Sect. 4.1 to re-evaluate these estimates.
One key excerpt from the original records translates as: "There were 49 shires affected by an earthquake. In the Beihai shire and Langya shire, ancestral temples, town walls and houses were destroyed. Landslides occurred in the mountains and water gushed out. More than 6000 persons were killed." (Xie et al. 1983, p. 8). This account is typical of the historical earthquake records of the time in China, with damage descriptions that are concise and relative simple. It implies the felt region is large (49 shires affected), but that only two shires suffered heavy damage (Beihai shire and Langya shire). Shire was a higher tier local government than Xian. Unfortunately, there was no name of a specific town or village in the shires in the original record for this event. According to the conventional practice at the time, the damage would have been assigned to the capital town of a shire. Hence, the C4 epicenter was put at the mid-point between the two capital towns, which could introduce a large error in the estimated epicenter location (Min et al. 1995;Wang and Lin 2020), represented by the large quoted uncertainty above. The distance between the capital towns is about 69 km, which is too far to guarantee both sites suffered an intensity as large as IX, because the equivalent radius for intensity IX in the magnitude 7.0-7.4 range (see below) is only about 16.7 km (Table 4, Wang and Lin 2020). This rules out the possibility of an epicenter at the mid-point.
According to "The Historical Atlas of China, Vol. 2-Qin Dynasty and Han Dynasty" (Tan 1996a), the administrative areas of the Beihai shire and Langya shire occupied the areas shown along with other shires of the time in Fig. 6. The figure also shows the locations of Han dynasty towns (purple triangles) and modern cities and towns (filled circles) for reference. The capital town of Beihai shire was named "Yingling" and the capital town of Langya shire was named "Dongwu." From this map we would expect the damage associated with an event of this size at the C4 location to be significant at the capital town of Gaomi, in the shire of the same name, to the east of the C4 epicentre, similar to that reported for the shires of Beihai and Langya assigned to Yingling and Dongwu, respectively. The original record did not mention Gaomi, despite it being a regional capital, implying the source is more likely to be further away and towards the west of Gaomi town than the C4 epicenter.
There is a seismic density anomaly to the west of the epicenter of C4 in Fig. 6, so we simulate the intensity pattern that would have arisen if the center of the anomaly is taken as a candidate epicenter instead of the C4 estimation. After the trial calculation, we find that magnitude of 7.0-7.4 is more reasonable (Table 4), and apply modern estimates of the intensity attenuation law to reconstruct the intensity contours. From the results of the simulation (Fig. 6) Wucheng town in the shire of Langya is located within the area of intensity X consistent with the estimation of C4 that epicentral intensity ≥ IX (Table 1, Min et al. 1995). The town of Zhuxu is just outside the circle of X intensity, also in Langya shire. There are other four towns, in the region of VIII intensity in this shire in Fig. 6. All of these reports of damage would also have been assigned to Dongwu, the capital town of Langya Shire. The simulation is consistent with the original record showing heavy damage in Langya Shire, but the simulated intensity in Dongwu is only VII.
Apart from Yingling, the only known Han dynasty town in Beihai shire was Anqiu, the closest of all to the C4 epicentre on Fig. 6, due east of the candidate epicenter. Its damage reports would have been assigned to Yingling, so it is not possible to constrain the candidate epicenter from the east on the map. Both Yingling and Anqui are located close to the line of IX intensity on the simulation, representing significant damage. The simulated IX intensity region covers a more mountainous region to the south of the candidate epicentre. The simulated VIII intensity circle just touches the western boundary of the shire of Gaomi, almost all of which is located within the intensity VII contour, and thus would not have experienced heavy damage. On the balance of evidence, we conclude that the candidate epicenter near Wucheng is more likely than the C4 estimate.
In this case, the only two intensity data points are available at Yingling and Dongwu, some 69 km apart. Accordingly, the uncertainty estimated in the ± 100 category in C4 seems too conservative. Hence we estimate the uncertainty estimated in the ± 50 category.

The 26 July 1668 M s 6¾ earthquake
This earthquake occurred as an aftershock, just one day after the great Tancheng M s 8½ earthquake. In the C4 catalog, the epicenter was put at Anqiu (blue circle in Fig. 7), with a reported error in location of ± 50 km. The magnitude was estimated at M s 6¾.The earthquake was felt in a large area: the distance from the southern point to the northern point is about 900 km, but distance from east to west is difficult to estimate, because the eastern part is in the sea. All descriptions of the level of destruction were assigned to Anqiu in the original records. The description of the destruction quoted in (Min et al. 1995, p. 221) is "Anqiu: the earthquakes occurred day after day, and houses collapsed intermittently. Shandong: total 25 Xian affected, including some in the damage area of the main shock; Jiangsu: 8 Xian affected; Anhui: 2 Xian affected; Zhejiang: 2 Xian affected; Hebei: at least 3 places affected." As a consequence, it is not possible to estimate the epicentral intensity I 0 . In contrast, the estimate of the felt area quoted above is well constrained by the list of towns affected. Anqiu is not located in the center of the felt region, so it is very unlikely to be the epicentre. The southernmost recording of felt effects was at Shaoxing about 600 km from Anqiu, and the northernmost at Tianjing, some 300 km from Anqiu. Moreover, no significant damage was recorded at the regional capitals of Weixian and Gaomi (Fig. 7), as might be expected from the C4 epicentre. Together, these observations imply the epicenter is more likely to be to the south of Anqiu.
In the south of Anqiu, there is a seismic density anomaly centered near the historical town Guanshuai, so we tested a candidate epicenter at that location. After the trial calculation, we find that magnitude of 6.5-6.9 is more reasonable (in this case the same as the C4 magnitude), as illustrated in Fig. 7. Almost all of the area of intensity VII and some parts of the area for intensity VI for this event occur within the area of intensity X of the M s 8½ main shock, whereas the regional capital of Anqiu was outside this area. Peak intensity X for the main shock implies all houses would have collapsed completely in this zone in the mainshock, which could explain the absence of any further destruction recorded for the later event there, except at Anqiu, which recorded intensity VI during the later event. The original epicentral location is biased towards the northern part of the felt region. While this is a possible solution, we conclude the revised location of 35.90°N and 119.06°E illustrated in Fig. 7 is a more likely epicenter than the C4 one also shown, but there is no need to revise the magnitude in this case. Unfortunately, there are no observations from outside Anqui to confirm this inference.

The 17 August 462 M s 6½ earthquake
This earthquake occurred in the Nanbeichao Dynasty. Its epicenter was placed near Weishan with a reported error in location of ± 100 km, and its magnitude was estimated at M s 6½.The epicentral intensity was estimated as VIII (Table 1, Min et al. 1995).
Excerpts from the original records translate as: "earthquake, noise came from north. In the Lu shire, the hills tottered and the ground shook. In Pengcheng (nowadays named as Xuzhou), the parapet of the town wall fell down and the houses toppled or collapsed. In Yanzhou, cracks appeared on the ground and water gushed. The shaking lasted on and off for more than two years". There were only two places mentioned in the records. One is Pengcheng (117.20°E, 34.26°N) and another is Yanzhou (116.83°E, 35.54°N). According to the conventional method in China, the epicenter would have been assigned at the midpoint between these two capital towns. Independently, the Atlas of the Historical Earthquakes in China (referred to as 'Atlas' below) provides another valuable catalog for reference (Institute of Geophysics in State Seismological Bureau and Institute of Chinese Historical Geography in Fudan University 1990). The original records quoted in C4 and in the Atlas are exactly the same for this earthquake, but the epicenter for this earthquake given in the Atlas is around 10 km further eastern than the C4 estimate (Table 5). In this case, we prefer the Atlas epicenter, because it is closer to the midpoint, but could see no reason to revise the maximum recorded or epicentral intensity, or the magnitude. In this case, Yanzhou and Xuzhou are the key intensity data points, separated by a distance of about 148 km. Given half of this distance is larger than 50 km, an estimated uncertainty in the ± 100 category seems reasonable.
To conclude this section, Table 5 summarizes where the catalog parameters have been updated (lower case label, e.g., 'a') or not (capital label, e.g., 'B' only) based on Table 5 Parameters of strong historical earthquakes (M s ≥ 6 in the C4 catalog) near the Tan  Southwest of Linqu Amended considering the new information presented above. The candidate zones identified by the seismic density anomalies in Fig. 8 are labeled by number (e.g., 'Z1') ordered with the oldest event first. The definition of the zones, based on the outlines of seismic density anomalies, is given in the next section.  Table 5. The dotted bracket indicates the total length of the contiguous slip segments of M s 8.5 earthquake in 1668 and that of the great event in 6280 ± 120 B.C 1 3

The relationship between modern seismicity, strong historical earthquakes and the Tan-Lu fault system
In this section, we compare the updated parameters of Table 5 with the map of seismic density anomalies from modern seismicity to examine the relationship between the two more broadly (Fig. 8). A seismic density anomaly is defined where there is at least one distinct peak, with a maximum seismic density ≥ 5, consist with the definition previously used in analyzing data from northern China (Wang et al. 2017). We delineate the boundary of a seismic density anomaly using the seismic density contour with a uniform value 3, traced in bold in Fig. 8 where the anomaly is associated with a historical event. Where two anomalies overlap, we distinguish between the two by manually tracing through the valleys on either side of the saddle point of the contours. The boundaries of the density anomalies defined in this way are traced in bold in Fig. 8 for zones Z1-Z7, labeled as in Table 5. The map also shows the off-fault aftershock on 18 September 1668 of the great Tangsheng earthquake for reference. The seismic density anomaly near this aftershock is more marginal, because its peak seismic density is less than 5, indicated by a label with a lowercase z and labeled z8. We conclude that all of the known strong historical earthquakes in this area are associated with a matching seismic density anomaly within these criteria and the uncertainties stated in Table 5. There are also many anomalies that are not associated with a historical earthquake, notably those near Fengxian and Yanzhou around Zone 6, already discussed in Sect. 3.2.

Spatial distribution of modern seismic density anomalies
The main feature of Fig. 8 is a set of three contiguous high-amplitude seismic density anomalies aligned with, and stretched out along, the strike of the Tan-Lu fault system, separated by two main saddle points. This contiguous zone spans from north of Tancheng in the south to Guanshuai in the north, and each anomaly is associated with a historical event.
A fourth anomaly centred on Xinyi further to the south aligns along the same trend within the fault zone, but is not associated with any historical earthquake. The epicenter of the M s 8½ main shock in 1668 is located in the southernmost edge of the contiguous anomaly in Zone 4, to the south of Linshu and the north of Tancheng. Its major aftershock of M s 6¾ occurred one day after the mainshock, and is located in the northernmost zone (Z2). In between there is a zone (Z3) associated with the epicenter of the 1672 M s 6 event near Juxian, filling in a 'gap' crated by the previous two events. Hence, this earthquake is also considered to be a major aftershock of the M8½ Tancheng earthquake in C4, despite the time lag of 4 years (Min et al. 1995). Its epicenter is associated with a more pronounced and higher-amplitude anomaly with a high peak seismic density of 24. These three zones of concentrated present day seismic energy release may therefore represent the low amplitude tail of a long-term aftershock sequence associated with these three events.
Another possible interpretation of the contiguous density anomalies Z4 is segmented rupture from the main shock. For example, with modern instruments it is possible to detect two sub-events of the M s 8 Wenchuan earthquake in 2008 (Zhang et al. 2008). Similarly, the M s 8 historical earthquake that occurred in 1679 near Beijing, has two modern-day seismic density anomalies with centers within the area of maximum recorded damage. One of these contained the macro-seismic epicenter on the ground surface and another was associated with modern day hypocenter where rupture started in the middle layer of the local crust (Zhang et al. 2008;Wang et al. 2017). A magnitude M8½ main shock would be expected to have a rupture length greater than 100 km, consistent with the combined length of zones Z2-Z4 (just over 1 degree of latitude on the map). Therefore, we cannot rule out the possibility that the whole contiguous anomaly was caused at least in part by a complex mainshock and its strong aftershocks, which might not be recorded completely.
The ≈M s 6 earthquake of 18 September 1668 with an epicentre near Tai'an in Table 5 is also regarded as an aftershock of M s 8½ Tancheng earthquake in the C4 catalogue, also shown in Fig. 8 for reference. Its epicentre is about two hundred km to the west of the Tan-Lu system, measured perpendicular to its strike. It is not associated with any major mapped fault. In C4, the magnitude of this aftershock is estimated as M s 6 from the felt region area only, indicating this estimate may have a large error. This aftershock is associated with a marginal seismic density anomaly z8 in the modern seismicity, which implies that the magnitude approximately is equal to M s 6.
The 70B.C. M s 7¼ earthquake has an estimated epicenter near Honghe, and occurs on one of the mapped faults on the Tan-Lu system. It is associated with a low-amplitude density anomaly (Z1), but there is no obvious elongation of contours along strike, as there is in Z2-Z4. If there is any trend, it is alignment with two smaller-scale faults on the map with NW strike, one of which ends near Linqu. It would be remarkable if these represented long-lived aftershocks of an event that occurred more than two millennia ago. This case illustrates the possibility that some historical earthquakes might be more likely caused by other small scale faults than the mapped Tan-Lu faults, albeit located on the Tan-Lu fault zone.
The M s 6¾ 1829 earthquake is also located in a low-amplitude local anomaly (Z7) off the main Tan-Lu fault system, near Wujing. Two small-scale faults aligned roughly parallel to the Tan-Lu strike are also located near Wujing, but again there is no elongation of the density contours to be confident of a relationship with the mapped active faults.
The M s 6½ earthquake of 462 AD is associated with a seismic density anomaly (Z6) at Weishan, and a nearby fault striking NNW to the SW edge of the anomaly. Again, there is no elongation of density contours to act as independent evidence of a possible relationship. Confirming this would require more information on the fault, such as its dip angle and the depth of the fault basement, and a higher resolution of epicenter location.
There is palaeoseismic evidence that the segment of Tan-Lu fault with evidence of slip during the M s 8½ Tancheng earthquake extends from Zhucheng to Suqian, a rupture length of about 220 km. The associated seismic density anomaly has a similar length to this inferred slip segment. There is also significant palaeoseismic evidence that another great event occurred around 6280 ± 120 B.C., which broke the Tan-Lu fault system from Linyi to Wuhe, to the south of Suqian . However, to the south of this segment, there is no significant modern seismicity or associated seismic density anomaly in this region in Fig. 8. We use this as a constraint to estimate the length of the slip segment for the Suqian great event (as named by ) is 210 km, very similar to that of the M s 8½ Tancheng earthquake, with a magnitude estimated in the range M s 8 to M s 8½.

Temporal evolution of seismicity in seismic density anomalies
In this section, we examine the possibility that the seismic density anomalies outlined in bold on Fig. 8 could arise from aftershocks of moderate-sized earthquakes (M s ≤ 5.3) in the modern catalogue. This includes the seven major earthquakes listed in Tables 1 and 5. This check is necessary, because for example Zone 5 is also associated with the largest earthquake since 1974 in the region shown, with magnitude M s 5.3 in 1995.
To look for evidence of temporal clustering due to aftershocks of modern events, we plot the annual number of earthquakes in all magnitude bins as a function of calendar year in Fig. 9, and list the relevant statistical parameters in Table 6. The annual rate of earthquakes M L ≥ 2 in each zone per 1000 km is between 0.77 and 1.65. Apart from Zone 5, the event rate in each zone is generally relatively stationary within Poisson counting errors, shown at 95% confidence on Fig. 9. This relatively stationary modernday event rate is typical of sequences associated with the long tail of aftershocks from strong historical earthquakes, for example in the eastern US (Page and Hough 2014) or near Beijing, China (Wang et al. 2017), representing the slowly-varying 'fat tail' of the Omori law for aftershocks. The event rate in Zone 5 is an exception. It has a clear outlier of 15 events in 1995, the year when the M L 5.3 earthquake occurred, compared to an average of 1,035 events per year per 1000 km 2 . This excursion is well outside the Poisson counting error, caused by aftershocks of the M L 5.3 earthquake. The density anomaly in this case is almost certainly due to clustering of energy release caused by this modern earthquake. However, this analysis shows that the effect of modern aftershock sequences (even this largest one) produces an anomaly that is much smaller than the amplitude of the other anomalies, lending more confidence to their association with the historical events.
In the marginal zone z8, there are only 5 earthquakes in the modern period, all of which occurred in 1985, the largest with M L 3.4. The local density anomaly is also very weak (only just above the density 3 contour). The strong aftershock M s 6 on 18 September 1668 near Tai'an is associated with this marginal zone. We conclude that the magnitude may be slightly less than 6, denoted M s ≤ 6 quoted in Table 6.

Conclusion
(1) By re-examining the original historical record and practices, and informed by maps of active faults, topography, clustering of modern seismicity, and understanding of macro-seismic attenuation laws, we have reassigned the epicenter and in some cases the associated magnitude and epicentral intensity for four out of seven strong historical events in the area with M s ≥ 6 in the C4 catalog. These more accurate parameters are a prerequisite for analyzing the relationship between strong historical earthquakes and modern seismicity. (2) There is a strong association of magnitude-weighted clustering of modern seismicity, measured by a seismic density metric, for all of these seven strong historical earthquakes. In one case the clustering can be ascribed to aftershocks of the largest event in the instrumental catalog (M L 5.3), but the others have more stationary event rates that are consistent with their being long-lived aftershock sequences of the historical events, or indicating otherwise persistent zones of lithospheric weakness. (3) A large number of anomalies in seismic density are not associated with any know major event in the historical, instrumental or palaeo-seismic record, any known mapped fault, or obviously with a short-term aftershock sequence of a modern event. This might imply that there are some strong events missing in the historical earthquake catalog from 'blind' faults, or swarm-like clustering due to another cause. (4) There is a particularly strong association of magnitude-weighted clustering of modern seismicity for the 26 July 1668 M s 8½ main shock on the Tan-Lu fault and its two major aftershocks. The first occurred on the main fault trend one day after the mainshock, and the second filled in a gap between the mainshock and this first aftershock in 1672, some four years later. These three events are all associated with contiguous local anomalies in seismic density that have contours stretched out along the fault strike, separated by local saddle points in the clustering metric. (5) A magnitude M s ≤ 6 event occurred some two months after the 26 July 1668 M s 8½ main shock, some 200 km away, and perpendicular to the main fault. Thus, it is possible that some strong aftershocks might occur far away from a main shock epicentre, and off the trend of the main fault. (6) There are no significant seismic density anomalies along the section of Tan-Lu fault south of Suqian, despite geological evidence for a M s 8-8½ palaeo-earthquake occurring around 6280 ± 120 B.C. Together with conclusion (4) this evidence of absence suggests long-lived aftershock sequences can have a duration of up to several thousands of years for the largest events. Our results are consistent with a recurring sequence of intraplate earthquakes in New Madrid Seismic Zone, with recurrence intervals of 250 to 4000 years. These intraplate earthquakes could be triggered by a perturbation deep in Earth's crust, when accompanied by viscous flow in Earth's lower crust (Kenner and Segall 2000), or possibly even deeper in the upper asthenosphere.