# New gravity data and 3-D density model constraints on the Ivrea Geophysical Body (Western Alps)

M Scarponi, G. Hetenyi, T Berthet, L. Baron, P. Manzotti, Benoît Petri, M.

Pistone, O Müntener

► To cite this version:

M Scarponi, G. Hetenyi, T Berthet, L. Baron, P. Manzotti, et al. New gravity data and 3-D density model constraints on the Ivrea Geophysical Body (Western Alps). Geophysical Journal International, 2020, 222 (3), pp.1977-1991. 10.1093/gji/ggaa263. hal-02912137

HAL Id: hal-02912137

https://hal.science/hal-02912137

Submitted on 20 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# New gravity data and 3-D density model constraints on the Ivrea Geophysical Body (Western Alps)

M. Scarponi,<fnref n="1" /> G. Hetényi,<fnref n="1" />,<fnref n="2" />,<fnref n="1" /> T. Berthet,<fnref n="2" /> L. Baron,<fnref n="1" /> P. Manzotti,<fnref n="1" />,<fnref n="3" /> B. Petri,<fnref n="1" />,<fnref n="4" />
M. Pistone<fnref n="1" />,<fnref n="5" /> and O. Müntener<fnref n="1" />

Accepted 2020 May 24. Received 2020 April 29; in original form 2020 February 4

## SUMMARY

We provide a high-resolution image of the Ivrea Geophysical Body (IGB) in the Western Alps with new gravity data and 3-D density modelling, integrated with surface geological observations and laboratory analyses of rock properties. The IGB is a sliver of Adriatic lower lithosphere that is located at shallow depths along the inner arc of the Western Alps, and associated with dense rocks that are exposed in the Ivrea-Verbano Zone (IVZ). The IGB is known for its high seismic velocity anomaly at shallow crustal depths and a pronounced positive gravity anomaly. Here, we investigate the IGB at a finer spatial scale, merging geophysical and geological observations. We compile existing gravity data and we add 207 new relative gravity measurements, approaching an optimal spatial coverage of 1 data point per 4–9 km<fnref n="2" /> across the IVZ. A compilation of tectonic maps and rock laboratory analyses together with a mineral properties database is used to produce a novel surface rock-density map of the IVZ. The density map is incorporated into the gravity anomaly computation routine, from which we defined the Niggli gravity anomaly. This accounts for Bouguer Plate and terrain correction, both considering the *in situ* surface rock densities, deviating from the 2670 kg m⁻<fnref n="3" /> value commonly used in such computations. We then develop a 3-D single-interface crustal density model, which represents the density distribution of the IGB, including the above Niggli-correction. We retrieve an optimal fit to the observations by using a 400 kg m⁻<fnref n="3" /> density contrast across the model interface, which reaches as shallow as 1 km depth below sea level. The model sensitivity tests suggest that the ~300–500 kg m⁻<fnref n="3" /> density contrast range is still plausible, and consequently locates the shallowest parts of the interface at 0 km and at 2 km depth below sea level, for the lowest and the highest density contrast, respectively. The former model requires a sharp density discontinuity, the latter may feature a vertical transition of densities on the order of few kilometres. Compared with previous studies, the model geometry reaches shallower depths and suggests that the width of the anomaly is larger, ~20 km in west–east direction and steeply E–SE dipping. Regarding the possible rock types composing the IGB, both regional geology and standard background crustal structure considerations are taken into account. These exclude both felsic rocks and high-pressure metamorphic rocks as suitable candidates, and point towards ultramafic or mantle peridotite type rocks composing the bulk of the IGB.

**Key words:** Composition and structure of the continental crust; Gravity anomalies and Earth structure; Europe; Numerical modelling; Continental tectonics: compressional; Crustal structure.

<footnote n="1">¹Institute of Earth Sciences, University of Lausanne (UNIL), Batiment Geopolis, Quartier UNIL-Mouline, 1015, Lausanne, Switzerland. E-mail: matteo.scarponi@unil.ch</footnote>

<footnote n="2">²Earth Sciences Department, University of Uppsala, Geocentrum, Villavägen 16, 752 36, Uppsala, Sweden.</footnote>

<footnote n="3">³Department of Geological Sciences, Stockholm University, SE-106 91, Stockholm, Sweden.</footnote>

<footnote n="4">⁴Université de Strasbourg, CNRS, IPGS UMR 7516, F-67000 Strasbourg, France.</footnote>

<footnote n="5">⁵Department of Geology, University of Georgia (UGA), Franklin College of Arts and Sciences, Geography-Geology Building, 210 Field Street, Athens, GA 30602-2501, USA</footnote>
# 1 INTRODUCTION

The present-day architecture of the Alps is the result of a succession of lithospheric-scale processes, involving the collision between the two major European (to the North) and African Plate (to the South), since ~35 Ma (Handy *et al.* 2015), following the opening and the subsequent subduction of the Piemonte-Liguria ocean (e.g. Handy *et al.* 2010). The Europe-Adria collision, with Adria as the largest among the microplates involved in the collision, provided most of the rock mass of the current Alps (Schmid *et al.* 2004). The resulting geometry of the Alps is structurally complex both at the surface and at depth, and has been continuously investigated by a broad spectrum of geological and geophysical methods for more than a century. One of the outstanding features of the Alps are the Ivrea-Verbano Zone (IVZ) and the Ivrea Geophysical Body (IGB).

The IGB is considered as a sliver of Adriatic lower lithosphere, which was emplaced in the upper crust along the inner arc of the Western Alpine domain and is located at unusually shallow depths (e.g. Schmid *et al.* 2017). This work focuses on its northeastern portion, where the IGB is partially outcropping at the surface and presents a variety of middle-to-lower crustal and mantle composition rocks, across the geologically well-known IVZ (Fig. 1). This almost-unique complex can be regarded as a nearly complete cross-section of the continental crust (Fountain 1976) and it has been the subject of a variety of laboratory sample analyses for rock physical properties (e.g. Fountain 1976; Khazanehdari *et al.* 2000), tectonic mapping (e.g. Schmid *et al.* 2004; Brack *et al.* 2010; Petri *et al.* 2019) and petrological studies aiming at deciphering magmatic processes (e.g. Rivalenti *et al.* 1995; Sinigoi *et al.* 2011; Karakas *et al.* 2019) and its emplacement (Quick & Denlinger 1993; Quick et al. 1994, 2003). This exposed cross-section is delimited on its northwestern side by the Insubric Line (Fig. 1), that is the western end of the largest fault system in the Alps, marking the boundary between the Adriatic Plate and the orogenic wedge. From a geological point of view, it separates the Southern Alps—with a very low grade of Alpine metamorphism—from the Western, Central and Eastern Alps—with moderate to high grade Alpine metamorphism (e.g. Schmid *et al.* 1989; Bousquet *et al.* 2012). Several contributions support the idea of this fault being a subvertical to vertical feature, cutting across our study area, both from structural evidences (e.g. Schmid *et al.* 1987, 1989; Berger *et al.* 2012) and geophysical investigations (e.g. Schmid & Kissling 2000).

In this paper, we focus on the IGB structure beneath the IVZ, which is characterized by two main geophysical anomalies: a pronounced positive gravity anomaly and a high seismic P-wave velocity anomaly. The IGB density structure and the associated gravity anomaly have been addressed by several contributions in the literature (e.g. Niggli 1947; Vecchia 1968; Marson *et al.* 1994), leading to the production of different crustal density models. These range from continuous bodies with a constant density and a sharp density-contrast interface (e.g. Berckhemer 1968; Kissling *et al.* 1984) to numerous and varying-size block-assembled models, using different density contrasts (e.g. Bürki 1990; Bayer *et al.* 1989; Rey 1990). The profile presented by Berckhemer (1968) shows an anomalous body similar in shape to a bird with a head, a neck and a beak, from which the term bird's head is sometimes used in the literature.

