We present a 3-D shear wave velocity model of the Mauléon and Arzacq Basins from the surface down to 10 km depth, inverted from phase velocity maps at periods between 2 and 9 s. These phase velocity maps were obtained by analyzing coherent surface wave fronts extracted from ambient seismic noise recorded by the large-N Maupasacq seismic array with a matched filtering approach. This new model is in good agreement with a local earthquake tomography study performed on the same acquisition dataset. Our passive imaging models reveal the upper crustal architecture of the Mauléon and Arzacq Basins, with new details on the basement and its relationship with the overlying sedimentary cover. Combining these new tomographic images with surface and subsurface geological information allows us to trace major orogenic structures from the surface down to the basement. In the basin, the models image the first-order basin architecture with a kilometric resolution. At depth, high velocity anomalies suggest the presence of dense deep crustal and mantle rocks in the hanging wall of north-vergent Pyrenean Thrusts. These high velocity anomalies spatially coincide with a positive gravity anomaly in the western Mauléon Basin. In addition, our models reveal major changes from the Chaînons Béarnais to the western Mauléon Basin across a set of orogen-perpendicular structures, the Saison and the Barlanès transfer zones. These changes reflect the along-strike variation of the orogenic evolution that led to the preservation of the former rifted domain and its underlying mantle in the orogenic wedge of the Western Pyrenees. We discuss the implications of these results for the 3-D architecture of the Mauléon Basin and its underlying basement.

Nous présentons un modèle 3-D de vitesse des ondes de cisaillement des bassins de Mauléon et d’Arzacq de la surface jusqu’à 10 km de profondeur inversé à partir de cartes de vitesse de phase pour des périodes entre 2 et 9 s. Ces cartes ont été obtenues à partir de l’analyse de fronts d’onde de surface cohérents extraits du bruit sismique ambiant enregistré par le réseau Maupasacq par filtrage adaptatif. Ce nouveau modèle est en bon accord avec la tomographie locale réalisée sur ce même jeu de données. Nos nouvelles images tomographiques révèlent l’architecture supra-crustale des bassins de Mauléon et d’Arzacq, avec des informations nouvelles sur la nature du socle et sa relation à la couverture sédimentaire. En combinant ces nouvelles images tomographiques aux informations géologiques, il est possible de tracer les principales structures orogéniques de la surface jusqu’au socle des bassins. Dans le bassin, les modèles nous fournissent une image de premier ordre des plis et chevauchements à l’échelle kilométrique. En profondeur, les anomalies rapides suggèrent la présence de roches de la croûte inférieure et du manteau dans le toit des chevauchements pyrénéens de pendage nord. Ces anomalies rapides coïncident spatialement avec l’anomalie gravimétrique positive dans la partie ouest du bassin de Mauléon. Nos modèles tomographiques documentent également des changements de structures majeurs entre les Chaînons Béarnais et la partie ouest du bassin de Mauléon à travers des structures perpendiculaires à l’axe de la chaîne, représentées par les structures transverses du Saison et du Barlanès. Ce changement structural reflète les variations latérales de l’évolution orogénique qui a conduit à la préservation des domaines de rift hyper-étirés et du manteau sous-jacent dans le prisme orogénique. Nous discutons les implications de ces résultats concernant l’architecture 3-D du bassin de Mauléon et du socle sous-jacent.

Imaging crustal structures with a fine spatial resolution is an important goal of modern seismology, with major implications in domains such as georesources and seismic hazard. In regions where strong crustal heterogeneity prevails, such as in orogens, passive imaging using local earthquakes remains challenging owing to the uneven distribution of seismicity, both in time and space, and to the sparsity of permanent seismic networks. With ambient seismic noise, seismologists can now reconstruct surface waves from pairs of seismic stations, thereby freeing themselves from relying on the occurrence of natural earthquakes (e.g., Shapiro et al., 2005). Ambient noise surface wave tomography (SWT) has gained a large interest in the academic world and have been applied to a very broad range of scales and contexts. Most studies focused on large scale applications like North America (e.g., Lin et al., 2008) or Europe (e.g., Lu et al., 2018). The extension of the method at a regional scale using lower periods (below 10 s) are often very efficient to image the contours of the sedimentary basins, usually characterized by lower velocities (e.g., Macquet et al., 2014). Fewer studies have focused on the use of surface waves from the ambient seismic noise to image structures inside a sedimentary basin (Jia and Clayton, 2021 being a recent exception). The advantages related to this passive source of seismic waves i.e., continuity of sources and reduced acquisition costs compared to active seismic methods, did however motivate focused applications on sedimentary basins for example for the characterization of an oil and gas field (e.g., Mordret et al., 2013), the subsurface imaging for the exploration of deep geothermal resources (e.g., Lehujeur et al., 2018; Planès et al., 2019) or the monitoring of CO2 underground storage sites (e.g., Gassenmeier et al., 2014). Another recent advance came from the recognition that large-N node deployments such as those commonly used in controlled-source acquisitions for the oil and gas industry provide rich and valuable datasets for passive imaging studies (e.g., Schmandt and Clayton, 2013). These two recent developments have opened important new perspectives not only for crustal-scale tomography, but also to get valuable insights on the 3-D geometry of sedimentary basins in structurally complex areas to complement active seismic reflection surveys.

