Post-LGM intensification of marine faulting: resolution-dependent hazard assessment

In marine environments, site-specific fault hazard assessments conventionally calculate average slip rates by measuring the displacement of the youngest mappable stratigraphic horizon. However, due to data resolution limits, this horizon is typically tens to hundreds of thousands of years old. Consequently, hazard surveys provide long-term averages that may mask recent, non-linear fault behavior. In this study, we examine how increasing stratigraphic resolution affects the calculated slip rates used for offshore seismic design. We investigate thin-skinned normal faults offshore Israel that pose hazards to major pipelines delivering gas to onshore power plants. Previous assessments measuring displacements of a 350 ka horizon yielded average slip rates of 0.25 mm yr⁻¹. However, by utilizing higher-resolution seismic data to measure displacements across a 14 ka horizon, we calculated a recent slip rate of 2.4 mm yr⁻¹. This tenfold increase indicates non-linear slip rates and raises the hypothesis that rapid post-LGM (Last Glacial Maximum) sea-level rise influenced the increased faulting. To examine this, we extended our time window to the latest Pleistocene, thereby demonstrating a correlation between sea-level fluctuations and variations in faulting. The finding that fault slip rates increased after the last glacial period has two main implications. First, it supports the hypothesis that rapid sea-level rise drives increased faulting, possibly due to changes in pore pressure along thin-skinned faults and detachment surfaces, contributing to the understanding of thin-skinned fault mechanics. Second, it suggests the need to use post-LGM stratigraphic horizons as seismic markers for geomarine hazard analysis, particularly in circum-Mediterranean margins where the Messinian salt giant propels salt tectonics.

1 Introduction

While building codes on land generally prohibit construction over active faults (e.g., Hart and Bryant, 1999), such regulatory constraints are largely absent in the marine environment. Consequently, when subsea infrastructures must cross active or potentially active faults, two measures are recommended. First, an optimal route crossing the least hazardous faults should be selected. A method for that was recently presented by Laor and Gvirtzman (2023), who classified all faults in their marine study area into relative hazard levels, allowing identification of the least hazardous faults. Second, site-specific surveys are needed to characterize faults that are crossed by the selected route. One way to characterize faults is to apply a probabilistic fault displacement hazard analysis (PFDHA; Youngs et al., 2003), typically used in onshore environments, following a probabilistic approach developed for ground-motion prediction (PSHA; Cornell, 1968; McGuire, 1976; Bender and Perkins, 1987). While the traditional “earthquake approach” of PFDHA calculates hazard curves based on earthquake recurrence intervals (magnitude-frequency relations), an alternative “displacement approach” can be applied in some marine environments without relying on seismic mechanisms (Youngs et al., 2003; Fiorini et al., 2023). However, implementing this method requires calculating the frequency of displacement events, which relies on precise determinations of the overall fault slip rate divided by the average displacement per discrete event (Youngs et al., 2003). In offshore settings, the lack of direct observational capabilities (e.g., paleoseismic trenching) and the dependency on geophysical surveys often limit the availability and accuracy of such data. In particular, for the shallow, thin-skinned faults studied here (offshore Israel), discrete slip events in geophysical data are uncertain to detect, making it difficult to determine whether these faults creep continuously or exhibit stick-slip behavior. Because separating discrete slip events from continuous deformation is currently impossible in this setting, as well as in other marine settings, probabilistic approaches such as PFDHA remain infeasible.

In practice, fault hazard surveys frequently rely on deterministic evaluations, making the direct calculation of fault slip rates the primary and most critical metric for structural design. In practice, site-specific surveys aimed at determining these recent slip rates for marine seismic design measure the displacement of the youngest possible stratigraphic horizon, which is typically 10⁴ or 10⁵ years old, depending on the available resolution of seismic and core/well data. Thus, hazard surveys commonly provide average slip rates that represent the recent tens to hundreds of thousands of years. In this work, we demonstrate how increasing stratigraphic resolution affects marine fault hazard assessment. We utilize the Levant continental margin, offshore Israel, as a natural laboratory to investigate thin-skinned normal faults propelled by salt tectonics (Fig. 1). These faults create high seabed steps and extend directly across Israel's national infrastructure corridor, raising the need for precise slip-rate quantification to ensure subsea asset resilience. Leveraging an extensive, high-quality seismic reflection dataset and ongoing deformation driven by the unstable Messinian salt giant, we aim to characterize fault kinematics and slip rates at a higher resolution than previously possible in this region.

Our results reveal that slip rates calculated after the last glacial period have increased by an order of magnitude. This finding is significant for two reasons: (1) It suggests that climate-driven water and sediment loading may increase fault activity. (2) The necessity of post-LGM stratigraphic markers in marine site-specific surveys aimed at determining fault slip rates.

Given that thin-skinned deformation is widespread across the Mediterranean, our findings here serve as a benchmark, providing a template for reexamining seismic hazards in similar salt-tectonic settings across additional Mediterranean sites.

Figure 1 (A) Regional elevation map. ESM- Eratosthenes Sea Mount. LB- Levant Basin. The study area (shown in B) is marked in a red polygon. Digital elevation model (DEM) from Hall (1994) and Gvirtzman et al. (2015). (B) Location map of the study area offshore northern Israel. High-resolution (20×20 m) colored bathymetric map derived from data collected for this study along the seismic survey track, overlaid on a background of white depth contours derived from the Gvirtzman et al. (2015) DEM. The bathymetric background map with the color scale shown in (A). Black and pink lines show faults from Groups 1 and 2, respectively (Group 2 is the focus of this study), after (Laor and Gvirtzman, 2023). The yellow line indicates the salt pinch-out. White lines mark the seismic sections featured in this paper. F1–F4 denote the faults analyzed herein.

