Calibrating a New Attenuation Curve for the Dead Sea Region Using Surface Wave Dispersion Surveys in Sites Damaged by the 1927 Jericho Earthquake

. Instrumental strong motion data is not common around the Dead Sea region. Therefore, calibrating a new attenuation equation is a considerable challenge. However, the Holy Land has a remarkable historical archive, attesting to numerous regional and local earthquakes. Combining the historical record with new seismic measurements will improve the regional equation. 10 On 11 July 1927, a rupture, in the crust in proximity to the Northern Dead Sea, generated a moderate 6.2M L earthquake. Up to five hundred people were killed, and extensive destruction was recorded, even as far as 150 kilometers from the focus. We consider local near-surface properties, in particular, the shear-wave velocity, as an amplification factor. Where the shear-wave velocity is low, the seismic intensity far from the focus would likely be greater than expected from a standard attenuation curve. In this work, we used the Multichannel Analysis of Surface Waves (MASW) method to estimate seismic 15 wave velocity at anomalous sites in Israel in order to calibrate a new attenuation equation for the Dead Sea region. Our new attenuation equation contains a term which quantifies only lithological effects, while factors such as building quality, foundation depth, topography, earthquake directivity, type of fault etc., remain out of our scope. Nonetheless, about 60% of the measured anomalous sites fit expectations; therefore, this new GMPE is statistically better than old ones. From our local point of view, this is the first time that integration of the 1927 historical data and modern shear-wave velocity 20 profile measurements improves the attenuation equation (sometimes referred to as the attenuation relation) for the Dead Sea region. In the wider context, regions of low-to-moderate seismicity should use macroseismic earthquake data, together with modern measurements, in order to better estimate the peak ground acceleration or the seismic


Introduction
Generating a modern and applicable attenuation equation is one of applied seismologists main interests. Considering the Dead Sea area, for which instrumental strong motion data are not available, this task is particularly challenging. Using the Holy Land's historically rich database, researchers had defined seismic intensities and estimated earthquake locations.
Investigating anomalous sites, with seismic intensities higher or lower than predicted from the basic regional attenuation 5 relation, may lead to a better attenuation equation. The local geological conditions can strongly influence the nature and severity of shaking at a given site. Assessing the local geological conditions by geophysical techniques at these anomalous sites, and adding a logarithmic term to a basic attenuation equation, should improve the equation.
This work focuses on the 1927 event, but it is part of wider research which extends to additional earthquakes. The 1927 event was chosen as it is the only devastating one recorded, albeit teleseismically, during the instrumented period. 10 Our main goal in this research is a tighter constraint on the attenuation equation derived from this event. This should allow us to examine whether this preliminary work coincides with our expectations of site amplification and de-amplification due to the local lithology.

