In the present study, we perform a gravity modelling at crustal scale along the trace of the CROP-04 (on-shore) and M-6B (off-shore) deep seismic reflection profiles, crossing the Southern Apennines and the Southern Tyrrhenian Sea (Italy). Along the 321 km-long modelled profile, we investigate the crustal-scale sources for the observed gravity anomalies through a simplified model of the crust and upper mantle across both onshore and offshore areas.

After a compelling review of the published Moho geometries in the area, that were retrieved from either active or passive seismic methods, we test them in the observed gravity field through forward modelling of the Bouguer gravity anomalies. The comparison between the different Moho interpretations highlights the major contributors to the observed Bouguer gravity at the crustal scale, defining a set of starting values of these parameters for our final model.

The proposed model locates the westward flexure of the Adriatic Moho, mimicking the subduction of the Adriatic lithosphere beneath the Peri-Tyrrhenian block and locates the step between the western (Tyrrhenian) and the eastern (Adriatic) Moho beneath the Apennines range providing a valuable geometrical and compositional model at the crustal scale. The model depicts a typical oceanic-to-continental crust transition in the Tyrrhenian domain and represents a solid starting base for further detailed modelling across the area.

In areas where the geodynamic and structural evolution results in a complex crustal setting, gravity forward modelling can provide significant contributions in constraining the deep tectonic structures at regional scale and the Moho topography (e.g. Reyet alii, 1990; Lefort & Agarwal, 2002; Riveroet alii, 2002, Düzgitet alii, 2006; Tassiset alii, 2013; Mancinelliet alii, 2015).

The Southern Apennines of Italy are the NW-SE trending, Neogene-Quaternary foreland fold-and-thrust belt, formed in the complex framework of the Africa-Eurasia convergence overlying the SW-dipping Apulian-Adriatic continental lithosphere (e.g. Scrocca, 2010 and references therein). The NE-verging fold-and-thrust belt and the associated foreland basin migrated eastward through time, incorporating a complex assemblage of tectonic units, involving Mesozoic-Cenozoic successions, mostly deposited on the Adriatic continental passive margin, in carbonate platform or pelagic basin environment (D’Argenioet alii, 1975). The emplacement of the compressional belt was synchronous with co-axial extension at its rear, opening the Tyrrhenian basin and affecting the internal (western) zone of the mountain belt. The intense seismic activity of the axial zone of the Southern Apennines, characterized by earthquakes with maximum magnitude up to 7.0, reflects the still on-going, WSW-ENE trending extension (Amato & Montone, 1997; Montone & Mariucci, 2016).

According to most Authors (following the pioneer works of Malinverno & Ryan, 1986; Patacca & Scandone, 1989), the development and space-time evolution of the compressional and extensional structures of the Apennines/Tyrrhenian system are driven by the roll-back of the passively subducting paleo-Ionian lithosphere, also driving the opening of a back-arc Tyrrhenian basin, formed at its rear (Amatoet alii, 1993; Piromallo & Morelli, 1997; Wortel & Spakman, 2000; Cimini & De Gori, 2001; Lucente & Speranza, 2001; Di Stefanoet alii, 2009). This process is also responsible for the present-day crustal structure of the region, where the Adriatic crust (25 to 50 km thick) progressively deepens westward from the foreland area to the mountain belt, over a lithospheric thickness ranging between 70 and 100 km (Calcagnile & Panza, 1980; Nicolich, 1989; Scarasciaet alii, 1994; Di Stefanoet alii, 2009). On the western side of the peninsula, a thinned peri-Tyrrhenian continental crust borders the oceanic crust of the Southern Tyrrhenian Sea (Nicolosiet alii, 2006), in a region also characterized by high heat flow (Nicolich, 1989; Scarasciaet alii, 1994; Cataldiet alii, 1995; Della Vedovaet alii, 2001). In contrast, low heat flow values are measured above the Adriatic lithosphere, which was already structured during the Mesozoic rifting stages (Scroccaet alii, 2005).

