# Sequential inversion of local earthquake traveltimes and gravity anomaly-the example of the western Alps

Philippe Vernant, Frédéric Masson, Roger Bayer, Anne Paul

► To cite this version:

Philippe Vernant, Frédéric Masson, Roger Bayer, Anne Paul. Sequential inversion of local earthquake traveltimes and gravity anomaly-the example of the western Alps. Geophysical Journal International, 2002, 150, pp.79-90. 10.1046/j.1365-246X.2002.01694.x. insu-03607072

HAL Id: insu-03607072

https://insu.hal.science/insu-03607072

Submitted on 13 Mar 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
# Sequential inversion of local earthquake traveltimes and gravity anomaly—the example of the western Alps

Philippe Vernant,<fnref n="1" /> Frédéric Masson,<fnref n="1" /> Roger Bayer<fnref n="1" /> and Anne Paul<fnref n="2" />

Accepted 2002 January 24. Received 2002 January 21; in original form 2001 May 31

## SUMMARY

We present a joint analysis of gravity anomaly and seismic arrival time data recorded in the western Alps. Seismological data were collected by a network of 126 permanent and temporary stations implemented in 1996. A set of ~550 local events has been recorded. Gravity data result from the addition of two new gravity surveys to an existing data base. A published velocity model obtained by local earthquake tomography (LET), was used to construct an initial 3-D gravity model, using a linear velocity–density relationship (Birch’s law). While the synthetic Bouguer anomaly field calculated for this model has the same shape and wavelength as the observed anomaly, its amplitude is strongly underestimated. To derive a crustal velocity–density model that accounts for both types of observations, we performed a sequential inversion of seismological and gravity data. The variance reduction of the arrival time data for the final sequential model was comparable to the variance reduction obtained by simple LET. Moreover, the sequential model explained ~90 per cent of the observed gravity anomaly. The main features of our model compared with the LET model are: (1) an important broadening of the high-velocity anomaly associated with the high-velocity high-density Ivrea Body, (2) a 10 km thick low-velocity zone beneath the nappes of Digne and Castellane and (3) a high-velocity zone at more than 25 km depth under the internal zone of the range.

**Key words:** Alps, gravity anomalies, inverse problem, *p* waves, tomography.

## 1 INTRODUCTION

In the western Alps, gravity interpretation has often been a fruitful complementary tool for seismic methods in studies of the deep structure of this orogenic belt. Since the early 1960s, detailed gravity profiles across this belt have been used to define 2-D models that have become more and more detailed as the geometrical resolution improves because of new reflection–refraction seismic surveys (Morelli 1963; Berckhemer 1968; Ménard & Thouvenot 1984; Rey 1989). The structural interpretations have focused on defining the geometry of the Moho (Waldhauser et al. 1998) and the characteristics of the Ivrea Body. This structure, located at the eastern border of the belt, corresponds to the famous Ivrea gravity anomaly (Niggli 1946; Coron 1963). Refraction and wide-angle reflection profiles recorded between 1958 and 1966 showed a high-velocity body interpreted as a wedge of upper mantle by Closs & Labrouste (1963) and Berckhemer (1968). Reinterpretation of seismic data and 2-D gravity modelling allowed Ménard & Thouvenot (1984) to propose a structural model based on a flaking of the European lithosphere. They divided the Ivrea Body into three units: the surface unit associated with the basic and ultrabasic rocks of the Sesia Zone and Lanzo massif (Fig. 1), the main unit at ~10 km depth connected to

the Frontal Penninic Thrust and the lower unit at 30 km depth. The results of the ECORS-CROP profile (ECORS-CROP DSS Group 1989b) showed a reflector at 30 km depth below the Brinçonnais zone in agreement with the hypothesis of Ménard & Thouvenot (1984). It was interpreted as being the top of a mantle wedge by the ECORS-CROP DSS Group (1989b) and Nicolas et al. (1990). More recently Roure et al. (1996) and Schmidt & Kissling (2000) have suggested that the reflector at 30 km could be the top of the lower crust instead of mantle material.

Joint interpretation of gravity and refraction/wide-angle reflection seismic data has been performed by a *posteriori* verifying the fit of the 2-D seismic models with gravity data using a classic density–velocity relationship (Birch 1961). In most cases, the agreement between seismic and gravity data was obtained without large modifications of the velocity models, except beneath the eastern part of the internal zone of the belt where velocity discontinuities are poorly defined (see, for example, the ECORS-CROP seismic cross-section in the western Alps ECORS-CROP 1989a,b). The 2-D models have mainly proved the continuity of the main crustal structures along the arcuate strike of the belt, but they failed to provide detailed knowledge of the 3-D geometry of the alpine crust (Rey 1989; Rey et al. 1990).

<footnote n="1">¹CNRS/LGTS, ISTEEM Université Montpellier II, 4 place E, Bataillon, 34095 Montpellier Cedex 05, France. E-mail: vernant@dstu.univ-montp2.fr</footnote>

