Gravity modeling of the Alpine lithosphere affected by magmatism based on seismic tomography

The Southern Alpine regions have been affected by several magmatic and volcanic events between the Paleozoic and the Tertiary. This activity has undoubtedly had an important effect on the density distribution and structural setting at lithosphere scale. Here the gravity field is used to create a 3D lithosphere density model on the base of a high-resolution seismic tomography model. The results of the gravity modelling demonstrate a highly complex density distribution in good 10 agreement with the different geological domains of the Alpine area represented by the European plate, the Adriatic plate and the Tyrrhenian basin. The Adriatic derived terrains (Southalpine and Austrolpine) of the Alps are typically denser (2850 kg m3), whilst the Alpine zone composed by European terrains provenance (Helvetic and Tauern window) presents lower density values (2750 kg m-3). Inside the Southalpine, south of the Dolomites, a well-known positive gravity anomaly is present and one of the aims of this work is to investigate the source of this anomaly that has not yet been explained. The modelled density 15 suggests that the anomaly is related to two different sources, the first involves the middle-crust below the gravity anomaly and is represented by localized mushroom-shaped bodies interpreted as a magmatic intrusion, while a second wider density anomaly affects the lower crust of the Southern Alpine realm and could correspond to a mafic and ultramafic magmatic underplating (gabbros and related cumulates) developed during Paleozoic extension.


Introduction 20
Very recently a 3D seismic tomography of the lithosphere in the Alps has been published (Kästle et al., 2018), which forms an ideal basis for a 3D density modeling of the globally available gravity field derived from the integration of satellite and terrestrial data. The 3D inversion of gravity data is an ill-posed problem, since more than one solution exists that reproduces the observations, and therefore constraints from a physical model or from other geophysical disciplines are necessary. The density model gives an added piece of information to the existing seismic velocity model, as systematic differences in the 25 relation between velocity and density give clues on the petrophysical rock composition. Moreover, the rheological rock parameters can be defined if both the seismic velocity and density are known. A particular focus of this work is on the lithospheric signatures of the magmatic activity that affected the Alps from Paleozoic to Paleocene, which conditioned the https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. mechanic characteristic of the crust during the Alpine phase, and in the phases preceding it Viganò et al., 2013;Zampieri, 1995). The presence of a gravity high located above the Venetian magmatic area comprising the Lessini 30 Mountains, Berici and Euganei hills, and Rovigo has been known for decades Ebbing et al., 2006;Zanolla et al., 2006), but an explanation of the high is lacking. The present model shall define the mass sources through a 3D density modeling that starts from the newly acquired seismic tomography with improved spatial resolution, thank's to the Alparray initiative.
The 3D modeling comprises the entire Alpine arc, and extends down to a depth of 200 km, the same depth of the seismic 35 tomography model. The conversion from seismic velocity to density is done with differentiated approaches, depending on depth and position in upper crust, lower crust or mantle. We use the most recent terrestrial-satellite derived gravity model, and apply the global correction for terrain, which leads to very different Bouguer gravity values compared to the classic local correction up to the 167 km Hayford radius. The global fields require special attention and we include a dedicated study in order to be used in the regional density modeling. The Alpine gravity field has been modelled before 40 Spooner et al., 2019), novel in our study is the use of the Bouguer field corrected for the global topography, which is more realistic than the 167 km classical Hayford radius reduction. A recent study (Spooner et al., 2019) used an approach in which density is constant for lateral crustal domains of typical extension of 80 km to 150 km, and the lithosphere is defined in layers of constant density. The layers have been defined with one layer for unconsolidated and consolidated sediments each, two layers for crystalline crust, one layer for subcrustal lithosphere, and one constant layer for the asthenosphere. In our approach 45 we do not pose any restrictions on the density variability, since we take full advantage of the high resolution seismic tomography model, and emphasis is put in our study to convert the seismic velocities to densities.
One of our central questions is whether we can demonstrate an alteration of the density in crust and upper mantle below the extended Venetian magmatic Alpine province and the older Athesian magmatism (Bellieni et al., 2010). Infact, during magmatic emplacement it can be expected that part of the molten crust is trapped in the lower crust. Only a portion (less than 50 10 wt%) of the entire molten rocks has erupted to the surface and less than 60 wt% of the whole melt can be emplaced as intrusions in the upper crust (De Min, personal comunication). Moreover, this crustal melting is often related to the emplacement of basic mantle derived magma at the Moho discontinuity (magmatic underplating), which provides the necessary temperature to induce a crustal melting.
The working hypothesis thus is, whether the above mentioned gravity high could be related to the magmatic events, and the 55 consequent crustal densification and/or juvenile crust formation in underplating.

