Classifying marine faults for hazard assessment offshore Israel: a new approach based on fault size and vertical displacement

For many countries, the methodology for offshore geohazard mitigation lags far behind the well-established onshore methodology. Particularly complicated is the assessment of fault hazard in the marine environment. The determination of whether a fault is active or not requires ultra-high-resolution seismic surveys and multiple coring and, unfortunately, frequently ends with uncertain results. Moreover, if a pipeline must cross a fault, it is not enough to determine whether the fault is active; slip rates are needed for resistant planning.

Here we suggest a new approach for fault hazard assessment for the master planning of infrastructure. We provide planners a way to choose a route that will cross the least hazardous faults; these faults will then be investigated in site-specific surveys for slip rates that will allow seismic design. Instead of following the onshore practice that is hard to implement in the marine environment, we suggest taking advantage of the marine environment where seismic data are commonly better in quantity and quality. Based on existing industrial 3D seismic surveys, we measure for each fault in the study area the amount of its recent (in our specific case, 350 ka) vertical displacement and the size of its plane. According to these two independently measured quantities, we classify the faults into three hazard levels. This allows planners to choose infrastructure routes that cross the least hazardous faults at an early stage of planning and direct them to sites that need further investigation.

Our case study is the Israeli continental slope, where numerous salt-related, thin-skinned, normal faults dissect the seabed, forming tens of meters high scarps. A particular hazardous zone is the upper slope south of the Dor disturbance, where a series of big listric faults rupture the seabed in an area where the sedimentation rate is 4 times faster than the vertical displacement rate. We suggest that this indicates exceptionally fast creep, seismic rupture, or rapid tremor and slip episodes.

1 Introduction

The need for geohazard assessment in the marine environment is increasing globally due to the growing number of infrastructures laid on the seafloor. To mitigate marine geohazards, numerous studies have been conducted in many world basins (Georgia Basin, Barrie et al., 2005; Sea of Marmara, Armijo et al., 2005; Gulf of Mexico, Prior and Hooper, 1999; offshore California, Clark et al., 1985, and references therein; Norwegian Sea, Shmatkova et al., 2015; Italian continental margins, Chiocci and Ridente, 2011; and more). Some of these studies focus on submarine landslides, and when faults are considered, they are sometimes treated as static seabed obstacles. Note, however, that even extremely accurate mapping of the seafloor does not provide the information needed to determine whether the fault is active or not.

Onshore, the criteria for fault activity is based on displacement of dated stratigraphic markers – 11 kyr (California, Bryant and Hart, 2007), 13 kyr (Israel, Sagy et al., 2013), and even older (see Styron and Pagani, 2020, compilation for examples). To determine if a specific fault is active, the continuation or displacement of the chosen marker is examined in outcrops or trenches.

In the marine environment, such an approach is much more complicated and requires high-resolution seismic surveys and core analyses (Posamentier, 2000; Kvalstad, 2007; Hough et al., 2011). Commonly a site-specific survey includes four steps (Prior and Hooper, 1999; Angell et al., 2003): (a) mapping the seafloor, (b) establishing a chronostratigraphic scheme by tying high-resolution seismic data to dated horizons in boreholes, (c) structural mapping of the fault and displacement measurements, and (d) geological interpretation and quantification. This commonly used approach is problematic because each survey requires months of work and frequently yields uncertain results.

One practical option is to define all faults rupturing the seabed as active faults in the absence of age information (On, 2016). This approach is based on the rationale that faults are likely active if they are identified at the seabed despite sedimentation. Note, however, that fault scarps can remain hundreds of thousands of years on the seabed without any additional jump if the sedimentation rate is lower than the displacement rate.

Another approach is to apply a probabilistic fault displacement hazard assessment (PFDHA, Wong and Stepp, 1998; Youngs et al., 2003; Angell et al., 2003), analogs to the practice developed for earthquake ground motion prediction (PSHA, Cornell, 1968, 1971). Such analysis provides a graph showing the annual frequency calculated for various displacement values. This probabilistic approach requires assumptions regarding (1) creep versus seismic slip and (2) the number of seismic events that had produced an observed displacement. In the case studied here (offshore Israel), we do not know if the studied faults produce earthquakes at all (maybe they only creep). Furthermore, if they produce earthquakes, we do not have any information about the magnitudes–frequency relationships. Thus, it seems that this approach will not yield robust results in our case.

The goal of this study is to provide a practical and relatively fast solution for early-stage planning of marine infrastructure that must cross a faulted zone. For instance, there is no choice in the case studied here, and planning requires a route that will cross the least hazardous faults. For this, we need criteria to determine the relative fault displacement hazard level. We base this determination on the amount of recent vertical displacement and the size of the fault plane. We assume that bigger faults with larger past displacements have a greater potential for larger future ruptures.

Our analysis takes advantage of the wealth of high-resolution seismic data frequently available offshore. Instead of investing in multiple coring to find out whether or not each specific fault in the study area displaces a recently dated horizon (a practice that frequently fails to provide an answer), we measure the vertical displacement of a 350 ka horizon (the most recent-dated horizon that covers all our research area) and the area of the fault plane. These two parameters are tough to measure in seismic data usually available on land and are easily measured in high-resolution surveys frequently available offshore.

