Probabilistic and Deterministic Seismic Hazard Assessments of the area comprised between west Gulf of Cádiz and east Alboran Sea

The increased deployment of subsea infrastructures and the exploration of marine resources have heightened the need to assess seismic hazard on the seabed, especially in tectonically active offshore regions. The area between the Gulf of Cádiz and the Alboran Sea, which is rich in coastal and submarine assets, is located within the Ibero-Maghrebian Region (IMR), a seismically active region. While previous research has addressed seismic hazard in adjacent inland areas using deterministic and probabilistic approaches, few studies have focused on offshore zones. Moreover, existing models often overlook the amplification effects introduced by bathymetry and seafloor conditions, and they rely on ground motion prediction equations (GMPEs) derived from inland data. Consequently, seismic hazard in marine environments remains poorly constrained. This study aims to evaluate the feasibility and applicability of both probabilistic (PSHA) and deterministic (DSHA) seismic hazard assessments for submarine areas, using updated seismogenic zonation, a high-resolution bathymetric model, and GMPEs suitable for soft marine soils. The objective is to produce reliable peak ground acceleration (PGA) estimates at the seabed in a complex tectonic context. The analysis was conducted using the OpenQuake (OQ) engine for PSHA and a custom MATLAB© script for DSHA. The seismogenic sources and parameters were taken from the European Seismic Hazard Model (ESHM20), and four GMPEs were selected for their relevance to offshore or soft soil conditions. These GMPEs were benchmarked against six regional offshore earthquakes recorded in the IMR as a plausibility check. Hazard estimates were computed over a EMODnet bathymetric refined grid and expressed as PGA maps for various return periods (100, 475, 2475, and 5000 years). The results show that significant PGA values occur over key submarine fault systems such as the Gorringe Ridge, Horseshoe Plain, and Arzew Fault. PSHA reference results for onshore site conditions and clipped to the study area was used as a methodological contrast to provide a reference for how the hazard results obtained amplify PGAs. The use of seabed-adapted GMPEs resulted in higher and more realistic PGA values compared to inland models. This study demonstrates the applicability of seismic hazard assessment methods to offshore environments and highlights the importance of incorporating seafloor conditions into hazard modelling. The findings offer a methodological basis for improving the seismic design of subsea infrastructures in tectonically complex marine regions such as the IMR.

1 Introduction

The majority of ocean infrastructures are anchored or seated on the seabed in some way (e.g., Guo et al., 2021; Rui et al., 2024). The stability of these assets may be compromised by the mechanical properties of the seafloor (e.g., Liu et al., 2021; Randolph and Gourvenec, 2017; Veiskarami et al., 2017). Earthquakes and the processes they trigger, such as tsunamis and landslides, affect the dynamic behaviour of the ground, with the potential to cause significant damage to infrastructures and facilities (Román-de la Sancha et al., 2022; Shan et al., 2022). Consequently, attention should be given to investigate local effects on offshore structures (e.g., Kee et al., 2023; Li et al., 2023) along with undersea seismic hazard and the influence of surrounding seismogenic sources (Ercilla et al., 2021; Karthigeyan, 2022). Coastal and offshore facilities should be designed to withstand earthquakes and their coseismic effects. In offshore settings, the water column, seabed morphology and compliant soft sediments rheology modify source-to-site path and site terms. The water layer alters impedance contrasts and incidence/refraction of P and SV waves at the seafloor, shifts spectral content, and can lengthen significant duration at the seabed relative to adjacent onshore sites (Chen et al., 2024a; Dhakal and Kunugi, 2023; Nakamura et al., 2015; Zhang and Zheng, 2020). These effects, together with basin geometry and bathymetric geography, motivate use of seabed-adapted GMPEs and 3-D distance metrics that include depth (e.g., R_RUP) when estimating seabed PGA (e.g., Chen et al., 2024a; Dhakal et al., 2021; Diao et al., 2014; Hu et al., 2020; Lan et al., 2021) . However, the strategies to evaluate undersea seismic hazard may not be entirely valid considering they rely on seismic response of emerged terrain.