Since the 1960’s, the complex crustal structure of the Southern Apennines-Southern Tyrrhenian system was explored using different geophysical surveys and techniques, comprising the seismic reflection profiles collected during the Deep Seismic Sounding (DSS) experiments in the 1960s–1990s (Scarasciaet alii, 1994; Nicolich, 2001; Cassiniset alii, 2003) and the deep seismic reflection profiles of the CROP Project in the 1980s–1990s (Scroccaet alii, 2003; Finetti, 2005; Cassiniset alii, 2007). In more recent periods, passive seismic methods, based on the natural earthquakes registered by both fixed and temporary seismic networks, offered further constraints to the crustal structure of the region (Steckleret alii, 2008; Di Stefanoet alii, 2009; Piana Agostinetti & Amato, 2009). Furthermore, significant contributes to unravelling the deep structure of the region were offered by other geophysical techniques, such as the inversion or forward modelling of gravity and magnetic data (e.g. Fediet alii, 2005; Tibertiet alii, 2005; Anelliet alii, 2007; Speranza & Chiappini, 2007) and the magnetotelluric profiles (Patellaet alii, 2005).

In this paper, we applied forward model techniques to a dense grid of gravity data, derived from the 1:500,000 Gravity Map of Italy (Bigiet alii, 1991), aimed to characterize the lateral extent and densities of the crustal-scale bodies contributing to the observed Bouguer anomaly and delineate depth and geometry of the Moho discontinuity. The final result is a simplified model of the crust and upper mantle structures, down to a depth of 70 km, along a regional transect across the Southern Apennines and the Southern Tyrrhenian Sea. The investigated transect is arranged along two WSW-ENE trending deep seismic reflection profiles, acquired in the framework of the CROP project (Scroccaet alii, 2003): the M-6B profile, crossing the Southern Tyrrhenian offshore, and the CROP-04 profile, crossing the onshore areas of the Southern Apennines fold-and-thrust belt, of the Bradanic foredeep and of the Apulian foreland (Fig.1).

The Bouguer anomaly map of the study area, presented in Fig. 2a, derives from the dataset of the 1:500,000 Gravity map of Italy that integrated data originally acquired by AGIP and CNR, consisting of ~270,000 station measurements. From this dataset, we derived a regular grid of 85952 points (1 x 1 km gridded), using a density of 2670 kg m-3 for the calculation of the Bouguer reduction.

Along the modelled profile, the Bouguer anomaly ranges between 20 and 200 mGal (Fig. 2b). Maxima are located at the western end of the transect, in the central part of the Southern Tyrrhenian Sea, at ~ 100 km from the coastline. From this point, an eastward decreasing regional trend is observed, reaching a minimum of 20 mGal in correspondence with the Southern Apennines mountain belt. Note that significantly lesser values of gravity are achieved along the strike of the mountain range, both towards NW and SE (about -20 and -40 mGal, respectively). The gravity values gently increase in the eastern part of the profile, towards the Adriatic Foreland.

The topography and bathymetry data (Fig. 1) were retrieved from the Global Relief Model (GEBCO Compilation Group, 2020), obtained by merging the Digital Elevation Model (DEM) onshore and the Digital Bathymetry Model (DBM) offshore. From west to east, the topographic profile (Fig. 2c) can be divided into 4 segments: (i) the flat and deep floor of the Southern Tyrrhenian Sea (at about 3500 m b.s.l.), corresponding to the Marsili basin, whose oceanic nature is proven by both air- and ship-borne magnetics (Nicolosiet alii, 2006); (ii) the continental slope, connecting the basin to the emerged area through a steep lower step (up to about 1300 m b.s.l.), followed by a relatively gentler upper slope; (iii) the mountain range of the Southern Apennines, about 100 km large and reaching a maximum elevation of about 1300 m a.s.l.; (iiii) the eastward sloping Adriatic foothills, merging into the almost flat foreland area.

We performed the gravity forward modelling along a WSW-ENE trending regional transect, extended down to a depth of 70 km. The total length of the modelled profile is 323 km and encompasses the 131 km-long M-6B and 154 km-long CROP-04 traces (Fig.2a). Gravity models were produced using the GM-SYS software (Geosoft, 2013).

Our modelling procedure consisted of two steps. The first step was mainly aimed to define the Moho geometry. To achieve this goal, we tested seven different models of the Moho depth along the studied profiles that were retrieved from literature, and considered representative of the possible different interpretations. This modelling was performed using a simplified structure, consisting only of two layers, i.e. crust and upper mantle, with density values of 2800 kg m-3 and 3200 kg m-3, respectively, separated by the Moho discontinuity.