A certain number of Bouguer anomaly maps have also been produced from field-data, highlighting a broad agreement on the order of magnitude of the anomaly's amplitude and on the main location of the pronounced positive anomaly, associated with the IGB (Fig. 2). Nevertheless, differences between these maps arise when it comes to compare iso-anomaly contour shapes and locations, crest

lines and local maxima, most likely because of the heterogeneous and sparse data coverage and interpolation artefacts (Fig. 2).

The main constraints on the IGB structure and geometry were first given by seismic refraction experiments across and along the Western Alps (e.g. Berckhemer 1968; Ansorge *et al.* 1979). These seismic investigations provided the first constraints on the gravity models of the shallow IGB structure, suggesting the existence of seismic interfaces as shallow as few kilometres. Later, the IGB has been detected with reflection seismics as well (e.g. Bayer *et al.* 1989; Thouvenot *et al.* 1990), highlighting shallow seismic features that served as guidance for the further development of 2-D density models with gravity data (e.g. Rey 1990; Bayer *et al.* 1989). More recently, important insights on the IGB structure and bulk properties came from local earthquake tomography, showing anomalous P-wave seismic velocities in the range of 7–8 km s⁻¹ at middle-crustal depth (e.g. Solarino *et al.* 2018; Diehl *et al.* 2009).

We also base our study on a seismically constrained Moho map over the whole region (Spada *et al.* 2013) which, together with the latest seismic tomography experiments, offers a suitable and homogeneously defined initial structure for the development of a higher-resolution 3-D crustal density model. Indeed, while tomographic images reveal a broad-scale connection between the top of the bird's head imaged by active seismics and the European Moho at more than 40 km depth and shallower on the Adriatic side at ca. 35 km, the size and geometry of the IGB at depth still remain uncertain. Here, we advocate for the need for higher-resolution imaging to investigate the shallow upper-crustal density-structure beneath the IVZ and to look at the IGB in more detail, using gravimetry and geological field surveys. This is especially important in the perspective of understanding how seemingly dense rocks are emplaced at shallow depth in an orogen, and in particular in light of the forthcoming scientific drillings in the area (Pistone *et al.* 2017). For the latter, progress towards an exhaustive investigation of the IGB and new insights on the IVZ crustal structure, by means of a complete and up to date gravity database of the area, are a strategic step.

For this purpose, we have collected new and compiled old gravity data, together with surface rock density information. We developed a new 3-D density model, defined by a single density-contrast interface within the crust, by iteratively minimizing the misfit between synthetic data and gravity anomaly observations. We define a new gravity anomaly correction, the Niggli correction, in order to account for the surface density deviations from the otherwise constant correction density 2670 kg m⁻³ and we incorporate this correction into the computations. The sensitivity of model geometry to the crustal density contrast is then investigated in order to explore the acceptable range of model geometries and density contrasts, which may add constraints on the nature of the rocks composing the IGB.

# 2 GRAVITY DATA ACQUISITION

A new gravity database was compiled by merging existing gravity data sets with our own newly acquired data set, aiming at a homogeneous spatial coverage of one gravity point per 4 to 9 km². This allows for a detailed mapping of the regional anomaly and brings the investigation scale closer to the available geological observations within the IVZ, as represented by the various tectonic units adopted from Petri *et al.* (2019) and from Manzotti *et al.* (2014) for what lies to the northwest of the Insubric Line (Fig. 4a). In the following, the data collection procedure and pre-processing practices (from raw data to absolute gravity values) are presented.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Topographic map of the Western European Alps (see inset for broader location). The red box is the region of the interest for this study, focusing on the geologically well-known Ivrea-Verbano Zone, indicated in yellow (from Petri *et al.* 2019), and bounded to the west by the Insubric Line (shown only within the red box, in cyan). Figs 2–5 and 7 show the same area than the red box.</figcaption></figure>

## 2.1 Existing data compilation