Geological introduction
The study area comprises the Alpine arc, with particular focus on the magmatism of the southern Alps (Fig. 1). The southern Alpine domain represents the undeformed retro-wedge of an orogeny of double vergence (Dal Piaz et al., 2003;Zanchetta et al., 2012). To the north it is delimited by the dextral Periadriatic line, also known as Insubric line. The South-Alpine terrains 60 https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License.
formed the passive margin of the Gondwana Paleozoic shelf (Muttoni et al., 2009;Stampfli and Borel, 2002). The central South-Alpine realms are characterized by three principal tectonic systems: the Val Sugana line, the Gudicarie line, and the Schio-Vicenza line, all active in the late Permian. The late Paleozoic-Meosozic rifting created a system of NS oriented faults, which induced the separation of the Lombardic basin from the Veneto platform, documented by the variable thickness of the Triassic and Cretaceous sediments (Bertotti et al., 1993). From the late Cretaceous to the first Eocene, the extensional regime 65 is inverted due to the convergence between Gondwana and Laurasia, which generates compressional systems which initiate to the west of the Gudicarie, and later extends to the eastern Dolomites (Zanchetta et al., 2012). Later, the Gudicarie is affected by two compressive phases; in mid-late Miocene the ENE-WSW oriented thrusts of the Valsugana develop, which exhume part of the Hercynian basement. In late Miocene and Pliocene the NE-SW oriented thrusts of Bassano and Grappa-Montello developed, together with reactivation of faults in the Gudicarie Castellarin and Cantelli, 2000). From 70 end Paleozoic up to Tertiary different magmatic events have affected the South Alpine realm with intrusive and effusive products, that range from basic to acid. According to (Castellarin and Cantelli, 2000) the volcanic and magmatic volumes have contributed to considerably increase the rigidity of the crust, with implications on the tectonics and geodynamics.

Permian Magmatism
The South Alpine units have been affected by large volumes of magmatism in the Paleozoic (285-260 Ma), with volcanic and 75 plutonic products, of acid and basic type (Rottura et al., 1998). The thickness of the deposits has an average value of 900 m, and reaches up to 2000 m (Rottura et al., 1998), and are represented by the ignimbrites and rare flows of the Athesian Volcanic District. These volcanic rocks are bordered by several outcropping granitic intrusion and the whole frame has some similarities with the coeval Ivrea-Verbano system (De Min, personal comunication), which has been interpreted as the remnants of a supervolcano (Quick et al., 2009). The units that belong to this age are the porpyhric Athesian platform, intrusions of Cima D'Asta, 80 Bressanone-Chiusa, Ivigna, Monte Croce, Monte Sabion, Monte Luco. These magmatic products indicate important crustal contamination which could be triggered by previous subduction (Variscan?), as indicated by their geochemichal properties (Rottura et al., 1998). The magmatism, probably has originated in a post-Hercynic phase in a trans-tensive tectonic regime.
According to (Rottura et al., 1998) the magmatism was generated due to a lithospheric thinning followed by: an asthenospheric upwelling, 85 -consequent partial melting of the asthenospheric mantle for adiabatic decompression, -production of gabbro which, stored in underplating, has thermally perturbed the fertile crust, -consequent production of crustal melts which mixed with most evolved mantle derived magma.

Triassic Magmatism
The Triassic magmatism affected a continental area which at the time confined with the western Tethys, and comprised the 90 Southalpine, Dinaric, Transdanubian and Austral-alpine domains (Abbas et al., 2018). In the Dolomites the magmatic rocks are rare pyroclastics, dykes, sills, subaerial and marine volcanics, next to intrusives found at Predazzo and at Monti Monzoni.
The origin of the magmatism has been hypothesized as due to a plume, subduction, or extensional events at lithospheric scale (Abbas et al., 2018;Casetta et al., 2018). The magmatism could be associated to activation of regional strike-slip faults during the Ladinian (Abbas et al., 2018). In comparison with the Paleozoic magmatic event, the Triassic magmatic volume appears 95 to be more modest and characterized by a stronger abundance of basic and intermediate products. In the Dolomitic area these intrusions and vulcanites outcrop exactly inside the Permian magmatic structure. According to (De Min et al., 2020) this could be related to the presence, in that time, of a restitic crust (depleted during the Paleozoic event) which scarcely interacted with the mantle derived melts stored in underplating. These latter could therefore have evolved without heavily contributing to a crustal fusion. Smaller outcrops of this event are also present along the Sugana valley and the Recoaro-Pasubio area. 100

Tertiary magmatism
From late Paleocene to Oligocene, during 25 Myrs, the terrains south of the Trento platform were affected by a volcanism that extended over an area of 2000 km 2 (Milani et al., 1999). The principal units form the Lessini mountains, the Berici and the Euganei Hills. The province is characterized by basalts, basanites, transitional basalts and trachytes, which are divided in three categories according to the chemical composition, that ranges from alkaline to tholeiitic. The volcanic activity is due to an 105 extensional regime that accompanied the formation of a foreland basin. The formation is connected to WNW-SSE compression originated from the Adria-Europe convergence during the Alpine orogeny (Macera et al., 2008;Milani et al., 1999;Zampieri, 1995). During this same period, the calc-alcaline intrusives of Adamello and Mesino-Bragaglia to the west of the southern Judicarie were emplaced. These batoliths are composed of several intrusive units that formed in Eocenic, in an interval between 31 and 35 Ma. Some Rb/Sr datings arrive up to 42 Ma for the location Corno Alto (Martin and Macera, 2014). 110 3 Geophysical Data