2 Scientific background

2.1 Geological history of the Levant Basin

The Levant Basin was formed in the late Paleozoic and early Mesozoic, alongside the opening of the Tethys Ocean that had separated Africa from Eurasia (Garfunkel and Almagor, 1984; Garfunkel, 1988, 1998; Robertson, 1998). At that time, several rifting phases created a system of horsts and grabens spreading from the northern Negev northwestwards into the Levant Basin (Bein and Gvirtzman, 1977; Garfunkel and Almagor, 1984; Garfunkel, 1988, 1998; Robertson, 1998). After the rifting stage, approximately at the end of the Early Jurassic (∼ 180 Ma), the Levant continental margins turned passive and continued to accumulate sediments for more than 100 million years (Gvirtzman and Garfunkel, 1997, 1998; Steinberg et al.,  2008; Bar et al., 2013).

At the end of the Turonian and the beginning of the Santonian (∼ 84 Ma), a change in the relative movement between Africa and Eurasia led to a change in the stress regime and folding along the Syrian Arc began (Krenkel, 1924; Henson, 1951; De-Sitter, 1962; Freund et al., 1975; Reches and Hoexter, 1981; Eyal and Reches, 1983; Sagy et al., 2018).

About 35 million years ago, a large area, including east Africa and northern Arabia, started rising above sea level. This process provided large amounts of clastic sediments to the Levant Basin, where the sedimentation rate increased significantly (Gvirtzman et al., 2008; Steinberg et al., 2011; Avni et al., 2012; Bar et al., 2016, 2013). These clastic sediments compose the Saqiye Group, which thickens from tens of hundreds of meters in the Israeli coasts to 1.5 km in the continental shelf area (Gvirtzman and Buchbinder, 1978) and 6 km in the deep Levant Basin (Steinberg et al., 2011).

About 6 million years ago, the connection between the Mediterranean Sea and the Atlantic Ocean was restricted during a short event termed the Messinian Salinity Crisis (MSC). During the crisis, the sea level dropped, and an evaporite sequence a few kilometers thick accumulated in the entire Mediterranean Sea (e.g., Ryan and Hsü, 1973; Hsü et al., 1973). The salt sequence offshore Israel is nearly 2 km thick in the deepest portion of the basin, thinning landwards and nearly pinching out to zero beneath the continental slope (Ryan and Cita, 1978; Mart and Gai, 1982; Gradmann et al., 2005; Bertoni and Cartwright, 2006; Netzeband et al., 2006; Gvirtzman et al., 2013, 2017).

In the Pliocene, the Nile, one of the largest rivers in the world, supplied a huge amount of sediments to the eastern Mediterranean that buried the Messinian salt and produced a giant delta with a well-developed deep-sea fan (Mascle et al., 2001). Alongshore currents transporting sediments from the Nile Delta through the Sinai coast to the Israeli coast gradually formed the continental shelf offshore Israel (Gvirtzman and Buchbinder, 1978; Goldsmith and Golik, 1980; Carmel et al., 1985; Stanley, 1989; Tibor et al.,  1992; Buchbinder et al., 1993; Golik, 1993, 2002; Buchbinder and Zilberman, 1997; Perlin and Kit, 1999; Ben-Gai et al., 2005; Zviely et al.,  2006, 2007; Klein et al., 2007; Schattner et al., 2015; Schattner and Lazar, 2016; Zucker et al., 2021). The slope of this continental shelf is currently faulted by faults, which are the target of this study.

2.2 Regional tectonic activity

It is generally agreed that the Levant continental margin is passive since the Middle Jurassic (Garfunkel, 1988; Gvirtzman et al., 2008) with no deep-seated active faults south of Mount Carmel (Fig. 1a). Neev et al. (1973) and Neev (1975) debated this consensus, suggesting that an active fault, which they named the Pelusium Line, runs all along the Israeli continental margin. On the other side, Garfunkel and Derin (1984) and Garfunkel (1988), argued that all faults crossing the Pliocene–Quaternary section offshore Israel are thin-skinned and salt-related. Nearly 30 years later, based on better seismic data, Gvirtzman et al. (2008) and Gvirtzman and Steinberg (2012) showed that a continental margin fault zone (CMFZ) does exist approximately at the same location marked by Neev et al. (1973), but this line stopped operating in the Miocene. The deep-seated faults along the CMFZ (∼ Pelusium Line) are sealed with Miocene strata and do not reach the surface.

Figure 1 Location maps. (a) Regional setting – bathymetry and topography from Hall (1994). Abbreviations: EEZ – Israel's exclusive economic zone; CTZ – Carmel–Tirza Fault; MC – Mount Carmel; DST – Dead Sea Transform; MLT – Mount Lebanon Thrust Fault. (b) The studied area with thin-skinned faults in black after Gvirtzman et al. (2015). Red polygons are borders of seismic surveys mentioned in the text.