To investigate the potential of passive imaging for the characterization of sedimentary basins, which has so far received little attention from the academic world, we have deployed a large-N array of 442 sensors, the Maupasacq experiment, in the western part of the Northern Pyrenees foothills, from March to September 2017. The inner part of the network formed a regular 50 × 30 km rectangular grid, extended to about 120 × 130 km, thanks to two additional circular deployments (Fig. 1). In the inner rectangular grid, the station spacing ranged from 1 km between the geophone nodes to 7 km between the broadband stations (details in Polychronopoulou et al., 2018; Lehujeur and Chevrot, 2020). The acquisition covers the Mauléon and Arzacq Basins with 191 geophone nodes, 197 short period instruments, and 54 broadband stations. Despite being widely covered by seismic reflection profiles and wells, thanks to decades of oil and gas exploration, the 3-D architecture of the Mauléon-Arzacq Rift System remains poorly constrained.

Different rift phases are regionally recognized during the Mesozoic, but it is well accepted that the north Pyrenean Rift Basins are mostly shaped by the Aptian-Cenomanian hyper-extension rifting stage before being shortened by the Pyrenean orogenesis from Late Cretaceous onwards. The structural complexity resulting from this polyphasic deformation history probably explains why active deep reflection seismic techniques failed to image even the first-order geometry of the basins (Daignières et al., 1994). This area is also characterized by a major positive gravity anomaly regionally known as the Labourd anomaly (Casas et al., 1997; Chevrot et al., 2014) (Figs. 2A and 2C), which has been recently related to the presence of a mantle body in the basement of the northern part of the Mauléon Basin (e.g., Jammes et al., 2009; Masini et al., 2014; Wang et al., 2016; Saspiturry et al., 2019b; Lescoutre and Manatschal, 2020) (Fig. 2D). In spite of its importance to understand the rifting and orogenic evolutions, the 3-D morphology of the continental Moho beneath the study area remains disputed (Daignières et al., 1994; Chevrot et al., 2015; Wang et al., 2016; Teixell et al., 2018).

In this work, we exploit phase velocity maps of fundamental mode Rayleigh waves at periods between 2 and 9 s in order to invert a 3-D shear wave velocity (VS) model of the area. These phase velocity maps were obtained by extracting coherent surface wave fronts from seismic noise originating from the Atlantic and Mediterranean sea using the Maupasacq array, as described in Lehujeur and Chevrot (2020). The method used an innovative matched-filtering technique to exploit the large number of seismological stations deployed simultaneously on the field and the strong directivity of ambient seismic noise. The present contribution can be seen as a continuation of that work as we invert the previously obtained phase velocity maps to image the 3-D variations of shear velocity beneath the Mauléon and Arzacq Basins. The inversion procedure involves a first point-wise Markov chain Monte-Carlo inversion (McMC) to determine a 1-D shear velocity model at each surface node of the 3-D tomographic grid. Then, these 1-D models are aggregated to get a preliminary 3-D model. This model is finally used to initiate a new 3-D iterative inversion, using a 3-D Gaussian model covariance matrix. This allows us to obtain a smooth 3-D model that fits the dispersion data.

The paper is organized as follows. In Section 2, we provide an overview of the geological context of the study area highlighting the main questions that remain concerning the architecture of the Mauléon Basin and his underlying basement. In Section 3, we describe the details of the two-step surface-wave inversion method. In Section 4, we confront and compare our VS model derived from surface waves to other tomographic studies covering the area, and in particular to a recent local travel-time tomography that exploited the Maupasacq dataset. We also compare these models to constrained geological information from geological maps and well data. Finally, in Section 5 after discussing the limits of passive imaging approaches, we propose a 3-D geological interpretation of the tomographic models.

Present-day Pyrenean structure

The Pyrenean Mountain Belt results from the north-south convergence between the Iberian and Eurasian plates from the Late Santonian to the Early Miocene (Puigdefàbregas and Souquet, 1986; Olivet, 1996; Rosenbaum et al., 2002; Macchiavelli et al., 2017; Van Hinsbergen et al., 2020). In its western part, the focus of this study, the belt is made of different structural units characterized by orogenic and pre-orogenic rift-related geological records (see Fig. 2B). The southern limit of the Aquitaine Foreland Basin is the North Pyrenean Frontal Thrust (NPFT), which corresponds to the Sainte-Suzanne Thrust in its westernmost part. Southward, the North Pyrenean Zone (NPZ) is separated from the Axial Zone (AZ) by the south-verging Lakhoura Thrust System. The AZ represents a large outcrop of exhumed pre-Mesozoic rocks that hold some of the highest altitudes of the Pyrenees. It is connected to the South Pyrenean Foreland Basin (SPFB) by a set of south-verging thrusts. While the AZ and the SPFB mainly result from the Eocene-Miocene antiformal nappe-stacking of pre-orogenic Iberian thick continental crust (proximal margin domains of Masini et al. (2014), Tugend et al. (2014) and Tugend et al. (2015)), the NPZ primarily results from the shortening of Cretaceous hyper-extended rift basins, with a basement-cover decoupling in the Upper Triassic evaporite layers. In the Western Pyrenees, the NPZ results from the inversion of the Mauléon Basin that has been thrusted on top of the Aquitaine Foreland Basin. At the front of the NPFT, and buried underneath syn-orogenic sediments, the Arzacq Rift Basin corresponds to the northern extension of the Cretaceous Pyrenean Rift System that escaped from most of the orogenic overprint (Masini et al., 2014; Tugend et al., 2015; Angrand et al., 2018; Issautier et al., 2020; Ducoux et al., 2021).

Tectono-stratigraphic evolution of the Western Pyrenees

