The Impact of Lake-Level Fluctuation on Earthquake Recurrence Interval over Historical and Prehistorical Timescales: The Case of the Dead Sea

Abstract

We review the impact of large historical lake water-level changes on seismicity via the stress field of the shallow crust where devastating earthquakes nucleate. A novel backward earthquake simulation presented in this chapter can be used to investigate the geological record for the past ten millennia (presented in this study) and even more. The simulation is based on a theoretical model, which explains the variability in the recurrence interval of strong earthquakes. We suggest that the water-level changes in ancient lakes located in tectonic depressions along the Dead Sea transform could contribute to the observed differences. It is found that the increase in the water level moderates the seismic recurrence interval. Based on this empirical correlation together with mechanical considerations, an additional indication is established regarding the water-level reconstruction and location of earthquakes in the Dead Sea area. This indication is based on simulated earthquakes, by superimposing the effective normal stress change due to the reconstructed water-level change on the estimated tectonic shear stress accumulated since the preceding seismic event.

1. Introduction

Earthquakes are one of the primary geological hazards, with historical and modern accounts serving as vivid testimony of the threats they pose. With catastrophic losses of human lives and health [1, 2], earthquakes may cause major economic losses and property damage [3, 4]. Earthquakes also play a significant role in sedimentary systems evolution, liquefying water-saturated clastic sediments by cyclic loading [5, 6, 7].

Despite the global, social, and scientific impact of earthquakes, in many cases, their triggering mechanisms remain poorly defined. Specifically, tectonic motion is the primary cause and consequence of earthquakes: For example., large earthquakes along transform faults are typically associated with slip rate and stress buildup and release (e.g., Dead Sea transform (DST)—[8, 9, 10, 11, 12, 13], San-Andreas fault—[14, 15] or North Anatolian fault—[16, 17, 18]).

However, earthquake triggers may have no direct relation to the local tectonic dynamics. Observations supported by theoretical calculations demonstrate that surface processes (e.g., glaciations, sedimentation, and water loading) are able to change significantly the stress field in the Earth’s crust [19, 20]. Field data from northern Europe and America suggest that melting of the ice sheets was connected to the paleo-seismologic evidence and to the development of late Pleistocene/early Holocene fault scarps [21, 22, 23, 24], whereas the seismicity in the tectonically stable Greenland and Antarctica may be solely due to the current ice unloading [25].

Modeling the slip-rate evolution of localized normal faults located below the glaciated Earth crust, Hetzel and Hampel [26], Hampel and Hetzel [27], and Hampel et al. [28] demonstrate that variations in the slip-rate in these faults depend on the changes in the differential stress in the crust, which are in turn connected to changes in the loading on the Earth surface. The glaciations induce a decrease in the differential stress in the crust and, thus, a decrease in the faulting rate. In contrast, during postglacial rebound and isostatic adjustment, the differential stress in the crust grows, resulting in increased rates of faulting, inducing in turn temporal clustering of the earthquakes.

Additionally, it has long been recognized that fluids can influence the occurrence of earthquakes. Injection of liquid into (or removal from) the subsurface is capable to induce seismicity [29, 30, 31, 32] even outside the region subjected to the direct fluid perturbations [33, 34, 35].

Moreover, seismic response in the field was observed to be connected to reservoir impoundments and their seasonal water level variations [36, 37, 38, 39, 40]. An increase in pore pressure may lead to two different kinds of seismic response: swarm-like rapid seismicity induced by immediate pore pressure buildup, and delayed seismicity associated with pore pressure diffusion from the reservoir to hypocentral depths [37, 41]. Delayed response, as observed in the field, is often associated with large magnitude earthquakes [37, 42]; it may extend significantly beyond the confines of the reservoir and may not show an immediate correlation with major changes in reservoir level.

Seismic activity change was associated with historical water-level changes in the Dead Sea lake [43, 44] and Salton Sea lake [45, 46, 47] at much longer timescales. The Salton Sea, a remnant of ancient Lake Cahuilla, lies on the southern San-Andreas Fault (SSAF) [48]. Luttrell et al. [47] conducted a study to examine how Coulomb stress perturbations influenced local faults in relation to changes in the level of Lake Cahuilla. Their objective was to assess whether the lake could potentially impact the timing of fault rupture. The researchers proposed that the loading resulting from the ancient Lake Cahuilla impeded fault failure along most of the SSAF, except for a potential small area within the lake boundary.

Subsequently, Brothers et al. [46] modeled changes in Coulomb stress due to the combined effects of lake loading caused by periodic flooding of Lake Cahuilla, movement of stepover faults, and increased pore pressure. They determined that these factors could elevate stress levels along the southern SSAF to the extent of triggering an earthquake. Over the past 1100 years, the periodic flooding of Lake Cahuilla has exhibited fluctuations of up to 100 meters. These variations have shown throughout this time period a significant correlation with the incidence of earthquakes along the southern, the least active segment, of the SAF.

Belferman et al. [44] took advantage of the high temporal resolution of the historical earthquake records for the Dead Sea Lake during the last two millennia. The relatively small water-level hikes of 15 m, characteristic for time intervals of centuries to millennia, were analyzed and shown to be able to exacerbate the seismicity pattern at the Dead Sea fault. Belferman et al. [43] demonstrated for the first time that plausible scenarios for lake-level history show striking correlations with the historical record of earthquake recurrence intervals. Additionally, the plausible correlation observed between these phenomena suggests a complementary indicator that could be used to constrain either the locations of historical earthquakes or lake-level fluctuations during periods of uncertainty in one of these variables.

