This contribution reviews the challenges of imaging collisional orogens, focusing on the example of the Pyrenean domain. Indeed, important progresses have been accomplished regarding our understanding of the architecture of this mountain range over the last decades, thanks to the development of innovative passive imaging techniques, relying on a more thorough exploitation of the information in seismic signals, as well as new seismic acquisitions. New tomographic images provide evidence for continental subduction of Iberian crust beneath the western and central Pyrénées, but not beneath the eastern Pyrénées. Relics of a Cretaceous hyper-extended and segmented rift are found within the North Pyrenean Zone, where the imaged crust is thinner (10–25 km). This zone of thinned crust coincides with a band of positive Bouguer anomalies that is absent in the Eastern Pyrénées. Overall, the new tomographic images provide further support to the idea that the Pyrénées result from the inversion of hyperextended segmented rift systems.

Cette contribution passe en revue les enjeux de l’imagerie des orogènes collisionnels, en se focalisant sur l’exemple du domaine pyrénéen. Au cours des dernières décennies, notre compréhension de l’architecture de cette chaîne de montagnes a connu des avancées importantes grâce au développement de techniques innovantes d’imagerie passive, s’appuyant sur une exploitation plus poussée de l’information dans les signaux sismiques, ainsi que de nouvelles acquisitions sismiques. Les nouvelles images tomographiques mettent en évidence la subduction continentale de la croûte ibérique sous les Pyrénées occidentales et centrales, mais pas sous les Pyrénées orientales. Les reliques d’un rift Crétacé hyper-étendu et segmenté se trouvent dans la zone nord-pyrénéenne, où la croûte imagée est plus mince (10–25 km). Cette zone de croûte amincie coïncide avec une bande d’anomalies de Bouguer positives qui est absente dans les Pyrénées orientales. Les nouvelles images tomographiques viennent globalement renforcer l’idée que les Pyrénées résultent de l’inversion de systèmes de rifts segmentés et hyper-étirés.

Whereas surface geology of most mountain belts located in Western Europe has been mapped in detail, their internal structure remains elusive. In principle, seismic tomography can allow us to image the deep architecture of continental orogens, but this is a particularly challenging problem owing to the strong heterogeneity of the crust and the great difficulty of field acquisitions in these regions. To get new insights into the deep architecture of mountain belts, geophysicists thus need to improve and develop new passive imaging methods that can handle all the complexities of seismic wave propagation in such environments, but also to acquire new high quality data with a sampling density that is sufficient to resolve geological structures.

The Pyrénées are a favorable place to assess the potential of seismic tomography for unravelling the deep roots of orogens. From a geological point of view, this mountain belt is one of the best-studied collisional orogenic systems worldwide after more than a century of geological investigations. Its global architecture is also notably simpler than the nearby Alps because oceanic subduction, with all the potential implications, never occurred along this convergent margin. Indeed, it is now well accepted that its Late Cretaceous and Tertiary orogenic evolution results from the shortening of narrow Cretaceous hyper-extended rifts and not from an oceanic basin sensu stricto (see paper for a synthesis on the rift-to-orogenic evolution recorded in the Pyrénées and also the review papers of Manatschal et al., 2021, Jolivet et al., 2021 and Mouthereau et al., 2021 in this special issue). The moderate convergence (<200 km) (e.g., Muñoz, 1992; Teixell et al., 2018; Mouthereau et al., 2014) also implies that the records of the pre-collisional history remain largely preserved, which also makes the Pyrénées a unique place to study the control of rift inheritance on the formation of collisional orogens (see Manatschal et al., 2021, this special issue). Another key advantage of the study area is that this mountain range has been also a very active area of geophysical studies since the early eighties (Fig. 1). It was notably one of the very first to be crossed by a deep reflection seismic survey, the ECORS-Pyrénées profile (ECORS Pyrenees Team, 1988).

During the last decade, the PYROPE and IberArray deployments have covered the Pyrenean domain with a regular 2-D backbone of broadband stations, which was complemented by six dense transects deployed across the Pyrénées (Fig. 1). The Pyrénées thus emerge as a mountain range that has been exceptionally well covered with seismic instruments. These recent acquisitions also motivated important efforts to develop new passive imaging techniques. In this study, we will review these methodological developments, which allowed us to obtain finely resolved images of the deep architecture of the Pyrénées. The new tomographic images will be confronted to those from early deep seismic sounding studies, which will allow us to discuss the advantages and shortcomings of passive versus active acquisitions. We will also illustrate how these new tomographic images have deeply impacted our understanding of the formation of this mountain belt, emphasizing the control of rift-inherited structures during the phase of convergence.

The paper is organized as follows. We start in Section 2 with a review of the results obtained with the deep seismic sounding studies of the ECORS program. In Section 3, we describe the PYROPE and OROGEN temporary deployments, which considerably improved our capacity to perform passive imaging studies in the Pyrénées. Recent methodological developments that are particularly relevant for imaging mountain ranges are summarized in Section 4. Several key tomographic results are then presented in Section 5. Finally, we discuss in Section 6 the implications of these results for our understanding of the formation and evolution of the Pyrénées and give some perspectives on the future of passive imaging of the deep architecture of mountain belts.