In the Western Pyrenees, the crustal basement is made of a complex lithological assemblage of pre-Mesozoic rocks. It includes metamorphic and sedimentary rocks that were variably affected by Variscan orogenic to post-orogenic deformations and intrusions (e.g., Saspiturry et al., 2019a). This complex inherited crust was eroded and peneplained before being unconformably covered by continental clastics (Permian and Lower Triassic). An ill-defined early phase of rifting (see Leleu et al., 2016, for a review) was recorded through the Triassic marine transgression that eventually led to the deposition of the Upper Triassic evaporites (Curnelle, 1983), which will represent the main décollement level during the subsequent rifting and orogenic deformations.

During the Jurassic, the Pyrenean domain recorded a quiet tectonic period through the development of a shallow marine carbonate platform (Canérot, 2017). Another ill-defined pulse of rifting, responsible for local emersion and erosion of the sedimentary cover, is regionally recorded during the Late Jurassic and Early Cretaceous in the area (e.g., Canérot et al., 1999). The last and most intense rifting phase, the so-called hyper-extension phase, then occurred between the Aptian and the Cenomanian. As the pre-orogenic crustal and basin architectures are mainly shaped by this hyper-extension phase of rifting, we will consistently refer to this latter phase by using a pre-, syn- and post-HE terminology (“HE” standing for “hyper-extension”) to avoid any confusion with preceding rifting phases. While the precise kinematics of hyper-extension is still debated, it corresponds to the main phase of crustal thinning and accelerated subsidence within the Mauléon and Arzacq Basins. Within the Maulón Basin, this hyper-extension phase eventually led to the exhumation of the subcontinental mantle at the seafloor during the Early Cenomanian (Fig. 2D; Jammes et al., 2009; Masini et al., 2014; Saspiturry et al., 2019b) and to the rework of mantle clasts within the syn-HE sediments (Fortané et al., 1986; Jammes et al., 2009; Lagabrielle et al., 2010; Debroas et al., 2010; Corre et al., 2016; Asti et al., 2019; Lagabrielle et al., 2019). Crustal thinning and tectonic exhumation of mantle rocks were proposed to be accommodated by a set of long-offset detachment faults dipping either to the north, to the south, or both, depending on the authors (Jammes et al., 2009; Masini et al., 2014; Gómez-Romeu et al., 2019; Saspiturry et al., 2019b, Lescoutre et al., 2019).

After a short post-rift relaxation until the Late Santonian, early orogenesis initiated with the inversion of the Mauléon Basin forming the early orogenic fold-and-thrust belt of the Western Pyrenees. During this stage, deformation was accommodated by the underthrusting of the former hyper-thinned continental and exhumed mantle basement of the Mauléon Basin underneath the European continental crust (Tugend et al., 2014; Mouthereau et al., 2014; Gómez-Romeu et al., 2019). At the surface, part of the pre-HE sedimentary cover together with the syn- and post-HE cover were detached from the hyper-extended basement (including Permian and Lower-Middle Triassic rocks) using the Keuper evaporites as a basal décollement level. The Lakhora and NPFT Fault Systems thrusted the deformed sedimentary cover on top of both the Iberian and European continents (Jammes et al., 2009; Teixell et al., 2016; Labaume and Teixell, 2020). Following a Paleocene transitional phase of more diffuse deformations (the “Pyrenean tectonic quiescence”) and a Lower Eocene proto-collision stage (for details, see Teixell et al., 2016; Waldner et al., 2021), shortening propagated further south from the Eocene to the Miocene and formed the Pyrenean Axial Zone and the Iberian Foreland fFold and Thrust Belt (the South Pyrenean Zone) up to the pro-foreland basin (the Ebro Basin). While the collisional flexure was also recorded on the European retro-foreland basin, shortening rates were more limited there in comparison to the southern side of the belt (e.g., Waldner et al., 2021; Teixell et al., 2016; Angrand et al., 2018). In the study area, syn-collisional deposits are lying underneath and northwards of the NPFT (i.e., in the Aquitaine/Arzacq Basins) as turbiditic deposits grading upsection to syn- and post-orogenic marine and continental clastics. At depth, shortening was accommodated by the north-directed underthrusting of the Iberian crust beneath Europe leading to the present-day imaged geometry of the crustal root (Wang et al., 2016; Chevrot et al., 2018) (Fig. 2D).

The non-cylindrical structure of the Western Pyrenees

