Research ArticleGEOLOGY

A 220,000-year-long continuous large earthquake record on a slow-slipping plate boundary

See allHide authors and affiliations

Science Advances  27 Nov 2020:
Vol. 6, no. 48, eaba4170
DOI: 10.1126/sciadv.aba4170


Large earthquakes (magnitude ≥ 7.0) are rare, especially along slow-slipping plate boundaries. Lack of large earthquakes in the instrumental record enlarges uncertainty of the recurrence time; the recurrence of large earthquakes is generally determined by extrapolation according to a magnitude-frequency relation. We enhance the seismological catalog of the Dead Sea Fault Zone by including a 220,000-year-long continuous large earthquake record based on seismites from the Dead Sea center. We constrain seismic shaking intensities via computational fluid dynamics modeling and invert them for earthquake magnitude. Our analysis shows that the recurrence time of large earthquakes follows a power-law distribution, with a mean of 1400 ± 160 years. This mean recurrence is notable shorter than the previous estimate of 11,000 years for the past 40,000 years. Our unique record confirms a clustered earthquake recurrence pattern and a group-fault temporal clustering model, and reveals an unexpectedly high seismicity rate on a slow-slipping plate boundary.


The understanding of earthquakes, in general, and seismic hazard, in particular, relies on our knowledge of past seismic history. The longer a time window that a record spans, the better that understanding can be. However, large earthquakes [moment magnitude (Mw) ≥ 7.0] usually have recurrence intervals longer than the time span of modern seismograph operation of about a century and infrequently occur on individual faults (1). Paleoseismological trenching at suitable sites can extend the record for surface rupturing earthquakes to the past few thousand years. Subaqueous paleoseismology exploits lacustrine and marine sediments to retrieve much longer records of paleoseismic shaking. These long records can probe our understanding of the physical behavior of fault systems in the wake of large seismic events and are therefore essential for improving seismic hazard assessment.

The earthquake recurrence pattern is a key issue regarding the understanding of fault behavior and seismic hazard assessment. Regular recurrence patterns exist in records of thousands of years on geometrically simple and fast-slipping strike-slip faults, such as the southern onshore section of the Alpine Fault in New Zealand (2) and Wrightwood Section of the San Andreas Fault (3). Earthquake recurrence patterns of slow-slipping faults (<5 mm year−1), e.g., the Dead Sea Fault (4, 5), are more difficult to determine because they usually have longer interseismic intervals. Long and continuous paleoseismic records of several 104 years on these slow-slipping faults will improve statistical analyses and hazard assessments because the length of the record can compensate for the uncertainty in event identification and dating (3). In addition, these paleoseismic records with extremely long time spans will be more representative of the complete fault history (6).

The Dead Sea Fault is a left lateral, >1000-km-long strike-slip plate boundary separating the African and Arabian plates (Fig. 1A). Short-term slip-rate estimates from GPS measurements of 4.2 to 5.8 mm year−1 along the central to the southern part of the fault (7, 8) are similar to the long-term geologic rate (9). The Dead Sea Basin is the deepest and largest continental tectonic structure along this plate boundary and has a width of 15 to 20 km and a length of ~150 km (Fig. 1). During the Quaternary, a sequence of terminal water bodies occupied the basin. The location of the Dead Sea Basin makes the soft and water-saturated sediments that accumulate in the basin excellent recorders of seismic shaking. Seilacher (10) and El-Isa and Mustafa (11) first hypothesized that the asymmetric folds of unlithified sediments exposed at the eastern Dead Sea margin are earthquake-triggered deformations. Along the western Dead Sea margin, layers of shattered folds in the form of intraclast breccia layers are common (12). The juxtaposition of these layers against syn-depositional faults (13) and the time correspondence of historical and archaeological earthquakes during the past 3 thousand years (ka) (1416) indicate that these layers are subaqueous seismites. These pioneering studies set the stage for our current subaqueous paleoseismic research.

Fig. 1 Tectonic setting of the Dead Sea Fault.

(A) Dead Sea Fault is a sinistral boundary between the African and Arabian plates (43). (B) Major active faults (43, 44) along the plate boundary, Dead Sea Transform; in this area, the fault is composed of four fault segments. The red star marks the drilling site; the black points mark places referred to in the study; the magenta triangles indicate historic and instrumental Mw ≥ 6.0 earthquakes near the drilling site (45). (C) The gray bars represent the fault rupture of historic Mw ≥ 7.0 earthquakes since 31 BCE (34, 36) that occurred along the focused part of the fault—the central Dead Sea Fault (up to 150 km north and south of the drilling site).