Moving northwards to offshore northern Israel and Lebanon, the current tectonic activity is different. The Carmel–Tirza Fault (CTF, Fig. 1a), a branch of the Dead Sea Transform (DST), disrupts the continental margin off the Galilee (Kafri and Folkman, 1981; Garfunkel and Almagor, 1984; Ginzburg and Ben-Avraham, 1987; Schattner et al., 2006; Sharon et al., 2020). In addition, the Mount Lebanon Thrust (Fig. 1a) disrupts the continental margin offshore Lebanon and produces earthquakes, such as the 551 CE M = 7.5 Beirut–Tripoli earthquake (Elias et al., 2007).

2.3 Thin-skinned, salt-related normal faulting along the Israeli continental slope

Unlike the deep-seated faults that stopped operating in the Miocene and do not reach the seabed, numerous thin-skinned normal faults rupture the seafloor all along the Israeli continental slope (Fig. 1b), creating steep steps that are tens of meters high (Almagor and Garfunkel, 1979; Garfunkel et al., 1979; Mart and Gai, 1982; Almagor, 1984; Garfunkel, 1984; Garfunkel and Almagor, 1984; Tibor et al., 1992; Gradmann et al., 2005; Martinez et al., 2005; Bertoni and Cartwright, 2005, 2006; Netzeband et al., 2006; Mart and Ryan, 2007; Cartwright and Jackson, 2008; Cartwright et al., 2012; Gvirtzman et al., 2013, 2015; Katz et al., 2015; Safadi et al., 2017; Gadol et al., 2019).

Recently, based on improved bathymetry data, the seabed traces of these faults were mapped in detail (Gvirtzman et al., 2015; Katz et al., 2015; Kanari et al., 2020), showing that their scarps are not buried by sediments. This apparently indicates that displacement rates are higher than burial rates. However, averaged over hundreds of thousands of years, displacement rates are roughly similar to sedimentation rates (Elfassi et al., 2019). This indicates that the fault scarps observed on the present seafloor may have formed by recent instantaneous seismic ruptures (Elfassi et al., 2019) or rapid episodic motions. In any case, these relatively shallow thin-skinned faults are incapable of producing large earthquakes (discussed below) because their fault planes are relatively small compared to crustal faults. The major hazard they pose is surface rupture, which may as well trigger slumps (Katz et al., 2015).

The recognition that the thin-skinned faults along the Levant continental margin are related to salt tectonics has been stated in many studies (Neev et al., 1976; Ben-Avraham, 1978; Almagor and Hall, 1979; Garfunkel et al.,  1979; Mart and Gai, 1982; Garfunkel, 1984; Garfunkel and Almagor, 1984; Tibor et al., 1992; Gradmann et al., 2005; Martinez et al., 2005; Bertoni and Cartwright, 2006, 2007; Loncke et al., 2006; Netzeband et al., 2006; Hübscher and Netzeband, 2017; Mart and Ryan, 2007; Hubscher et al., 2008; Cartwright and Jackson, 2008; Clark and Cartwright, 2009; Cartwright et al., 2012; Gvirtzman et al., 2013; Gadol et al., 2019; Ben Zeev and Gvirtzman, 2020; Hamdani et al., 2021). In particular, it has been suggested that faulting was initiated by basinward salt flow (Gradmann et al., 2005; Bertoni and Cartwright, 2006, 2015; Allen et al., 2016; Cartwright et al., 2018; Kirkham et al., 2019) triggered by basinward tilting of the continental margin, as a result of coastal uplift (Cartwright and Jackson, 2008; Elfassi et al., 2019; Hamdani et al., 2021).

The beginning of faulting was initially dated to a relatively broad time interval between the late Pliocene and the early Pleistocene (e.g., Garfunkel et al., 1979; Almagor, 1984; Gradmann et al., 2005; Netzeband et al., 2006). Later, based on 3D high-resolution seismic surveys, Cartwright and Jackson (2008) showed that offshore central Israel faulting began in the mid-Pliocene. Then, in the late Pliocene, it spread northward, and in the early Pleistocene, southward. Elfassi et al. (2019) established a new chronostratigraphic scheme for the Pliocene–Quaternary section offshore Israel that allows better fault dating. By combining seismic and bio-stratigraphic data, they divided the Pliocene–Quaternary sequence into four units (Fig. 2): unit 1 – Pliocene (5.33–2.6 Ma), unit 2 – Gelasian (2.6–1.8 Ma), unit 3 – Calabrian–Ionian (1.8–0.35 Ma), and unit 4 – Ionian–Holocene (< 0.35 Ma). Based on the improved chronostratigraphy, Elfassi et al. (2019) measured displacement rates on several faults offshore central Israel (in the Sara–Myra survey, Fig. 1b) and concluded that during the Pliocene faulting activity was minor (< 4 m Myr⁻¹); then, in the Gelasian, it peaked to rates of > 100 m Myr⁻¹ (10 cm kyr⁻¹). Later it decreased to rates of ∼ 50 m Myr⁻¹ (5 cm kyr⁻¹).

Figure 2 Uninterpreted (a) and interpreted (b) seismic section across the Levant continental margin offshore Israel (location in Fig. 1b). Chronostratigraphic and seismostratigraphic of the Pliocene–Quaternary section after Elfassi et al. (2019). Green – seafloor; light blue – base unit 4; purple – base unit 3; blue – base unit 2; red – base unit 1 (and top evaporites); yellow – base evaporates. Thin-skinned faults in black lines.