As introduced above, the structure of the NPZ and southern Aquitaine Basin is mostly related to the inversion of the Mauléon hyper-extended rift system. 2-D geological models generally state that the thin-skinned shortening style recorded within the NPZ results from the decoupling created by the evaporitic décollement, above which the cover got passively squeezed in a pop-up structure while the basement was thickenned by thrusts (Tugend et al., 2014; Mouthereau et al., 2014; Dumont et al., 2015; Teixell et al., 2016; Labaume and Teixell, 2020). If this evolution fits with first order observations in the eastern part of the Mauléon Basin, sections and models proposed for the western side of the study area rather advocate a thick-skin crustal pop-up that sampled the autochthonous basement of the Mauléon Basin within the orogenic wedge (e.g., Jammes et al., 2009; Masini et al., 2014; Lescoutre and Manatschal, 2020; Saspiturry et al., 2020a). This pop-up structure is edged by the lateral continuations of the Lakhoura and NPFT systems in the south and north, respectively (Teixell et al., 2016; Lescoutre et al., 2021; Lescoutre and Manatschal, 2020; Saspiturry et al., 2020a; Fig. 2D). It should be noticed that the sampled basement is outcropping within the Basque Massif forming the south-western margin of the Mauléon Basin (Fig. 2B). There, it is made of Paleozoic rocks (with its pre-Keuper cover) to the south and of the Ursuya granulites to the north, the latter bearing a small body of serpentinized mantle (Boissonnas et al., 1974; Vielzeuf and Kornprobst, 1984; Masini et al., 2014; Lescoutre, 2019; Saspiturry et al., 2019a, 2019b). The pre-orogenic sediment draping of these basement rocks suggests that it corresponds to the autochtonous basement of the Mauléon Basin, made of thinned continental crust and mantle rocks (Jammes et al., 2009; Masini et al., 2014; Wang et al., 2016; Chevrot et al., 2018; Lescoutre et al., 2021; Saspiturry et al., 2020a). The Labourd positive gravity anomaly in the western part of the Mauléon Basin gives further support to this hypothesis (Figs. 2A and 2C). Although this anomaly was initially attributed to the presence of lower crustal rocks at shallow depth (Grandjean, 1994; Vacher and Souriau, 2001; Pedreira et al., 2007), it has also been interpreted as a shallow piece of dense lithospheric mantle (Casas et al., 1997; Jammes et al., 2010). This interpretation has gained support from recent works documenting P-wave velocities of 7.3 km/s below 10 km depth (Fig. 2D; Wang et al., 2016; Chevrot et al., 2018). It was also observed that this gravity anomaly is progressively attenuating from west to east (Gottis, 1972; Boillot et al., 1973, Fig. 2A), suggesting a structural and/or a lithological change accross the basin (e.g., Masini et al., 2014; Fig. 2A). Lescoutre and Manatschal (2020) recently explained this structural change as reflecting the transition between shifted Mauléon-Arzacq and Basque-Cantabrian Rift Axes located north-east and south-west, respectively, of the Basque Massif (Fig. 2B). In this model the N20° orogen-perpendicular Saison and Barlanès transverse structures (TS; Rat, 1988; Razin, 1989; Masini et al., 2014; Saspiturry, 2019; Lescoutre et al., 2021; Fig. 2B) played a significant role on the non-cylindrical structure of the Mauléon Basin. These structures divide the study area in 3 segments from west to east: the western and eastern Mauléon segments, and the Chaînons Béarnais segment (Fig. 2B). The Saison structure is characterized by a steep dipping dextral shear zone (Zolnaï, 1975; Schoeffler, 1982; Richard, 1986; Saspiturry, 2019), whereas the Barlanès structure corresponds to the lateral termination of the Chaînons Béarnais folds. Even though they are assumed to be of fundamental importance for the closure of rifted domains during orogenesis (e.g., Lescoutre and Manatschal, 2020; Saspiturry, 2019), their precise geometry and role at depth remain to be precised, which was one of the main motivations in the choice of the emplacement of the Maupasacq seismic acquisition used in this study.

Surface wave tomography

The inversion of surface wave dispersion curves to determine the vertical variations of elastic parameters (density, shear and compressional wave velocities) is a classical non-linear inverse problem commonly encountered in earthquake or ambient noise tomography. The forward problem can be efficiently computed using the propagator matrix method, which considers a superposition of homogeneous layers over a half-space (e.g., Thomson, 1950; Haskell, 1953; Knopoff, 1964; Gilbert and Backus, 1966; Aki and Richards, 2002). The linearization of the problem has been extensively used to invert surface wave dispersion measurements because it is computationally efficient and can rapidly converge to a solution (e.g., Dorman and Ewing, 1962; Xia et al., 1999). However, the inversion results may depend on the initial model. Due to the non-uniqueness of the solution, the lack of prior knowledge about the structures can therefore lead to very different models that fit the observations equally well (e.g., Bodet et al., 2005). Global search approaches such as grid search (e.g., Macquet et al., 2014), Monte-Carlo methods (e.g., Shapiro and Ritzwoller, 2002; Socco and Boiero, 2008; Maraschini and Foti, 2010), genetic algorithms (e.g., Lomax and Snieder, 1994), or neighborhood algorithms (Sambridge, 1999; Mordret et al., 2014) are thus often preferred since they only require solving a forward problem. In addition, they can also provide the posterior covariance matrix of the model.

Building a 3-D S-wave velocity model based on the phase velocity of surface waves is classically done by first estimating the local dispersion curves at each geographical location and then solving many independent 1-D inverse problems (e.g., Shapiro and Ritzwoller, 2002). Other studies have proposed to integrate these two inversion problems into one to determine directly the 3-D structure either by linearized or Monte-Carlo methods (Fang et al., 2015; Zhang et al., 2018). These techniques still rely on the 1-D approximation for the forward computation of dispersion curves, but they allow controlling the lateral coherence of the velocity parameters.

In this work, we employ another two-step inversion scheme to invert a set of surface wave dispersion maps at several periods, and we apply it to the Rayleigh waves phase velocity dispersion maps from Lehujeur and Chevrot (2020). In the first step, we invert the surface wave dispersion curves at each geographical location independently using a 1-D Markov chain Monte-Carlo inversion. This step provides us a probability density function of the S-wave velocity as a function of depth at each location. The median profiles are combined to form a preliminary 3-D S-wave velocity model. In the second step, we perform a 3-D linearized inversion starting from the result of the point-wise inversion, and we regularize the problem using a non-diagonal covariance matrix over the model space to smooth the model laterally and vertically, while preserving the fit to the dispersion data.

Preliminary point-wise depth inversion