2 Background

2.1 The effect of glacial and interglacial periods on thin-skinned deformation

Over the past decades, submarine mass failures have been widely investigated along continental margins and are now recognized as one of the most significant geological processes shaping continental slopes (e.g., Masson et al., 2006; Owen et al., 2007; Lee, 2009; Tappin, 2010; Katz et al., 2015; Gadol et al., 2020). The mechanisms of slope failure are generally explained by a combination of preconditioning factors, such as rapid sediment accumulation, the presence of methane hydrates, the build-up of high pore-water pressures, and triggers, such as seismic activity or changes in hydrostatic loading. In this context, fluctuations in sea level are increasingly considered an important factor: sea-level fall reduces hydrostatic pressure and may destabilize gas hydrates, while rapid sea-level rise enhances pore pressure and reduces slope stability (Paull et al., 1996; Trincardi et al., 2003; Maslin et al., 2004; Wien et al., 2006; Sammartini et al., 2019; Li et al., 2025). Nevertheless, the global correlation between the frequency of submarine landslides and sea-level changes remains debated and has not yet been unequivocally demonstrated (Urlaub et al., 2013).

The present study touches on another process that may be affected by sea-level fluctuations – thin-skinned faulting. It has already been suggested that the water load imposed on continental shelves during sea-level rise induces lithospheric flexure that steepens slopes and enhances shear stresses along fault systems (Mitrovica and Peltier, 1991). Modeling results further indicate that rapid Holocene sea-level rise was sufficient to increase Coulomb failure stress by more than 1 MPa, a magnitude large enough to reactivate existing faults (Luttrell and Sandwell, 2010; Brothers et al., 2013).

A prominent case study for the destabilization of a continental margin is the Holocene Storegga Slide offshore Norway, one of the largest documented submarine landslides. This event occurred during a period of accelerated sea-level rise, amplified by meltwater discharges from proglacial lakes Agassiz–Ojibway in North America. The combination of rapid relative sea-level rise, increased pore pressure in slope sediments, and synchronous seismic activity is thought to have triggered the catastrophic failure (Bryn et al., 2005; Haflidason et al., 2005; Bondevik et al., 2012). The Storegga case thus illustrates how sea-level forcing, coupled with fault reactivation and slope instability, can jointly produce large-scale submarine failures with far-reaching geohazard implications (Smith et al., 2013).

Despite these indications, direct observational evidence linking sea-level variations to thin-skinned marine faults remains sparse. It is unclear whether the same processes that destabilize continental slopes with a vertical scale of tens or a few hundreds of meters operate at larger scales of a few thousand meters to influence thin-skinned faults. Hydrostatic pressure perturbations generated by water loading-unloading cycles are known to decay rapidly at depths shallower than 1 km (Locat and Lee, 2002), while seismic monitoring typically targets deeper crustal levels. In light of this knowledge gap, here we apply high-resolution geophysical data capable of examining whether or not variations in fault slip rates correlate with sea-level variations.

2.2 Marine construction standards

Modern offshore construction standards reflect the challenges and complexities associated with installing infrastructure on the seabed, particularly in geologically dynamic settings. Unlike onshore practice, where hazard mapping and fault characterization rely on decades of established protocols, the marine environment lacks universally accepted procedures for hazard assessment and infrastructure design (Laor and Gvirtzman, 2023). This gap is largely due to the intrinsic difficulties in accessing, imaging, and dating marine faults, as well as the higher sedimentation and erosion rates characteristic of many continental margins (Kvalstad, 2007; Hough et al., 2011). Consequently, decisions regarding the routing and protection of seafloor infrastructure (such as pipelines or cables) often require integrating multiple high-resolution seismic and bathymetric datasets with site-specific geological and geotechnical investigations (Prior and Hooper, 1999; Angell et al., 2003).

A key difference in offshore geohazard assessment is the reliance on extensive seismic mapping, which is typically higher in both quality and quantity than comparable surveys on land. This abundance enables regional-to-local-scale classification of fault activity and displacement hazards without the need for exhaustive coring and stratigraphic dating for every fault encountered (Prior and Hooper, 1999; Angell et al., 2003). In practice, the inability to definitively date fault activity in the marine environment has led some standards to treat all seabed-rupturing faults as potentially active, adopting a precautionary, engineering-first approach (On, 2016). Alternatively, probabilistic analyses comparable to those used in earthquake hazard assessment have been adapted to estimate the expected frequency and size of fault displacement events, guiding both the early planning and the detailed engineering phases of offshore projects (Wong and Stepp, 1998; Angell et al., 2003; Youngs et al., 2003).

Recent approaches emphasize pragmatism and risk minimization at early planning stages. Laor and Gvirtzman (2023) propose classifying offshore faults by measurable parameters such as recent (e.g., 350 ky) vertical displacement and fault-plane size, using broad seismic data coverage. This method enables infrastructure planners, working under uncertainty, to route pipelines or cables across zones of lowest apparent hazard, focusing costly site-specific studies only where unavoidable crossings occur. Such pragmatic classification provides a functional hazard map explicitly designed for marine conditions, while acknowledging the current limits of fault activity characterization at sea.