Heifetz et al. (17) proposed that earthquake-forced shear, leading to sediment turbulence called the Kelvin-Helmholtz instability, is a plausible driver for folding and shattering in the soft sediments in the Dead Sea Basin. The laminated Dead Sea sediments consist of stably stratified water-saturated mud in which the density increases with depth, and thus are not susceptible to Rayleigh-Taylor instability, which requires inverted densities (17). Further analysis shows that the geometry of folded layer textures in the Dead Sea sediments obeys a power-law of −1.9, similar to Kelvin-Helmholtz turbulence in other environments (18). The shear seismic energy provided by earthquakes may exceed the gravitational potential energy, thus allowing the laminated and stably stratified layers to move horizontally in the same direction but with different velocities, creating shear that localizes at the layer interfaces (17, 18). Previous numerical simulations (17, 18) confirm the “intraclast breccia layer” as the final stage of soft-sediment deformation. Moreover, these simulations confirm that the observed textures are a proxy for shaking intensity by associating the stage of the deformation (in different thicknesses) with minimum ground accelerations.

The interpretation of deformed sediment layers has been the primary means for recovering paleoseismic records of different magnitudes and time spans along the central Dead Sea Fault for up to 60 ka ago (4, 16, 19). However, previous paleoseismic records along the central Dead Sea Fault (4, 5, 1921) are either short or incomplete in their record of moderate earthquakes (5.0 < Mw < 7.0), and the constraints of local intensities and magnitudes are poor. Here, we investigate earthquakes recorded in the sedimentary sequence of the deep ICDP (International Continental Scientific Drilling Program) Core 5017-1 from the Dead Sea depocenter (Fig. 1B; Materials and Methods). The previous dating constrains the age of the 457-m-long composite core spanning from ~220 ka ago to the present (fig. S1 and table S1) (22, 23). We use this record to (i) identify and measure all folded and brecciated layers, (ii) constrain the shaking intensities of individual events via computational fluid dynamics modeling, and (iii) assess the recurrence pattern of large earthquakes and fault behavior model during the past 220 ka.


Earthquake indicators and paleoevents

As an ultimate repository for mass wasting (24) with an average sedimentation rate of 2 mm year−1 (25), the Dead Sea depocenter provides the most complete record of earthquake shaking along the plate boundary. Alternating laminae of white aragonite and dark detritus that characterize the sedimentary sequence of the ICDP Core 5017-1 serve as sensitive markers for identifying earthquake-induced deformation (24, 26). We identify in situ folded layers and intraclast breccia layer as earthquake indicators (seismites), based on their resemblance to outcrop observations of seismites that are known to be earthquake induced (4, 1216).

In situ folded layers. Similar to the structures preserved in the Dead Sea margin outcrops, the folded layers in the drilling core from the Dead Sea depocenter also appear as folded aragonite-detritus laminae in various forms of (i) linear waves (Fig. 2, A and B), (ii) asymmetric billows (Fig. 2, C to F), and (iii) coherent vortices (Fig. 2, G to J) (18). These delicate aragonite laminae are well preserved and can be traced in the strata, indicating that the layers are deformed in situ and have not undergone any notable transportation. That is, notable transportation would disaggregate and destroy the delicate submillimeter-thick aragonite laminae. Layer-parallel displacements characterize these in situ folded layers. The ICDP Core 5017-1 is positioned in the center of the Dead Sea abyssal plain. This makes improbable postdepositional causes for layer-parallel shears such as sloping substrates or downhill water flow above the sediments. In total, we identify 367 in situ folded layers in the ICDP Core 5017-1 (table S2). Figure S2 shows more examples of in situ folded layers in the drilling.

Fig. 2 Paleoearthquake indicators in the ICDP Core 5017-1.

(A to J) In situ folded layers; (A and B) linear waves, (C to F) asymmetric billows, and (G to J) coherent vortices. (K to M) Intraclast breccia layers. The vertical light blue bars indicate the position of events. Core depth: (A) 11,010.0 to 11,012.0 cm; (B) 16,604.0 to 16,608.0 cm; (C) 10,929.9 to 10,932.4 cm; (D) 26,582.7 to 26,585.2 cm; (E) 32,861.0 to 32,862.5 cm; (F) 35,921.8 to 35,923.8 cm; (G) 13,754.4 to 13,758.0 cm; (H) 10,605.4 to 10,606.9 cm; (I) 36,425.9 to 36,427.9 cm; (J) 12,528.0 to 12,532.0 cm; (K) 14,492.5 to 14,500.0 cm; (L) 39,206.4 to 39,210.4 cm; and (M) 10,772.0 to 10,787.0 cm.

Intraclast breccia layer. Similar to the structures preserved in the Dead Sea margin outcrops, this type of layer from the Dead Sea depocenter consists of mixed aragonite-detritus laminae fragments (Fig. 2, K to M). The in situ deformation process of this type of layer is recognized by (i) the lack of erosion processes at the base of the layer and (ii) some remaining parts of in situ folded layer can be observed directly in the base of the intraclast breccia layer (Fig. 2, K to M, and fig. S3). Together, these features provide evidence for the evolution of the intraclast breccia texture under seismic shaking and differentiate it from any other detrital layers with laminae fragments that formed by secondary sedimentary processes such as mass transport. In total, we identify 46 intraclast breccia layers in the ICDP Core 5017-1 (table S2). Figure S3 shows more examples of intraclast breccia layers in the core.