We first determine a preliminary model following the inversion approach described in Lehujeur et al. (2018) except for some aspects that are detailed below. We invert the phase velocity dispersion curves of fundamental mode Rayleigh waves sampled logarithmically at 9 periods between 2 and 9 s on the surface nodes of the 3-D tomographic grid (Lehujeur and Chevrot, 2020). We assume that the phase velocity measurements at the 9 periods are independent. The 1-D models are described by the shear velocity (VS) and thickness of 8 homogeneous layers overlying a half-space. The VP and density inside each layer, needed for the resolution of the forward problem are derived from VS, using the VP-to-VS and VP-to-density scaling relationships given by Brocher (2005):
VP=0.9409+2.0947VS0.8206VS2+0.2683   VS30.0251VS4,
graphic
(1)
ρ=1.6612VP0.4721VP2+0.0671VP30.0043VP4+0.000106VP5.
graphic
(2)

These relations are valid for VS ≤ 4.5 km/s (Eq. (1)), and for 1.5 ≤ VP ≤ 8.5 km/s (Eq. (2)). Here, we used them mostly to reduce the size of the model space to be explored by the Markov chain Monte-Carlo algorithm, which reduces significantly the computational effort required by the first step of our inversion procedure and does not impact the solution in a significant way as the Rayleigh waves are mostly sensitive to VS (e.g., Xia et al., 1999). The dispersion curves are computed with the codes from Herrmann (2013). All the inverted parameters are restricted to a prescribed range using a uniform probability density law, which reduces even more the size of the model space to explore and avoids testing unrealistic models (see the area delimited with black dotted lines in Fig. 3). In addition, we impose that the VS variation between two consecutive layers are between −0.5 and +1.0 km/s. This additional constraint forces the tested models to remain relatively smooth, while still allowing the velocity to decrease with depth.

The inversions are run using 10 independent Markov chains launched in parallel. The proposal distribution is adjusted during the inversion to maintain an acceptance rate of the Metropolis algorithm around 25%, and the Markov chains run until 2500 models are kept, so that approximately 10 000 models are tested. We retain the median of the 1000 best models found as the solution of the inversion (Fig. 3). After combining all the median 1-D models, we obtained the 3-D S-wave velocity model shown in Figures 4 and 5. The solutions at nearby locations are very consistent although the inverse problems are solved independently, which emphasizes the stability of the point-wise approach. However, the small-scale irregularities suggest that the model is noisy. A simple cure would be to smooth the model, but this would degrade the fit of observed data, as discussed later.

3-D inversion

To obtain the final 3-D VS model, we perform a new linearized 3-D inversion, starting from the solution obtained with the point-wise inversion described in the previous section, following the method described in Montagner and Tanimoto (1990) for global surface-wave tomography. We regularize the problem using a non-diagonal Gaussian 3-D covariance matrix to smooth the model, while preserving the fit to the observations. The model is parameterized into a regular 34 × 23 × 30 Cartesian grid of 1.612 × 1.636 × 0.345 km size cells. As in the preliminary inversion, VP and density are scaled to VS using the relations from Brocher (2005) (Eqs. (1) and (2)). The least squares cost function to minimize is defined by Tarantola and Valette (1982):
χ2(mn)=χd2(mn)+χm2(mn),
graphic
(3)
where,
χd2(mn)=12[dobsg(mn)]T·Cd1·[dobsg(mn)],
graphic
(4)
χm2(mn)=12[mnm0]T·Cm1·[mnm0],
graphic
(5)
where dobs is the data vector that contains the dispersion measurements, g is the operator that predicts the Rayleigh wave dispersion curves, mn is the model vector that contains the values of VS inside each tomographic grid cell, n is the iteration number of the linearized inversion, and m0 is the starting VS model. We assume that all the measurements are independent and thus we use a diagonal data covariance matrix Cd. We also assume that the spatial variations of VS can be described by a Gaussian covariance matrix.
Cm=σm2S,
graphic
(6)
where σmforumla is the prior uncertainty on the VS model assumed homogeneous and tuned to regularize the inverse problem, and S the Gaussian smoothing kernel:
Sij=exp{12[(rirj)T·L1·(rirj)]}.
graphic
(7)
where ri is the location of cell i. Because we expect a finer vertical resolution, we impose two distinct correlation lengths along the horizontal and vertical dimensions:
L=[Lh2Lh2Lv2],
graphic
(8)
where the horizontal (Lh) and vertical (Lv) correlation lengths are set to 2.5 km and 1 km, respectively.
This non-linear inverse problem is solved iteratively following Tarantola and Valette (1982), where the model is updated at each iteration according to:
mn+1=m0+Cm·GnT(Cd+Gn·Cm·GnT)1·[dobsg(mn)+Gn·(mnm0)],
graphic
(9)
with Gn is the matrix containing the Fréchet derivatives computed in the current model mn:
Gn(T,z)=c(T)VS(z)|mn,
graphic
(10)
where c denotes the phase velocity, VS is the S-wave velocity and T and z are the period and depth, respectively.
The inversion is initiated with m0, a smoothed version of the model mpriorobtained in the preliminary inversion:
m0=S¯·mprior,
graphic
(11)
where S¯ij=Sij/kSikforumla is the normalized smoothing kernel. We run the second inversion step for several values of the prior uncertainty σmforumla, used as a damping parameter. The inversion usually converges after 3 to 10 iterations depending on the damping coefficient used. The retained prior uncertainty is σmforumla = 0.4 km/s, slightly below the maximum curvature of the L-curve computed following Hansen and O’Leary (1993) (see the red dashed curve in Fig. 6).