More recently, Hill et al. [45] present a direct relationship between changes in Lake Cahuilla water level and earthquakes based on computed time-dependent Coulomb stress changes for the last 1100 years. The authors present coupled 3D finite element model that takes into account time-varying surface water loads, fault geometry, crustal poroelasticity, and viscoelastic relaxation in the ductile substrate. The simulation results indicate that the past six major earthquakes on the southern San-Andreas Fault probably occurred during high stands of Lake Cahuilla. It was shown that increases in lake-level result in positive Coulomb stress changes on most of the southern San-Andreas Fault, with stressing rates as high as two to three times the tectonic loading rate.

In this chapter, we summarize the results and methods to identify the correlation between the historic water-level reconstructions at the Dead Sea and seismicity patterns in the area over the past two millennia. We also offer an additional analytical model to cope with dependence on initial conditions, which can be used in the future for interpreting these phenomena for earlier extended periods.

2. The study area: Dead Sea

2.1 Tectonics

The Dead Sea Lake fills tectonic depressions along the Arabian and African plate boundaries. The Arabian Plate is rifting away from the African and Nubian plates along the active spreading axis of the Red Sea, with the Sinai-Levant Block lagging behind, sheared along the Dead Sea fault (Figure 1) [49, 50, 51].

Figure 1.

Tectonic plates in the Middle East. The Dead Sea fault is one of the Dead Sea transform (DST) fault system, which transfers the opening at the Red Sea to the east Anatolian fault. Jordan Valley (JV), Dead Sea (DS), Arava/Araba Valley, and Gulf of Aqaba (GoA) located along the DST.

The Dead Sea fault is the major fault system along the plate boundary, and one of the best studied continental transform faults in the world [8, 11, 12, 49, 52, 53, 54, 55]. The tectonic motion at the DST is characterized predominantly by a left-lateral strike-slip regime with a velocity of up to ∼5 mm/yr. along various segments [51, 56, 57, 58]. Motion is transferred from the opening Red Sea to the escaping Anatolian Block via sinistral shear motion along the Dead Sea transform (e.g., [49, 51, 59]). This sinistral strike-slip plate boundary comprises the basins of the Gulf of Aqaba/Elat, Araba/Arava, the Dead Sea, the Jordan, Hula, Beqa’a, Ghab, and Karasu, a distance of about 1000 km [60].

Many factors and processes operating on different scales may have a profound influence on the dynamics of the deformations in the Dead Sea region. One of them—the influence of the water bodies deep underneath the Jordan-Arava valley—is explored in depth in this study.

2.2 Seismicity

The Dead Sea fault system is a major source of prehistoric, historical, and contemporary earthquakes [8, 11, 61, 62, 63, 64, 65]. Tectonic motion on the DST has always been addressed as a primary cause of earthquakes in the region. Large earthquakes (M > 5.5) are typically associated with slip rate and stress buildup at the strike-slip faults of the transform [9, 10, 11, 13, 66, 67, 68, 69]. Yet, normal faults might release M > 6 earthquakes [64].

Due to paleo-seismological research, we have access to 220 kyr of earthquake records, obtained based on different dating methods (e.g., uranium-thorium dating, radiocarbon dating, infrared stimulated luminescence dating) of deformation or changes that occur in rock structures due to earthquake shaking (namely seismites) [65]. A promising environment for paleoseismological investigations is provided by caves because they serve well as repositories for geological information. Kagan et al. [70] constrained dates of 38 seismites from the Soreq and Har-Tuv caves, 40 km of the DST. This was based on interpretation and amalgamation of individual seismite ages for bracketing the chronology of the strongest (M7.5-8) earthquakes covering the time interval since 185 ka.

Kagan et al. [70] found in their study a correlation between the paleoseismic events that occurred during the past 75 kyr identified with speleothems with other independent paleoseismic records in the region. Some of them were recorded in Dead Sea cores [71] and correlate with damage at archeological sites [8, 72]. Others correlate with seismites in the Lisan Formation [73, 74] or in the pre-Lisan Formations [65, 75].

Earthquake imprints formed in laminated sediments of the Lisan Formation [76] such as intraclast breccias [7], represent an uninterrupted geological record of the earthquake in the Dead Sea region for the past 14 to 70 kyr. Porat et al. [77] inferred that injection of clastic dikes is one of the most impressive liquefaction features generated by strong, M ≥ 6.5 earthquakes. The authors constrained the liquefaction date to shortly after the deposition of the Lisan Formation (15 kyr) by using optically stimulated luminescence ages. The dikes ages vary between 15 and 10 ka, thus suggesting that the period following the draw-down of Lake Lisan period was accompanied by strong earthquakes.

Migowski et al. [71] reconstructed an earthquake record using the dated Ein Gedi and Ze’elim cores seismites from the Dead Sea shores for the past 10 kyr. Later, for the last 4000 years, significantly more seismites were detected from the Ein Feshkha sedimentary sections in the Dead Sea area [64].

