The first Bouguer gravity map for the entire Mediterranean region, published by Makris et al. (1999), was a milestone.These data have been acquired starting from the beginning of the 60 s, and from today’s perspective, the numerical methods used for data processing did not meet modern methodological requirements. The old Bouguer gravity map by Makris et al. (1999) is shown in Figure 2A. It has a large scale, is of lower resolution and unsuitable for resolving short wavelength anomalies, which are important for modeling purposes.

The new gravity field compilations published by the “AlpArray Gravity Research Group, AAGRG” (Zahorec et al., 2021) on a 4 km × 4 km grid for both Bouguer and Free air anomalies provide a database suitable for 3D modeling (Figures 2B, 3). These compilations follow the most recent gravity processing standards. After removing long wavelength regional fields, the gravity map shows local structures valuable for modeling shallow density variations in the subsurface. The Bouguer anomalies are complete (CBA) and generated using the latest criteria and reference frames (both positional and gravity reference systems). Atmospheric corrections were also applied. Error statistics (Zahorec et al., 2021) based on cross-validations and interpolation residuals indicate high data accuracy. For example, the Austrian dataset shows interpolation residuals between −8 and +8 mGal, cross-validation residuals between −14 and +10 mGal, with standard deviations well below 1 mGal. The database accuracy is approximately ± 5 mGal for most areas. This compilation shows negative gravity values in the Alps and the Po Basin in northern Italy. The Ivrea high in the western Italian Alps and the dominant high in the Bouguer anomaly in the Liguro-Provençal basin are clearly visible. Detailed descriptions of individual anomalies can be found in Zahorec et al. (2021).

The Bouguer anomaly map (Figure 3A) features positive anomalies offshore and moderate to negative values on land. Bluish colors indicate strong negative gravity in the Apennines, the Po Basin, and the Southern Alps, reflecting crustal root density deficits. A significant advancement of this compilation is its resolution, showing short wavelengths even in the offshore region of the Liguro-Provençal Basin, Corsica, and northern Sardinia. The effects of the Ivrea body and previously unknown positive anomalies in the Liguro-Provençal remain in the residual field. In Corsica, the anomalies are significantly smaller than in the offshore Mediterranean Sea. Detailed information on calculating the CBA by filtering/smoothing is provided in the text.