However, the efficacy of these approaches relies on our ability to constrain the activity of such faults. Understanding recent fault displacement rates is therefore essential for developing region-specific offshore construction standards and for more accurately assessing the seismic and surface rupture hazards posed by marine faults. Improved knowledge in this area will ultimately inform the design of infrastructure and guide appropriate engineering adaptations to ensure both safety and resilience in offshore developments.

2.3 Thin-skinned tectonics offshore Israel

The Israeli continental shelf and slope comprise Plio-Quaternary siliciclastic sediments derived from the Nile Delta through longshore currents (Emery et al., 1960; Gvirtzman and Buchbinder, 1978; Goldsmith and Golik, 1980; Carmel et al., 1985; Tibor et al., 1992; Buchbinder et al., 1993; Buchbinder and Zilberman, 1997; Ben-Gai et al., 2005; Frey-Martínez et al., 2006; Segev et al., 2006; Zviely et al., 2006, 2007; Klein et al., 2007; Schattner et al., 2015; Schattner and Lazar, 2016; Gadol et al., 2020; Zucker et al., 2021). The stratigraphy of this Plio-Quaternary section has been defined through several frameworks. Initially, the Plio-Quaternary section was broadly divided into two primary units: the Pliocene (unit-1) and the Pleistocene (unit-2), separated by regional unconformity (Schattner et al., 2017). The Pleistocene section is further subdivided internally by a prominent mid-Pleistocene unconformity (UC1 in Schattner et al., 2012) into a lower sequence (SU2a), characterized by sigmoidal progradation and aggradational depositional patterns, and an upper sequence (SU2b) marked by oblique tangential clinoforms (Schattner et al., 2017; Schattner and Lazar, 2016). Later, higher-resolution studies using integrated bio- and seismo-stratigraphic data subdivided the Plio-Quaternary section into four main depositional units: Unit 1 corresponds to the Pliocene (5.33–2.6 Ma), Unit 2 to the Gelasian (2.6–1.8 Ma), Unit 3 to the Calabrian-Ionian (1.8–0.35 Ma), and Unit 4 to the Ionian-Holocene (< 0.35 Ma) (Elfassi et al., 2019).

The Plio-Quaternary sedimentary package overlies the Messinian evaporite sequence, which lies approximately 1 km below the seafloor. The Messinian salt layer, deposited during the Messinian Salinity Crisis (5.97–5.3 Ma), reaches a maximum thickness of approximately 1.8 km in the deep basin and thins dramatically beneath the continental slope to a few tens of meters (Netzeband et al., 2006; Cartwright and Jackson, 2008; Steinberg et al., 2011 Gvirtzman et al., 2017; Manzi et al., 2018; Meilijson et al., 2018). The line at which the evaporitic sequence thins below seismic resolution is termed the salt pinch-out line (Fig. 1b; Hsü, 1977; Rowan et al., 1999; Ryan, 2009).

The post-Messinian structural evolution of the Levant Basin is a multiphase process characterized by two distinct thin-skinned tectonic domains that operate simultaneously but are driven by different mechanisms (Ben Zeev and Gvirtzman, 2020). Approximately 2.7 million years after cessation of salt deposition (Pliocene-Gelasian transition), continental uplift and margin tilting initiated westward salt flow, which in turn carries the overlying rock column basinward (Cartwright and Jackson, 2008; Elfassi et al., 2019; Hamdani et al., 2021; Moneron and Gvirtzman, 2025). This westward motion of rocks overriding the salt relative to the stable salt-less upper slope produces margin-parallel (NNE-trending) extension and normal faulting offshore Israel. Subsequently, a contractional domain developed in the deep basin, featuring folds and thrust faults trending NW to WNW, forming the circum-Nile deformation belt (CNDB; Netzeband et al., 2006; Cartwright and Jackson, 2008; Gvirtzman et al., 2015; Ben Zeev and Gvirtzman, 2020; Zucker et al., 2020; Moneron and Gvirtzman, 2025). This deep-basin shortening is propelled by a pressure gradient, causing outward salt expulsion (out-squeezing) due to the massive sedimentary loading of the advancing Nile Delta (e.g., Loncke et al., 2006; Gvirtzman et al., 2013; Hübscher and Netzeband, 2017; Zucker et al., 2020). The CNDB deformation began at approximately 1.8 Ma, coinciding with the establishment of a continuous continental shelf offshore Israel, which serves as a marker for the Nile Delta's maturation (Ben Zeev and Gvirtzman, 2020). Where these two systems converge in the southern Levant Basin, their orthogonal motions – the westward gliding of the continental slope and the northeastward advancing of the CNDB – interact and mutually restrain each other. This convergence is facilitated by strike-slip faults that allow the deep basin overburden to slide north-northeastward relative to the continental margin (Ben Zeev and Gvirtzman, 2020; Laor and Gvirtzman, 2026).