Recently, Liu et al. [65] reconstructed the longest archive of earthquakes from a 250 m core from the Deed Sea bottom, covering 220 kyr. Many of the seismites in the studies mentioned above, dated for the last 2 millennia, are matched by independent historical records [8, 68, 73, 78, 79, 80, 81, 82, 83, 84]. Using geological, archeological, biblical, historical, and seismological data, researchers infer part of the historical records for the Dead Sea area. Ben-Menahem [85] surmised that about 110 major earthquakes affected the area during the past 2500. Of these, 42 originated along the Dead Sea fault system itself, while 68 were imported from the Helenic-Cyprian arcs and the Anatolian-Elburz-Zagros fault systems.

Numerous publications list earthquakes that hit the Dead Sea and its surroundings during the last 2 millennia [8, 61, 62, 63, 79, 80]. Belferman et al. [44] filtered from the scores of listed events only the most destructive ones in the Dead Sea basin, typically causing local intensities of VII or higher in Jerusalem. For a minimal epicentral distance of 30 km, this would translate to a magnitude M > 5.7, according to the attenuation relations of Hough and Avni [86]. Considering site effects in the rocky terrain typical of the Jerusalem vicinity, where the near-surface shear velocity is relatively high, the modeled magnitudes exceed M6 (see Eq. (6) and Table 2: Beit HaKerem and Motza, which can be found in the reference [87]). In Belferman et al. [43], most earthquakes from the earlier list were amended. All earthquakes in Belferman et al. [43] were selected by simultaneously satisfying two criteria: (1) The acceptance regularizes the relation between recurrence intervals and lake level, and (2) they are corroborated by evidence from reassessed historical archives and/or seismites in the northern basin of the Dead Sea (Ein Feshkha and Ein Gedi sites).

2.3 Lake-level fluctuation

Several water bodies occupied the tectonic depression along the Dead Sea transform during the Neogene-Quaternary. The size and composition of these water bodies have changed through time, reflecting the changing climatic [88, 89, 90, 91, 92, 93, 94, 95] and physiographic [90, 96] conditions in the region. Layers of sediments accumulated in the tectonic depression reflecting several lake phases: the Neogene Sedom lagoon [97, 98, 99, 100, 101], the mid- to late-Pleistocene Lake Amora [93, 101, 102, 103], the last interglacial Lake Samra ∼130 kyr [104], the last glacial Lake Lisan 18 ∼ 70 kyr [101, 105, 106], and the Holocene to modern Dead Sea ∼15 kyr [92, 107].

The history of saline water bodies in the Dead Sea basin starts with the Sedom lagoon which has been related to the Neogene ingression of seawater into the Jordan-Arava valley. The reasons for the ingression and the exact age of the lagoon are not clear, and its disconnection from the open sea may be associated with various events that include eustatic changes in the Sea, the tectonic uplift of the Judea-Samaria Anticline, the formation of morphological sill in the Jezreel Valley, etc. [93].

Mid- to late-Pleistocene Lake Amora probably extended over a large part of the Dead Sea basin-Jordan Valley [93, 95]. Zak [101] suggested that the age of the Pleistocene sedimentary sequence at Mt. Sedom deposited by the Amora Formation lies between ∼1 Ma and ∼ 100 ka. The transition between Amora and Samra Formations corresponds to ∼130–140 ka, when the lake stood lower than ∼380 m below mean sea level (bmsl). The lithology and chronology of the Samra Formation suggest that the lake fluctuated mostly between ∼310 and ∼ 350 m bmsl [104], being significantly lower than the last glacial Lake Lisan (mostly ∼280 ± 40 m bmsl) [95, 96, 108, 109]) but higher than the Holocene Dead Sea (mostly ∼400 ± 20 m bmsl) [92, 107, 110, 111]. The lake-level curve for the past 150 kyr, presented in Figure 2, was digitized from Waldmann et al. [104].

Figure 2.

Water-level change in the last 150 kyr, for Samra, Lisan and the Dead Sea (DS) Lake. Derived from the work of Waldman et al. [104].

The most widespread sediments in the rift valley are those deposited at the bottom of Lake Lisan (Figure 2) that occupied the Dead Sea basin and adjacent basins between ∼70 and 18 ka., during the last glacial period in the northern latitudes [105, 112, 113, 114]. The age of the highest lake level of Lake Lisan ∼27–23 ka, [93, 95, 96, 108] or ∼ 26–24 ka [115, 116] corresponds to the time of a cold episode in the global climate, while the relatively low lake level of ∼55–30 ka. Refs. [95, 108, 115] corresponds to the warmer climatic conditions.

The curve of Holocene Dead Sea level fluctuations presented for past 10 kyr is based on the recognition and dating of shoreline deposits in the fan delta outcrops [88, 111, 117] or facies of drilled sediments [92, 107, 118, 119]. The Dead Sea fluctuated within the range of 370 to 430 mbsl during the past 10 kyr [92]. For the last centuries, a precise reconstruction of the Dead Sea level has been available [90].

During the past 4000 yr., the levels of the Dead Sea fluctuated within the range of 390 to 415 mbsl [111, 120]. Most of the time, the lake level was below 402 mbsl, and the Dead Sea was confined to its deep northern basin [111].