For the gravity field interpretation presented in Section 6, it is still necessary to eliminate a regional field from the CBA. The following procedure was used: 1) subtraction of the XGM 2019e_2159 satellite field (ICGEM website, http://icgem.gfz-potsdam.de/tom_longtime), 2) slight spline smoothing to eliminate remaining small local artifacts. The resulting residual field at the end of this procedure is referred to here and later as the “residual field” (Figure 3B).

In addition to the new Bouguer map, a Free air anomaly map has been compiled by the AAGRG (Zahorec et al., 2021). This new map provides more details than previous compilations for the region (e.g., Sandwell et al., 2014).

Although most gravimetric models use Bouguer gravity fields as the reference, we chose free air anomalies (Figure 3C) to be compatible with the previous 3D density model for the Alps and their forelands - 3D ALPS (Spooner et al., 2019b; Spooner et al., 2019a), which uses Free air anomalies and was available for our studies. This approach allows the inclusion of rock-related density inhomogeneities in the crust between the surface and the model reference level. Our goal was to extend the 3D ALPS model to the offshore Ligurian Sea, ensuring compatibility with the existing density model of the Alps.

As expected, the image of the Alpine Free air anomalies is strongly influenced by the topography of the study area: gravity highs largely correlate with high topographic elevations (Western Alps, Apennines, Corsica, etc.), while extended lows indicate lower density in the Po basin of northern Italy. Deviations suggest density changes in the Earth’s crust and lithosphere. For the 3D model calculations, also the Free air anomaly of the AAGRG was smoothed, as a detailed 3D gravity model using the original Free air anomaly would have been too large in terms of model size and the number of gravity stations to process for interactive modeling (both storage and computation time).

Considerable interest has focused on studying the lithosphere of the western Mediterranean using seismic methods. An overview of campaigns in the Liguro-Provençal Basin up to 2002 is provided by Rollet et al. (2002). Recent surveys include the CROP deep seismic profiles (De Voogd et al., 1991; Finetti and Morelli, 1973), TGS-NOPEC and SARDINIA profiles (Gailler et al., 2009; Jolivet et al., 2015), the GROSMarin 3D seismic refraction project (Dessa et al., 2011), the LOBSTER campaign (Dannowski et al., 2020; Kopp et al., 2023), and the SEFASILS cruise (Dessa et al., 2020). These studies, along with other gravimetric and geological data, serve as important constraints for our gravity modeling. See Figure 4 for the location of seismic lines used to constrain our density model.

Other important constraints for our 3D density modeling come from earlier shipboard seismic measurements by Makris et al. (1998) and measurement campaigns by Morelli and colleagues (e.g., Morelli and Nicolich, 1990; Morelli et al., 1975), which provide insights into the lithosphere structure in the western Mediterranean. Due to the early stage of density modeling development at the time, the density models of these authors in the Ligurian Sea are relatively simple.

Figure 5 shows earthquake distribution in the Ligurian Basin, based on the USGS earthquake catalog (U.S. Geological Survey, 2021). Events between January 1900 and February 2021 with a minimum magnitude of Mw 2.5 are included, with the largest earthquakes reaching magnitudes up to Mw 6.5 at the Ligurian margin. Earthquakes are divided into shallower, 5–20 km deep (yellow/orange), and deeper than 50 km (dark violet) focal depths in Figure 5.

Within the Liguro-Provençal basin, no events deeper than 50 km were observed. On the Italian mainland, especially southeast of Genoa, many earthquakes occur at depths of 10–20 km due to the subduction of the Adriatic plate Doglioni et al. (1999). A cluster of earthquakes of magnitude 2.5 is located in the NW Liguro-Provençal Basin. A notable feature at 42°N and 7.0° to 7.5°E is the accumulation of earthquakes with epicenters in depths between 10 and 50 km (orange and pink). Thorwart et al. (2021) interpreted this cluster and the focal mechanisms as an indication of “re-activated pre-existing rifting-related” structures. This cluster corresponds spatially with steep gradients in the gravity field of the upper mantle, as shown in Figures 5, 15. Several amphibious ambient noise studies have been conducted in the region (Guerin et al., 2019; Wolf et al. (2021); Nouibat et al. (2021)). From these observations the 3D group velocities were used to calculate 1D depth inversion for S-wave velocities. The results show that the Moho drops from 12 km below the southwestern center of the Liguro-Provençal basin to 20–25 km below the Ligurian coast in the northeast and to over 30 km depth near the coast of Provence.

Curvature computation has become prominent in geophysics since Roberts (2001). In potential field interpretation, curvature is used to enhance the identification and delineation of subsurface geological structures such as edges of density (or magnetization) contrasts. Its attributes in both seismic and potential methods offer deeper insight into the geometry of the corresponding distributions (e.g., velocity, gravity and magnetic fields). Ebbing et al. (2018) and Li (2015) demonstrated how gravity gradients from recent satellite missions support the interpretation of regional to global gravity fields. However, for our local study, satellite data resolution was insufficient for interpreting laterally rapidly changing structures. Instead, we used the in-house software curvature by Dr. Sabine Schmidt to analyze residual fields processed from Bouguer anomalies (Figure 3C). Free-air and Bouguer anomalies can be used equivalently if the topographic reduction in the Bouguer anomaly accounts for the actual density distribution of the masses between the surface and the reference level. In the case of the free-air anomaly, the effect of these inhomogeneities is still present in the field and must therefore be included in the modeling. Figure 6 shows the “indices of shape curvature” of the residual field (Figure 3B), distinguishing between convex, concave, flat, and ridge manifestations and correlating them with geological structures. Different colors indicate typical shape characteristics: valleys (blue to green), flat areas (yellow), ridges (orange to red), and dome-like convexities (magenta).

The overall shape index picture is complex, but a narrow band with “valley characteristics” (green-blue colors) in the center of the Ligurian Basin can be identified between high elevations off the French Mediterranean coast and Corsica (yellowish and red colors). This contrast may indicate a narrow rift zone, a hypothesis explored further in Section 6.

Terracing and clustering identify patterns, anomalies, and structures in the Earth’s subsurface based on gravity field variations. Both methods were prepared for use in the Ligurian Sea by Strehlau et al. (2022). Terracing partitions a gravity anomaly dataset into distinct amplitude or gradient intervals, or “terraces,” highlighting variations and aiding in identifying subsurface features, geological boundaries, and density anomalies. Cooper (2020) presented the method, which uses the Laplacian function and has recently been replaced by the potential field shape index. The filtering process iterates until the result provides unique homogeneous segments and the “terracing effect” is visible. For the terraced residual field in Figure 3B, the gravity field in Figure 7A shows a simplified but clearer picture. Reddish colors indicate high residual gravity and green-blue colors indicate low values. The result is similar to the shape curvature analysis, revealing a narrow strip of low gravity in the study area center. Clustering applied to terraced gravity field data uses mathematical algorithms, such as k-means or hierarchical clustering, to group similar gravity anomaly measurements into clusters. These clusters represent regions with analogous gravity anomaly patterns, indicating underlying geological structures or subsurface phenomena with common density distributions. This method helps to reveal spatial relationships and coherent structures within the data. The number of clusters, k, is specified by the user. Initially, the cluster centers’ positions are randomly determined and then optimized using the method of least squares. Further methodological details are described in Florio and Lo Re (2018) and Kodinariya and Makwana (2013). In this study, three clusters (k = 3) were used and applied to the residual field computed from the CBA in Figure 3A. The clustering result is shown in Figure 7B. It highlights a narrow zone of low gravity (cluster 1) in the center of the basin, surrounded by a broader cluster (cluster 2). This result illustrates, like the terracing results, the advantage of cluster analysis in representing terraced distributions, effectively reducing the map to its essential content.

Our forward modeling concept involves interactive fitting of potential fields (Free air anomalies) using IGMAS + (Interactive Gravity and Magnetic Application System), a free software tool with nearly 40 years of development (Anikiev et al., 2021; Anikiev et al., 2023b; Götze et al., 2023). IGMAS + uses an analytical solution of the volume integral for gravity and magnetic effects of homogeneous bodies bounded by polyhedrons of triangulated model boundaries (Götze and Lahmeyer, 1988; Schmidt et al., 2011; Anikiev et al., 2023a). The backbone model is constrained by multidisciplinary data such as geological maps, well data, seismic reflection and refraction profiles, structural signatures from seismic receiver functions, and local surveys. The software supports spherical geometries and optimized memory for fast inversion of material parameters and changes in model geometry. Its highly interactive technique makes IGMAS + user-friendly, operating in real time while preserving model topology. Due to its triangular model structure, IGMAS + handles complex structures (multi-Z surfaces) like overhangs of salt domes well. An inversion tool for automated modeling was recently published (Alvers et al., 2023). This software was used to model the gravity field of the lithospheric subsurface to a depth of 300 km in the Liguro-Provençal Basin. Gravity data and constraints were provided by the datasets described earlier.

The 3D model’s basis was the density distribution in the 3D ALPS model (Spooner et al., 2019b; Spooner et al., 2019a). Model structures from 3D ALPS remained unchanged in the eastern and northern parts and were interactively modified in the southern Liguro-Provençal Basin. 3D views of the model are shown in Figures 8A,B. Due to the complex Alpine density distribution between the land surface and the Bouguer anomaly reference level (0 m), we modeled the Free air anomaly (Figure 3C). The disadvantage is that the model’s spatial resolution is sparse, and some details of mountains and model surfaces cannot be represented, especially for Corsica and parts of the French-Italian Alps. Therefore, we truncated the drawing of the Free air anomaly isolines at a certain value. All larger anomaly values appear as pink patches in the maps (e.g., Figures 8A,B, 10A,B). The model is constructed with 33 vertical planes: 31 central planes for modeling and two peripheral planes in the west and east (not shown in Figure 8) to avoid edge effects (e.g., Alvers et al., 2015). The 2.911 model stations were placed on a regular grid with about 7 km spacing (Figure 8A). The reference gravity (Free air anomalies) appears at the top. The lower part of the model (between 50–300 km depth) is represented as a voxel cube, where each voxel is assigned a density value.

In addition to Figure 8A, more model details are shown in Figure 8B. The top surface of the model is hidden, and the stations are turned off. Inside the model, you can see the vertical planes that define the interior and the georeferenced position of the seismic profiles from Figure 4. The perspective view also shows the subsurface structures, indicated by different colors corresponding to model densities in Figure 8C below the model views. Model densities and corresponding colors were used in our model of the Liguro-Provençal Basin which are shown in Figures 8, 9, 11 (after Spooner et al. (2019b); Spooner et al. (2019a)). Note: Density values with six decimal places are calculated by inverting the DENSITY parameter; only the first two digits are significant. The vertical planes (12, 23, 26) will be shown later as examples of the model structure (Figure 9), as well as the seismic profiles P01 and P02 + Makris (Dannowski et al., 2020).

Reliable seismic constraints were available for much of the study area, so the model geometry was not significantly changed, and densities were kept as constant as possible compared to Spooner et al. (2019b)). Free air anomalies were propagated up to 6 km to ensure model stations were always above the model surface, which is up to 3 km for the southern French Alps and 2.7 km for Corsica.