Seismic tomography database
The seismic tomography database is the starting point for the 3D density modeling. The recent model of (Kästle et al., 2018) uses seismic stations distributed over the entire Alpine Arc, some stations in part of central Europe and many stations in Italy.
The tomography analyzes surface waves generated by earthquakes with periods between 8 s up to 250 s and ambient noise 115 observations in a period range of 4 s to 60 s. Given the sensitivity of the seismic periods, the surface waves resolve the deep structures, the ambient noise records the crustal structure and sedimentary basins, with an overlap of the two records for periods between 10 s and 35 s. The tomography inversion uses the (Ekström, 2011) model as a starting model for the inversion. The checkerboard test shows that the Alps and Italy are well resolved, except for the French territory, the Adriatic and Tyrrhenian seas. The Dinarides and Pannonian basin belong to the areas with less resolution. Some problems exist also in the resolution 120 of superficial structures in the eastern Po-basin, due to scarce availability of seismic stations. The model resolves the upper mantle well in the Alpine foreland of the Adria plate, of central Europe, of the Pannonian basin and the Ligurian Sea, with reduced resolution below the Alpine region. The model reports shear (Vs) velocities with a resolution in latitude of 0.1°, in https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. longitude of 0.146°, and a vertical resolution of 100, m for the upper 600 m, of 500 m between 1 km and 2.5 km, of 1 km between 3 km and 12 km, of 2 km between 14 km and 84 km, and 5 km between 85 km and 200 km. The present study selects 125 a portion of the model limited to longitudes 8 -14 ° and latitudes 43 -49 °.

Gravity field and topographic reduction
The gravity data are calculated from the global gravity model EIGEN6c4 (European Improved Gravity model of the Earth by New techniques, (Förste et al., 2014)). The model combines observations from different geodetic satellites, as LAGEOS-1/2 SLR, GRACE GPS-SST, GOCE and incorporates existing models for the high frequency part. It uses the DTU 2'x2' global 130 gravity anomaly grid from satellite altimetry (Andersen, 2010;Andersen et al., 2010) over the oceans and the EGM2008 model (Pavlis et al., 2012) over the continents. The model consists in a table of spherical harmonic coefficients complete up to degree and order 2190. For degrees above N=235 in continental areas, the model reproduces the field of EGM08, over oceanic areas the model DTU12.
It is customary to reduce the gravity field for the effect of the topography, for which models in spherical harmonic expansion 135 are available, as the model of Hirt and Rexer, 2015). The layers building topography have been divided into glaciers, oceans, lakes and solid crust, with their respective homogeneous densities equal to 917 kg m -3 , 1000 kg m -3 , 1030 kg m -3 , 2670 kg m -3 . We have used the particular model dv_ELL_RET2014PLUSGRS80, which includes the gravity of the ellipsoid. We refer to  for details on how the model is obtained.
Due to the fact that several options exist in how to define the functional to obtain the gravity value from the spherical harmonic 140 models of the earth gravity and the topographic effect, we briefly give here the fundamental definitions. Further details can be found for instance in (Torge, 2011;Torge and Müller, 2012) and (Barthelmes, 2013).
The Bouguer anomaly is defined as: Here the gravity effect of the topography ( ) is subtracted from the free-air anomaly ( ). The Free air anomaly according 145 to Molodensky (Torge, 2011) is defined as: with ( ) the observed value in point P, and ( ) is the gravity value of the reference ellipsoid calculated at the height of the telluroid. The gravity disturbance is defined as: Here the gravity field of the ellipsoid is subtracted from the observed gravity field, both calculated in the observation point P.
Here W is the earth gravity potential and U is the potential of the ellipsoid, and h is the geometric height. (Vaníček et al., 2004) show that due to the second term which enters the definition of the gravity anomaly, the gravity disturbance must be used when 155 deriving the gravity potential of the topography for the Bouguer correction. Infact, otherwise the second term of Eq. (5) would be applied twice. In spectral domain the derivatives translate to the following: where and are the Stokes coefficients of the gravity potential of either earth or topography, to which the Stokes 160 coefficients of the ellipsoid have been subtracted. The calculated gravity anomaly and disturbance are shown in Fig. 2, and it is seen that the gravity disturbance is systematically higher than the gravity anomaly, confirming (Vaníček et al., 2004).
Subtracting the gravity effect of topography calculated with the model RET2014  from the gravity disturbance derived from EIGEN6c4, the "Bouguer gravity disturbance" is obtained (Tenzer et al., 2019) (Fig. 3). We fulfill all calculations (Barthelmes, 2013) with the publicly available gravity calculation service of ICGEM (International Center for 165 Global Gravity Field Models; www.icgem.com). The result is a Bouguer map which is quite different from the one obtained from a local correction of the topography up to the Hayford radius of 167 km. Infact, the Bouguer values are -60 mGal over the Alpine arc, whereas the regional Bouguer values are classically below -180 mGal e.g. Spooner et al., 2019;Zanolla et al., 2006). Another area with clear negative values is close to the front of the Tuscan-Emilian Appennines, presumably due to the combined effect of sediments of the Po plain and the crustal thickening due to 170 stacking of sediment layers of the Appennines. The markedly different values obtained with the global reduction are inherent to the distant masses that are neglected in the regional reduction, but which are largely compensated by the isostatic crustal thickness variations (Szwillus et al., 2016). The 3D density modeling we intend to fulfill will be made on a window of less than 6° degrees latitude by 6° degrees in longitude, and a maximum depth of 200 km. The limit in lateral and depth extension of the model produces a bandlimited gravity field that can be modelled. It is therefore necessary to reduce the lowest degrees 175 of the gravity field of topography, since they introduce a field which is due to distant masses and are uncorrelated to the regional properties of the gravity field which is going to be modelled in the present study. We find that starting with degree and order N>10 the field of topography starts resembling the regional topography, and use this value as lower limit for the spherical harmonic expansion of the fields. This value agrees with the findings of (Steinberger et al., 2010), who define degree 10 the limit of lithospheric contributions, lower degrees having their origin in the deep mantle or lower. The band-limited 180 Bouguer map obtained for degrees 10 < N ≤ 2190 (EIGEN 6c4 and RET2014) has the same features as the non-band-limited map, with the difference that now the Boguer anomalies are more similar to those of the regional topographic reduction. Infact now the Bouguer disturbance is -160 mGal below the Alps, and -130 mGal below the Appennines.