In what follows, we use the chronostratigraphy of Elfassi et al. (2019) to map the most recent horizon (350 ka) in the entire study area (light blue – base unit 4 in Fig. 2b) and identify the zones with the most potent recent activity.

2.4 The Dor and the Palmahim disturbances

Two huge (10s of km) thin-skinned, rootless structures were observed in the 1970s along the Israeli slope – the Palmahim and the Dor disturbances (Fig. 1b; Garfunkel et al., 1979; Garfunkel, 1984; Almagor, 1984). Some studies described these disturbances as deep-seated tectonic structures (e.g., Neev et al., 1976), while others suggested gravitational instability structures induced by basinward sliding over late Miocene (i.e., Messinian evaporites) detachment surface (Garfunkel et al., 1977). According to Garfunkel (1984), these disturbances are similar to other gravitationally induced rootless structures, which are related to the flowage of underlying salt or shale, that are known in deltas and continental margins in other parts of the world (Bruce, 1973; Evamy et al., 1978; Harding and Lowell, 1979; Crans et al., 1980). The Palmahim structure has been described as a rotational slide, bounded between two translational faults (Mart et al., 1978; Garfunkel et al., 1979), while the Dor disturbance seems to be the focus of a regional deformation zone (Garfunkel, 1984; Gadol et al., 2019).

3 Data and methods

This study aims to map seabed fault scarps and their subsurface continuation. Accordingly, we start with bathymetry analysis quantifying fault scarps, and then we use seismic data to map faults in the subsurface and to measure displacements.

Table 1 Seismic data.

3.1 Bathymetry

The Israel national bathymetry survey provides pixel resolution of 15 m until a water depth of ∼ 700 m (Sade et al.,  2006, 2007) and 50 m between isobaths 700 and 1700 m (Tibor et al.,  2013). In addition, we used bathymetric grids with ∼ 10 m cell size, derived from four 3D seismic surveys listed in Table 1 (Aviya; Dalit; Yam Hadera; and Sara–Myra).

To quantify the height of fault scarps at the present seafloor, we developed an algorithm that uses the fault map prepared by Gvirtzman et al. (2015) and automatically calculates elevation differences from both sides of the fault segment every 50 m. A fault segment is a visually mappable lineament in a bathymetric map, regardless of its connection to other segments in the sub-seabed. The algorithm begins with manually moving fault segments, marked by Gvirtzman et al. (2015), to their most accurate location, that is, along the maximal slope of the seabed fault scarp. Then, for each point along each fault, the algorithm measures the dip angle and the true fault direction (dip direction), ignoring the possibility that the fault scarp may have changed by erosion and/or sedimentation. For each point along the faults, the algorithm searches the two closest points from both sides of the fault according to the true dip direction and the dip angle. The calculated output includes the three components of the fault movement: throw, heave, and displacement. This algorithm was applied to all grids described in Table 1, and the measurements were used for the throw analysis.

3.2 Seismic reflection data

The seismic data used here include 2D and 3D industrial seismic reflection surveys processed in the time domain (TWT, two-way time) and 3D seismic cubes that were pre-stack depth migrated (Table 1). All surveys were loaded and interpreted using the Kingdom IHS software (IHS Markit, 2022). Preliminary mapping of the four seismic units described above was done by Elfassi et al. (2019). Ben-Zeev and Gvirtzman (2020) expanded this mapping to cover Israel's exclusive economic zone (EEZ; Fig. 1a). Here, we recheck and remap these horizons in detail along the continental slope where faults are common and map these four units in 3D seismic volumes (Fig. 1b, Table 1).

Subsurface faults mapping

Each fault with a seabed expression was mapped in the subsurface. When a seismic volume was available, we mapped the faults in 3D manually and frequently showed that separate seabed segments connect in the subsurface. If only the 2D seismic lines were available, connectivity between segments would sometimes remain uncertain. In addition, we measured the cumulative displacement of the ∼ 350 ka horizon (the theoretical vertical resolution is a quarter of a wavelength, in our case 4–5 m).

The subsurface mapping of faults adds several layers of information on top of seabed mapping: (1) it allows measuring the displacement of dated horizons and thus indicates the rate of motion; (2) it allows distinguishing between small surface faults that are minor and small surface faults that connect at the subsurface to large faults; (3) it allows identifying hidden fault segments, which do not appear on the bathymetry but may rupture it in the future; and (4) it provides a 3D view of the fault plane, which is essential for structural analysis (as well as estimation of potential earthquake magnitudes, if these faults rupture seismically; see discussion).

4 Results

4.1 Measurements

4.1.1 Scarp height

Figure 3a shows the heights of seabed scarps measured from both sides of all faults every 50 m. The map shows that between the Palmahim and the Dor disturbances, fault scarps are relatively low (<20 m), whereas from the Dor disturbance northwards, they are significantly higher (20–90 m). This observation is consistent with extension measurements that also increases northwards (Cartwright and Jackson, 2008; Ben Zeev and Gvirtzman, 2020).