The extensional zone offshore Israel exhibits hundreds of normal faults that offset the seafloor, creating tens of meters high seabed steps (Neev et al., 1976; Ben-Avraham, 1978; 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; Loncke et al., 2006; Netzeband et al., 2006; Mart and Ryan, 2007; Cartwright and Jackson, 2008; Clark and Cartwright, 2009; Cartwright et al., 2012; Gvirtzman et al., 2013, 2015; Katz et al., 2015; Hübscher and Netzeband, 2017; Gadol et al., 2020, 2025; Ben Zeev and Gvirtzman, 2020; Hamdani et al., 2021; Laor and Gvirtzman, 2023, 2026). These faults form primarily above the salt wedge, as documented in numerous studies. However, in the Dor structure area offshore northern Israel (Fig. 1b), a distinct group of listric faults is located landward of the salt pinch-out line. These faults (Group 2 in Laor and Gvirtzman, 2023) have been classified within the high-hazard category because they are relatively big with relatively large recent displacements, compared with surrounding faults. For comparison, recent (350 kyr) slip rates of faults located above the salt wedge are approximately 0.05 mm yr⁻¹ (Elfassi et al., 2019; Laor and Gvirtzman, 2023), while recent (350 ky) slip rates of the big listric faults near Dor structure (Group 2; Fig. 1b) are approximately 0.25 mm yr⁻¹ (Laor and Gvirtzman, 2023). Moreover, while slip rates of the salt-wedge faults have decreased from 0.11 mm yr⁻¹ in the Gelasian (2.6–1.8 Ma) to 0.05 mm yr⁻¹ recently (350 kyr; Elfassi et al., 2019), the faults near the Dor structure reached high slip rates in the past 350 ky. Given their elevated activity levels and strategic location within Israel's primary national infrastructure corridor, these Group 2 faults (pink line in Fig. 1B) constitute a critical focus for seismic hazard assessment, prompting this study to concentrate specifically on these structures.

3 Methods

3.1 High-resolution seismic reflection acquisition and processing

A sub-meter high-resolution multi-channel seismic survey was conducted offshore northern central Israel in June 2021 (Fig. 1). The survey consisted of eight seismic lines designed to cover the shelf edge and upper slope area (200–1200 m water depths). The seismic data were acquired using a 0.5 to 3 kHz Geo Marine Survey Systems Sparker source (Geo-Source 400) deployed with a 48-channel Geo Marine Survey Systems streamer, with a group interval of 3.125 m. To ensure optimal spatial sampling while profiling to a depth of up to ∼ 200 m, a multi-pinging, continuous-recording approach was applied. The acquired ultra-high-resolution (UHRS) sparker seismic data display a broad frequency spectrum ranging from ∼ 300 to 1500 Hz, with the dominant peak energy concentrated between 400 and 600 Hz. Assuming typical interval velocities of 1550–1650 m s⁻¹ for the shallow unconsolidated marine sediments (Hamilton, 1972) and applying the quarter-wavelength criterion (Sheriff and Geldart, 1995), the estimated vertical resolution ranges between ∼ 0.25 and 0.7 m. This decimeter scale resolution enables detailed imaging of shallow faulting and small-scale deformation structures.

Raw seismic data were loaded into Emerson-Paradigm (2022; now AspenTech SSE) software for processing and interpretation. Geometry and common midpoint (CMP) binning at a 1.56 m spacing were computed using custom MATLAB scripts and assigned to the data. Next, seismic traces were filtered to enhance signal quality, followed by amplitude scaling to balance signal strength across the dataset while preserving relative reflection amplitudes. Deconvolution steps, including de-signature and de-ghosting corrections, were performed using the ECHOS DeconQ module to reduce source and receiver ghost effects. Finally, a pre-stack time migration was applied to focus reflection events and improve structural interpretation. To accurately measure fault displacement in depth (meters), the two-way travel times (TWT) of the targeted reflectors were converted using specific interval velocities. While regional, deep-penetrating seismic data from hydrocarbon exploration provide a general velocity framework for the Plio-Quaternary section in the area (Sagy et al., 2015, 2020), the specific velocities for our shallow target interval (the upper ∼ 200 m) were determined based on typical interval velocities for shallow, unconsolidated marine sediments in this setting (Hamilton, 1972), applying a linear velocity-depth gradient. Accordingly, the time-to-depth conversion for the last 250 kyr was performed using the following interval velocities: 1550 m s⁻¹ for the uppermost Unit 4D (0–14 ka), 1600 m s⁻¹ for the underlying Unit 4C (14–140 ka), and 1650 m s⁻¹ for the deeper Unit 4B (140–250 ka). This approach accounts for the gradual downhole compaction of the shallow marine sediments and ensures that depth conversions are anchored in the locally measured seismic properties.

3.2 Horizon dating using seismo- and sequence stratigraphy

Sequence stratigraphy deals with repetitive units where the depositional conditions that prevailed during their sedimentation were similar, and they are bounded by erosional or non-depositional unconformity surfaces (Galloway, 1989). These surfaces, which separate the units, formed during sea-level fall or rise (Vail et al., 1977).

The physical formation of these erosional unconformities during a relative sea-level fall is primarily driven by the basinward shift of the shoreline and the lowering of the base level (Posamentier et al., 1992). During sea level drops, the previously submerged continental shelf is sometimes exposed to subaerial conditions. This exposure leads to the development of a subaerial unconformity through terrestrial weathering processes, including fluvial erosion, sediment bypass, and pedogenesis (Sloss et al., 1949). Concurrently, near the shelf edge and upper slope, the drop in sea level physically lowers the wave base. This downward shift forces the wave energy to actively scour the progressively shallowing seafloor, stripping away previously deposited strata and creating a regressive surface of marine erosion (Plint, 1988). Together, these subaerial and shallow-marine erosional mechanisms remove older sequences and generate the prominent unconformity surfaces observed at the basin margins (Sloss et al., 1949; Plint, 1988; Posamentier and Allen, 1999; Catuneanu et al., 2006, 2011; Catuneanu, 2017).