The drop from the late nineteenth century high stand of 390 mbsl to the present (2023) elevation of 428 m bmsl (www.water.gov.il) resulted mainly due to artificial water diversion of runoff water for the basin, superimposed on the climatic trend [111]. The current rate of lake-level decline is significantly high. However, drops of similar magnitude are not unusual in the late Holocene record, which yields evidence of long periods of droughts in the region [111].

3. Poroelasticity and its role in induced seismicity

Water-level-induced seismicity was addressed in several studies of water reservoirs [41, 121, 122] and was also analyzed using the theory of linear poroelasticity [123, 124]. For the first time, a relationship of this kind was considered in the pull-apart basin of tectonic depressions (Dead Sea lake) for earthquake intervals of timescale 100–200 years in Belferman et al. [43, 44]. The increase in water level contributes to an immediate increase in vertical stress due to the weight of the water mass (the loading effect) and subsequently to the increase in horizontal stress due to elastic coupling controlled by Poisson’s ratio [44]. Resulting stress changes are then superimposed on the background tectonic stress field [44]. Further, generation or absence of post-impounding seismicity depends on several additional factors, presented in Figure 3, and discussed in Belferman et al. [44].

Figure 3.

Factors that have the potential to influence the seismic activity prompted by alterations in water levels. Please refer to the text for more information.

Based on the real-time observation of 30 cases of post-impounding seismicity it was indicated that the complex interaction of the water-induced stress with the state of pre-existing stress near the reservoir, together with geological and hydrologic conditions at the site, determines whether these stress changes are sufficient to generate the seismic activity [125]. It was indicated that in areas of high strain accumulation and high levels of natural seismicity, stress changes induced by the reservoir may be small compared to natural variations in contrast to regions of moderate strain accumulation (low to moderate natural seismicity). Thus, in the low to moderate natural seismicity regions, the post-impounding can usually be accompanied by earthquakes. In aseismic areas with slow strain accumulation, the reservoir-induced stresses may be insufficient to raise the stress level to a state of failure [125].

Seismic activity induced by surface water-level fluctuations is also affected by the faulting regime (Figure 3), determined in turn by the respective orientations of the three principal stresses [126]. In regions where the vertical compressive stress is not minimal (normal and strike-slip faulting), seismic activity is more sensitive to the change in the effective stress due to water-level change than in regions where it is (thrust faulting) [123, 125, 127, 128, 129].

In many cases (for instance, the Koyna, Aswan, and Oroville reservoirs), the hypocenter of earthquakes, which are associated with changes in water levels, does not occur directly under the reservoir. In these cases, the faults passing through the reservoir serve as a conduit for pore pressure diffusion [41, 121, 130]. This allows the influence of the reservoir to extend to greater distances and depths, where larger magnitude earthquakes take place.

In addition to the conditions mentioned above, based on observations that some reservoirs continue to be active, whereas others show seismic quiescence under similar water-level fluctuations [37], it was indicated in Belferman et al. [44] that seismicity can strongly depend on the hydromechanical properties of faults (Figure 3).

4. Methods

4.1 Modeling setup

We consider a vertical fault underneath the lake (or reservoir) bed, extending along the outplane axis y (Figure 4), embedded in 2D (plain strain) geometry of the upper crust. We assume that fault failure will result from a change in the effective stress normal to the fault [123, 124, 131, 132, 133, 134]. A step change in the water load, ${\sigma }_w,$ (Figure 4a) is applied on the entire top boundary of the half-infinite structure modeled, representing a valley filled by the lake, where uniaxial strain conditions apply.

Figure 4.

This figure, derived from the work of Belferman et al. [44], illustrates the model used in this study. (a) A step-function water load, applied on a top boundary of the structure, representing a poro-elastic seabed; (b) schematic of the structure below the lake with a vertical fault oriented along the outplane axis y. the calculations are based on uniaxial strain conditions.

4.2 The impact of Lake-level fluctuation on pore pressure

As a result of the water load on the seabed during water-level (WL) change, considering the diffusion, the pore pressure ($p$) is a function of time ($t$) and sub-bottom depth ($z$). Belferman et al. [44] consider a step function load on the half-space, for which, following Roeloffs [123], the analytical solution is:

$$p\left(z,t\right)=\gamma p_s+\left(1-\gamma \right)p_s\mathit{erfc}\left(z/\sqrt{4\mathit{ct}}\right)$$E1

where $p_s$ is the pore pressure on the top boundary of the half-space, $\gamma$ is the loading efficiency, and $c$ is the hydraulic diffusivity. For the post-diffusion stage $\mathit{erfc}\left(\frac{z}{\sqrt{4\mathit{ct}}}\right)\to 1$ at $t\to \infty$ (where $t$ is the diffusion time). For any earthquake hypocentral depth, z, the assumption of the post-diffusion stage is valid, if the diffusion time scale, indicating the travel time of excess pore pressure signal from seabed to hypocentral depth, $t_{\mathit{\text{diff}}}=z^2/c$, is small compared to the earthquake’s recurrence interval (RI) ($t_{\mathit{\text{diff}}}\ll t_{\mathit{RI}}$) [44]. For an average earthquake hypocentral depth underneath the DSF, $z\cong 20\mathit{km}$ [135, 136, 137] and $c=4m^2/s$ justified for the faulted rock, the diffusion time scale indicating a travel time of excess pore pressure from lake bed to the hypocentral depth $t_{\mathit{\text{diff}}}\cong 3\mathit{yr}$. It is negligible compared to the RI of moderate-to-large earthquakes (M > 5.5) in this area, estimated as $t_{\mathit{RI}}=1600\mathit{yr}$ for Lake Lisan and as $t_{\mathit{RI}}=100\mathit{yr}$ for the Dead Sea, based on seismite analyses from the Peratzim Creek and Ein Gedi cores [7, 8, 71, 74].