The data misfit (first term of Eq. (3)) obtained with the prior model resulting from point-wise inversion (step 1) is χd2(mprior)=169forumla in arbitrary units (Fig. 6, blue dashed line). The starting model from equation (11) has a much larger data misfit χd2(m0)=465forumla, which demonstrates that simply smoothing the output of a point-wise inversion deteriorates the fit to dispersion data. After convergence of the second inversion, the data misfit for σmforumla = 0.4 km/s is χd2(msol)=124forumla (Fig. 6), which is lower than the misfit obtained with the prior model. Examples of observed and modeled dispersion curves for the prior, starting, and optimized models are shown in Figure 7.

Local earthquake tomography

In order to provide independent constraints on the S-wave velocity structure, we also show the results of a VS model obtained from local earthquake tomography (Villaseñor et al., 2019), using arrival times picked on records from the stations of the Maupasacq experiment. During the 6 months of operation of this temporary network, a total of 1980 local earthquakes were detected and located (Polychronopoulou et al., 2018). From this catalog, we selected for the tomography well recorded earthquakes with at least 8 P arrival times, 2 S arrival times and an azimuthal gap smaller than 200°. This resulted in a dataset consisting of 996 earthquakes with a total of 87 122 P-wave and 72 445 S-wave arrival times.

The tomographic inversion method used is based on Benz et al. (1996) as modified by Tryggvason et al. (2002) in order to incorporate P- and S-wave arrival times. The travel times (forward problem) are calculated using the finite-difference code of Podvin and Lecomte (1991), which provides accurate results even for models with large lateral variations in velocity structure. The tomographic method inverts simultaneously for P- and S-wave velocity structure and earthquake relocation. This is a non-linear inverse problem so in order to find its solution it is first linearized and then solved iteratively. The model is parameterized in terms of constant velocity cells of 1.2 × 1.2 × 1.2 km in size, and the dimensions of the model region are 130 km and 115 km along the E-W and N-S directions respectively.

The starting model used in the inversion is the 1-D velocity model previously determined for earthquake location (Polychronopoulou et al., 2018). The final model shown here was obtained after 10 iterations. The initial root-mean square (RMS) error of the used arrival times was 0.107 s for P-waves and 0.101 s for S-waves, and the final RMS was 0.082 s for P-waves and 0.081 s for S-waves. This represents a variance reduction of 23% for P-waves and 20% for S-waves, which is typical for travel time tomography. We estimated the resolving power of the dataset performing reconstruction tests of synthetic checkerboard models (Supplementary Material S3). Due to the high density of stations, we were able to resolve reliably anomalies of 5 km in size and larger in the center of the LET model (southern part of the Maupasacq deployment), where the high velocity anomalies of interest are located.

The S-wave velocity model obtained using local earthquake tomography shows finer details than the one obtained from SWT in the regions that are well sampled by crossing rays. However large regions of the model that are not illuminated by rays are unconstrained. This problem is less significant in surface wave tomography because of the characteristics of the sensitivity of surface waves. Similarly the actual value of the velocity anomalies obtained from local earthquake tomography is typically overestimated, while the values obtained from surface waves are more realistic. Therefore, combining the complementary characteristics of both methods we can better describe the features of the seismic velocity structure.

The final VS model obtained after inversion of the surface wave dispersion curves is shown in Figures 8 and 9. As expected, this model is smoother than the prior model shown in Figure 4 but the dominant structures are very similar. Figure 10 displays VS at 1 km depth in our surface-wave tomography (SWT) model and in the local earthquake tomography (LET) VS model described here and in Villaseñor et al. (2019). Whereas these two models were derived from completely independent datasets, they show very similar and coherent velocity structures. For example, both models image higher than average velocities in the Mauléon Basin (SW part in the map shown in Fig. 10). However, the average velocity in the LET model is clearly on average faster than in the SWT model. Figures 1113 show cross-sections in the SWT and LET models. These vertical cross-sections are also in good agreement, but the vertical resolution is finer in the SWT model than in the LET model, especially in the area of the Arzacq Basin (Fig. 11). For example, the expected top of the basement of the Arzacq Basin (ca. 3 km in depth on the northern part of the Fig. 11) is more sharply defined in the SWT model whereas the LET model shows a more gradual transition. In contrast, the LET model has a finer lateral resolution compared to the SWT model, especially at depths larger than 6 km.

Both models show anomalies that are in good agreement with surface geology. For example, the elliptical low velocity anomaly observed in both models in the central part of the Mauléon Basin and elongated along the WNW-ESE direction (y1 in Figs. 10, 12 and 13) i.e., parallel to the strike of the main orogenic structures, corresponds to the core of a mapped syncline. The core of this syncline is filled by Upper Cretaceous turbidites with Lower Cretaceous and Jurassic marls and limestones underneath. Conversely, the Saint-Palais anticline in the central part of the C1-C3 section is characterized by higher seismic velocities (y2 in Figs. 10 and 13). The low velocity anomaly y3 in Figures 10 and 11, seen in both models, is related to the Cenozoic sediments accumulated in the Aquitaine Basin. On map view, the regional shape of the NPFT, which juxtaposes Mesozoic rocks onto younger sediments and makes a lateral ramp north of the Saison structure, has a clear signature at shallow depth in both models (y4 in Fig. 10). Finally, the higher velocities in the southern part of the models can be related to the metamorphic rocks of the Arbailles and Labourd Paleozoic Massifs (Fig. 10).