Existing gravity data have been compiled from several sources. The Swiss Federal Office of Topography (Swisstopo: https://www.swisstopo.admin.ch/) provided one data set covering some of the North-Italian valleys and a data set covering the fraction of Swiss territory of our interest. Additional data has been obtained from the Istituto Nazionale di Oceanografia e di Geofisica Sperimentale (OGS: https://www.inogs.it/).

This preliminary compilation served as basis for field-campaign planning across our study area, which presents several areas of limited accessibility due to their rough topography and limited road network in the high mountain areas. Therefore, the measurement locations have been defined in order to refine the previously available data set and to obtain a data coverage as homogeneous as possible. The new data is available upon request and will be publicly released by end of 2020.

## 2.2 New gravity data collection

We reached the majority of the locations by road vehicle, on foot (0.5–3 hr hiking to the target location for a dozen of points) and by helicopter, the latter to reach the most remote areas. By the end of our field efforts, requiring *ca.* 36 d in the field, in September 2017, June 2018 and September 2019, respectively, we have made measurements at 207 new gravity points, which has drastically improved the spatial coverage map (Fig. 3). The average distance of a newly measured gravity point to their respective nearest neighbour—in the previous data set—is 2.1 km, while it is 2.9 km to nearest 4 existing points, on average.

All measurements were organized into daily loops, with each day or half-day of field campaign starting and finishing at the exact same location. This practice allows to estimate the gravimeter instrumental drift, a mechanically inherent property of the instrument which

needs to be estimated and then subtracted from each measurement, assuming that the drift is linear during the loop.

At each site, a Scintrex CG-5 relative gravimeter was used to measure the local gravity variation, with respect to the previous site or the starting point of the associated daily loop. In this way, relative gravity variations are recorded at each step along the loop.

As a result, we obtained a relative but self-consistent network, with fixed gravity variations within the network, that needs to be referenced to an absolute gravity point. We connected our network to the established gravity reference points at the Simplon Pass (Switzerland) and at the Polytechnic of Milan (Italy), where the very first loop of the campaign started. After measuring gravity in the absolute reference point in Milan, additional reference points for further loops were defined progressively, allowing for an expansion of the network loop by loop.

## 2.3 Point location and elevation estimate

For each gravity point, the elevation estimate at the measurement site is of fundamental importance for the subsequent data processing and modelling, as each metre of vertical displacement upward in free atmosphere corresponds to *ca.* 0.3086 mGal decrease in gravity. For the purpose of estimating elevation, a stand-alone GNSS (Global Navigation Satellite System) receiver antenna (TopCon, www.topconpositioning.com) has been used at each site to record satellite signals for *ca.* 20 min at 1 Hz sampling rate. Repeated measurements at the same location show that this setup allows a 0.5 m precision on the elevation estimate, obtained via PPP (Precise Point Positioning) processing from Natural Resources Canada (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php). This open-access service provides longitude and latitude coordinates in the WGS84 reference system as well. The ellipsoidal height, given by the PPP procedure, is then converted to the physical height, by
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Compilation of previous gravimetric works, as referenced in the legend, having mapped Bouguer anomalies in the area. All document a pronounced positive anomaly following the geological Ivrea-Verbano Zone (IVZ, in yellow). (a) Location of the crest lines and the local maxima. (b) The 0 mGal anomaly contour lines. Panels (a) and (b) highlight discrepancies between the maps.</figcaption></figure>

interpolating the geoid map of Italy (Corchete 2010) at our measurement locations. To ensure the best precision, we measured the gravimeter and GPS antenna elevations from the ground with 1 cm precision, and referenced our gravity data points homogeneously. The elevation of the newly acquired points ranges from 117 to 2009 m.

## 2.4 Gravity data pre-processing

Processed GNSS coordinates and raw gravity data, organized into ordered loops, were processed via the GRAVPROCESS software

(Cattin *et al.* 2015). This tool calculates point-wise absolute gravity values, free-air gravity anomaly and complete Bouguer gravity anomaly values, by applying corrections for tides, atmospheric pressure, instrumental drift, latitude, elevation as well as the Bouguer Plate and terrain corrections.

The terrain correction is especially important in mountain areas and contributes to the signal by several mGals (1 mGal = $10^{-5}$ m s$^{-2}$); it is thus fundamental for a correct estimation of the Bouguer gravity anomaly. In our study, we used the SRTM digital elevation model (Shuttle Radar Satellite Mission, Farr *et al.* 2007) for the terrain correction in the whole region, together with the measured site elevations. The gravitational effect of topography has
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Gravity data coverage around the IVZ (green contour), obtained by merging previous gravity databases and newly collected gravity measurements. The black crosses represent the 207 new relative gravity points from this project. The colour map indicates the distance to the nearest gravity point, white-to-grey representing our intended optimal coverage.</figcaption></figure>

been computed within a 166.917 km radius circle (second Hayford
zone) around each gravity point, by interpolating the digital eleva-
tion model over a 2-D mesh grid, whose resolution is decreasing
with the distance from the gravity data point. The mesh geometry
parameters in GRAVPROCESS are user-defined (see the *Mesh* box
in the software interface) and were fixed at 0.1 km, 166.917 km and
0.3 during our calculations, representing respectively the maximum
grid resolution, the maximum distance from the gravity point for
topographic effects and the resolution decay with the distance from
the gravity point. We refer to the GRAVPROCESS user manual for
further details. Statistics on the quality of the newly processed data
indicate data point uncertainties well below 1 mGal: the mean value
of standard deviation of repeated measurements at a given point is
0.015 mGal (min. 0.005, max. 0.14), the daily instrument drift mean
is 0.007 mGal (min 0.001, max 0.13) and uncertainties on the final
anomaly products range between 0.18 and 0.27 mGal.

We then processed uniformly a unique gravity database, made
up of both compiled and newly measured data. The gravity data
from various sources match well, as data represented along various
profiles crossing the entire study area do not show visible systematic
shifts among the different data sources. We cannot exclude the
existence of smaller-amplitude systematic errors or fluctuations, on
the order of one or few mGals; nevertheless, their effect would be
lower than the variations introduced by the 5-km projection of data
along our profiles, necessary for the modelling.

3 GRAVITY DATA PRODUCTS

The IVZ presents a variety of outcropping dense rocks at the surface, typical of middle-to-lower crust (e.g. Fountain *et al.* 1990). Their respective densities clearly differ from the 2670 kg m⁻³ correction density used in standard gravity anomaly computations (Hinze

2003). Hence, we define and apply the *Niggli correction* during the
gravity anomaly computation. We named the corrected anomaly the
*Niggli anomaly* in honour of one of the first authors who studied the
relationship between gravity anomaly and dense-rock observations
in the IVZ (Niggli 1947).

This correction accounts for surface rock-density deviations from
the 2670 kg m⁻³ reference value, and hence includes geological
surface observations which can be local but of significant amplitude.
As a first step in applying this correction, we compile a surface
density map, which is in turn incorporated into the GRAVPROCESS
terrain correction computation, as described below.

3.1 Surface rock-density map

Our study area encompasses the western end of the Southern Alps
and the adjacent units of the Western Alps, separated by the Insubric
Line (Fig. 1). The studied portion of the Western Alps is located
in an area where the nappe stack involves, from top to bottom,
(i) continental units deriving from the Adria palaeomargin (i.e.
the Sesia-Dent Blanche nappes); (ii) remnants of the Mesozoic
Piemonte-Liguria ocean (i.e. the Zermatt-Saas and Combin Zone)
and (iii) slices of the Iberian-European basement (e.g. Gran San
Bernard, Gran Paradiso and Monte Rosa massifs, e.g. Schmid *et al.*
2004; Manzotti *et al.* 2014). Southeast of the nappe stack (i.e.
southeast of the Insubric Line boundary, Fig. 4a), the Southalpine
domain is represented by the IVZ and the Serie dei Laghi (e.g. Zingg
1983; Boriani 1971). The former is a piece of Permian middle to
lower Adriatic crust whereas the latter is a piece of Permian upper
Adriatic crust covered by Mesozoic sediments; both belonged to the
upper plate during the Alpine orogeny.

In order to produce a surface density map of this area, we com-
piled information from different sources and followed a systematic
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Geological units we compiled for our study area, following Bigi *et al.* (1983), Manzotti *et al.* (2014) and Petri *et al.* (2019). Each compiled unit has been associated with the dominant rock composition. The legend on the right refers to each unit type as listed in Table 1. (b) Surface rock-density map of the Ivrea-Verbano Zone, associating a constant value with every tectonic unit as mapped in Petri *et al.* (2019) and Bigi *et al.* (1983). The density of each unit was chosen as a mean value, on the basis of their mineralogical composition, field observations and laboratory rock-sample analyses (Khazanehdari *et al.* 2000; Brack *et al.* 2010; Petri *et al.* 2019). Mean unit density is shown in Table 1 and obtained as a result of the associated unit composition.</figcaption></figure>

and consistent approach, both for what lies Southeast of the Insubric Line (the IVZ, its subunits and the Serie dei Laghi, Petri *et al.* 2019) and for what lies to the northwest (Manzotti *et al.* 2014). Finally, the compiled units have been merged into a single map (Fig. 4a).

Concerning the IVZ, we follow the tectonic unit boundaries from Petri *et al.* 2019. We have associated each unit with a reference petrological signature (Fig. 4a) by comparing with previous maps (Brack *et al.* 2010) and with laboratory rock analyses of samples collected within the units themselves (Khazanehdari *et al.* 2000; Table 1).

In order to follow a consistent procedure for the density estimation of these and subsequent units, we compiled the published petrophysical data of IVZ rock samples with specific mineral model proportions (Table 1) from Khazanehdari *et al.* (2000) and we referred to a mineral property database (www.mindat.org) to assign individual mineral densities (Table 1). For minerals with end-member compositions, we considered weighted averaged densities according to expected mineral abundances in the area. Finally, we computed rock densities according to the mineral densities and the given mineral abundances (Table 1, Fig. 4b).
Evaluating the rock densities from mineral density and compositions is likely to be an upper bound as they are not accounting for possible presence of distributed fractures and porosity, yet rock sample porosity measurements appeared to be negligible in this area (<1 per cent for those reported by Khazanehdari et al. 2000).

Some remarkable deviations from the reference 2670 kg m⁻³ arise in the IVZ, especially in the areas of exposed lower crust, including upper mantle rocks that were integrated as lenses in the lower crust (e.g. ~3370 kg m⁻³ for the Balmuccia and Premosello peridotites and ~3320 kg m⁻³ for the Finero peridotite).

Concerning the tectonic units on the northwestern side of the Insubric Line (Fig. 4a), we mainly referred to the previous tectonic maps of Bigi et al. (1983), Manzotti et al. (2014) and Brack et al. (2010). For each unit, we have considered a representative and simplified composition in terms of rocks and mineral percentages (Table 1). The continental Monte Rosa and Gran Paradiso units display high percentage of orthogneiss (mainly dominated by quartz, K-feldspar and biotite). In the units derived from the Piemonte-Liguria ocean, calcschists (mainly containing calcite, quartz and white mica) are dominant in the Combin Zone, whereas the Zermatt-Saas Zone largely consists of serpentinite. A high percentage of weakly retrogressed granulite is found in some units (e.g. Valpelline Series and IIDK, Table 1) of the Sesia-Dent Blanche nappes.

As a result, we obtain a surface rock-density map (Fig. 4b), composed of different units and each characterized by a representative mean density. Despite being an approximation in petrological and geological terms, this map serves well the purpose of accounting for major surface rock-density variations and deviations from 2670 kg m⁻³, to be included in the gravity anomaly computation. Some units are partially extending outside the study area, still we compile them in order to account for distant contributions on the gravity measured at each site.

To further consider the existing geological knowledge at smaller spatial scales, we have included additional information from the SAPHYR rock property catalogue (Zappone & Bruijn 2012, Zappone et al. 2015), containing point-wise bulk-density information from laboratory analyses on rock samples from 164 locations across the IVZ. This information has been included in the density map, by smoothly interpolating between the added SAPHYR points and the assigned background unit densities (explained above). In practice, each bulk-density value from SAPHYR was assigned to the nearest point in our digital elevation model (90 m resolution) and the smooth interpolation was performed following a Gaussian-shaped bell of 1 km half-width, which is comparable to the size of the smallest mapped unit considered here. As a result, we obtain a surface density map (Fig. 4b), showing the density associated with each mapped unit and containing local heterogeneities from point-wise information.

## 3.2 Bouguer gravity anomaly and Niggli correction

The point-wise complete Bouguer gravity anomaly map shows a pronounced positive feature, that extends throughout the whole study area trending northeast to southwest (Fig. 5a). The maximum anomaly amplitude varies along the crest line along strike, while the cross-strike extent above 0 mGal is ca. 30 km, decreasing gently to the east and steeply to the west.

In the calculation of the complete Bouguer anomaly, the following classical equation is used:

$$BA = G_{\text{obs}} - G_0 - FA - BP - TC, \quad (1)$$

where $G_{\text{obs}}$ is the observed gravity value, $G_0$ is the ellipsoidal gravity at the station location, FA is the free-air correction, BP is the Bouguer Plate correction, and TC is the topographic correction. In this equation, a constant and standard density value of 2670 kg m⁻³ is considered both in the Bouguer Plate and in the terrain correction terms. In order to consider the in situ rock densities as introduced and compiled above, both terms need to be corrected. The terrain correction will consider a 3-D density model between topography and station elevation, and the Bouguer Plate correction will consider the same 3-D density model between the station elevation and sea level. The amount of correction depends on the difference of local densities with respect to the standard 2670 kg m⁻³. The values of the density-dependent corrections BP* and TC* compared to the standard corrections BP and TC are shown on Figs 6(a) and (b), and can reach ca. 40 and ca. 5 mGal, respectively.

We then define the Niggli correction (NC) as a term that considers the rock density-dependence between topography and sea level:

$$NC = (BP - BP^*) + (TC - TC^*) \quad (2)$$

Adding this correction term NC to the equation of the Bouguer anomaly, we obtain the Niggli anomaly (NA)

$$NA = BA + NC = G_{\text{obs}} - G_0 - FA - BP - TC + NC \quad (3)$$

In other words, the Niggli anomaly is a type of Bouguer anomaly where local rock densities have been taken into account for the Bouguer Plate and terrain correction terms:

$$NA = G_{\text{obs}} - G_0 - FA - BP^* - TC^* \quad (4)$$

The Niggli correction is obtained by incorporating the surface density map in the same GRAVPROCESS routine that has been previously used to compute the Bouguer gravity anomaly. At the terrain correction stage, each triangle face of the surface mesh geometry, vertically extended downward until sea level, defines a volume unit. The topographic masses between surface and sea level are a compact ensemble of such volume units. The Bouguer anomaly is computed via summation of line integrals along the perimeter of every single volume unit and, thanks to the linearity of the problem, the Niggli correction can be defined by the same calculation, by multiplying each unit integral with the associated local density deviation from the reference value of 2670 kg m⁻³. As a consequence, the Niggli correction is positive in the areas with a lower surface rock density while it is negative in the areas of higher surface rock density.

While developing the 3-D crustal density model to explain the field observations, we will refer to and aim at reproducing the Niggli anomaly.

This allows us to simplify the definition of the model geometry in our approach described below, as we account for the shallow and complex structures mapped at the IVZ surface directly with the Niggli correction term. We consider the compiled surface rock-units as a first valid representation of the crustal complexity between the surface and the sea level in our study area. Hence, we eventually model the density distribution at depth only up to sea level, assuming in our model that the layer above sea level has already been accounted for based on the compiled surface rock-density map (Fig. 4b).

We note that the main features of the Niggli anomaly map are similar to that of the Bouguer anomaly map, but the local maxima reported by earlier studies (Fig. 2) are either diminished in amplitude or vanished.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Point-wise Bouguer gravity anomaly (BA), computed from the currently available merged gravity dataset. (b) Point-wise Niggli correction (NC), computed by taking into account mapped surface rock-densities in the terrain correction, as described in the text. This Niggli correction is negative in areas of denser rocks at the surface, and positive in areas of lighter rocks at the surface.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Effect of considering *in situ* rock densities on gravity corrections, plotted versus elevation. (a) Bouguer Plate correction BP, with the standard constant value of 2670 kg m⁻³ everywhere (red), and when considering the 3-D variation of densities (green). Yellow lines indicate the orientation of higher and lower density effects. (b) Terrain correction TC, with the standard approach (red), and when 3-D density variations are taken into account (green). Yellow legend indicates the trend of the effects caused by higher and lower density variations with respect to the standard value.</figcaption></figure>

## 4 3-D DENSITY MODEL CONSTRUCTION

The crustal density model has been developed using the IG-MAS + software (Schmidt *et al.* 2010). This tool allows the creation of 3-D models using triangulated polyhedrons and/or triangulated grids, with the possibility of manipulating geometries and physical parameters such as density, susceptibility and magnetization. For this work, density and model geometry have been iteratively modified to reach an optimal fit to the data.

Here we present a 3-D crustal density model, defined by one single density-contrast interface with respect to a constant-density

background, from 60 km depth up to sea level. Despite the reality being most likely closer to spatial density gradients, the sharp contrast of our model serves efficiently in fitting the pronounced gravity anomaly we observe, which would be increasingly complex and non-unique with broad density gradients. We therefore model the geometry of the interface defining the density contrast and the amplitude of the density jump itself, with respect to a given background. Therefore, the modelling is not constraining the absolute density values across the interface geometry, but only their difference. The initial model geometry is based on the regional Moho discontinuity map of Spada *et al.* (2013), locally modified to match
the 7 km s⁻¹ P-wave velocity contour surface from the local earthquake tomography model of Diehl *et al.* (2009), in the vicinity of the IGB gravity anomaly. The software package GeomIO (Bauville & Baumann 2019) was used to interpolate the initial model interface on a 3 km by 3 km grid to serve as the input in IGMAS + . In order to develop this model to fit the gravity anomalies, we have assigned a target density contrast of 400 kg m⁻³, which is consistent with what has been suggested by previous authors: 400 kg m⁻³ by Berckhemer (1968), 350 kg m⁻³ by Kissling *et al.* (1984), from 360 to 410 kg m⁻³ from Bürki (1990).

Within the IGMAS + workflow the input model geometry is in-
terpolated and displayed over an ensemble of user-defined west–east
cross-sections (10 km spaced from each other in north–south direc-
tion). We have left the far-field structure untouched and coincident
with the Moho discontinuity map, constrained by previous studies
(Spada et al. 2013). On the contrary, our geometry modifications
focused on the IGB gravity anomaly and where the IGB is sup-
posed to be more pronounced and extending towards the surface,
presenting the iconic bird's head shape as suggested by previous
investigations (e.g. Berckhemer 1968).

At each modelling step, the user is free to modify either the
geometry and/or the value of the associated model density contrast;
the software displays the updated comparison between synthetic
and observed data, allowing an immediate visual feedback.

In order to investigate the feasibility of a 400 kg m⁻³ density–
contrast model, we iteratively modified the geometry along each
cross-section within our study area. At regular intervals between
geometry modifications, we ran an automated inversion (available
within IGMAS + tools) which provides, at a given geometry, the
density–contrast value minimizing the misfit between observed and
synthetic data. We observed that a minimum of 5 rounds of geometry
change and density contrast inversions were necessary to reach a
stable solution with *ca.* 10 mGal RMS misfit, and beyond a few
more iterations the user effect became important without obvious
model updates.

For each small geometry modification, the inverted density con-
trast varied of small amounts: *ca.* ± 10 kg m⁻³ with respect to
the previous iteration. This practice served as a step-by-step guide
during the model development. In order to fit the positive gravity
anomalies, our model geometry rapidly converged towards a gently
E–SE dipping structure, similar to the ones reported in the literature
(e.g. Kissling *et al.* 1984).

5 RESULTS

5.1 Model geometry

The proposed reference model for the IGB presents a 400 kg m⁻³
density–contrast interface, which is significantly closer to the sur-
face than the initial model. The modelled structure is more than
10 km shallower in the areas where the anomaly is most pronounced,
reaching as shallow as 1 km depth below sea level, which is also
shallower than other proposals in the literature (e.g. 3 km by Berck-
hemer 1968). Given this density contrast, we require the IGB to be
not only shallower but also to have a *ca.* 1.5–2 times wider neck
than the earlier models (Berckhemer 1968). However, we have no
resolution to resolve the very thin beak *of the bird’s head* proposed
earlier. The neck width was not a well-constrained feature as only
one side of it was seismically imaged (Fig. 9d).

With this reference model, we reproduced the main features of the
gravity anomaly across the IVZ in terms of maximum amplitude,

east-west extension and along-strike variations (Fig. 7). In the next
section, we discuss the model sensitivity to the density contrast,
showing how the geometry is modulated by the density variations
with respect to the reference model.

Six east-west profiles were extracted from the model's target area for a more detailed comparison between observations and synthetics (Fig. 8), presenting a satisfactory fit to the data in the areas where the anomaly is most pronounced. The far-field effects are less well-explained than the central area, most likely because of the thick sedimentary cover effects on the Eastern side of the study area (Po Plain) and other possible but much smaller density anomalies different from the IGB itself.

Along these profiles, the model misfit is defined as the RMS of
the anomaly residuals, within an interval focused around the main
peak for each west–east profile Fig. 8. This misfit definition omits
the far-field residuals, which are not directly related to the IGB and
are beyond the scope of this work. For this reference model, we
have a RMS of ~7.5 mGal over the selected points and profiles of
interest (Fig. 8).

5.2 Model sensitivity

The 400 kg m⁻³ density contrast of our reference model provides a
satisfactory fit to the data (Fig. 9a) and suggests a shallow structure
that explain the gravity anomaly (Fig. 9b). However, very small
modifications of individual vertices of the model geometry (on
the order of a few 100 s of metres) and subsequent density–contrast
adaptations are suitable to keep the fit to the data almost unchanged.

Hence, we have investigated the model geometry sensitivity with
respect to the density–contrast values, in order to see how the model
behaves when we deviate from the reference model 400 kg m⁻³. We
present end-member models, whose geometry have been developed
starting from the reference model at two different density contrasts:
a lower (300 kg m⁻³) and a higher (500 kg m⁻³) density–contrast
model.

In both cases, we have modified the IGB structure in 3-D until a
reasonable fit the data was recovered again (Fig. 9c). While visually
the fits seem to be acceptable, both models are slightly worse in
terms of RMS misfit than our reference model. For example, for
the profile shown in Fig. 9: 5.23 mGal for the high, and 6.70 mGal
for the low density–contrast model, compared to 4.23 mGal of the
reference model.

Changing the density–contrast has consequences on the model
geometry in both the lower and higher density–contrast cases
(Fig. 9d). The former requires the IGB to have an even wider neck
and a shallower and wider head: the IGB head reaches and is lim-
ited by the virtual sea level, implying that 300 kg m⁻³ represents
the lower limit for the IGB density contrast. In this case, a broader
density gradient instead of a single interface does not seem to be a
realistic solution, as it would have difficulties explaining the sharp
gravity anomaly shape. The high density–contrast model presents a
narrower neck and a deeper top of the interface, located at *ca.* 2 km
below sea level. In this case, the higher density–contrast is providing
a sharp peaked gravity anomaly, suitable to fit the observations: we
may obtain a plausible solution by distributing this density contrast
across a few kilometres thick interface by using a spatial gradient,
instead of a sharp contrast located at the interface boundary. We
estimate this spatial gradient to be up to 4 km thick.

Nevertheless, all three tested models require a shallow-lying
anomalous structure with a somewhat sharp density contrast to
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Gravity data and model. (a) Observed Niggli gravity anomaly, interpolated over the study area. Selected gravity data points along six west-east oriented profiles serving for the comparison with the model across the most pronounced anomaly are highlighted. (b) Interpolated synthetic Niggli gravity anomaly, obtained from 3-D density modelling (see text for details). The same west-east profile locations are given as a reference.</figcaption></figure>

explain the observed gravity data. This range of possible density contrasts is directly usable for petrological modelling of the IVZ (Pistone et al. 2019). Furthermore, in order to test the user-uncertainty and the effect of non-uniqueness in gravity modelling on our work, we developed two other end-member models with higher and lower density contrasts (500 and 300 kg m⁻³, respectively), this time starting directly from the initial structure, and not from the reference model solution. In both cases, we reached the same order-of-magnitude RMS of previous models and the same main results: the IGB being thinner and further away from surface in the higher density case and thicker and closer to surface in the lower density case. During these tests, it became clear that near-surface points are better constrained from the gravity measurements, and the constraints are looser at depth. Given our data set's resolution ability (about 1 gravity point per 4–9 km² and at ca. 1 mGal uncertainty), the horizontal position of the IGB neck is resolved at <5 km on either side, and the vertical position of the IGB head-top is resolved at ca. 1 km, depending on the chosen density contrast.

# 6 DISCUSSION

The effect of the Niggli gravity anomaly correction (Fig. 5b) is not negligible when it comes to study areas with such complex surface geology as the rocks of the IVZ (Fig. 4). In our study area,

the Niggli correction flattens certain local maxima of the Bouguer gravity anomaly (in particular the ones located close to denser rock units) and it shifts the main gravity anomaly crest-line to the east, as an effect of accounting for denser and high elevation (i.e. thicker) rock units on the western side of our domain (e.g. units 1, 3 and 13).

Still, the model presents some limitations beyond the inherent non-uniqueness of gravity modelling and potential field methods in general, which could be addressed in future works, using other geophysical data such as receiver functions analyses, allowing for seismic discontinuity mapping and additional geometry constraints.

Here below we address the limitations of the proposed 3-D density model.

Regarding the amplitude of the Niggli correction, some uncertainties affect the computation itself due to the uncertainty in the single mineral-density estimates as solid solution end-members which differ substantially for pure Mg and Fe end-members. This is the case for olivine, pyroxenes and garnet (Table 1), where mineral density variations can change by few hundreds of kg m⁻³ depending on the chemical composition. However, on average for a given rock, and ranging through different mineral end-members, typical density variations are on the order of ten to a few tens (and up to a maximum of 100) kg m⁻³ of the final rock density. The density extrapolation to whole unit density may be affected by the same order-of-magnitude uncertainty, if one was to explore
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between observed and modelled Niggli gravity anomaly along the six west-east oriented profiles presented in Fig. 6. Profile locations are given on the left-hand side. For each profile, gravity points within 5 km distance on both sides have been considered. The mGal scale is valid for all six profiles, shown in an offset way. The RMS misfit has been computed around the main peak, between the two vertical black dashed lines.</figcaption></figure>

all the possible combinations of mineral abundances and rock proportions. Nevertheless, 100 kg m⁻³ is estimated as a conservative and safe upper bound for the uncertainty related to unit density association.

As an example to quantify the effect of mineral density uncertainties on the observed gravity anomaly, we computed the Niggli correction of a middle-sized, middle-elevation unit in the centre of our study area (i.e. the IIDK unit, identified by number 3 in Table 1 and Fig. 4a), with a 100 kg m⁻³ density change of the whole unit with respect to its own reference density. The associated Niggli correction variation at a gravity point nearby is ~2.7 mGal, which is compatible with the scopes and below uncertainties of our modelling work.

Another limitation stems from the model geometry itself. The 3-D crustal structure is built by working on parallel cross-sections which are 10 km apart from each other, not allowing for much smaller details or for single unit bodies to be modelled in detail

(e.g. Balmuccia peridotite unit, extending N-NE to S-SW with ~4 km length and 1 km width). Also, the single-interface model is not able to capture possible density variations along strike of our study area, as only one density contrast value is used for the IGB. However, the choice of the single-interface modelling approach captures the main signal related to the IGB. Finally, the user effect in the construction of the model involves spatial ambiguities and flexibility below the kilometre scale. The main features such as the shallow-lying anomaly, the steeply E-SE dipping body, and the arcuate shape in map view of the IGB are well constrained.

The proposed model suggests 300 kg m⁻³ as a possible lower bound for the density contrast associated with the IGB with respect to the surrounding lithologies, which would not fit the observations if distributed across a broad vertical gradient instead of a sharp interface. Hence, assuming a standard (upper crustal) background density of 2670 kg m⁻³ our model requires the IGB composing rocks to be at least 3000 kg m⁻³. Such a density may point to both mafic
<figure><img src="image_11.png" /><figcaption>**Figure 9.** West-east gravity profile and corresponding density model cross-section at *y* = 50 km across the study area. (a) Observed and modelled Niggli gravity anomaly. The RMS misfit has been computed around the main peak, between the black dashed lines. (b) Respective crustal density model cross-section using a 400 kg m⁻³ density-contrast interface, in comparison with the initial model defined by the 7 km s⁻¹ *Vp* contour from local earthquake tomography (dashed line). (c) and (d) Same profile but use a lower (300 kg m⁻³) and a higher (500 kg m⁻³) density contrast to explore the sensitivity of the model geometry. See main text for more details.</figcaption></figure>

and high-grade felsic metamorphic rocks commonly documented in middle to lower continental crust (e.g. gabbro, stronalite, kinzigite; Hacker *et al. 2015*). However, the limit of ~3000 kg m⁻³ should be regarded as a minimum density for the rocks composing the IGB. Indeed, the density model would equally fit the observed data if there was the same crustal background density increase with depth in the whole model instead of a uniform constant one. The assumption of such a gradient is closer to reality: the present-day lower continental crust is likely made of mafic rocks and/or high-grade metamorphic rocks (Bois *et al. 1989*; Fountain *1989*; Rey *1993*; Schmid & Kissling *2000*; Hacker *et al. 2015*). This implies that we can expect the background crustal density to reach ~3000 kg m⁻³ above the Moho. Keeping the same density contrast of the IGB of 300 kg m⁻³ along the model interface would imply that, in turn, the density of the IGB is 3300 kg m⁻³ at its root, which is at ~30 to 40 km depth below IGB head. Such a density for the IGB would exclude felsic rocks (Hacker *et al. 2015*) and rather point either to mafic rocks that underwent high pressure metamorphism (eclogites; Aoki & Takahashi *2004*), ultramafic plutonic rocks such as pyroxenites or garnet hornblendites (e.g. Müntener & Ulmer *2006*) or mantle peridotites.

Geological arguments can help to sharpen the lithological com-
position of the IGB. The IGB gravity anomaly is bounded to its
western side by the Insubric Line (Fig. 4), as underlined by with the
0 mGal anomaly contour line (Fig. 5b). Both geological (Schmid
et al. 1987, 1989; Berger et al. 2012) and near-surface geophysi-
cal arguments (Greenwood et al. 2018) indicate that the Insubric
Line is a crustal scale, subvertical to vertical structure: this implies
that the IGB structurally belongs to the Southern Alps. While com-
mon in the Western Alps, Alpine high-pressure metamorphic rocks

are unknown in the Southern Alps, as these are mostly dominated
by high-temperature metamorphic rocks. This makes high-pressure
metamorphic rocks (eclogites) a rather unlikely lithology for the
IGB. In consequence, plagioclase poor, ultramafic plutonic rocks or
mantle peridotite appear as the main candidate for the composition
of the IGB.

These considerations remain valid for higher values of IGB den-
sity contrast: the preferred 400 kg m⁻³ model, and also our model’s
upper limit at 500 kg m⁻³. Such scenarios would require to have
mafic IGB rocks at shallower depth, and to have mantle density at
the root of the IGB.

Additional field observations and a wider study area with more
detailed geological considerations would be necessary to investi-
gate in further detail the broader surroundings of the IVZ, such as
the Sesia Zone and the Po Plain sedimentary basin. However, the
conclusions on the dominant gravity anomaly of the IGB in relation
to the rocks exposed in the IVZ remain valid. These results are
also mutually valuable in the perspective of forthcoming projects
in the same area such as the DIVE project (Drilling the Ivrea Ver-
bano zonE; Pistone et al. 2017), which would provide additional
joint geophysical and petrological investigations across scales in
the IVZ.

Finally, the new constraints on the density, shape and total volume of dense IGB rocks in the IVZ call for investigating how such rocks could have been emplaced at shallow depth. Detailed geological observations coupled to 3-D thermomechanical modelling, with appropriately selected physical parameters (density, rheology) may shed light on the involved mechanisms, and may help to understand why other deep-crustal section are exposed but remain rare on the continents.
# 7 CONCLUSIONS

We have integrated a newly compiled gravity database with 207 new gravity points we measured across the IVZ, in order to investigate the IGB at higher resolution and in a novel way. By defining the Niggli anomaly, we have incorporated the geological surface observations from the field into the geophysical investigations in form of an adapted density correction, bridging across different spatial scales, from regional gravity anomaly trends to laboratory rock sample analyses.

The result is a 3-D single-interface crustal density model of the IGB, which presents an optimal fit to the observed data at 400 kg m⁻³ density contrast with respect to its environment. The model result supports a dense body with a structure located close to surface, with new constraints suggesting 0–2 km depth below sea level.

The modelling procedure highlights the need for a shallow and pronounced density contrast or spatially sharp (<4 km) density gradient in order to fit the observations. The range of plausible density contrasts is between 300 and 500 kg m⁻³, with the preferred model at 400 kg m⁻³.

Considering a natural increase of densities with depth, and the metamorphic grade of rocks in the area, felsic and high-pressure metamorphic rocks (eclogites) were excluded as main composing rocks if the IGB. Instead, mafic or ultramafic plutonic rocks are most likely for the lithology of the IGB.

Further efforts will be spent to investigate the structure beneath the IVZ, especially by jointly analysing gravity and seismological data we have collected in the same region. Starting in 2017, 10 broad-band seismic stations have been installed along a west-east profile across the IVZ and have continuously recorded seismic data for 2 yr (https://doi.org/10.5281/zenodo.1038209). The seismic data will be processed by the receiver function method in order to investigate seismic discontinuity locations and sharpness, and the seismic crustal structure. The joint use of gravity and seismic data will put new constraints on the IGB geometry and on its physical properties. This, in turn, will provide useful input for geological and numerical models exploring the emplacement of deep-crustal bodies at shallow depth, which may explain the mechanism and relative rareness of Ivrea-type outcrops across the globe.

# ACKNOWLEDGEMENTS

We would like to acknowledge many colleagues for the precious contributions and various inputs that made this work possible: U. Marti, R. Barzaghi, the OGS institution, A. Zappone, T. Diehl, R. Cattin and S. Mazzotti for the data sharing, the internal discussions and the GRAVPROCESS software sharing and editing, H.J. Götze and C. Spooner for IGMAS + introduction; C. Alvizuri, G. Moradi and S. Subedi for the fieldwork efforts together with all the local people in the IVZ and its surrounding for the precious help during the measurement campaigns. We would like to acknowledge the Swiss National Science Foundation (SNF) which is supporting this research (grant numbers PP00P2_157627 and PP00P2_187199). MP acknowledges the support of the SNF-sponsored Fellowship (PZ00P2_168166 grant). Author contribution statement: MS and GH designed the study, MS, GH, TB and LB collected the new field data and performed the processing, PM provided the database for density estimations, MS performed the modelling, MS, GH, PM, BP, MP and OM contributed to the interpretation, MS and GH drafted and wrote the manuscript, all authors revised the manuscript and approved the final version. The final IGB model geometry is available at the following doi: 10.5281/zenodo.3897507

# REFERENCES

Ansorge, J., Mueller, S., Kissling, E., Guerra, I., Morelli, C. & Scarascia, S., 1979. Crustal section across the zone of Ivrea-Verbano from the Valais to the Lago Maggiore, *Boll. Geofis. Teor. Appl.*, **21**(83), 149–157.

Aoki, I. & Takahashi, E., 2004. Density of MORB eclogite in the upper mantle, *Phys. Earth Planet. Inter.*, **143–144**, 129–143.

Bauville, A. & Baumann, T.S., 2019. geomIO: an open-source MATLAB toolbox to create the initial configuration of 2-D/3-D thermo-mechanical simulations from 2-D vector drawings, *Geochem. Geophys. Geosyst.*, **20**(3), 1665–1675.

Bayer, R., Carozzo, M.T., Lanza, R., Miletto, M. & Rey, D., 1989. Gravity modelling along the ECORS-CROP vertical seismic reflection profile through the Western Alps, *Tectonophysics*, **162**(3-4), 203–218.

Berckhemer, H., 1968. Topographie des "Ivrea-Körpers" abgeleitet aus seismischen und gravimetrischen Daten, *Schweiz. Mineral. Petrogr. Mitt.*, **48**(1), 235–246.

Berger, A., Mercolli, I., Kapferer, N. & Fügenschuh, B., 2012. Single and double exhumation of fault blocks in the internal Sesia-Lanzo Zone and the Ivrea-Verbano Zone (Biella, Italy), *Int. J. Earth Sci.*, **101**, 1877–1894.

Bigi, G., Cosentino, D., Coli, M., Parotto, M., Sartori, R. & Scandone, P., 1983. Structural model of Italy, sheet no. 1, 1: 500,000. CNR, Italy.

Bois, C., Pinet, B. & Roure, F., 1989. Dating lower crustal features in France and adjacent areas from deep seismic profiles, in *Properties and Processes of Earth's Lower Crust*, AGU, pp. 17–31. doi:10.1029/GM051p0017.

Boriani, A., 1971. Geology of the Junction Between the Ivrea-Verbano and Strona-Ceneri Zones (Southern Alps), in *Memorie Dell'Istituto Di Geologia e Mineralogia, Università Di Padova*, **28**, 1–36.

Bousquet, R. et al., 2012. Metamorphic framework of the Alps-Carte metamorphique des Alpes CCGM/CGMW. CCGM/CGMW.

Brack, P., Ulmer, P. & Schmid, S.M., 2010. A crustal-scale magmatic system from the Earth's mantle to the Permian surface: field trip to the area of lower Valsesia and Val d'Ossola (Massiccio dei Laghi, Southern Alps, Northern Italy), *Swiss Bull. für angewandte Geol.*, **15**(2), 3–21.

Bürki, B., 1990. Geophysical interpretation of astrogravimetric data in the Ivrea Zone, in *Exposed Cross-Sections of the Continental Crust, Series C: Mathematical and Physical Sciences*, Vol. **317**, pp. 545–561, eds Salisbury, M.H. & Fountain, D.M., Springer.

Cattin, R., Mazzotti, S. & Baratin, L.M., 2015. GravProcess: an easy-to-use MATLAB software to process campaign gravity data and evaluate the associated uncertainties, *Comput. Geosci.*, **81**, 20–27.

Corchete, V., 2010. The high-resolution gravimetric geoid of Italy: ITG2009, *J. Afr. Earth Sci.*, **58**(3), 580–584.

Diehl, T., Husen, S., Kissling, E. & Deichmann, N., 2009. High-resolution 3-DP-wave model of the Alpine crust, *Geophys. J. Int.*, **179**(2), 1133–1147.

Farr, T.G. et al., 2007. The shuttle radar topography mission, *Rev. Geophys.*, **45**(2), doi:10.1029/2005RG000183.

Fountain, D.M., 1976. The Ivrea—Verbano and Strona-Ceneri Zones, Northern Italy: a cross-section of the continental crust—new evidence from seismic velocities of rock samples, *Tectonophysics*, **33**(1-2), 145–165.

Fountain, D.M., 1989. Growth and modification of lower continental crust in extended terrains: the role of extension and magmatic underplating, in *Properties and Processes of Earth's Lower Crust, Geophysical Monograph Series*, Vol. **51**, pp. 287–299, eds Mereu, R.F, Mueller, S. & Fountain, D.M., AGU, doi:10.1029/GM051p0287.

Fountain, D.M., Percival, J. & Salisbury, M.H., 1990, *Exposed Cross-Sections of the Continental Crust*, Springer.

Greenwood, A. et al., 2018. High-resolution seismic reflection survey across the Insubric Line, Italian Alps, in *Proceedings of the 8th International Conference on Environmental and Engineering Geophysics*, 10–14 June 2018, Hangzhou, China.

Hacker, B., Kelemen, P.B. & Behn, M.D., 2015. Continental Lower Crust, *Annual Review of Earth and Planetary Sciences*, **43**:167–205.

Handy, M.R., Schmid, S.M., Bousquet, R., Kissling, E. & Bernoulli, D., 2010. Reconciling plate-tectonic reconstructions of Alpine Tethys with the geological-geophysical record of spreading and subduction in the Alps, *Earth Sci. Rev.*, **102**, 121–168.
Handy, M.R., Ustaszewski, K. & Kissling, E., 2015. Reconstructing the Alps–Carpathians–Dinarides as a key to understanding switches in subduction polarity, slab gaps and surface motion, *Int. J. Earth Sci.*, **104**(1), 1–26.

Hinze, W.J., 2003. Bouguer reduction density, why 2.67?, *Geophysics*, **68**(5), 1559–1560.

Karakas, O. et al., 2019. The pace of crustal-scale magma accretion and differentiation beneath silicic caldera volcanoes, *Geology*, **47**(8), 719–723.

Khazanehdari, J., Rutter, E.H. & Brodie, K.H., 2000. High-pressure-high-temperature seismic velocity structure of the miderrustal and lower crustal rocks of the Ivrea-Verbano zone and Serie dei Laghi, NW Italy, *J. geophys. Res.*, **105**(B6), 13 843–13 858.

Kissling, E., Wagner, J.J. & Mueller, S., 1984. Three-dimensional gravity model of the northern Ivrea-Verbano Zone, Geomagnetic and Gravitometric Studies of the Ivrea Zone, *Matér. Géol. Suisse. Géophys.*, **21**, 55–61.

Manzotti, P., Zucali, M., Balle`vre, M., Robyr, M. & Engi, M., 2014. The tectonometamorphic evolution of the Sesia—Dent Blanche Nappes (internal Western Alps): review and synthesis, *Swiss J. Geosci.*, **107**, 309–336.

Marson, I., Orlando, L. & Stoka, M., 1994. Gravity model on the CROP profile, *Proc. Symp. CROP Alpi centrali. Quaderni di Geodinamica Alpina e Quaternaria*, **2**, 161–169.

Müntener, O. & Ulmer, P., 2006. Experimentally derived high-pressure cumulates from hydrous arc magmas and consequences for the seismic velocity structure of lower arc crust, *Geophys. Res. Lett.*, **33**, L21308, doi:10.1029/2006GL027629.

Niggli, E., 1947. Über den Zusammenhang zwischen der positiven Schwereanomalie am Südfuß der Westalpen und der Gesteinszone von Ivrea.

Petri, B., Duretz, T., Mohn, G., Schmalholz, S.M., Karner, G.D. & Müntener, O., 2019. Thinning mechanisms of heterogeneous continental lithosphere, *Earth planet. Sci. Lett.*, **512**, 147–162.

Pistone, M., Müntener, O., Ziberna, L., Hetényi, G. & Zanetti, A., 2017. Report on the ICDP workshop DIVE (Drilling the Ivrea–Verbano zonE), *Scientific Drilling*, **23**, 47–56.

Pistone, M., Ziberna, L., Hetényi, G., Scarponi, M., Zanetti, A. & Müntener, O., 2019. Deciphering the petrophysics and petrology of the Ivrea Geophysical Body (Alps) by combining laboratory data and thermodynamic modelling, in *Proceedings of the AGU Fall Meeting, S23G-0712*.

Quick, J.E. & Denlinger, R.P., 1993. Ductile deformation and the origin of layered gabbro in ophiolites, *J. geophys. Res.*, **98**(B8), 14 015–14 027.

Quick, J.E., Sinigoi, S. & Mayer, A., 1994. Emplacement dynamics of a large mafic intrusion in the lower crust, Ivrea-Verbano Zone, northern Italy, *J. geophys. Res.*, **99**(B11), 21 559–21 573.

Quick, J.E., Sinigoi, S., Snoke, A.W., Kalakay, T.J., Mayer, A. & Peressini, G., 2003. Geologic map of the southern Ivrea-Verbano Zone, northwestern Italy, US Geological Survey.

Rey, D., 1990. Gravity and aeromagnetic maps of the western Alps: contribution to the knowledge of the deep structures along the ECORS-CROP seismic profile, *Mém. Soc. Géol. France*, **156**, 107–121.

Rey, P., 1993. Seismic and tectono-metamorphic characters of the lower continental crust in Phanerozoic areas: a consequence of post-thickening extension, *Tectonics*, **12**, 580–590.

Rivalenti, G., Mazzucchelli, M., Vannucci, R., Hofmann, A.W., Ottolini, L., Bottazzi, P. & Obermiller, W., 1995. The relationship between websterite and peridotite in the Balmuccia peridotite massif (NW Italy) as revealed by trace element variations in clinopyroxene, *Contrib. Mineral. Petrol.*, **121**(3), 275–288.

Schmid, S.M., Aebli, H.R., Heller, F. & Zingg, A., 1989. The role of the Periadiatic Line in the tectonic evolution of the Alps, *Geol. Soc. Lond., Spec. Publ.*, **45**, 153–171.

Schmid, S.M., Fügenschuh, B., Kissling, E. & Schuster, R., 2004. Tectonic map and overall architecture of the Alpine orogen, *Eclogae Geol. Helv.*, **97**(1), 93–117.

Schmid, S.M. & Kissling, E., 2000. The arc of the western Alps in the light of geophysical data on deep crustal structure, *Tectonics*, **19**, 62–85.

Schmid, S.M., Kissling, E., Diehl, T., van Hinsbergen, D.J. & Molli, G., 2017. Ivrea mantle wedge, arc of the Western Alps, and kinematic evolution of the Alps–Apennines orogenic system, *Swiss J. Geosci.*, **110**(2), 581–612.

Schmid, S.M., Zingg, A. & Handy, M., 1987. The kinematics of movements along the Insubric Line and the emplacement of the Ivrea Zone, *Tectonophysics*, **135**, 47–66.

Schmidt, S., Götze, H.J., Fichler, C. & Alvers, M., 2010. IGMAS+–a new 3D gravity, FTG and magnetic modeling software, in *GEO-INFORMATIK Die Welt im Netz*, pp. 57–63.

Sinigoi, S., Quick, J.E., Demarchi, G. & Klötzli, U., 2011. The role of crustal fertility in the generation of large silicic magmatic systems triggered by intrusion of mantle magma in the deep crust, *Contrib. Mineral. Petrol.*, **162**(4), 691–707.

Solarino, S. et al., 2018. Mantle wedge exhumation beneath the Dora-Maira (U) HP dome unravelled by local earthquake tomography (Western Alps), *Lithos*, **296**, 623–636.

Spada, M., Bianchi, I., Kissling, E., Agostinetti, N.P. & Wiemer, S., 2013. Combining controlled-source seismology and receiver function information to derive 3-D Moho topography for Italy, *Geophys. J. Int.*, **194**, 1050–1068.

Thouvenot, F., Paul, A., Senechal, G., Hirn, A. & Nicolich, R., 1990. ECORS-CROP wide-angle reflection seismics: constraints on deep interfaces beneath the Alps, *Mémoires de la Société géologique de France*, **156**, 97–106.

Vecchia, O., 1968. La zone Cuneo-Ivrea-Locarno, élément fondamental des Alpes. Géophysique et géologie, *Schweizerische Mineralogische Petrographische Mitteilungen*, **48**, 215–226.

Zappone A., & Zappone, R.H.C., 2012. The Swiss Atlas of physical properties of rocks (SAPHYR), *Swiss Geophysical Commission Report*, http://sgpk.ethz.ch/static/jahresbericht/2012/zappone2012.pdf.

Zappone, A., Wenning, Q. & Kissling, E. 2015 SAPHYR: The Swiss Atlas of physical properties of rocks, *Swiss Geophysical Commission Report*, http://sgpk.ethz.ch/static/jahresbericht/2015/SAPHYRreport2015.pdf.

Zingg, A., 1983. The Ivrea and Strona-Ceneri zones (Southern Alps, Ticino and N-Italy)—a review, *Schweiz. Mineral. Petrogr. Mitt.*, **63**(2-3), 361–392.