Method
The aim of this work is to develop a three-dimensional density model of the Alpine region constrained by seismic tomography. 185 The method is developed in four steps: -Moho discontinuity investigation; -Seismic wave velocity-density conversion for mantle and crustal layers; -Gravimetric forward modelling and comparison between modelled results and the observed gravity data; -Inversion of the gravimetric residual and model density correction 190

Moho discontinuity depth definition
The Moho discontinuity depth was defined by studying the vertical variation of the tomographic seismic velocity. First of all, a preliminary velocity field analysis was fulfilled using other Moho models proposed by (Grad et al., 2009) (MohoGrad), (Tesauro et al., 2008) (EuCrust07) and (Laske et al., 2013) (Crust1.0). From this field analysis it was discovered that the mean S-wave velocity of the crust-mantle transition is 4.1 km s -1 with a standard deviation of 0.2 km s -1 (Table 1). This velocity 195 value is in agreement with the velocity change at Moho between 3.85 km s -1 and 4.48 km s -1 proposed in the 1D global model AK135 (Kennett et al., 1995). Therefore, on the basis of this statistics, a preliminary result was obtained researching the maximum vertical gradient of the velocity profile, within a velocity range between 3.8 km s -1 and 4.4 km s -1 , that represents the mean value plus and minus one standard deviation. To improve the stability and the reliability of this approach, another condition has been imposed on the velocity surrounding the Moho surface obtained from the maximum gradient analysis. The 200 S-wave velocity reported 4 km above and below the Moho surface must be respectively smaller and larger than the typical velocity of the lower crust and uppermost mantle. The mean P-wave velocity of the European lower crust is estimated to about 6.5 km s -1 , while in the upper mantle, the P-wave velocity of the peridotite is larger than 8 km s -1 (Christensen and Mooney, 1995). From this observation, the limits for the P-wave velocity above and below the Moho have been assigned to be of 6.8 km s -1 and 7.4 km s -1 . 205 The results show a good agreement with the Moho models proposed by the other authors (Fig. 4). Some differences are represented by an important crustal thickening located in the Ivrea-Verbano (D1) area and below the Istrian region (D2). Also the geometry and the depth found for the two main orogens (D3) are significantly different from the reference models.

Seismic wave velocity to density conversion for mantle and crustal layers
The Moho discontinuity allowed to divide the entire tomographic volume in crustal and mantle layer. For each layer, a different 210 approach and relations were used, in order to define the density values starting from the tomographic seismic velocity.