In Figures 1416, we compare our tomographic models with seismic reflection profiles, surface geology and well data from some key boreholes of the Mauléon Basin (details in Supplementary Materials S1 and S2). In order to facilitate these comparisons, we converted our models from depth to two-way travel times using
twt(x,z)=2[ztopozdyVP(x,y)+ztopozsrdVrep],
graphic
(12)
where twt(x,z) is the two-way travel time corresponding to depth z and position x along the profile, VP the P-wave velocity estimated by multiplying the VS by a VP/VS ratio of 1.73, zsrd the reference level taken at 500 m above the see level, ztopo the topography level, and Vrep the replacement velocity chosen at 2.5 km/s. The time adjustment related to zsrd and Vrep is used to avoid propagating the irregular surface topography in the twt profile. It consists in filling the volume between the topography and the reference level by an imaginary layer of constant velocity. The values chosen for Vrep and zsrd have little impact on the resulting figures.

In the central part of the western section (Fig. 15), a velocity inversion is observed at ~ 2 s TWT in the SWT model (Fig. 15A). This inversion is in good agreement with the Bellevue well, which documents Jurassic marls and limestones over the imbricated Lower Cretaceous limestones and Upper Triassic evaporites (Fig. S1). This tectonic contact probably corresponds to the Bellevue Thrust that is cropping out at the surface further north. A similar velocity inversion can be recognized in the northern part of the SWT model on the western section at 2.5 s TWT. This depth corresponds to the occurrence of the Sainte-Suzanne Thrust in the Orthez well that brings shallow high-velocity rocks (Lower Cretaceous and Jurassic carbonates) on top of younger clastics sediments (Fig. S1). Further south, a shallow and smooth vertical velocity inversion can also be observed north of the Arbailles Massif, where the nearby Ainhice well documented a duplication of the stratigraphy at depth (Lescoutre et al., 2019, 2021; Saspiturry et al., 2019b; Fig. S1). This structure could correspond to a thin-skin thrust at the front of the Arbailles Massif (Fig. 15).

Limits of passive imaging methods

Part of the resolution discrepancy observed between SWT and LET may be explained from the different parameterization and regularization schemes used in the two types of inversions. However, because velocity anomalies tend to be smeared along the propagation directions, surface waves better constrain vertical structural variations compared to horizontal ones, whereas for body waves it is the opposite. Vertical smearing is indeed strongly expressed in the LET model, particularly at shallow depth. This may partly explain the smaller amplitudes of velocity anomalies in the upper crust in this model. The comparison of LET and SWT also reveals that whereas the shallow velocity anomalies in the two models are rather well correlated, the LET model is on average faster than the SW model. This suggests that the starting 1-D model used in the LET inversion is on average too fast, a bias that remains in the final 3-D model. This can simply be understood by the trade-off between the average velocity model and the origin time of the earthquakes. Therefore, the initial model used in LET still has a strong imprint on the final 3-D model, and it is thus crucial to build this model very carefully, using all the a priori information that is available. SWT could provide such key constraints on absolute velocities, in particular at shallow levels.

Other important differences come from the intrinsic limitations of both approaches. In LET, the resolution is controlled by ray coverage, and thus by the distribution of earthquakes. In the Western Pyrenees, the seismicity is concentrated inside a narrow band that approximately follows the limit between the Axial Zone and the North Pyrenean Zone, with hypocentral depths rarely exceeding 20 km (Souriau and Pauchet, 1998b; Chevrot et al., 2011). Owing to the distribution of earthquakes, the resolution is thus limited to the top 15 km, but it degrades notably beneath the Arzacq Basin at depths below 10 km. In contrast, the main limitation of SWT comes from the amplitude and distribution of noise sources, the geometry of the array, and the sensor types. In this study, we exploited surface waves excited by energetic oceanic sources at periods from 2 to 9 s (microseismic band). At longer period, the sensitivities of geophones and short period sensors become very low, and the records are dominated by instrumental noise. Therefore, at periods longer than ~ 6 s only the broadband sensors can be exploited. In addition, the wavelengths of Rayleigh waves at these periods are larger than 20 km i.e., of the order of the size of the region that we want to image. This may explain the degradation of the lateral resolution in the deeper part of the model, especially beneath the NW-SE Mauléon transect shown in Figure 12. On the other hand, at shorter period the station spacing (≥ 1 km) becomes too coarse to sample the wavefield with at least 2 samples per wavelength, and cycle skipping issues appear. For these reasons, the resolution is more uniform in the SWT model but limited to the top 10 km of the crust.

In any case, the similarity of LET and SWT models, in addition to the different but complementary sensitivity of body and surface waves to shallow structures, suggest that joint inversions would be a natural way to further improve the resolution and robustness of crustal tomographic models. Because these inverse problems are still solved separately with different model parameterization and regularization schemes, this will require developing a new generation of inversion codes.

Geological interpretations and integration into a 3-D structural model

We now attempt to build a coherent 3-D model of the Mauléon Basin. We base our interpretation mostly on the shape of the velocity anomalies observed in the new tomographic models and how these anomalies agree or not with the available geological information, gravity data and seismic profiles, keeping in mind the strengths and weaknesses of each approach discussed above. As already pointed out in the previous section, the velocity inversions can be associated with major thrusts (see Sect. 4; Fig. 15). The inversion seen on the eastern section (Fig. 14) can thus be interpreted as the southward continuation of the south-dipping NPFT beneath the bottom of the Les-Cassières-2 borehole (Fig. S2). The western profile is located to the west of the Saison structure (Fig. 2A), where the positive Bouguer gravity anomaly interpreted as related to a shallow body of subcontinental mantle has been imaged by Wang et al. (2016). This piece of mantle is located in the hanging-wall of the Sainte-Suzanne Thrust (Fig. 15). East of the Saison structure, the high velocity mantle body is laterally shifted to the south and replaced by lower velocity rocks in the hanging wall of the steeper NPFT. These rocks might correspond to middle/lower crustal (metamorphic) rocks, serpentinized mantle or sedimentary rocks. The geometry of the iso-velocity lines at depth, in both the eastern and western NS sections, suggests that the high-velocity rocks (VS> 3.8 km/s) get shallower towards the north underneath the Mauléon Basin, and thus that upper crustal rocks thin at the expense of deeper lithospheric rocks. Note that similar relationships between basement velocities and crustal tapering were reported from offshore Ocean-Continent transitions such as the Iberian (Dean et al., 2000) or the Newfoundland (Lau et al., 2006) margins. This interpretation is in line with previous models from the Mauléon Basin involving north-dipping detachment faults exhuming granulite and subcontinental mantle rocks (Jammes et al., 2009; Masini et al., 2014; Gómez-Romeu et al., 2019; Saspiturry et al., 2019b; Lescoutre et al., 2019).