Early deep seismic sounding studies documented a vertical Moho offset of about 15 km in the central Pyrénées (Daignières et al., 1982) and of about 5 km beneath the eastern Pyrénées (Gallart et al., 1980). The localization of the Moho offset, close to the limit between the Axial Zone (AZ) and the North Pyrenean Zone (NPZ), suggested that the North Pyrenean Fault (NPF) was a major structure, interpreted as the former plate boundary between Iberia and Europe. A few years later, the ECORS-Pyrénées profile (Fig. 1), the first seismic survey to cross an entire orogenic belt (McCaig, 1988; ECORS Pyrenees Team, 1988) revealed the fan-shape structure of the Pyrenean crust, with the NPF passively advected above north-vergent thrusts, and confirmed the existence of a thick crustal root made of the shortened Iberian crust beneath the AZ. Deep reflectors detected beneath the NPZ suggested that the Iberian crust is underthrusted beneath the southern edge of the European crust (Choukroune and The ECORS Team, 1989; Choukroune, 1992) but the extent and nature of the subducted material remained poorly constrained, leaving the grounds for various contradictory interpretations (Fig. 2). For example, Muñoz (1992) promoted continental subduction involving the Iberian lower crust with a nappe-stacking of the Iberian Upper crust, whereas according to Roure et al. (1989) the thick crustal root was produced by stacking fragments of lower crust, with limited northward extent of underthrusted material. These two models involve very different amount of convergence, with far-reaching consequences for the modes of formation and history of the Pyrénées.

The subsequent ECORS-Arzacq profile (Fig. 1), deployed in the western Pyrénées, provided little constraints on deep reflectors especially beneath the Mauléon Basin (Daignières et al., 1994). The crustal section proposed by Teixell (1998) is not compatible with the strong positive Bouguer anomaly measured in the western Pyrenean foothills. A more recent interpretation (Teixell et al., 2018) does incorporate an uplifted piece of lower crust or mantle rocks beneath the Mauléon Basin, but its size is incompatible with the magnitude of the positive Bouguer anomaly and with the strong (∼1 s) teleseismic P wave travel time difference between the rays that hit the surface inside the gravity anomaly and the rays outside (Chevrot, unpublished result). To summarize, whereas the ECORS program provided important insights into the crustal structure of the central Pyrénées, the deep architecture of this mountain range remained elusive.