Our observations reveal that some in situ deformed layers have been re-deformed by a latter deformation within a short core section. Using geophysical and chemical datasets and careful sedimentary structure analysis, we recognize these redeformed seismites. The potential uncertainties in thickness measurements (for redeformed seismites) that are induced by the redeformation range from a few millimeters to several centimeters and therefore have no notable effects on shaking intensity estimation. In addition, regarding the latter in situ deformations, the shape and thickness of total folded sediments constrain the intensity of seismic shaking regardless of the type and thickness of the redeformed seismites. In total, we identify 413 independent seismic shaking markers from the Dead Sea center.

Ground acceleration constraint for deformed layers

We update the previous computational fluid dynamics modeling of Wetzler et al. (18) by extending the upper limit of layer thickness and ground acceleration from 0.5 m and 0.6g to 1.0 m and 1.0g, respectively (Fig. 3; Materials and Methods). We run a series of two-dimensional numerical simulations using the physical properties of the soft sediments at the bottom of the Dead Sea based on the Kelvin-Helmholtz instability mechanism. According to the numerical simulations, formation of a layer of (i) linear waves, (ii) asymmetric billows, (iii) coherent vortices, and (iv) intraclast breccia requires a minimum acceleration of 0.13, 0.18, 0.34, and 0.50g, respectively (Fig. 3). Also, our numerical simulations indicate that the thickness of the deformed layer scales with acceleration (Fig. 3 and fig. S4). Because we could not directly determine the Reynolds number of an individual unstable layer at the time of seismic shaking, we constrain only the lower boundary of acceleration needed to initiate deformation of a layer with a certain thickness. By considering both the shape and thickness of the deformed layers, we identified 18, 67, 139, 141, 240, and 413 events with acceleration of ≥0.65g, ≥0.50g, ≥0.34g, ≥0.26g, ≥0.18g, and ≥0.13g, respectively (table S2).

Fig. 3 Numerical simulation on in situ folded layer and intraclast breccia structures in the Dead Sea sedimentary sequences.

(A) Typical structures from Dead Sea depocenter Core 5017-1. (B) Typical structures from Dead Sea onshore outcrops (Fig. 1B). (C) Schematic diagrams based on snapshots from the numerical simulations demonstrating the four structures. (D) Quantitative estimation of the accelerations that are needed to initiate the four structures with different thicknesses; the deformations normally occurred when Richardson number ≤ 0.125.

Return time statistics

We calculate the timing of the seismites between dated horizons by linear interpolation between dated points of the cored section. Previous computational fluid dynamics modeling indicates the sediment turbulence develops at the interface of water and water-saturated sediments, and one earthquake produces only one in situ deformation (17, 18). Therefore, the age of the first undisturbed lamina overlying the deformed layer constrains the time of formation of each seismite. The 413 acceleration ≥0.13g events that Core 5017-1 records have a mean return time of 530 ± 40 (SEM) years and an SD of 900 years (Fig. 4A and Table 1). The strong seismic shaking events (acceleration, ≥0.34g) that Core 5017-1 records have a much longer mean return time of 1500 ± 190 years and an SD of 2200 years. The acceleration ≥0.13, ≥0.18, ≥0.26, ≥0.34, ≥0.50, and ≥0.65g events have the coefficient of variation (COV; SD divided by mean) of 1.6, 1.7, 1.5, 1.5, 1.6, and 1.5, respectively. Earthquake distributions with COV ≤ 0.7 are “quasi-periodic,” distributions with COV around 1 are random, while distributions with COV > 1 are “clustered” (aperiodic) (2, 27). Thus, the events that Core 5017-1 records at different shaking intensity levels are clustered in time, which is in line with some previous paleoseismic investigations in the region (4, 5). Besides, the acceleration ≥0.13g and ≥0.34g events present a power-law–like probability density of the recurrence interval (Fig. 4, B and C, and Table 1).

Fig. 4 Return time statistics of seismites and magnitude constraint for strong seismic shaking events during the past 220 ka.

(A) Temporal distribution of moderate (PGA ≥ 0.13g) and strong (PGA ≥ 0.34g) seismic shaking events. (B and C) Histograms for return times of PGA ≥ 0.13g and PGA ≥ 0.34g events. We plot two distribution types (exponential and power-law) for each dataset. (D) Normalized return time data to show return time distribution of moderate and strong seismic shaking events. (E) Magnitude constraint for strong seismic shaking events by applying the three regional empirical attenuation relations (2831), taking the past 2-ka earthquake scenario as an analogy for the paleoseismic record, and assuming that most Mw ≥ 6.0 earthquakes occurred with D ≥ 30 km from the drilling site (see the text for details); D, epicentral distance. (F) Comparison of different temporal distributions of large earthquakes on the central Dead Sea Fault Zone derived from three different geological records. (G) Magnitude-frequency distribution of modern (33) (olive colored) and paleoearthquakes (pink colored) on the central Dead Sea Fault during the past 220 ka; the number of modern earthquakes (fig. S5) is extrapolated to 220 ka.

Table 1 Statistical analysis of recurrence times for the referred records.

View this table:


Lower-bound magnitude determination for the paleoevents