Crustal conversion
A huge set of relations were proposed by different authors in order to convert the seismic wave velocity in density (Brocher, 2005) and most of these equations were found to transform principally the compressional wave velocity (Vp) to density. For this reason, it has been necessary to convert the tomographic Vs in Vp through the "Brocher's regression fit" (Brocher, 2005), 215 a relation obtained from seismic, borehole and laboratory data for different kinds of lithology and valid for Vs between 0 and 4.5 km s -1 . Given the extreme variability of the crustal rocks present in the study area, a density estimation was carried out through two equations: the "Gardner's rule" by (Gardner et al., 1974) and the "Linear regression fit" proposed by (Christensen and Mooney, 1995). The first is an empirically derived equation that relates seismic P-wave velocity to the bulk density of the lithology. It is very popular in petroleum exploration because it can provide information about the sedimentary lithology 220 from interval velocities obtained from seismic data. (Christensen and Mooney, 1995) proposed a widely used linear density-Vp relation for crystalline rocks obtained from seismic refraction and laboratory data, in order to define the continental crust features as a function of the different geodynamic context. In this work, for distinguishing the sedimentary rocks velocity domain from the crystalline domain, and therefore which of the two relations to use, the velocity value at which the two curves intersect has been chosen (Fig. 5). The velocity is Vp = 6 km s -1 , as seen in (Brocher, 2005) and Fig. 5. 225

Mantle conversion
A different approach has been used for the velocity conversion in the upper-most mantle; in this case, the tomographic velocity model was compared with a Synthetic Upper Mantle model (SUM). The SUM was modelled through the open-source software Perple_X (Connolly, 2005), which allows the computation and mapping of the phase equilibria, rheological properties and density. 230 For this modelling the TECTON bulk rock composition model (Griffin et al., 2009) of the mantle has been used, which represents an estimate of the mean composition of a Neoproterozoic-Phanerozoic Sub-Continental Lithospheric Mantle (SCLM), gained from xenolith suites and peridotite massifs data. In Table A-1, the major elements of the TECTON composition are shown, of which we used the most abundant elements, up to a total content of 98.9 %, leaving out the minor elements with percentages less than 1 %. The software calculates the stable phases for a two-dimensional grid of temperature 235 (T) and pressure (P) conditions, typical for the mantle up to a depth of 200 km (Tmin=350 °C to Tmax=1600 °C, Pmin=0.1 GPa to Pmax=7 GPa). The software requires setting some calculation parameters, which are defined here. The thermodynamic model is the (Holland and Powell, 1998)  The result consists in the equilibrium phase diagram, which shows the fields of stability of different equilibrium mineral 245 assemblages for the adopted bulk-rock composition, depending on the P-T combinations. For each assemblage a number of parameters are given as well, from which we select Vs and density, shown in Fig. 6.
The conversion of the seismic velocity to density was made comparing the tomographic velocity at its given depth, with the synthetic velocity at the pressure corresponding to the depth. The pressure was calculated integrating the lithospheric density column from top to the corresponding depth. From the topography to the bottom of the crust, the crustal density model 250 calculated from the crustal seismic velocities was used. The resulting pressure at the base of the crust was used for converting the seismic velocity to density in the resolution cell below the Moho, and then iteratively downwards to the base of the model.
The temperature was left unconstrained and corresponds to the value that allows to match the tomographic and synthetic velocity at the given pressure. The above described density values define the starting model which is corrected through the gravity inversion discussed in the next chapters. 255

Gravimetric forward modelling and comparison between modelled results and the observed gravity data
In the previous section it has been argued how the wave velocity was transformed into density values. The next step was to discretize the entire model with tesseroids, mass elements defined in spherical geocentric coordinates. This task has been done building these elements as function of the original three-dimensional tomographic grid, from which we defined the density value of the tesseroids. In order to make the modelled gravity effect comparable with the observed gravity anomaly, the 260 obtained density model was reduced by the reference density model PREM ("Preliminary Reference Earth model"). This onedimensional reference model represents the mean physical features of the earth, including elastic parameters, seismic attenuation, density, pressure and the earth-radius dependent gravitational value. Being a global model, it was modified in order to obtain a one-dimensional model comparable with a continental crust, eliminating the first 3 km of seawater and lowering the Moho depth. All changes to the original PREM model have been done to maintain the pressure value at the base 265 of the model constant and we call this modified model co-PREM ( Fig. A1 in Appendix).
The Fig. 7 (a) shows the result of the gravitational forward modelling, developed through the software Tesseroids (Uieda et al., 2016). The gravity anomaly map presents two well defined negative anomalies induced by the crustal thickening below the orogens and the extensive positive anomaly produced by the Tyrrhenian crust. At large scale, the modelled anomaly reproduces in broad outline the observed data (Fig. 3) and from the difference between observed-modelled anomalies, we 270 obtained the residual gravity map (Fig. 7). The residual map shows that the modelled anomaly differs from the observed map by values that range from -40 mGal to +220 mGal. Furthermore, the residual outlines three regions characterized by different trends: the Tyrrhenian basin shows a positive residual around +80 mGal, Germany and part of Switzerland is characterized by negative values reaching peaks around -40 mGal, while the residuals are generally positive in all of the remainder study area, reaching values of +200 mGal. 275 https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License.