All three model curves in Figure 9 along vertical Sections 12, 23, and 26 in Figure 8B (green/white stippled lines) and 10 (red lines) fit the reference anomaly. This is evident in the residual map below (Figure 10C). Since the 3D ALPS model of Spooner et al. (2019b) was compiled for the Central Alps and the Ligurian Sea margins, we slightly modified the geometry of the geological bodies. Our workflow was to fit the model Moho such that it was consistent with the Moho compilations from seismic campaigns (profiles in Figure 4) and in areas not covered by seismic lines with the Moho of Grad et al. (2009). Then, the geometries of the geological bodies were adjusted, and their densities were changed (inverted) if necessary. Differences between reference and model gravity are mostly less or equal than 5 × 1 0 − 5 m / s 2 (5 mGal).

All densities and corresponding colors in Figure 9 match those in Figure 8C. Section 12, located at the far west of the area, intersects a gravity high offshore of France, followed by another gravity high further south. The first high is attributed to near-surface structures (consolidated sediments), while the second is due to an elevated Moho position. Subsurface structures along vertical Section 23 differ significantly further east. Here, a pronounced, localized gravity high on land (green body, “u crust Apennine,” density 2,720 kg/m³) is bordered by a relative gravity minimum, caused by increasing water depth, the tapering of the 2,720 kg/m³ body, and the rise in 2,400 kg/m³ “consolidated sediments.” A similar pattern appears in Section 26’s gravity profile, east of the previous section. Initially dominated by the high-density 2,720 kg/m³ body, the profile shows a uniform trend until gravity increases again in the south, due to the higher densities of the Corsica-Sardinia block. Refer to Figure 11 for comparison.