Figure 3 Vertical offset measurements along faults. (a) Height of seabed scarps derived from bathymetry analysis. (b) Vertical offsets at the base unit 4 horizon measured from seismic data. Assigning 350 ka to the base unit 4 horizon (Elfassi et al., 2019), its vertical offset is transformed into vertical displacement rate (the left-hand side of the scale bar in panel b). (c) Vertical offset measured at the base of unit 4 (gray dots) and scarps height at the seafloor (orange dots). Note that vertical offsets in bathymetry increase northwards, whereas vertical offsets at the base of unit 4 increase in the vicinity of the Dor disturbance. Bathymetry from Tibor et al. (2013).

4.1.2 Vertical displacement rate

The problem with analyzing bathymetry alone is that faults scarps are reduced by sedimentation and erosion and do not correctly represent fault displacement. Therefore, we also measure fault throw along the youngest regionally mappable horizon (base unit 4, Fig. 3b), which yield vertical displacement rates averaged for the past 350 kyr (the best possible representation of “recent” in the study area). These measurements highlight an exceptionally active zone in the vicinity of the Dor disturbance with vertical displacement rates reaching 40–50 cm kyr⁻¹ (Fig. 3b); this anomaly is not detected in the bathymetric analysis (Fig. 3a), emphasizing the need for subsurface measurements. To further illustrate the Dor anomaly, Fig. 3c shows a projection of all seabed and subsurface offset measurements along a south–north section, emphasizing peak throws near the Dor disturbance (at X ∼ 32^∘38^′0^′′ N), nearly 2 times larger than in surrounding areas.

4.1.3 Sedimentation rate

Considering the 350 ka age of base unit 4 (Elfassi et al., 2019), recent sedimentation rates (thickness of unit 4 divided by 350 kyr) can be calculated for the entire study area (Fig. 4c). Results indicate relatively low (< 60 cm kyr⁻¹) values in the deep basin, increasing to ∼ 90 cm kyr⁻¹ in the mid-slope and > 150 cm kyr⁻¹ along the basinward propagating shelf edge (Ben Zeev and Gvirtzman, 2020). Particularly interesting is the off-shelf area near the Dor disturbance reaching > 200 cm kyr⁻¹ (the impact of this observation on fault interpretation is discussed below).

Figure 4 (a) Faults on bathymetry background (after Tibor et al.,  2013). (b) Base unit 4 structure map. (c) Unit 4 sedimentation rate. Half grabens separated the Dor disturbance from the shelf edge and emphasized its dome shape seen in panel (b). These half grabens are filled with a thick section of unit 4 with a sedimentation rate exceeding ∼ 1.8 m kyr⁻¹ (c). A high sedimentation rate is also observed along the shelf edge, expressing shelf progradation during the past 350 kyr.

In addition to the shelf edge belt, a large thickness of unit 4 is observed in a deep half graben separating a prominent dome at the center of the Dor disturbance from the shelf edge (Fig. 4b). The accommodation space created by this half graben is quickly filled by sediments arriving from the nearby shelf edge. South of the Dor disturbance, the half graben is separated from the shelf edge (Fig. 5b). North of the disturbance, the two features create a continuous sedimentary package (Fig. 5a). Noteworthy, the listric faults east of the half graben are different from all the other faults as will be discussed below.

Figure 5 Cross section (a) north and (b) south of the Dor disturbance (seismic location in Fig. 4a). Normal faults in black lines. Seismic reflectors as in Fig. 2.

4.2 Fault classification

4.2.1 Vertical displacement categories

To classify faults according to their vertical displacement, we assign each fault segment a single value of maximum throw measured anywhere along it (a) at the seabed (height of scarp) and (b) at the base of unit 4 (vertical offset). Results are presented in Fig. 6 in three colors – white represents faults producing seabed scarps < 30 m, red represents faults with seabed scarps of 30–60 m, and turquoise represents faults with seabed scarps > 60 m. This illustration is consistent with Fig. 3, showing that seabed fault scarps higher than 30 m (red) are more common near Dor and northwards (Fig. 6a). In contrast, fault scarps higher than 60 m (turquoise) are observed only north of Dor (Fig. 6b) with the exception of one outlier near the Palmahim disturbance. Noteworthy, the threshold values of 30 and 60 m were defined for convenience, such that all three groups will include a reasonable number of faults, and the third group with exceptionally high values will be smaller. If needed, these threshold values can be changed.

Figure 6 Fault classification by vertical throw after assigning each fault segment a single value, representing the maximum vertical displacement measured anywhere along it. (a) Faults forming seabed scarps smaller (white) and higher (red) than 30 m. (b) Faults forming seabed scarps larger than 60 m (turquoise). (c) Faults displacing base unit 4 by more than 60 m (blue) and smaller than 60 m (black). (d) Faults displacing base unit 4 by more than 100 m (yellow). Note that faults with the largest vertical throw are concentrated around the Dor disturbance. Background in all maps is shaded relief of bathymetry (Tibor et al., 2013).

