A three-dimensional model of the Pyrenean deep structure based on gravity modelling, seismic images and petrological constraints

Summary

The present study aims at building a 3-D density model of a whole collisional mountain belt at the lithospheric scale. The geometry of internal structures is constrained by seismic studies, while density values are imposed by petrological considerations. The model must fit both the Bouguer and geoid anomalies. This approach is applied to the Pyrenees, for which numerous geophysical data are available. The crustal structure is constrained from the results of seismic profiles. The most important crustal feature is the large Moho jump which coincides with the surface expression of the North Pyrenean Fault, considered as the suture between the Iberian and Eurasian plates. The Eurasian crust is about 30 km thick, whereas the Iberian crust is 50 km thick in the central part of the range. Thick sedimentary basins at the Piedmont, and two large heavy intracrustal bodies located in the western and central parts of the range, are also important features of the crustal model. At the lithospheric scale, a low-velocity strip under the North Pyrenean Fault down to 80–100 km depth, interpreted as the subduction of the Iberian lower crust, has been revealed by tomographic models. The density of this structure has been constrained using petrological arguments, which assume that the crustal material is transformed into eclogite at about 12–15 kbar (40–50 km depth). An eclogite-type density value, much higher than the mantle density, is thus imposed from 50 to 100 km depth.

Theoretical Bouguer and geoid anomalies are computed for this 3-D density model using a summation of density block anomalies. A good agreement is found with observations, giving strong support to the proposed model. The crustal structure, as seen by seismic profiles, explains most of the observed long-wavelength negative Bouguer anomaly, as well the local anomalies related to the intracrustal bodies. A possible mantle origin of these bodies is ruled out. On the other hand, the transformation of the subducted lower crust into dense eclogite provides a successful fit to the long-wavelength positive geoid anomalies. In addition, the water released during eclogite formation may help to explain the high conductivity zone detected down to 80 km south of the North Pyrenean Fault. The 3-D model proposed here thus reconciles petrological considerations and geophysical data.

Bouguer anomalies, density, eclogite, geoid, Pyrenees, tomography

Introduction

The Pyrenean range results from the north–south convergence of the Iberian and Eurasian plates, which collided about 65 Ma (see Choukroune 1992; for a review). The suture between the two plates, the North Pyrenean Fault (NPF), is a major tectonic accident running north of the range from the Atlantic Ocean to the Mediterranean Sea (Fig. 1). The convergence episode which lead to the rise-up of the range followed an extensive episode related to the opening of the Bay of Biscay: the left-lateral strike slip motion induced crustal thinning and the opening of pull-apart basins (see, e.g. Beaumont et al. 2000). The transition from extensive to compressive regimes was made through a transpressional episode (75–100 Myr ago) that resulted in the rising up of small lherzolite massifs, observed in the North Pyrenean Zone all along the NPF (Vielzeuf & Kornprobst 1984). To the north, the North Pyrenean Zone overrides the Aquitaine Basin along the North Pyrenean Frontal Thrust. South of the NPF, the Palaeozoic Axial Zone is in a large part inherited from a previous orogeny. It includes large granitic massifs whose uplift has been reactivated during the last 65 Ma, leading to summits up to 3400 m. A few allochthonous Palaeozoic massifs (Souquet & Peybernès 1987) are also present in the North Pyrenean Zone, in particular the North Pyrenean Massifs in the central part of the range, and the Labourd Massif to the west. South of the Palaeozoic Axial Zone, the South Pyrenean Zone, made of Mesozoic and Cenozoic sediments, overthrusts the molasse of the Ebro Basin to the South.

Figure 1

Main structural units of the Pyrenean range (redrawn after Choukroune 1992). NPFT: North Pyrenean Frontal Thrust; NPM: North Pyrenean Massifs; SPT: South Pyrenean Thrust. Light grey indicates zones of Palaeozoic material, dark grey indicates zones of Quaternary volcanic rocks. Also shown with dots are two regions exhibiting prominent anomalies in gravity and seismic data (see text). The locations of the four cross-sections of Fig. 2 are shown.