Because of the intrinsic geometrical spreading, dissipation, and possibly dispersion of the seismic energy, ground motion effects of a moderate earthquake nearby generate similar shaking intensities to those generated by a large earthquake farther away. Therefore, determining the location of an earthquake along the length of the Dead Sea Fault or on other nearby faults is the key for the magnitude constraint. It follows that magnitude estimation based on a single station is difficult and dependent on the location of possible source faults. For this sinistral boundary between the African and Arabian plates, we model the potential source region of earthquakes as having a fixed width and a length of a few hundred kilometers (Fig. 1). The nearest faults are ~5 km from the Core 5017-1 drilling site.

We apply three empirical attenuation relations (2831) developed for the Dead Sea region to constrain the lower magnitude limit of paleoseismic events that the ICDP Core 5017-1 records (Materials and Methods). Among them, two attenuation relations (2830) are described by macroseismic intensity, magnitude, and epicentral distance (D), and one relation (31) is described by peak ground acceleration (PGA), magnitude, and epicentral distance. To apply the two attenuation relations that are described by seismic intensity, we convert accelerations into seismic intensity via the linear relationships between PGA and the modified Mercalli intensity scale (MMI) that Wald et al. (32) proposed (table S2). According to the three regional empirical attenuation relations, PGA ≥ 0.13g or MMI ≥ VI½ corresponds to Mw ≥ 5.5, 5.3, and 5.7, by taking Dmin = 5 km (table S3). Therefore, we interpret the lower-bound magnitude of the recorded PGA ≥ 0.13g (MMI ≥ VI½) events as Mw ≥ 5.3. Considering that the incompleteness of moderate earthquakes in the paleoseismic record, we infer that the mean recurrence of Mw ≥ 5.3 earthquake on the central Dead Sea Fault Zone is <530 ± 40 years. This value is much shorter than the previously obtained mean recurrence of ~1600 years for the same magnitude, based on the paleoseismic record between 60 and 14 ka ago (4).

Magnitude constraint for strong seismic shaking events

The ICDP Core 5017-1 records 139 strong seismic shaking events with PGA ≥ 0.34g (MMI ≥ VIII) (table S2). According to the regional empirical attenuation relations (2831), one felt intensity, for example, MMI ≥ VIII (PGA ≥ 0.34g) at the drill site, will require a moderate earthquake (6.0 < Mw < 7.0) with a D between 5 and 30 km, an Mw ≥ 7.0 earthquake with a D of 30 km, or an Mw ≥ 8.0 earthquake with a D of up to 150 km. We adopt a maximum magnitude in the region of Mw 8.0, which will require a rupture of ~300 km along a major fault. Therefore, we consider only large earthquakes with D ≤ 150 km (the central Dead Sea Fault Zone) as triggers of the strong seismic shaking events (MMI ≥ VIII) that Core 5017-1 records. On this slow-slipping plate boundary, the Dead Sea Fault (zone) is the only real contributing fault because most of the transform margin slip rate is on the Dead Sea Fault, and the other faults were either far away or had low slip rates. Within the Dead Sea region, no other major faults have sufficient length and slip rate to materially contribute to the rate of Mw ≥ 7 earthquakes. The historic earthquake catalog from the region shows that the 1822, 1712, 1408, 1170, 1139/1140, and 859/860 CE Mw ≥ 7.0 earthquakes occurred with distance to the drill site ≥300 km north of the Dead Sea (the northern Dead Sea Fault Zone). However, none of them have an expression in the Dead Sea Core 5017-1 record.

The spatial distribution of instrumental and historic moderate and large earthquakes on the central Dead Sea Fault Zone during the past 2 ka do supply additional clues for magnitude constraint for these strong seismic shaking events. The instrumental (33) and historical (5, 3436) earthquake catalogs reveal that during the past 2 ka, all major earthquakes (Mw ≥ 6.0) occurred with D ≥ 30 km from the drilling site (Fig. 1B). By taking the past 2 ka earthquake scenario as an analogy for the paleoseismic record, we assume that most Mw ≥ 6.0 earthquakes occurred with D ≥ 30 km from the drilling site. Under this basic assumption and the three regional empirical attenuation relations, (i) an intensity of MMI ≥ VIII (PGA ≥ 0.34g) requires an earthquake with Mw ≥7.0, ≥7.0, and ≥7.3; (ii) an intensity of MMI ≥ VIII½ (PGA ≥ 0.50g) requires an earthquake with Mw ≥7.4, ≥7.3, and ≥7.6; and (iii) an intensity of MMI ≥ IX (PGA ≥ 0.65g) requires an earthquake with Mw ≥ 7.8, ≥7.6, and ≥7.8 (Fig. 4E and table S3). Therefore, we interpret the corresponding lower-bound magnitudes of strong seismic shaking events in the ICDP Core 5017-1 to be Mw ≥ 7.0, 7.3, and 7.6, respectively.