<footnote n="2">²LGIT, CNRS and Université Joseph Fourier, BP53, 38041 Grenoble Cedex 9, France</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Tectonic sketch map of the western Alps. (a) Valais, Piemont and Briançon basement, (b) Austro-Alpine basement, (c) Piemont zone (ophiolites in black), (d) Subbriançon zone, (e) Briançon zone, (f) Flyschs, (g) External crystalline massifs. Numbers indicate (1a) Dora Maira, (1b) Lanza, (1c) Sesia zone, (2a) Pelvoux, (2b) Mercantour-Argentera, (2c) Simplon-Tessin, (2d) Mont-Blanc, (2e) Aar, (2f) Belledonne-Beaufortain, (3) Sedimentary nappes of Digne and Castellane. The main geologic structures of the region are shown by the small sketch map. These features will be reported on all the maps of this study. The dark grey line is the border. FPT, Frontal Penninic Thrust; FBT, Frontal Briançonnais Thrust. Seismic stations are plotted as black triangles.</figcaption></figure>

During the previous decade, new geophysical data have been collected in order to obtain refined 3-D images of the Alpine lithosphere. Within the framework of the GéoFrance 3-D Alps project (1997), a passive seismological experiment was carried out by Paul *et al.* (2001) in the south-western Alps to record local earthquakes and compute a local earthquake tomographic (LET) velocity model. Simultaneously, new gravity data have been collected and added to previous surveys resulting in a new detailed and precise Bouguer anomaly map of the western Alps (Masson *et al.* 1999).

These two data sets provide an opportunity for joint analysis of gravity anomalies and seismic traveltimes in a region where anomalous high-density and *P*-wave velocity bodies have already been observed, albeit with poor resolution. In this paper, we show that the 3-D LET model obtained by Paul *et al.* (2001) does not explain the Bouguer anomaly, caused by the lack of resolution in some parts of the seismic model caused by a heterogeneous distribution of seismic rays. Taking into account the integrating property of the gravity field and the complementary resolving power contained in gravity
and seismic tomography, we derived a co-operative inversion of the two data sets (Lines *et al.* 1988). We first present the sequential inversion process and apply it to the western Alps. Our preferred density–velocity model and its structural implications are discussed in the final section of the paper.

## 2 FAILURE OF THE DIRECT MODELLING OF GRAVITY DATA FROM THE INITIAL LET MODEL

We first attempt to model the gravity anomaly observed in the studied area from the tomographic P-wave velocity model. The classic direct approach (forward modelling) is used to verify the goodness of the $V_p$ model by conversion of velocity into density. The associated theoretical gravity field is computed and compared with the observed one. In that aim, we use a new gravity map and seismological data resulting from surveys carried out in the south-western Alps within the framework of the GéoFrance 3-D Alps project (Masson *et al.* 1999; Paul *et al.* 2001). From west to east, the studied area (Fig. 1) overlaps the external zone including the Digne and Castellane nappes, the internal flysch nappes thrusting on the external domain, the external crystalline Pelvoux and Argentera massifs, the Briançonnais and Piemont zones (Blue Schists Piemont zone, eclogitic Viso domain and High Pressure Dora Maira massif) separated from the Apulian domain by the Insubric line. This internal zone is characterized by nappe piles of lower crustal and upper mantle slices (Nicolas *et al.* 1990; Schmidt & Kissling 2000) leading to a very complex Moho topography overlain (Waldhauser *et al.* 1998) by the dense and high-velocity Ivrea Body. It is therefore an important zone to test the consistency of the density and velocity models.

### 2.1 The seismological data and the local earthquake tomography

From 1996 August to December, 67 temporary stations were installed in the southernmost part of the western Alps. They complemented a network of 59 permanent stations managed in this region by the universities of Grenoble, Nice and Genova (Fig. 1). The average interstation distance was 10–15 km. From permanent and temporary stations, 347 local earthquakes were selected. Finally, 104 quarry blasts recorded by Italian network and 99 complementary earthquakes located a great depth and recorded by the permanent network of the University of Genova were added to the previous data set. Paul *et al.* (2001) inverted arrival times simultaneously for velocity ($V_p$ and $V_p/V_s$) and hypocentre parameters using the classic program SIMULPS (Thurber 1983; Eberhart-Philips 1993; Evans *et al.* 1994). Depth slices of the resulting 3-D $V_p$ model are shown in Fig. 2 with locations of relocated earthquake foci. The thick black line delineates the area with correct resolution according to criteria discussed in detail by Paul *et al.* (2001). Note that the surface of the well-resolved area strongly diminishes with depth and shifts to the east caused by the lack of hypocentres at depths larger than 10 km under the western half of the studied region. This model displays strong velocity contrasts. Its main features are the following (Paul *et al.* 2001).

(1) A low-velocity anomaly beneath the Digne and Castellane nappes in the external domain. This anomaly is visible from the surface up to 5 km depth.

(2) High velocities at shallow depths (0–5 km) beneath the Piemont zone, south-west of the Dora Maira Massif.

(3) A north–south high-velocity anomaly (7.4–7.5 km s⁻¹) under the Dora Maira massif and the westernmost Po plain at depths greater than 8 km corresponding to the Ivrea Body. This high-velocity structure, well defined at 12 and 30 km depth, is not well imaged between 16 and 20 km depth. Synthetic tests indicate that this result could be an artefact of the inversion owing to the rather inhomogeneous distribution of rays (Paul *et al.* 2001). The Ivrea Body is certainly a continuous unit of high-velocity material from 12 km depth down to the upper mantle.

### 2.2 The new Bouguer anomaly map of the western Alps