In this work, we identified unconformities that separate units with different seismic facies using seismic interpretation and horizon picking. Our picking strategy was designed to systematically track true stratigraphic boundaries rather than strictly following a specific phase (e.g., a peak or trough) of the seismic wavelet. Accordingly, a mapped horizon may coincide with a distinct seismic reflection in certain locations, while in others, it represents a continuous transition zone between different reflectivity packages (seismic facies). This facies-guided approach ensures that our structural measurements accurately capture the displacement of physical geological boundaries, accounting for local variations in the seismic wavelet signature. We defined sequence boundaries (SB) according to Hunt and Tucker (1992) approach, so that they represent formation at the end of the regression. The geometry of the SB reflections truncated the top set of the former reflections by Toplap or Truncation geometry. The Truncation geometry is characterized by a high angle between the truncated reflection and the erosional reflection, and usually forms erosional topography. The Toplap geometry is characterized by a moderate angle and topography. Other criteria for SB identification are the Onlap relationship in the foreset and topset regions on the erosional surface.

Our stratigraphic framework is calibrated by two independent chronological anchors: 14 ka (late Pleistocene) and 350 ka (mid-Pleistocene), providing a foundation for our age model. While previous studies in the Levant margin have utilized seismic and sequence stratigraphy to delineate the late Quaternary succession (Schattner and Lazar, 2016), our interpretation aligns with these regional stratigraphic schemes. Although the intermediate horizons identified in this study cannot be directly correlated with every surface proposed in earlier works, our framework shows no significant contradictions. Rather, it complements regional observations, emphasizing that while sequence stratigraphy provides an essential relative temporal framework, the precision of our displacement measurements relies primarily on these well-constrained anchor horizons (14 and 350 ka).

3.3 Reconstructing the history of displacement rates

To calculate the displacement accommodated by the fault during each period, we reconstructed the displacement of the oldest horizon over time. The vertical displacement across faults was measured for each geological horizon. The relative movement of the hanging wall against the footwall was reconstructed by calculating the cumulative offset of the earliest horizon and representing it as a function of time (Fig. 2).

Figure 2 The methodology for backward reconstruction of the oldest horizon movement along the faults. The vertical offset was measured for each horizon, and the relative movement of the hanging-wall block relative to the footwall block was reconstructed over time.

Slip rates for discrete time periods were determined by dividing the incremental displacement of the basal reflection by the duration of each chronostratigraphic interval. This analytical procedure was applied to four faults (F1–F4) that exhibit clear displacement of all traced horizons. These structures were selected because they are representative of the “Group 2” fault population (Laor and Gvirtzman, 2023), which is characterized by the highest average displacement rates observed over the last 350 ka. Furthermore, because they lie within the national infrastructure corridor of offshore Israel, these faults are a primary focus of regional seismic hazard assessment. Preliminary mapping of this fault population confirms that these structures (F1–F4) exhibit displacement patterns and slip-rate trends comparable to other faults in the vicinity, supporting their suitability as representative benchmarks for the area's kinematics.

4 Results

4.1 Chrono-Seismo-Sequence stratigraphy

Figure 3 shows the division of the Plio-Quaternary section into four units based on the seismo- and bio-stratigraphy of Elfassi et al. (2019), with the base of Unit 4 dated to 350 ky (light blue).

Figure 3 Seismic section from industry seismic data crossing the Group 2 faults (marked as black lines; Laor and Gvirtzman, 2023). Division of the Plio-Quaternary section into four dated units is based on Elfassi et al. (2019). Base unit 2 (blue horizon) is dated to 2.6 Ma; Base unit 3 (pink horizon) is dated to 1.8 Ma; Base unit 4 (light blue horizon) is dated to 350 ka. For the location of the line, see Fig. 1.

Figure 4a further subdivides Unit 4 into four subunits (4A–4D), based on the new high-resolution survey conducted in this study (details in the Methods section above). The base of the youngest subunit 4D (yellow reflection) is interpreted as “Surface A” identified by Neev et al. (1976), Almagor et al. (2000), Schattner et al. (2010), and Gadol et al. (2025), marking the end of the Last Glacial Maximum and the onset of sea-level rise. Recently, this horizon was further dated to approximately 14 ky (Hamann et al., 2008; Gadol et al., 2025).

Figure 4 (A) Subdivision of Unit 4 into four sub-units (4A–4D), based on a new high-resolution seismic survey conducted in this study. Light blue horizon is dated to 350 ka (following Elfassi et al., 2019); green horizon is dated to 250 ka (based on C); purple horizon is dated to 140 ka (based on C); and the yellow horizon is dated to 14 ka (Neev et al., 1976; Almagor et al., 2000; Schattner et al., 2010; Gadol et al., 2025), (B) Detailed view of the shelf edge zone showing geometries that were used to divide the section into sequences according to sequence stratigraphic theory. TWTT – two-way travel time. TST – transgressive system tract. HST – highstand system tract. LST – lowstand system tract. (C) Global sea level curve showing glacial and interglacial cycles (modified after Waelbroeck et al., 2002). Note the four sea-level drops at approximately 350, 250, 140, and 20 ka, which are used as sequence boundaries.