In the along-strike section (Fig. 16), the 3.2–3.6 km/s iso-velocity lines become shallower at the western edge of the LET profile where the basement rocks of the Labourd Massif are cropping out (Fig. 2B). The trend of these iso-velocity lines indicate the geometry of the top basement at depth. At the eastern edge of the E-W section, the high velocities at shallow depth are in good agreement with the nappe-stack of the Chaînons Béarnais that consists of mixed sedimentary cover (including metamorphosed carbonates) and basement rocks (including serpentinized mantle bodies) (Labaume and Teixell, 2020). At depth, in the central part of the model, a steep velocity contrast between high velocities to the west and lower velocities to the east can be observed on the LET model (Fig. 16). In the easternmost part of the E-W profile, a similar observation can be made with an abrupt change towards lower velocities to the east. These steep velocity gradients correspond to the location at depth of the Saison and Barlanès transverse structures mapped at the surface (Fig. 2B). Therefore, these structures are segmenting the Mauléon Basin and the Chaînons Béarnais in a roughly orogen-perpendicular direction. As seismic velocities are to first-order positively correlated to densities, this observation is consistent with the eastward attenuation of the positive gravity anomaly attributed to shallow mantle rocks across the Saison and Barlanès structures (Figs. 2A and 12). The exact nature of the basement rocks and the depth of the Moho on both sides of the Saison structure remain poorly constrained. However, as stated above, the tomographic model of Wang et al. (2016) suggests the presence of subcontinental mantle at shallow depth (~ 10 km) below the western Mauléon segment. Because the LET model shows relatively lower velocities from the west to the east of the Saison structure, and the gravity anomaly shows a similar transition (Fig. 12), the basement of the eastern Mauléon segment is most likely made of less dense material more altered/serpentinized (e.g., mantle rocks or mafic/metamorphic lower crust). This interpretation is in good agreement with previous structural models that interpreted the Saison structure as a crustal transfer zone active during the Pyrenean orogeny (Masini et al., 2014; Lescoutre and Manatschal, 2020; Saspiturry, 2019). The Saison and Barlanès structures accommodated deformation between the western Mauléon segment, where the hyperextended rift domain together with the subcontinental mantle rocks have been transported in the hanging wall of a north-vergent thrust, and the Chaînons Béarnais segment, where most of this hyperextended domain has been underthrusted beneath the European crust (Fig. 17). This distributed deformation also suggests a change in the structural level of indentation during the orogeny. While the European crust anomalously indented the Iberian mantle at an early stage of convergence in the west leading to its shallow sampling (i.e., thick-skinned style, see Lescoutre and Manatschal, 2020), it is likely that the indentation rather used the basement-sediment interface in the east (Fig. 17). This scenario explains the apparent increase of thin-skinned shortening eastward, because of a significant accommodation of shortening within the basement in the west. Note that these complex 3-D structures are restricted to the hanging-wall of the north-dipping “slab” previously imaged by Wang et al. (2016) that should be consistently made of the Iberian basement formerly located in the southern border of the basin.

Our study demonstrates that, using dense large-N deployments, it is possible to obtain finely resolved images of fold and thrust belts from the exploitation of surface and body waves with passive imaging approaches. Obviously, the level of details in our tomographic images is not on par with a typical seismic reflection survey and resolving the different sedimentary horizons inside a sedimentary basin is clearly beyond the reach of surface wave tomography. Nevertheless, our study provides robust first-order constraints on the deep architecture of the Mauléon and Arzacq Basins, where previous controlled source seismic reflection studies actually gave rather poor results, and for a fraction of the cost of such acquisitions. In particular, our tomographic models confirm the presence of orogen-perpendicular structures in the study area that controlled the along-strike change in the orogenic basement structure and composition. They also control the local preservation of hyper-extended rifted crust and mantle at shallow depth going along the gravimetric anomaly. We are thus convinced that passive imaging represents a valuable source of information that should be considered in the future, especially in environments where controlled source acquisitions are challenging or are impeded by legislation.

This work was supported by OROGEN, a tripartite research project between the CNRS, TOTAL and BRGM, and by the ANR AAPG program (project CLEARVIEW, ANR-17-CE23-0022). AV received funding from the Spanish government through the “Severo Ochoa Centre of Excellence” accreditation (CEX2019-000928-S). We thank associate editor David Pedreira, as well as Josep Antón Muñoz and an anonymous reviewer for their constructive comments.

Cite this article as: Lehujeur M, Chevrot S, Villaseñor A, Masini E, Saspiturry N, Lescoutre R, Sylvander M. 2021. Three-dimensional shear velocity structure of the Mauléon and Arzacq Basins (Western Pyrenees), BSGF - Earth Sciences Bulletin 192: 47.

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.