The water-level variation in Lake Cahuilla has a different tendency compared to the Dead Sea. Subsequently, Hill et al. [45] consider the periodic loading on the poroelastic half-space for which the analytical solution following Roeloffs [123] is:

$$\stackrel{\sim }{p}\left(z\right)=\gamma p_s+\left(1-\gamma \right)p_s{e}^{-z\sqrt{\omega /2c}}{e}^{-\mathit{iz}\sqrt{\omega /2c}}$$E2

were $\omega$ is the angular frequency of water loading ($\omega =\frac{2\pi }{T}$, where T is period time). Similar to a step function load for the post-diffusion stage, for periodic loading, we obtain the long-term limit: $e^{-z\sqrt{\frac{\omega }{2c}}}{e}^{-\mathit{iz}\sqrt{\frac{\omega }{2c}}}=e^{-z\sqrt{\frac{\pi }{\mathit{Tc}}}}{e}^{-\mathit{iz}\sqrt{\frac{\pi }{\mathit{Tc}}}}\to 1$ at $T\to \infty$.

Here, T is defined by the periodic change of water level. Hence, for a long period of time, the diffusion effect is negligible because the pore pressure at the depth reached the value at the top boundary of the half-space. This is expressed in Hill et al. [45] in the “memory” effect of pore pressure at depth, whereby subsequent lakes can contribute to higher pore pressure owing to the diffusive time lag of a previous lake superimposing on the next. This memory effect is pronounced only in one case where the time interval (T) between maximal lake stands is about 41 years (Lake F and Lake E in Figure 2, [45]).

Dead Sea lake WL change is distinctly different from the periodic trend that is characteristic of Lake Cahuilla. This may be due to the pulverization of fault rock on the southern San-Andreas Fault [138].

4.3 Analytical formulation

A poroelastic model with brittle yielding was formulated by Belferman et al. [44], utilizing the Mohr-Coulomb failure criterion. This model was used to calculate the changes in earthquake recurrence interval in response to the increase of water level at basins overlying strike-slip faults.

Within the poroelastic component of the model, produced by the WL change, horizontal stress changes normal to the strike-slip fault calculated under a uniaxial (vertical) strain condition (Eq. (10b) in [44]). This applies to a post-diffusion stage, i.e., when pore pressure at hypocentral depth equilibrates with the lake’s bed ($P_f={\sigma }_w)$. In the case discussed, it is appropriate to skip the diffusion stage since the diffusion time is much shorter than the interseismic period. In this case, the effective vertical stress change becomes zero, while the effective horizontal stress (${\mathit{\Delta \sigma }}^{\prime }$) change becomes:

$$\Delta {\sigma }^{\prime }=\frac{1-2\nu }{1-\nu }\left(\beta -1\right){\sigma }_w$$E3

where $\beta$ is Biot’s coefficient and $\nu$ is Poisson’s ratio, ${\sigma }_w=\mathit{\text{ρgΔh}}$, where $\rho$ is density of water and $g$ is the acceleration of gravity.

Further, under the general loading conditions, the critical value of shear stress (${\tau }_c$) responsible for yielding at the fault is connected to the effective stress normal to the fault (${\sigma }_n^{\prime }$), according to the Mohr-Coulomb failure criterion (e.g., [139]):

$${\tau }_c=C+\mathit{tan}\left(\phi \right){\sigma }_n^{\prime }$$E4

where ${\sigma }_n^{\prime }={\sigma }_n-P_f$ is effective normal stress affected by pore pressure ($P_f)$, $C$ is cohesion, and $\phi$ is the angle of internal friction at the fault.

In order to calculate the change of recurrence interval (∆t) induced by the water level change, we use the Coulomb failure envelope and Mohr circle [140] with the following initial characteristics:

$${\sigma }_0=\frac{{\sigma }_1+{\sigma }_3}{2},{\tau }_0=0,R_0=\frac{{\sigma }_1-{\sigma }_3}{2}$$E5

where ${\sigma }_1$ and ${\sigma }_3$ is the minimal and maximal principal stress, which is horizontal for the strike-slip faults. In this case, if during the interseismic period, the water level increased by $\mathit{\Delta h}$, the horizontal stress change can be specified by ${\mathit{\Delta \sigma }}^{\prime }=\Delta {\sigma }_{\mathit{xx}}^{\prime }=\Delta {\sigma }_{\mathit{yy}}^{\prime }$, following the definitions of principal stresses in 2D [140]. This change manifested in the shift of the Mohr circle toward the origin by ${\mathit{\Delta \sigma }}^{\prime }$ while for the pre-seismic stress state, the Mohr circle center is: ${\sigma }_0^{\prime }={\sigma }_0+{\mathit{\Delta \sigma }}^{\prime }$ (where ${\mathit{\Delta \sigma }}^{\prime }$ is negative at the post-diffusion stage, see Figure 5b, [44]).

Figure 5.