The crustal domain in the model is supported by extensive seismic survey data (Figure 4). This is not the case for the upper mantle, whose architecture has recently become better known through the AlpArray Consortium tomography and the 4DMB SPP, but still lacks clarity on whether interpreted plates are consistent with observed surface deformation and topography. For mantle modeling, we relied on Kumar et al. (2022), combining shear wave tomography models into a statistical ensemble to derive three scenarios for plate lithospheric thickness and geometry. These scenarios were used for geodynamic simulations to calculate topography, surface velocities, and mantle flux. The scenario representing a detached plate in the Alps and a connected plate in the northern Apennines, consistent with observed seismicity at intermediate depths, was adopted. It captures the study area’s topography and vertical surface velocities, serving as a first approximation to a constrained upper mantle.

The gravity field caused by the density inhomogeneities in the upper mantle between 50 and 300 km depth is shown below in Figure 15. This figure shows that neglecting the gravity anomalies in the upper mantle significantly effects on the shape and magnitude of the gravity in the study area: it varies from − 30 × 1 0 − 5 m / s 2 (−30 mGal) in the southwest to 70 × 1 0 − 5 m / s 2 (70 mGal, maximum) in the Apennines in the northeast.

In addition to the north-south profiles through the 3D model shown in Figures 9, 11 presents 3D model slices along arbitrary orientations—specifically along profiles P01 and P02-Makris (refer to Figure 4). This allows for a direct comparison between the 2D modeling by Dannowski et al. (2019) and our recent 3D modeling (Figures 11A,B), confirming overall agreement between the two. Notably, the free-air gravity rise in the Corsica region (Figure 11A, southeast) is primarily driven by topographic elevation. The significant gravity low around profile km 290 in the 3D model is attributed to thicker sediment layers, consolidated sediments (depicted in yellow and blue), and the submersion of mantle structures (gray and orange). Both models reveal similar general structures and density distributions; however, the 3D model provides a better fit to the gravity field in specific areas, such as west of Corsica.