Site Response
Ground motion is controlled by a number of variables, including source characteristics, source distance, propagation 15 directivity, near-surface geology, etc. The elastic properties of near-surface materials, and their effect on seismic wave propagation, are crucial to earthquake and civil engineering, and environmental and earth science studies.
Seismic surface waves are initiated at the moment that the earthquake wave front impinges on the surface. These waves spread out, and the surface shakes as they pass. Surface wave amplitude at the surface is controlled by the mechanical properties of the rocks below. These often consist of low velocity weathered rock over bedrock with much higher velocities. 20 When seismic waves pass from a high-velocity layer to a low-velocity layer, their amplitudes and duration typically increase.
The phenomena of site amplification, as a result of soft sediments overlying hard bedrock, is well known since the early days of seismology (Milne, 1898). Site-effects are also well known and were investigated after several major earthquakes: Mexico City 1985 (Singh et al., 1988), Yerevan 1989 (Borcherdt et al., 1989), San Francisco 1989 (Hough et al., 1990), Los Angeles 1994 (Hall et al., 1994) and Kobe, 1995 (Aguirre andIrikura, 1997). Therefore, local lithology is a crucial factor for 25 estimating site amplification, defined as the amplitude ratio between the surface layer and the underlying bedrock. Site amplification at a specific site can be attributed to many factors, such as basin effects, focusing effects, topography, and reverberation of the seismic waves in the upper layers due to acoustic impedance differences ( Figure 1).
The amplification, A, is proportional to the reciprocal square root of the product of the shear-wave velocity, Vs. (Eq. (1)) (Aki and Richards, 2002): where ρ is the density of the investigated soil. As shear-wave velocity decreases by a given fraction the amplification increases by half that fraction (for a constant density). Since density plays a minor role (Dal Moro, 20152014;Xia et al., 5 1999) the Vs value can be used to represent site conditions.
The most widely used index in the classification of the soil response is the average shear-wave velocity in the uppermost 30 meters, the Vs30. This index was accepted for site classification in the USA -National Earthquake Hazards Reduction Program -NEHRP (Building Seismic Safety Council, 2001). In Europe by the new provisions of Eurocode 8 (BSI, 2011), and in Israel it is accepted by the design provisions for earthquake resistance of structures -SI 413 (The Standards Institution  10 of Israel, 2013). The value of 30 meters comes from the USA and European building codes, where it was found empirically that this depth is directly proportional to deeper and shallower values (Boore et al., 2011). Zaslavsky et al. (2012) argued that the use of Vs30 is a simplification that cannot be justified in the complex geological conditions in Israel, yet no alternatives have thus far been proposed. Therefore, in this scenario, the Israel Standards Institute still adopts the Vs30 index.
In modern attenuation equations, also known as ground motion prediction equations (GMPE), coefficients are derived from 15 strong motion data, namely from ground acceleration measurements. In the past, and in areas lacking the technology to record earthquakes, it was impossible to measure the peak ground acceleration (PGA) directly. Therefore, it is common to categorize historical earthquakes with seismic intensity scales that describe the damage at each site or area (Ambraseys, 2009;Guidoboni and Comastri, 2005)

The M6.2 1927 Jericho earthquake 20
The left-lateral Dead Sea transform separates the Sinai-Levant Block from the Arabian Plate ( Figure 2). The 6. 2ML July 11, 1927, Jericho earthquake (Ben-Menahem et al., 1976Shapira, 1979) was the strongest and most destructive earthquake to hit the Holy Land during that century. Furthermore, it was the first to be instrumentally recorded by seismographs. The epicentral location was originally estimated at a few kilometers south of the Damia Bridge, which is 30 kilometers north of Jericho (International Seismological Summary -ISS Bulletin of 1927). In the following decades new estimates have been 25 published: Shapira et al. (1993) calculated the epicenter to be near Mitzpe Shalem. Zohar and Marco (2012)  The damage from the earthquake was heavy, especially in places near the source, but not only there: In Nablus, located 70 km from the epicenter (Figure 2), 60 people were killed, 474 were injured, and more than 700 structures were destroyed, most of which were built on soft sediments (Blankenhorn, 1927;Willis, 1928). By comparison, Jerusalem is only about 30 kilometers from the source and the damage there was much smaller, especially in property. However, in Mount Scopus and the Mount of Olives (eastern neighborhoods in Jerusalem), the damage exceeded that in other parts of Jerusalem (Abel, 5 1927;Brawer, 1928). Other cities also suffered from this earthquake: Tens of people were injured and even died, and hundreds of houses were ruined in Ramleh and Lod (Brawer, 1928). Jericho in the Jordan Valley also suffered significant damage, especially in terms of buildings collapsing ( Figure 3). The total number of victims was about 350-500 (Ambraseys and Melville, 1988;Amiran, 1951;Arieh, 1967;Ben-Menahem, 1991). Beyond the casualties, several environmental effects were reported: The Jordan river flow ceased near the Damia bridge for about 21.5 hours (Willis, 1928) and a one-meter 10 seiche wave was observed in the Dead Sea (Abel, 1927;Blankenhorn, 1927). Some evidence suggests that the earthquake was felt up to 700 kilometers from the epicenter (Ben-Menahem, 1991), although a different interpretation suggests this distance was only 300 kilometers (Ambraseys and Melville, 1988).
Compiling historical evidence, Avni (1999) where MMI is the Modified Mercalli Intensity (assumed to be equivalent to MSK), M is the magnitude and d is the distance from the epicenter. Raphael and Agnon (2018) note that damage from the 1927 Jericho earthquake was higher east of the transform (on the Arabian Plate) than on the west (Sinai-Levant Block). This observation, consistent with their archaeoseismic findings for 25 earthquakes in antiquity, requires further study.