Consistent with our hypothesis that fault scarps are decreased by sedimentation and erosion, classification according to vertical offsets at the base of unit 4 (Fig. 6c and d) portrays a different picture with peak vertical displacements in the vicinity of the Dor disturbance without increasing northward trend (again, one outlier near Palmahim). In particular, we highlight the faults bounding the Dor disturbance from the east (Fig. 6d), where large throws (> 100 m) at the base of unit 4 are observed. These faults coincide with the listric faults mentioned above (Fig. 5b). Uncommonly, these faults form seabed scarps higher than 60 m (Fig. 6b) despite the exceptionally high sedimentation rate observed at that location (Fig. 4c).

4.2.2 Hidden fault segments

To map fault planes in the subsurface and measure their area, we use high-resolution 3D seismic volumes. Figure 7 illustrates that 35 fault segments rupturing the seabed on the eastern side of the Sara–Myra survey converge at depth to seven major faults. Noteworthily, a part of the fault marked by red (Fig. 7b) has no surface expression (Fig. 7a). This hidden fault segment ruptures the three lower horizons (Fig. 7d–f), reaching base unit 4 in several locations (Fig. 7c), and is unseen at the seabed (Fig. 7a). Namely, the partial seabed segments do not represent the actual fault size. A similar analysis conducted for the Yam Hadera seismic survey illustrates that several major fault segments (marked green, pink, and yellow) are hidden (Fig. 8a and b). The sub-seabed fault mapping in 3D highlights the segments “missing” in the bathymetry.

Figure 7 Subsurface mapping of fault planes. (a) 35 faults segments rupturing the seabed in the eastern part of the Sara–Myra survey (location in Fig. 1). (b) A 3D view of fault planes illustrating that the 35 fault segments at the seabed belong to seven major faults (each fault represented by one color). An example of some fault plane area measurements. (c–f) Structural maps of four subsurface horizons (base units 4–1), each with faults crossing it (same colors as in panel b). Note the hidden fault segment (dashed red line in panel a), which does not disrupt the seabed but may rupture it in the future.

Figure 8 (a) Seabed faults in the Yam Hadera seismic survey with hidden fault segments marked in the same colors as the fault to which it is connected in panel (b). (b) All faults are displacing base unit 4. Each fault is represented by one color. (c) Three-dimensional illustration of five faults with their measured fault plane area. Note that the yellow and the pink faults are not detected at the seabed in some parts (hidden fault segments) despite their large plane area (101 km² and 37 km², respectively), and their colors are the same as in panel (b).

Figure 8c presents an example of five 3D-mapped faults with their measured plane area. The red and blue faults are two sides of one graben rooted in the Messinian salt layer with a relatively minor fault plane area, whereas the pink, yellow, and green faults have a significantly larger fault plane area despite their shallower penetration only to the top of the Pliocene (base unit 2) horizon (Fig. 8c).

4.2.3 Geometry and location relative to the salt wedge

Another way for fault classification is based on their geometry and location relative to the underlying salt layer (Fig. 9). Group I produces horsts and grabens (marked blue) mostly along the base of the continental slope, west of the salt wedge-out line. The faults of Group I displace the entire Pliocene–Quaternary section down to the top salt horizon (Fig. 10a, cross section aa^′), and their dip angle varies around 45^∘ (Fig. 11).

Figure 9 Classification of faults according to structure and location. The hanging wall in red, the footwall in green, and the grabens in blue. Group I consists of horsts and grabens, running along the base of the continental slope west of the salt wedge-out boundary (black line). Group II consists of down-stepping normal faults with hanging walls always on the basinward side, mostly located east of the salt wedge-out line. Group III is strike-slip faults. Bathymetry from Tibor et al. (2013).

Figure 10 Seismic cross sections illustrating the difference between Group I (a-a^′, Sara–Mira survey) located west of the salt wedge and Group II (b-b^′, Yam–Hadera survey) located on the continental slope east of the salt wedge (Messinian evaporites are missing). Note that the faults of Group II do not displace the top unit 2 horizon (unit 2 in blue). Location in Fig. 9.

Figure 11 Panel (a) shows the relationship between the average faults dip angle, the fault width, the location of the fault along the slope (water depth), and the offset of base unit 4 horizon. Group II, located on the upper slope, is characterized by larger fault width (1.5–4.5 km), gently dipping (∼ 30^∘) fault planes and larger fault offsets. Panel (b) shows the difference between the offsets as they were measured on the seabed and base unit 4 (cumulative offset). The differences in faults offsets are much larger in Group II.

Group II consists of seaward-dipping faults producing a series of down-stepping stairs (growth faults, rotated blocks, and half grabens) mainly in the upper slope, east of the salt wedge-out line (Fig. 9). These faults are highly listric (Fig. 10b, cross section bb^′) as already described above (Fig. 5). They are characterized by smaller dip angles of about 30^∘ (Fig. 11) and do not displace unit 1 (Fig. 10b, cross section bb^′).

Group III are relatively long strike-slip faults with a few hundred meters of lateral displacement (Ben Zeev and Gvirtzman, 2020). Their vertical throw is relatively small, and its direction changes along the strike (Fig. 9).