Before 2010, the small number of permanent seismological stations installed in the Pyrénées was a limitation for structural studies (Fig. 3A). Stations were not only few but also poorly distributed, most of them being installed in the central elevated part of the range to monitor the local seismicity. In addition, at that time the data from the permanent IGN and IGC networks in Spain were not easily accessible. In fact, most of the available data came from the travel time picks in the bulletins from the different French and Spanish organisms that operated seismological networks in the Pyrénées. For these reasons, early regional tomography studies of the Pyrénées (e.g., Souriau and Granet, 1995; Souriau et al., 2008) should be considered with great care. During the last decade, great efforts have been made to acquire new data in the Pyrenean domain and to share in real time seismic data from all the operating networks (http://www.sispyr.eu), making the Pyrénées a mountain range among the most densely covered with seismic stations (Fig. 3B). As a result, this belt is thus nowadays one of the most favorable location to explore the potential of passive imaging approaches in a convergent margin setting.

The IberArray and PYROPE backbone deployments (2010–2013)

The seismological component of the TOPOIBERIA project (2006–2010), the IberArray seismic network (http://iberarray.ictja.csic.es/), covered the Iberian Peninsula in three successive deployments, from south to north. This initiative was a strong impulse to motivate the PYROPE project, a temporary deployment of broadband stations in the French side, synchronized with the beginning of the third IberArray deployment, north of Iberia. The main objective of the PYROPE project was to improve the seismic ray coverage in the Pyrenean domain by deploying a dense temporary array of seismological stations to collect high quality broadband data (Chevrot et al., 2014).

The Pyrenean transects

The 2-D backbone PYROPE array was complemented by two dense profiles that crossed the Pyrénées, approximately following the ECORS-Pyrénées and ECORS-Arzacq lines (Fig. 1). At the end of the PYROPE project, once the interest of such dense transects for studying deep lithospheric structures had been firmly established (Chevrot et al., 2015), additional deployments were planned within the framework of the OROGEN project to fully cover the mountain belt, from the Bay of Biscay to the Mediterranean Sea (Fig. 3B).

The large-N Maupasacq acquisition in the Mauléon Basin

From March to September 2017, a dense array of ∼450 sensors, the Maupasacq array, was deployed in the western Pyrénées foothills, partially overlapping the strong gravity anomaly of the Mauléon Basin (Polychronopoulou et al., 2018). The array was composed of about 200 10 Hz SERCEL nodes, 200 short-period sensors, and 50 broadband sensors (Fig. 4). The main purpose of this deployment was to explore the potential of passive approaches for imaging a basin in a foothill context, in a place where controlled-source experiments provided poor images of the structures in the basement (Daignières et al., 1994). Another original aspect of this acquisition is that it involved different types of sensors with geophone nodes, short-period, and broadband sensors, thereby allowing us to investigate their potential for the different types of passive imaging studies.

The recent passive seismic acquisitions have transformed the Pyrénées into an outstanding natural laboratory to experiment new tomographic approaches. The main challenge is to reach a kilometric spatial resolution, sufficient to characterize the main intracrustal geological units. This may at first sight appear as a rather modest objective. However, because passive approaches exploit seismic waves at frequencies typically below 1 Hz, this target resolution represents a fraction of seismic wavelengths. Such a resolution is well beyond the limits of classical asymptotic (ray based) tomography approaches that rely on the assumption that the medium is smooth at the scale of the wavelengths. Improving the resolution of tomographic models of continental orogens thus implies both a densification of seismic networks and a better exploitation of complete seismic wavefields, which raises important methodological as well as numerical challenges.

3-D modelling and inversion of teleseismic body waves

A first strategy to improve the resolution of passive imaging methods at the regional scale is to develop efficient hybrid numerical methods to model and invert short period teleseismic body wave fields, in the 1–20 s period range. The main idea is to inject an incident wave field computed in a spherical Earth model inside a 3-D regional spectral-element grid. In that case, the time consuming 3-D computations are limited to the regional grid and accurate synthetic seismograms, which account for free surface topography, 3-D heterogeneities (in VP, VS, density, anisotropy, and even attenuation) can be obtained with reasonable computational resources. The principles of the injection method are summarized in Figure 5 (Monteiller et al., 2021).

A first step came from the development of a method that coupled the Direct Solution Method (DSM) (Kawai et al., 2006) for the global propagation with the spectral-element method (SPECFEM3D) (Komatitsch and Tromp, 1999) in the regional domain (Monteiller et al., 2013). Whereas this early implementation was important to firmly establish the potential of hybrid methods to compute synthetic seismograms in complex 3-D regional models, DSM is a cumbersome method to use for the computation of the incident wavefield because it involves storing and manipulating a large number of spherical harmonic coefficients, a problem that becomes particularly severe for shallow sources. This motivated the development of a new injection technique relying on AxiSEM (Nissen-Meyer et al., 2007, 2014), which was first used in a tomographic study of the western Alps (Beller et al., 2018). A final development has been to also implement the injection of an incident plane wave computed with the frequency-wavenumber (FK) method (Monteiller et al., 2021). When the regional domain is small (less than several hundreds of kilometres), the effects of the curvature of the wavefront and of the Earth can be neglected, and FK provides an acceptable approximation of the incident wavefield. All these numerical tools are now gathered in RegHyM, a complete package to compute synthetic seismograms in 3-D regional models that is now available and distributed to the scientific community (https://gitlab.com/Seismic_Imaging/RegHyM).

Eikonal/Helmholtz surface wave imaging

Surface waves represent another important source of information to constrain lithospheric structures. Using records from dense local or regional seismic arrays, phase velocities can be directly mapped from the amplitude and phase of surface-wave records with the Helmholtz equation. Because this approach does not require the resolution of an inverse problem, it avoids degrading the resolution by introducing arbitrary damping and smoothing regularization constraints into the inversion. Following the booming of dense regional deployments, Helmholtz tomography rapidly gained in popularity and was extensively used for example to image the shear velocity structure of the crust and uppermost mantle beneath North America using USArray Transportable Array (TA) data (e.g., Pollitz, 2008; Lin et al., 2009; Jin and Gaherty, 2015; Shen and Ritzwoller, 2016). In spite of these remarkable successes, Helmholtz tomography has so far received little interest for imaging the lithospheric roots of continental orogens. One can thus reasonably expect major advances coming from this approach in the near future.

The main challenge in Helmholtz tomography is the estimation of the Laplacian of the amplitude field, which can be problematic from noisy and sparsely sampled seismic data. Neglecting the amplitude term in the Helmholtz equation leads to the simpler eikonal equation, which only involves the gradient of the phase field.

In a seminal study, Wielandt (1993) argued that using the eikonal equation to derive phase velocities results in models contaminated by strong artefacts. However, Lehujeur and Chevrot (2020b) derived analytical expressions of the errors arising from neglecting the amplitude in eikonal tomography. They showed that, in general, these errors are quite strong but vary sinusoidally with the incoming wave direction. Therefore, provided good azimuthal coverage, these errors will tend to cancel out, and eikonal tomography will give phase velocity maps with a spatial resolution on par with Helmholtz tomography or even full waveform inversion (Lehujeur and Chevrot, 2020b).

Extraction of coherent surface wave fronts from ambient seismic noise

Following the pioneering work of Shapiro and Campillo (2004) and Shapiro et al. (2005), seismic noise correlations have been extensively used to derive shear waves velocity models of the Earth’s crust in various regions such as the Alps and the Pyrénées (Yang et al., 2007; Stehly et al., 2009; Palomeras et al., 2017; Kästle et al., 2018; Lu et al., 2018, 2020; Zhao et al., 2020; Macquet et al., 2014). The idea is to exploit the permanent excitations of the Earth by the ocean waves, which generate seismic waves that are recorded on distant continents (e.g., Ardhuin et al., 2011). By correlating the seismic ambient noise recorded at a pair of stations, it is then possible to reconstruct the surface wave signal that would have been produced by a source located at one station and recorded at the other (e.g., Shapiro and Campillo, 2004), often referred to as “the empirical Green’s function”. After measuring the dispersion of surface waves between all the pairs of stations belonging to a regional array, one can then invert the resulting measurements to obtain a 3-D shear wave velocity model of the lithosphere.

Array analysis of ambient seismic noise has shown that it is mainly composed of coherent surface wave fronts produced by the action of swells and oceanic waves on the sea floor (e.g., Chevrot et al., 2007). Lehujeur and Chevrot (2020a) developed a matched-filtering method to isolate and extract these coherent surface-wave fronts and applied it to the ambient seismic noise recorded by the Maupasacq array. Figure 6 shows examples of isolated Rayleigh wave fronts at 5.2 s period coming from two different azimuths. These surface waves can then be exploited for seismic tomography in the same manner as those emitted by large earthquakes. The 3-D model derived from these dispersion measurements is presented in Lehujeur et al. (2021), this issue. The method opens new perspectives for applying Eikonal/Helmholtz tomography on surface waves extracted from ambient seismic noise at periods lower than 20 s, for which it is extremely difficult to get reliable measurements from surface waves produced by earthquakes, owing to their very long propagation through heterogeneous crusts.

Improving ray coverage in seismic ambient noise tomography with higher-order correlations

To improve the ray coverage and thus potentially the spatial resolution of ambient noise tomography, the possibility of exploiting higher-order noise correlations has been investigated (Brives, 2020). The main idea is to turn a network of seismic stations into a set of virtual sources that can be used to measure the velocity of Rayleigh waves between any pair of stations, even if these stations were not operational simultaneously (Stehly et al., 2008; Ma and Beroza, 2012; Brives, 2020). This method, called C2 for “Correlations of seismic noise Correlations” has several benefits. First, it allows us to measure the phase velocity of Rayleigh waves between the stations of the two PYROPE transects, which were deployed during different period of times, thereby considerably improving the ray coverage in the internal Pyrenean domain (Fig. 7b). Second, by improving the quality and robustness of reconstructed Green’s functions, it is possible to increase dramatically the number of paths that provide robust dispersion measurements exploitable in a tomographic inversion. With the C2 method, between 6,000 and 14,000 Rayleigh wave velocity measurements were collected, depending on the period range considered, i.e. more than twice as many measurements compared to the initial study performed by Macquet et al. (2014). These measurements have been used to compute 2-D group velocity maps at period ranging from 3 s to 70 s (Fig. 7c). As shown on Figure 7c, the spatial resolution of C2 method group velocity map is improved, in particular in the central part of the Pyrénées, and is characterized by a much denser ray coverage in the Bay of Biscay, which was poorly constrained with the C1 method.

This study demonstrates that in addition to new acquisitions, improved analysis techniques can also lead to a significant improvement of the number and quality of dispersion measurements, which is crucial to obtain crustal images with a finer resolution.

Improving ray coverage in local earthquake tomography with deep learning approaches

Another important source of information comes from the regional Pyrenean seismicity (Sylvander et al., 2021) which produces P and S waves that are recorded by nearby seismic stations. Whereas these waves are routinely picked on permanent stations to build catalogues of seismicity, records from temporary deployments are rarely exploited and most of the time absent of seismic bulletins. The manual picking of P and S waves is indeed very time consuming, even for a seasoned analyst. However, recent studies have found that deep neural networks can be trained to perform this routine analysis, with results that are on par with human picks (e.g., Ross et al., 2018; Zhu and Beroza, 2019). This opens important new perspectives for building complete and extensive catalogues of P and S travel time picks from the exploitation of PYROPE and IberArray datasets that could be exploited for local earthquake tomography. This work is currently underway.

3-D imaging of seismic anisotropy

Seismic anisotropy is a powerful tool for characterizing and mapping fossil lithospheric deformations that accumulate over geological history (Silver, 1996).

Perhaps the most simple and unambiguous manifestation of seismic anisotropy is shear wave splitting, which splits shear waves into two orthogonally polarized quasi-shear waves propagating with different velocities (Savage, 1999). Shear wave splitting is quantified by the splitting delay accumulated between the two quasi-shear waves and the fast polarization directions, which follow mantle fabrics and the preferential orientation of olivine crystals.

Seismic anisotropy studies can therefore potentially provide first-order constraints on mountain building processes and is often used as a valuable indicator for plate kinematics (examplified by Jolivet et al., 2021 in this special issue), but they still suffer from a lack of resolution, which prevents us from characterizing the relative contributions of the lithosphere and of the asthenosphere. To address this important problem, the teleseismic full waveform inversion method has been generalized to image seismic anisotropy (Beller and Chevrot, 2020). Moving from isotropic to anisotropic inversions is a major challenge because a general 4th-order elasticity tensor is described by 21 elastic coefficients. To make the problem tractable, some simplifying assumptions are usually made on the symmetry class of the medium and on the orientation of the planes and axes of symmetry. For example, anisotropic media are often described in terms of Horizontal Transverse Isotropy (HTI), an hexagonal media with a symmetry axis lying in the horizontal plane. Whereas hexagonal symmetry is in general a good approximation of mantle anisotropy (e.g., Becker et al., 2006), the symmetry axis may not necessarily be horizontal. The main idea of this new approach is to invert for the 21 elastic coefficients, in addition to density, and to perform a post-inversion analysis of the resulting elasticity tensors at each node of the tomographic grid. This analysis involves first determining the principal axes of the elasticity tensor and then to decompose the elasticity tensor rotated into its intrinsic system of coordinates, the system of coordinates described by the main elements of symmetry. For example, in the case of hexagonal symmetry, the symmetry axis will correspond to one of the axis of the reference frame, and the two other orthogonal axes can be chosen arbitrarily in the plane perpendicular to the symmetry axis. This approach thus allows us to image any type of anisotropic medium without making any a priori assumption on the nature of seismic anisotropy. Numerical synthetic tests have shown that with this method it is not only possible to constrain the dip of the symmetry axis, but also to resolve and discriminate seismic anisotropy in the lithosphere and asthenosphere (Beller and Chevrot, 2020).

3-D Modelling of regional gravity anomalies

The nature and geometry of the dense bodies that produce the positive Bouguer gravity anomalies in the Pyrénées has been a long-standing controversy (Torné et al., 1989; Grandjean, 1994; Casas et al., 1997; Vacher and Souriau, 2001; Pedreira et al., 2007; Jammes et al., 2009). Early interpretation of the positive Bouguer anomaly observed on the ECORS-Pyrénées profile invoked a slice of either mantle or lower crustal material emplaced in the upper crust beneath the NPZ (Torné et al., 1989). A later study ascribed the Bouguer anomalies in the western and central Pyrénées to slices of upper mantle exhumed during the transtensional motion of Iberia (Casas et al., 1997). This study also emphasized that the positive Bouguer anomalies are observed in the central and western Pyrénées, but not in the eastern Pyrénées, which was interpreted as indication of a smaller amount of Cretaceous extension in the east before the convergence. With the exception of Vacher and Souriau (2001) and Pedreira et al. (2007), these early studies relied on 2-D modelling, which has the detrimental effect of overpredicting the magnitude of density anomalies.

Recent efforts focused on developing a new spectral-element method to compute and invert gravity anomalies at the regional scale (Martin et al., 2017). The computation of the volumetric integration of density models with spectral-element quadratures is both more precise and much faster than the classical semi-analytical approaches. In addition, by imposing surface topography on the upper edge of the 3-D spectral element mesh, topographic or terrain corrections are implicitly taken into account in the computations of gravity anomalies. Because the computations are performed on the same grid as the one used in the wave propagation problem, this method opens the way to the joint inversion of seismic and gravity data (Martin et al., 2021; Wang et al., 2022). Precise Bouguer corrections can be computed by integrating the contribution of masses between the sea level and surface topography with the spectral-element method (Martin et al., 2017). The map of Bouguer gravity anomalies derived with this approach using ETOPO1 surface topography and a crustal density of 2.67 is shown in Figure 8. The deep crustal roots beneath the Axial Zone in the central Pyrénées produce a broad negative Bouguer gravity anomaly. The other salient features are the two strong positive Bouguer anomalies beneath the Mauléon Basin (MB) and the Saint-Gaudens region (SG).

The RF sections

The migration of receiver functions along the six PYROPE and OROGEN transects (Fig. 1) provided the first global vision of the deep crustal architecture of the whole Pyrénées. Figure 9 shows the compilation of migrated sections for the 5 N–S Pyrenean transects, sorted from west to east (Chevrot et al., 2018).

The westernmost transect (A–A’) shows the continental subduction of the Iberian crust and the presence of a shallow strong discontinuity beneath the Mauléon Basin at around 10 km depth. Further east, the second profile (B–B’) shows very similar structures, with a north-dipping continental subduction and the presence of a discontinuity at about 20 km depth beneath the NPZ. Beneath the central transect (C–C’), which approximately follows the ECORS-Pyrénées profile, the north-dipping subduction of the Iberian plate is still observed, with a Moho located around 20 km depth. The first three profiles thus all suggest a thinned European crust beneath the NPZ. The reverberations in the shallow sedimentary layers of the Aquitaine Basin are responsible for the characteristic stripe patterns observed in the northern parts of transects A–C. Further east, the D–D’ section shows a completely different crustal architecture. The northward underthrusting of the Iberian plate is no longer observed and the orogenic crust is laterally segmented into three distinct units, delimited by small but clear Moho offsets. In the central unit, the polarity of crustal thickening is reversed compared to the western sections, with a Moho that is tilted towards the south.

The easternmost transect (E–E’), close to the Mediterranean coastline, is also segmented into three crustal domains and characterized by the absence of underthrusted Iberian crust. Crustal thickness is close to normal along the whole profile, with little Moho topography. The absence of a crustal root in the eastern Pyrénées suggests that the present-day topography cannot be explained by isostasy alone and requires some form of dynamic support in this area. Finally, the E–W section (F–F’) shown in Figure 10 samples the transition from the thickened Pyrenean domain to the Mediterranean Sea and shows a progressive eastward thinning of the crust (see also Díaz et al., 2018).

The first applications of teleseismic P wave FWI

Numerical hybrid methods opened (Monteiller et al., 2013) the way to the full waveform inversion of teleseismic P waves. The Gauss-Newton L-BFGS inversion algorithm was first tested in synthetic tomographic experiments (Monteiller et al., 2015) and subsequently applied for the very first time on the data recorded by the westernmost Pyrenean transect (Wang et al., 2016). In that study, vertical and radial component P wave records from 5 teleseismic sources were inverted to determine models of both compressional (VP) and shear (VS) wave velocities (Fig. 11). The VS model is visually characterized by a finer resolution, a direct consequence of the shorter wavelength of shear waves compared to compressional waves (the shear waves are produced by P-to-S conversions and therefore have the same frequency content as the P waves). Overall, the VS model shows an excellent agreement with the receiver function migrated section (Fig. 9, section A–A’), in particular for the geometry of the imbricated Iberian and European Moho. However, FWI tomography provides images that are more prone to geological interpretation. For example, the combination of VP and VS provides tighter constraints on the nature of crustal rocks (e.g., Watanabe, 1993; Christensen, 1996; Brocher, 2005). Constraining density (see Sect. 5.7) should allow us to further constrain the composition of deep crustal materials. The low velocity anomaly dipping north beneath the European mantle in both models provides further evidence for the underthrusting of continental material (i.e. “continental subduction”) beneath the western Pyrénées. Scaling the VP model to density using a standard linear Birch-type relationship shows that this model predicts the shape and amplitude of the positive Bouguer anomaly in the Mauléon Basin (Wang et al., 2016). Because the fast anomaly beneath the Mauléon Basin (see also Lehujeur et al., 2021, this special issue) is connected to the European mantle, it probably corresponds to a mantle body that has been formerly exhumed during the Cretaceous episode of rifting. Therefore, the tomographic model suggests that we image a piece of preserved hyper-extended rift, which demonstrates the important role of rift-inherited structures in the formation and evolution of the Pyrénées.

Figure 12 shows the model obtained with full waveform inversion (Wang, 2017) on the central profile (C–C’) in Fig. 1), which approximately followed the ECORS-Pyrénées transect (black line in Fig. 1). The most salient features are the north-dipping low-velocity anomaly beneath the southern border of the European plate, which corresponds to underthrusted Iberian crust, and the thinner European crust beneath the North Pyrenean Zone. This image is in good agreement with both the RF section (section C–C’, Fig. 9) and the ECORS-Pyrénées section (Roure et al., 1989).

Ambient noise tomography of the Pyrénées with the C2 method

Using one year of noise recorded by permanent and temporary PYROPE stations, Macquet et al. (2014) obtained the first tomographic model of the Pyrenean domain with ambient noise tomography. The model successfully imaged the geometry of sedimentary basins and the Moho, the thicker Iberian crust and the thinner European crust, and the structural dissymmetry between the South Pyrenean Zone and the North Pyrenean Zone at shallow crustal levels. However, only a fraction of station pairs provided exploitable dispersion curves, owing to the uneven azimuthal distribution of noise sources, which are mainly located in the Atlantic Ocean close to Galicia in the West and in the Mediterranean Sea in the East (Chevrot et al., 2007). In addition, the data from the Pyrenean transects were not available at the time of this early study.

Brives (2020) used the group velocity measurements shown in Figure 7 to build a probabilistic 3-D VS model of the Pyrénées with a Bayesian inversion method (Lu et al., 2018). The final VS model obtained provides at each pixel the probability distribution of VS and the probability of the presence of a seismic discontinuity (e.g. the Moho). It was built assuming that at each location the medium can be described by a 4-layer lithospheric template (sediment, upper crust, lower crust, upper mantle). Figure 13 shows the VS model at 4 km and 35 km depth. At 4 km depth, low velocity anomalies correspond to thick sedimentary basins such as the Aquitaine Basin and the Bay of Biscay. The 35 km depth slice (Fig. 13, right) highlights variations in crustal thickness, with low velocities (VS < 3.7 km/s) in the central Pyrénées and the Armorican and Iberian Massifs, where the crust is particularly thick. High velocity anomalies (VS > 4 km/s) are found in the Bay of Biscay and in the Gulf of Lion, where the mantle is shallower.

Figure 14 presents two cross-sections across the 3-D VS model along the I–I’ and J–J’ transect shown in Figure 13. For each cross-section, the upper panel shows the most probable shear wave velocity model and the lower panel the probability of having an interface at each depth. In cross-section I–I’, which is located close to the receiver function section B–B’ (Fig. 9), we observe that the Moho depth increases as we move toward the north up to the North Pyrenean Fault. By contrast, beneath the North Pyrenean Fault, we observe fast S-velocities of the order of 3.8 km.s−1 at 15 km depth, and the probability of having an interface is maximum at 20 km depth, which indicates a shallow Moho depth. This is consistent with the observations made on the receiver functions in this area. As we move further east along the J–J’ cross-section, the Moho shows little topography. Below the Ebro Basin, we observe a flat Moho at a depth of 30 km, and a slightly shallower Moho below the Aquitain Basin. Similarly, the VS model does not show clearly an increase of the Moho depth below the axial zone. This is consistent with the absence of crustal root shown by the receiver functions of the eastern profiles (Fig. 9). For more detailed comparisons of crustal thickness models derived from ambient noise tomography, the reader is referred to Brives (2020).

These results highlight the benefits of using iterative seismic noise correlation to image crustal structure over a large-scale region.

Eikonal tomography

The iterative extraction of coherent waves from Maupasacq data identified a variety of noise sources in a broad range of azimuths that are usually buried beneath the predominant Atlantic and Mediterranean noise sources (Lehujeur and Chevrot, 2020a). In their models, phase velocities in the Arzacq and Mauléon basins nicely correlate with surface geology at periods shorter than 3 s. These dispersion measurements were later exploited to constrain the 3-D VS structure of the Mauléon and Arzacq basins (see Lehujeur et al., 2021, this special issue) and address the question of mantle sampling within orogens. Note that this latter question is far from being restricted to the Pyrénées (e.g. Ivrea body in the Italian Alps or Ronda peridotites in the Western Betics as discussed in Bessière et al., 2021, this special issue).

Local earthquake tomography

The Pyrénées are characterized by a moderate seismic activity that can potentially be exploited for tomographic imaging. The only local earthquake tomography that covers the whole Pyrénées dates back to the 1990s (Souriau and Granet, 1995). This early study imaged in the upper crust the fast velocities in the basement rocks of the axial zone, the low velocities in the Aquitaine and Ebro basins, and a band of high velocity anomaly beneath the NPZ interpreted as evidence for uplifted blocks of upper mantle material. At a deeper level, it also imaged the thick crustal root beneath the Axial Zone. Since then, the permanent networks on both sides of the Pyrénées have been greatly improved both in terms of number and quality of sensors. There is thus a wealth of information basically untapped in the recent catalogues of seismicity produced by the different organisms that operate seismological networks in the Pyrénées.

Here we present new but preliminary results of local earthquake tomography in the Pyrénées: first a local-scale study for the region covered by the large-N deployment described in Section 3.3, and second a more regional study of the western Pyrénées that includes the local one.

The data used to obtain the 3-D P- and S-wave velocity models of the local-scale model (Mauléon area) consist of arrival times of P and S waves from local earthquakes recorded inside the large-N experiment. The P-wave picks where obtained with the approach described in Tselentis et al. (2011), consisting in filtering the record in the S-transform domain for denoising purposes, and then applying the Kurtosis criterion. The S waves arrival times, which are often more difficult to pick, were obtained following Lois et al. (2013). After picking all the stations of the network and associating the arrival time picks that were compatible with local events, a catalog of 1,624 earthquakes was obtained. From this catalog, we selected well recorded earthquakes for the tomographic inversion that met the following criteria: number of P arrival times ≥8, number of S arrival times ≥2, and azimuthal gap ≤200°. This resulted in a dataset of 996 earthquakes with 87,122 P arrival times, and 72,445 S arrival times. The size of the model region was a rectangle of 130 km × 115 km in the E–W and the N–S directions, respectively (see Fig. 15).

This dataset produces a high-resolution model inside the region covered by the dense station deployment (e.g., Fig. 15a), but the ray-path sampling falls off quickly. To increase the data coverage both laterally and in depth, we also considered a model for a larger region, using arrival time data from the large-N experiment and from a regional dataset of earthquakes in the western Pyrénées recorded by permanent stations of the monitoring networks and also by the temporary stations described in Sections 3.1 and 3.2. Using the same selection criteria, we obtained a dataset consisting of 4,175 earthquakes recorded on 553 stations, for a total of 103,314 P arrival times, and 94,512 S arrival times. The size of the regional model is a rectangle of 255 km × 220 km in the E–W and the N–S directory respectively (between 2.5°W to 0.5°E longitude and 42°N to 44°N in latitude, see Fig. 16a).

For the tomographic inversion we used the method of Benz et al. (1996) as modified by Tryggvason et al. (2002) in order to include also S arrival times. This method performs an iterative joint inversion for velocity structure and earthquake relocations.

Images of the resulting local and regional models are shown in Figures 15 and 16, respectively. The local model images in detail the shallow structure of the region, with main structures trending WNW–ESE (Fig. 15a). The southern boundary of the Arzacq basin is also clearly defined. At larger depths (>8 km) the regional model shows a high velocity anomaly coincident in location with the strong positive Labourd-Mauléon Bouguer gravity anomaly. The geometry of this fast anomaly is in good agreement with the results of full waveform inversion shown in Figure 11. The eastern limit of this high velocity body is sharply defined at 0.8°W (Fig. 15b). Because of the relatively small size of the local model, there are few rays that cross the parts of the model deeper than the deepest seismicity (∼15 km). On the other hand, the regional model does not have the resolution of the local one, but can image larger depths, and also a larger area. This is illustrated in Figure 16, which shows that structures down to 35 km depth are still resolved.

Seismic anisotropy

Previous analyses of SKS splitting in the Pyrenean domain revealed predominant fast polarizations oriented parallel to the strike of the belt (Barruol and Souriau, 1995; Barruol et al., 1998). This orientation has been variably ascribed to the translation of Iberia with respect to Europe (alternatively Permian, Cretaceous or both), to the pure shear deformation caused by the N/S Pyrenean collision, or to the mantle counter-flow caused by the Alpine slab retreat and the Gulf of Lion back-arc rifting initiated in the Oligocene (∼35 Ma) (Barruol et al., 2004; Jolivet et al., 2021). From the exploitation of IberArray data, Díaz et al. (2015) found similar fast directions throughout the whole Iberian Peninsula, which seems to favor large-scale asthenospheric flow rather than a Pyrenean-related process. The recent analysis of shear wave splitting on the stations from the permanent French and Spanish networks and from the PYROPE and IberArray temporary experiments improved considerably the spatial coverage of splitting measurements in the Pyrenean domain (Bonnin et al., 2017). This analysis confirmed the large-scale coherent pattern of seismic anisotropy in that region, but also documented abrupt changes of splitting parameters, along the Pyrenean transects, in particular between the Pyrénées and the Aquitaine basin (Fig. 17). These new data suggest that both a spatially coherent asthenospheric flow related to the opening of the western Mediterranean and fossilized lithospheric fabrics are required to explain the apparent shear wave splitting patterns.

Density structure of the lithosphere

From a compilation of geological and geophysical information a 3-D model of crustal density was built with the Geomodeller software from the BRGM (Wehr et al., 2018). The final model fits the observed Bouguer anomalies within 3.5 mGal. Figure 18 shows the Bouguer anomalies corrected for the contribution of the variations of crustal thickness and of sediments. Strong positive residual anomalies are observed beneath the North Pyrenean Zone in the western and central Pyrénées. In the eastern Pyrénées, where the crustal thickness is around 30 km, a broad and/or deep low-density anomaly is required to explain gravity data, suggesting that this part of the range is underlaid with a hot and buoyant lithosphere that provides some dynamic support to topography. Overall, the gravity anomalies are in rather good agreement with the results of tomographic studies.

Going one step further, focusing on the western Pyrénées, Martin et al. (2021) have demonstrated that we can obtain 3-D density models from gravity inversions constrained by tomographic models. The initial density model is derived by scaling the VP model obtained by full waveform inversion of teleseismic waves (Wang et al., 2016). The computations of Bouguer gravity anomalies are performed with the 3-D spectral-element method described in Martin et al. (2017). The inversion is regularized by enforcing a structural similarity between the VP and density models by minimizing the norm of their cross-gradient. The final density model (Fig. 19) provides additional constraints on the thermal state and composition of the lithosphere in the western Pyrénées.

The recent progresses on our understanding of the deep architecture of the Pyrénées can be ascribed to the development of innovative passive imaging approaches as well as the new dense high-quality seismic acquisitions. FWI of teleseismic body waves has been applied for the very first time in the Pyrénées, firmly establishing it as a sensitive approach to image continental orogens at the lithospheric scale, with results in agreement with those of deep seismic reflection acquisitions but for a fraction of their cost. FWI leads to 3-D models of seismic velocities (and of density) and not only of seismic interfaces as with receiver functions. These models thus carry additional information that is particularly important to constrain the nature and composition of crustal rocks.

We have identified two important sources of information that remain largely untapped. First, surface waves should allow us to improve dramatically the resolution of tomographic models in the lithosphere and in particular to constrain the depth of the lithosphere/asthenosphere boundary. The eikonal/Helmholtz tomography approach, which is rapidly developing, should play a pivotal role to better exploit surface waves recorded by dense regional arrays such as PYROPE and IberArray for the Pyrénées or AlpArray for the Alps.

In addition, from the azimuthal anisotropy of surface waves, its depth distribution, and the waveforms of teleseismic P and S waves, we are now opening the doors for imaging the 3-D distribution of seismic anisotropy in the near future. Since the distribution of lithospheric deformations are mainly controlled by the mechanical behavior of the lithosphere, this should provide key insights into mountain building processes in particular and more generally into plate kinematics. Second, the local and regional seismicity recorded by temporary stations should allow us to improve dramatically the catalogues of P and S picks for local earthquake tomography. Since picking manually all the P and S arrivals on these datasets would be extremely time consuming, we will need to turn to semi-supervised approaches relying on deep neural networks. With these new complete catalogs, we should be able to improve significantly the quality of the images at upper and mid crustal levels.

Additionally, the Maupasacq deployment allowed us to demonstrate the potential of passive approaches to image the deep architecture of sedimentary basins (see Lehujeur et al., 2021, this special issue). For example, the velocity lineaments seen with passive methods (Fig. 15a) are also in very good agreement with the kilometric-scale folds in the Mauléon Basin and with the major thrusts such as the North Pyrenean Front (Lehujeur et al., 2021). Such structural details are getting close to the usual targets of industry in foothills domains.

Finally, to further improve our tomographic images of continental orogens, we will need to develop approaches that will take advantage of the complementary of various types of data such as surface-wave and body-wave seismic records, or seismic and gravity data. This will imply developing a new generation of inversion codes, with efficient interfacing to perform large-scale 3-D inversions of tomographic models defined with a much larger number and types of parameters, such as the 21 elastic coefficients and the density. There is thus still a long way to go, but at least we now have a much clearer vision of the directions to follow.

Thanks to a decade of new seismic acquisitions and geophysical studies in the Pyrénées from which the OROGEN project took a significant role, important progresses have been achieved in our understanding of the deep architecture of this mountain range. The recent tomographic results challenge the view of a Pyrenean range shaped by the transtensional motion of Iberia with respect to Europe (e.g., Roure et al., 1989; Choukroune, 1992). In contrast, they strongly support the idea that the Pyrénées result from the inversion of a former hyper-extended rift system. They also reveal a strong asymmetry between the eastern and western-central Pyrénées, with a continental subduction of the Iberian plate associated to shallow exhumed mantle beneath the northern Pyrenean foothills observed in the Western and Central Pyrénées but absent in the eastern Pyrénées. They also provide a plausible response on the source of the large positive Bouguer anomaly measured in the NPZ, in the western and central Pyrénées. We suggest that these anomalies result from the crustal thinning and mantle exhumation that occurred during the Cretaceous episode of rifting, and that they have been locally preserved from orogenic overprints. The different crustal architecture between the eastern and western-central Pyrénées is most likely inherited from the Cretaceous geometry of the rift system, which controlled the distribution of deformation when the system was inverted and shortened (see discussion in Chevrot et al., 2018), even though at this point the influence of the opening of the western Mediterranean basin (Jolivet et al., 20202021) cannot be ruled out.

We thank Emanuel Kästle, Gary Pavlis, and the associate editor David Pedreira for their insightful and constructive reviews. This work is part of the OROGEN research project, a tripartite partnership between the CNRS, TOTAL and BRGM.

Cite this article as: Chevrot S, Sylvander M, Villaseñor A, Díaz J, Stehly L, Boué P, Monteiller V, Martin R, Lehujeur M, Beller S, Brives J, Bitri A, Calassou S, Collin M, Ford M, Jolivet L, Manatschal G, Masini E, Mouthereau F, Vidal O. 2022. Passive imaging of collisional orogens: a review of a decade of geophysical studies in the Pyrénées, BSGF - Earth Sciences Bulletin 193: 1.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.