Approximately 1600 new data points were surveyed between 1997 and 1999 in the western French Alps to increase the existing data density in the high mountainous areas (Masson *et al.* 1999). These data and older surveys were merged and tied to the IGSN71 system. A new high-resolution Bouguer map was then constructed with topographic corrections up to 167 km distance and a density reduction of 2600 kg m⁻³ (Masson *et al.* 1999). A maximum error of 2 mgal is estimated for the zones of highest elevation. The regional trend of the map is a decrease of –150 mgal from the external domain to the internal zone caused by crustal thickening (Fig. 3a). Taking into account recent 3-D interface modelling of the Alpine crust–mantle boundary (Waldhauser *et al.* 1998), we estimated the gravity contribution of Moho depth variations. Its effect has been removed from the Bouguer anomaly map to compute a residual map shown in Fig. 3(b) and reproduced from Masson *et al.* (1999). Various density contrasts for the lower-crust–upper-mantle boundary were tested to minimize the long-wavelength anomalies of the residual map. The residual anomalies are mainly related to crustal heterogeneities. However, as reliable information on the Moho remains scarce in some portions of the studied area, the possibility of some remaining Moho effects cannot be rejected. Nevertheless, these large-wavelength uncorrected Moho effects do not blur the gravity anomalies caused by the local crustal structures. The main feature of the residual map of Fig. 3(c) is the north–south elongated positive anomaly (100 mgal) of the Ivrea Body that ends at the latitude of Cuneo. This map also shows a broad negative anomaly (–40 mgal) beneath the external nappes of Digne and Castellane. A comparison with the 3-D $V_p$ model of Fig. 2 shows that the anomalies of the external domain and the Ivrea anomaly are correlated, respectively, to low-velocity regions in the shallowest layers (0–5 km) and very high-velocity anomalies in the deepest layers (10–30 km).

### 2.3 Limitations of the forward approach

The empirical Birch’s law (1961) is commonly used to define the linear relationship between density and P-wave velocity, depending on the mean atomic weight. As the slope ($\Delta\rho/\Delta V_p$) of Birch’s law is constant independently of the mean atomic weight used, the density contrasts deduced from the LET model do not depend on the mean atomic weight, whereas absolute densities do. Therefore, we will only show absolute velocity models in this study.

In the forward approach, the $V_p$ model of Fig. 2 is converted into an *a priori* density contrast model to calculate the gravity effect of crustal heterogeneities at nodes of the regular grid (2 × 2 km²) of the residual anomaly. We adopt the space discretization of the LET model of Paul *et al.* (2001). The N–S and E–W horizontal extent of the grid is 160 × 160 km². The horizontal node spacing is 10 km, except in the region of the Dora Maira massif where the spacing reduces to 5 km. The depths of the node layers are –2.5, 0,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Smoothed P-wave velocity maps obtained by local earthquake tomography (Paul *et al.* 2001). The black lines show the limit of the resolved part of the model. The black dots are the nodes of the seismic model. The main features of this model are a north-south high-velocity anomaly under the Dora Maira massif and the westernmost Po plain at depths greater than 10 km corresponding to the Ivrea Body and a low-velocity anomaly beneath the Digne and Castellane nappes visible from the surface up to 5 km depth.</figcaption></figure>

5, 10, 15, 20, 25, 30 and 50 km. In order to calculate the theoretical gravity field, the crust is partitioned into prismatic elementary cells with uniform density contrast. From 0 to 40 km depth, the density grid is divided in eight layers, each one containing 420 cells. The velocity nodes are located at the centre of the top boundary of the density cells.

$p_0$ is the *a priori* density contrast model determined from the LET model. We define **G** as the matrix of the forward problem, where element $G_{ij}$ is the residual gravity field at point i induced by cell j with unit density. The theoretical field vector $d_0$ is calculated by solving the matrix form (Van de Meulebrouck *et al.* 1984; Richard *et al.* 1984),
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Bouguer anomaly map of the western Alps and the neighbouring regions (Masson *et al.* 1999). (b) Residual Bouguer anomaly map obtained subtracting the Moho effect. This anomaly is mainly caused by crustal heterogeneities (Masson *et al.* 1999). (c) Enlargement of the residual Bouguer anomaly map of the region under study. (d) Synthetic Bouguer anomaly map computed from the LET model of Fig. 2 converted to density using Birch's law.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Diagram showing the velocity–density relationships used to compute the gravity anomaly produced by the LET model. (1) Birch's law (Birch 1961), (2) Nafe & Drake (1957; Ludwing *et al.* 1970), (3) a Barton's simplification of the Nafe–Drake relationship (Barton 1986), (4) the minimum bound of the Nafe–Drake curve, (5) the maximum bound of the Nafe–Drake curve, (6) Glaznev *et al.* (1996). (b) Comparison along an E–W profile (latitude: Gap) between the observed residual Bouguer anomaly (Fig. 3) and the computed one using the LET model and the velocity–density relationships.</figcaption></figure>

$$Gp_0 = d_0. \qquad (1)$$

While the resulting anomaly (Fig. 3d) has the same shape and wavelength as the observed anomaly (Fig. 3c) it shows a weaker total amplitude anomaly. Three possible explanations may be given for such a discrepancy.