Inversion of the gravimetric residual and model density correction
Theoretically, the gravity residual map represents the extent to which the starting density model diverges from the one of the infinite possible density distributions that produced the measured gravity anomaly. In order to improve the reliability of the starting density model an inversion algorithm was developed based on the residual between observed and modelled data. As seen in the Fig. 8, the first step is to separate the long wavelength residual component from the short wavelengths, and to 280 estimate the different contributions of the mass sources that generate the residual signal.
The distinction between crustal and mantle contribution, which corresponds to the long and short wavelengths respectively, has been done by fitting a third-order polynomial surface with the values of the residual map (Fig. 9), and the choice of the order was made by a qualitative evaluation depending on the features and the large scale trend of the residual anomaly map.
The evaluation consisted in the observation that increasing the degree starting from zero added long wavelength signal to the 285 mantle component, which is uncorrelated to both topography and geology. Starting with degrees higher than 3 the polynomial started to show smaller scale features which correlated with the geometrical and physical properties of the crust. We therefore think that the third order polynome is qualitatively the correct choice to separate the mantle from the crustal component. The mantle component has a long period trend oriented NW-SE. The third order polynome allows to separate this very longwavelength signal from the crustal contribution. 290 Then, once we have obtained the two components, the algorithm provides the "linear iterative inversion" section, in which a rough first density correction is computed with the "Infinite slabs gravimetric effect" equation: where ∆ is the density correction, ∆ is the short or long gravimetric residual component, represents the thickness of the considered layer (crustal or mantle layer) and G is the gravitational constant. Then the density correction is applied to the 295 starting density model, the vertical distribution of the total mass correction has been assigned through the following equation: where _ and _ are respectively the new and the old density values, ∆ is the density correction, is the height of the considered layer (crust or mantle layer) and ∆ represent the height of the mass element belong the single density profile.
We apply the density correction, obtain a new mass distribution, and compute the forward modelling of the corrected density 300 model in the next step. The results of the latter operation is compared with the observed data producing, for each iteration, a new gravimetric residual map. This algorithm allows to define a density correction, which is applied to the three-dimensional starting model and allows to compare the gravity effect of this new mass distribution with the observed data.
https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. Figure 10 shows horizontal slices of the final density model at increasing depths ranging between 5 km and 40 km. Here we describe the main features.

Alpine region density distribution 305
In the first kilometres, up to 15 km depth ( Fig. 10 (a),(b),(c)), the low density values are due to the two main sedimentary basins situated in the fore-deeps of Alps and Apennines. In the European plate, the Molasse basin develops in W-E direction along the northern Alps front (LMB), and this sedimentary layer is wider and thinner in the western part whilst it straitens and 310 thickens in the eastern Bavarian part. In the Adriatic plate, the Po basin (PBL) is characterized by a thicker and less dense sedimentary layer than the Molasse basin. The Molasse basin is shallower (average thickness 5 km, maximum thickness 9 km), than the Po basin (average thickness 9 km, maximum thickness 13 km). As shown in Fig. 10  At greatest depth ( Fig. 10(g) and 10(h)), the modelled density features the Alps and Apennines crustal roots and the inner part (LIP) of the two belts is characterized by lower densities (2800-2850 kg m -3 ). Beyond the orogens crustal roots, the mantle is 335 marked by a density average around 3300 kg m -3 .
We have computed also the weighted average density for the crustal layer (Fig. 11). The average density distribution reproduces and synthesizes all the features seen above. In the European plate, different domains are present for the area close (2870 kg m -3 ) in the eastern part. The previously discussed W-E high-density ridge which develops from the Ivrea-Verbano anomaly to Schio-Vicenza fault system (RH), and that affects the Trento platform, is also seen in the average crustal density map.

