1 Introduction

Many archaeological sites in the Levant exhibit traces of earthquake damage. The affected buildings are commonly used to acquire information about the location, the time, and the local intensity of past earthquakes (e.g., Stiros and Jones 1996). However, retrieving information from archaeoseismic observations about the coseismic kinematics remains challenging and this is the goal of our research. It was not by chance that the term “archaeoseismological” was used for the first time in connection with earthquakes in the Levant. The Palestine Exploration Fund, founded in 1865, initiated surveys for the exploration of the Levant. In “The Survey of Western Palestine” from their expedition 1871–1878, Condor and Kitchner (1882) describe several archaeological sites in detail, of which numerous show damage now associated with seismic activity (Karcz and Kafri 1978). The earliest archaeological evidence for earthquake damage is reported by Alfonsi et al. (2012), who identify two Neolithic (7500–6000 BCE) earthquakes at Tell es-Sultan in the ancient city of Jericho, within the Dead Sea Fault (DSF) zone.

Information about palaeoseismic sources was obtained directly from archaeological structures that were bisected by fault ruptures. Examples from the DSF include Qasir el Telah (Niemi et al. 2001), Misyaf (Meghraoui et al. 2003), Tell Siçantarla (Altunel et al. 2009), and Tell Ateret (Marco et al. 1997; Ellenblum et al. 1998, 2015). The research at these sites focused on the earthquake history and estimates of the offset of the ruptured structures. In order to improve our knowledge of the properties of pre-instrumental earthquakes in this study, we further analyze the offset walls at Tell Ateret, located in the Jordan Gorge north of the Sea of Galilee (Fig. 1). Tell Ateret provides an opportunity to test whether, in addition to the measurement of the lateral displacement, it is also possible to deduce dynamic source parameters such as (a) the slip velocity of the fault during the earthquakes and (b) estimates of which side of the fault moved during an earthquake. Damaged structures at the site indicate the total relative slip at the site, however, which side moved and by which amount is not evident a priori. And (c) discriminate coseismic slip from possible pre- or post-seismic creep that contributed to the total displacement. Because the magnitude of the causing earthquake is mainly dependent on the offset and length of the ruptured fault segment, any creep contribution to the offset would result in overestimated magnitude. Therefore, determining the extent of potential creep would improve the knowledge base of the fault behaviour and thus contribute to future seismic hazard analysis.

Fig. 1
figure 1

a Digital terrain map of the border between the Arabian Plate and the Sinai Sub-Plate. The black arrows indicate the left-lateral plate movement along the Dead Sea Transform Fault. The blue lines indicate the fault strands which ruptured during the 1202 and 1759 earthquakes. The area within the black rectangle at Tell Ateret is shown enlarged in b. The black line show the main strands of the DSF; the numbers indicate the average slip rates (Gomez et al. 2007; Nemer and Meghraui 2020)

Full size image

This publication uses data collected within the framework of the PhD thesis “The Discrete Element Method in Archaeoseismological Research—Two Case Studies in Israel” of Gregor Schweppe published in 2019 at the faculty of mathematics and natural science of the University of Cologne.

2 Study area

The archaeological ruins of the Crusader fortress on Tell Ateret are 13 km north of the Sea of Galilee (SG), on the main trace of the DSF. The DSF is a left-lateral boundary fault between the Arabian plate and the Sinai subplate (Fig. 1). The entire fault system, which forms the tectonic pull-apart basins of Gulf of Aqaba (GoA), the Dead Sea, and the Hula Basin (HB), transfers the opening of the Red Sea to the Bitlis-Zagros collision zone (Quennell 1956; Freund 1965). The section between the GoA and the SG is mainly characterized by a single strike-slip fault, but in the Hula Basin, the fault splays into several non-parallel strands, the most prominent branches from east to west are the Rachaya Fault (RaF), Yammouneh Fault (YF), and the Roum Fault (RF) (see caption Fig. 1).

The DSF has been active since the Miocene (Bartov et al. 1980). The current seismicity is characterized by infrequent large earthquakes with periods of small to moderate earthquakes in between (Hamiel et al. 2009). The instrumental observation of the earthquakes of the Levant began in 1898 with the installation of the first seismological station (Helwan HLW) south of Cairo, Egypt. Since then the number of seismic stations and networks in the area increased. The catalogs of the larger networks in the region have been merged to a single catalog, which is continuously updated and can be accessed via http://seis.gii.co.il/ (last accessed 02.2020).