We test our magnitude conversion versus known historic earthquakes on the central Dead Sea Fault Zone. Six seismites in Core 5017-1 dated at −2 ± 44 years before the present (yr B.P.), 42 ± 44 yr B.P., 148 ± 44 yr B.P., 1248 ± 44 yr B.P., 1555 ± 47 yr B.P., and 1626 ± 47 yr B.P. correspond to the 1956 CE (Mw 5.5; D, ~5 km), 1927 CE (Mw 6.25; D, ~30 km), 1834 CE (Mw ~ 6; D, ~60 km), middle 8th century (Mw > 7; D, ~100 km), 419 CE (Mw ~ 6; D, ~40 km), and 363 CE (Mw ~ 6.8; D, ~70 km) earthquakes, respectively (table S4) (36). According to the regional empirical attenuation relations, we constrain the magnitudes of the six paleoearthquakes (seismites) with intensities of VII (0.18g), VI½ (0.13g), VI (0.09g), VII (0.18g), VI (0.09g), and VII (0.18g) as Mw 5.6, Mw 6.1, Mw 6.2, Mw 7.1, Mw 6.0, and Mw 6.9, respectively, which are in line with recorded historic magnitudes (table S4). This test supports our magnitude conversion based on the regional empirical ground motion attenuation relations.

Integrate a 220-ka-long large earthquake record

A few short gaps may exist in the paleoseismic record based on the ICDP Core 5017-1, due to low coring recovery rate and halite deposition at some depth. In addition, there is the possibility that some remote Mw ~ 8.0 earthquakes with D > 150 km induced the PGA < 0.34g (MMI < VIII) events. We incorporate previously established large earthquake records from the central Dead Sea Fault Zone into our present 220-ka paleoseismic series to develop two constraining limitations. Begin et al. (19) developed a 60-ka-long Mw ≥ 7.0 paleoseismic record based on seismites identified from Dead Sea western onshore outcrops (Peratzim Valley; Figs. 1B and 4F). In addition, Kagan et al. (20) and Kagan (21) developed a 185-ka-long Mw ~8.0 paleoseismic record based on damaged cave deposits in the Soreq and Har-Tuv Caves, located 40 km due west of the Dead Sea Fault (Figs. 1B and 4F). The long record comprises 26 events that occurred in the past ~185 ka.

Because the paleoseismic record of Begin et al. (19) is too short, we incorporate the paleoseismic record of Kagan et al. (21) into our present 220-ka paleoseismic series. Among the 26 events, four speleoseismites correspond to PGA < 0.34g (MMI < VIII) seismites (in situ deformed layers) in Core 5017-1, and four speleoseismites have no corresponding seismites in the drilling (table S5). We, therefore, incorporate these eight speleoseismites into the 220-ka paleoseismic series. In addition, we incorporate four Mw ≥ 7.0 historic earthquakes that occurred near the drilling site (D < 150 km) during the past 2 ka (Fig. 1C) into the 220-ka paleoseismic series, because these four historic earthquakes either correspond to PGA < 0.34g (MMI < VIII) events or do not show discernible responses in Core 5017-1 (table S5). In total, we incorporate 12 Mw ≥ 7.0 events (including 8 Mw ~ 8.0 events) into the 220-ka paleoseismic series (Fig. 4F). Thus, the integrated 220-ka-long large earthquake record comprises 151 Mw ≥ 7.0 events, 75 Mw ≥ 7.3 events, and 26 MMI ≥ IX (PGA ≥ 0.65g) events. The integrated record yields a mean recurrence of 6200 ± 1800 years for the 26 largest events, similar to Kagan et al.’s (21) mean recurrence of 6900 ± 1000 years for Mw ~8 earthquakes derived from the speleoseismite record. We therefore shift the magnitude corresponding to MMI ≥ IX (PGA ≥ 0.65g) events from Mw ≥ 7.6 to Mw ≥ 7.8 (table S3), which is also in line with the magnitudes converted from the other two regional empirical attenuation relations (28, 31).

A previous compilation of instrumental, historic, and paleoseismic studies inferred that the Gutenberg-Richter distribution has been stable in the Dead Sea region during the past 60 ka, with a b value of ~0.95 (37), and similar to a b value of 0.97 deduced from instrumental dataset (Fig. 4G). The integrated 220-ka-long large earthquake record shows a b value of 0.95 (Fig. 4G, magenta color), assuming a Gutenberg-Richter distribution for the Dead Sea region as a whole. This similar b value of the Gutenberg-Richter distribution and approximate match in a value after matching the catalog time span to the 220-ka total length implies that our magnitude constraint for strong seismic shaking events is appropriate, and our record of large earthquakes is relatively complete.

Earthquake recurrence pattern and fault behavior model and their implications for seismic hazard

The paleoseismic record of Begin et al. (19) yields a mean recurrence of 4600 ± 1500 years for Mw ≥ 7.0 earthquakes and a COV of 1.0 during the past 60 ka (Fig. 4F and Table 1). The exponential and power-law distributions, however, are not fit well to the distribution of these recurrence times (Fig. 5A). In contrast to the paleoseismic record of Begin et al. (19), the paleoseismic record of Kagan et al. (21) yields a mean recurrence of 6900 ± 1000 years for Mw ~8.0 earthquakes and a COV of 0.7 (Fig. 4F and Table 1). The earthquake recurrence times approximately follow a Weibull distribution (Fig. 5B). This type of distribution, together with the small COV, suggests a quasi-periodic earthquake recurrence pattern for the largest earthquakes.