Figure 11B highlights a dipping Moho interface toward the Italian mainland. Both models accurately represent the crust-mantle transition zone in the Liguro-Provencal Basin. The low-velocity zone (LVZ) in the upper serpentinized mantle may also be interpreted as a high-velocity zone (HVZ), potentially indicating rift-related underplating in the crust. The modeling confirms higher average densities and elevated gravity values in the central Ligurian Basin. Velocity and density modeling results are consistent with the seismic findings of Dannowski et al. (2020), which report a 6–8 km thick sedimentary cover and a seismic Moho at depths of 11–15 km below the sea surface.

In summary, the new constrained 3D modeling shows a satisfactory fit between the Free air anomaly and the modeled Free air anomaly. Deviations in the residual gravity map are within ± 5 × 1 0 − 5 m / s 2 ( ± 5 mGal) in Figure 10. The error we report for the modeled field is also ± 5 mGal. Larger deviations occur only in the southwestern Alps and Corsica, primarily due to the need for more vertical planes to approximate the topography/Free air anomaly in mountainous areas. However, adding more planes for extensive model areas, such as the entire offshore area, would unnecessarily complicate the model.

Gravitational Potential Energy (GPE) is a key concept in geophysics, representing the energy stored in an object due to its position in the gravity field. It is directly proportional to the object’s mass and its height relative to a reference point. The modeling software calculates the GPE at the gravity stations’ positions, resulting in a 1D calculation. The 2D distribution shown in the maps in Figure 12 is derived through interpolation.

Background: 1) Dependence on Mass: GPE is directly proportional to an object’s mass. Heavier objects have more GPE at the same height. 2) Dependence on Altitude: GPE is directly proportional to the height above a reference point. Higher elevations result in greater GPE. 3) Gravitational Field Strength: On Earth, g is approximately 9.81 m/s² but varies with location and altitude. 4) Tectonic Reference Level (TRL): The choice of reference level is crucial for GPE calculations and vertical stress. Following Coblentz et al. (1994), we assume an average value of 2.373 × 1,014 N/m, equivalent to the potential energy of continental plates and basins with topography below sea level.

GPE helps explain various geophysical phenomena. Applying current processing methods to our gravity data set provides valuable insights, as demonstrated in previous studies (e.g., Ghosh et al., 2009; Flesch and Kreemer, 2010; Schmalholz et al., 2014; Neres et al., 2018).

The GPE per unit area (A) is defined as:

G P E = m g h / A , given m = ρ h A we get:

G P E = ρ g h 2

where m is mass, g is gravity, h is height, and ρ is rock density. For inhomogeneous density columns, we use piecewise constant densities:

G P E = g ∑ ρ i h i

The unit of GPE is k g / s 2 ( J o u l e / m 2 or N / m ). Comparing GPE with stress, we see:

G P E = ∑ σ i h i

where σ i is the stress contribution of layer i with constant density ρ i and thickness h i . In a constant density environment, stress increases linearly with depth, while GPE increases quadratically.

The Ligurian Sea’s stress regime results from the convergence of the African and European plates. This collision shapes the region’s tectonics. Additional factors influencing the Ligurian margin include the rollback of the Ionian-Adriatic subduction, gravitational collapse of the Apenninic lithosphere, and lateral extrusion of the southwestern Alps (Morelli et al., 2022). Due to this complexity Eva and Solarino (1998) suggested that the Ligurian Basin should be analyzed separately from the adjacent Alpine chain. Compression within the Ligurian Sea may reactivate fault systems in the Nice Arc. The GPE distribution in the Liguro-Provençal Basin reveals a complex interplay of tectonic forces. Positive GPE values in the central region and around Corsica indicate thicker or denser crust and this stored energy translates into compressional forces as the crust attempts to reach gravitational equilibrium. Lower GPE near the French coast suggests thinner or less dense crust. These areas are often associated with compressional forces rather than extensional forces, as the crust is under less gravitational potential energy stress. This complex stress environment results from the Ligurian Basin’s opening or extension, influenced by broader plate interactions. The coexistence of these varying GPE values highlights the dynamic interplay of tectonic forces in the region, which is corroborated by independent studies that describe similar heterogeneous stress distributions and their effects on geological structures (e.g., Baroux et al., 2001).