Based on sequence-stratigraphic principles, we identified the four horizons picked within Unit 4 as SBs underlain by truncations (red arrows; Fig. 4b) of high-stand systems tracts (HST) and overlain by onlap geometry (black arrows; Fig. 4b) of transgressive systems tracts (TST). Since all four horizons represent low sea level, and the yellow and light blue horizons are dated to 14 and 350 ka, respectively, we interpret the intermediate horizons to 250 ka (green) and 140 ka (purple), respectively, according to the global sea-level curve (Fig. 4c; Waelbroeck et al., 2002). In summary, we now have four horizons with estimated ages corresponding to four major sea-level drops within the last 350 ky (Figs. 4a–5).

Figure 5 High-resolution seismic section conducted in this study. Faults analyzed in this study are marked as black lines (same faults presented in Fig. 3). Subdivision of Unit 4 into four sub-units (4A–4D). Base unit 4A (light blue horizon) is dated to 350 ka; Base unit 4B (green horizon) is dated to 250 ka; Base unit 4C (purple horizon) is dated to 140 ka; Base unit 4D (yellow horizon) is dated to 14 ka. The inset provides a close-up view of the seafloor and ∼ 14 ka reflectors, illustrating the high resolution of the data and the similarity in offset magnitude across the reflectors. TWTT – two-way travel time. For the location of the line, see Fig. 1.

4.2 Displacement rates

First, we performed backward restoration of fault movement since 1.8 Ma on four faults identified in the industrial, relatively low-resolution seismic survey (e.g., Fig. 3); results of the displacement history and displacement rates are shown in Fig. 6a–b. Then, we performed the same analysis on the same faults (Fig. 5) utilizing the new stratigraphic subdivision of Unit 4 to 4A–4D (Fig. 6c–d). These results demonstrate that when slip rates are averaged over 350 ky, the obtained values are 0.1–0.25 m ky⁻¹, whereas higher-resolution analysis reveals that these average values can be separated into much lower values for Units 4B–4C and much higher values for Unit 4D. For practical purposes, considering the most recent displacements, slip rates have increased 10-fold.

Figure 6 (A) Displacement history of the 1.8 Ma horizon identified in the industrial, relatively low-resolution seismic survey. (B) Displacement rates based on A. (C) Detailed presentation of Unit 4 subdivided into three subunits (4B–4D), based on high-resolution data collected in this study. (D) Detailed displacement rates of the subdivided Unit 4. Note the relatively high rates averaged over the last 14 ky (Unit 4D).

In a wider view (Fig. 7), downslope faults located above the salt wedge (Group 1 in Laor and Gvirtzman, 2023), indicate the highest slip rates during deposition of Unit 2 (2.6–1.8 Ma) and a decrease ever since (Elfassi et al., 2019). However, upslope faults located landward of the salt wedge (Group 2 in Laor and Gvirtzman, 2023) illustrate increased activity in the last 350 ky and especially in the last 14 ky. This result highlights the hazard posed by Group 2 even more than estimated by Laor and Gvirtzman (2023).

Figure 7 Summary graph of displacement rates throughout the entire history of salt tectonics in the Levant Basin. Classification of faults into two groups follows Laor and Gvirtzman (2023). Displacement rates for Group I were calculated by Elfassi et al. (2019), reached 0.11 m ky⁻¹, averaged over the Unit 2 period (2.6–1.8 Ma). Displacement rates for Group II were calculated in this study based on industrial data, reaching 0.25 m ky⁻¹ averaged over the Unit 4 period (350 ka-present). When recalculated using the new high-resolution data from this study, Group II displacement rates average 2.4 m ky⁻¹ over the last ∼ 14 ky, a 10-fold increase compared to calculations based on industrial data.

5 Discussion

5.1 Displacement Rates, Sedimentation, and Erosion

Considering Unit 4 as a whole, its average sedimentation rate in the area of the analyzed faults is 2 m ky⁻¹ (Laor and Gvirtzman, 2023), which is an order of magnitude faster than the average slip rates (0.25 m ky⁻¹). This means that on a 350 ky scale, fault scarps are expected to be completely buried. Therefore, the observed noticeable fault scarps (tens of meters high; Fig. 8a) must indicate that displacement rates are not constant and that a very recent rapid slip has created a fault step that has not yet been buried (Fig. 8a). However, a simple calculation shows that even the high slip rates inferred for the last 14 ky are insufficient to explain the observed fault scarps unless these scarps are much younger than 14 ky. If faults move in a constant velocity (creep motion) of 2.5 m ky⁻¹, with a burial rate of 2 m ky⁻¹ (Laor and Gvirtzman, 2023), the expected height of the fault scarp after 14 ky should be ∼ 7 m, whereas in reality the observed fault scarps are a few times larger. The striking similarity between the displacement patterns of the 14 ky horizon and the present-day 30 m bathymetric offset (Fig. 8a) suggests that the prominent fault scarp formed during a very recent event, likely triggered by sudden seismic slip or rapid creep (Fig. 8b). This temporal uncertainty carries significant implications for seismic hazard assessment, as the preservation of such a sharp scarp implies a recent activity that remains poorly constrained within our current chronological framework.

Figure 8 (A) Location of fault F2 (highlighted in black). TWTT – two-way travel time. (B) Displacement history of fault F2, showing significant displacement at 14 ka. Based on displacement-to-sedimentation rate ratios calculated for various scenarios described in the text, the most recent displacement on this fault likely occurred within the last 14 ky (indicated by dashed pink line).