(i) The Birch's law is not appropriate for the crustal rocks of the Alpine belt. The definition of a velocity–density relationship for a strongly heterogeneous model that covers a large area is not a simple problem, and Birch's law may be an oversimplification. We use a strategy inspired by Barton (1986) in order to check the validity of Birch's law for the region studied. We compute the gravity anomaly deduced from the LET model using (1) Birch's law (1961), (2) the Nafe–Drake relationship (Nafe & Drake 1970; Ludwig *et al.* 1970), (3) Barton's simplification of the Nafe–Drake relationship (Barton 1986) which associates the typical velocities of the crystalline crust (5.7–7.0 km s⁻¹) with a constant density of 2.8 g cm⁻³, (4) the maximum and the minimum bounds of the Nafe–Drake curve and (5) a relationship defined by Glaznev *et al.* (1996). Fig. 4(a) presents the five Vp/ρ relationships. Whatever the relationship used, it is not possible to explain the observed gravity anomaly by a simple velocity–density conversion of the LET model (Fig. 4b). As we compute the density contrasts, we only use the Δρ/ΔVp slopes of the relationship and as they are roughly similar, the gravity profiles obtained are similar too (Fig. 4a). Therefore, the choice of the law is not crucial. The Nafe–Drake relationship was originally established for marine sediments (Nafe & Drake 1957) so we decided to use Birch's law (1961).

(ii) The velocities of the LET model are underestimated. The maximum velocity obtained at the top of the Ivrea Body by LET (~7.4 km s⁻¹) is consistent with the velocities previously obtained in the 1960s and 1970s by several refraction and wide-angle seismic experiments (Closs & Labrouste 1963; Berckhemer 1968; Choudhury *et al.* 1971; Perrier 1973; Giese & Prodehl 1976; Ansorge *et al.* 1979; Thouvenot & Perrier 1981). Producing agree-

ment between observed and synthetic residual anomalies requires an increase in the velocity anomaly of the Ivrea Body by a factor of 2. This hypothesis is also unrealistic.

(iii) The LET model minimizes the lateral extent of structures and consequently reduces the amplitude and extent of the synthetic Bouguer anomaly. This observation is supported by synthetic tests, which show that, because of the uneven distribution of hypocentres, some parts of the LET model are totally unresolved, particularly at depths greater than 15 km (Paul *et al.* 2001). Moreover, synthetic LET tests also document that the Vp anomaly associated with the Ivrea Body at the 16 and 20 km depth node layers is strongly underestimated. This hypothesis is certainly the most suitable for explaining most of the discrepancy between the observed and computed residual anomalies.

To check this hypothesis, we decided to use a cooperative inversion (Lines *et al.* 1988) of the two data sets in order to compute a 3-D velocity–density model consistent with both gravity and seismic data.

# 3 COOPERATIVE INVERSION: THE SEQUENTIAL METHOD

## 3.1 The method

Lines *et al.* (1988) propose two main kinds of cooperative inversion of geophysical data. The first one is the joint inversion where all data are inverted simultaneously. This strategy is often used for data sets combining gravity and seismics (Oppenheimer & Herkenhoff 1981; Lees & VanDecar 1991; Kaufmann & Long 1996; Zeyen & Achauer 1997). The respective weight of the two data sets is the main difficulty of the method (Lines *et al.* 1988). The second one is the sequential inversion where each data set is inverted successively. The *a posteriori* information resulting from the previous inversion of the first data set is transformed into *a priori* information to invert the second data set.
<figure><img src="image_6.png" /><figcaption>Figure 5. Flow-chart of the sequential inversion procedure.</figcaption></figure>

To avoid the weighting problem of the joint inversion and to facilitate reuse of the seismic method (Thurber 1983) previously applied by Paul *et al.* (2001) for the LET, we adopt a sequential strategy.

Our approach consists of the reiteration of a set of *n* iterations of the seismological inversion, leading to a new velocity model and new event locations, followed by one inversion of the gravity data and the computation of a new density contrast model. This procedure is repeated until the convergent criterion is satisfied. Using this method, it is possible to estimate qualitatively what information is brought to the model by seismic and gravity data.

Fig. 5 illustrates the organization of the sequential method. As already noted for the LET inversion (see Section 2.1), the starting velocity model used for the first inversion of traveltime data is the initial 1-D velocity model estimated by Paul *et al.* (2001) from 1-D inversion of their data set (Kissling *et al.* 1994). Then, the iterative SIMULPS program calculates the 3-D $V_p$ model and new earthquake locations from the arrival times of local earthquakes. This process is stopped after *n* iterations. This 3-D absolute $V_p$ model is then converted into a 3-D relative density contrast model using Birch's law. At this stage, the linear inverse gravity problem is solved leading to a new density contrast model that is transformed back to a new $V_p$ model. This completes the first loop of the sequential inversion. The following loops use the same procedure defining as input the final velocity model of the previous loop. The sequential process is stopped when the standard deviations between observed data and theoretical values calculated from the models stop decreasing significantly between two loops.

The density contrast model is computed using the linear stochastic method (Van de Meulebrouck *et al.* 1984). After dividing the crustal model into prisms (see Section 2.3), the linear equation of the forward problem has the matrix form of eq. (1). We assume that the measurement errors in the residual gravity field are independent and define a diagonal covariance matrix on the experimental gravity errors. We also assume an uncertainty on the *a priori* density contrast solution by a covariance matrix, the terms of which are given by