The high-resolution seismic volumes Sara–Myra and Yam Hadera allow detailed investigation of the difference between Group I and Group II according to three parameters: displacement (at the seabed and at the base unit 4), fault plane dip, and fault width. The upper panel in Fig. 11 illustrates that the fault width negatively correlates with the dip angle (larger width for gently dipping faults). The lower panel of Fig. 11 further illustrates that displacement (seabed and base unit 4) negatively correlates with water depth (faults in shallower waters have larger displacements).

Group I, located in the deeper waters (blue), is characterized by relatively small vertical offsets (better seen in the base unit 4 horizon), high dip angles (> 45^∘), and relatively short fault width (0.5–2 km). Group II (the listric faults), located in shallower waters (yellow), is characterized by larger (base unit 4) vertical offsets, lower dip angles (∼ 30^∘), and larger faults widths (1.5–4.5 km). These observations highlight the listric faults (Group II), located east of the salt wedge-out line (Fig. 10b), which are big in size and in vertical offsets.

5 Discussion

5.1 Seabed versus subsurface mapping of faults

Detailed mapping of the seafloor has become standard practice in marine geohazard assessment, and the demand for improved resolution is continuously growing. Here we show that bathymetry is not enough for faults investigation, even if it is extremely detailed, because fault scarps are strongly affected by sedimentation and erosion; hence their heights do not represent the real offsets. In fact, the subsurface mapping may be more informative even if its resolution is lower. For example, peak vertical displacements of faults near the Dor disturbance are twice the size of those measured along nearby faults; yet this is not observed on the bathymetry because the scarps are quickly buried. Sedimentation rates averaged on 350 kyr indicate > 200 cm kyr⁻¹ near the Dor disturbance (Fig. 4c). Moreover, a 6 m long core retrieved nearby (location in Fig. 4c) with a sedimentation rate of > 850 cm kyr⁻¹ (Ashkenazi et al., 2022) indicates that the sedimentation rate may have increased in the last couple of thousands of years. Note that the sedimentation rate calculation includes all sources of material accumulated due to the downslope transport of materials.

The drawback of these measurements is their dependency on the quality of the seismic data. Where only 2D lines are available, the measured value represents the throws at the survey–fault intersection, which may represent the tip of the fault; moreover, some faults may not be crossed by any seismic profile.

Additional support for the advantage of subsurface mapping is the structural map of the 350 ka horizon (Fig. 4b) and the calculated sedimentation rate map (Fig. 4c). These maps show that the most active regions in the study area are the half grabens surrounding the Dor disturbance from the east (Fig. 5). These half grabens are rapidly subsiding (thick unit 4), and the faults bounding them are the most active. We suggest that while the faults of Group I are driven by the salt flow that produces extension above it (Hamdani et al., 2021), the faults of Group II are also affected by the gravitational collapse of the continental slope.

5.2 Fault classification

Based on the maximal vertical displacement of base unit 4 (Fig. 6c and d), we classify all fault segments mapped on the seabed (rupturing the seabed) to three vertical offset levels. Vertical offset smaller than 60 m is considered low, 60–100 m is considered moderate, and > 100 m is considered high (Fig. 12a).

Figure 12 Fault classification by vertical displacement (a) and size (b). Each seabed fault segment is assigned a value based on its subsurface structure, i.e., the maximal vertical displacement measured along the fault segment at the base unit 4 horizon and the total area of all segments, connected at the subsurface. When the 3D mapping of a fault is unavailable, fault size is expressed by its length in a map view. (c, d) Vertical displacement at base unit 4 versus fault size (length/area). Pink, yellow, and blue present three levels of vertical displacement and size, which are proxies for surface rupture and potential earthquake magnitudes, respectively. While classification by the two criteria correlates for most faults, black circles mark faults whose vertical displacement is high relative to their size, and red circles mark faults that are big relative to their vertical displacement. Bathymetry from Tibor et al. (2013).

Based on the size (area of fault plane or its length on surface projection), we classify all faults mapped at the subsurface to three levels. Fault planes smaller than 10 km² or shorter than 5 km are considered small, an area of 10–20 km² or length of 5–10 km is considered moderate, and an area larger than 20 km² or length longer than 10 km is considered big (Fig. 12b).

It should be noted that unlike the classification by vertical displacement, which is performed on seabed segments, the classification of faults by size is performed on fault planes, and a single fault plain frequently combines many seabed segments (i.e., the number of fault planes in our database is significantly smaller than the number of seabed segments).

Though the two classification criteria are independently measured, and despite a certain degree of arbitrariness in choosing the cutoff values (60 and 100 m of vertical displacement; 10 and 20 km² for fault plane area), it is interesting to compare the resulting maps. For most faults in the study area, the two criteria yield a similar category (Fig. 12c and d). That is, fault segments with high vertical displacement levels are usually a part of a big fault and vice versa, similar to observations related to deep-seated tectonic faults (Wells and Coppersmith, 1994). Exceptions, marked in Fig. 12c and d by black circles (moderate vertical displacement and small faults), mainly belong to Group II, which is exceptional in many ways, as shown above. Conversely, exceptions marked by red circles (big faults with small vertical displacement) belong to Group III, which are strike-slip faults whose vertical displacement is not expected to correlate with their dimensions.