Starting from the suggestions about the Moho geometry, derived from the first step, the second step was aimed to the construction of a synthetic, 4-layer model, describing the structure and the lateral variations of density occurring of the upper crust, lower crust and upper mantle, locally covered by lighter, recent sediments. This model is based on the compilation of the density values available in the literature, along with some observations derived from the topographic/bathymetric profile.

In order to test the Moho geometry along the studied section, we went out a careful review of previously published works, concerning the investigation of the crustal structure of the Southern Apennines through active or passive seismic techniques. Among the many different proposed interpretations, we selected seven models, representative of the variability of the proposed solutions, and deriving from active wide-angle (DSS) seismic profiles (Scarasciaet alii, 1994; Nicolich, 2001; Cassiniset alii, 2003), deep crustal seismic reflection (CROP) profiles (Finetti, 2005), passive teleseismic receiver functions (Steckleret alii, 2008; Piana Agostinetti & Amato, 2009) or deriving from a combination of active (refraction and reflection seismic) and passive (receiver function) data (Di Stefanoet alii, 2011).

A quick comparative view of the Moho depth and geometry proposed by the selected papers is shown in Fig. 3. Almost all the models, with the exception of Steckleret alii (2008), who did not analyse the Tyrrhenian side of the section, share some common traits, highlighting the presence of two distinct Moho surfaces, characterising the western and the eastern part of the transect. Namely these are the shallower, gently east-dipping Tyrrhenian Moho beneath the Tyrrhenian Sea and the western part of the Italian peninsula; and a deeper and steeper, west-dipping Adriatic Moho in the eastern side of the profile. Moreover, the studies based on active seismic data, suggest that in the central part of the transect, the Tyrrhenian crust is superposed to the Adriatic crust (Fig. 3). On the other side, the presence of Moho doubling in the central part of the section is only hypothesized in the studies by Piana Agostinetti & Amato (2009) and Di Stefanoet alii (2011), that in any case clearly image a sharp step between the shallower Tyrrhenian and the deeper Adriatic Mohos.

However, the seven models show other significant similarities and differences when observed in detail (Fig. 3). At the western end of the section, beneath the Tyrrhenian basin, the depth of the Tyrrhenian Moho ranges from 8-10 km (Nicolich, 2001, Finetti, 2005) to more than 20 km (Cassiniset alii, 2003). Moving NE-ward, the Tyrrhenian Moho gently dips beneath the Apennines, to converge to a depth of 25-30 km beneath the mountain belt. At the eastern end of the profile, the Adriatic Moho shows a reasonably restricted depth range (28-33 km) and dips SW-ward, reaching a maximum depth of ~50 km beneath the main mountain ridge (Scarasciaet alii, 1994; Cassiniset alii, 2003; Piana Agostinetti & Amato, 2009). The central part of the profile shows the largest differences among the proposed crustal interpretations: the studies based on active seismic data (Scarasciaet alii, 1994; Nicolich, 2001; Cassiniset alii, 2003; Finetti, 2005) image a region of Moho doubling, extending from the Tyrrhenian coastline to the axial region of the Apennines, where the Tyrrhenian Moho is superposed to the Adriatic Moho. The presence of such Moho doubling was also explicitly hypothesised in the study by Piana Agostinetti & Amato (2009). Finally, a more complex geometry is imaged by Finetti (2005), including an intermediate Moho step, just beneath the Tyrrhenian coastline.

In order to test their degree of correlation with the observed gravity signals, we calculated the Bouguer gravity anomalies for the seven crustal models, described above (Fig. 4). For this test, we used an extremely simplified two-layer structure: i.e. crust and lithospheric mantle, separated by the Moho discontinuity, with density values of 2800 kg m-3, 3200 kg m-3, respectively. For the models made after Steckleret alii (2008) and Piana Agostinetti & Amato (2009), which are limited to the onshore section, the Moho surface was extended westward with an almost flat geometry. In the Southern Tyrrhenian Sea, the lithosphere is thin (about 30 km, Calcagnile & Panza, 1980; Panzaet alii, 2003), as also confirmed by the high heat flow (Della Vedovaet alii, 2001), due to the presence of a shallower asthenosphere. Therefore, in the western part of all the models we added a third layer depicting the eastward-deepening asthenosphere, with a minimum depth of 30 km at the western end of the section, and we assigned to this layer a density of 3130 kg m-3, slightly lighter than that of the overlying mantle.