Fig. 5 Recurrence pattern of large earthquakes (Mw ≥ 7.0) on the slow-slipping plate boundary during the past 220 ka.

(A) Histograms for recurrence times of Mw ≥ 7.0 events in the paleoseismic record of Begin et al. (19). (B) Histograms for recurrence times of Mw ~ 8.0 events in the paleoseismic record of Kagan et al. (21). (C) Histograms for recurrence times in the integrated 220-ka-long Mw ≥ 7.0 record. (D) Normalized recurrence data for the three datasets for comparison.

The integrated 220-ka-long Mw ≥ 7.0 earthquake record yields a mean recurrence of ≤ 1400 ± 160 years and a COV of 1.4 (Fig. 4F and Table 1). Recurrence times follow a power-law distribution (Fig. 5C). In contrast to the recurrence patterns revealed by the abovementioned two large earthquake records, the integrated 220-ka-long Mw ≥ 7.0 record shows a clustered recurrence pattern (Fig. 5D and Table 1). The distribution of earthquake recurrence times shows a maximum probability density at the shortest recurrence times (< 0.7 ka) and exhibits a long tail at the longest recurrence times (> 5.5 ka). The high probability density at the shortest recurrence times may represent intra-cluster recurrence, and the long tail indicates intercluster recurrence in the group-fault temporal clustering model of McCalpin (6). The group of faults in this case may be the bounding faults of the Dead Sea Basin, plus ruptures within the basin itself. For example, ruptures entering the basin from the south would cause strong shaking at the drilling site and cluster by promoting rupture to the north on the west-bounding fault. Additional Mw ≥ 6 events could be required in the overlap zone to maintain the net plate boundary slip rate. Lower apparent rates of large events outside the basin may correspond to the longer intercluster times of the ICDP Core 5017-1 record.

The specific locations of the three study sites (Peratzim Valley, Soreq and Har-Tuv Caves, and ICDP Core 5017-1; Fig. 1B) could partially explain the different earthquake recurrence patterns revealed by the three large earthquake records. The Peratzim Valley is situated west of the southern Dead Sea Basin. The seismites preserved at the site may primarily record earthquakes originating from the southernmost end of the central Dead Sea Fault and thus lead to the aperiodic recurrence pattern with a much longer mean recurrence. The Soreq and Har-Tuv Caves located on the northwestern side of the Dead Sea Basin tend to be more frequently affected by the movement of the northern part of the central Dead Sea Fault. This location may not record some Mw ~ 8.0 earthquakes, which extend south from the southern Dead Sea Basin. The speleoseismite record, therefore, may correspond more closely with nearby faults and in a quasi-periodic recurrence pattern. The ICDP Core 5017-1, with a central location and much higher sediment accumulation rate, may have recorded earthquakes from all faults that surround the drilling site. The group-fault temporal clustering will appear if each quasi-periodic behavioral fault is out of phase with each other. Alternatively, the group-fault temporal clustering might arise from earthquake rupturing and stress transfer between different faults near the drilling site over short time intervals and then followed by long quiescence.

The unique location of the ICDP Core 5017-1 makes the integrated 220-ka-long Mw ≥ 7.0 earthquake record the most complete on the Dead Sea Fault and the most representative of the fault history. The integrated record yields a mean recurrence of ≤ 1400 ± 160 years for Mw ≥ 7.0 earthquakes, which is significantly shorter than the mean recurrence of 4600 ± 1500 years for the same magnitude based on Begin et al.’s 60-ka-long paleoseismic record (19). Moreover, our integrated record shows that Mw ≥ 7.0 earthquakes are clustered during the past 40 ka (COV = 1.4), with a mean recurrence of ≤ 700 ± 110 years. Thus, the new results do not support the previous conclusion that the seismic regime in the Dead Sea Basin has been stationary for the past 40 ka and characterized by periodic recurrence pattern with a mean recurrence of ~11 ka during that time period (19). Therefore, our integrated 220-ka-long Mw ≥ 7.0 record indicates a higher seismic hazard than previously appreciated for the slow-slipping Dead Sea Fault.

Different from the periodic recurrence of earthquakes on fast-slipping and geometrically simple strike-slip faults, e.g., the Alpine Fault in New Zealand (2), we infer aperiodic earthquake behavior on the slow-slipping and the geometrically complex sinistral boundary between the African and Arabian plates. Our results suggest that researchers may underestimate the seismic hazard potential of similar slow-slipping faults with irregular rupture. Our study highlights the potential of in situ deformed sediment layers in a subaqueous environment as a strong-motion paleoseismometer to record long seismic sequences covering multiple recurrence intervals of large earthquakes. Long records are vital for accurate hazard assessment. Our quantitative method of seismic record reconstruction, with paleoearthquake intensity (ground acceleration) and magnitude estimation, may also prove suitable for similar subaqueous environments along other faults.