Two additional factors complicate the simple comparison between sedimentation and displacement. First, several of the horizons used for offset measurements show truncation and may represent erosional unconformities. Any removal of a section across these surfaces would alter the apparent offset and could either exaggerate or reduce the reconstructed slip. Second, the assumption of constant sedimentation rates within subunits is an approximation. Rates may have varied within Unit 4, and such variations would influence both burial and the interpreted timing of scarp formation.

5.2 Correlation with sea level changes

The recent increase in displacement rates and the remarkable preservation of the ∼ 30 m fault scarp on the present-day seafloor despite high sedimentation rates, raises a fundamental question: why did the activity of these faults rise so much during the last 14 ky? Our hypothesis is that this acceleration is intrinsically linked to the steep post-LGM 120 m sea-level rise. To further examine this hypothesis, we extended our analysis backward in time and illustrated a correlation between sea-level changes and fault motion along three upslope faults (Fig. 9). The colored dashed lines in Fig. 9b represent reconstructed motion across the faults, with age anchors at 250, 140, and 14 ky, corresponding to peak glacial periods, and linear age interpolation between these points based on relative unit thicknesses (Fig. 9a), assuming constant sedimentation rates within subunits. Comparison of the fault motions to the sea-level curve indicates a noticeable correlation between post-glacial sea-level rise and increased fault activity.

Figure 9 (A) Example of subdivision into undated sub-units for one of the faults (F4). TWTT – two-way travel time. (B) Correlation between displacement histories of F2, F3, and F4 faults (pink, light blue, and green dashed lines) and sea level change (the lower undulating blue line). Black dots indicate age anchor points. Between these points, the interpolated ages were calculated based on the ratio of unit thickness to time. Fault activity increases during periods of sea level change. Note that the displacements do not account for portions removed by erosion; therefore, original displacements may have been originally larger (for example, at 140 ka). Sea level curve modified after Waelbroeck et al. (2002).

Nonetheless, it remains uncertain whether the observed fault acceleration responds to the onset of sea-level rise or legs sometime after the quick rise. This ambiguity stems from both the resolution limits of the data and the unknown lag time between climate forcing and structural response. In any case, it seems that fault activity follows a cyclic pattern coinciding with glacial-interglacial transitions.

Because sea level has been relatively stable over the last ∼ 7 ky, any younger slip that contributes to the 30 m high present-day seabed scarp may represent a delayed response to the post-LGM sea level rise, or alternatively, require an additional cause unrelated to eustatic forcing; such causes may be related to localized pore-pressure changes, internal salt-tectonic dynamics, or other short-term processes not resolved here. In light of our results and the questions they raise, we propose that coupling between sea-level fluctuations and fault activity should be systematically tested along additional continental margins.

Figure 10 Sensitivity analysis of the calculated fault slip rates to interval velocity variations. The baseline slip rates (solid black line) were calculated using our standard velocity model (dashed black line). The upper (pink) and lower (light blue) bounds were derived by applying a ±100 m s⁻¹ perturbation to the interval velocities across all stratigraphic units. The convergence of these bounds with the baseline model demonstrates that the observed temporal kinematic trend, particularly the acceleration in slip rates since the Last Glacial Maximum (LGM), remains robust regardless of velocity uncertainties.

5.3 Fault Mechanism and response to loading

The observed acceleration in fault slip rates after the Last Glacial Maximum raises fundamental questions regarding the mechanical behavior of these faults, and particularly how they respond to water and sediment loading. The transition from terrestrial to marine fault analysis necessitates accounting for the reduction in effective normal stress on the fault plane, driven by the combined effects of the water column and post-glacial sediment loading. When loading increases and pore fluids cannot dissipate fast enough within fine-grained shales, pore pressure increases and frictional resistance along the fault plane is reduced (e.g., Dugan and Flemings, 2000; Flemings et al., 2003). This mechanical process is recognized in marine settings as a primary driver for fault instability and submarine mass movements (e.g., Sultan et al., 2004; Lafuerza et al., 2009). Although a direct correlation between measured pore pressure fluctuations and episodes of rapid fault slip has not yet been documented in the literature, such fluctuations may result from the sea-level variations documented here.

Another complexity related to the thinned-skinned faults studied here is the inability to determine whether these faults constantly creep or behave in a stick-slip manner with seismic rupture events (or a combination of the two mechanisms). The late Pleistocene reconstruction presented in Fig. 9 indicates that fault motion is far from being constant. However, whether the documented variations in slip rates express changes in the intensity of the stick-slip motion or changes in the velocity of the creep motion is currently unresolved.

5.4 Uncertainty in Fault Displacement Rates Measurements

The calculation of fault slip rates from seismic reflection data inherently involves uncertainties associated with time-to-depth conversion (seismic velocity model) and sediment compaction during burial. The conversion of vertical fault displacement measured in the two-way travel time (TWT) domain to the depth domain relies on the assumed interval velocities for the sedimentary units sequence. To illustrate the uncertainty introduced by this conversion, we tested the sensitivity of the calculated slip rates to velocity variations of ±100 m s⁻¹, which are considered an upper bound for velocity errors. As illustrated in Fig. 10, the uncertainty resulting from those variations is minor and negligible for our interpretation. The absolute slip rates remain well-constrained, confirming the validity of the observed kinematic trend regardless of the specific velocity bounds.