In the last 20 years, a large number of refraction and reflection seismic profiles have been performed in or around the Pyrenees (e.g. Hirn et al. 1980, 1984; Gallart et al. 1980, 1981; Daignières et al. 1981, 1982, 1994; Séguret & Daignières 1986; Pinet et al. 1987; ECORS Pyrenees Team 1988; Fernandez-Viejo et al. 1998). A large amount of data is therefore available concerning the thickness and seismic structure of the Pyrenean crust. These data have revealed the existence of an important Moho jump, located below the surface expression of the NPF. On the Eurasian side, the crust is about 28–30 km thick, whereas it reaches 50–55 km on the Iberian side in the central part of the range. The thickness of the Iberian crust decreases at both ends of the range, it is of the order of 40 km at the Atlantic side, and it does not exceed 35–40 km near the Mediterranean. More recently, a tomographic image of the Pyrenees has been obtained down to 150 km from the inversion of the propagation times of both local and teleseismic events (Souriau & Granet 1995). At crustal level, some bodies of high P- and S-wave velocities have been detected in the North Pyrenean Zone. Their locations coincide with those of prominent positive Bouguer anomalies (de Cabissole 1989; dotted zones in Fig. 1). At the lithospheric scale, a 30 km wide strip of negative P velocity (− 2 per cent) is observed south of the NPF down to a depth of 80–100 km. It has been interpreted as Iberian lower crust material subducting under the Eurasian plate. This anomaly has also been detected from a magnetotelluric survey in the central Pyrenees (Pous et al. 1995), and is possibly associated with partial melting (Ledo et al. 2000).

Seismic reflectors have been used as constraints to build density structures that are consistent with observed Bouguer anomalies. It has been shown that the crustal thickening may explain the very large negative Bouguer anomaly observed south of the NPF (Torné. 1989). On the other hand, the large anomalous body located in the upper crust in the western part of the North Pyrenean Zone has been modelled as a lower crustal block, in order to fit the local Bouguer anomaly (Grandjean 1994).

The problem of the global equilibrium of the range has rarely been considered. There is, however, an apparent inconsistency between the results of gravity and tomography: density modelling can explain the observed Bouguer anomalies as resulting from a thickened crust of 50 km, whereas seismic tomography observes a subduction of this crust down to 80 or 100 km depth. On the other hand, the important Moho jump below the NPF results in an overcompensation of the elevations south of this fault, according to Airy's model. A subcrustal structure is thus necessary to ensure compensation (e.g. Zeyen & Fernandez 1994; Ledo et al. 2000).

The goal of the present paper is to propose a new density model able to reconcile the tomographic images and the gravity anomalies, using petrological arguments. Some previous models have been either limited to local anomalies (Grandjean 1994) or restricted to 2-D models (e.g. Brunet 1986; de Cabissole 1989; Torné. 1989; Ledo et al. 2000). The present paper proposes a 3-D density model of the whole Pyrenean range, including most of the available geophysical data. No gravity model has yet taken advantage of the 3-D tomographic models of Souriau & Granet (1995), which provide important geometrical and geophysical constraints. Furthermore we aim to fit both observed Bouguer and geoid anomalies: Bouguer anomalies provide information mostly about superficial structures, whereas geoid anomalies may give insight into deep, large-scale features, thus in particular into subcrustal structures.

In the following sections, we will give a description of the 3-D density model. Results for Bouguer anomalies and geoid heights are then given and compared to the observed values. A comparison with previous density models will be addressed and finally petrological interpretations will be discussed, as well as the consequences of the presence of partial melt and free water under the Pyrenees.

3-d density modelling of the pyrenees

The 3-D density structure of the Pyrenees is modelled as an assemblage of blocks, each with a constant density. This section details the geophysical data used to define block geometries. 2-D views of the blocks can be found in Fig. 2. Theoretical Bouguer anomalies and geoid heights are then computed, using a classical algorithm which transforms the volume integral for each constant-density block into a surface integral, and sums over all the blocks (e.g. Talwani & Ewing 1960). With such a procedure, only a density contrast of the block with respect to the surrounding medium is needed. We define a 1-D reference medium with only three layers: upper crust, lower crust and mantle, with densities of 2.75, 2.93 and 3.28 g cm⁻³, respectively. These values come either from direct measurements, or from density–seismic velocity systematics (de Cabissole 1989). The boundary between upper and lower crust is fixed at 15 km depth. This choice has absolutely no influence on the final results, since no lateral density variation is imposed within the lower crust (between 15 and 30 km depth). The gravity and geoid anomalies computed in the next sections will result from density perturbations with respect to this 1-D model. Observed Bouguer anomalies are taken from Bayer et al. (1996) and observed geoid heights are from Denker et al. (1997).

Figure 2

Observed (solid lines) and computed (dashed lines) anomalies as a function of latitude (°) along the four south to north cross-sections located in Fig. 1. For each cross-section, the top graph gives Bouguer anomalies Δg (mGal), the middle graph gives geoid anomalies ΔN (m) and the bottom graph shows the intersections of the modelled 3-D blocks with the corresponding cross-section. Numbers correspond to the different densities. The dark arrow on top of each panel gives the location of the NPF.

Topography