Finally, we provide a simplified map that combines the two measured parameters into a single hazard level (Fig. 13). In this map, a high level is assigned to a fault segment, which either is characterized by high vertical displacement or large planes; low means low vertical displacement and small plane area; moderate is all the rest. This map simplifies the use of our analysis for early planning of new infrastructures on the seabed, which is the aim of this study.

Figure 13 Final simplified fault displacement hazard map classified into three hazard levels according to a combination of the two criteria presented in Fig. 12 (i.e., vertical fault displacement and size). The combination is conservative; i.e., a high level is assigned to a fault segment, which either is characterized by high vertical displacement or belongs to a big fault; low means low vertical displacement and small size; moderate are all the rest. Bathymetry from Tibor et al. (2013).

5.3 Listric faults south of the Dor disturbance

The listric faults south of the Dor disturbance (part of Group II) are particularly exceptional. Their planes dip gently with lower angles; they have a bigger width but do not penetrate unit 1; they are located on the steep slope, east of the salt wedge; and, particularly important, they produce large seabed scarps despite their location in a high sedimentation zone. In fact, the sedimentation rate at that location is 4 times larger than the vertical displacement rate (∼ 200 vs. ∼ 50 cm kyr⁻¹, respectively, Figs. 3b and 4c). Allegedly, this observation indicates that these faults are creeping faster than the sedimentation rate, they slip seismically, or they operate in rapid episodes of tremor and slip (ETS) or “slow earthquakes” (Ito and Obara, 2006; Ikari et al., 2011). The possibility of seismic rupture was already raised by Elfassi et al. (2019) for the deep basin faults of Group I in the Sara–Myra survey, where sedimentation rates are similar or slightly higher than vertical displacement rates. In that case, continuous creep seems unlikely because its rate is similar to the burial rate and cannot produce significant seabed scarps, unless the scarps are produced by a recent rapid creep event.

5.4 Earthquakes and faults

If the thin-skinned faults offshore Israel are seismically active, they might produce earthquakes and ground shaking in addition to surface rupture. This possibility is apparently supported by the many epicenters located near the faults and particularly around the Dor disturbance (Wetzler and Kurzon, 2016). The problem is that the depths of these earthquakes are much deeper (10–30 km) than the shallow thin-skinned faults (1–2 km). Katz and Hamiel (2019) argued that many hypocenters coincide with the CMFZ at a depth of about 18 km, but this is inconsistent with Gvirtzman and Steinberg (2012), who showed that the CMFZ stopped operating in the Miocene.

The accuracy of hypocenter depths offshore Israel is highly uncertain, as stated by Wetzler and Kurzon (2016), because of the lack of seismic stations at sea and because of the simplified velocity model they extended from the onshore area. Therefore, at this stage, we cannot determine whether the recorded earthquakes offshore Israel are produced by thin-skinned faults or by deeper sources.

Another complication is the relations between fault plane area, which commonly exceeds 10 and even 20 km² (Figs. 8 and 12b), and earthquake magnitudes, which are mostly 2 < M < 4 (Wetzler and Kurzon, 2016). These values are inconsistent with the empiric relations measured in deep-seated faults (Wells and Coppersmith, 1994), where fault planes of 10–20 km² are typically associated with M ∼ 5 earthquakes (which have not been recorded yet). However, deep-seated faults are different from thin-skinned faults in many ways, leaving us with unclearness. This short discussion indicates that the seismicity of the thin-skinned faults needs more research, which is crucial for hazard assessment. At this stage, we cannot tell if the thin-skinned faults creep very fast, rupture seismically, or produce episodes of tremor and slip (ETS).

6 Summary and conclusions

1. The need for geohazard assessment in the marine environment is increasing globally. Yet, in the field of hazard maps for planning and building, the offshore regions are commonly lagging decades behind the onshore practice.

2. Mapping active faults in the marine environment is particularly complicated. If the onshore practice is followed, a Holocene horizon needs to be detected in the subsurface; then, for each fault, the question whether this horizon is displaced or not needs to be answered. This requires high-resolution seismic surveys and numerous coring and thus cannot be done for large regions.

3. In site-specific surveys, detailed bathymetry has become the main tool for mapping faults. Yet, we demonstrate that this is insufficient because fault scarps are decreased by sedimentation and erosion, particularly in sediment-rich environments such as continental margins.

4. Here, we take advantage of the marine environment (wealth of seismic data) to produce maps that cannot be produced onshore. First, we map a subsurface horizon dated to 350 ka in the entire study area. Second, we measure vertical fault displacements along this horizon. Third, we map fault planes combining several fault segments and measure their size.

5. By classifying all faults according to their vertical displacement and size, we prepare two hazard maps, which are further combined into a single simplified fault hazard map.

6. Our maps are particularly useful for master planning. The sedimentation rate map alone immediately reveals tectonically active grabens, and the hazard maps help define the most hazardous zones.

7. Using our maps, we revealed a particularly problematic zone in the upper slope south of the Dor disturbance. In this area, a series of big listric faults are characterized by large vertical displacements. The sedimentation rate in this location is also exceptional – 4 times faster than the vertical displacement rate – and, still, fault scarps are prominent. We suggest that this indicates rapid creep, seismic rupture, or episodic motions.