The Dead Sea WL reconstructions for the last two millennia. The dashed curves are suggested by the literature sources. Turquoise anchor points follow Bookman et al. [111] used in WL interpretation, while one point (in dark blue) is shifted to left in error interval of ±45 yr. solid, black line water curve is a plausible scenario suggested by this study.

We assume that the far-field tectonic strain rate, ${{\epsilon }^{\dot{}}}_{\mathit{xy}}$, stays constant over any interseismic time period, $\mathit{\Delta t}$, [44] while the elastic strain, $\mathit{\Delta \epsilon }$, accumulated during this period can be presented by: $\mathit{\Delta \epsilon }={{\epsilon }^{\dot{}}}_{\mathit{xy}}\mathit{\Delta t}$. Based on Hooke’s law and considering only horizontal shear strain, shear stress accumulation at the fault during $\mathit{\Delta t}$ is:

$$\Delta {\tau }_{\mathit{xy}}=2G{{\epsilon }^{\dot{}}}_{\mathit{xy}}\mathit{\Delta t}=\frac{\mathit{\text{Ccos}}\left(\phi \right)}{t_{\mathit{RI}}}\mathit{\Delta t}$$E6

where $t_{\mathit{RI}}$ is a tectonic recurrence interval and C ≠ 0 at the healed fault (e.g., [141]). Hence, a permanent tectonic strike-slip stress rate can be defined as $\frac{\mathit{\text{Ccos}}\left(\phi \right)}{t_{\mathit{RI}}}$ [44]. The change of the shear stress during the interseismic period increases the radius of the Mohr circle while for the pre-seismic stress state, the radius is: $R{\prime }_0=R_0+\Delta {\tau }_{\mathit{xy}}$.

Byerlee’s law envelope [142] was used to define the strength of a seismogenic zone at the fault immediately after the earthquake. An initial stress state, ${\sigma }_{\mathit{\text{init}}}$, is defined as a Mohr circle with a radius $R_0$ centered at the point $\left({\sigma }_0,0\right)$ and restricted by the Byerlee’s law envelope. The larger hypocentral depth is associated with a larger ${\sigma }_0$. Based on the results of Byerlee [142] and on additional laboratory experiments (for a review, see [143]), we adopt, for simplicity’s sake, the friction angle $\phi =0.54$ and initial cohesion $C=0$.

The pre-seismic Coulomb failure envelope is defined by the nonzero cohesion coefficient, $C\ne 0$, specific for the healed fault zone (e.g., [141] and references therein) and friction angle $\phi =0.54$. When the circle reaches the failure envelope, the rock fails at the fault oriented most favorably for sliding (in our case, it is the pre-defined strike-slip fault). The stress will then drop again to ${\sigma }_{\mathit{\text{init}}}$.

After stress release, the time to the next earthquake, $\mathit{\Delta t}$, is calculated from the solution of the Mohr-Coulomb failure criterion for a strike-slip tectonic regime and a water level change, $\mathit{\Delta h}$, characteristic to the Dead Sea lake [44]:

$$\left\{\begin{array}{c}{\left(\tau -{\tau }_0\right)}^2+(\sigma -{\left({\sigma }_0+∆\sigma \prime \right)}^2={\left(R_0+\frac{\mathit{\text{Ccos}}\left(\phi \right)}{t_{\mathit{RI}}}\mathit{\Delta t}\right)}^2\\ \tau =C+\mathit{tan}\left(\phi \right)\sigma \end{array}\right.$$E7

Solving this problem for single solution, we get $\mathit{\Delta t}$ us a linear function of the water level change $\mathit{\Delta h}$:

$$∆t=\left(C+\mathit{tan}\left(\phi \right)\mathit{\text{ρgΔh}}\right)\frac{t_{\mathit{RI}}}{C}$$E8

4.4 Best fit random method of WL curve prediction

In the interpretation of water level over the past two millennia [92, 111, 144], there are several uncertainty. Specifically, the water-level dating method (Radiocarbon dating) could have an error of about ±45 yr., as estimated from the radiocarbon dating of shoreline deposits in fan delta outcrops [111]. The entire past bi-millennial Dead Sea level record is constrained by less than twenty “anchor points” (the data obtained by the dating collected from surveyed paleo-shorelines, [111]). However, historical water level records are quite precise elevation-wise, as they are obtained from different points around the lake [92, 111].

Consequently, the interpretations of the curves (Figure 5) are not identical and are not unambiguous, except for some limitations [43]. Accounting for these limitations, which include the “anchor points” (the data obtained by the dating collected from surveyed paleo-shorelines, [111]) of water level, the 10 million WL curves were generated [44] for the last bi-millennial interval, using a uniformly distributed random number generator. For the simulation, we are setting a ten-year time step.

The linear correlation between RI of widely recorded medium to large historical earthquakes (M > 5.5) available from the literature and WL interpolations was tested on the basis of an estimate of the value of the Pearson product-moment correlation coefficient, R. These statistics were used to assess the suitability of each randomly interpolated WL curve and for identifying out-of-correlation earthquakes.

4.5 Earthquake simulation

The analytical model presented in this section has been discretized with the time step of 1 year and applied to the sequence of WL change samples [44]. From the starting point (AD 33, see Table 1 which can be found in the Appendix of the reference [44]), the simulator moves forward with time, along the WL curve. For each step, the WL change, $∆h_i$, and the amount of accumulated tectonic stress are calculated. After each stress release, the time to the next earthquake, $∆t$, is calculated from the solution of the Mohr–Coulomb failure criterion for a strike-slip tectonic regime. Eq. (9) takes the following discrete form:

$$\left\{\begin{array}{c}{\left({\tau }_i-{\tau }_0\right)}^2+(\sigma -{\left({\sigma }_0+\frac{1-2\nu }{1-\nu }\left(\beta -1\right)\mathit{\rho g\Delta }{h}_i\right)}^2={\left(R_0+\frac{\mathit{\text{Ccos}}\left(\phi \right)}{t_{\mathit{RI}}}\mathit{\Delta t}\right)}^2\\ ={\tau }_i=C+\mathit{tan}\left(\phi \right){\sigma }_i\end{array}\right.$$E9

For each time step, the algorithm determines whether there is a single solution, or two, or nil. A case of no solutions means that the Mohr circle is yet to reach the failure envelope, as the accumulating tectonic stress and the WL increase are still insufficient. At this case, the $\mathit{\Delta t}$ increase with the time steps. The system of Eq. (7) may have a single solution when the failure criterion is met at the end of some timestep or two solutions when it is met before the end of the timestep. A case of two solutions is rounded down to a case of a single solution if a time step (1 year) is small compared to the earthquake RI (about hundreds of years), where the Mohr circle reaches the failure envelope and the earthquake accrues. In this case, the calculated $\mathit{\Delta t}$ saved and resets to zero. At the final stage, we get an array of RI.

4.6 The backward earthquake simulation approach

This formulation is based on the same basic assumptions that were proposed in Belferman et al. [43, 44], as well as set out in the “Analytical formulation of the direct model” section. Similar to the previous simulation, we assume that after an earthquake, the Mohr’s circle returns to the initial state of stress. But in contrast to the previous case, our starting point of reference is the last earthquake, as it is much better constrained. Considering that at the moment before the specific earthquake occurred, Mohr’s circle reaches the failure envelope, the set of Eq. (7) gets one solution. The horizontal effective stress change (${\mathit{\Delta \sigma }}^{\prime }$) calculated from Eq. (3). While ${\sigma }_w=\mathit{\text{ρgΔh}}$, calculated for the WL change $\mathit{\Delta h}$ at the moment of this specific earthquake.

The recurrence interval, $\mathit{\Delta t}$, calculated from Eq. (8), and based on this value, we calculate the time of the previous earthquake $t_p=t_c-\mathit{\Delta t}$. Running this algorithm on the WL change during the period 2–10 kyr (based on the level history of [92]), combined with the WL change received from the best fit random method [44] including the anchor points [111] for the past two millennia, we get a new catalog of simulated earthquakes presented in Figure 6.

This formulation, like the previous one, is based on the solution of a system of equations (Eq. (7)) in the domain determined by positive definite values of $\mathit{\Delta t}$. This is only possible when the WL change ($\mathit{\Delta h}$) does not exceed the value $\frac{C}{\mathit{tan}\left(\phi \right)\mathit{\rho g}}$. When the WL change exceeds the defined value, the model jumps to the next point where the WL change is within the definition domain.

5. Results and discussion

Even though historical water-level records exhibit a significant degree of elevation accuracy, given that they are derived from diverse surveyed locations around the lake [92, 111], the dating of water levels might carry an approximate error of ±45 years and even more, as deduced from radiocarbon dating of sediment deposits in fan delta outcrops [111] or sediment cores [92]. The substantial differences in potential interpolations arise from the uncertainties in dating and resolution aspects of WL reconstructions. This was shown in Belferman et al. [43], using best fit random method WL curve generation. Ten million different curves were generated, constrained with anchor points established by field surveys and radiocarbon dating [111].

Dating historical earthquakes can exhibit a high level of precision, and the validation of accuracy occurs when various historical references align in agreement [61, 62, 63]. But in many cases, the epicenter of even a historical earthquake can be imprecise or not known.

This chapter extends our simulation to prehistoric time (10 ka), for which further challenges include identifying and dating shoreline and earthquake markers [71, 92]. The water-level curve is combined from the one indicated in Belferman et al. [43] (for 2 ka) with that provided by Migowski et al. [92] (for 2–10 ka). Our model for that level curve generates 50 earthquakes.

The recurrence intervals are simulated using the reverse earthquake simulation algorithm, while the starting point of reference is a simulated earthquake that occurred in 1907 CE. We start from this point because it is the last earthquake received from the earthquake simulation presented in Belferman et al. [43], which correlates with the documented earthquake from the field of 1927 CE. The recurrence intervals are depicted in Figure 6 as black points plotted against the water level (orange curve). The calculated years of simulated earthquakes are presented along the time axis by blue stars (Figure 6).

Figure 6.

Orange curve represents the best fit random water level [43] combined with Migowski et al. [92] vs. simulated and historic RIs, correspondingly. The blue dots mark the dates of the seismic events, while the black dots indicate the recurrence interval between these events.

Two types of spatial segments have been identified (gray bars in Figure 5). Light gray bar (∼3490–4020, 4370–4730, and 9380–10,000 years) indicates the time segments where the WL change exceeds the definition domain specified in our model. The dark gray bars indicate the time segments of earthquake sequence (4024–4065, 4724–5068, and 9141–9196) received for periods of height WL stand, with recurrence interval smaller than 35 yrs. The domain excludes the time segments where the WL change almost reaches the critical value of $\frac{C}{\mathit{tan}\left(\phi \right)\mathit{\rho g}}$ (marked in dark gray in Figure 6).

Simulated earthquakes are presented vs. earthquakes from the literature in Figure 7. The blue points (Figure 7) indicate the simulated events supported by literature. Orange points (Figure 7) indicate simulated events, for which there is no evidence in the literature. A cluster of earthquakes was generated in these segments. In Figure 7, we left only the first and the last earthquake of the cluster.

Figure 7.

A comparison dates of historic vs. simulated earthquakes based on the suggested best-fit WL curve [43] combined with WL curve presented by Migowski et al. [71].

In the segments marked in light gray, our model did not perform the calculation. The earthquakes received in these segments were calculated in the previous step. Therefore, our model may not have generated earthquakes that occurred in these segments, such as the one indicated by Ben-Menahem [85] to 1560 B.C. or the seismite estimated by Migowski et al. [71] to have occurred around 7700 B.C.

Of the 50 simulated earthquakes during a period of 10,000 years, only 15 events were not confirmed in the literature (Figure 7). All these 15 events occurred before 3000 BCE. In other words, for all simulated earthquakes by our model within a span of 5000 years, there is evidence in the literature of a significant earthquake occurring in the Dead Sea area. Five earthquakes not confirmed by literature generated at 3068–3533 BCE. In Ben-Menachem’s 1991 study, the archeological excavations at Tleilat el-Rassul of Kenyon were referenced, indicating the possibility of destruction caused by earthquakes that occurred between 3300 and 3600 BCE. In addition, Migowski et al. [71] dated the possible earthquake to 3300–3367 BCE from Ein Gedi core. This evidence can satisfy the simulated events 339 3 BCE, 3533 BCE. Another prehistoric cluster (four events) of earthquakes was simulated in the period 7196–7410 BCE due to a relatively high water level. This period is not well represented by the archeological record due to scarcity of masonry structures.

The advantage of using the backward simulation method, presented for the first time in this chapter, is that it starts from the last earthquake along the Dead Sea fault, for which ample information is available. By contrast, the previous simulation [43] starts from the last earthquake of the studied period. Therefore, it is suitable for a period of 2000 years where earthquakes have historical evidence. The backward method enables us to simulate earthquakes that occurred in prehistoric times, including during the Lake Lisan period. For this period, the water level was significantly higher, and the data regarding both phenomena are subject to even greater uncertainty. Nevertheless, the paleoseismicity studies [8, 12] indicate that despite such high water levels, the seismic activity may have been subdued compared with activity for the Holocene Dead Sea.

However, it is worth noting that our model is limited by the maximum change in water level, and the water level during the Lisan period exceeds this limit. On this occasion, in future studies, it is worth integrating additional physical properties of Earth’s crust. Field data from Scandinavia and North America suggest that the melting of Pleistocene ice sheets was connected to paleoseismic activity [19, 21, 22, 23, 24, 145, 146]. Conversely, the low level of seismicity in Greenland and Antarctica may be due to the current ice load [25], which in turn induces changes in the differential stress in the Earth crust [26, 28].

Hampel et al. [20], based on numerical models that include one or several fault planes embedded in a rheologically layered lithosphere, indicate how loading and flexure of the lithosphere by ice or water decreases the slip rates of nearby faults, whereas unloading and rebound of the lithosphere accelerates the slip (see also [47]). Slip rate variations are caused by changes in the differential stress, which is the difference between the maximum (${\sigma }_1$) and minimum (${\sigma }_3$) principal stresses. When a fault is under a load, the differential stress decreases during loading and flexure of the lithosphere and increases during unloading and rebound. This leads to a decrease and subsequent increase in the fault slip rate.

This configuration needs to be considered when developing the earthquake simulation model for the Lisan period. Additional considerations for further studies are the mechanical condition of the recorder (namely the sediment bed) and the effect of water column height on the dynamics of seismite generation [147, 148].

6. Conclusions

The correlation model, developed in [44] and supported by the comparability of a simulated earthquake with a historical catalog for 2,000 years [43] imitates reaction of pore-elasto-plastic crust. This model is validated by the results obtained in this chapter for 10 millennia of geological records of earthquakes and shows a good correlation for the past five millennia.

Based on the findings presented in the Results-Discussion chapter, we indicate that flexure of the lithosphere, suggested by Hampel et al. [20], may not be significant over two or five millennia, but it should be considered when analyzing longer time periods as Lake Lisan.

Considering this in our future model, we expect to receive an explanation of the low seismicity pattern during the Lisan high stand and solve the model limitations where the WL change exceeds the definition domain, specified in our model.

Nevertheless, the new back-counting model presented in this work enables us to extend the relationship between WL changes and seismicity over a period span of ten millennia.

Acknowledgments

This project was supported by grants from the Ministry of Energy (grant no. 213-17-002) and the German–Israeli Foundation for Scientific Research and Development (GIF; grant no. I-1280-301.8). We would like to express our gratitude to Dr. Regina Katsman and Prof. Zvi Ben-Avraham for their valuable contributions in initiating this chapter.