$$C_{ij} = \Lambda_{ij} \exp\left(-\frac{d_{ij}}{\lambda}\right), \qquad (2)$$

where

$$\Lambda_{ij} = c\sigma_i\sigma_j. \qquad (3)$$

From a stochastic point of view, this *a priori* covariance matrix ($C_{ij}$) could be deduced from the *a posteriori* covariance matrix of the velocity model. If the *a posteriori* uncertainty on $V_p$ is high, $\sigma_i$ is also high, so the density value in cell *i* can vary strongly during the gravity inversion. Conversely, if the *a posteriori* uncertainty on $V_p$ is low, indicating that the $V_p$ value is well constrained in the cell *i*, $\sigma_i$ is low, which prohibits strong variations of the density contrast in the cell. The data sets are complementary because gravity data provide information where seismological data are sparse.

$\Lambda_{ij}$ is the product of the standard deviations $\sigma_i$ and $\sigma_j$ of the density contrast in prisms *i* and *j* deduced from the *a posteriori* standard deviation of the $V_p$ model and *c* is a damping factor. $d_{ij}$ is the distance between cells *i* and *j*, and $\lambda$ is the correlation radius of the model. In this paper, the parameters *c* and $\lambda$ are chosen in order to significantly decrease the standard deviations between observed gravity and traveltime data and theoretical values computed from density and $V_p$ models in the iterative process. The factor *c* is assumed to range from 0 to 1. When it is close to 0, the weight of the gravity data is weak in the cooperative inversion and the final model will be the initial model given by Paul *et al.* (2001). When *c* is close to 1, the $V_p$ models become unstable in the iterative process. This factor must be chosen in order to obtain a trend of gravity variance reduction parallel to the seismic one, which is controlled by the damping factor of the seismic inversion.

The correlation length $\lambda$ introduces a forced correlation between the physical properties of two density prisms. If it equals to zero, the *a priori* information between two prisms is not correlated (this is the case in the seismic inversion), and when $\lambda$ is different from 0 the prisms are correlated and the correlation increases with $\lambda$. The greater the correlation length $\lambda$ is, the larger the structures and the anomaly wavelengths are.

## 3.2 Testing several parameters and sequences

In a first step, we search for the best sequence of inversions that explains the two data set and leads to density and $V_p$ models compatible with previous regional studies. We tested several sequences by varying the number of iterations *n* for each seismological inversion. As *n* = 6 is the number of iterations used by Paul *et al.* (2001) to obtain their final model, using this value amounts to performing the first gravity inversion directly from the final model. Whatever the values of the parameters *c* and $\lambda$, the final model for *n* = 6 cannot explain the gravity data without considering unreliable velocity values. The consistency between the two data set is obtained with a low *c* value that leads to velocities greater than 9.0 km s⁻¹.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top: variance of the seismic data versus the number of iteration. Crosses correspond to the LET performed by Paul *et al.* (2001) Triangles and dots corresponds to this study (triangle: variance after the seismic inversion, dot: variance after the gravity inversion). Paul *et al.* (2001) and this study have parallel decreasing of the variance indicating that the cooperative inversion explains the seismic data and the simple seismic inversion. Bottom: same thing for the standard deviation of the gravity data. The decrease of the variance of the seismic data and of the standard deviation of the gravity data are parallel, indicating that the value of *c*, the attenuation factor, is adapted to the studied case.</figcaption></figure>

After testing many sequences, it appears that the best one is $n = 1$. The optimal model is obtained for an attenuation factor $c = 0.10$ and a correlation length $\lambda = 10$ km. Fig. 6 documents how the two sets of data are progressively explained by the sequential inversion. The curves showing variance reductions in traveltime and gravity data have similar shapes, indicating that the chosen value for $c$ is well adapted to the case studied. The standard deviations for gravity and seismological data do not change significantly after the sixth iteration. The variance of the time delays decreases from 0.052 to 0.021 s². The standard deviation for the gravity field reduces from 58 to 10 mgal. Note that, after six iterations of LET, Paul *et al.* (2001)'s final model variance is 0.019 s²,

which is smaller than what we obtain here. However, the difference is small and our model explains gravity data, whereas their model does not.

Fig. 7 shows the computed residual gravity anomaly for the final model (left) and the difference between observations and computations (right). It documents that our model fits the observed gravity field very closely. Only very short-wavelength anomalies are not explained because of the large horizontal size of the cells, and perhaps to the use of an erroneous uniform value of the density for topographic correction.

## 4 DISCUSSION

It is well known that interpretation of a gravity residual field by a linear inversion method leads to a non-unique solution, which may be described by its average and its covariance matrix from a stochastic point of view. In our case, the fundamental non-uniqueness of the density solution is reduced by integrating a priori velocity information and Birch's linear relationship between velocity and density.