ICDP Core 5017-1 and age model

We retrieved the ICDP Core 5017-1 from the Dead Sea depocenter, with a water depth of 297 m, from November 2010 to March 2011. The 457-m-long composite core penetrated Holocene and late Pleistocene sediments. Mud, alternations of aragonite-detritus laminae, and halite comprised the sedimentary sequence. The millimeter- to submillimeter-scaled detritus laminae are deposited during the rainy seasons and are composed of allogenic quartz, calcite, and clay. The delicate aragonite laminae are composed of authigenic aragonite crystals (38). The core was dated with methods of 14C, U-Th, tuning, and δ18O stratigraphy correlation (22, 23). Samples for 14C dating focused on terrestrial plant remains to avoid a potential inheritance effect from saline water. The measured 14C ages were calibrated by using OxCal 4.3 and IntCal13 with 1σ error. The U-Th dating was carried out on primary aragonite laminae and reported with 2σ error. The age of Core 5017-1 ranges from ~220 ka ago to the present. In total, we use 53 age points for the age model (fig. S1 and table S1). We use linear interpolation between dated layers to calculate the timing of each seismite. We did not consider it necessary to propagate uncertainties of the age determination to each seismite. Previous studies (2, 3) have shown that when analytical uncertainties from dating are small relative to the time between events, dating uncertainty has little effect on parameter estimates such as COV and virtually no effect on mean recurrence time.

Computational fluid dynamics modeling

Heifetz et al. (17) first modeled the physical process of soft-sediment deformation preserved in the Dead Sea sediments. This linear fluid dynamics modeling found that the Kelvin-Helmholtz instability is a plausible mechanism for observed soft-sediment deformations. Subsequently, Wetzler et al. (18) improved the modeling by considering nonlinear processes of soft-sediment deformation. Wetzler et al. (18) examined more than 300 soft-deformation structures, which include the linear waves, asymmetric billows, coherent vortices, and intraclast breccia textures preserved in the late Quaternary Dead Sea sedimentary sequences, and found that the geometry closely follows a power-law (with an exponent of −1.9), compatible with Kelvin-Helmholtz turbulence in other environments. The viscosity and density of modern Dead Sea sediments in their water-saturated condition are taken as a reasonable analog for late Quaternary Dead Sea sediments and measured for modeling purposes. Boundary conditions of the model were designed to reproduce different values of ground acceleration. Minimum ground accelerations needed to induce the layers of (i) linear waves, (ii) asymmetric billows, (iii) coherent vortices, and (iv) intraclast breccia layers with different thicknesses are determined.