The results of the modelling are illustrated synthetically in Fig. 4a-4g, labelled with reference to the seven considered models. In all the models, the long-wavelength gravity anomalies are mimicking the Moho geometries, respecting a general pattern, where regional gravity-highs correlates with Moho shallowing while regional gravity-lows are related to Moho deepening. At the same time, all the models show significant discrepancies between the calculated and observed gravity anomalies, due to the lateral and vertical density variations within the crust, which are not considered in the oversimplified, 3-layer modeled structure.

The Moho geometries proposed by models in Fig. 4a and 4b show the better overall fitting, whilst the model in Fig. 4c results in a marked westward shifting of the calculated anomaly in respect to the observed anomaly, where the minima are located beneath the Peri-Tyrrhenian on-shore area. The shallow Tyrrhenian Moho in the south-western end of the modeled profile of Fig. 4d generates a pronounced misfit, where the calculated anomaly is significantly lower than the observed anomaly while in the central part it is higher than the observed regional minima. The model in Fig.4e shows an overall good fit along the transect. An exception is found in the NE part of the model, where the calculated anomaly is significantly lower than the observed. The model in Fig. 4f provides a general good fit in the eastern part, but shows a misfit in the central part of the transect, where regional minima is observed. Similarly, in the western part of the model a misfit is observed, but this is not related to the tested reference model because the passive seismic data in that study does not cover the offshore region. Also the model in Fig. 4g, generally fits the observed gravity minima in the onshore part of the profile, where the reference model provides reference Moho depths.

In general, the thickness of the Tyrrhenian and Adriatic Moho, the steepness of the westward-dipping Adriatic Moho and the location and extent of the Moho doubling in the central part of the profile are the most important parameters affecting the gravity signature. The combination of these parameters governs the calculated Bouguer anomalies at the modeled scale and has strong influence on the location of the minima and the regional trend along the entire profile.

Keeping in count the results of these preliminary gravity models, along with the topographic/bathymetric profile of Fig. 2c, we adopted the Moho geometry, shown by the dashed line in Fig. 3. This model is generally consistent with the Moho depth extracted from the past studies. As in most previous interpretations, we distinguish two different Moho discontinuities: a newly-formed Tyrrhenian Moho and an older Adriatic Moho. In the SW end of the profile, the relatively shallow Tyrrhenian Moho flattening at about 12 km depth, reflects the presence of a deep marine basin (about 3500 m b.s.l.). In the NE end (i.e. the Adriatic coastline) the Adriatic Moho is at a depth of ~30 km, flattening at 28 km beneath the foreland area. The central part of the profile, beneath the mountain ridge, is characterized by the superposition of a gently east-dipping Tyrrhenian Moho at about 27 km and a deeper and steeper, west-dipping Adriatic Moho, reaching a maximum depth of ~50 km. The adopted Moho geometry was tested with the simplified 3-layer model (Fig. 4h) and then used as a base for our final forward model of the Bouguer anomalies along the transect.

After defining a reliable Moho geometry and structure, in order to fit and properly model the vertical and lateral density contrasts between the main geological bodies, it is necessary to assign them appropriate density values, considering the stratigraphic and structural setting of the study area.

Table 1 offers a comprehensive view of the main lithological units, occurring along the considered profile, as retrieved from previously published works (Tibertiet alii, 2005; Improtaet alii, 2003; Mostardini & Merlini, 1986; Anelliet alii, 2007; Tassiset alii, 2013; Biella, et alii, 2007; Menardi Noguera & Rea, 2000; Mancinelliet alii, 2019). With the aim of providing a preliminary crustal model along the considered transect, in this study we decided to test a simplified structure, consisting of five layers, top to bottom: Pliocene-Quaternary sediments; Upper Crust; Lower Crust, Lithosperic Mantle and Asthenosphere, disregarding the internal complexity of the continental upper crust.

We are conscious that this is a gross simplification. The relatively shallow structures, formed by the tectonostratigraphic units of the Southern Apennines fold-and-thrust belt, are object of a wide literature, proposing alternative models (e.g. thin-vs. thick-skinned) for the tectonic style of the fold-and-thrust belt (e.g. Menardi & Nouguera 2000; Vezzaniet alii, 2010; Butleret alii, 2004; Shineret alii, 2004; Finetti, 2005; Brozzetti, 2011; Cippitelli, 2007; Scroccaet alii, 2007; Anelliet alii, 2007; Speranza & Chiappini, 2007; Scrocca, 2010; Mazzoli, 2013). Other studies envisage the widespread presence of Early Tertiary igneous rocks within the upper crust (e.g. Improtaet alii, 2014). The complex internal structure and the overall thickness of the upper crust may significantly affect the medium-scale gravity signals. However, here we are mainly interested in modelling the long-wavelength gravity signals, reflecting the crustal structure of the region.