Modelled density in the Venetian gravity anomaly area
Three representative transects across the main gravity anomaly structures of the central-eastern Southern Alps are illustrated in Fig. 12, and discussed in the following.
The A section (Fig. 12) shows the density distribution along the main direction of the positive gravity anomaly, parallel to the Schio-Vicenza fault. The low density in the southern part of the profile, represents the sediment layer that reaches a thickness 355 of around 13 km in the deepest part (longitude 12.5° E in the profile). As seen before, between the South-Giudicarie fault (SGF) and the buried Appennines thrust front, there is again the high-density plateau (RH) situated between 25km and the Moho discontinuity, with a density of 2900-3000 kg m -3 . In the middle crust, two limited mushroom-shaped bodies are lying between 10 and 25 km, the first (MB1) is located below the Venetian Volcanic Province (VVP) and above the high-density plateau (RH) close to latitude 11.5° N, while the second (MB2) is situated below an area between the south Giudicarie fault 360 (SGF) and the Periadriatic Line (PL) corresponding to the outcropping Adamello batholith (AB).
The NE-SW oriented section (Fig. 12 -B section), is orthogonal to the well-developed sedimentary wedge of the Alpine and Apennines foredeep (PB), which is 10 km thick in the deepest part and has density values between 2350 up to 2550 kg m -3 .
Using the 2700 kg m -3 contour line as a marker of the top of the Mesozoic carbonate platform, it is possible to note that the contour line depth increases between latitude 46° N and latitude 44° N. The 2700 kg m -3 contour line reaches the surface in 365 https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. front of the Valsugana thrust system (at around 46° N latitude) (Fantoni and Franciosi, 2010). North of the Valsugana thrust front (VTF), corresponding to the southern alpine belt, the model shows density values of 2750 kg m -3 for the upper crust and 2800 up to 2900 kg m -3 for the middle-lower crust. In the central part of the section the model reports a very high-density plateau (RH) situated in the lower crust, that develops from the Apenninic buried thrust front (ATF) up to the outcropping south-vergent thrust (VTF). 370 The last section (close to W-E) crosses the Ivrea-Verbano area (Fig. 12 -section C), and then runs along the Po sediment basin. Lateral density variations of the upper crust are well correlated with the different basin domains represented from west to east by the Lombardian basin (LB), the Trento platform (TP) and the Belluno basin (BB). The mid-crust shows high-density (e.g. greater than 2800 kg m -3 ) values already at a depth of 20 km, and this value distribution corresponds to the previously discussed high-density ridge (RH) seen in the average density map (Fig. 11). 375