Euler Deconvolution (ED) is a mathematical method used in geophysics to estimate the depths and locations of underground sources in potential fields (Reid et al., 1990; Paŝteka, 2006; Saleh and Pašteka, 2012). ED analyzes gradients or derivatives of observed potential field data to infer subsurface sources, which are modeled as simple geometric bodies (spheres, cylinders). ED provides estimates for:

• Depth (Z) of the source below the Earth’s surface.

ED assumes simplified source shapes and may not always represent complex geological structures accurately. It is often used with other geophysical and geological data to refine interpretations and improve subsurface exploration accuracy. In our analysis, ED is applied to both Bouguer and Free air gravity using the REGDER software (Pašteka, pers. comm).

The main parameters for the Euler Deconvolution (ED) are specified above the Figure A1 in the supplementary materials: for the Bouguer anomaly analysis, a window size (WS) of 7 units is used, and for the Free air anomaly (FA), the WS is 5 units. The structure index is set to 2 for both fields, indicating a three-dimensional borehole. In the Liguro-Provençal Basin, the analyses differ for the two fields. For the Bouguer Anomaly (BA), only a few source points are identified in the Ligurian Sea, with most sources located around Corsica and the southwestern Alpine arc. The concentration of source points in the southern area aligns with the anomaly in the shape curvature (Figure 6). The FA analysis also shows a concentration of Euler source points, though strong gradients in the Free air anomaly elsewhere limit its utility. Regarding the question of whether the gravity field indicates a possible rift structure, the ED results alone do not provide a conclusive answer. Detailed ED results are provided in the Appendix.

Rollet et al. (2002) discussed back-arc extension, tectonic inheritance, and volcanic activity in the Ligurian Sea, identifying three crustal domains: 1) continental thinned margins, 2) transition areas to the basin, and 3) a narrower atypical oceanic area. Margin structures feature tilted blocks forming syn-rift sedimentation and segmentation. Using Rollet et al.’s (2002) nomenclature, the terraced Bouguer anomaly (Figure 13B) shows: “orange” colors represent “margins,” dark pink/reddish colors indicate the “transitional domain,” and lighter pink areas in the central part represent “atypical oceanic crust.” Structurally, this area is at a relatively high bathymetric position in the 3D model (Figure 9). This fit is expected as Rollet et al. (2002) also based their interpretation on gravimetric and magnetic data.

With the results of Canva et al. (2021) we identify a correlation between areas of possibly exhumed mantle and the central positive anomaly in the Free air gravity anomaly. According to their interpretation, the bright pink stripe in the Bouguer anomaly (Figure 13B) could also represent areas of exhumed mantle. This correlation is supported by travel time tomography results from the SEFASILS cruise (Canva et al., 2021) and the southeastern part of the P02-Makris profile (Dannowski et al., 2020). In the neighboring Tyrrhenian Sea, seismic studies reveal exhumed mantle that is linked to the opening of the basin (Prada et al., 2016). There, the authors could constrain in space and time that processes of mantle exhumation, and later magmatic intrusions of MOR-type and intraplate basaltic rocks were involved in the opening of the basin. Consistent with Rollet et al. (2002), we interpret the darker pink area in the Bouguer anomaly as hyper-extended crust formed during the rifting phase.

In areas with continental thinned margins, the Moho depth is around 20–25 km due to significant crustal thinning from extensional tectonics. Transitional zones between continental margins and oceanic crust typically exhibit Moho depths of 15–20 km. In the narrower atypical oceanic areas, the Moho is shallower, around 10–15 km, indicating more significant crustal thinning and extension compared to continental margins. Overall, the Moho depth in the Ligurian basin ranges from 10 to 25 km, varying with the specific tectonic and geological context of each sub-region.