With our sequential approach, we found a 3-D velocity model that nearly explains the traveltimes and the LET model. Moreover, this model also explains the gravity anomalies while the LET one does not. Map-view slices of this final velocity model are shown in Fig. 8 and four cross-sections along lines A, B, C and D are shown in Fig. 9 (left). A comparison with the result of the LET (Figs 2 and 9, right) documents that the LET model does not explain the gravity data because of its underestimation of the lateral extent of the velocity anomalies. This is clearly the case for the high-velocity anomaly corresponding to the Ivrea Body. Some structures imaged by the sequential inversion are not found by the LET. This is, for example, the case for the deep high-velocity body located beneath the internal zone west of the Ivrea Body at 25–30 km depth (Fig. 9, cross-section A). The resolution of the local earthquake tomography changes dramatically across the model (Paul *et al.* 2001). To first order, it depends on the number of rays that cross the area. Conversely, the resolution of the gravity inverse problem decreases with depth and does not change within a given layer. During the sequential inversion, the gravity inversion steps only slightly modify the well-resolved cells of the $V_p$ model, whereas the poorly resolved cells are modified significantly. This is the way the sequential inversion proceeds to explain the gravity observations without degrading the fit to the traveltime observations. As the distribution of the earthquakes is very heterogeneous (all events with foci deeper than 10 km concentrate in the eastern part of the model under the Dora Maira Massif and Po plain) the resolution of the LET is low for a large part of the model. In the map-view at depths equal to or larger than 15 km (Fig. 2), the well-resolved areas are only those inside the thick black lines that correspond to the central part of the Ivrea Body.

Several structures of the LET model are strongly modified by the sequential inversion process. The east–west extension of the Ivrea Body at depths larger than 20 km is doubled, particularly in the northern part of the model (sections A and B in Fig. 9). The model explains the large E–W extension of the Ivrea gravity anomaly and cannot efficiently change the density estimated *a priori* from the ‘well-resolved’ velocity of the upper layers. In the southern part, which is characterized by a higher seismicity level, the Ivrea body is only slightly modified. This high-density, high-velocity structure appears to be a more or less vertically continuous body, between 10 and 30 km depth. Its maximum velocity is ~7.4 km s⁻¹ at 10 km depth and 7.9 km s⁻¹ at 30 km depth, indicating that it might be a
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left: synthetic Bouguer anomaly computed from the final sequential model. Right: difference between the observed and synthetic Bouguer anomaly. The average difference is close to zero for the entire study area.</figcaption></figure>

slice of upper-mantle wedging westward into the European Alpine crust (Ménard & Thouvenot 1984). A high-density–velocity zone is also observed beneath Dora Maira from the top of the Ivrea Body to the surface (Fig. 9, sections A and B). These shallow large densities are not associated with the Dora Maira massif, as its mean density is 2700 kg m⁻³ (Rey et al. 1990). Such a high-density zone was already demonstrated north of the region under study beneath the Sesia–Lanzo massif along the ECORS-CROP profile (Rey et al. 1990). It might be associated with high-grade metamorphic mafic rocks or with an upward extension of the upper mantle Ivrea Body.

Beneath the internal zone, at depths greater than 20 km, a high-velocity zone is revealed by the sequential inversion. This zone could be associated with slice of high metamorphic lower crust as it is documented on the cross-section of the Western Alps proposed by Schmidt & Kissling (2000). In section A (Fig. 9) we show the location of the high *V*<sub>p</sub> contrast reflector revealed by wide-angle seismic experiments near the section (ECORS-CROP 1989a,b). This reflector is located at the top of the high-velocity zone. Analysing the geophysical results along the ECORS-CROP profile, Nicolas et al. (1990) interpreted this zone as a large upper-mantle unit. This is a new evidence for stacking of upper mantle and (or) crustal slices beneath the internal Alps.

In the external domain and beneath the Digne and Castellane nappes, the upper crust is characterized by a 10 km thick low-velocity zone (Fig. 9, section D1). From the LET to the sequential models, the thickness of the low-velocity layer has been doubled. According to the geological cross-sections proposed by Ritz (1991), the thickness of the nappe piles is ~5 km. We suggest that the lower 5 km of low-velocity material could correspond to a hidden Permo–Carboniferous basin proposed by Ménard & Molnar (1988). We do not discuss the near-surface velocities values in the south-western part of this study further, because there are no other seismic data, apart from that of Paul *et al.* (2001) in this region.

## 5 CONCLUSION

Computation of a theoretical gravity anomaly from the LET seismic model of Paul *et al.* (2001) points out that this model is not able to explain the observed gravity field. To solve this problem, *P*-wave arrival times and gravity data are inverted sequentially five times. The least-squares optimal solution jointly satisfies the two sets of data and shows that the LET model underestimates the extent of the structures. In the ‘resolved’ parts of the model, the seismic tomography defines the location and the extreme values of the velocity anomalies well. However, inversion of the gravity field completes the *P*-wave tomography in crustal domains where rays are sparse.

The main features of the sequential inversion model are:

1. the Ivrea Body is larger than the LET one and vertically continuous under 10 km depth. It is roofed by relatively high-velocity material up to the surface;

2. a high-density-high-velocity zone is observed in the internal zone at depths greater than 25 km. Its top coincides with a wide-angle reflector observed from the ECORS-CROP seismic experiment;

3. a 10 km thick low-velocity zone under the Digne and Castellane nappes may include a 5 km deep Permo–Carboniferous basin that is hidden by the nappes.