Topography is required for geoid calculations, to include the effects of masses above the sea level. To obtain the topography of the Pyrenees, we use the ETOPO5 data file with a resolution of 5 min×5 min (NOAA 1988). From this map the shape of the mountain belt is approximated using nine truncated pyramidal blocks up to 2200 m high. Higher altitudes form only local summits that do not significantly contribute to the geoid. A constant density of 2.67 g cm⁻³ is assumed inside these blocks. This value is slightly lower than the assumed density of the upper crust (2.75 g cm⁻³), but it corresponds to the value used to build the observed Bouguer anomalies (Bayer et al. 1996). It is therefore required in the geoid computation to obtain a consistent picture.

Sediments

The Pyrenees are surrounded by sedimentary basins—the Aquitaine Basin to the north and the Ebro basin to the south. The sediment thickness varies greatly, with a maximum of more than 10 000 m near Pau, to the northwest of the range. The geometry and thickness of these sediments are taken from compilations provided by BRGM et al. (1974) from seismic and well log data. Seven blocks form this sediment cover, with different densities depending on their depth and age: we follow here an approach similar to that of Grandjean (1994), giving higher densities to older and deeper sediments, which have undergone more compaction. We assume Tertiary sediments for the entire Ebro basin and the most superficial part of the Aquitaine basin (sediments down to a depth of 5000 m) to have a density of 2.6 g cm⁻³. Upper cretaceous sediments (between 2000 and 5000–7000 m depth) are given a density of 2.7 g cm⁻³. Finally, sediments from 7000 to 10 000 m near Pau are given a density of 2.75 g cm⁻³, thus they do not exhibit any density contrast with respect to the average upper crust.

Crustal structure

The crustal structure on the European plate (north of the NPF) is relatively simple: seismic profiles show its thickness to be constant at 30 km (Gallart et al. 1980, 1981; Daignières et al. 1982; Pinet et al. 1987; Ecors Pyrenees Team 1988; Daignières et al. 1994; see also Daignières et al. 1981 for a longitudinal profile). The crustal thinning due to the opening of the Bay of Biscay appears further west (e.g. Fernandez-Viejo et al. 1998). Some crustal thinning has been included on the very northwestern edge of our model, with a Moho depth decreasing northwards to a value of 25 km (Pinet et al. 1987; see Lefort 1993 for a synthesis of the Bay of Biscay).

A synthetic view of the Iberian crustal thickness can be found in Verges et al. (1995). The Iberian Moho away from the Pyrenees is 33 km deep, and it plunges down towards the NPF, where it reaches 40 km in the western Pyrenees (Daignières et al. 1982), 50–55 km in the central part (Hirn et al. 1980; ECORS Pyrenees Team 1988) and 30–40 km in the eastern part (Gallart et al. 1980; Daignières et al. 1982). All these depths are confirmed by the longitudinal profile performed south of the NPF (Gallart et al. 1981). To account for these lateral differences in Moho depths in a convenient way, we divided the Iberian crust into three blocks (with equal densities), bounded by the NPF northwards and by the South Pyrenean Thrust southwards, with laterally varying depth as described above. The western block can be seen in cross-section #1, the central block is cross-cut by profiles #2 and #3, and the eastern block is visible along cross-section #4 (Fig. 2).

Two different layered structures for the Iberian crust have been considered in previous gravity models: Torné. (1989) included a mid-crustal layer with a density of 2.80 g cm⁻³, whereas Grandjean (1994) proposed the same upper/lower crustal density structure (densities 2.75/2.90 g cm⁻³) for both the European and the Iberian sides, the lower crust being much thinner on the European side (less than 10 km thickness, compared to an almost 20 km thick Iberian lower crust). Such differences are poorly constrained by geophysical data, so it was decided to assume the same lower crustal geometry on both sides of the range, as did Ledo et al. (2000) in their 2-D model. We also assumed that the discontinuity between the upper crust and lower crust does not exhibit any significant topography.

Crustal bodies

Two main anomalous regions within the Pyrenean crust have been recognised from geophysical data. The first one is located below the Mauléon basin and Labourd massif, the second one south of the city of Saint-Gaudens (see dotted regions in Fig. 1). Hereafter they are referred to as the ‘Labourd’ and ‘Saint-Gaudens’ anomalies. These anomalies can be seen in the Bouguer anomaly map (Fig. 3), and are also observed in the tomographic model of Souriau & Granet (1995). The Labourd anomaly has been extensively studied by Grandjean (1994) from a joint analysis of gravity data and teleseismic residual times. He proposed a density of 2.9 g cm⁻³ together with a P-wave velocity of 7.0 km s⁻¹ (16 per cent faster than his upper crust velocity of 6.0 km s⁻¹). In the large-scale tomographic model of Souriau & Granet (1995), it appears slightly smaller and weaker, with only a +6 per cent P-wave anomaly (which means an absolute velocity of 6.3–6.6 km s⁻¹). This diminishing is possibly due to the fact that the anomaly lies on the edge of the tomographic model and therefore suffers from poor ray sampling. To define the geometry of this anomaly, we made a compromise between the geometries proposed by Grandjean (1994) and Souriau & Granet (1995): it extends from 1 to 12 km depth, becomes wider below 7 km depth (below sediments) and then narrows down to 12 km depth.