In the late 1970s, Giese et al. (1978) provided early insights into the complex offshore conditions of the Ligurian Sea, Western Alps, and Apennines. In “Alps, Apennines, Hellenides” (Closs et al., 1978), they noted the static and dynamic interdependencies of the Adriatic microplate with the Western Alps and the Apennine Arc, assuming an oceanic crust with a mantle depth of 15 km in the Ligurian Sea. Our 3D modeling (Section 5.3.1) shows a Moho depth of about 12–16 km in the central Ligurian Sea. These results, combined with structural comparisons (Figure 14), indicate that the Ligurian Sea’s opening and boundary structures are controlled by the interplay between inherited crustal heterogeneity and external tectonic forces, such as subducting lithosphere rollback.

According to Rollet et al. (2002), various crustal domains can be identified, which are also reflected in our interpretation of gravity field data. The first domain, demarcated by a reddish-pink line and labeled “TD,” corresponds to the transitional domain. The white line, labeled “AOD,” delineates the extent of the (inner) “atypical oceanic domain,” characterized by an elevated Moho position and a central gravity high in both the free-air and Bouguer anomalies (Figure 3). The bright yellow line in the center, labeled “ODG,” represents the “oceanic domain” as defined by Gueguen et al. (1998), marked by the central maximum in the gravity field. Light and dark blue lines, based on the interpretation of Asch (2005) and previously illustrated in Figure 13A, differentiate between oceanic crust in the center and a “rifted, thinned continental crust.” The latter is distinguished by lower free-air anomaly gravity values and significantly deeper Moho levels compared to the central high. Additionally, a red dashed line off the west coast of Corsica in the east indicates a proposed fracture zone, characterized by a steep gradient in the gravity field. No significant differences were found between our results and those from the referenced studies.

Dessa et al. (2020) characterized the northern Ligurian Sea as a geologically complex offshore region, influenced by multiple phases of tectonic deformation. They emphasized the role of salt structures, particularly along a narrow zone approximately between coordinates 43°20′N/8°00′E and 43°50′N/8°30′E, as well as deeper crustal tectonics, in controlling sedimentation patterns. Due to the regional scale of our 3D model, these localized features are too small to be resolved at high spatial resolution. Nonetheless, a local minimum in the free-air gravity anomaly (Figure 3C) and the interpretation of gravitational potential energy (GPE) in conjunction with the vertical stress map (Figure 12) may indicate morphological effects and accumulation of salt deposit accumulation within the crust.

Morelli et al. (2022) provided a detailed reconstruction of the morpho-structural setting and tectonic evolution of the Alpine and Apennine margins and the Ligurian Sea. Their analysis of bathymetric and geophysical data reveals a close relationship between recent active tectonics and the morpho-dynamic evolution of the Ligurian basin margin. This relationship is particularly evident in the Alpine margin, where uplift and inherited structures under compressive regimes are associated with widespread seismicity, showing spatial correlations with the GPE and stress maps in Figure 12. Figure 15 presents a map of the voxel cube gravity effect superimposed on earthquake distribution, indicating an accumulation of earthquakes in areas with higher gravity field gradients. The voxel cube contains densities converted from velocities at depths between 50 and 300 km, highlighting the notable coincidence of hypocenters and gravity field shapes.

Our investigations reveal highly segmented gravity fields in the Liguro-Provençal Basin, as evidenced by the residual Bouguer gravity (Figure 3B), Free air anomaly (Figure 3C), shape curvature (Figure 6), and clustering of the terraced gravity fields (Figures 7, 13B). This level of resolution was previously unattainable. However, it remains uncertain whether these anomalies are definitively part of a rift structure. Comparisons with other regions, such as the South China Sea (e.g., Hagen, 2021), show much clearer gravity signals associated with rift structures. Seismic data from the LOBSTER experiment and structural evidence from the sefasils experiment suggest the deep reflection horizon in the south corresponds to the Moho, while shallower reflections in the north are attributed to the base of the sediment (Canva et al., 2021). The extensive salt complex observed may account for the narrow minimum in our residual gravity maps in the central Liguro-Provençal Basin. However, these local structures were not the primary focus of our 3D modeling due to the required resolution.

From the perspective of gravity field analysis (Bouguer, Free air anomaly, and residual fields), terracing and clustering, curvature, GPS, and 3D modeling indicate that rifting has likely influenced the basin’s geological evolution. To summarize these findings, we provide a hierarchical processing approach that could be beneficial for similar studies in other regions:

• AAGRG and smoothed Free air anomaly (chapter 3.3) ⇒ strong indication