Pliocene-Quaternary sediments are present on top of the modelled section, but only locally they reach significant thickness, affecting the gravity field at regional scale. This is mostly the case of the Bradanic trough, representing the youngest foredeep basin of the Southern Apennines, with a maximum thickness of 1500 m (Patacca & Scandone, 2007), and of the Tyrrhenian basin seabed, where a thickness of sediments is ~600 m. The latter value is constrained by the ODP well #650, which recovered ~600 m of Pleistocene–uppermost Pliocene clastics, mudstones, and nannofossil oozes (Kastens & Mascle, 1990) and was plotted by Finetti, (2005) along the M-6B deep seismic profile. We modelled these shallow sediment sources, assigning them a laterally homogeneous density value of 2300 kg m-3 (Table 1).

For the other major layers of our model (i.e. upper crust, lower crust and upper mantle), we assigned different density values for the different structural domains of the considered transect.

In the western side of the transect, under the Tyrrhenian Sea, we hypothesized the presence of oceanic crust, whose stratigraphy consists of two layers: a ~2 km thick upper layer, with a density of 2800 kg m-3 (upper oceanic crust) and a ~6 km thick lower layer, with a density of 2850 kg m-3 (lower oceanic crust). This vertical superposition mimics the typical stratigraphy of the oceanic crust, with a sedimentary layer superposed to an upper crust of vesicular basalts, whose uppermost section (32 m) was penetrated by the ODP well #650 (Kastens & Mascle, 1990), overlying a lower crust heavily intruded by gabbros (Pontevivo & Panza, 2006; Savelli & Ligi, 2017).

In the onshore area, within the continental upper crust of the Southern Apennines, we have considered two different units, separated by a tectonic boundary, corresponding to the frontal thrust of the Southern Apennines: an eastern, Adriatic upper crust, consisting entirely of platform carbonates with average density of 2650 kg m-3; and a western, Apenninic upper crust, consisting of a mixed lithological assemblage, containing also lighter rocks of pelagic basin environment, with an average density taken as 2600 kg m-3. For the continental lower crust, below both the Apenninic and the Adriatic upper crust, we considered an homogeneous density of 2850 kg m-3.

Further lateral variations were introduced in the peri-Tyrrhenian region, corresponding to the continental platform and slope, representing the transition between the higher-density oceanic Tyrrhenian crust and the lower-density continental Apenninic/Adriatic crust. In this transitional zone, between km 50 and 100 along the modelled profile, we assigned a value of 2700 kg m-3 for the upper crust and a value of 2850 kg m-3 for the lower crust.

Beneath the Tyrrhenian and the Adriatic Moho, the Lithospheric Mantle was also differentiated, hypothesizing a different age and thermal state. An higher density value (3300 kg m-3) was assigned to the older and colder Adriatic Lithospheric Mantle, in respect to the presumably hotter and lighter, newly-formed Tyrrhenian Lithosperic Mantle (3200 kg m-3). A similar scenario was suggested along the CROP-03 profile where different physical properties have been envisaged for the two different Tyrrhenian and Adriatic domains, as suggested by both geophysical and geological data (Barchiet alii, 1998; Pauselliet alii, 2006).

Considering the previously defined Moho geometry and crustal stratigraphy, with the associated density values, we used gravity forward modelling techniques for refining the geometry of the different crustal blocks, composing our final “best-fit” crustal model. The curve in the upper panel of the Fig. 5 shows and compares the calculated and the observed gravity patterns along the transect.