The sequential method is a straightforward and efficient method of coupling tomography and gravity interpretations in the crustal continental domain with a high seismicity level and where velocity and density contrasts are high. The interest of the method is the complementarity of the two data sets. Indeed, this experiment shows that the well-resolved area of the LET model does not change and the gravity data bring fruitful information for the other parts of the model. The main difficulty consists in estimating the parameters introduced in the algorithm to obtain a convergence of the iterative process and ‘realistic’ models. In this study, this is performed during the gravity inversion by way of a priori information on the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Smoothed P-wave velocity maps obtained by sequential inversion. This model has to be compared with the model of Fig. 2. The main features of this model are: (1) the broadening of the north–south high-velocity anomaly corresponding to the Ivrea Body; (2) the existence of a high-velocity zone west of the Ivrea Body at depths greater than 20 km; and (3) the thickening of the low-velocity anomaly beneath the Digne and Castellane nappes visible from the surface up to 10 km depth (5 km in the LET model). This structure could indicate the existence of a hidden Permo-Carboniferous basin under the nappes.</figcaption></figure>

density contrasts deduced from density–velocity relationship and the confidence of the LET model. The sequential method could be applied to many other geophysical problems; for example, by combining wide-angle seismic data and gravity data. For instance, the inverse method developed by Zelt & Smith (1992) for the interpretation of the wide-angle seismic experiments could be coupled with gravity inversion within the framework of the sequential approach proposed in this paper.

## ACKNOWLEDGMENTS

This work was sponsored by the GéoFrance 3-D project (CNRS, MENRT and BRGM). We would like to thank all the participants of the gravity and seismological experiment who helped during the fieldwork to make these experiments successful. We also gratefully acknowledge U. Achauer and R. Keller for many helpful suggestions.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Vertical cross-sections along lines A, B, C and D of the map. A1, B1, C1 and D1 correspond to the sequential model (gravity + seismology) while A2, B2, C2 and D2 correspond to the LET model (seismology). The white line indicates the location of a deep reflector imaged by the ECORS-CROP wide-angle experiment. FPT, Frontal Penninic Thrust; DM, Dora Maira; ARG, Argentera; IB, Ivrea Body; CN, Castellane Nappe.</figcaption></figure>
## 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. Theor. Appl.*, **21**, 83, 149-157.

Barton, P.J., 1986. The relationship between seismic velocity and density in the continental crust—a useful constraint?, *Geophys. J. R. astr. Soc.*, **87**, 195-208.

Berckhemer, H., 1968. Topographie des Ivrea-Körpers, abgeleitet aus seismischen une gravimetrischen Daten, *Schweiz. Mineral. Petrogr. Mitt.*, **48**, 235-246.

Birch, F., 1961. The velocity of compressional waves in rocks to 10 kilobars, part 2, *J. geophys. Res.*, **66**, 2199-2224.

Choudhury, M., Giese, P. & de Visintini, G., 1971. Crustal structure of the Alps: some general features from explosion seismology, *Bull. Geofis. Theor. Appl.*, **13**, 51-52, 211-240.

Closs, H. & Labrouste, Y., 1963. Recherches sismologiques dans les Alpes Occidentales au moyen de grandes explosions en 1956, 1958 et 1960, in *Mém. Coll. Année Géophysique Int.*, Vol. 12, CNRS, Paris.

Coron, S., 1963. Aperçu gravimétrique sur les Alpes Occidentales: Année géophysique internationale, in *Mém. Coll. Année Géophysique Int.*, Vol. 12, CNRS, Paris.

Eberhart-Philips, D., 1993. Local earthquake tomography: earthquake source regions, in *Seismic Tomography: Theory and Practice*, pp. 613-643, eds Iyer, H.M. & Hirahara, K., Chapman & Hall, London.

ECORS-CROP Gravity Group (Bayer, R., Carozzo, M.T., Lanza, R., Miletto, M. & Rey, D.), 1989a. Gravity modelling along the ECORS-CROP vertical seismic reflexion profile through the western Alps, *Tectonophysics*, **162**, 203-218.

ECORS-CROP Deep Seismic Sounding Group, 1989b. A new picture of the Moho under the western Alps, *Nature*, **337**, 249-251.

Evans, J.R., Eberhart-Philips, D. & Thurber, C.H., 1994. User's manual for SIMULPS12 for imaging $V_p$ and $V_p/V_s$: a derivative of the 'Thurber' tomographic inversion SIMUL3 for local earthquakes and explosions, *US Geol. Surv. Open File Rep.*, 94-431.

Giese, P. & Prodehl, C., 1976. Main features of crustal structures in the Alps, in *Explosion Seismology in Central Europe*, pp. 347-375, eds Giese, P., Prodehl, C. & Stein, A., Springer, Heidelberg.

Glaznev, V.N., Raevsky, A.B. & Skopenko, G.B., 1996. A three-dimensional integrated density and thermal model of the Fennoscandian lithosphere, *Tectonophysics*, **258**, 15-33.

Groupe de Recherche Géofrance 3D, 1997. Géofrance 3D: l'imagerie géologique et géophysique du sous-sol de la France, *Mém. Soc. Géol. France*, **172**, 53-71.

Kaufmann, R.D. & Long, L.T., 1996. Velocity structure and seismicity of southeastern Tennessee, *J. geophys. Res.*, **101**, 8531-8542.