Figure 3

Maps of Bouguer anomalies (in mGal, top panels), and geoid anomalies (in m, bottom panels). In each case the top map shows observed anomalies and the bottom map shows computed anomalies. Thick lines: main faults; dashed lines: border. The four cross-sections of Fig. 2 are located.

For the Saint-Gaudens anomaly, we have used the geometry of the seismic anomalies obtained in the tomographic model: this region is in the centre of the model, and one can be much more confident of the geometry proposed by seismic tomography than for the Labourd anomaly. The proposed block is about 60 km long, 25 km wide and extends from 7 to 15 km depth, where it is much smaller.

The petrological nature of these anomalies is still unknown, although they are generally attributed to lower crust or mantle bodies. Their densities have been considered as free parameters in the model, in order to gain insight into their composition.

Iberian lower crust subduction

Seismic tomography (Souriau & Granet 1995) revealed a 30 km wide strip of low-velocity anomalies (− 2 per cent) extending down to 80–100 km depth in the central and eastern Pyrenees. This structure has been ascribed to the subduction of the Iberian lower crust during convergence. Such a feature has indeed been reproduced by thermomechanical models (Chéry et al. 1991; Beaumont & Quinlan 1994), and is also consistent with magnetotelluric data (Pous et al. 1995). Such a vertical feature could not be seen in seismic reflection, and former density models did not include it if we except the recent 2-D model of Ledo et al. (2000). Here it is accounted for by using the geometry found in tomography: a first constant-density block makes the transition from the thickened lower crust to the vertical subducting strip, at depths of 50–60 km, and a second block represents the vertical strip itself down to 100 km depth. This strip is about 200 km long and 30 km wide (Souriau & Granet 1995). To constrain the densities of these two blocks, we use arguments based on the petrological evolution of lower crust material when depth (hence pressure) increases. The lower crust is generally composed of interstratified layers of granulitic and mafic rocks. At about 12–15 kbar (around 40–50 km depth), along a modelled geotherm for the Pyrenees (Zeyen & Fernandez 1994), mafic rocks gradually transform into eclogite, a very dense assemblage of garnet and pyroxenes (Fig. 4). The transition to eclogite facies must be very smooth and spans a large depth interval (white zones in Fig. 4), since many reactions are involved: plagioclases and amphiboles disappear whereas garnets and pyroxenes form. Eclogites are very dense rocks with a density even higher than that of ultramafic (mantle) rocks: we use a value of 3.48 g cm⁻³, which represents an average value of different eclogites from different tectonic settings (Christensen & Mooney 1995). We propose a model where the first block, extending from 50 to 60 km depth, can be understood as a transitional region from crustal mafic rocks to pure eclogite: we impose an arbitrary intermediate density of 3.35 g cm⁻³ for this block (+ 0.07 g cm⁻³ with respect to the surrounding mantle). The density value and geometry of this intermediate block are largely speculative, and are just intended to make a smooth transition from the light horizontally thickening crust to the heavy vertically subducting eclogite strip (see Fig. 2, cross-sections #2 and #3). Note that the densities used here are not corrected for the pressure increase with depth. This effect would influence the 1-D reference model and the density anomalies in the same way, so that the density contrast would remain unchanged.

Figure 4

Main metamorphic facies of mafic crustal rocks (redrawn after Bucher & Frey 1994). A geotherm computed from heat flow measurements in the Pyrenees (Zeyen & Fernandez 1994) is also shown (dashed at upper crust conditions). The geotherm enters the eclogite facies at a pressure of about 12–15 kbar, corresponding to a depth of about 40–50 km.

Results

Bouguer anomalies

Fig. 3 (top panels) shows maps of observed and computed Bouguer anomalies. The most evident feature in the observed Bouguer anomalies is the long-wavelength negative anomaly south of the NPF (in blue), with a maximum amplitude of around −120 mGal. Overprinted to the north of this large-scale negative anomaly are two zones of smaller positive anomalies above the Labourd and Saint-Gaudens anomalies. In the description of the computed anomalies, the long-wavelength negative anomaly is examined first, then smaller-scale features are addressed. Four cross-sections through the 3-D model along the lines shown in Fig. 3 are displayed in Fig. 2.

The geometry of the negative anomaly is correctly fitted by our 3-D model, in both its north–south and east–west directions. On the edges of the mountain range, cross-sections #1 and #4 show that the geometry imposed by seismic reflection profiles explains correctly the Bouguer anomalies. On the western side (cross-section #1), the wedge-shaped lower crust (Daignières et al. 1994), becoming less and less thick westwards (Daignières et al. 1981), explains the smaller amplitudes, around −70 mGal. On the eastern side (cross-section #4), the thickened lower crust flattens horizontally towards the south and becomes thinner eastwards (Gallart et al. 1980; Hirn et al. 1980). This structure explains why the negative anomaly is wider than to the west, still with an amplitude larger than −100 mGal.

In the central part of the range, the negative anomaly in our model is made up of the sum of three contributions. Sediments have a small negative contribution on both sides of the axial zone, generally less than −30 mGal. The lower crust on the Iberian side gives a very strong negative contribution reaching −130 mGal. Finally, these strong negative contributions are counteracted by the positive contribution of the dense subducted lower crust, with a maximum anomaly around +34 mGal. Along profiles #3 and #4, the computed Bouguer anomaly matches the observed one with an excellent accuracy of better than 10 mGal. The fit is not as good for profiles #1 and #2. For this last profile, the northern flank of the negative anomaly is shifted southwards with respect to the observed profile by about 10 km. This feature is due to the slightly out-of-phase addition of the respective contributions of the thickened lower crust and the heavy subducting strip. The 3-D geometry of the proposed model exhibits a slight difference between the location of the thickened lower crust, deduced from seismic reflection profiles, and the location of the subducting strip, deduced from seismic tomography. However, a 10 km shift is on the edge of the resolution achievable by the different geophysical methods. Apart from these details, the results show that a good fit is obtained from the available crustal models and that the subduction of the Iberian lower crust down to 100 km depth, with eclogite transformation at a depth of around 50–60 km, is compatible with the large-scale Bouguer anomalies.

The long-wavelength picture of the Bouguer anomalies being correctly explained by the model, we can focus now on the density of small-scale crustal features. The two positive anomalies corresponding to the Labourd and Saint-Gaudens regions are correctly fitted by the 3-D density model, with input densities of 2.93 g cm⁻³ and 2.95 g cm⁻³, respectively. In the case of the Labourd anomaly, the observed anomaly of +50 mGal is matched by the model with a body much closer to the surface than in Grandjean (1994)'s geometry. The fit is obtained here thanks to the 3-D contribution of the heavy subducted material, which is only 20 km away to the east. Petrological implications of the densities imposed for the crustal anomalous blocks are addressed in the Discussion section.

Geoid heights

The most important feature of the present 3-D model is the presence of heavy material (eclogite) at great depth, between 50 and 100 km. We showed in the previous section that this heavy material is so deep within the Earth that its signature in Bouguer anomalies is relatively weak. Therefore, the good fit to observed Bouguer anomalies that has been shown here cannot be seen as a strong argument for the eclogite hypothesis. The great size of this eclogite strip (200 km long, 30 km wide, 50 km thick) results in a long-wavelength high-density anomaly. It should therefore give a strong signal in geoid anomalies, which are very sensitive to deep structures.

The bottom panels in Fig. 3 show maps of observed and computed geoid anomalies. Computed values are relative heights with respect to an arbitrary zero level. The observed anomalies are adapted from the quasigeoid of Denker et al. (1997). To extract the signal due to the Pyrenean structures, we have removed a 3-D, long-wavelength trend of the geoid undulations. It has been approximated by the weighted average of the geoid defined at each edge of the grid of computation (see Fig. 3). It mostly consists of a smooth upwards slope of 0.7 m km⁻¹ in the SSW direction.

The computed geoid heights globally match the observed values. The positive anomalies along the range axis, with maximum amplitudes of 3–5 m south of the NPF, are recovered. However, we note that the zone with positive anomalies, if defined by the region with ΔN greater than +1 m, is too narrow by about 20 km along profiles #1 and #2, and by about 10 km along profiles #3 and #4 (Fig. 2). Maximum amplitudes are correctly fitted, except along cross-section #3, where the observed maximum anomaly of 4.6 m is underestimated by 0.6 m. As for Bouguer anomalies, the computed positive geoid anomalies come from the summation of three main contributions: the topography has a strong positive incidence, of the order of +11.5 m; the thickened lower crust contributes negatively, with values around −10.5 m; and finally the heavy subducting eclogite adds a positive contribution, around +3.6 m. It is clear that a small shift of one of these contributions with respect to the others can easily affect the resulting global anomaly. In this respect, the fit obtained by the 3-D model is satisfactory.

Discussion

Comparison with previous density models

The most important feature of our 3-D model of density under the Pyrenees lies in the narrow strip of Iberian lower crust subducted along the NPF and transformed into dense eclogite. The geometry of this feature comes from results of seismic tomography, and is consistent with magnetotelluric data (Ledo et al. 2000) and also with some geodynamical models of the regional shortening. Proposed values vary between 100 and 165 km of north–south shortening (see Beaumont et al. 2000 and references therein). One point, however, needs to be clarified here: previous models of density under the Pyrenees did not include such a feature, but they did fit the Bouguer anomalies (de Cabissole 1989; Torné. 1989; Zeyen & Fernandez 1994). The fact is that the positive contribution of the subducted eclogite to the bulk Bouguer anomaly is relatively small. It represents about 20 per cent of the negative contributions of sediments and lower crust together. This means that although eclogite is much denser than the surrounding mantle, it is much too deep to create high Bouguer anomalies. From this point of view, our model provides an alternative for fitting Bouguer anomalies. However, fitting geoid anomalies is a much more powerful tool, since they are very sensitive to deep, long-wavelength features.

Very recently, Ledo et al. (2000) also obtained a good fit to both Bouguer and geoid anomalies along a 2-D north–south cross-section through the Pyrenees with a completely different model. They also included a strip of lower crust subducting into the mantle and extended this feature down to 80 km depth on the basis of their magnetotelluric measurements. However, they did not invoke a mineralogical transformation, and allocated a constant density (2.93 g cm⁻³) to the lower crust down to 80 km. They included a wide mantle wedge above the subducting lower crust, a feature which is made possible, but not required, by the seismic observations (Daignières et al. 1989). They also introduced a jump at the lithosphere–asthenosphere discontinuity. These two features give the positive contribution required to fit the positive geoid anomalies. Their study and ours nicely illustrate how two completely different density models can fit Bouguer and geoid anomalies. We mention here two lines of argument supporting the present 3-D density model. First, from a petrological point of view, there is no reason to hamper the metamorphism of lower crustal material brought to pressures in excess of 10 kbar. Transformation into very dense rocks seems unavoidable. Second, if the subducting strip was made of untransformed crustal material, then it would create a very large seismic contrast of about −12 per cent with respect to the surrounding mantle [taking the ‘mafic garnet granulite’ values of Christensen & Mooney (1995) for the lower crust]. An even higher velocity contrast would be obtained if partial melt is present, as proposed by Ledo et al. (2000) to account for the observed high conductivities. Such a huge seismic anomaly would have been detected in the tomographic analysis, which exhibits a good resolution in the corresponding depth range (Souriau & Granet 1995). An anomaly of only −2 per cent has been reported. Thus, although it is possible to explain Bouguer and geoid anomalies by a subduction of untransformed lower crust, it is inconsistent with seismic tomography and petrological constraints. The 3-D model of the geoid also confirms that the dipping structure is not present to the west of the range, as observed from the tomographic models.

Thermal considerations

Throughout the present paper a thermal origin of the seismic anomaly has been neglected. We outline here some arguments for the petrological interpretation. As the Iberian lower crust subducts under the European crust, cold material is penetrating into the warmer mantle. This subduction probably started soon after the collision between the Iberian and Eurasian plates 65 Myr ago (e.g. Choukroune 1992). We can roughly estimate the time needed for the subducting lower crust to re-equilibrate thermally with the surrounding mantle using a characteristic thermal diffusion length formula, δ = √κt, where κ is the thermal diffusivity (e.g. Turcotte & Schubert 1982). With δ = 15 km (half the thickness of the subducted lower crust in our model) and κ = 10⁻⁶ m² s⁻¹, a time of 7 Myr is required for the centre of the subducted crust to start warming up. One can therefore expect a nearly complete thermal re-equilibration after 65 Myr of convergence.

More elaborate models taking into account ongoing subduction have been created to investigate the thermal and elastic states of the Pyrenean crust and mantle (Chéry et al. 1991). In the latter study, the lower crust before compression is about 200 °C colder than the mantle. After 12 per cent shortening of the range corresponding to the present state, the conductive heat transfer has warmed up the subducting crust by about 75–100 °C, whereas the surrounding mantle is cooled down. From this model, it thus appears that the Iberian lower crust might be about 100 °C colder than the mantle, but not much more. Thermal computations based on heat flow measurements (Zeyen & Fernandez 1994) confirm this last conclusion. Converting this 100 °C into P-wave seismic velocity (with a conversion factor ∂ lnV_P/∂T = −7 × 10⁻⁵ °C⁻¹; see Christensen & Mooney 1995), one obtains an anomaly of +0.7 per cent. Therefore, the thermal contribution to the observed seismic contrast must be small.

From another point of view, one can imagine a combined effect of composition and temperature inducing the observed −2 per cent seismic contrast. For instance, we have advocated the necessity of eclogite transformation, because crustal metamorphic rocks have velocities about −12 per cent lower than mantle rocks. The argument was built without any thermal variation. The thermal contrast needed to decrease the −12 per cent seismic anomaly down to the observed value of −2 per cent must be of the order of −1000 °C, a value that can definitely be ruled out in the present context. Therefore, an alternative model to the presence of eclogite, based on thermal variations, seems unrealistic. The presence of eclogite, colder than the mantle by 100 °C (Chéry et al. 1991), will satisfactorily reconcile the models with the observed anomalies.

Petrological nature of the intracrustal bodies

Heavy bodies in the Pyrenean upper crust (Labourd Massif and Saint-Gaudens anomalies) have been attributed either to lower crustal material or to mantle material (e.g. Souriau & Granet 1995). This interpretation comes from the fact that both types of rocks crop out in the Pyrenees. We have reported in Fig. 5 the constraints we have from geophysical data and from modelling of these crustal bodies (grey area) using densities and P-wave velocities. P-wave velocities observed by tomography are of the order of 6.5–6.7 km s⁻¹, and Grandjean (1994) proposed a value of 7 km s⁻¹. Densities obtained by modelling and fitting of Bouguer anomalies are either 2.9 (Grandjean 1994), 2.93 or 2.95 g cm⁻³ (this study). To be conservative, we have extended the proposed density range to between 2.8 and 3.1 g cm⁻³. Clearly, geophysical results are inconsistent with a mantle nature for these bodies. Dunite has a P-wave velocity of 8.21 km s⁻¹ for a density of 3.31 g cm⁻³, values for pyroxenite are 7.69 km s⁻¹ and 3.27 g cm⁻³, and values for eclogite are 7.88 km s⁻¹ and 3.48 g cm⁻³, respectively (Christensen & Mooney 1995). These three points lie outside the range proposed for the anomalous bodies. Candidate rocks are metamorphic rocks (mafic granulite, amphibolite, greenschist facies basalt), but one should be cautious interpreting these results. Many parameters, such as the presence of water, a local thermal anomaly and the presence of anisotropy, can shift the position of the points in Fig. 5. However, the conclusion that crustal bodies in the Pyrenees cannot be of mantle nature seems definitely reliable. Thus they do not seem to be related to the outcropping lherzolites, which form very small-scale bodies all along the range, and which are not detectable in gravity or geoid anomalies with the large grid considered in the present study.

Figure 5

Density as a function of P-wave velocity (km s⁻¹) for different crustal and mantle rocks. Data points are from the compilation of Christensen & Mooney (1995), along an averaged geotherm, and for a depth of 10 km. ECL: eclogite; DUN: dunite; PYR: pyroxenite. Other points are for typical crustal rocks. Error bars are standard deviations around the average of experimental values. The grey area corresponds to the values of density and P-wave velocity derived by geophysical methods for the crustal anomalous bodies (Labourd and Saint-Gaudens anomalies, see text), indicating that these bodies cannot be of mantle origin.

Isostatic compensation of the range

A rough estimate of the mass balance shows that the Pyrenean relief is overcompensated at crustal level. For Airy's model, a mean topography of about 1000 m would be locally compensated by a lower crust root of about 8 km, a value significantly smaller than the 20 km observed for the Moho offset in the central part of the range. A model of flexure of a semi-infinite elastic plate loaded by the Pyrenean topography has been proposed to explain the equilibrium of the range (Brunet 1986). This model requires, however, an additional forcing whose origin was assumed to be a dense lithospheric root.

The model we have obtained shows that an Airy-type compensation is not ruled out, thanks to the heavy eclogite material which compensates the excess in crustal asymmetry. A 1 km topography with a 20 km Moho offset (mass deficit 0.2 kg m⁻³) is isostatically compensated by a 34 km denser root (mass excess 0.1 kg m⁻³). Such an equilibrium is compatible with the low uplift rate of the range, inferred from comparison of levellings over the past century (Fourniguet & Lenôtre 1986).

A small present-day convergence rate of 0.5–1 cm yr⁻¹ between the African and Eurasian plates at the longitude of the Pyrenees has been measured from space geodesy (Crétaux et al. 1998). It implies that the subduction continues at present, thus favouring a dynamical equilibrium in which the dense root may contribute to sustain the collision process (e.g. Houseman et al. 1981; Fleitout & Froidevaux 1982).

Water and partial melt under the Pyrenees

Recently, magnetolluric (MT) measurements have been conducted in the Pyrenees (Pous et al. 1995; Ledo et al. 2000). A high conductivity zone has been found below the NPF to depths of at least 80 km, in good agreement with the low-velocity strip observed by seismic tomography. High conductivity anomalies are generally attributed to the presence of free fluids, partial melt or graphite inclusions (Jones 1992). Of these, the presence of partial melt has been the favoured interpretation of MT measurements (Ledo et al. 2000; Partzsch et al. 2000). We propose in this section a review of the conditions required for the presence of partial melting and connected free water.

Fig. 6 shows the Pyrenean geotherm as computed from heat flow measurements (Zeyen & Fernandez 1994). Overprinted on the geotherm is the region where eclogite forms from the mafic lower crust, around depths of 45–55 km. The formation of eclogite is a dehydration process: hydrated minerals are not stable any more, and decompose into dry minerals plus free water (e.g. Bucher & Frey 1994). Therefore, in this depth range free water will be released into the system. The distribution of this fluid on the grain surface is determined by the dihedral angle Θ between the solid grains and the melt. If Θ > 60°, water does not form an interconnected network and stays trapped at grain boundaries. We have plotted the experimentally determined Θ = 60° isopleths for an aqueous fluid–forsterite–forsterite triple junction (Mibe et al. 1999 and references therein). The 50° and 70° isopleths can be taken as a confidence domain for the occurrence of fluid interconnectivity. The Pyrenean geotherm crosses the 60° isopleth at a depth of 60–70 km, meaning that the free water released by the dehydration is trapped within the subducting eclogite down to that depth. At greater depths, percolation of the fluid along grain boundaries is possible, allowing for a lower resistivity. Also plotted in Fig. 6 is the solidus of H₂O-saturated mantle peridotite (Kawamoto & Holloway 1997). The Pyrenean geotherm crosses the wet mantle solidus at a depth of around 70–80 km. Since water coming from eclogite formation forms a connected network at this depth, it can reach the surrounding mantle and trigger its partial melting. Note, however, that if the Pyrenean geotherm is correct, this partial melting can occur only at a depth of 70–80 km, which corresponds to the bottom end of the observed high conductivity region. Therefore, partial melt might not be the cause of high conductivities all the way down from 50 to 80 km depth. We suggest that eclogite formation can be a process inducing high conductivities. The experiments of Popp & Kern (1993) on different minerals (gypsum, carnallite and serpentinite) showed that dehydration reactions during prograde metamorphism are associated with discontinuous behaviour in electrical conductivity. Extrapolating these experiments to the dehydration of the lower crust material, one can suspect an abrupt increase of conductivity around 45–55 km depth. We therefore suggest that the observed high conductivity zone may have two different origins: water released by eclogite formation at depths of 45–70 km, and partial melting of the surrounding mantle at greater depths. These fluids are probably concentrated in thin regions (cf. the 10 km thick high conductivity zone, Ledo et al. 2000), and therefore are barely detectable by tomographic models.

Figure 6

Pyrenean geotherm based on computations from heat flow measurements (Zeyen & Fernandez 1994) as a function of depth and pressure. The thick dashed line gives the 60° isopleth of the dihedral angle at aqueous fluid–forsterite–forsterite triple junctions (Mibe et al. 1999). The 50° and 70° isopleths can be taken as error bars around the 60° curve. When the dihedral angle is smaller than 60°, water can form an interconnected network. Also shown is the solidus of mantle H₂O-saturated peridotite (Kawamoto & Holloway 1997), above which wet mantle may begin to melt. Shown in grey along the Pyrenean geotherm are zones of dehydration due to the transformation of crustal material into eclogite and water interconnection.

Conclusions

Available geophysical data have been used to build a 3-D density model of the deep Pyrenean structure. Many transverse and longitudinal seismic profiles give strong constraints on the crustal structure, in particular on its thickness. The tomographic model of Souriau & Granet (1995) shows a low-velocity strip lying vertically under the Pyrenees, down to 80–100 km depth. It has been interpreted as the subduction of the Iberian lower crust under the Eurasian crust. To constrain the densities of these features, we use petrological arguments: lower crustal material subducting below 50–60 km must be metamorphosed in the eclogite facies, much denser than the surrounding mantle. The resulting 3-D density model gives a good fit to both observed Bouguer and geoid anomalies. Therefore, the crustal structure as imaged by seismic profiles is consistent with gravity data, and the presence of dense material at great depths explains the high positive geoid anomalies. This study shows that, despite the non-uniqueness of the gravity models, the joint interpretation of Bouguer and geoid anomalies may give important constraints on the deep structure of a mountain belt as soon as geometrical constraints are available.

A recent interpretation of magnetotelluric data proposed a different model, with lower crustal material plunging down to 80 km without being transformed. This 2-D model also fits Bouguer and geoid anomalies but is inconsistent with petrological arguments and seismic tomography, as untransformed crustal material would give a much larger velocity contrast than observed in seismic tomography. Partial melt can be an explanation for the observed high conductivities, but only at the deepest level (80 km depth). We propose that the transformation of lower crust into the eclogite facies, and the subsequent water release, can be another cause for the high conductivities observed by magnetotelluric measurements.

The metamorphism of lower crustal material is well documented and commonly invoked in geological studies. It seems to be of great help to build a consistent 3-D density model able to reconcile all the geophysical observations. Its consideration is promising for future models in any context of continental collision where subduction or large-scale crustal thickening appears.

Acknowledgments

We thank M. Daignières and M. de Saint-Blanquat for their comments on the manuscript.