Figure 11 Sensitivity analysis of fault slip rates to sediment compaction. (A) Conceptual model showing the geometric relationship between measured displacement (D_i) and sediment compaction. The porosity-depth curve (right) illustrates the reduction in stratigraphic thickness (compaction) as a function of burial depth (d_r). (B) Comparison between original slip rates derived directly from seismic data (solid black line) and “decompacted” slip rates (dashed orange line). The latter were calculated by artificially increasing the measured displacement of pre-14 ka units by 20 % to account for potential stratigraphic thickness loss due to compaction. The results demonstrate that this conservative correction has a negligible impact on the calculated slip rates, reinforcing the robustness of the observed kinematic trend.

A secondary source of potential error in our calculations is sediment compaction that affects reconstructed unit thicknesses (Fig. 11). Note, however, that the average displacement rate calculated from 350 ky to the present and from 14 ky to the present do not depend on compaction. By definition, these rates are the measured displacement divided by the time interval. This means that the main conclusion, stating that post-LGM slip rates are an order of magnitude greater than the 350 kyr average, is independent of the compaction process.

The error associated with neglecting decompaction is related to intermediate stages, during which displacement d_r is reconstructed (Fig. 11A). This error is proportional to the reduction in porosity with depth, and considering that our measurements were taken at a depth range of 100–300 m. Standard porosity-depth curves for shales suggest a thickness reduction of approximately 10 % (Sclater and Christie, 1980; Fig. 11A).

To rigorously test the resilience of our findings, we modeled an extreme scenario by artificially increasing the measured displacement of buried pre-14 ka units by 20 %. Although an observed porosity reduction of ∼ 10 % across the 100–300 m interval suggests a more moderate thickness adjustment, we adopted this 20 % value as a conservative sensitivity bound to account for potential uncertainties in sedimentary compaction and lateral lithological variations (Sclater and Christie, 1980). This upper-limit approach shows that our results are minimally affected by compaction-related errors (Fig. 11B).

5.5 Implications for Hazard Assessment

Our findings demonstrate that displacement rates averaged over a 350 ky interval underestimate post-LGM rates, suggesting that such long-term averages mask higher-frequency accelerations in fault activity. Continuing this line of thinking, we further raise the possibility that an additional increase in stratigraphic resolution is needed to separate the rapid sea-level rise at 18–8 ka from the relatively stable sea-level in the last 7 ky. Figure 12 illustrates two possible scenarios related to the question of quick or delayed response of faults to sea level rise. Option 1 is that the displacement rate peaked in the early Holocene, exactly coinciding with sea-level rise, while Option 2 is that the displacement rate peaked in the late Holocene, lagging after the quick sea level rise by a few thousand years. To further constrain the very recent displacement rates in the late Holocene, we call for integration of ultra-high-resolution sub-seabed seismic with core dating that will allow measuring the displacement of early and late Holocene horizons to compare their slip rates. This will narrow the gap between offshore and onshore paleoseismic studies.

Figure 12 Illustrated plot showing increased displacement rates for each time interval based on average displacement rates. The green and red curves show two alternative scenarios for variations in displacement rate over the last 14 ky, both yielding the same average rate. Note that Option 1 (green) would provide relief in geohazard planning due to sea level stability in the late Holocene, while Option 2 (red) would worsen it. This graph highlights the importance of enhanced temporal resolution for geohazard assessment.

6 Conclusions

This study demonstrates that the stratigraphic resolution utilized in seismic interpretations has a fundamental impact on fault hazard assessment for marine infrastructure. By increasing the temporal resolution from long-term averages (350 ka) to a more recent interval (14 ka) using high-resolution seismic reflection data, we identified fault slip rates more than 10 times higher than previously reported. These findings indicate that conventional hazard surveys, which rely on older, deeply buried stratigraphic horizons, may significantly underestimate recent fault activity and the associated risks to marine infrastructure. The preservation of 30 m fault scarps despite high sedimentation rates indicates very recent fault activity, potentially within historical times. However, we still cannot determine whether these faults creep or rupture seismically.

Our analysis shows a correlation between glacial-interglacial sea-level cycles and fault displacement rates. We present the hypothesis that elevated pore pressure within low-permeability sediments, driven by rapid sea-level rise and associated sediment loading, acts as a primary mechanical driver for this instability. While this pore-pressure mechanism provides a plausible working hypothesis for the observed acceleration, further geotechnical in-situ measurements are required to confirm this causal link and to distinguish it from other potential forcing factors. For infrastructure design, these results highlight that long-term averaged slip rates of ∼ 0.25 mm yr⁻¹ may underestimate hazard potential compared to recent rates exceeding 2.4 mm yr⁻¹. If these accelerated rates are indeed a response to sea-level-driven changes in effective stress, they may represent a transient peak related to the post-glacial environment. Consequently, sub-meter-resolution stratigraphic analysis is essential to refine hazard estimates and account for non-steady-state fault behavior in the marine environment. The coupling between environmental forcing, specifically sea-level fluctuations, and thin-skinned fault activity introduces new considerations for geohazard forecasting. While we distinguish between specific sea-level changes and broader environmental forcing, our findings suggest that fault systems in these settings may exhibit a degree of sensitivity to external environmental triggers. This sensitivity has direct implications for infrastructure planning in the eastern Mediterranean, particularly as climate-driven sea-level rise continues. Future research should extend this high-resolution stratigraphic approach to additional Mediterranean sites and focus on subdividing Holocene sediments to further quantify the relationship between sea-level changes and fault kinematics.