Kissling, E., Ellsworth, W.L., Eberhart-Philips, D. & Kradofler, U., 1994. Initial reference models in local earthquake tomography, *J. geophys. Res.*, **99**, 19635-19646.

Lees, J.M. & VanDecar, J.C., 1991. Seismic tomography constrained by Bouguer gravity anomalies: applications in western Washington, *Pure appl. Geophys.*, **135**, 31-52.

Lines, L.R., Schultz, A.K. & Treitel, S., 1988. Cooperative inversion of geophysical data, *Geophysics*, **53**, 8-20.

Ludwig, J.W., Nafe, J.E. & Drake, C.L., 1970. Seismic refraction, in *The Sea*, Vol. 4, pp. 53-84, ed. Waxwell, A.E., Willey, New York.

Masson, F., Verdun, J., Bayer, R. & Debeglia, N., 1999. Une nouvelle carte gravimétrique des Alpes Occidentales et ses conséquences structurales et tectonique, *C.R. Acad. Sci. Paris*, **329**, 865-871.

Ménard, G. & Molnar, P., 1988. Collapse of a Hercynian Tibetan Plateau

into a late Paleozoic European Basin and Range province, *Nature*, **334**, 235-237.

Ménard, G. & Thouvenot, F., 1984. Ecaillage de la lithosphère européenne sous les Alpes Occidentales: arguments gravimétriques et sismiques liés à l'anomalie d'Ivrea, *Bull. Soc. Géol. Fr.*, **26**, 875-884.

Morelli, G., 1963. Groupe d'études des explosions alpines: Mémoire collectif, Closs and Labrouste editors, année géophysique internationale, XII, 2, CNRS, Paris.

Nafe, J.E. & Drake, C.L., 1957. Variation with depth in shallow and deep water marine sediments of porosity, density and the velocities of compressional and shear waves, *Geophysics*, **22**, 523-552.

Nicolas, A., Polino, R., Hirn, A., Nicolich, R. & ECORS-CROP Working Group, 1990. ECORS-CROP traverse and deep structure of the western Alps: a synthesis, *Mém. Soc. Géol. Fr.*, **156**, 15-27.

Niggli, E., 1946. Über den Zusammenhang zwischen der positiven Schwereanomilie am Südfuss der Westalpen und der Geisteinzone von Ivrea, *Eclog. Geol. Helv.*, **39**, 211-220.

Oppenheimer, D.H. & Herkenhoff, K.E., 1981. Velocity density properties of the lithosphere from three-dimensional modeling at the Geysers-Clear Lake region, California, *J. geophys. Res.*, **86**, 6057-6065.

Paul, A., Cattaneo, M., Thouvenot, F., Spallarossa, D., Béthoux, N. & Fréchet, J., 2001. A three-dimensional crustal velocity model of the southwestern Alps from local earthquake tomography, *J. geophys. Res.*, **106**, 19367-19389.

Perrier, G., 1973. Structure profondes des Alpes Occidentales et du Massif Central Français, *PhD thesis*, University Paris.

Rey, D., 1989. Structure crustale des Alpes occidentales le long du profil 'ECORS-CROP' d'après la sismique réflexion et le champ de pesanteur, *PhD thesis*, University Montpellier.

Rey, D. *et al.*, 1990. Gravity and aeromagnetic maps on the Western Alps: contribution to the knowledge on the deep structures along the ECORS-CROP seismic profile, in *Deep Structure of the Alps, Mem. Soc. Geol. Fr.*, **156**, 107-121.

Richard, V., Bayer, R. & Cuer, M., 1984. An attempt to formulate well-posed questions in gravity: application of linear inverse techniques to mining exploration, *Geophysics*, **49**, 1781-1793.

Ritz, J.F., 1991. Evolution du champ de contraintes dans les Alpes du Sud depuis la fin le l'Oligocène: Implications sismotectoniques, *PhD thesis*, University Montpellier.

Roure, P., Choukroune, P. & Polino, R., 1996. Deep seismic reflection data and new insights on the bulk geometry of the mountain ranges, *C. R. Acad. Sci. Paris*, **322**, 345-359.

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.

Thouvenot, F. & Perrier, G., 1981. Seismic evidence of a crustal overthrust in the Western Alps, *Pure appl. Geophys.*, **119**, 163-184.

Thurber, C.H., 1983. Earthquake location and three-dimensional crustal structure in Coyote Lake area, Central California, *J. geophys. Res.*, **88**, 8226-8236.

Van de Meulebrouck, J., Bayer, R. & Burg, J.P., 1984. Density and magnetic tomography of the upper continental crust: application to a thrust area of the french Massif Central, *Ann. Geophys.*, **2**, 5, 579-592.

Waldhauser, F., Kissling, E., Ansorge, J. & Mueller, S.T., 1998. Three-dimensional interface modelling with two-dimensional seismic data: the Alpine crust-mantle boundary, *Geophys. J. Int.*, **135**, 264-278.

Zelt, C.A. & Smith, R.B., 1992. Seismic traveltime inversion for 2D crustal velocity structure, *Geophys. J. Int.*, **108**, 16-34.

Zeyen, H. & Achauer, U., 1997. Joint inversion of teleseismic delay times and gravity anomaly data for regional structures: theory and synthetic examples, in *NATO Science Series, Partnership Subseries 1, Disarmament Technologies*, **17**, 155-168.