Discussion
Several geophysical models have been developed for the Alpine lithosphere in the past, especially using seismic (Gebrande et al., 2006;Kästle et al., 2018) and gravity data as modelling constraints Spooner et al., 2019). In this work for the first time, a high resolution three-dimensional seismic tomography model has been adopted as starting model for the computation of a 3D density model using the algorithm discussed before. The results show a highly complex density 380 distribution in good agreement with the different geological domains of the Alpine area represented by the European plate, the Adriatic plate and the Tyrrhenian basin. Notwithstanding the offsets that can occur between the outcropping geology, and the intracrustal tectonic domains, we find that the average crustal density map correlates with the different tectonic realms.
The Adriatic derived terrains (Southalpine and Austrolpine) of the Alps are typically denser (2850 kg m -3 ), whilst the Alpine zone composed by European terrains provenance (Helvetic and Tauern window) presents lower density values (2750 kg m -385 3 ). This density relation has been also noted by (Spooner et al., 2019), indicating that this physical parameter could be used as a marker for the different tectonic domains.
The South-Alpine and Austroalpine provinces, were the locus of sporadic, yet occasionally widespread magmatism. The magmatic occurrences are associated with strongly variable tectonic regimes, which range from clearly collision-related (the Austroalpine Paleogene plutonism, i.e., Vedrette di Ries/Riesen Ferner), to late-or slightly post-orogenic associated with basin 390 formations shortly after the main Variscan orogeny (Permian magmatism, e.g., the Athesian Volcanism and Ivrea Verbano and Serie dei Laghi Group), to clearly post-orogenic (Triassic, e.g., the Predazzo-Monzoni complexes), and finally to an almost within-plate environment, peripheral to the Alpine belt (the Oligocene Veneto Volcanic province) (Bellieni et al., 2010).
These magmatic activities have undoubtedly changed the thermal setting at the time and, after cooling, increased the crustal rigidity especially in the Venetian Alps (Dolomites) which appears today the least deformed and shortened in the whole 395 Southern Alps domain . https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License.
In this domain, south of the Dolomites, the well-known positive Bouguer gravity anomaly is present, located between the towns of Verona and Vicenza, which covers the Venetian Volcanic Province (VVP). The short-wavelength content of this gravity high can reasonably be explained by the presence of a shallow Mesozoic carbonate platform (Trento Platform) and basement which hosted several volcanic bodies, as demonstrated by field evidence and the intra-sedimentary magnetic bodies 400 proposed by (Cassano et al., 1986) on the base of magnetic and well data. According to this interpretation, the isolated mushroom-shaped body modelled in the middle-upper crust could be related to the Triassic and Cenozoic magmatic activity.
These intrusions develop below the Venetian Volcanic Province and could correspond to the Predazzo-Monzoni complex and the Adamello and Mesino-Bragaglia plutons emplaced respectively in the Dolomites and westward of the South-Giudicarie fault. However, the remainder part of this signal corresponds to signal generated by the lower crustal source component, which 405 gives origin to the long wavelength of the positive Bouguer anomaly, as also proposed by (Zanolla et al., 2006). During Permian, the Austroalpine and Southalpine units were affected by magmatism, high-temperature/low-pressure metamorphism (HT-LP) and extensional tectonics. Intra-continental basins hosting Permian volcanic products have been interpreted as developed either in a late-collisional strike-slip or in a continental rifting setting (Spalla et al., 2014;Spiess et al., 2010) and references therein. As regards the geodynamic context, all evidence suggests that this magmatic pulse was induced by post-410 Variscan lithospheric attenuation that produced crust-mantle detachment and large ultra-mafic crustal underplating, together with the development of a Permian basin (Dal Piaz and Martin, 1998). For a fuller understanding of this process, we may refer to the modern model tested in the Western Southern Alps, in the Ivrea-Verbano zone, where the entire crustal section of an underplated and metamorphosed Hercynian basement crops out, overlain by volcanic and sediment formations (Quick et al., 1995(Quick et al., , 2009Sinigoi et al., 2011). In the central Southern Alps the Permian activity is represented by the porphyry plateau of 415 the Bolzano province (Athesian plateau) and the numerous granitic-granodioritic intrusions. Considering the large scale of this process, it is reasonable to suppose the scattered presence of mafic underplating in the whole Alpine realm, confirmed by the Permian-Triassic gabbro bodies and related HT-LP metamorphism emplaced in the Austro-Alpine and South-Alpine continental crust (Spalla et al., 2014). According to this hypothesis, the high-density values modelled for the Adriatic lower crust (high-density ridge HR, Fig. 10 (d) -(f)) are consistent with the densities of a basaltic mantle-derived rock, hence this 420 anomalous lower crust could be interpreted as mafic underplating emplaced between the crustal-mantle transition during Permo-Triassic post-collisional lithospheric extension.
Historically, the deep source of the Vicenza-Verona positive gravity anomaly has been explained by thinned crust of 25-27 km ( (Cassinis, 2006;Finetti, 2005;Gebrande et al., 2006;Scarascia and Cassinis, 1997). This approach has been promoting the development of several Moho models of the Alps that invoke a relatively shallow mantle situated in the central-southern 425 alps, especially for the gravity-modelling derived Moho models (Braitenberg et al., 2002;Ebbing et al., 2001). Recently, other seismic methodologies such as local earthquake tomography, receiver functions and ambient noise tomography, have allowed to develop 3-D high-quality crustal models of the Alps (Kästle et al., 2018;Spada et al., 2013), and this family of Moho models point to a crustal thickness of 30-35 km in the Vicenza-Verona positive gravity anomaly region, corroborating the possibility that the gravimetric anomaly is related to an intra-crustal denser source (mafic underplating). Also (Mueller and Talwani,430 https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. 1971) in their Alpine crustal model, proposed high density distribution for the South Alpine lower crust, that connects the Venetian gravity high with the Ivrea-Verbano zone.
Unfortunately, the kinematics of magmatic underplating has been described mainly by petrophysical models explaining the genesis of the outcropping rocks, while the available geophysical images of underplated material remain relatively sparse and confined to specific tectonic environments (Thybo and Artemieva, 2013). Our work adds a geophysical model based on the 435 gravity field to demonstrate the existence of underplating to an area related to an orogen. Underplating from a gravimetrc analysis had been postulated for instance also in the context of the Parana' large igneous province  and the West Siberian large igneous province (Braitenberg and Ebbing, 2009).

Conclusion
The scope of this work is to study the density distribution of the Alpine lithosphere, with the specific target of characterising 440 the sources of the well-known Venetian gravity anomaly high. A new approach based on three-dimensional seismic tomography as starting model for the computation of a 3D density model has been developed. We propose an innovative approach which includes several numerical algorithms and that provides a Moho definition, a velocity-density conversion, gravimetrical forward modelling and finally a linear inversion based on the gravimetric residual map. The results for the Alpine structure highlight that all terrains of Adriatic provenance (Austro-Alpine and South-Alpine) are denser (2850 kg m -3 ) than the 445 European terrains (2750 kg m -3 ), suggesting that this physical parameter could be used as a marker to characterize the different tectonic domains involved in the Alpine regions (Spooner et al., 2019). We also considered the alternative source models for this gravity anomaly, in which different authors provide a stretched and thinned Adriatic crust located between the Giudicarie thrust front and the Schio-Vicenza fault, but today several Moho models of the Alps exist, that refute this model.
The approach we use, could represent a generally applicable method to achieve insights about large-scale crustal density 460 distributions, starting from seismic tomography, and it allowed us to carry out a high-resolution density model for the Alpine https://doi.org/10.5194/se-2020-101 Preprint. Discussion started: 30 June 2020 c Author(s) 2020. CC BY 4.0 License. region. The expected improvement of more robust and high-resolution seismic tomography models, as well as the develop of a non-linear inversion algorithm, could enhance the reliability of the density models in the future.   Figure A1 Density of the PREM model and the modified coPREM model used as reference in the modeling procedure.

Data Availability
The 3D density model is available in the publicly accessible repository zenodo.org at link https://www.zenodo.org/record/3885404#.Xt5IHzozbIV.

Author Contributions
Both authors contributed to the conception and design of the study, manuscript writing, and approved the submitted version. 495 Davide Tadiello was in charge of all computations.    (Pfiffner, 2015) and main magmatic and volcanic bodies (Bigi et al., 1990a(Bigi et al., , 1990bSchuster and Stüwe, 2008)