The pre-instrumental earthquakes in the area have been the subject of interdisciplinary research. The wealth of information from historical, geological, and archaeological archives were used to compile historical earthquake catalogs (Guidoboni and Comastri 2005; Sbeinati et al. 2005; Ambraseys 2009). Location and fault strands activated during the 1202 and 1759 earthquakes which affected Tell Ateret (Fig. 1) have been discussed in several papers (Marco et al. 1997, 2005; Ellenblum et al. 1998; Daëron et al. 2004, 2005, 2007; Nemer et al. 2008).

3 History of the Tell Ateret

The archaeological strata of the Tell Ateret include (early to late) an Iron Age II fortification, a Hellenistic complex with houses and fortifications, a medieval Crusader castle called Vadum Iacob (meaning Jacob’s Ford), and a Mamluk mosque. Vadum Iacob constitutes the most conspicuous remnant and includes a formidable 4.4 m wide defense wall and inner structures. The name (Latin: Vadum Iacob; Hebrew: Metzad Ateret; Arabic: Qasr al-'Atra) commemorates the legendary abode of the biblical Patriarch Jacob (Gen. 32:10) during the disappearance of his son Joseph (Ellenblum 2003). The most recent structure on the Tell is a Mamluk and Ottoman pilgrimage site with a mosque, which smoking paraphernalia determines its age of late Ottoman period.

During the Crusader period, after the victory of the battle of Montgisard in 1177, the Christian forces wanted to consolidate their power in the Muslim territory  north of the SG. The purpose of the construction of the Crusader Fortress at Ateret in 1178 was to take control of the only crossing of the Jordan River between its source in the north and the SG. The crossing was part of an important trading route from Damascus to Acre (Ellenblum 2007). In 1179, the construction site was conquered by Saladin’s Muslim forces (Barber 1998; Ellenblum 1998, 2007). The unfinished fort was consequently abandoned with the exception of a small mosque constructed within the ruin of the fortress in the Ottoman period and also occasional squatter settlements.

Table 1 in the Appendix summarizes various intensities and magnitudes of past earthquakes of the area. Two of these earthquakes offset the crusader fortress, first on 20 May 1202, and again on 30 October 1759 (Marco et al. 1997; Ellenblum et al. 1998). Ambraseys and Melville (1988) estimated the magnitude of the 1202 earthquake at MS 7.6 with intensities from VIII to IX over a wide area, ranging from Nablus, south of the SG, to Arqa, in current day Lebanon. For the 30 October 1759 earthquake, Ambraseys and Barazangi (1989) estimated MS 6.6 based on the intensity and the radius of the isoseismal of the area in which the earthquake was felt. From the historical records, Sbeinati et al. (2005) derived intensities up to VIII. The earthquake history of the fortification walls is summarized in Fig. 2.

Fig. 2
figure 2

Schematic illustration of the stages of accumulated slip (values are rounded) in the Ateret structures since 1179, timeline from bottom to top (after Ellenblum et al. 2015). The black arrow points north. Left (south) of crusader castle are older man-made structures that have also been affected by past earthquakes (pre 1179). The black horizontal line indicates the course of the DSTF

Full size image

Palaeoseismological research revealed the surface rupture of the 1202 earthquake along the YF in Lebanon (Daëron et al. 2005), Ateret (Ellenblum et al. 1998), and Bet Saida (Marco et al. 2005), with a probable fault length of ~200 km. The October 1759 rupture was found preserved in the Ottoman mosque at Ateret and in palaeoseismic trenches at Bet Saida (Marco et al. 2005), and the November 1759 rupture was located on the Rachaya Fault (RaF), which has a potential rupture length of about ~50 km (Nemer et al. 2008).

4 Methods

We use a numerical model of the fortress wall to back calculate what type of ground movement caused the current deformation. This approach requires an accurate documentation of the current state of the deformation of the northern fortification walls and a numerical model representing the state of the ruin before it was deformed. Based on estimates of the local tectonic setting and the historical record, different scenarios have been developed to represent diverse potential ground movements. Finally, the results of the simulation are compared with the observed deformation.

4.1 Data acquisition

We documented the site of Tell Ateret with two 3D laser scanning surveys. In an initial survey in 2013, we focused on the deformed northern fortification wall (Fig. 3). We combined 17 individual 3D laser scans to a single-point cloud that provides a detailed image of the current state of the deformation as it has stood since its excavation in 2007 (in the following referenced to as point cloud, PCL). In the second survey, the entire site and its surroundings were extensively surveyed with 3D laser-scans (Hinzen et al. 2017), resulting in a detailed digital terrain model with 5 cm grid resolution (in the following referenced to as DTM). For the development of a numerical model in which the deformation of the fortification wall can be simulated, a thorough reconstruction of the original orientation of the fortification wall is vital. Both the PCL and the DTM provide the foundation for the development of the model of the northern fortification wall. In addition, the PCL was used for comparing the results of the numerical simulations with the in-situ situation.

Fig. 3
figure 3

a Digital terrain model (DTM) of the ruin of the Crusader fortress at Tel Ateret and its surroundings. The white rectangle shows the section of the northern wall of the fortress. The grid resolution of the DTM is 5 cm (after Hinzen et al. 2017). b Clear view from a bird’s eye perspective from inside the fortress on the laser scan point cloud (PCL) of the ruptured section of the northern wall which 4.4 m wide; scale varies with perspective

Full size image

4.2 Model development

Model construction and all simulations were carried out with 3DEC, a discrete element (DE) software for geotechnical analysis of discontinuous media such as jointed rocks and masonry, developed by Itasca Consulting Group Inc. (itascacg.com, last visit 11.2019). All elements are represented by blocks. The blocks can be subdivided into smaller blocks allowing designation of different material properties. The 3DEC code has been used in several previous archaeoseismological applications and its physical correctness verified (e.g. Psycharis et al. 2003; Alexandris et al. 2004; Oliveira et al. 2012; Lemos et al. 2015).

The geometry of the model is based on the reconstruction of the laser-scan data. The original angle of some corners in the fortification wall could be measured on well-worked ashlars which span the corners and which are still intact (Fig. 4). The ashlars of the inner and outer shell have a width of 0.5 m. The overall width of the fortification wall is 4.4 m. In the enlarged view to the DTM, the current course of the fortification walls is visible (Fig. 4). Two assumptions were made for the reconstruction: (1) the width of the fortification walls is constant within the range of the model and (2) the angles of the inner and outer shell at common corners of the wall are equal. The second assumption is necessary since some parts of the outer shell have not yet been excavated and are still covered by loose material from construction ramps. The latter were still in place as the fortification was under ongoing construction when it was sacked by Saladin on 30. August 1179 (Ellenblum et al. 1998).

Fig. 4
figure 4

Reconstruction of the original wall course. The orange and green lines mark the current and the reconstructed course of the inner and outer shell, respectively. The green marker on the north-eastern corner indicates a 1.75 m slip between the current wall line and the assumed original position. The upper photo shows the offset of the inner shell. The red highlighted block was fitted by the stonemason to extend over the corner. The blue rectangle marks the section of the fortification wall which is part of the DE model. The red circle highlights the offset of the Ottoman Mosque, which is also shown in the lower photo (photo taken in 1995 by S. Marco). The black line indicates the course of the DSF; the arrows indicate the relative direction of movement (DTM by Hinzen et al. 2017) (after Schweppe 2019)

Full size image

Figure 5 illustrates a top view of the final model, which was used as initial state for all numerical simulations. The model is composed of (a) two shells with 0.5 m wide rectangular blocks; (b) a 3.4 m wide filling between the two shells, realized by a convolute of tetrahedral blocks (edge length of 0.5 m) representing the basalt cobbles of the filling; (c) construction ramps confining the fortification walls, also realized by tetrahedral blocks (edge length of 1 m); and (d) a baseplate composed of two separate blocks to enable horizontal movement, incorporating two levels and a slope with a height difference of 1.5 m to represent the natural topography of the site. In total, the final model was composed of 52,864 rigid blocks.

Fig. 5
figure 5

Top view to the discrete element model of Ateret walls. The ramps inside and outside the fortress are shown in sandy-yellow. Blocks of the inner and outer shells are shown in light grey and the basalt-mortar filling in orange, respectively. The eastern (Arabian Plate) and western (Sinai Plate) baseplates are shown in dark and light grey, respectively. The black arrow points north (after Schweppe 2019)

Full size image

Appropriate material properties must be estimated for all structural elements of the model. As the archaeological remains could not be directly sampled for laboratory tests, published data for similar materials were collected. The limestone of the well-worked ashlars is similar to limestone of the Hazera Formation, material properties of which were published by Kamai and Hatzor (2008). The material parameters of the mortar were taken from Krausz (2002). The density for the basalt is in range of the values estimated by Bourbie et al. (1987). The material properties of the ramps were those of densified soil following Bowles (1996). Table 2 in the Appendix summarizes the material properties of all components.

Table 1 List of reliable historical earthquakes (Zohar et al. 2017) from ca. 760 BCE until 1927 CE that damaged at least one locality in Israel and its close surroundings (see Fig. 1). Date—time of occurrence in year and if known—also the month, day and hour; reported damaged localities—localities reported to have been damaged within the research area (Fig. 1) that we consider of moderate (MR) or higher degree of reliability. Asterisk denotes earthquakes that caused damage also beyond our area of interest; estimated magnitude in previous studies—list of studies and the magnitudes they estimated for that earthquake.; Avg. mag.—the average value of the estimated magnitudes (Zohar et al. 2017); size deg.—In terms suggested by Ambraseys and Jackson, (1998). N-S extent (km)—the distance between the northernmost and southernmost damaged localities, in km and Casu.—estimated number of casualties according to the historical reports: “-”—no casualties or not mentioned or not known, F—few (10 or less), M—many (N10).
Full size table

The properties for the wall joints, normal-stiffness (kn), and shear-stiffness (ks) were estimated following Bui et al. (2017) and Lourenco et al. (2005).

$$ {k}_n=\frac{1}{\left(\frac{1}{E_m}-\frac{1}{E_b}\right)\ast {h}_b}=\frac{E_b\ast {E}_m}{\left({E}_b-{E}_m\right)\ast {h}_b} $$
(1)

and

$$ {k}_s=\frac{k_n}{2\ast \left(1+v\right)} $$
(2)

where Em and Eb are the Young’s moduli of the mortar and blocks respectively, v is the Poisson ratio of the wall blocks, and hb is the block height, which is 0.67 m for the outer and 0.54 m for the inner shell ashlars. In order to create a stable starting model for the dynamic simulations, the complete model was set to static load conditions for a model time of 5 s to allow the model to settle; the state at the end of the 5 s period was then the starting model for all calculations.

The time history for the ground motion consisted of an analytical cycloidal pulse (Zhang and Makris, 2001) and was used in all simulations. The pulse fulfils the boundary conditions of zero acceleration, velocity, and displacement at the start and zero acceleration and velocity at the end of the signal, while the final displacements match those between the two plates. This signal type was previously successfully applied in archaeoseismological research; the displacement history of the pulse is similar to the horizontal ground movement close to a strike slip event (Hinzen 2009; Hinzen and Montabert 2017).

In the high-resolution PCL of the northern fortification wall, a lateral offset of 1.75 m has been estimated (Fig. 6); this value was used to reconstruct the original position of the wall shells.

Fig. 6
figure 6

Laser scanning points are shown in grey. The blue solid lines mark the current position of the ruptured wall (outer perimeter of the ashlars). The blue dashed lines extrapolate the current wall position which due to the deformation are no longer perfectly parallel; the dashed green lines show the inferred original position (Schweppe 2019)

Full size image

4.3 Simulations

We ran a total of 58 simulations for two scenarios. The first scenario focused on the influence of different movement directions and slip velocities of the plates on the resulting deformation pattern. In this scenario, only one movement with a total offset of 1.75 m was considered. The left lateral movement, which is characteristic for the DSF can be the result of different absolute movements of the eastern and western plates. These displacements are not a priori obvious, and therefore, one of the modelling tasks was to decide which movements were the most probable. Figure 7b–e illustrate the potential absolute movements that result in a net left lateral movement. The four movement possibilities were tested in the simulations of the first scenario.

Fig. 7
figure 7

Perspective view of a digital terrain model of Tel Ateret (1 m grid spacing) indicating the plate movements considered in the first scenario. a Location of Tel Ateret on the west bank of the Jordan River. Scale varies with perspective; the fortress measures 140 × 50 m; the legend on the bottom left indicates the elevation above sea level. b The black line shows the trend of the DSF cutting through the fortress. The red arrow (not to scale) indicates the northward movement of the Arabian plate; the Sinai plate is locked. c The eastern and western sides of the DSF move in opposite directions with the same amount of slip. d The Sinai plate moves south, while the Arabian plate is locked. e Both plates move north with the larger slip of the Arabian plate causing a net left lateral motion

Full size image

A movement scenario in which both plates are displaced southwards (opposite of Fig. 7e) was not considered because it contradicts the tectonic regime at the DSF. The slip velocities of the first scenario range from 0.1 to 5 m/s. The slow movement is of particular interest in order to estimate the effect of a potential creep movement to the fortification wall. If a significant fraction of the observed wall offset was due to creep, this would reduce the magnitude estimate of the earthquake. The first scenario consisted of a total of 22 simulations.

In the second scenario, the effect of two consecutive movements deforming the northern fortification wall was simulated, to test the current hypothesis that two earthquakes dislocated the fortress walls. According to Ellenblum et al. (2015), the younger (second) event displaced the site by 0.5 m. To achieve a total offset of 1.75 m, the preceding earthquake must have displaced the fortification walls by 1.25 m. Therefore, the simulations of the second scenario were separated into two stages. The first stage offsets the baseplates 1.25 m and in the second stage, the resulting positions of the first stage are further offset by 0.5 m. The first stage included six simulations with different offset velocities. In the second stage, each end status of the first stage simulation was displaced by 0.5 m, with five different offset velocities, resulting in a total of 30 simulations for this stage.

4.4 Comparison of the simulations with measured data

The following procedure was used to objectively compare the results of each test with the in situ observed deformation pattern as shown in the PCL of the northern fortification wall. Due to the lack of common fixed points between the PCL and the numerical model, the best fit between both datasets has been estimated with the Iterative Closest Point (ICP) technique as described by Rusinkiewicz and Levoy (2001). This technique estimates the best fit between the observed and calculated results by reducing the root-mean-square (RMS) of the residual distances in an iterative process. In 3DEC the coordinates of the starting and end position of the vertices of each block were recorded. From the simulation results, a dataset of points representing the outside of the ashlars of the shells was exported and compared with the measured values.

5 Results

5.1 First scenario

The first scenario focused on testing the effect of the direction of plate motions and displacement velocity on the resulting deformation. In nine simulations, the eastern plate moved 1.75 m north (Fig. 7b) with slip velocities of 0.1 m/s, 0.25 m/s, 0.5 m/s, 1 m/s, 2 m/s, 2.5 m/s, 3 m/s, 4 m/s, and 5 m/s.

In the simulation with slip velocities above 2 m/s, blocks which were not covered by the construction ramps dropped out of the wall due to inertial forces. The inertial forces occur during the starting and stopping phases of movement, and with increasing slip velocities, inertial forces contribute progressively to the deformation of the structure. Schweppe et al. (2017) described two failure mechanisms for free-standing walls utilizing cycloidal pulses related to the inertial forces based on the relation between the direction of ground motion due to plate movement and the collapse direction. The so-called first failure mechanism describes the effect when the direction of collapse is in the direction of the ground motion that occurs during the stopping phase. In the second failure mechanism, the direction of collapse is opposite to the direction of ground motion and occurs during the starting phase. The failure effects were present in the calculation results when the displacement velocities were 2 m/s and higher. However, due to the construction of the fortification wall, the blocks of the inner and outer shell did not show the same falling characteristic. For the first failure mechanism, the blocks of one shell were hampered by the filling between the two shells; the blocks of the other shell can drop out. In other words, in a northward movement, blocks fall out of the northern shell, but not from the southern shell. The opposite effect was observed for the second failure mechanism. This result also occurred during simulations with high velocities such that blocks are hampered from falling by the filling and then still drop out during the stopping phase of the movement.

Inertial forces increase proportional to the acceleration of ground motion, the effect of which can be seen in the first nine simulations where higher slip velocities cause stronger inertial forces and thus larger deformation. In the four simulations with velocities of 1 m/s and lower, no effect of the inertial forces on the deformation pattern is observed. To quantify the effect of inertial forces, it is necessary to analyze the displacements of the blocks relative to the movement of the baseplate. A null or small relative displacement indicates that the block only shifted with the baseplate, and inertial forces did not add additional displacement. Conversely, a large vector of relative displacement indicates a considerable influence of inertial forces. This result only applies to the blocks which rest on a moving plate and thus might help differentiate which plate(s) moved during the earthquake(s).

Figure 8 summarizes the vectors of relative displacement between blocks and the baseplate for the low slip velocities. It is evident that the blocks to the right (east) of the fault line show only small vectors of relative displacement. The largest of these occur in the southern shell close to the fault line, where the shell is pulled apart. The largest overall displacement occurs left of the fault line on the locked plate. During the movement of the western plate, the fortification wall is pulled apart and the blocks on the western plate are dragged along with the movement of the eastern plate as long as the friction between interlocking ashlars is large enough. It is important to note that this pulling does not affect the outer shell of the wall section based on the moving plate because its movement is hampered by filling which is backing the shell. Therefore, only the blocks of the inner shell show relative displacement to the baseplate that does not result from inertial forces, but can be attributed to the pulling. The relative displacement vectors of the four simulations shown in Fig. 8 are similar, a strong indicator that no significant inertial forces contributed to the deformation pattern.

Fig. 8
figure 8

Top view to the model of the fortification wall. Vectors of the relative displacement between the baseplate and the ashlars of the shells at the end of the four simulations with low slip velocities are plotted on top of each other; slip velocities are indicated by the vector colour as shown in the legend. The grey line marks the course of the fault, and the grey arrow points in the direction of the ground motion. All scales are in metres, and the origin (0/0) is the centre of the model (m) (Schweppe 2019)

Full size image

In simulations with slip velocities above 2 m/s, both first and second failure mechanisms occurred. Figure 9 shows the absolute displacement vectors (including the displacement of the baseplate). The total deformation in the simulation with 1 m/s slip velocity is the result of the plate movement and a pulling effect (Fig. 8). At a slip velocity of 3 m/s, numerous blocks drop out of the outer shell in the range of the moving plate in the direction of ground motion (first failure mechanism) and out of the inner shell in opposite direction of the ground motion (second failure mechanism). This is the effect of inertial forces during the stop and start phase of movement, respectively. In the simulation with the highest slip velocity, the fortification wall tilts in the direction opposite to the ground motion due to strong inertial forces. In this case, the fortification wall would suffer total destruction if there were no construction ramps buffering the blocks. The strength of the pulling is also influenced by slip velocity. At low slip velocities, the pulling effect reached blocks at larger distance from the fault line than at higher slip velocities (Fig. 9), with the latter causing a “sharper cut” of the wall.

Fig. 9
figure 9

Top view to the model of the fortification wall. Vectors show the absolute displacement of the ashlars at the end of the simulations with slip velocities of 1 m/s, 3 m/s, and 5 m/s (plotted on top of each other); slip velocities are indicated by the vector colour as shown in the legend. The grey line marks the course of the fault, and the grey arrow points in direction of the ground motion. All scales are in metres, and the origin (0/0) is the centre of the model (m) (Schweppe 2019)

Full size image

In the other simulations of the first scenario, further movement directions were examined. In six simulations, both plates were moved in opposite directions (Fig. 7c) with velocities of 0.5 m/s, 1 m/s, 2 m/s, 3 m/s, 4 m/s, and 5 m/s, and subsequently, the western plate moved southwards (Fig. 7d) with 1 m/s, 2 m/s, 3 m/s, 4 m/s, and 5 m/s. In further two simulations of this scenario, both plates were moved northward (Fig. 7e) with different displacements and slip velocities.

The quality of the fit between simulation results and the in-situ deformation was quantified with the ICP method in the form of the residual distance between both datasets (Rusinkiewicz and Levoy 2001). The median of the residual distances is a measure of how well the simulations fit the observations; the mean deviation of the median (MD) is a measure for the spread of the residual distances. Where both the median and the MD values are small, the calculations are judged to agree well with the observations (cf. Table 3). Figure 10 shows the MD values in relation to the median of the RMS values for all simulations of the first scenario. It is evident that both the displacement velocity and movement direction strongly influence the deformation in the simulation and thus change the quality of the fit. This result positively clarifies that deformation pattern is measurably influenced by slip velocity. Those cases where the eastern plate was locked and the western plate moved in a southern direction show the largest misfit. The simulations where the western plate is locked and the eastern plate moves north show the best fit, and the simulation with a slip velocity of 3 m/s constitutes the preferred scenario.

Fig. 10
figure 10

Mean deviation from the median (MD) plotted versus the median of the residual distances between simulation results and in situ observations for all simulations of the first scenario. Results with the same direction of ground motion are outlined by polygonal splines with the same colour. The labels show the velocity (m/s) of the according simulation; (EP = Eastern Plate and WP = Western Plate). The arrow points to the solution where the EP moves north at 3 m/s with a small MD and median considered the preferred scenario (after Schweppe 2019). Detailed values are given in Table 3 in the Appendix

Full size image

5.2 Second scenario

Based on the results of the first scenario, in the second scenario simulations, the western plate was kept locked and velocities above 1 m/s were considered. In this way, the current hypothesis assuming two earthquakes shifted the fortress wall was tested. The total displacement of 1.75 m was separated into two movements. In the first stage, the eastern plate was shifted 1.25 m north with 1 m/s, 2 m/s, 2.5 m/s, 3 m/s, 4 m/s, and 5 m/s slip velocity in six different simulations. With the exception of the somewhat smaller dislocations, the results of these tests were similar to the results of the first scenario with the same slip velocities. In the second stage, each end state of the simulations of the first stage was offset by another 0.5 m with slip velocities of 1 m/s, 2 m/s, 3 m/s, 4 m/s. to 5 m/s. The median and MD of these calculations are listed in Table 4.

Figure 11 compares the results of the second scenario in the same manner as in Fig. 10. The different combinations of the two movements have a discernible effect on deformation. The simulation with the best fit to the observations has a slip velocity of 3 m/s for the displacement of 1.25 m in the first stage and a slip velocity of 1 m/s for the second stage offset of 0.5 m. Comparing the best fit of results from both scenarios shows an overall better fit for the second scenario and strongly supports the two-earthquake hypothesis.

Fig. 11
figure 11

Mean deviation from the median (MD) plotted versus the median of the residual distance of all simulations of the second scenario. Simulations with the same first stage are highlighted with the same colour. The labels show the velocity (m/s) of the according simulation. The arrow points at the preferred solution with a small median and low MD with slip velocities of 3 and 1 m/s for the first and second earthquakes, respectively (after Schweppe 2019). Detailed values are given in Table 4 in the Appendix

Full size image

6 Discussion

We examine the deformation of the northern fortification wall of the ruin of the crusader fortress of Tell Ateret with DE models. The historical record and the geological observations support the hypothesis that two earthquakes, the first on 20 May 1202 and the second on 30 October 1759, successively deformed the fortification walls of the fortress (Ellenblum et al. 2015). We show that slow creep movement is unlikely to have caused the observed damage because such motion would not have produced the observed deformations by inertial forces. Based on the results of our simulations and a quantitative comparison to observations from the laser scan reconstruction, the most likely scenario is that two distinct plate movements, where the first offset is 1.25 m at a slip velocity 3 m/s and the second offset 0.5 m at 1 m/s.

Intensities and magnitudes for both the two historical earthquakes have been previously derived from the historical information. The results of our simulations can also be compared to the macroseismic studies. Ambraseys and Melville (1988) estimated an intensity of IX for the earthquake in 1202. Based on the relationship between MM intensity and PGA in California after Wald et al. (1999), an intensity IX would correspond to a PGA between 6.38 and 12.16 m/s2. The cycloidal pulse for the first movement with an offset of 1.25 m and a slip velocity of 3 m/s had a PGA of 11.3 m/s2 and lies within the range given by Wald et al. (1999). The second cycloidal pulse with an offset of 0.5 m and a slip velocity of 1 m/s had a PGA of 3.14 m/s2. This translates to an intensity VII for the 1759 event, which agrees with the results given by Sbeinati et al. (2005).

To estimate the magnitude of the simulated events, a relationship for earthquakes in the Middle East from Ambraseys and Melville (1982) can be used:

$$ {M}_s=1.1+0.4\log \left({L}^{1.58}{D}^2\right) $$
(3)

where L is the rupture length and D is the relative displacement in cm. Daëron et al. (2005) suggest the potential rupture lengths of 200 km and 50 km for the 1202 and 1759 events, respectively. With 1.25 m offset, the 1202 event corresponds to a magnitude of MS 7.4 and the second event in 1759 with an offset of 0.5 m results in a magnitude of MS 6.7. Both magnitude values lie in the range of published estimate for these earthquakes (Ambraseys and Melville 1988; Ambraseys and Barazangi 1989; Marco et al. 1997).

Previous archaeoseismic studies at the southern end of Tel Ateret (Ellenblum et al. 2015) revealed repeated rupturing of constructions since the post–Iron Age II era and the authors estimated an average slip rate at the site of 2.6 mm/year. This slip rate is significantly less than estimated from GPS surveys with a lower margin of 4.9 ± 1.4 mm/year (Le Beon et al. 2008; Sadeh et al. 2012; Masson et al. 2015; Wechsler et al. 2018) and the longer-term geologic slip rates of ~5 mm/year (e.g. Marco and Klinger 2014) and 7 mm/year reported for farther north in Syria (Meghraoui et al. 2003). Even though these disparate estimates indicate a certain unevenness in the spatial and temporal distribution of slip along the DSF, the archaeologically determined value of 2.6 mm/year can be compared to the results of this study. Assuming that the whole stress on the fault segment at Tel Ateret was released during the 1202 earthquake, the 818 years since the event would accumulate to a 2.1 m deficit in displacement. If we subtract the 0.5 m slip of the 1759 earthquake, a deficit of about 1.6 m remains. According to Ambraseys and Melville (1988), the deficit is equivalent to a gap of a magnitude MS 7.5 earthquake in this part of the Jordan valley (assuming the geodetic slip rate of 5 mm/year the magnitude would be MS 7.8.).

The uncertainties in our approach using a discrete element model are inherent given that the model is only an approximation to the real geometry of the fortification wall and cannot be an exact copy. The material parameters of the model are based on literature values which are, at best, in the range of the used building materials. Further, only a simple ground motion model (cyclic pulse) was used. In all scenarios, the full displacement was assigned to the main earthquake; the assumption that aftershocks also contributed to displacement were not considered. However, we show that inertial forces contribute to the deformation. Therefore, splitting the offset to a main shock and aftershocks would decrease these inertial forces. In addition, the precise slip markers provided by the archaeological walls was considered sufficient justification for using one horizontal movement for each event and neglecting the possible effect of vertical components.

Both terrestrial laser scans that were made at the site (DTM and PCL) provide a possibility for comparison with future surveys. The PCL was made in 2013 and with the estimated minimum slip rate at the site of 2.6 mm/year, a re-scanning after a period of some 10 years should be sufficient to detect creep effects by comparing the two point clouds. And in the event of an earthquake, there is, of course, the opportunity to compare pre- and post-event deformation. A precondition for future comparisons is that the site is sufficiently protected from further deterioration by weather and/or looting.

7 Conclusions

This archaeoseismological study presents the first time that a DEM has been used to deduce dynamic source parameters of an earthquake from the observed deformation of an ancient structure. The DEM results provided a quantifiable differentiation in the deformation pattern of the wall of Tel Ateret based on the variation of key earthquake source parameters: direction of motion, displacement, and slip velocity. Based on the modelling results, post-seismic creep has been excluded as the main reason for the offset. In addition, best fit estimates between model deformations and laser scans of the contemporary ruins support the hypothesis that two rapid movements of the Arabian plate are more likely to have caused the observed deformations than a single movement of the Arabian plate or a major movement of the Sinai Plate. For the best-fitting simulation, the earthquake of May 20th 1202 was associated with an offset of 1.25 m at a slip velocity of 3 m/s. The consequent earthquake, of 30 October 1759, resulted in an offset of 0.5 m at a slip velocity of 1 m/s. In context of the historical record, these values agree well with the previous investigations of these earthquakes. The 1.25 m and 0.5 m slip, with no significant creep contribution, indicates a slip deficit at the Ateret of about 1.6 m and corresponds to an earthquake of magnitude 7.4.

The number of simulations that are possible such as in this study has practical limits (at this time) due to computation times. Particularly, some of the simulations with slow slip velocities took more than a week of computation time on an average performing workstation. Future advances in software code, advanced macro models of brick and mortar systems, and increased hardware capabilities will allow an increased the number of simulations and thus allow more parameter variations, paving the way for probabilistic modelling in quantitative archaeoseismological studies.