The area encompassing the Strait of Gibraltar, between the Alboran Sea in the Mediterranean Sea and the Gulf of Cádiz in the Atlantic Ocean, is also known as the Ibero-Maghrebian Region (IMR) (Buforn et al., 1995). A substantial number of communication and energy cables (TeleGeography, 2026) are present in the IMR. Major infrastructures are also located here, including oil and gas pipelines connecting Algeria and Morocco with Spain (Medgaz; https://www.medgaz.com/, last access: 31 July 2026) and several offshore wind farms. Moreover, the European-African tunnel is planned to be developed in the IMR (SECEGSA; https://www.secegsa.gob.es/proyecto-del-enlace-fijo, last access: 31 July 2026). These assets, among other issues, underscore the area's substantial geostrategic importance (Báez et al., 2021). The IMR is at the boundary between the Eurasian and Nubian tectonic plates (Buforn et al., 2016). It comprises the southern Iberian Peninsula and the Maghreb (Buforn et al., 1995), and it is considered to have a medium to moderate tectonic activity, with frequent seismicity due to the convergence of the tectonic plates. Significant and devastating earthquakes have been documented within the region (Buforn et al., 2015), some of which have triggered tsunamis and submarine landslides (Table 1). Consequently, a comprehensive Seismic Hazard Assessment (SHA) is of utmost relevance within the IMR.

Nowadays, there are two broad standard approaches to SHA (Hamidatou et al., 2024; Kumar et al., 2024): Deterministic Seismic Hazard Analysis (DSHA) and Probabilistic Seismic Hazard Analysis (PSHA; Reiter, 1990). These methods have also been criticised (Kossobokov and Panza, 2022). Currently, PSHA has been broadly adopted by the scientific community and for mitigation plans, whilst DSHA is still considered useful in worst-case scenario modelling situations (Grasso and Maugeri, 2012; Mostafa et al., 2019), commonly used by civil protection agencies. DSHA was the first methodology proposed in the late 1960s for the design of nuclear power plants and was later applied to other large industrial infrastructures (AEIS-IGN, 1979; NRC, 1973). DSHA and PSHA can complement each other, providing further information about the existing seismic hazard (Wang et al., 2012).

Table 1 Location, date and moment magnitude (M_w) of notable earthquakes with casualties in the IMR.

Despite the significant record of earthquakes, submarine landslides and tsunamis in the IMR (Rodriguez et al., 2017; Vázquez et al., 2022), the area has not been updated in full SHA. The first attempt of SHA in IMR was made by Molina Palacios (1998). This work evaluates the seismic hazard in terms of intensities and PGAs of the Iberia-Africa contact zone by a probabilistic approach. This involved a seismic zonation using areal sources, applying available estimation tools to flat terrain, using ground motion prediction equations (GMPEs) from rocky inland areas, and relying on seismotectonic knowledge from the 1990s. Conversely, the emerged land around the IMR has exhaustive seismic hazard studies, using both deterministic and probabilistic approaches. For example, the PSHA method has been the most widely implemented (e.g., Crespo et al., 2014; Danciu et al., 2024; IGN-UPM, 2017; Salgado Gálvez et al., 2015). Neo-deterministic Seismic Hazard Assessment (NDSHA) has been applied for the Iberian Peninsula (García-Fernández et al., 2022) and DSHA for Sevilla city (Sá et al., 2021). In the North African part, Poggi et al. (2020) performed a PSHA, and Mourabit et al. (2014) also conducted an NDSHA approach in this region. At the European continental scale, the most recent work covers the Euro-Mediterranean region, where the IMR is located. A fully probabilistic framework was adopted in the implementation of the hazard model, using harmonised input datasets across national boundaries. The seismogenic zone mosaic uses the 2020 European Seismic Hazards Model (ESHM20) (Danciu et al., 2024) for PSHA hazard assessment. The ESHM20 model uses the same principles as the 2013 Seismic Hazard Harmonisation in Europe (SHARE) project model (ESHM13) (Wössner et al., 2015), with state-of-the-art procedures applied homogeneously to the pan-European region. However, the SHA gap in the marine environment persists, hence further assessment of large-scale seismic hazard under seabed conditions within the marine IMR is needed, taking into consideration earthquake characteristics and spatial distribution based on comprehensive undersea data. In order fill this gap, this study has used seismic zoning and parameters from the ESHM20 (Danciu et al., 2021a, 2024). The novelty of this work also involves the use of attenuation functions derived from analogous marine regions regarding their seismotectonic context, considering there are no estimates of GMPEs for the IMR at this time. The DSHA here has been carried out using an in-house MATLAB© code, while for the PSHA the OpenQuake engine software was used (OQ v3.21) (Pagani et al., 2014; Silva et al., 2014).

2 Study area

The study area is characterized by a complex geodynamic context resulting from the interaction between the Nubia and the Iberian tectonic plates (Fig. 1) (e.g., Custódio et al., 2016; Neres et al., 2016). The plate boundary of the study area is usually divided into three main areas according to their stress regime: the SW as simple shear zone, the pure shear Gulf of Cádiz as the Eastern Betic, and the South Moroccan compressional arc. Tectonic models of NW-SE to WNW-ESE show oblique convergence between the Nubia and Iberian plates (Herraiz et al., 2000; Reilly et al., 1992) with displacements of 2 to 5 mm yr⁻¹ (Fernandes et al., 2007; Nocquet, 2012) of the western Betic Cordilleras with respect to the Iberian Massif (Gonzalez-Castillo et al., 2015; Palano et al., 2013).

The Gulf of Cádiz (western sector of IMR) is located on the contact boundary between the Eurasian and African plates. It extends from the Azores Islands to the Strait of Gibraltar, through the Iberian massif and the Algarve region to the north, and, in the eastern part, through the orogenic arc of the Betic-Rifeña mountain range. It is characterised by a NW-SE oblique convergence regime, mainly controlled by the formation of the Betic-Rifeño Orogen and by the accommodation of post-orogenic compressional tectonic activity (IGN, 2025a). This tectonic configuration makes seismic activity in the Gulf of Cádiz significant, with moderately deep earthquakes (h<40–60 km) (Custódio et al., 2016). Some moderate magnitude earthquakes have also been recorded (Martín-Dávila and Pazos, 2003). Most of the focal mechanisms in this area are reverse and rifting type. There is also historical and instrumental evidence of several high-magnitude earthquakes in the region. These include the earthquakes of 1755 (M_w7.7–8.5; Chester, 2001; Fonseca, 2020; Mendes-Victor et al., 2009) and 1969 (M_w7.8; López Arroyo and Udías, 1972), both causing tsunamis that affected the coasts of Morocco, Portugal, and Spain. The Gorringe Bank, the São Vicente submarine canyon (São Vicente fault), Horseshoe, Marquês de Pombal, and the Southwest Iberian Margin (SWIM) faults show the greatest magnitudes in the region (M_w7.2–8.4; Grevemeyer et al., 2017; Martínez-Loriente et al., 2013) surrounded by the Gorringe Ridge, the Coral Patch, the Horseshoe and Seine abyssal plains.

Figure 1 Location of the area of interest of this work in relation to the tectonic setting. Abbreviations are as follows: WAB: Western Alboran Basin, EAB: Eastern Alboran Basin, AOB: Algerian Oceanic Basin, SVF: São Vicente Fault, MPF: Marquês de Pombal fault, SWIM: Southwest Iberian Margin lineaments (1, 2 and 3). AIF: Al Idrisi fault, YSF: Yusuf fault. Background: shadowed Digital Terrain Model (altimetry and coloured bathymetry) in the study area (extracted from EMODnet; https://emodnet.ec.europa.eu/en/bathymetry, last access: 24 August 2026).

The Alboran Sea is the westernmost part of the Mediterranean Sea and the oriental sector of IMR. It is bordered by the Alpine Mountain ranges of the Betic in the Iberian Peninsula and the Rif in the north of Africa. It is a contact zone between the converging Eurasia and Nubia tectonic plates. Among the tectonic structures, the submarine mountain ranges stand out, with a length of more than 50 km, oriented NE-SW and limited to the north and south by reverse faults with opposite dip. Several sets of conjugate directional faults are also found in this area: the NNE-SSW sinistral faults, such as the Al-Idrisi; the NE-SW sinistral faults, such as the Carboneras; and NW-ESE dextral faults, such as the Yusuf. Most of the seismic events recorded in the Alboran Sea have lower magnitudes (M_w<4; IGN, 2025b). The seismic activity in the Alboran Sea is interconnected with the southern Iberian Peninsula through the tectonic system running E-W from the Rif and Alboran Ridge, Eastern Alboran Basin, to the Arzew faults and Tell Atlas (e.g., Mourabit et al., 2014; Leprêtre et al., 2018). Some relevant events in the Alboran Sea are the Adra earthquake in 1910 (M_w6.1; Stich et al., 2003), and the Al Hoceima earthquakes of 1994 (M_w6.0), 2004 (M_w6.3) and 2016 (M_w6.4) (Kariche et al., 2018).

The majority of recorded earthquakes in the IMR are between 0 and 30 km, although a considerable number of earthquakes are of intermediate depth ($\mathrm{30}\le h\le \mathrm{180}$ km), mostly located in the Gulf of Cádiz (IGN, 2025c). Deeper earthquakes (h≥600 km) are less frequent and located in the province of Granada (Buforn and Udías, 2007). About 90 % of the registered seismicity in the IMR is located in the Gorringe Bank, High Atlas and the Granada-Málaga-W of Alboran Sea area. In the Alboran Sea, there is also significant shallow seismicity at depths of less than 30 km, especially in the active Betico-Alboran-Rif shear zone; to the W, earthquakes go deeper.

3 Materials

For seismic hazard estimation on the seabed, using the DSHA and PSHA methods, the following input parameters are required: a digital terrain model (DTM); a seismogenic sources and their associated seismic parameters; and GMPEs.

3.1 Ocean and land digital terrain models

The DTM used in this study was obtained from the EMODnet Bathymetry World Base Layer (EBWBL), composed of the bathymetric grid for Europe from the European Marine Observation and Data Network (EMODnet; https://emodnet.ec.europa.eu/en/bathymetry, last access: 24 August 2026) (Martín Míguez et al., 2019), with spatial resolution of 4 arcsec, and the raster for oceans and land from the General Bathymetric Chart of the Oceans (GEBCO; https://www.gebco.net/, last access: 24 August 2026, GEBCO Bathymetric Compilation Group 2024, 2024; Mayer et al., 2018), with a spatial resolution of 15 arcsec.

3.2 Seismogenic source models and seismicity parameters

Considering the SHAs carried out in the emerged IMR zone, several seismogenic source models have been identified (Table 2). The one implemented in the ESHM20 (Danciu et al., 2021a, 2024) has been used for this work given it is the latest update. This zonation includes areal seismogenic sources (shallow and deep), active faults and background smoothed seismicity (Fig. 2). To achieve our SHA, all seismic sources influencing our study area were incorporated, together with their corresponding characteristic seismic parameters.

Table 2 Some seismogenic zonation models covering the study area.

Figure 2 Seismogenic source models of ESHM20: (A) areal seismogenic sources (shallow and deep); (B) active faults and (C) smoothed seismicity; in function of (1) Depth (km), (2) Predominant rupture mechanism for each seismogenic source, and (3) Magnitude of maximum credible earthquake (MCE) (Data from Basili et al., 2024; Danciu et al., 2021b, 2024).

The seismic parameters used for the DSHA (Fig. 2), were extracted from the ESHM20 model (Danciu et al., 2021a, b). The maximum magnitude < MM_max> from the maximum magnitude distribution function, was adopted as the maximum credible earthquake (MCE) along with the predominant rupture mechanism assigned to each planar source ID at its depth (Tables S1 and S2 in Supplement). The areal seismogenic sources are located at depths ranging from 12 to 60 km (Fig. 2A1). Seismic features derived from Active Shallow Crustal Region (ASCR) seismogenic sources have magnitudes ranging from 6.6 to 8.8 (Fig. 2A3). The predominant rupture mechanism is inverse in the band from the Gulf of Cádiz towards the Algerian Mediterranean coast, through the southern Rif and parallel to the Atlas chain (Fig. 2A2, B2 and C2). This rupture mechanism transitions into strike-slip toward the Alboran Sea and the adjacent Iberian coastal zone, while evolving into normal faulting further inland within the Iberian Peninsula. The deep source beneath the Strait of Gibraltar (∼150 km) was initially discarded in line with the original DSHA approach, since it is outperformed by shallower ASCR sources. However, this work includes the distribution of all seismic sources. The active faults are located at depths ranging from 6 to 25 km (Fig. 2B1) and magnitudes ranging from 6.9 to 8.5 (Fig. 2B3). For the background smoothed seismicity, point sources have depths ranging from 12 to 60 km (Fig. 2C1) and magnitudes ranging from 6.5 to 8.6 (Fig. 2C3). For the PSHA, the same dataset used in the ESHM20 model has been used (Danciu et al., 2021b, c). These files include all the seismogenic sources and their characteristic parameters, required to carry out the probabilistic calculation in our study area. Since we adopted the same calculation software (OQ Engine) used in the ESHM20 model for the seismic hazard calculation, it was not necessary to modify the seismogenic source files.

3.3 Ground motion prediction equations (GMPEs)

GMPEs are essential in seismic hazard analysis, as they model how ground motion attenuates with distance. These equations estimate shaking intensity at a site based on factors like earthquake magnitude, distance to the source, faulting style, and local geological conditions (Barani et al., 2016). GMPEs are commonly derived from strong ground-motion records for an earthquake (e.g., Campbell and Bozorgnia, 2003). Functions relate the magnitude M and the distance R of the seismic scenario (M, R). With an appropriate set of quality records, it is possible to fit a parametric model f(M, R, ϕ_k) to estimate the intensity of strong motion (Y) for a given seismic scenario. The fitted parameters reflect the characteristics of elastic and anelastic energy attenuation through distance scaling, but also near and far-source, faulting mechanism, and hanging wall effects. The heterogeneous geological media causes a dispersion in the (M, R) sample pairs. Therefore, the correlation between variables can be strong, but never perfect, meaning there is an unavoidable uncertainty ε. If the expected seismic intensity at the project site Y is interpreted as a conditional random variable on the pair (M, R), the fitted model (Al Atik et al., 2010):

$$\begin{array}{}\text{(1)}& Y\sim f\left(M,R,{\mathit{\varphi }}_{\mathrm{k}}\right)+\mathit{\epsilon }\end{array}$$

provide an estimate of the median ground motion and the residual ε identified by a Gaussian distributed random variable.

GMPEs are commonly defined in terms of relative fault distance to the SHA site assuming a planar rupture geometry. The rupture trace is defined as the projection of the upper edge of the rupture onto the ground surface. The rupture plane, trace, and surface projection allow the definition of four finite fault distance metrics: R_RUP, R_JB, R_X, R_Y, Z_TOR (Fig. 3). R_RUP is the closest distance between the site and the rupture surface. R_JB is the Joyner and Boore distance or the closest distance between the site and the surface projection of the rupture. R_X is the horizontal and R_Y the vertical distance between the site and the top edge of the rupture. Finally, Z_TOR is the depth to the top edge of the rupture. These metrics are usually incorporated in GMPEs, as R in the mathematical formulation of Eq. (1) (e.g., Douglas, 2024). Some of the GMPEs formulation Eq. (1) consider the distance R as a function of depth (e.g., R_RUP). In these cases, the estimation of the intensity of the motion is affected by the terrain relief where the site is located. Thus, the estimated shake motion Y in elevated areas will be lower than in marine abyssal plains because they are closer to seismogenic sources, and the results of Y estimation in regions with significant unevenness will be sensitive to the use of a DTM.

The wealth of available strong ground motion data has enabled the development of numerous GMPEs, the vast majority on the rocky ground inland (e.g., Douglas, 2024). However, in offshore regions, traditional rock-based models may not adequately capture the unique conditions of the seafloor (Chen et al., 2024a; Tan and Hu, 2023; Wang et al., 2024). The growing importance of subsea structures and the associated risks of seabed instability have prompted the development of efforts to obtain GMPEs for offshore areas. Nevertheless, the limited availability of data and the scarcity of strong motion records have restricted their collection to very few offshore areas. A significant proportion of recent research has focused on the active area of the Japan Trench using the S-net and KiK-net stations data (Kanazawa et al., 2016), where several GMPEs have been recently developed (e.g., Dhakal et al., 2024; Hu et al., 2020, 2023; Nakanishi and Takemura, 2024).

Figure 3 Earthquake source and distance metrics from an offshore evaluation site in and out a seismogenic source area. Site A outside the ground projection of the seismogenic source, with (H≠0) and without surface altimetry or bathymetry (H=0). Site B inside the ground projection of the seismogenic source, without surface altimetry or bathymetry (H=0).

For soft soils, a GMPE model was already described by Idriss (1991) in the early 1990s (IDR91). The IDR91 model is suggested for the entire moment magnitude range and a distance of 0–100 km or more. The Joyner-Boore metric is used as the source-to-site distance (R_JB) to obtain Ln (Y), where Y is in g, and the uncertainty (ϵ) in Ln (Y) depends on the earthquake magnitude in the function described by Idriss (1991) if M_w < 7.25, otherwise the standard error term is constant (ϵ=0.38). This model is suited for crustal earthquakes and it has been implemented in SHA for submarine zones in the Mediterranean Sea (e.g., D'Acremont et al., 2022; Lafuerza et al., 2012; Sultan et al., 2007). Therefore, IDR91 was also considered in this work to compare it with other recent approaches to offshore GMPEs. Recently, the HZ23 model (Hu et al., 2023), has been developed for offshore earthquakes in the Japan Trench area and derived from onshore model ZH06 (Zhao et al., 2006). HZ23 was developed for M_w between 4.0 and 7.1, and hypocentral distances of 20–300 km. This model uses R_RUP as distance metric in Eq. (1) for crustal, interface, and slab events to obtain Ln (Y) with a constant standard error (ϵ=0.7784), where Y is in cm s⁻². Additionally, the averaged shear-wave velocity to 30 m depth (Vs30) was incorporated as a parameter in Eq. (1) for the definition of the site class. The NT24 model (Nakanishi and Takemura, 2024), was also developed for the Japan Trench area and derived from the onshore model MF13 (Morikawa and Fujiwara, 2013). It can assess PGAs for M_w between 5.5 and 7.4 for crustal, intraplate, and slab earthquakes (M_w>5.0 for crustal events) and epicenter distances of 1 to 300 km. This model uses R_RUP to obtain log₁₀ (Y) with a constant standard error (ϵ=0.326). Finally, the DKK24 model (Dhakal et al., 2024) was developed for M_w between 5.5 and 7.4, and hypocentral distances of 1 to 300 km. This model uses R_RUP as distance metric in Eq. (1) for crustal, interplate and intraslab events to obtain log₁₀ (Y) with a constant standard error (ϵ=0.25183). Therefore, four GMPEs have been used in this work: IDR91, HZ23, NT24, and DKK24 (Fig. 4), considering both the seminal studies on soft soils and recent works proposing new GMPEs for submarine environments. In this study, the PGA (g) is the Y seismic intensity used in all these selected attenuation functions.

The predominantly offshore setting of the study area (at the intersection of oceanic, active continental, and stable continental crust) poses significant challenges for selecting appropriate GMPEs. Recent tectonic classifications (Chen et al., 2018; Poggi et al., 2020) describe the study region as a Variscan-Hercynian zone in an Active Shallow Crustal Region (ASCR). In the case of the HZ23, NT24 and DKK24 models and their implementation in our study area, the IMR and the Japan Trench area share some geo-tectonical similarities. In the IMR, the African Plate is subducting beneath the Eurasian Plate (Gutscher et al., 2002), while in the Japan Trench area, the Pacific Plate is subducting beneath the Okhotsk Plate (Nishikawa et al., 2023). Both areas experience significant seismic activity due to the tectonic interactions. The IMR has a history of large earthquakes (Buforn et al., 2015), while the Japan Trench area is known for its frequent and intense earthquakes, including the devastating 2011 T$\overline{o}$hoku earthquake (Ritsema et al., 2012).

The selected GMPEs were benchmarked against horizontal PGAs for six regional offshore earthquakes in the IMR (Fig. 4), including five crustal and one intraslab event, using onshore/near-coastal accelerograms. This benchmark is a statistical plausibility check within overlapping M_w-R ranges, but it does not constitute a validation nor recalibration of seabed GMPEs in the IMR. The primary evidence supporting seabed GMPEs (HZ23, NT24, DKK24) remains their calibration on ocean-bottom data (S-net/KiK-net) and the documented land-seafloor differences reported in the literature (e.g., Chen et al., 2015, 2024b; Dhakal et al., 2021; Hu et al., 2023). This approach follows previous IMR hazard studies (e.g., Molina Palacios, 1998; Olaiz et al., 2025; Terrinha et al., 2009; Vilanova and Fonseca, 2004), which also relied on inland proxies due to the absence of seabed recordings. Therefore, this comparison should be interpreted as an acceptability check within the GMPEs magnitude-distance validity ranges, rather than a direct GMPE re-calibration. The scarcity of offshore strong-motion data underscores the need for OBS deployments to enable future seabed-specific GMPE validation in the IMR. The comparative analysis of IDR91, HZ23, NT24, and DKK24 was conducted using the following events:

1. M_w 4.5 Alboran Sea earthquake, 31 January 2016 (Fig. 4A),

2. M_w 4.9 Alboran Sea earthquake, 4 January 1994 (Fig. 4B),

3. M_w 5.5 Oran earthquake, 6 June 2008 (Benfedda et al., 2020) (Fig. 4C),

4. M_w 6.1 SW Cape St. Vincent earthquake, 12 February 2007 (Fig. 4D)

5. M_w 6.3 Alboran Sea strike-slip earthquake, 25 January 2016 (Fig. 4E),

6. M_w 6.8 Boumerdès earthquake, 21 May 2003 (Khellafi et al., 2013) (Fig. 4F).

The acceleration data for earthquakes in Fig. 4A, B, D, and E are sourced from IGN (https://www.ign.es/web/sis-catalogo-acelerogramas, last access: 8 July 2025).

Figure 4 Relationship plots between PGA (g) and hypocentral distances (R_Hyp) of the selected GMPEs. Solid lines denote the median prediction curves for each earthquake, whereas the shadowed region between the dashed lines denote one standard deviation range (±s_ε).

The predicted PGAs, through the HZ23, NT24 and DKK24 models and their respective standard deviation bands around each curve are consistent with the observed ones from the six tested earthquakes (Fig. 4). Most of the observed PGAs are within ±s_ε bands. However, IDR91 behaves as an upper limit boundary of models in the majority of cases, and, for M_w6.3 Alboran Sea earthquake (Fig. 4E), the observed PGAs are usually lower than the predicted values of the GMPEs for seabed. To evaluate the model's accuracy, residuals were obtained as the difference between the observed PGAs and the estimated PGA from the selected GMPE. As Fig. 5 shows, a near-zero tendency of residual distribution is found, with mean residual values presented in Table S3. Despite the oscillations discernible in the graphs, residuals are situated between the curves (dashed lines in the graphs) that delineate the uncertainty bands (±s_ε) (Fig. 5A, C, E, G). The residuals demonstrate bounded variance homoscedasticity, with RMSE values of about 0.077, 0.236, 0.047, and 0.040 g for IDR91, HZ23, NT24, and DKK24, respectively (Table S3), with the M_w6.8 earthquake exception. All four GMPEs benchmarked against nearby observed ground motions demonstrate adequate goodness of fit. However, it is the NT24 and DKK24 models (Fig. 5E–H) that most effectively capture the variability of the observed values, even at the extremes of the distribution of the set of observations. Specifically, IDR91 (Fig. 5A, B) emerges as the most conservative GMPE in terms of providing seismic intensity values that exceed the observed ones, with the exception of the closest recorded PGAs. The GMPE that exhibits a reduced capacity to capture this variability is HZ23 (see Fig. 5C, D), with the highest RMSE relative to the others GMPEs (Table S3), but predicting a higher seismic intensity, according to the observed residuals, primarily for the largest earthquake (M_w6.8). Its prediction is, in this instance, limited by the lower boundary curve of the uncertainty band of one σ_s.

The final step in the diagnosis undertaken here evaluates the properties of the prediction errors in the logarithms of the PGA (Fig. S1), as defined in Eq. (1). The averages of the epsilon values should be zero, where NT24 and DKK24 best fit the theoretical value (Fig. 6A). Their standard deviations are either below or very close to the values that the authors of the GMPEs provide for their models. Furthermore, it cannot be discarded that residuals in the logPGA have a normal distribution, as shown by the histogram for and the Q-Q plot (Fig. 6B, C) and indicated by the p-value for the Kolmogorov-Smirnov and the Shapiro-Wilk tests in NT24 (bold p-values in Table S3) at 95 % confidence.

Figure 5 Residual plots of the goodness of fit for the predicted values by the selected GMPEs: IDR91 (A, B), HZ23 (C, D), NT24 (E, F), and DKK24 (G, H), against the registered data in the earthquakes used for diagnosis.

4 Seismic hazard assessment

Seismic hazard refers to the probability of exceeding a certain earthquake intensity in a specific area within a given time period (Kijko, 2019; McGuire, 2004). Seismic hazard can be assessed using two main methods: Deterministic Seismic Hazard Analysis (DSHA) and Probabilistic Seismic Hazard Analysis (PSHA). In DSHA, hazard is defined as a specific percentile (such as the 50th, 84th, or 98th, as used in this study) of ground shaking intensity measured as peak ground acceleration (PGA) caused by one or several earthquakes. This method uses basic earthquake data and ground motion models based on a chosen earthquake scenario (Krinitzsky, 1995, 2002). Despite their deterministic character, the results obtained have frequently been used as a proxy for the upper limit of seismic intensity. Seismic hazard in PSHA is determined by calculating ground motion frequencies or exceedance rates using a mathematical model based on earthquake statistics and ground motion functions (McGuire, 2004, 2008). This model takes into account uncertainties related to earthquake size and moment (Kramer, 1996). PSHA methods can be categorized into parametric, using the total probability theorem (Cornell, 1968; Kramer, 1996), and nonparametric methods that utilize extreme value distribution functions (Epstein and Lomnitz, 1966). Other approaches, such as NDSHA, are often discussed together with DSHA and PSHA in relation to their representativeness for the largest events. Although both SHA methods use seismological and geological information, they define and calculate seismic hazard differently, which improves the understanding of seismic hazard forecasting at the studied site or region (Orozova and Suhadolc, 1999; Wang and Cobb, 2012).

4.1 Deterministic Seismic Hazard Assessment (DSHA)

In DSHA, the capable controlling source at a site is identified, typically a mapped fault segment and, where faults are diffuse, an areal or background sources. The maximum credible earthquake (MCE) for that source and an appropriate GMPE are then used to compute a chosen percentile (e.g., 50th or 84th) of ground motion at the shortest source-to-site distance under the adopted distance metric (Krinitzsky, 2002; Reiter, 1990). DSHA does not explicitly incorporate probability and recurrence functions, the seismic hazard is obtained as a certain percentile of ground motion. The 50th percentile is mostly used in preliminary DSHA studies, as well as the 84th percentile when it is used for inland critical structures such as nuclear or hydraulic power plants. Thus, some DSHA use the second quartile PGA plus a standard deviation, which is equivalent to an exceedance probability of 16 % (Ben-Zion et al., 2003), calculated from simple earthquake and ground motion statistics (Krinitzsky, 1995, 2002). To perform these calculations, this study adopts a DSHA approach, following the standard algorithmic workflow as used in Wang et al. (2012), Vipin (2013), Candia et al. (2019), and Ramkrishnan et al. (2021) which implement the classical deterministic hazard assessment procedure (Reiter, 1990). Deterministic ground motions are computed for areal seismic sources in regions lacking well-defined fault geometries. This method has precedent in offshore hazard studies and is particularly suitable for submarine environments where seismicity is diffuse and fault structures are unresolved. Similar approaches have been discussed (i.e., Krinitzsky, 2002; Kramer, 1996) and are widely used in engineering applications where deterministic estimates are required for preliminary infrastructure design.

Computationally, the DSHA is implemented through the following steps:

•  

  (i) Build a catalogue of seismogenic sources (faults and, where needed, areal/background zones) {SS_j} ${}_{j=\mathrm{1}\mathrm{\dots }{N}_{\mathrm{s}}}$.

•  

  (ii) Assign the seismic potential to each source SS_j with its MCE_j and the prevailing focal mechanism (normal, reverse or strike-slip).

•  

  (iii) Select the set {fⁱ} of the empirical GMPEs or attenuation relationships Yⁱ=fⁱ (M, R, ${\mathrm{\Phi }}_{\mathrm{k}}^{\mathrm{i}}$, ${\mathit{\sigma }}_{\mathrm{Y}}^i$), with their corresponding parameters ${\mathrm{\Phi }}_{\mathrm{k}}^i$ and the random uncertainties ${\mathit{\sigma }}_{\mathrm{Y}}^i$ for each i-th prediction equation.

•  

  (iv) Arbitrarily select the desired probability of exceedance PR_exed= Pr[$Y\mathit{>=}{y}_{\max }/M,R$], for a seismic scenario given by the (M,R) pair.

•  

  (v) Calculate the p-th percentile equivalent to the probability of exceedance:

  $$\begin{array}{}\text{(2)}& {\mathrm{PC}}_{\mathrm{Pr}}=\mathrm{Pr}[Y<y_p/M,R]=\mathrm{1}-\mathrm{Pr}[Y>=y_p/M,R]\end{array}$$
•  

  (vi) Calculate the standard normal random variable Z₀₁, of mean 0 and variance 1, that matches the percentile PC_Pr.

•  

  (vii) Loop for each site P located at the geographic coordinate position ($x,y,h$) with latitude, longitude and hypsometry (height or depth above mean sea level) do,

•  

  (viii.a) Loop for each j-th seismic source SS_j (j=1, .... N_s), assuming that the worst-case scenario (M,R, focal mechanism) is selected, defined as the occurrence of an earthquake of magnitude M= MCEj from a point of the j-th seismogenic source at the shortest possible distance R=R_min =min{d(P,SS_j)}, based on the distance to be handled (R_JB, R_Epi, R_HYP, etc.) in the attenuation function (Fig. 3).

•  

  (viii.a.1) Calculate the median ground motion (50th percentile or GMPE treated as deterministic) of the seismic parameter $Y_{\mathrm{50},j}^i$ at the site, with each i-th GMPE: $Y_{\mathrm{50},j}^i$ = fⁱ($M,R,{\mathrm{\Phi }}_{\mathrm{k}}^i$).

•  

  (viii.a.2) Calculate the p-th percentile of the seismic parameter at the site with the i-th equation of motion:

  $$\begin{array}{}\text{(3)}& \log Y_{p,j}^i=\log Y_{\mathrm{50},j}^i+Z_{\mathrm{01}}{\mathit{\sigma }}_Y^i\end{array}$$
•  

  (viii.a.3) Obtain the on-site seismic ground shaking parameter Y_p,j on site P produced by each j-th earthquake source.

•  

  (viii.b) Deterministic evaluation of the th-percentile of the seismic hazard $Y_p(x,y,h)$ at the site located at ($x,y,h$) as the largest parameter of the intensity of ground shake obtained from each seismogenic source: $Y_p(x,y,h)=$ max{Y_p,j}

•  

  (ix) Write Y_p($x,y,h$) output to change to a new location of the P site to assess.

However, since in each GMPE the logY value is distributed as a normal random variable of mean logY with standard deviation σ (Kramer, 1996), and $Y_p(x,y,h)=$ max{Y_p,j}, both Y_p and Y_p,j are random variables. The p-th percentile Y_p is accurately calculated as follows:

$$\begin{array}{}\text{(4)}& \mathrm{Pr}\left[Y_p<={\displaystyle \frac{y_p}{M,R}}\right]=F\left({\displaystyle \frac{\log Y_p-\log y_{\max /M,R}}{\mathit{\sigma }}}\right)={\displaystyle \frac{p}{\mathrm{100}}}\end{array}$$

in accordance with the classical scheme, with a single dominant seismic source, and where F denotes the cumulative density function of the standard normal distribution (e.g., mean μ=0, and variance =1). Now, considering the distribution of $Y_{p,{\mathrm{N}}_{\mathrm{s}}}$ as the extreme value of a set of values (Ang and Tang, 2007; Coles, 2001), the above approximation is computed as:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathrm{Pr}\left[Y_{p,N_{\mathrm{s}}}<={\displaystyle \frac{y_p}{M,R}}\right]=\mathrm{Pr}\left[Y_{\mathrm{1}}<={\displaystyle \frac{y_p}{M,R}}\right]\\ & \times \mathrm{Pr}\left[Y_{\mathrm{2}}<={\displaystyle \frac{y_p}{M,R}}\right]\times \mathrm{\dots }\times \mathrm{Pr}\left[Y_{N_{\mathrm{s}}}\le {\displaystyle \frac{y_p}{M,R}}\right]={\displaystyle \frac{p}{\mathrm{100}}}\end{array}\end{array}$$

This implies incorporation of the effect of the remaining N_s seismic sources in the hazard percentile estimation in the hazard assessment. The estimation is more complicated than in the case of a single source, as it is now necessary to solve the nonlinear equation:

$$\begin{array}{}\text{(6)}& \prod _{j=\mathrm{1}}^{N_{\mathrm{s}}}F\left({\displaystyle \frac{\log Y_{p,N_{\mathrm{s}}}-\log y_{j/M,R}}{\mathit{\sigma }}}\right)={\displaystyle \frac{p}{\mathrm{100}}}\end{array}$$

As this paper uses an approach based on a logic tree scheme that weights the GMPEs to capture the epistemic uncertainties in the hazard estimation, the above equation is rewritten to solve for each i-th GMPE fⁱ ($M,R,{\mathit{\sigma }}_{\mathrm{k}}^i$) used and its corresponding i-th random uncertainty σⁱ:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\phi }}^i\left(\log Y_{p,N_{\mathrm{s}}}^i\right)& =\prod _{j=\mathrm{1}}^{N_{\mathrm{s}}}F\left({\displaystyle \frac{\log Y_{p,N_{\mathrm{s}}}^i-\log y_{j/M,R}}{{\mathit{\sigma }}^i}}\right)\\ & -{\displaystyle \frac{p}{\mathrm{100}}}=\mathrm{0}\end{array}\end{array}$$

from which each $Y_{p,N_{\mathrm{s}}}^i$ is obtained.

Figure 6 Diagnosis of the residuals in logarithms of PGA. (A) Box-Whisker plot of residuals for each GMPE model used in this study. (B) Stacked histogram in relative frequencies for the residuals in logarithms of PGA. (C) Q-Q plot of the standardised residuals in logarithms of PGA.

4.2 Probabilistic Seismic Hazard Assessment (PSHA)

The principles of PSHA find their roots in the work of Cornell (1968), Cornell et al. (1971), and McGuire (1976). Cornell uses the total probability theorem (Kramer, 1996) to consolidate information on various variables (magnitude, distance, source location) into a single measure of seismic hazard for a specific site. Once probability relationships are derived from the observed seismicity (locations, sizes, recurrence) in the area, exceedance probabilities of ground motion can be quantified (Kijko, 2019). It takes into consideration the seismicity of the area from the distribution of locations, sizes, recurrence times and energy attenuation (Benito and Jiménez, 1999), to integrate the probabilities of exceeding a certain level of ground motion from various seismic sources, at different distances and with different magnitudes, that may occur at a site.The basis for carrying out PSHA is described in several steps (Kramer, 1996; McGuire, 2004; Reiter, 1990):

1. identification and characterization of seismic sources;

2. determination of the recurrence distribution of each source, in time or magnitude;

3. establishing the response parameter based on magnitude, distance, and site conditions with its uncertainty; and

4. calculating the probability that the response parameter is exceeded over a time period and estimating the hazard at the studied site, resulting from the sum of the hazard probabilities of all seismic sources influencing it.

This method allows for obtaining the mean annual exceedance rate (λ_IM>x) for a given site and exposure time, defined as the probability that a specified level x of intensity measure (IM) of ground motion is exceeded at this site (McGuire, 2004):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathit{\lambda }}_{\mathrm{IM}>x}& =\sum _{i=\mathrm{1}}^{N_{\mathrm{s}}}{\mathrm{v}}_i\underset{r=\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}.\underset{m=M_{\min }}{\overset{M_{\max }i}{\int }}P\left[\mathrm{IM}>x\right|m,r]f_{\mathrm{Mi}}\left(m\right)\\ & f_{\mathrm{Ri}}\left(r\right)\mathrm{d}m\mathrm{d}r\end{array}\end{array}$$

Where:

• λ_IM>x is the annual exceedance rate of a soft soil movement intensity IM that exceeds a threshold value x.

• $P\left[\mathrm{IM}>x\right|m,r]$ is the probability that IM exceeds x given an earthquake of magnitude M=m at a distance R=r.

• f_Mi(m) is the probability distribution of the earthquake magnitude.

• f_Ri(r) is the probability distribution of the distance to the site.

• N_s is the total number of seismic sources considered.

And the mean annual occurrence rate of earthquakes from source i:

$$\begin{array}{}\text{(9)}& v_i=\exp \left({\propto }_i-{\mathit{\beta }}_i{m}_{\mathrm{0}}\right)\end{array}$$

with α=aln(10) and β=bln(10) being the parameters of the Gutenberg-Richter recurrence law:

$$\begin{array}{}\text{(10)}& {\mathit{\lambda }}_{\mathrm{m}}=\exp \left(\mathit{\alpha }-\mathit{\beta }m\right)={\mathrm{10}}^{a-bm},\end{array}$$

adjusted from the seismic catalog of the project using statistical technique of maximum likelihood. The probability function in Eq. (8) is:

$$\begin{array}{}\text{(11)}& P\left[\mathrm{IM}>x\right|m,r]=\mathrm{1}-F_{\mathrm{IM}}\left(x\right)=\mathrm{1}-\mathrm{\Phi }{}^{\prime }\left({\displaystyle \frac{\ln x-{\mathrm{lnIM}}_{\mathrm{average}}}{{\mathit{\sigma }}_{{\ln }_{\mathrm{IM}}}}}\right)\end{array}$$

where the function $P[\mathrm{IM}>x|m,r]$ indicates the IM ground motion parameter generated by a source i of magnitude (m), distance (r), and uncertainty epsilon (ϵ) in relation to the threshold value x. Finally, ƒ_Mi(m) and ƒ_Ri(r) are the parametric probability density functions for which a functional model is assumed, estimated from the magnitude and distance catalogue.

Results are expressed as seismic hazard curves (SHC), which can be obtained for different intensity measures, the most common being Peak Ground Acceleration (PGA), representing the exceedance probabilities over usually 50 years, and the ordinates of the acceleration spectrum S_a (T). Another result is the response spectra, commonly referred to as uniform hazard spectra (UHS), derived from the SHC and obtained for each chosen return period.

4.3 Epistemic uncertainty treatment: logic tree scheme in DSHA and PSHA

There are numerous uncertainties in the DSHA and PSHA methodologies that arise from limited comprehension of the seismic process. The random nature of earthquakes and their occurrence generate variations in the spatial, temporal, and magnitude distribution when defining and characterizing seismic sources. In addition, the way seismic waves attenuate as they propagate from the source to the site of interest can be uncertain, generating uncertainties about the attenuation models or GMPEs, whose statistical adjustment of their coefficients incorporates an estimation error. Moreover, the arbitrary choice of GMPEs by users is itself an unsure process. All these uncertainties can be classified as random and epistemic (Kiureghian and Ditlevsen, 2009; McGuire and Shedlock, 1981).

In this work, the logic tree method has been used to incorporate and analyse the inherent uncertainties both for DSHA and PSHA analysis. Its formulation uses weighting factors, assigned to seismic source models and GMPEs based on their likelihood (Coppersmith and Youngs, 1986; EPRI, 1987). The logic tree is constructed from a series of branches connected through nodes, from which the computational process twigs according to a possibility or appropriate weighting factors {wⁱ}. These factors are assigned based on expert judgment (Budnitz et al., 1997), depending on the suitability and importance of each seismic source model and GMPE used. The treatment of uncertainties by the logic tree scheme is made mathematically explicit in the calculation process in the equations to obtain Y_p,j at site P produced by each j-th seismic source:

$$\begin{array}{}\text{(12)}& Y_{p,j}={\sum }_i{w}^i{Y}_{p,j}^i\end{array}$$

And also, to obtain Y_p,Ns at each site P:

$$\begin{array}{}\text{(13)}& Y_{p,N_{\mathrm{s}}}={\sum }_i{w}^i{Y}_{p,N_{\mathrm{s}}}^i\end{array}$$

which includes the effects of the N_s seismic sources. The logic tree scheme shown here is used both for the estimation of the PGA based on a single control source, with the most unfavourable seismic scenario, and for the PGA estimation based on all seismogenic sources with their particular seismic scenario.

The logic trees used in this study (Fig. 7) incorporate a combination of three seismic source models and four GMPEs, selected to estimate the PGA in our study area. In the logic trees used in the DSHA (Fig. 7a) and PSHA (Fig. 7b), the weights assigned to each seismic source model were those used in the ESHM20 (Danciu et al., 2024). The GMPEs weights in this study are based on statistical performance and regional relevance, in order to properly address epistemic uncertainty. At each node, branch weights are normalized so that $\sum w=\mathrm{1}$.

4.4 DSHA and PSHA calculation specifications

The DTM from EMODnet (Fig. 1) has been used to generate the computational grid in the seismic hazard calculations of this work. The DTM was resampled to generate grid points spaced at 0.1° (approximately the standard error in earthquake epicentre determination; Panza et al., 1990). In addition to its geographical coordinates, each point was assigned its respective bathymetric or topographic elevation. This mesh size, commonly used in the region (e.g., Danciu et al., 2024; IGN-UPM, 2017; Poggi et al., 2020), was chosen for its balance between computational efficiency and spatial resolution. This allowed capturing the variations in seismic hazard in this study while keeping a manageable computational load. Thus, computationally, SHA calculations, using DSHA and PSHA, were performed on a grid of 3854 sites (Fig. S2). The elevation of the DTM is used in the estimation of the relative geographical distances, from the site to the seismogenic source, used by the GMPEs.

Figure 7 Schematic representation of the logical tree structure of the seismic hazard models: (A) DSHA and (B) PSHA. The weights assigned for each branch are shown in red (modified from Danciu et al., 2024).

For DSHA calculations, a computational routine was developed to calculate the seismic intensity with a selected GMPE, from different seismic source types, incorporating the extreme value distribution into the DSHA as the scenario in which all seismic sources contribute (maximum Y distribution in Eq. 5), and to account for the logical tree with them. Since the extreme value distribution must be determined at each site where the hazard is estimated, the convergence characteristics of this algorithm allow it to be repeated for up to 100 iterations. This number is sufficient to achieve an admissible error, $|Xo-X*|$ in the root X* solution of Eq. (7) (Burden et al., 2015). To reduce the computational time compared with earlier DSHA implementations (e.g., Wang et al., 2012), the routine was optimized through parallel processing in MATLAB© (R2024a, Update 6, available at https://www.mathworks.com, last access: 8 July 2025). The routine is publicly available at https://es.mathworks.com/matlabcentral/fileexchange/183034-dsha (last access: 24 August 2026).

PSHA calculations were performed for the mean values of PGA (in g) and four annual probabilities of exceedance (APEs), i.e., 0.01, 0.002105, 0.000404, and 0.0002, which correspond to the return periods of 100, 475, 2475, and 5000 years, respectively. These hazard maps are prescribed in the European seismic-resistant standards. Here, the OQ Engine configuration and the input files used are publicly available from ESHM20 (Danciu et al., 2021b, c), but the calculations were performed under soft soil conditions, with a Vs30 = 200 m s⁻¹, corresponding to class E according to the Eurocode 8 classification (CEN, 2004) and the National Earthquake Hazards Reduction Program (NEHRP) (BSSC, 2004). The DTM was incorporated into the calculations. For epistemic uncertainty representation in this study, a total of 100 000 logic trees end-branches were sampled using the OQ random technique. The programming of the new GMPEs, as well as all PSHA calculations in this study, was carried out using OpenQuake (OQ) engine (Pagani et al., 2014) version 3.21.0 in development mode (available at https://docs.openquake.org/oq-engine/manual/3.19/getting-started/installation-instructions/universal.html, last access: 8 July 2025).

5 Seismic hazard assessment results and discussion

The findings of this study demonstrate the effectiveness of using both DSHA and PSHA methodologies to generate seabed seismic hazard maps in the IMR. The spatial hazard patterns in these maps are mainly controlled by the seismogenic source model and its parameters. The selected offshore oriented GMPEs are statistically consistent with the onshore/near coastal data used for benchmarking. This advancement represents a significant improvement over previous results in this field. In particular, this study extends and updates the only published maps (i.e., Molina Palacios, 1998) in the IMR.

The seismic hazard maps presented in the results obtained with either DSHA or PSHA suggest that the highest PGA values are concentrated in deep seafloor areas, where areal seismogenic sources of medium and high activity and faults are found. It has been demonstrated in previous studies that PGA maps depict the areas of highest seismic intensity as being limited to coastal and inland regions of the Iberian Peninsula and northern Africa or zones characterized by high seismic activity along major tectonic fault line (e.g., Danciu et al., 2024; IGN-UPM, 2017; Poggi et al., 2020). However, the present results expand the geographical coverage of previous studies by highlighting the overlapping influence of the Gorringe and Horseshoe active zones with the SVF, MPF, and SWIM faults (distributed across the western IMR) and by emphasizing the significant role of the Arzew, AIF, and YSF faults in shaping the diffuse seismicity sources in the Alboran Sea. This is of great interest in the design of offshore infrastructures and their current and future deployment on the ocean floor, as vulnerable infrastructures may be affected by significant seismic effects not previously considered (Karthigeyan, 2022; Kee et al., 2023).

5.1 Probabilistic approach

The seismic hazard maps for the area between the western Gulf of Cádiz and the eastern Alboran Sea, with return periods of 100, 475, 2475, and 5000 years are shown in Figure 8. The expected PGA values vary according to the analysed return period, generating four different seismic hazard scenarios in the IMR. For the 100-year return period (Fig. 8A), the mean horizontal PGA values are between 0.17 and 0.60 g. For 475-year return period (Fig. 8B), these values are between 0.26 and 1.06 g. In the case of 2475-year (Fig. 8C), the PGA values range from 0.38 to 2.59 g. Finally, for a 5000-year return period (Fig. 8D), the PGA values are around 0.44–3 g, respectively. For each map shown, the highest PGA values are found in geographic areas over deep submarine zones of the IMR that exhibit the highest levels of seismic activity and intensity. Consequently, the highest seismic hazard values in the IMR are observed across most of the studied extent of the Gulf of Cádiz. This pattern is primarily attributed to deep seismic sources, located at approximately 60 km depth, which are associated with more intense reverse or strike-slip faulting ruptures. Notable examples include the areal sources PTAS390, ESAS129, and MASS344 within the Gorringe and Horseshoe regions. Shallower right-lateral (RL) and reverse (RR) faults, including the SVF, MPF, and SWIM, also contribute to the observed hazard levels. In the Alboran Sea, the highest PGA values are mainly organised in accordance with the most active shallow fault systems, Arzew, AIF and YSF, forming a pattern of strike-slip oriented regions. It is evident that these values exhibit a reduced sensitivity to the selected GMPEs, in comparison to those obtained in the western sector of the Gulf of Cádiz. The PGA values, attributable to the AIF surface faults, are located over a SW-NE band that traverses the Alboran Basin and extends northwards, accompanied by a band of more attenuated intensities NW-SE, as a consequence of the YSF faults' orientation. At the easternmost point of the IMR, the Arzew Faults, which are characterised by elevated levels of activity, delineate a high-hazard zone along the Algerian coast that extends in an east-west direction, reaching almost the centre of the Alboran Basin. This phenomenon is attributable to the geometrical configuration of these semi-deep faults.

Figure 8 Probabilistic seismic hazard maps for the mean horizontal PGA (g) at APEs of (A) $\mathrm{1}/\mathrm{100}$, (B) $\mathrm{1}/\mathrm{475}$, (C) $\mathrm{1}/\mathrm{2475}$, and (D) $\mathrm{1}/\mathrm{5000}$. All values are estimated under soft soil conditions with Vs30 = 200 m s⁻¹.

Regarding the locally high PGA values (approaching the upper bound of ∼3 g) that appear in the longest return period map (Fig. 8D) and, to a lesser extent, in the 2475-year return period map (Fig. 8C) should be interpreted with particular caution because long return periods are increasingly controlled by rare, large-magnitude scenarios and by the upper tail of the ground-motion variability. It is likely that these results stem from extrapolating the HZ23 model beyond its intended range of applicability. Although HZ23 contributes only 10 % of the total weighting scheme, it has comparatively large aleatory variability (ϵ=0.7784). Therefore, when the IMR hazard calculation requires evaluating scenarios that may exceed the GMPE calibration ranges, the HZ23 branch can disproportionately influence the high-end tail of the combined hazard estimate. This effect is amplified at very low annual exceedance probabilities (2475–5000 years), where extrapolation beyond validated magnitude ranges and model variability have greater leverage on the final PGA levels. Consequently, the highest PGA patches in Figs. 8C–D should be regarded as reflecting increased epistemic uncertainty associated with GMPEs extrapolation and variability. The overprediction observed in the results of the probabilistic seismic hazard maps suggests that, while the HZ23 model may offer valuable insights within its original context, its use in the IMR should be approached with caution. Although the hazard results remain consistent with the prevailing seismicity of the study region, the highest PGA values obtained for the 2475- and 5000-year return periods should be interpreted with caution. Because these hazard levels are increasingly controlled by rare large-magnitude scenarios and by extrapolation beyond the calibration ranges of some of the GMPEs, a portion of the largest PGA estimates may reflect methodological artefacts associated with the upper tail of the hazard calculation, rather than fully constrained physical ground motion levels. This limitation should be taken into account, and priority should be given to the future development of regionally calibrated offshore GMPEs for the tectonic context of the IMR to reduce uncertainty in the large magnitude range.

Results show significant differences compared to previous PSHA studies in the IMR (e.g., Molina Palacios, 1998), particularly in marine zones. These results not only refine previous assessments in IMR but also improve the understanding of the tectonic influence on offshore seismicity. More recent SHAs typically highlight the higher seismic intensity in the coastal areas of the Iberian Peninsula and North Africa (e.g., Danciu et al., 2024; IGN-UPM, 2017; Poggi et al., 2020). However, these new results extend the geographical scope of seismic hazard studies in the region by demonstrating considerable seismic activity in the marine areas of the IMR. This emphasizes the need to consider seabed seismic dynamics in regional SHA studies.

The resulting seismic hazard maps could be highly valuable for offshore infrastructure design, as they highlight the potential vulnerability of current and future subsea facilities in the region. Hazard results are presented in terms of seabed PGAs, which enables first-order screening of offshore shaking demand. Spectral acceleration Sa(T) hazard maps, which are more directly linked to structural response at different periods, are not provided because the selected offshore GMPEs that have been implemented are currently developed for PGA (and, in some cases, PGV) (e.g., Dhakal et al., 2024; Nakanishi and Takemura, 2024).

5.2 Deterministic approach

The obtained mean PGA (g) seismic hazard map (Fig. 9), considering the distribution of maximum control scenarios at each of the 3854 study sites (Fig. S2), was generated by incorporating the DTM into the calculation, as well as the specific rupture mode of each fault and seismogenic source (Tables S1 and S2) and the GMPEs considered in this study: IDR91, HZ23, NT24, and DKK24. As illustrated in the applied logic tree (Fig. 7A), the weighted implementation of each GMPE has resulted in a consolidated seismic hazard scenario, which has been weighted according to Eq. (13). The expected PGA values in the deterministic seismic hazard scenario are between 0.58 and 2.55 g. The spatial distribution of accelerations shows that the highest PGA values are found in the Gorringe seafloor and extend to the Horseshoe Plain. These abyssal areas are influenced by seismogenic sources and faults with moderate to strong seismicity, generating PGA outliers in the area, ranging approximately between 1.2 and 2.55 g. To the east, PGA values range from 0.65 to 1.8 g. This relates to a sequence of moderate seismicity sources along the Algerian coastal zone, together with the Arzew faults. Notably, the use of a DTM to determine the three-dimensional location of each calculation site, may result in amplified PGA values in areas where the site depth is closer to the seismic source, particularly in the deeper oceanic regions of the IMR. In general terms, the PGA (g) map obtained show a geographical distribution of the maximum values over the Atlantic areas of Gorringe, Horseshoe and Coral Patch, as well as off the Algerian coast, in the Mediterranean Basin of the Alboran Sea. This coincides with the areas of highest seismogenic potential in the IMR.

Figure 9 Deterministic seismic hazard map in terms of mean PGA (g) under soft soil conditions.

As Reiter (1990) points out, the DSHA does not contemplate the inherent uncertainty in the estimation of seismic hazard, since the frequency of occurrence and magnitudes are not explicitly considered. So, a less accurate version, distributing the uncertainty due to randomness in the attenuation functions over a large region, is preferable, even though it provides larger values as in this work. These PGA values can be used as preliminary upper bounds of those that can be achieved, for different exceedance probabilities, rather than working with more precise, but less accurate PGAs.

The decision for the SHA method is based on its suitability for assessing earthquake hazards, as earthquakes are inherently random and unpredictable. It is important to acknowledge that hazard estimates cannot be fully verified due to the nature of seismic events (Musson, 2012a, b; Wang, 2012). Therefore, the use of DSHA could be recommended for critical industrial facilities of socioeconomic relevance.

5.3 Distribution of PGAs amplification relative to a baseline PSHA using the ESHM20 model (rock, onshore).

Baseline PSHA seismic hazard maps were computed using the published ESHM20 source model and GMPE logic tree exactly as released (rock, onshore site conditions) (Danciu et al., 2021a, b, c, 2024), clipped to our study zone in the IMR, to provide a reference for how the hazard results obtained, by incorporating seabed-specific GMPEs and bathymetry, alter PGA values. These baseline results are not part of the official ESHM20 products and are used here strictly as a methodological contrast. The distribution of PGA values obtained by baseline PSHA (ESHM20 model; rock) is increased throughout all its range (Fig. 10A and B). Incorporating the relief of the ocean floor through bathymetry has an effect, although not as significant as the amplification of the attenuation curves. The deeper the seismic intensity increases, the closer the control point is to the seismogenic sources in relation to inshore territory or mean sea level.

Figure 10 Distributions and relative-change maps of PSHA (in terms of PGA) at the seabed. (A–B) Histograms (in relative frequency) and box–whisker plots of estimated PGA values for return periods of 475 years (A) and 2475 years (B), comparing baseline PSHA results under rock/onshore site conditions (PGArock; blue) with the seabed PSHA results obtained in this work (PGAseabed; green) within the study area in the IMR. Baseline PSHA (ESHM20 model; rock) results were computed by the authors using the published ESHM20 seismogenic source model and GMPE logic tree exactly as released, clipped to the study area. They are not official ESHM20 products and must not be interpreted as an offshore extension of ESHM20. (C–D) Relative-change maps are shown expressed as the ratio ΔPGA = PGAseabed/PGArock for return periods of 475 years (C) and 2475 years (D).

The combination of GMPEs, bathymetry and the spatial distribution of seismogenic sources in this study has provided a new geographical pattern of amplification rates PGA (see Figs. 10C and D) very similar for different return periods. As demonstrated in the literature (see Dhakal et al., 2021; Hu et al., 2020), offshore sites exhibit differing distance attenuation when compared to inshore sites. This discrepancy is due to the incorporation of modifications to offshore zones as specified by the authors of the GMPEs within the model. It is acknowledged that some of the offshore GMPEs used (HZ23 and NT24) are derived from onshore GMPEs. The maximum amplification reaches approximately about 7 to 9 times in deepest areas of the IMR eastern limits, all very specifically located in AOB. In the remaining zones, the rate of increase of the PGA is within the medium ranges in comparison to the reference PGA_rockvalues in rock and emerged terrain (PGA calculated using the ESHM20 seismic model with the PSHA approach). It is observed that, in most of the studied area, the amplification rates were found to be less than five times in all seismic hazard scenarios. The site effect is primarily attributable to the attenuation characteristics of soft seafloor soils and their depth found in WAB and AOB (Chen et al., 2023; Diao et al., 2014), which are incorporated through the GMPEs used.

Regarding the spatial distribution of PGAs amplification rates, the highest amplification rates (represented by yellow, orange and red colors), exceeding five to nine times max the PGA values from the baseline PSHA (ESHM20 model; rock), are concentrated in regions where the expected PGA values are relatively low, such as the Algero-Balearic Basin. Mid-range amplification rates (blue, green and olive colors) are distributed across the entire study area. In this case, regions where the PGA values are already high, such as the Gorringe Ridge, Horseshoe Abyssal Plain, and the Algerian margin, the amplification rates tend to be moderate, typically ranging between 2 and 4 times the PGA_rock. On the other hand, the lowest amplification rates (blue color) are observed in the West Alboran Basin (WAB), following the Alboran Ridge and Channel, as well as in the Gulf of Cádiz. As illustrated in these amplification maps, the phenomenon is most evident in the West Alboran WAB and the Algerian Ocean Basin (AOB), extending to the west in the Gorringe Ridge and Tagus Plain. In such instances, the site effect assumes a heightened significance, particularly in the context of basin bottoms characterised by the accumulation of substantial sedimentary deposits.

Within IMR geographical framework, no discernible spatial correlation is identified between the amplification outcomes resulting from the GMPEs used in this study. Areas exhibiting low seismic intensity, such as AOB, demonstrate a substantial amplification effect, while regions characterised by higher PGA values exhibit average amplifications, as observed in the Gorringe Ridge area, or even attenuations to the west of Gorringe Bank and in the AIF and YSF fields. Excluding these areas where the PGA values are highest or lowest, the study area shows an amplification rate of less than three times, despite the variations in the PGA values obtained with the GMPEs used and with the baseline ESHM20.

Due to the absence of data concerning strong motion events within the IMR, there exists inadequate empirical data within the magnitude-distance range to directly develop a GMPE using empirical regression techniques. The authors of this study argue that there is a significant need for underwater seismic observation stations. This need arises from the considerable uncertainty and limited understanding of attenuation processes within the IMR. This is regardless of the extensive knowledge already available on amplification and attenuation processes caused by the substantial water column and the rheological characteristics of the ocean floor. To address this gap, the authors recommend implementing underwater seismic stations similar to those used in the United States (Boore and Smith, 1999) and the S-Net network in Japan (Kanazawa, 2013).

This approach would enable two key advancements: (1) calibration and empirical modification of existing GMPE, as demonstrated in HZ23; (2) reformulation of a regionally adjustable backbone attenuation function, or a system of ergodic functions, to better capture ground motion variability. Such methods have already been successfully applied in seismotectonic settings, as shown in the ESHM20 model (Kotha et al., 2020; Weatherill and Cotton, 2020). These improvements would be tailored to the specific Vs30 and D14 distribution characteristics of the IMR as recently detailed by Hu et al. (2024).

6 Conclusions

This study addresses the critical need for a comprehensive seismic hazard assessment in the IMR, particularly in the marine domain between the Gulf of Cádiz and the Alboran Sea. Despite the well-documented historical and instrumental seismicity of the region, previous hazard studies have largely focused on terrestrial zones, leaving a considerable gap in our understanding of undersea seismic hazard, especially in terms of PGA at the seabed level. This work contributes to filling that gap by applying both probabilistic (PSHA) and deterministic (DSHA) methodologies, adapted to offshore conditions and integrated with recent developments in seismogenic zonation (ESHM20) and GMPEs. This study has demonstrated that PGA estimation for seabed conditions is both feasible and meaningful, providing crucial insight for hazard mitigation planning in marine environments.

Key findings show that significant PGA values are present across deep marine zones of the IMR, particularly over the Gorringe Ridge, Horseshoe Plain, and along active fault zones near the Algerian margin. The results confirm a consistent spatial distribution of hazard for both DSHA and PSHA approaches, with maximum values occurring over regions characterised by dense seismogenic activity and structural complexity. The incorporation of four GMPEs (IDR91, HZ23, NT24, and DKK24) allowed a comparative analysis of offshore attenuation behaviour. Among these, NT24 and DKK24 provided the best statistical fit to observed regional ground motion data, while HZ23 exhibited overestimation tendencies beyond its intended magnitude range, indicating potential constraints in its application to this region. The logic-tree framework adopted in both DSHA and PSHA enabled the epistemic uncertainty to be systematically incorporated into the hazard modelling. The overestimation observed with HZ23 highlights the importance of developing or calibrating GMPEs specifically for the IMR, particularly for large magnitude, short distance offshore events. Given the extrapolation beyond the validated model range, we recommend that future SHA applications in the IMR down-weight HZ23 in the logic tree or exclude it of the analyses to avoid potential overestimation of PGAs in the region.

The findings provide a regional-scale hazard characterization and a first-order screening framework for assessing potential seismic exposure of subsea infrastructures, such as pipelines, submarine telecommunication cables, offshore energy facilities, and potential transcontinental tunnels. Current hazard assessments based solely on terrestrial models may significantly underestimate seabed accelerations and associated risks in the IMR. Moreover, the results underscore the importance of incorporating seabed-specific GMPEs and digital terrain models in offshore SHA, as bathymetry demonstrably affects source-to-site distances and thus ground motion predictions. In this context, the PGA-based hazard maps provide first-order screening of seabed shaking for offshore infrastructures. However, because subsea assets exhibit a wide range of fundamental periods, a complete design-oriented hazard characterization ideally includes spectral acceleration Sa(T) at representative periods. The Sa(T) characterization remains a necessary next step for detailed structural design of subsea structures with longer natural periods. Future work will extend the framework to Sa(T) hazard and uniform hazard spectra once suitable offshore spectral GMPEs and/or regionally calibrated models for the IMR become available.

However, some limitations remain to be addressed. The scarcity of direct offshore strong-motion records and the use of a uniform Vs30 proxy to represent soft soils in the IMR required the adoption of GMPEs derived from analogue seismotectonic contexts, such as the Japan Trench. While this approach was statistically validated, the use of offshore GMPEs developed for the Japan Trench region represents one of the main sources of uncertainty in the study. Although these models represent the most suitable currently available formulations for offshore applications and capture several physical processes specific to seafloor ground motions, the tectonic analogy remains imperfect. The IMR is a complex plate boundary region characterized by diffuse deformation, whereas the Japan Trench represents a more clearly defined subduction setting. These uncertainties should be addressed in future studies through the expansion of OBS monitoring networks in the region and region-specific GMPE calibration. Similarly, the DTM resolution, though consistent with regional studies, may benefit from refinement in areas of steep bathymetric gradients to improve hazard detail, accordingly, to include possible site-effects. Moreover, the potential interaction between earthquakes and coseismic hazards (such as landslides and tsunamis) requires further research. This could be approached through integrated, multihazard assessment frameworks to better understand and manage these complex interdependencies.

The integration of Neo-deterministic Seismic Hazard Assessment (NDSHA) method and machine-learning approaches for source characterisation are promising avenues. Under the adopted modelling choices, both DSHA and PSHA yield consistent spatial patterns of high hazard over major offshore structures in the IMR.