MASW Theory
The MASW method is environmentally friendly, non-invasive, low-cost, rapid, robust, and provides reliable Vs30 data (Miller et al., 2002). Multichannel records make it possible to separate different wavefields in the frequency and velocity domains. Fundamental and higher modes can be analyzed simultaneously, but generally, only the fundamental mode is used 5 because it has the highest energy (Park et al., 1998).
The MASW method consists of three main steps: (A) Acquisition of experimental data, (B) signal processing to obtain the experimental dispersion curve, and (C) inversion to estimate Vs30 ( Figure 5). The inverse problem consists of estimating a set of parameters that describe the soil deposit, based on an experimental dispersion curve. Inversion problems based on wave propagation theory cannot be solved in a direct way due to their non-linearity. Thus, iterative methods must be used where a 10 theoretical dispersion curve is determined for a given layer model and compared to the previously obtained experimental dispersion curve (Ryden et al., 2004). Vs30 typically does not converge to one stable value. In other words, for the same dispersion curve, one will get slightly different Vs30 depending on the initial model. Yet, the ensuing uncertainty is of little consequence since Vs30 enters logarithmically into the attenuation equation (see below). 15 We carried out the surveys with a linear array of 24 vertical geophones (R.T. Clark's geophones with natural frequency of 4.5 Hz) at equal intervals of 2-3 meters over a total length of 46-69 meters. For the survey sound source we used a fivekilogram sledgehammer striking a twenty-centimeter square aluminum plate at variable offsets of 5, 10, 15, 20, 25 and 30 meters (both forward and reversed). The seismic data were recorded on a Geometrics Geode seismograph at a sampling rate mostly of 8 kHz for 0.5-2 seconds (Table 1). For an acceptable Signal to Noise Ratio, we used the so-called "vertical 20 stacking" approach, which is a summation of multiple synchronized repetitions of the test (usually five times).
Rayleigh wave dispersion curves are obtained by the MASW module of the RadExPro® software, whose calculation procedure is based on a paper by Park et al. (1998), and also by the WinMASW® software. From all the dispersion images that we calculated from each offset shot (Figure 6), we choose the smoothest and clearest one ( Figure 6) to compute the site's Vs30 profile. An inversion process then finds the shear-wave velocity profile whose theoretical dispersion curve is as close as 25 possible to the experimental curve ( Figure 6). The data and coefficients are automatically inverted via genetic algorithms which represent an optimization procedure belonging to the classification of global-search methods. Genetic algorithms are commonly used to generate high-quality solutions to optimization and search problems by relying on bio-inspired operators such as mutation, crossover and selection Compared to traditional linear inversion methods based on gradient methods (Jacobian matrix) these inversion techniques produce a very reliable result in terms of precision and completeness (Moro et 30 al., 2007).

Velocity model
All models were considered to be a stack of homogeneous linear elastic layers, neglecting lateral variations in soil properties.
The number of unknowns for a layered model, when considering only shear-wave velocity, is three for each layer: density, thickness, and one elastic constant. Therefore, the number of unknowns is 3n-1 (where n represents the number of layers).
The change in density with depth is usually small in comparison to the change in shear modulus and is normally neglected 5 (Park et al., 1997).

Number of layers & layer thicknesses
The resolution of surface wave surveys decreases with depth. Thin layers are well resolved when they are close to the surface, whereas at great depth, the resolution is limited and only large changes can be detected (Foti et al., 2014).
Regardless of the number of the layers of the site, Vs30 is almost the same in each case (Figure 7). For those reasons, as well 10 as the lack of density information, we did not restrict each model to a specific number of layers. Without boreholes and lithostratigraphic data, which is the case in our work, a useful rule of thumb is to assume layer thicknesses increasing with depth, to compensate for the decreased resolution with depth, an intrinsic shortcoming of surface wave testing (Foti et al., 2014).

Depth of investigation 15
We used a five-kilogram sledgehammer and summed up five strikes. For some sites, this type of source is insufficient to determinate a shear-wave profile down to 30 meters. In addition, at some sites, we were not able to spread the geophones at intervals of more than two meters, which limited the length of the seismic line. This length probably excludes longer wavelengths which limits the depth of investigation. Lastly, as the shear-wave velocity of the lowest frequency is highermore data is available for deeper layers. Therefore, the penetration depth will decrease in areas with low shear-wave 20 velocity. For instance, if we can clearly detect a phase velocity of about 300 m/sec at 5 Hz, we can roughly estimate a depth of investigation of approximately 20-30 meters according to the following equation: where n equals 2-3 (Foti et al., 2014;Moro, 2015). In other words, this equation emphasizes that the depth of investigation is about a half to a third of the largest wavelength observed. 25 We carried out the surveys with a linear array of 24 vertical geophones (R.T. Clark's geophones with a natural frequency of 4.5 Hz) at equal intervals of 2-3 meters over a total length of 46-69 meters. For the survey sound source we used a fivekilogram sledgehammer striking a twenty-centimeter square aluminum plate at variable offsets of 5, 10, 15, 20, 25 and 30 meters (both forward and reversed) ( Figure 6A). The seismic data were recorded on a Geometrics Geode seismograph at a 5 sampling rate mostly of 8 kHz for 0.5-2 seconds (Table 1). For an acceptable Signal to Noise Ratio, we used the so-called "vertical stacking" approach, which is a summation of multiple synchronized repetitions of the test (usually five times).
Rayleigh wave dispersion curves are obtained by the MASW module of the RadExPro® software, whose calculation procedure is based on a paper by Park et al. (1998), and also by the WinMASW® software. From all the dispersion images that we calculated from each offset shot ( Figure 6B), we choose the smoothest and clearest one ( Figure 6C) to compute the 10 site's Vs30 profile. An inversion process then finds the shear-wave velocity profile whose theoretical dispersion curve is as close as possible to the experimental curve ( Figure 6D). The data and coefficients are automatically inverted via genetic algorithms which represent an optimization procedure belonging to the classification of global-search methods. Genetic algorithms are commonly used to generate high-quality solutions to optimization and search problems by relying on bioinspired operators such as mutation, crossover and selection compared to traditional linear inversion methods based on 15 gradient methods (Jacobian matrix) these inversion techniques produce a very reliable result in terms of precision and completeness (Dal Moro et al., 2007).
From 24 surveys, we succeeded in extracting Vs30 for 19 of the 20 sites studied (the Hartuv data were too noisy for interpretation) ( Table 2 and Supplementary). These would be used to recalibrate the attenuation equation arrived at by previous investigators at 133 sites19 of the 133 sites 20

Velocity model
All ground models were considered to be a stack of horizontal homogeneous elastic layers, neglecting lateral variations in soil properties. The number of unknowns for a layered model, when considering only shear-wave velocity, is three for each layer: density, thickness, and one elastic constant. Therefore, the number of unknowns is 3n-1 (where n represents the number of layers). The change in density with depth is usually small in comparison to the change in shear modulus and is 25 normally neglected (Park et al., 1997).

Number of layers & layer thicknesses
The resolution of surface wave surveys decreases with depth. Thin layers are well resolved when they are close to the surface, whereas at great depth, the resolution is limited and only large changes can be detected (Foti et al., 2014).
Regardless of the number of the layers at the site, Vs30 is almost the same in each case (Figure 7). For these reasons, as well as the lack of density information, we did not restrict each model to a specific number of layers. Without boreholes or other direct lithostratigraphic constraint, which is the case in our work, a useful rule of thumb is to assume layer thicknesses increasing with depth, to compensate for the decreased resolution with depth, an intrinsic shortcoming of surface wave testing (Foti et al., 2014).

Depth of investigation 5
We used a five-kilogram sledgehammer and summed up five strikes. For some sites, this type of source is insufficient to determinate a shear-wave profile down to 30 meters. In addition, at some sites, we were not able to spread the geophones at intervals of more than two meters, which limited the length of the seismic line. This length probably excludes longer wavelengths which limits the depth of investigation. Lastly, as the shear-wave velocity of the lowest frequency is highermore data is available for deeper layers. Therefore, the penetration depth will decrease in areas with low shear-wave 10 velocity. For instance, if we can clearly detect a phase velocity of about 300 m/sec at 5 Hz, we can roughly estimate a depth of investigation of approximately 20-30 meters according to the following equation: where n ranges between 2 and 3 (Foti et al., 2014;Dal Moro, 2014). In other words, this equation emphasizes that the depth of investigation is about a half to a third of the largest wavelength observed. 15 Zohar and Marco (2012) relocated the 1927 epicenter to a point near the Almog settlement. We used this most recently published epicenter to calculate new epicentral distances for the 133 sites. Since Equation 2 above is dependent upon d, we checked the variable scatter in the points, but found that the changes in the best-fit coefficients were very minor, so that we assumed for all purposes to use the original. 20 Figure 8 shows a scatter plot of the original MMI (assumed equivalent to MSK) vs. new d for their 133 sites. Hough and Avni (2011) fit this data with a curve whose equation best describes the attenuation equation for this event. Using the mathemathical form of their curve, we calculated upper and lower limits such that 60% of the points are enclosed. This we call the 60% prediction boundary. We consider that the lithological effects probably account for much of the scatter beyond this boundry, due to amplification and de-amplification. 25 Zohar and Marco (2012) relocated the epicenter to a point near the Almog settlement. We used this most recently published epicenter for calculating the new epicentral distances, d. Figure 8 shows a scatter plot of MMI vs. d for their 133 sites. Hough and Avni (2011) fit this data with a curve which best describes the attenuation relation for this event. Using the mathemathical form of their curve, we calculated upper and lower limits such that 60% of the points are enclosed. This we call the 60% prediction boundary. We consider that the lithological effects probably account for much of the scatter beyond 30 this boundry, due to amplification and de-amplification.

MASW surveys
From 24 surveys, we succeeded in extracting Vs30 for 19 of 20 sites (the Hartuv data were too noisy for interpretation) ( Table   2 and Supplementary).

Discussion
A number of researchers have studied the 1927 event. Avni (1999) tried to reduce the impact of local geology and attempted 5 to generate basic attenuation curves for specific azimuths. Zohar and Marco (2012)  measurements. An attenuation equation with a term that depends on the Vs30 index should lead to a better understanding of past events, and to more useful predictions of future earthquakes.

Survey locations and validation 10
The decision as to where exactly each survey should take place was based on Avni's thesis (Avni, 1999). Where the location was not sufficiently known, we rechecked the reference given by Avni. In most cases, there was evidence of specific damaged buildings. We tried to locate these buildings on historical maps (1927)(1928)(1929)(1930)(1931)(1932)(1933)(1934)(1935)(1936)(1937)(1938)(1939)(1940)(1941)(1942)(1943)(1944)(1945). Unfortunately, most sites were located inside urban areas, where we could not carry out the seismic surveys. Therefore, we surveyed in nearby open areas as close as possible to the referenced damage zones. 15 To validate our results, we compared them with a summary of thousands of seismic evaluations around Israel carried out over the years by the Geophysical Institute of Israel (GII), and compiled in a report by (Aksinenko and Hofstetter, 2012).
These evaluations were based upon refraction and borehole velocity measurements yielding Vs and/or Vp values, as well as the effects of topography and geology. The spacing of their data was such that often a number of GII sites had to be averaged to provide a value within several kilometers for comparison with our MASW values. However, Figure 9 shows that the GII-20 based values are in consistent agreement with those of the MASW. However, this comparison is a bit tricky because Vs30 results for two sites 5 km or much less distant could be significantly different, as shown in Figure 10. Remembering that Vs30 enters a logarithmic term, we find our approach potentially useful. Table 2 lists our 24 sites alphabetically, with their respective Vs30 values, the computed errors, and epicentral distances, d. 25 The Vs30 values vary from a low of 232 m/sec in Beit Alfa, -85 m.s.l. (Figure 11

A new attenuation equation
In the present case of the 1927 earthquake, the sources of the data are mostly historical documents and not strong data 5 measurements. This makes it difficult to quantify site response into a single equation. In the practical modern attenuation relationship, Vs30 is a crucial index. A term that depends on Vs30 has previously been constrained for several large data sets (Abrahamson et al., 2014;Boore et al., 1997;Campbell and Bozorgnia, 2008). We chose the Boore et al. (1997) attenuation equation (Eq. (4)) in order to emphasize site response:.
where Y is the ground-motion variable (peak horizontal acceleration or pseudo-acceleration response in g), M is the moment magnitude, r is the epicentral distance in kilometers, VA, and all b terms are frequency dependent coefficients to be determined. By adding Boore et al.'s (1997) Vs term to Hough and Avni (2011) where VA and C4 are adjustable coefficients. The first four coefficients remain the same as we assert that the magnitude, attenuation, geometrical spreading and site response are all independent. We adopt the value of VA from Boore's (1987) equation (Eq. (4)), as it represents a single value independent of the frequency. We took formerly derived GMPE, with its coefficients, and added another term, by regressing only for the new coefficient, then optimizing C4 and VA by Least Squares Fitting (LSF), as shown in Figure 11

The performance of the new attenuation equation
With these coefficients, 58% or 11 of 19 sites, were amplified or de-amplified as we expected. For the entire distance range (up to 250 km) the Vs30 corrections leave 42% sites out of the prediction boundary (eight of nineteen sites). Seismic intensities at all these eight sites are overpredicted by the attenuation equation (Eq. (2)) ( Figure 12Figure 13). We expect that 25 Vs30 at these sites will be higher than 655 m/sec in order to obtain de-amplification. However, our results show the opposite effect -these eight sites are characterized by lower Vs30 which drives amplification. This can be caused by the fact that measurements were taken over agricultural fields, of which the upper layers (the first few meters) are characterized by low shear-wave velocity, decreasing the average Vs. Another reasonable explanation is that we did not succeed in extracting the average shear-wave velocity down to 30 meters and perhaps we missed some high-velocity shear-wave layers at deeper layers. In such cases, we constrain the last layer to be thicker in order to estimate Vs30 for all our surveys.

Conclusions
In this research, we investigate site amplification and de-amplification around Israel. According to previous studies (Aki, 5 1988;Boore, 2003;Borcherdt, 1994;Field and Jacob, 1995;Joyner and Boore, 1988), the local lithology can amplify or deamplify wave amplitude. The commonly used modern seismic method -MASW -allowed the extraction of Vs profiles at 19 sites reportedly damaged by the 1927 ML6.2 earthquake. We use these profiles to update the attenuation equation for the Dead Sea region by including the Vs30 term.
According to this new equation, 11 sites, which constitute 58% of our measured samples, move into the 60% prediction 10 boundary. This suggests that the prediction boundary actually encompasses over 80% of the macroseismic observations. This fit is better than any available attenuation equation for the Dead Sea region. However, as we have used only 19 sites, we should consider further research and provide wider results. Although our final equation (Eq. (6)) shows amplification and deamplification depending on Vs30, it does not take into consideration any other factor, such as building quality, foundation depth, topography, earthquake directivity, type of fault, etc. Obviously, for better results, we must use additional methods 15 and jointly invert some other seismic data such as: refraction (S and P waves), Horizontal to Vertical Spectral Ratio (HVSR), MASW of the transverse component of Love waves, MASW of the radial component of Rayleigh wave, Extended Spatial Auto-Correlation (ESAC), etc. Also, with these data in hand, a full inversion for the epicenter will be in order.
Despite the scarcity of data, this is the first time that an integration of historical data with shear-wave velocity profile measurements improves the attenuation relation. In order to better estimate the peak ground acceleration or the seismic 20 intensities that will be caused by future earthquakes, attenuation relations are necessary for areas characterized by high seismicity. Some of the regions of low to moderate seismicity have rich sources of historical earthquake data. The integration of historical data with modern shear-wave velocity profile measurements will lead to a better understanding of future earthquakes.