As already observed in the previous paragraph, the long-wavelength (about 200 km) undulations of the observed Bouguer anomaly, mostly derive from the Moho geometry, but are also affected by the lateral variations of thickness and density of the three major layers, i.e. upper crust, lower crust and upper mantle. At the north-eastern end of the profile, the positive trend of the observed anomaly is related to the shallowing of the Adriatic Moho to depths of ~30 km near the Adriatic coastline. In the central part of the transect, the gravity minimum corresponds to the crustal thickening (doubling) underneath the Southern Apennines, where the Tyrrhenian Moho (~27 km depth) overlies the deeper Adriatic Moho (~50 km depth). In the south-western part of the section, the gravity maximum is related to the thinning of the Tyrrhenian crust, reaching here its minimum thickness of ~12 km, in correspondence of the abyssal plain of the Southern Tyrrhenian Sea.

Superposed to this regional trend, local undulations of the observed Bouguer anomaly, with wavelengths between 10 and 20 km, represent the effect of shallower sources. These are related to the relatively low-density Pliocene-Quaternary deposits, infilling the basins in the Tyrrhenian side and the Bradanic trough (yellow blocks in Fig. 5).

Keeping in mind the effects of the deep (i.e. crustal structure) and shallow (recent basins) sources described above, the calculated gravity curve mostly fits the observed values. The residual short-wavelength discrepancies between the calculated and the observed curves, are here interpreted as possibly related to the heterogeneous upper crustal structure of the continental crust of the Southern Apennines thrust belt, which have not been addressed in this work.

Starting from original measurement stations, we created the Bouguer anomaly map of a wide region of Southern Italy, encompassing the Southern Apennines and the adjacent Tyrrhenian offshore, with a reduction density of 2670 kg m-3. To investigate the crustal geometries of this region, we performed a forward gravity modelling along the M-6B and CROP-04 deep seismic profiles, from the Central Tyrrhenian Sea to the Adriatic shoreline.

The proposed final model (Fig. 5), reasonably fits the observed gravity field, and is coherent with the first-order geological and geophysical constraints:

  • - the geometry of the Moho fits and averages the results of previous works, performed in the last decades across the study area;

  • - the density values adopted for the crustal scale geological bodies are also fully consistent with previous studies;

  • - considering the geometries and densities of the modeled crustal bodies, the proposed model represents a typical oceanic-to-continental crust transition and is consistent with the topographic/bathymetric data, depicting the boundaries between the abyssal plain, the continental slope and platform and the emerged area.

The Moho doubling imaged the central part of the section does not necessarily imply the occurrence of a subduction process by itself (see e.g. Lavecchiaet alii, 1988; Lavecchiaet alii, 2003: Pauselliet alii, 2006). However, the strong reduction of lithospheric thickness detected by passive seismology (Calcagnile & Panza, 1980), reflected in the high surface heat flow, west to the Moho doubling, is consistent with the subduction retreat hypothesis. In our interpretation, the SW-dipping Adriatic crust envisaged in our model can be interpreted as a remnant of a lithospheric subduction of the thinned continental lithosphere, active until Middle Pleistocene times that is laterally continuous with the southernmost and still active subduction of the Ionian slab (Patacca & Scandone, 1989; Dellonget alii, 2019 and references therein). Above this retreating slab, the lithospheric mantle underlying the Tyrrhenian crust was progressively replaced by asthenospheric materials forming a new and shallower Tyrrhenian Moho (Scroccaet alii, 2005).

Possible developments of this research project may include:

  • - considering and incorporating a more detailed interpretation of the foreland fold-and-thrust belt, i.e. the upper continental crust of the Apennines, by integrating available geological, structural and geophysical data;

  • - analyze the distribution of instrumental seismicity along the transect, in order to define the rheological properties of the crust and of the upper mantle, in relation to the present-day geodynamic scenario.

In conclusion, our study highlights that gravity forward modelling can be used as a relatively quick and inexpensive method to refine and constrain the crustal structure of other sectors of the Apennine belt, as well as of other, similarly complex orogenic belts, where a range of possible solutions has been previously offered by different types of geophysical surveys.

This research is part of the MUSE-4D project, financed by MUR in the framework of PRIN-2017 (project #2017KT2MKE_003): U.R Responsible M.R. Barchi, National P.I. G. Lavecchia. We gratefully thank the Editor, C. Chiarabba, and the Reviewers, F. Speranza and N. Piana Agostinetti, whose constructive and thoughtful comments and suggestions greatly improved the manuscript. We thank Eni for providing the gravity data used in this work: the data were provided under a specific agreement with the Dipartimento di Fisica e Geologia of Perugia University and cannot be redistributed to third parties. The modelling was performed using the Oasis Montaj software, © Geosoft.