In this study, we update the previous computational fluid dynamics modeling of Wetzler et al. (18) by extending the upper limit of layer thickness and ground acceleration from 0.5 m and 0.6g to 1.0 m and 1.0g, respectively. We model the flow using the commercial software Fluent (, a computational fluid dynamics modeling tool. Wetzler et al. (18) provide more information on the computational fluid dynamics modeling.

The effects of water depth and secondary deformation were also considered. The water depth itself does not directly affect the dynamics of Kelvin-Helmholtz instability between two mud layers in the lakebed because they are in near hydrostatic balance before the deformation. Thus, only the difference in densities between the layers bears on the dynamics. The relative motion of the water layers with respect to the mud delivers the energy via an internal gravity wave. The dispersion relation for such waves is controlled by the thickness of the deforming muddy layer h, up to an order of a meter. The eigenfrequency (f) isf=sqrt[g(ρm+ρb)/h(ρmρb)](1)where g is the effective acceleration of gravity, ρm is the density of the minerals, and ρb is the density of the brine. With ρmb of ~2 and brine content of around a third, f is approximately 1 Hz.

Earthquake ground motion attenuation in the Dead Sea region

Three representative empirical ground motion attenuation relations for the central Dead Sea Fault are available from the literature.

Intensity (I)–magnitude (M)–D relations

(i) Malkawi and Fahmi (28) developed empirical intensity (I)–magnitude (M)–D relations characterizing earthquake ground motion attenuation in Jordan and conterminous areas, based on instrumental and historic datasets from the region. We chose to use their historic attenuation relation instead of the instrumental relation because the instrumental dataset does not include any large earthquakes, while the historic dataset includes abundant large earthquakes. The equation is expressed asI=5.76+1.52MS2.09 ln(D+25)(2)

When applying the relation, we convert the surface-wave magnitude (MS) scale into Mw scale according to (39, 40).

(ii) Hough and Avni (29) published a large-magnitude attenuation equation for the Dead Sea region calibrated using 133 seismic intensity (I) determinations of the 1927 Mw 6.3 Jericho earthquakeI=0.64+1.7M0.004D1.67log(D)(3)

Darvasi and Agnon (30) improved the attenuation relation of Hough and Avni (29) by considering local near-surface property (shear wave velocity; Vs30) as an amplification factor. The calibrated attenuation relation is expressed asI=0.64+1.7M0.004D1.67log(D)2.1ln(Vs30/655)(4)

In this study, we apply this calibrated attenuation relation (Eq. 4) instead of equation 3 of Hough and Avni (29). According to previous shear wave measurements on the Dead Sea sediments (41), we use a value of ~900 m/s for the Vs30 index. Following the attenuation relations (Eqs. 2 and 4), we estimate the change in M with D by fixing I.

PGA-M-D relation

On the basis of 57 PGA recordings of 30 instrumental earthquakes that occurred between 1979 and 2001 on the central and southern Dead Sea Fault, Al-Qaryouti (31) proposed an attenuation relation of PGA (g)log(PGA)=3.45+0.50M0.38log(D)0.003D±0.31P(5)in which P represents the normal distribution and P = 0 and 1 for the 50th and 84th percentiles, respectively. Following the attenuation relation, we estimate the change in M with D by fixing PGA.

These empirical attenuation relations provide the opportunity to constrain magnitudes for our paleoseismic events recorded in the ICDP Core 5017-1. However, we note that all three empirical attenuation relations have some shortcomings when applying them to our 220-ka paleoseismic record. In the attenuation relation of Malkawi and Fahmi (28), 46 historic earthquakes were used to constrain the magnitude scaling, far more than the other two attenuation relations. However, only 71 data points of maximum intensity were included, and some magnitude uncertainties might exist in the historic earthquakes (between 19 and 1896 CE). The attenuation relation of Hough and Avni (29) has robust intensity-distance scaling and fit for the 1927 Mw 6.3 Jericho earthquake because 133 intensity points were included. However, the attenuation relation was based on only a single instrumental moderate earthquake. The attenuation relation of Al-Qaryouti (31) is based on instrumental PGA values, which provide high measurement accuracy. However, the number of PGA points is limited because the observation time period is relatively short (1979–2001). In addition, much of the data come from single-site recordings. In the same year, Meirova et al. (42) published an alternative PGA-M-D relation but limited to data west of the Dead Sea fault. It predicted much higher magnitudes than the relation of Al-Qaryouti (31) and was not applied to our present study.

We applied all these three empirical attenuation relations (28, 30, 31) to constrain the lower magnitude limit of paleoseismic events recorded in the ICDP Core 5017-1. The attenuation relation of Malkawi and Fahmi (28) is more suitable for the inverting point intensity of strong seismic shaking events (MMI ≥ IX) to magnitude than the others. This is because the catalogs of Hough and Avni (29) and Al-Qaryouti (31) used to reconstruct empirical attenuation relations lack large earthquakes (Mw ≥ 7.0). We estimate that magnitude uncertainties of historic earthquakes in the catalog of Malkawi and Fahmi (28) do not have much impact on the strong seismic shaking intensity range relevant to our present study.

b value and magnitude of completeness for the modern earthquake catalog

We calculate the magnitude of completeness (Mc) for a rectangle of an area 150 km north and south to the drilling site, between 1990 and 2015 (fig. S5). Mc is computed iteratively from the goodness of fit using a Kolmogorov-Smirnov test between observed and theoretical Gutenberg-Richter distributions with b value from maximum likelihood estimation. We select the smallest magnitude for which the difference between the distributions is negligible. More specifically, we use a threshold of ΔD = 0.05 as the definition of negligible when ΔD is the Kolmogorov-Smirnov metric distance between the two cumulative distribution functions, resulting in Mc = 3.1 and b value = 0.97.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We appreciate the three anonymous reviewers for thorough, constructive, and encouraging review and their strong recommendation. We thank J. Hall for English editing. We also thank Z. Guo, Q. S. Liu, and Y. John Chen for supporting Y.L.’s 3-week visit during September 2017 at the Department of Ocean Science and Engineering, Southern University of Science and Technology, China, to continue working on the project. Funding: This research was supported by the University of Liege under Special Funds for Research, IPD-STEMA Program (R.DIVE.0899-J-F-G to Y.L. between 2019 and 2020), Post-Doctoral Fellowship of the Faculty of Exact Sciences at Tel Aviv University (to Y.L. between 2016 and 2017), the Israel Science Foundation (Center of Excellence grant #1436/14 and grant #1645/19 to S.M.), the International Continental Scientific Drilling Program (ICDP), the DESERVE (Dead Sea Research Venue) Virtual Institute under the auspices of the Helmholtz Association ( (to A.A.), and the Israel Science Foundation (ISF 363/20 to N.We.). Author contributions: Y.L. put forward the path to the goal, developed the geological observations and interpretations, and prepared the manuscript, which was edited by all the coauthors; N.We. contributed to the numerical simulation; N.Wa. participated in the ICDP campaign and the associated sampling and logging party together with A.A., who also contributed to the interpretation of historic earthquake catalog and the inverting point intensity for magnitude-frequency relations; G.P.B. contributed to the inverting point intensity for magnitude-frequency relations and English editing; S.M. initiated and contributed to all aspects of the project. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Correspondence and requests for additional data should be addressed to Y.L. (yinlusedimentology{at};{at}

Stay Connected to Science Advances

Navigate This Article